香京julia种子在线播放

We used a Wilson-Cowan model of the mean activity of seven basal ganglia nuclei, namely, D1-MSN (D1-medium spiny neuron), D2-MSN, FSI (fast spiking interneuron), GPe-TA (globus pallidus externus—tonically active), GPe-TI (GPe-tonically inactive), STN (subthalamic nucleus) and GPi (globus pallidus internus). The dynamics of the system is given by

with the population-firing-rate vector

and the sigmoidal activation function

Here, λ_CTX denotes the firing rate of cortical inputs, Ẏ=dYdt=[dλD1dt,dλD2dt,dλFSIdt,dλTAdt,dλTIdt,dλSTNdt,dλGPidt]T the temporal derivative of the rate vector Y(t), and ^T the matrix transpose. For the sake of simplicity, the time constant τ = 15 ms is kept constant for all nuclei. To test the sensitivity of the results to this choice, parameter search and simulations are also performed for τ = 1ms. This produces qualitatively similar results, but shifts the frequency of beta oscillations (peaks in spectra in Figure 6D) out of the experimentally observed range (data not shown).

The values of θ and λ_max (Table 1) were chosen in order to get realistic distributions of instantaneous firing rates under different input conditions (Supplementary Figure 2). For a specific nucleus, λ_max and θ are fixed across all network configurations and for both physiological and pathological conditions. Activity is modeled assuming zero delays (i.e., instantaneous update), as this is convenient for the overall ease of analysis of the system. Moreover, oscillations in a circuit can be purely delay driven. As we are interested in the role of effective connectivities on oscillations, a choice of zero delay removes the confounding factor.

The coupling matrices

and

denote the recurrent and input connection matrices, respectively. An element J_i,j of A or B represents the effective strength of the connection from population j to population i. The colors of the parameters in Equations (3) and (4) indicate whether a particular projection is considered as a free (red) or fixed (blue) parameter. The calculation of the fixed parameters is described in detail in Section 2.1.1. The code for this model and analysis scripts are shared in the git repository https://github.com/jyotikab/Homology_BG.git.

We used a genetic algorithm to find valid network configurations that reproduce experimental observed phase and activity relationships in physiological and Parkinsonian conditions. The pseudocode is given in Algorithm 1.

The crossover number X (number of elements among the 20 free parameters to be swapped) and the random mutation probability Y (probability of a random mutation) were tuned in order to cover a reasonable amount of parameter space in a time span of 3 days (≈3,000 iterations). A value of X = 2 and Y = 0.1 gave the best results (data not shown). The range for all the negative parameters was (0, −6.) and positive parameters was (0, 13). In the case of a random mutation, the new value was drawn from the corresponding positive or negative uniform distribution. The parameter search was run for three different seeds to ensure that the results didn't depend on the initial conditions. The distributions are shown in Supplementary Figure 1.

The cost function used in the parameter search is binary rather than real-valued. We derive a series of phase and activity relationships based on experimental findings on the response of the GPe-nuclei (arkypallidal and prototypical) to cortical slow wave activity (SWA) and beta activation (Beta) (Mallet et al., 2008; Abdi et al., 2015). The SWA was recorded under the influence of anesthesia, whereas the cortical activity showed beta activation during hindpaw stimulation of the animal. For each candidate network configuration, we set λ_CTX in Equation (1) to provide SWA and beta stimulus as defined in Table 3. We then test whether the network's activity fulfills the relationships set out in Table 4. If it fulfills all the criteria in the left column, the network configuration is considered a valid example of a physiological network. If it fulfills all the criteria in the right column, it is considered a valid example of a parkinsonian network. If neither, than the configuration is discarded (see Algorithm 1).

The (Pearson's) cross-correlation coefficients Corr(λ_x,λ_y) in Table 4 are used to test phase relationships between the activities of the nuclei x and y. A positive cross-correlation between the activities of STN (λ_STN) and GPe-TA (λ_TA), for example, imply an in-phase relationship. Although cross-correlation is not a direct test for phase relationship, it is sufficient in this case since the input is sinusoidal. The phase differences between λ_TA/TI and λ_STN in Figures 2C,D are calculated from the Fourier transformed signals (using Fast Fourier Transform of the NumPy library) at stimulus frequency. For physiological networks, the experimental data suggests that GPe-TA and GPe-TI are largely non-modulated by the cortical activity (Mallet et al., 2008). This was imposed by selecting the networks with low Fano factor

for the GPe-TA and GPe-TI activity. Here, Var(·) and Mean(·) correspond to the variance across time and the time average, respectively. For parkinsonian conditions, the experimental data shows that λ_TA and λ_TI are strongly modulated by cortical activity which is reflected in a significant (positive/negative) cross-correlation and large Fano factors.

GPi suppression (GS) measures the ability of the network dynamics to effectively suppress the GPi activity and is defined as the ratio of the change in average GPi activity after the stimulus onset t^* to the activity before stimulus onset, i.e.,

where Δ₁ = Δ₂ = 500 ms.

Susceptibility to oscillations (SO) denotes the tendency of the system to oscillate under a transient stimulus (here a square wave). It is calculated as one minus the mean spectral energy, i.e.,

where i ∈ {GPi, GPe-TA, STN, GPe-TI} and SE_i denotes the spectral entropy for the amplitude spectrum Λ_i(ω) of nucleus i, calculated as the Fourier transform of the activity, λ_i.

We measure these dynamical properties of the system following a perturbation with a single square pulse of amplitude 4 spikes/s. Analogous results were obtained for amplitudes in the range 1–8 spikes/s (data not shown). Note that the square pulse input contains components in all frequency bands and is therefore an appropriate test stimulus for susceptibility to oscillations.

The distributions of values generated for the free effective connectivity parameters are manipulated in two ways, “replace by mean” and “shuffle.” The manipulations are carried out either individually (one parameter at a time, restoring the original values of a parameter after each manipulation) or cumulatively (manipulating each parameter in succession, without restoring the original values of previously considered parameters).

All network simulations and the genetic algorithm are implemented in Python (http://www.python.org). Results were analyzed using SciPy and NumPy libraries and visualizations were carried out using Matplotlib (Hunter, 2007).

We developed a firing rate model of the basal ganglia consisting of seven neuronal populations, namely, D1-MSN(D1-medium spiny neuron), D2-MSN, FSI (fast spiking interneuron), GPe-TA (globus pallidus externus - tonically active, or arkypallidal), GPe-TI (GPe-tonically inactive, or prototypical), STN (subthalamic nucleus) and GPi (globus pallidus internus) as shown in Figure 1. These nuclei were chosen with the purpose of modeling a minimal basal ganglia circuit implementing the three main functional pathways, i.e., the direct pathway (D1 → GPi), the indirect pathway (D2 → GPe → GPi) and the hyperdirect pathway (STN → GPi). For the sake of simplicity, we limit our analysis to the basal ganglia dynamics and its output. Consequently, other connected nuclei such as thalamus are not represented, and cortex is modeled as a feedforward excitatory input. Consistent with the recent experimental findings (Mallet et al., 2012; Abdi et al., 2015) we modeled GPe as a network of two subpopulations. These subpopulations are distinct in terms of their response to cortical activity and may be part of different functional pathways. Since the strengths of the majority of the projections emanating from and projecting to GPe-TA and GPe-TI are unknown, they are considered to be free parameters (marked as red dashed lines in Figure 1). These include projections to striatum (D2-MSNs, D1-MSNs and FSIs), projections from D2-MSNs to GPe and recurrent projections between GPe and STN. Although it has been suggested that GPe-TA populations (arkypallidal) projects upstream to striatum much more than GPe-TI population (prototypical), here both are included as free parameters in order to include the possibility of few but strong projections from prototypical neurons to striatum. Similarly, the projections from both arkypallidal and prototypical neurons to STN are considered, even though it is assumed that the prototypical population projects downstream. All the cortical projections, i.e., to the striatum as well as to STN are also considered as free parameters, since their relative effective strengths are also unknown. The model dynamics and structure, and the estimation of the fixed connectivity parameters (blue solid lines in Figure 1), are defined in Section 2.1 in Materials and Methods.

The firing rate model described in Sections 2.1 and 3.1 and is incompletely specified due to the large number of unknown connectivity strengths. To address this issue, rather than making assumptions to specify these, we generate a large ensemble of model configurations by performing a parameter search across the 20 free parameters (see matrices A and B in Section 2.1). A genetic algorithm is used to find configurations of the free parameters shown in Figure 1 that conform to a series of criteria based on their activity and phase relationships (see Table 4 in Materials and Methods). A configuration is considered valid if the firing activity and phase relationships of arkypallidal and prototypical GPe subpopulations fit the experimental observations made for healthy and 6OHDA-lesioned rats (i.e., modeling Parkinson's disease) presented in Abdi et al. (2015) and Mallet et al. (2008).

The method is described in Section 2.2 and resulted in 1,214 network configurations classified as “physiological” and 1,265 classified as “parkinsonian.” The corresponding activity and phase relationships are shown in Figure 2 (c.f. Figures 6 and 10 in Abdi et al., 2015). In the physiological condition, the beta activation in cortex tends to increase the average activity in both arkypallidal and prototypical subpopulations (Figure 2A) and both populations mostly fire in-phase with the STN and cortical activity (Figure 2B). In the parkinsonian condition, however, beta activation has an opposite effect on the subpopulations. The activity in prototypical neurons decreases, whereas the activity in arkypallidal neurons increases, when the cortical activity switches from SWA (slow wave activity) to beta activation (Figure 2C). The subpopulations also differ in their phase relationship with STN activity in the parkinsonian condition, with arkypallidal neurons firing in-phase and prototypical neurons firing anti-phase with STN activity (Figure 2D).

It should be noted that every data point in Figures 2A,C indicates the activity of the arkypallidal and prototypical populations in one example network (i.e., one valid parameter configuration that fits the criteria). In contrast, each data point shown in Figures 6 and 10 in Abdi et al. (2015) corresponds to a different recorded neuron. Since this experimental data was gathered from many animals (20 healthy and 16 6OHDA-lesioned), these neurons may well belong to different homologous networks. It is also noteworthy that the model results are consistent with experimental results even in an aspect not used as a validity criterion (Table 4), namely that the phase differences in physiological networks show a wider distribution (with a peak ≈ 0°, but some points showing a phase difference of 270° or 180°), as compared to the narrow distribution of phase differences in parkinsonian networks (with peak ≈ 180° for prototypical and peak ≈ 0° for arkypallidal neurons). This might be due to relatively looser constraints on the phase difference in the heathy condition, where only weak modulation of arkypallidal and prototypical neurons by cortical activity is required (see Section 2.2).

Figure 3 shows the response of the seven neuronal populations in the model to slow wave and beta activity in the cortex. The activity is shown for four randomly selected networks from both the physiological and parkinsonian classes. In the physiological condition (Figures 3A,C), both GPe-TA and GPe-TI show less modulation by cortical activity than in parkinsonian networks (Figures 3B,D), in which the opposite phase relationships between STN-arkypallidal (GPe-TA) and STN-prototypical (GPe-TI) can clearly be seen. Moreover, the parkinsonian networks exhibit an increase in GPi activity when cortical activity is high, whereas GPi mostly shows a decrease in its activity during high cortical activity in physiological networks, suggesting an anti-phase relationship between GPi and cortex (CTX). The parkinsonian networks also show a high frequency oscillation over the imposed slow wave activity unlike the physiological networks (Figures 3A,B), indicating an overall higher susceptibility to oscillations.

Table 5 gives the medians and quartiles of the effective connectivities generated for the physiological and parkinsonian ensembles. Although relationships between effective connectivities did not form part of the classification criteria, there is a good fit between many of the structural relationships observed in the generated ensembles and experimentally observed data. Some of the matches to well explored structural relationships include:

Corticostriatal projections to D1-MSNs are stronger than those to D2-MSNs in physiological networks, but the situation reverses in parkinsonian networks. This is in agreement with the general trend for Parkinson's disease, i.e., the indirect pathway becomes stronger than the direct pathway due to dopamine depletion (Flores-Barrera, 2010; Escande et al., 2016)

(physiological: J_D1,CTX = 7.8, J_D2,CTX = 2.6; parkinsonian: J_D1,CTX = 2.7, J_D2,CTX = 9.7)

Our results also predict that cortical projections to FSIs are stronger in parkinsonian networks (physiological: J_FSI,CTX = 1.9, parkinsonian: J_FSI,CTX = 8.9). This prediction, however, is in contradiction with the studies like Wiltschko et al. (2010), Clarke and Adermark (2015) which show that dopamine agonists increase FSI activity due to presence of D1 like D5 receptors.

Experimental data on the effective strengths of projections emanating from and projecting to GPe-TA-TI has been largely missing, until recently. In one such study, Glajch et al. (2016) measured IPSCs in ex vivo slices from two types of GPe neurons, one expressing NPas1⁺ and other expressing PV⁺. Neurons expressing NPas1⁺ correspond to arkypallidal (GPe-TA) neurons and those expressing PV⁺ correspond to prototypical (GPe-TI) neurons, as also shown by Abdi et al. (2015). It turns out that several relationships between the effective strengths of GPe-TA-TI connections in physiological and parkinsonian conditions observed by Glajch et al. (2016) can also be found in our generated ensembles:

The projections from GPe-TA to D2-MSNs are stronger than those to D1-MSNs. Moreover, both of these projections are stronger in in 6-OHDA conditions

(physiological: J_D2,TA = −1.2, J_D1,TA = −0.83; parkinsonian: J_D2,TA = −1.4, J_D1,TA = −0.93).

Only one prediction is in contradiction with Glajch et al. (2016), which predicts that GPe-TI projections to D1-MSNs remain unchanged when measured in naive and 6-OHDA lesioned mice, whereas our distributions show a slight decrease in GPe-TI projections to D1-MSNs in parkinsonian conditions (J_D1,TI = −0.18) as compared to physiological conditions (J_D1,TI = −0.3). Predictions of the relationships that go beyond current experimental observations are listed and discussed in Section 4.2. In total, the good fit of the generated structural relationships in our network ensembles to experimentally observed relationships from a variety of sources allow us to conclude that these structural relationships are crucial to evoke the different dynamics in the physiological and parkinsonian conditions.

The effective connection strengths generated by the genetic algorithm show a considerable variation across the ensemble of 1,214 physiological and 1,265 parkinsonian network configurations, respectively (Figure 4A and Table 5). Whereas some parameters have a narrow distribution (e.g., J_D1,TI for both categories, J_TI,TI and J_TI,STN for parkinsonian networks), others are broadly distributed (e.g., J_TA,D2 for parkinsonian networks, J_D1,CTX and J_FSI,CTX for both categories). Moreover, pairs of effective connection strengths are typically correlated (Figures 4B,C). Both in physiological and parkinsonian networks, a number of effective-weight pairs are positively correlated (e.g., {J_D1,TA, J_D1,TI}, {J_FSI,TA, J_FSI,TI}, {J_TA,STN, J_TI,STN},{J_D1,CTX, J_D2,CTX}), others exhibit negative correlations (e.g., {J_TA,TI, J_TA,TA}, {J_TI,TA, J_TI,TI},{J_TI,TA, J_TI,STN}). For most connection pairs, the magnitude of correlations is slightly smaller in parkinsonian networks. Exceptions are, for example, strong positive correlations for {J_TI,TI, J_D2,TI} and {J_TI,TI, J_D2,TA} pairs in parkinsonian networks, which are not observed in the physiological case. Further, {J_TI,STN, J_TI,D2} and {J_TI,STN, J_TI,TA} pairs are more negatively correlated in parkinsonian networks as compared to the physiological case.

For some projections, the difference between the physiological and parkinsonian distributions is striking, for example the much reduced cortical projection strength to D1-MSNs J_D1,CTX. The distributions of other connections, e.g., inhibition from GPe-TA to GPe-TI, J_TI,TA, are more similar.

To evaluate the functional significance of the observed variability and covariability in the effective connectivity, we test whether the classification of networks into “physiological” and “parkinsonian” networks according to the criteria described in Table 4 is robust with respect to (i) replacing effective connection weights by the mean for that connection across the ensemble (Figure 5A, top row), or (ii) shuffling the weights across network configurations in the same ensemble (Figure 5A, bottom row). The classification robustness is defined as the fraction of networks that retain their original classification after a given manipulation. The two manipulations are described in further detail in Section 2.4.

The replace by mean manipulation tests the assumption that the distribution of a parameter is just fluctuations around a mean and is thus not important for the dynamics of the network. The shuffle manipulation allows us to investigate whether it is sufficient for the network dynamics to preserve the marginal weight distributions across network configurations. Both types of assumptions are very common in modeling studies integrating experimentally obtained parameters, and neglect correlations between the effective weights.

For the majority of connections, replacing the weight by the corresponding ensemble average preserves the original classification for more than 50% of both physiological and parkinsonian networks. For some cases the classication robustness is very high, for instance, replacing the weight J_D1,TA of the GPe-TA inhibition to D1-MSNs by the ensemble average preserves the classification of about 95% of networks. However, for a number of connections (e.g., J_TI,STN), replacing the weight by the ensemble average causes a substantial fraction of networks (more than 50%) to no longer exhibit the dynamics described in Table 4 for their original classifications. Similar results are obtained by shuffling weights across the ensemble of network configurations (bottom panel in Figure 5A). Overall, the classification robustness is lower in this case.

The above analysis is based on replacing the weight of individual connections, leaving the other connections unaltered. It is not clear whether manipulations of one parameter affect the robustness with respect to other parameters. To address this, we repeat the analysis in a cumulative fashion, i.e., performing the manipulations for each parameter in turn without resetting the values changed in the previous manipulation. The parameters are manipulated in order of decreasing classification robustness in the individual case, as shown in Figure 5A. For the replace by mean manipulation, this order is fixed, and shown as annotations to the data points (Figure 5B, solid curves). The shuffle manipulation is carried out ten times to test the robustness of the results with respect to different random weight realizations. Potentially, each of these ten realizations has a distinct classification-robustness order. The dashed curves in Figure 5B depict averages across the ensemble of ten realizations.

With this ordering, a cumulative replacement of weights by their respective ensemble averages leads to a maximum reduction of classification robustness to about 0.3 for physiological networks and 0.5 for parkinsonian networks, although initially the physiological networks are more robust with respect to multiple manipulations. Counterintuitively, if all weights are replaced by their ensemble averages (rightmost circles in Figure 5B), the fraction of unchanged networks jumps back to 1. This suggests that the correlations observed between the parameters in Figures 4B,C are not artifactual, but are necessary for producing the dynamics of the two categories. By restricting some parameters to their means, but allowing others to maintain their distributions, some of the resulting network configurations fall outside of the volume in 20-dimensional parameter space that contains the original ensemble, and no longer exhibit the correct dynamics. As more parameters are restricted, the effect increases. However, once all parameters are set to the means of their respective distributions, the resulting network configuration is within the original volume, and thus retains the original classification. A similar effect can be seen with the shuffle manipulation, which transforms an n-dimensional ellipsoid of network configurations into an n-dimensional sphere, thus also causing many network configurations to be shifted to points outside the original volume. Here, the classification robustness decays monotonically.

To quantify the relationship between correlation and classification robustness, in Figure 5C we plot the fraction of networks that retained their classification for each parameter against the average cross correlation coefficient between that parameter and the others, calculated from the matrices shown in Figures 4B,C. A strong dependency for parkinsonian and a weak dependency for physiological networks is evident, with a line fit yielding an R-value of −0.75 for parkinsonian and −0.33 for physiological networks.

In total, these results demonstrate that there is a wide variability in the parameter configurations that pass the selection criteria and can reproduce the network activity shown in Figure 2. Our results show the dangers of restricting the parameter search to produce a single set of parameters, especially by taking the mean of distributions, unless this is carried out systematically for all distributions. Similarly, our results illustrate that the variability in parameter configurations cannot be boiled down to the marginal distributions of the individual parameters, due to presence of correlations between them.

In order to investigate the functional consequences of the variability of the generated solutions shown in Figures 4, 5, we evaluated the network configurations on the basis of their dynamical properties, identified from experimental observations on basal ganglia activity in physiological conditions and Parkinson's disease. It has been shown that insufficient GPi suppression is usually associated with stymied movement. For example, Boraud et al. (2000) showed that the ratio of inhibited-to-activated GPi neurons is significantly reduced during movement in a MPTP-treated monkey. In an another experiment, the delayed suppression of GPi neurons was shown to be correlated with slowing down of movement (Leblois et al., 2006). Akinesia, which is the loss or impairment of voluntary movements, is also associated with oscillations in the basal ganglia system. The akinetic symptoms of the Parkinson's patients have been shown to grow worse with an imposed low-frequency stimulation of 10–20 Hz in STN (Timmermann et al., 2004; Chen et al., 2007). A tACS (transcranial alternating stimulation) of 20 Hz in cortex leads to slowing of movements in healthy individuals (Pogosyan et al., 2009).

Another parkinsonian symptom associated with the oscillation frequency of basal ganglia, specifically the GPe, STN, GPi and thalamus, is tremor. Electrophysiological studies have found single units oscillating at tremor frequencies in STN and pallidum (Raz et al., 2000; Levy et al., 2002) as well as oscillations of large populations in STN oscillating at 8 − 20Hz in tremor-dominant Parkinson's patients (Moran et al., 2008). These dynamical features also reflect the effect of striatal bias downstream in the basal ganglia. A lack of striatal bias toward “Go” can contribute to insufficient suppression of GPi rates. Moreover, we previously showed in a numerical study that an excess of striatal bias toward “No-Go” can initiate oscillations in the GPe-STN circuit (Kumar et al., 2011).

We therefore evaluate the networks' responses to a transient square pulse with respect to two dynamical features, namely the effective suppression of GPi rates to a transient cortical activity (GS) and the susceptibility of basal ganglia circuit to oscillations (SO). It should be noted that this input was not used for parameter fitting in genetic algorithms to distinguish between parkinsonian and physiological networks; the response to this input is an emergent feature of the effective connectivities generated to fulfill the parkinsonian or physiological criteria given in Table 4. The GPi suppression (GS) is calculated as the normalized difference of GPi rates in presence and absence of the cortical activity, such that a value of one for GS represents efficient suppression of GPi and zero or negative values indicate ineffective suppression of GPi rates. The susceptibility to oscillations (SO) is calculated from the spectral entropy (SE). A white noise signal, for example, will show a power spectrum with no clear peak and hence high entropy, whereas for an oscillatory signal, the power will be concentrated in a certain frequency band, yielding a lower entropy. We define the susceptibility to oscillations as (1 − SE), such that values toward zero indicate low SO (high SE) and values toward one indicate high SO (low SE). A full definition of both measures and the details of the transient cortical input used as stimulation are given in Section 2.3.

An illustration of the relationship of these measures to the activity evoked in the GPi by a transient square pulse is given for some example networks in Figures 6A,B. The physiological networks (Figure 6A) show no oscillatory activity (SO = 0.03 and SO = 0.09) with either effective (brown; GS = 1.0) or ineffective suppression of GPi activity (orange; GS = −0.93). The two examples of parkinsonian networks in Figure 6B both show an increase in the GPi activity (hence ineffective GPi suppression) during the cortical transient input (shown in black), the corresponding values for GS are therefore negative (−0.9 and −0.93). One example shows strong oscillatory activity (brown; SO = 0.93) whereas the other shows no oscillations (orange; SO = 0.03). The absence and presence of oscillatory activity in physiological and parkinsonian networks, respectively, can also be seen in the amplitude spectrum plotted in Figures 6C,D. The mean normalized spectrum (thick curves) for parkinsonian networks (Figure 6D) shows a clear peak in the beta range (15–25 Hz) for all nuclei, whereas the mean spectrum for physiological networks shows no such peak (Figure 6C).

The results of projecting the network configurations classified as physiological and parkinsonian onto the dynamical feature space defined by SO and GS is given as 2D histograms in Figure 7. The physiological networks displayed in Figure 7A have a mean at SO = 0.06, GS = 0.97 (indicated by a cyan marker), which indicates low susceptibility to oscillations and sufficient suppression of GPi activity. The distribution of SO values is well fit by a gamma distribution with shape parameter k = 1.48, scale parameter θ = 0.045 and a mean of 0.06. The distribution of GS values is well fit by a power law, with parameter α = 45 with a mean of 0.97. The parkinsonian networks (Figure 7B) show a mean at SO = 0.53, GS = −0.86, indicating a much higher susceptibility to oscillations and insufficient suppression of GPi activity. Unlike the physiological networks, in this case a closer look at the 1D histograms reveals distributions of values that are well fit by a triangular distribution for SO with lower limit a = 0.07, upper limit b = 1.07 and mode c = 0.3, and a uniform distribution for GS between −1.0 and −0.75 but also exhibiting a small scattering of outliers reaching up to much higher values.

These results demonstrate clearly that although the physiological and parkinsonian networks generated by the genetic algorithm exhibit large variability in the 20-dimensional structural space, they cluster together in the 2-dimensional dynamical feature space and exhibit an unambiguous separation between networks classified as physiological and those classified as parkinsonian. Most of the parkinsonian networks show low GS and high SO (Figure 7B, bottom), whereas most of the physiological networks show high GS and low SO (Figure 7A, top left).

However, there is some overlap between the physiological and parkinsonian networks (e.g., bottom left corner, SO ≈ 0.1 and GS ≈ −0.9). This could indicate that these two dynamical features are not enough to disambiguate these parkinsonian and physiological networks. However, it may also turn out that some networks are barely distinguishable even when projected onto a space of many such functional features. Such networks might be particularly useful for giving insight into the transition from physiological state to parkinsonian and vice versa.

It should also be noted, that although the homologies were constrained based on the characteristics of GPe(TA/TI)-STN subnetwork, they show expected dynamical features that depend on the activity of the entire BG network. That is, 70% of parkinsonian networks show a GS < −0.5 and a SO > 0.35, whereas 95% of physiological networks show GS > 0.8 and SO < 0.2. This indicates that the network characteristics of the BG subcircuit (GPe-STN) such as firing rate and phase relationships as shown in Mallet et al. (2008) and Abdi et al. (2015) are predictive of the network state of the entire BG in GS-SO space.

The model used here is a firing rate model, which is a highly simplified formulation that captures solely the evolution of the mean activity of a population. This kind of model inherently limits the insights that can be gained about the dynamics of the system, for example it can give no information about the structure of individual spike trains or higher order correlations within a population. However, it has the advantage of being computationally less expensive and with a restricted parameter space with respect to higher resolution simulation methods (e.g., spiking neuronal networks). It should be noted, that an another modeling work has described a spiking neural network model of the basal ganglia by including the GPe-TA/TI subpopulations (Lindahl and Hellgren Kotaleski, 2016). Although, the model description of the aforementioned work is more detailed (spiking neural networks) as compared to mean field models in our case, the former constrains the system to exactly one solution of effective weights. It would be interesting to check if these effective weights falls under the solution space spanned by the homologies in this work.

In this study we considered those parameters, for which experimental data sets exist, to be fixed. However, in Section 3.3 we discovered that replacing the generated values for the free parameters with the mean of the corresponding distribution does not necessarily preserve the dynamical properties of the networks, at least in part due to correlations between the free parameters. This argument could be equally well applied to the “fixed” parameters. A better approach might be to consider the fixed parameters as free, albeit within the biological range suggested by the experimental data. We plan to incorporate this extension in the future work. This extension will enables us to address changes in effective connectivities in the fixed parameters such as increase in the inhibition from FSI to D2 (Gittis et al., 2011; Corbit et al., 2016) and decrease in striatal lateral inhibition during parkinsonian conditions (Taverna et al., 2008) which are expected to be emergent trends in the weight distributions or can be used as additional post-hoc constraints to reduce the number of homologies.

Ideally one would choose the parameters of the activation function according to F-I curves obtained in in-vitro experiments. However, the activation function S(·) used in the model corresponds to the stationary firing rate response of an entire nucleus, rather than a single neuron. One could consider the parameters of the activation function as free parameters which are adjusted during the optimization process together with the synaptic weights. This would lead to a scenario where not only the synaptic weights but also the activation function parameters are widely distributed across network configurations. Future work can determine whether the additional insight gained by such an approach is sufficient to justify the increased computational complexity.

The cortical input in this model is considered as a feedforward input; the feedback connections from basal ganglia to cortex such as the pallidocortical (Chen et al., 2015) and thalamocortical projections are not included. Hence, the model cannot make predictions with respect to experimental observations such as pharmacological blocking of FSIs cause large negative suppressions of cortical LFPs which are time locked to fast involuntary movements (Klaus and Plenz, 2016). Only excitatory corticostriatal projections are considered, whereas recent data shows that somatostatin interneurons in the cortex also affect the striatal MSNs (Rock et al., 2016). The framework is sufficiently flexible to allow the inclusion of other projections and nuclei, however this would inevitably lead to an increase in number of free parameters.

The oscillations measured in our model in response to a square pulse of cortical input arise due to the change in the effective connectivities between the nuclei. However, oscillations might arise due to other causes, e.g., they may be purely delay driven. Since this model does not consider any delays, the results cannot provide any insight on the role of delays on the oscillations. Secondly, it has also been shown that the beta oscillations arise in striatum by increasing the activity of the striatal cholinergic interneurons pharmacologically (McCarthy et al., 2011). This model does not account for striatal cholinergic interneurons and make an alternative suggestion for the origin of oscillations in the absence of striatal cholinergic interneurons.

In this study we collated experimental data on the strengths of connections within the basal ganglia to determine the fixed parameters of the model (see Section 2.1.1), and showed in Section 3.2 that the relationships between the distributions generated for the free effective connectivity parameters are a good match to experimental observations. On the basis of the connectivity strengths generated by our parameter search (see Table 5), we can also predict structural relationships for physiological and parkinsonian networks which are yet to be verified by experimental results:

Striatopallidal projections are stronger in parkinsonian networks; this is consistent with the hypothesis that the indirect pathway becomes stronger in parkinsonian networks. Parkinsonian networks in our results also show a higher susceptibility to oscillations. These oscillations could initiate in the GPe-STN circuit and could be initiated by a strong indirect pathway, (i.e., a strong inhibition from striatum to GPe), as demonstrated in our recent numerical study (Kumar et al., 2011).

(physiological: J_TA,D2 = −0.4, J_TI,D2 = −0.45; parkinsonian: J_TA,D2 = −2.1, J_TI,D2 = −1.6).

It is also noteworthy that these predictions are not a result of one model specification, but around thousand models that meet the firing rate and phase relationships to be considered physiological or parkinsonian. Therefore, firstly, we can state these predictions more confidently than if they had emerged from a fitting to a unique model specification, and secondly, as no limits on the structural relationships were placed on the parameter search, it seems reasonable to conclude that these emergent differences in these relationships for the physiological and parkinsonian ensemble are relevant for generating the differential dynamics of the two conditions.

We also checked these predictions under three different conditions: (a) increasing the inhibition from D2 to FSI (J_D2,FSI) for parkinsonian networks; (b) keeping the projections from GPe-TA to STN (J_STN,TA) fixed and set to zero; (c) including the projections from D1-MSNs to GPe (J_TA,D1, J_TI,D1) as free parameters. The parameters distributions for first two conditions are largely similar to the original distributions with no or few differences as described below.

In this work, we demonstrated that a particular limited set of features (Table 4) from the GPe-STN sub-circuit can predict dynamical characteristics (GS and SO) relevant for the behavior of the entire BG network under physiological and parkinsonian conditions. It remains to be investigated whether other feature sets have a similar predictive power and what the minimal set of predictive features is.

To test whether the mapping from the feature set in Table 4 to the dynamical characteristics (GS and SO) is one to one (injective), we performed a new set of numerical experiments where networks were constrained by the dynamical features GS and SO rather than by the original firing rates and phase relationships. We then checked whether the firing rate and phase relationships listed in Table 4 emerged as a result. Here, a transient square pulse was used as an input signal in the genetic algorithm. All the networks that showed GS ≤ −0.7 and SO ≥ 0.45 were classified as parkinsonian networks whereas the ones that showed GS ≥ 0.85 and SO ≤ 0.2 were classified as physiological networks. Only a fraction of the resulting networks are consistent with the criteria in Table 4 (Supplementary Figure 7). Hence, we can infer the dynamical feature GS and SO from the firing rates and phase relationships of GPe-TA-TI-STN subnetwork, but not vice versa. For the networks that fulfill the criteria in Table 4, the marginal weight distributions are very similar to the ones shown in Figure 4 (Supplementary Figure 6).

Our parameter search led to over a thousand valid solutions (in 20 dimensional space) for both the physiological and the parkinsonian classification criteria. Although this could simply be a result of lack of constraints on the model, this could also indicate the presence of homologous networks in basal ganglia, in the sense that quite different parameter combinations can give rise to essentially the same dynamics. Homology has not yet been shown to occur in basal ganglia although other neuronal networks have shown occurrences of this phenomenon. However, it seems rather likely due to the complexity of the connectivity (Prinz et al., 2004).

What does the presence of homologous networks in our results imply? One possibility is that it is simply a consequence of the fact that we are mapping a high-dimensional parameter space (here network connections) onto a low dimensional space of network activity features (firing rates and phase relationships). Another possibility is that some parameters have only a small effect on the network activity features that we are interested in studying and can therefore be defined in a “sloppy” manner, i.e., with values are distributed over a large area in the parameter space (Gutenkunst et al., 2007). Our results show that indeed some parameters can be replaced by the mean or shuffled, whilst retaining the original classification for more than 95% of the ensemble. However, other parameters are much more sensitive to such manipulations. Thus, the existence of homologous networks could imply that that rather than absolute values of the parameters, it is their relative values that determine the dynamical state.

This interpretation is supported by the presence of correlations between the generated connectivity values, and the greater sensitivity of the dynamical system to those parameters that exhibit stronger cross-correlations (Figures 4B,C, 5C). The presence of correlations between the free parameters is also suggested as a possible reason for the presence of multiple solutions by Achard and De Schutter (2006). We propose that the approach of generating ensembles of homologous networks and analyzing them in terms of dynamical features is advantageous, as it gives us a framework for investigating complex interacting biological systems such as the basal ganglia which exhibit variability in structural parameters either due to nature (genetic tendencies) or nurture (plasticity in response to environmental factors).

In particular, it will enable us to examine whether the diversity shown by Parkinson's disease as a pathology is rooted in structural variability, to characterize the dynamic and structural properties of the transitions between different dynamical regimes, and uncover possible compensatory mechanisms. For example, we note that whereas the vast majority of physiological network configurations are located in a cluster characterized by effective suppression of GPi activity (high GS) and low susceptibility to oscillations (low SO), and parkinsonian networks typically display low GS and high SO, some network configurations do not fall into these clusters. Specifically, some networks classified as physiological with respect to their rates and phase relationships exhibit low GS, and some parkinsonian networks have low SO. These networks may simply be outliers, but they could also indicate the existence of transitional states between physiological and parkinsonian conditions.

Apart from the two suggested dynamical features (GS and SO), many others could be used to decompose a parkinsonian or physiological dynamical regime into more easily understandable dynamical features. Recent work by Escande et al. (2016) shows that the activity of D1-MSNs in response to optogenetic stimulation in cortex showed a gradient decrease from sham to partially lesioned 6-OHDA (and asymptomatic animals) to heavily lesioned (symptomatic animals). Hence, the gain of D1-MSNs could be yet another dynamical feature that could used for further functional classification.

On the basis of our results, we conclude that it may ultimately be more fruitful to study Parkinson's disease in terms of its dynamical features rather than its structural changes. Not only are dynamical features relatively low-dimensional, they also provide symptomatic targets for therapy (e.g., reduce susceptibility to oscillations) without presupposing specific underlying structures which may hold only in the mean, and not be applicable to a given patient.

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.