Thalamocortical dynamics underlying spontaneous transitions in beta power in Parkinsonism

&NA; Parkinson's disease (PD) is a neurodegenerative condition in which aberrant oscillatory synchronization of neuronal activity at beta frequencies (15–35 Hz) across the cortico‐basal ganglia‐thalamocortical circuit is associated with debilitating motor symptoms, such as bradykinesia and rigidity. Mounting evidence suggests that the magnitude of beta synchrony in the parkinsonian state fluctuates over time, but the mechanisms by which thalamocortical circuitry regulates the dynamic properties of cortical beta in PD are poorly understood. Using the recently developed generic Dynamic Causal Modelling (DCM) framework, we recursively optimized a set of plausible models of the thalamocortical circuit (n = 144) to infer the neural mechanisms that best explain the transitions between low and high beta power states observed in recordings of field potentials made in the motor cortex of anesthetized Parkinsonian rats. Bayesian model comparison suggests that upregulation of cortical rhythmic activity in the beta‐frequency band results from changes in the coupling strength both between and within the thalamus and motor cortex. Specifically, our model indicates that high levels of cortical beta synchrony are mainly achieved by a delayed (extrinsic) input from thalamic relay cells to deep pyramidal cells and a fast (intrinsic) input from middle pyramidal cells to superficial pyramidal cells. From a clinical perspective, our study provides insights into potential therapeutic strategies that could be utilized to modulate the network mechanisms responsible for the enhancement of cortical beta in PD. Specifically, we speculate that cortical stimulation aimed to reduce the enhanced excitatory inputs to either the superficial or deep pyramidal cells could be a potential non‐invasive therapeutic strategy for PD. HighlightsCoupling changes within and between circuit nodes lead to cortical beta enhancement.Input propagation delays play a crucial role in the up‐regulation of cortical beta.Beta power could be modulated by altering lamina specific inputs.


Keywords:
Parkinson's disease Dynamic causal modelling Beta oscillations Thalamocortical interactions Effective connectivity A B S T R A C T Parkinson's disease (PD) is a neurodegenerative condition in which aberrant oscillatory synchronization of neuronal activity at beta frequencies (15-35 Hz) across the cortico-basal ganglia-thalamocortical circuit is associated with debilitating motor symptoms, such as bradykinesia and rigidity. Mounting evidence suggests that the magnitude of beta synchrony in the parkinsonian state fluctuates over time, but the mechanisms by which thalamocortical circuitry regulates the dynamic properties of cortical beta in PD are poorly understood. Using the recently developed generic Dynamic Causal Modelling (DCM) framework, we recursively optimized a set of plausible models of the thalamocortical circuit (n ¼ 144) to infer the neural mechanisms that best explain the transitions between low and high beta power states observed in recordings of field potentials made in the motor cortex of anesthetized Parkinsonian rats. Bayesian model comparison suggests that upregulation of cortical rhythmic activity in the beta-frequency band results from changes in the coupling strength both between and within the thalamus and motor cortex. Specifically, our model indicates that high levels of cortical beta synchrony are mainly achieved by a delayed (extrinsic) input from thalamic relay cells to deep pyramidal cells and a fast (intrinsic) input from middle pyramidal cells to superficial pyramidal cells. From a clinical perspective, our study provides insights into potential therapeutic strategies that could be utilized to modulate the network mechanisms responsible for the enhancement of cortical beta in PD. Specifically, we speculate that cortical stimulation aimed to reduce the enhanced excitatory inputs to either the superficial or deep pyramidal cells could be a potential noninvasive therapeutic strategy for PD.
Increased oscillations in the beta band have been observed both in Parkinson's disease (PD) patients (Brown et al., 2001) and experimental animal models of the disease (Bergman et al., 1994;Sharott et al., 2005). However, it remains unknown why dopamine depletion leads to excessive synchronization across the CBGTC circuit during PD (Jenkinson and Brown, 2011;Hammond et al., 2007, Leblois, 2006. The onset of measurable oscillations in experimental Parkinsonism takes several days post dopaminergic cell loss (Mallet et al., 2008a). One potential explanation for this observation is that the reduction in dopaminergic drive may lead to plastic changes and give rise to abnormal synchronization in neural activity within and between different nodes of the CBGTC circuit. Regardless of the exact mechanism, a positive correlation between excessive beta activity and motor deficits has been reported by several clinical studies (Eusebio et al., 2009;Kuhn et al., 2008). When Parkinsonian motor deficits are attenuated with pharmacological (Levodopa) or neuromodulatory interventions (deep brain stimulation or optogenetics), a reduction in synchronization is observed in the beta-frequency band across different species, including humans (Brown et al., 2001;Eusebio et al., 2009;Kuhn et al., 2008;Levy et al., 2002;Priori et al., 2004;Silberstein et al., 2005), 1-methyl-4-phenyl-1,2, 3,6-tetrahydropyridine treated non-human primate models of PD (Heimer et al., 2006;Nambu and Tachibana, 2014) and a 6-hydroxydopamine (6-OHDA)-lesioned rat model of PD (Gradinaru et al., 2009;Sharott et al., 2005).
Although excessive synchrony in the beta band (i.e. beta power) is traditionally described as a sustained event when averaged over seconds (Brittain and Brown, 2014;Brown, 2007;Lopez-Azcarate et al., 2010), it primarily manifests as intermittent events of high beta power or "beta bursts" (Feingold et al., 2015;Sherman et al., 2016;Tinkhauser et al., 2017;Little et al., 2012;Leventhal et al., 2012). Beta bursts have been defined operationally as epochs of beta oscillations that surpass a certain thresholdand their presence has been quantified in physiological (Sherman et al., 2016;Feingold et al., 2015) and pathological neural activity (Tinkhauser et al., 2017;Little et al., 2012). More specifically, in Parkinson's disease, the probability of long beta bursts has been positively correlated with PD motor symptom severity (Tinkhauser et al., 2017;Little et al., 2012).
Adaptive Deep Brain Stimulation (aDBS) is an intervention that has been developed to account for the transient nature of pathological neural synchrony in the beta band. In contrast to conventional DBS (cDBS), which continuously delivers high-frequency stimulation, aDBS adapts stimulation delivery according to the level of beta power (Little et al., 2013(Little et al., , 2012Rosa et al., 2015), showing greater clinical efficiency (higher motor symptom relief and fewer secondary effects) than cDBS and random stimulation (Little et al., 2013).
From a neuronal network perspective, several studies have proposed that altered basal-ganglia output leads to excessive beta synchrony and motor impairments in PD (Bevan et al., 2002;Holgado et al., 2010;McCarthy et al., 2011;Terman et al., 2002). Employing Dynamic Causal Modelling (DCM) (Friston et al., 2003), a framework for specifying, fitting and comparing mathematical models of neural circuitry, Moran et al. (2011) and Marreiros et al. (2013) indicated modulation of the hyperdirect pathway and the projection from the subthalamic nucleus and globus pallidus externus as potential mechanisms for beta power enhancement in dopamine-depleted states.
Some experimental studies, on the other hand, support the role of cerebral cortex in the generation and modulation of beta oscillations (Jensen et al., 2005;Yamawaki et al., 2008). This perspective has motivated the consideration of cortical interlaminar interactions in the regulation of beta power. In the healthy state, the generation and modulation of beta oscillations has been investigated using DCM (Bhatt et al., 2016), revealing a link between a set of laminar specific interactions within the primary motor cortex and the enhancement/suppression of beta power evoked by movement. Using a theoretical model, Sherman et al. (2016) suggested that high beta power events in the physiological state emerge through cortical laminar interactions conditioned by temporal characteristics of the distal and proximal synaptic drives in the neocortex.
Motivated by these studies, we hypothesized thatin the Parkinsonian statean alteration of interlaminar and laminar-specific connectivity in the Thalamocortical (TC) loop contributes to the mechanisms generating the parkinsonian spectral profile. We focused on the TC loop due to the anatomical and functional characteristics of this network: 1) the cortex is an optimal target for non-invasive therapeutic techniques such as TMS and TACS (Barker et al., 1985;Cantello, 2002;Kobayashi and Pascual-Leone, 2003;Herrmann et al., 2013); 2) the thalamus is the only CBGTC node projecting directly to cortex, allowing for the integration of information from subcortical structures to the motor cortex (Wise and Donoghue, 1986;Brazhnik et al., 2016) and 3) cortex and thalamus establish a reciprocal relationship (Hooks et al., 2013), which is thought to play a key role in physiological and pathological sensory and motor computations (Sherman & Guillery, 2001). In PD, where motor impairments are the cardinal symptoms, understanding the synaptic dynamics and organization of the thalamocortical (TC) circuit could potentially shed light on pathophysiological mechanisms. Accordingly, we used cross spectral density (CSD) -DCM (Moran et al., 2009(Moran et al., , 2011 with a neural mass model of the TC loop (van Wijk et al., 2018) to characterize its contribution to spontaneous beta power fluctuations observed in the motor cortex of 6-OHDA-lesioned Parkinsonian rats.

Electrophysiological recordings in parkinsonian rats
The spectral data used in this study was based on motor cortex field potentials (electrocorticograms) recorded in 36 urethane-anesthetized rats rendered Parkinsonian by unilateral 6-OHDA lesions of midbrain dopaminergic neurons. To record electrocorticogram (ECoG) data, a steel screw electrode was implanted over the right somatosensory-motor cortex ipsilateral to the 6-OHDA lesion and referenced to a steel screw electrode implanted over the ipsilateral cerebellar hemisphere. Electrophysiological recordings were carried out 21-42 days after surgery for the induction of 6-OHDA lesions, thus allowing for changes in the CBGTC circuit to stabilize. For detailed descriptions of electrode implantation, anaesthesia, surgical induction of 6-OHDA lesions and related procedures, please refer to (Mallet et al., 2008a(Mallet et al., , 2008bSharott et al., 2017). Only ECoG recordings made during periods of spontaneous 'cortical activation' were considered in this study (Mallet et al., 2008a(Mallet et al., , 2008bSharott et al., 2017). All experimental procedures were carried out on adult male Sprague-Dawley rats (Charles River, Margate, UK) and were conducted in accordance with the Animals (Scientific Procedures) Act, 1986 (UK).

Data processing
All operations described in this section were performed in Matlab 2017a/2018a. Data and code that support the findings of this study are available from the corresponding author (HC; hayriye.cagnan@ndcn .ox.ac.uk) upon request. Recordings were down-sampled to 1000 Hz from 16,000 Hz. To characterize the spontaneous beta power fluctuations typically observed in PD, we defined two conditions based on instantaneous beta powercondition one being Low Beta (LB) power and condition two being High Beta (HB) power. These conditions were based on fluctuations in beta power that enabled us to select data-features (i.e., timeseries) for subsequent dynamic causal modelling that were representative of the two conditions.
To extract the beta power envelope, we applied a second order bandpass Butterworth filter with cut-off frequencies at 15-35 Hz to the ECoG recording and subsequently employed the Hilbert transform to compute the envelope of the ECoG in the beta frequency band. Each envelope was then divided into non-overlapping epochs of 500 msec The two conditions were subsequently derived from a relative threshold applied to the area under the envelope across the 500msec epochs: (1) LB epochs consisted of segments whose envelope area fell below the 5th percentile of the envelope area observed across all epochs, and (2) HB epochs consisted of segments whose envelope area was above the 95th percentile of the envelope area observed across all epochs (Fig. 1). 5th and 95th percentile thresholds were determined per dataset. From each recording, we randomly selected 5 epochs per condition (n ¼ 5). This number corresponds to the minimum number of epochs found in either of the two conditions across all recordings. A detailed comparison between the above threshold and more conventional ones can be found in the Supplementary material (Fig. S1).

Dynamic causal modelling (DCM)
DCM for cross spectral density is used to infer the hidden (neuronal) states (z) and synaptic parameters (θ) that generate spectral features of observed data (u) (Moran et al., 2009(Moran et al., , 2011. Hidden states and unknown parameters cannot be observed directly but can be estimated under a generative or forward model. This model comprises a biophysical neural mass model and the spectral composition of neural and channel noise (Moran et al., 2008). The neural mass model is expressed in terms of a differential equation with the following form: The neural mass model (f )together with a likelihood model mapping hidden states to observed measurementsconstitutes a generative model; namely, a probabilistic mapping between neural fluctuations and the spectral content of observed activity. Using a Bayesian framework, DCM estimates the (posterior) probability density over the synaptic parameters, which are the most likely value of the hidden parameters, given the observed data (Moran et al., 2011). The generative (neural mass) model calls on its biophysical parameters to describe the evolution of voltages (v) and currents ðiÞ in each subpopulation of neurons (Jansen and Rit, 1995). In addition to estimating the posterior density over model parameters (e.g., synaptic connection strengths and the amplitude of neuronal fluctuations), DCM also provides an estimate of the evidence for a particular model or network architecture implicit in the generative model. This allows one to compare different models or hypotheses using Bayesian model comparison. A complete description of the mathematical framework that underwrites DCM can be found in .

Neural mass model of the thalamocortical circuit
A neural mass model of the Thalamocortical circuit was created comprising two formally distinct neural mass models of the motor cortex and the thalamus using the new generic framework for Dynamic Causal Modelling (van Wijk et al., 2018) (Fig. 2). Here, we adopted the motor cortex microcircuit (MMC) model developed by Bhatt et al. (2016) and coupled it to a model of the thalamus, based on thalamic anatomical literature (Shepherd and Grillner, 2010;Douglas and Martin, 2004). As Fig. 1. Extraction of low beta and high beta power features isolated from ECoG data. Panel A. shows the segmentation of the envelope into 500 msec epochs (5 s as an example). Panel B. depicts the area under the beta band envelope for each epoch. If an epoch had an area under the curve below the 5th percentile of the area under the envelope observed across all epochs (blue line), it was classified as low beta (*1); if an epoch had an area under the curve above the 95th percentile of the area observed under the envelope across all epochs (red line), it was classified as high beta (*2). 5th and 95th percentile thresholds were determined per dataset. Epochs in between the two percentiles were not considered. Panel C. shows the corresponding low beta (dark blue) and high beta (red) epochs in the ECoG signal filtered at 15-35 Hz. superficial pyramidal cells (SP) in the supragranular layer, middle pyramidal cells (MP) in the granular layer, deep pyramidal cells (DP) in the infragranular layer and, inhibitory interneurons (II) as a common inhibitory subpopulation to the 3 cortical laminae. Intrinsic synaptic connections among the above subpopulations comprise a reciprocal connection between superficial and middle pyramidal cells, a reciprocal connection between superficial and deep pyramidal cells, a reciprocal connection between each of the three pyramidal subpopulations and the inhibitory subpopulation and finally, a self-inhibitory connection to each cortical node. The bottom part of the diagram depicts the thalamus and its subpopulations: reticular thalamic cells (RET) as the inhibitory subpopulation of the thalamus and relay cells (REL) as the excitatory subpopulation of the motor thalamus. Intrinsic synaptic connectivity of the thalamus comprises a reciprocal connection between relay and reticular cells and self-inhibitory connection of reticular cells. As corticothalamic extrinsic connections, deep pyramidal cells were considered to send afferents to both relay and reticular subpopulations, while the model space for thalamocortical projections is described in section 2.3.4 and illustrated in Fig. 3 (top panel).
with previous models of the sensory cortexand incorporating the work of Yamawaki et al., 2014-Bhatt et al. (2016 used 3 excitatory subpopulations (neuronal ensembles consisting of "superficial", "middle" and "deep" pyramidal cells located in the supragranular, granular and infragranular cortical layers, respectively) and one common inhibitory subpopulation (inhibitory interneurons) to model the primary motor cortex. In the MMC model, the coupling between these subpopulations (GABAergic or glutamatergic synapses) is tailored according to synaptic characteristics of the primary motor cortex: a reciprocal connection between superficial and middle pyramidal cells (Yamawaki et al., 2014), a reciprocal connection between superficial and deep pyramidal cells (Hooks et al., 2013;Yamawaki and Shepherd, 2015;Anderson et al., 2010;Weiler et al., 2008), a reciprocal connection between each of the three pyramidal subpopulations and the common inhibitory subpopulation (Fino et al., 2013), and a cell type specific self-inhibitory connection Yoshimura and Callaway, 2005). The self-inhibitory connections aim to capture laminar-specific inhibition, mediated by local inhibitory neurons (K€ atzel et al., 2011). For further discussion on recurrent inhibitory connections in the context of the canonical microcircuit model, please refer to Auksztulewicz and Friston (2015).
In this study, the thalamus was modelled using an excitatory subpopulation (neuronal group of thalamic relay cells) and an inhibitory subpopulation (neuronal group of thalamic reticular cells) (Shepherd and Grillner, 2010) that were connected as follows: a reciprocal connection between relay and reticular cells (Harris, 1987;Cox et al., 1997) and a self-inhibitory connection of reticular cells (Shu and McCormick, 2002). Although we acknowledge that there are distinct thalamic nuclei (i.e. neuronal ensembles receiving afferents from different brain regions (Sherman & Guillery, 2001) the thalamus was modelled here as a single neuronal mass model. An important extension of the current work would be to subdivide the motor thalamus (ventral anterior, ventral lateral and ventral medial nuclei in rodents) into input zones that receive GABAergic drive from the Basal-Ganglia and glutamatergic drive from the cerebellum (Kuramoto et al., 2011;Nakamura et al., 2014).
To model the extrinsic synaptic interactions between the motor cortex and thalamus we used two corticothalamic projections from deep pyramidal cells to thalamic relay cells and thalamic reticular cells (Bourassa et al., 1995;Jones, 2001). Although the thalamus is thought to project to all layers of the cortex (Hooks et al., 2013) and the ventromedial nucleus (VM) of the motor thalamus has been shown via immunochemistry studies to project mainly to layers I and II of the motor and anterior cingulate cortices (Arbuthnott et al., 1990;Clasc a et al., 2012;Kuramoto et al., 2015), it is not clear which thalamocortical projections are important in modulating beta power. To resolve this, we considered different models to test the impact of including different connections on model evidence (section 2.3.4.). Operationally, the difference between within-structure (intrinsic) connections and between-structures connections (extrinsic), relies on their propagation delay parameters (Table.1). Aiming to bring bio-plausibility to the TC neural mass model architecture, here the propagation delays of extrinsic connections (thalamocortical and corticothalamic projections) were set to 8msec while the propagation delays of intrinsic dynamics (intracortical and intrathalamic connections) were set to 1 msec

Neural dynamics
In DCM, neural dynamics (i.e., fluctuations in voltages and currents) at the subpopulation level is described by two key operations (Eq. (2)): a convolution operator and an output operator (Moran et al., 2007). The convolution operator transforms presynaptic inputs (firing rate) into postsynaptic membrane potentials based on a synaptic impulse response function, which considers the nature of the synapse (i.e. excitatory or inhibitory).
The output operator consists of a non-linear function that converts the postsynaptic membrane potentials into a firing rate to be relayed to another subpopulation. This is conveyed through a sigmoid function S which captures the membrane sensitivity and firing threshold of each subpopulation. Furthermore, the shape of the sigmoid function (slope) measures the efficacy of a presynaptic ensemble to generate output. This output is additionally scaled by the synaptic coupling strength as illustrated by the following generic second order differential equation: Here, the averaged membrane potential v of the subpopulation j in the source a is influenced by subpopulations of the same source with synaptic strength γ and subpopulations from different sources with synaptic strength λ. Intrinsic synapses γ show a positive synaptic strength if glutamatergic and negative synaptic strength if GABAergic. S denotes the sigmoid function above and T the subpopulation-specific membrane time constant. Endogenous fluctuations or input, I is modelled as a mixture of white and pink noise and drives middle pyramidal cells and thalamic relay cells. The rationale for modelling afferent input to both the motor cortex and thalamus (as opposed to restricting the model to thalamic input) rests on the fact that both motor cortex and thalamus receives input from other (unmodelled) components of the motor system. Examples here include inputs from supplementary motor areas and premotor cortex to primary motor cortex (Jones et al., 1975) and inputs from basal ganglia and cerebellum to motor thalamus (Kuramoto et al., 2011;Nakamura et al., 2014).
In this study, we used DCM for cross spectral density (Moran et al., 2009(Moran et al., , 2011 where the data generated by a model of neural hidden states Table 1 Prior expectations set for the parameters of the baseline condition (Low beta). CT-corticothalamic projections; TC-thalamocortical projections; MMCmotor microcircuit; THAL-thalamus; SP-superficial pyramidal cells; MP-middle pyramidal cells; DPdeep pyramidal cells.  (7)). A summary of the parameters described in this section and their prior values are shown in Table.1. Prior values were based on previous DCM studies (Bhatt et al. 2016) and optimized for our study.

Model inversion
In DCM, model inversion iteratively tunes the model's parameters to optimize the fit of the predicted electrophysiological data to the observed data. Using a standard (variational) Bayesian scheme, model inversion uses priors to constrain the search of parameter space to explain the observed spectral features of electrophysiological data. When fitting the data (i.e., inverting the model), the optimization of model parameters uses a variational Laplace scheme to minimize a (free energy) bound on (negative) log model evidence. This free energy approximation to model evidence is subsequently used for model comparison Friston and Stephan, 2007).
In brief, model evidence is the (marginal) likelihood of observing data given a model, pðyjmÞ. It reflects a balance between accuracy (goodness of fit between predicted and observed spectral densities) and complexity Fig. 3. Competing models of the Thalamocortical circuit as described by factors 1 and 2 (9 architectures times 16 modulatory configurations). The diagram on the top (factor 1: architecture) describes the 9 families of models constructed to elucidate which thalamocortical projections are the most plausible explanation for the generation of beta oscillations (1.-3.) accounts for a singular projection from thalamus to motor cortex via superficial pyramidal cells, middle pyramidal cells and deep pyramidal cells; (4.-6.) accounts for two afferents to two excitatory subpopulations of the motor cortex via superficial and middle pyramidal cells, middle and deep pyramidal cells and superficial plus deep pyramidal cells, and (7.-9.) accounts for projections to the superficial pyramidal subpopulation and inhibitory interneurons, the middle pyramidal subpopulation and inhibitory interneurons and deep pyramidal cells and inhibitory interneurons. To disclose the synaptic modulation (intrinsic and/or extrinsic) responsible for an enhancement of beta power, the models on the bottom (factor 2: modulatory configuration) feature 16 different modulatory configurations, under each of the 9 architectures described above. The first eight set of connections (1.-8.) entail extrinsic and intrinsic synaptic modulation (except for model 8, with no intrinsic modulation) and the second eight set of connections (9.-16.) considers intrinsic modulation only. The intrinsic modulatory connections in family 2 were: (1. and 9.) cortical modulation via reciprocal connection between superficial and deep pyramidal subpopulations; (2. and 10.) cortical modulation via reciprocal connection between superficial and middle pyramidal subpopulations; (3. and 11.) cortical modulation via self-inhibitory connection of the inhibitory interneurons subpopulation; (4. and 12.) thalamic modulation via reciprocal connection between reticular cells and relay cells; (5. and 13.) cortical and thalamic modulation via reciprocal connection between superficial and deep pyramidal subpopulations plus reciprocal connection between reticular cells and relay cells; (6. and 14.) cortical modulation via reciprocal connection between superficial and middle pyramidal subpopulations plus reciprocal connection between reticular cells and relay cells; (7. and 8.) self-inhibitory connection of the inhibitory interneuron subpopulation plus reciprocal connection between reticular cells and relay cells and lastly, (8. and 16.) the null hypothesis that neither extrinsic nor intrinsic connections change to explain condition specific changes in cortical beta power (i.e. enhancement of beta).
(divergence between prior and posterior parameter estimates) (Stephan et al., 2010). This balance depends upon the expected precision of the observed data. Given the high signal to noise ratio in the data obtained using the electrocorticographic recording method, the expected precision of observed data was assumed to be high (with a log precision of 12).
At this stage, low beta power was set as our baseline condition (with prior expectations optimized to best explain its spectral features). Condition-specific effects (B parameters) on both extrinsic (between-regions) and intrinsic (within-regions) coupling strengths (Moran et al., 2007) were used to explain periods of high beta power. In other words, we estimated the changes in synaptic efficacy required to move from a low beta power condition to a high beta power condition.

Bayesian Model Comparison and parameters analysis
A set of models were implemented which varied according to 2 factors: i) the laminar-specificity of thalamocortical projections that generate beta oscillations, and ii) the changes in synaptic connectivity within the TC loop (intrinsic and/or extrinsic) required to induce a transition from a low beta power condition to a high beta power condition.
The first factor comprised 9 families (types) of models. These models had identical intrinsic and corticothalamic connections as described in section 2.3.1 and illustrated in Fig. 2 (Fig. 3, top panel).

but differed in the laminar targets of thalamocortical afferents: 1) superficial pyramidal cells; 2) middle pyramidal cells; 3) deep pyramidal cells; 4) superficial plus middle pyramidal cells; 5) middle plus deep pyramidal cells; 6) superficial plus deep pyramidal cells; 7) superficial pyramidal cells plus inhibitory interneurons; 8) middle pyramidal cells plus inhibitory interneurons and 9) deep pyramidal cells plus inhibitory interneurons
The second factor comprised 16 families of models that varied in the set of connections that could show condition specific effects. For each one of the 9 architectures in the first factor, we explored condition specific effects by including or not the following features: intracortical modulatory synapses; intrathalamic modulatory synapses and extrinsic (between cortex and thalamus) modulatory synapses (Fig. 3, bottom panel). There were therefore 9 Â 16 ¼ 144 candidate models in total.
Bayesian Model Comparison (BMC) was used to determine the model with the highest log-model evidence among the models described above (Stephan et al., 2010). We then characterised the parameters of the winning model at the group level using Parametric Empirical Bayes (PEB) . In brief, PEB uses a hierarchical model, where the parameters from each subject's DCM are summarized by a posterior density (the expected connectivity strength and posterior covariance). Random effects at the between subject level are similarly inferred to inform the group-level parameter estimates. This means that PEB allows for estimation of fixed and random effects in an optimal fashion; implicitly reducing the influence of "outlier subjects" on the posterior density at the group level. Note that in Bayesian inference, (unstandardized) effect sizes are available explicitly in terms of posterior expectations and Bayesian credible intervals (Standardized effect sizes such as correlations are replaced by differences in model evidenceimplicit in Bayesian model comparison). Here, only a subset of parameters was taken to the group level and assumed to exhibit random effects: synaptic coupling strength of intrinsic and extrinsic connections (G and A parameters in the DCM respectively); condition-specific effects on coupling strength (B parameters) and the time constants of subpopulations (T).

Model selection
The model with the highest evidence for the transition from low beta epochs to high beta epochs (Fig. 4) was that with i) an architecture featuring thalamocortical projections from the REL-DP and REL-II in cortex and ii) modulatory changes in: intrinsic connections at the cortical level between SP-MP; intrinsic connections at the thalamic level, between REL-RET; corticothalamic extrinsic connections from DP-REL and DP-RET and thalamocortical connections from REL-DP and REL-II. The set of differential equations explaining the neural dynamics of the winning model can be found in the supplementary materials (Fig.S5). The winning model shows a free energy difference (i.e., log Bayes factor) of approximately 6 from the next closest model (Fig.S2). This corresponds to very high evidence for the winning model, in relation to alternative explanations.
Using fixed-effects Bayesian Model Comparison (FFX-BMC) to make inferences at the family level, the architecture with thalamic projections to DP and II showed the highest evidence across subjects (Fig. 5A). Similarly, the condition specific effects in the reciprocal connection between SP -MP, REL-RET and DP-REL plus connections from DP-RET and REL-II had the highest posterior probability (Fig. 5B). These results confirm our hypothesis that both the laminar-specificity of extrinsic connectivity and intrinsic connections are key elements underlying the modulation of oscillatory activity in the beta band.

Parameter analysis
The results from our second level analysis (PEB modelling of A,G,B and T parameters at the group level) suggest that the transition from low beta state to high beta state is induced by i) an increase in synaptic strength in connections from relay cells to inhibitory interneurons, relay cells to deep pyramidal cells, middle pyramidal cells to superficial pyramidal cells and relay to reticular cells; plus ii) a reduction of synaptic strength in connections from superficial to middle pyramidal cells, deep pyramidal to both relay and reticular cells and from reticular to relay Fig. 4. Observed and expected power spectral densities (PSD). (A) Spectral features to be explained by a DCM: Red lines depict the mean high beta spectral densities and blue lines the mean low beta spectral densities from each of the 36 rats. (B) Group mean of HB spectral densities in red and LB spectral densities in blue. Respective variabilities (75th and 25th percentiles of the mean spectra) denoted in light red and light blue. (C) Goodness of the fits between mean data spectral densities and spectral densities generated by the winning model. The full red line shows the mean of high beta data and the dark red dashed line the high beta spectra estimated by the winning model (correlation coefficient, r ¼ 0.9997). The full dark blue line refers to the mean of low beta data and the dark blue dashed line to the low beta spectra produced by the winning model (correlation coefficient, r ¼ 0.9957). cells (Fig. 6). Additionally, from the posterior distribution of our B parameters we assessed the effect size of each modulatory connection on the enhancement of betaillustrated in the bar plot below (Fig. 6).
Considering the connections that showed the greatest change to explain beta enhancement: REL-DP, MP-SP and REL-RET; we further analysed, via forward modelling, the impact of simultaneous alteration of the above connection strengths on the magnitude of beta power. As such, we aimed to characterize the contribution of these three parameters to  (Fig. 3, Factor 1: Architecture, number 9). These results suggest that thalamocortical projection to the deep pyramidal cells and cortical inhibitory subpopulation (in thick lines) were crucial for the generation of beta oscillations and that this effect was consistently observed across subjects (posterior probability of 1). Diagram and bar plot (B) indicate the modulatory connections of our winning model (Fig. 3, Factor 2: Modulatory configuration, number 6). The diagram shows the set of connections as thick lines to have a higher likelihood (compared to the homologous 15) of inducing the power spectral changes observed (beta enhancement). These being: a reciprocal connection between superficial and middle pyramidal subpopulations, reciprocal connection between thalamic relay and reticular subpopulation and a reciprocal extrinsic connection between deep pyramidal cells and thalamic relay cells, an extrinsic connection from deep pyramidal cells to thalamic reciprocal cells and from thalamic relay cells to cortical inhibitory interneurons. The bar plot shows a posterior probability greater than 0.99 for the modulatory configuration described above and a negligible posterior probability of approximately 0.004 for a modulatory configuration which assumed the same modulatory characteristics as the winning model except for the intrinsic synaptic mechanisms of the motor cortex; i.e., presenting an intracortical modulation via reciprocal connections between superficial and deep pyramidal cells instead of reciprocal connections between superficial and middle pyramidal cells.   6. Average modulatory effect of B parameters (condition-specific parameters) obtained via Parametric empirical Bayes analysis . The two bar plots on the left-hand side illustrate the absolute connection strength of each modulatory connection in the low and high beta conditions. The bar plot in the centre shows how connectivity strength of B parameter changed at the group level in order to induce an increase of beta power. Negative values of change indicate a reduction in connectivity strength and positive values an increase. The anatomy of these connections is illustrated in the diagram on the right-hand side.
the gradual transition between the two stateslow and high beta power.
To do so, we have effectively replaced the MP-SP and REL-RET posteriors with values between À1 and 1 and plotted the summed beta spectral output (15-35 Hz) normalised by the summed beta spectral output  in the baseline condition (i.e. low beta power). This post hoc simulation allows us to characterize the selective effect of specific connections on the expression of cortical beta power level. Fig. 7A, suggests that when the intrinsic thalamic connection from relay to reticular cells have low levels of coupling strength, high beta power appears abruptly as the coupling strength from middle to superficial pyramidal cells is increased. This emergence of beta is modulated by the gradual change in coupling strength from thalamic relay cells to deep pyramidal cells, where reduced coupling leads to higher levels of beta. On the other hand, Fig. 7B, suggests that when the connectivity from relay to reticular cells is high, a concurrent increase in the coupling from both relay to deep pyramidal cells and middle to superficial pyramidal cells is required to achieve beta power enhancement. In short, the level of cortical beta depends on the magnitude of excitatory inputs to both superficial and deep pyramidal.
We additionally explored the spectral output of each subpopulation of our TC neural mass model (Fig. 8) and observed that all nodes generated spectral curves within the beta band during both conditions and increased their power from condition one (LB) to two (HB) as expected. Please note that the plotted curves represent the source-space activity, which has been estimated separately for each data set.

Discussion
In this study, we aimed to identify the network mechanisms that contribute to the dynamic regulation of beta synchrony in the parkinsonian motor cortex. In-vivo studies of the basal ganglia thalamocortical (BGTC) circuit suggest that alterations in the firing rate across the direct and indirect pathways are responsible for the motor impairments observed in PD (Albin et al., 1989;DeLong, 1990;Nambu, 2004). Similarly, in-silico simulations of the BGTC circuit propose that an altered coupling from the subthalamic nucleus to globus pallidus externus, and strengthening of the hyperdirect pathway play an important role in the enhancement of beta synchrony following chronic dopamine depletion (Marreiros et al., 2013;Moran et al., 2011).
Our study complements the literature on PD, while exploring two novel concepts: i) laminar-specific dynamics within the motor circuit as a putative mechanism for the spontaneous modulation of beta power and ii) short-term synaptic processes, i.e. transient alterations in effective connectivity to be responsible for the spontaneous and intermittent nature of beta power observed in Parkinsonian time-series.
Focusing on the Thalamocortical loop of the BGTC circuit, we have employed DCM to identify a model of TC interactions that offers plausible substrates for the transient enhancement of cortical beta. Our study suggests two core features of the thalamocortical circuit that may underwrite the genesis of beta oscillations in the parkinsonian state: 1) laminar specific thalamocortical projections; and 2) modulation of synaptic strength across all network levels (i.e. within and between structures).
To model different levels of beta synchrony, we extracted low beta epochs and high beta power epochs from motor cortex ECoG recordings acquired from urethane-anesthetized rodents rendered Parkinsonian by 6-OHDA lesions. This rodent model is useful in capturing the chronic dopamine depletion that is common to 'late stage' PD, and has been widely used for studies of the mechanisms by which excessive beta synchrony arises and propagates within the BGTC circuit in Parkinsonism. Although urethane has been shown to alter the function of multiple neurotransmitter receptors -which will dictate various aspects of neural activity such as firing rateits effects occur to a much smaller extent when compared to other anaesthetics (Hara and Harris, 2002). The dose of urethane allows clear differentiation of slow wave and activated states; only the latter of which resembles the awake brain state and contains abnormally sustained beta oscillations in the dopamine depleted animal (Mallet et al., 2008a, b;Sharott et al., 2017). Finally, the abnormal beta oscillations present in the BGTC circuit in anesthetized and behaving 6-OHDA lesioned are similar in many respects to those present in unmedicated people with PD (Avila et al., 2010;Brazhnik et al., 2016;Degos et al., 2009;Mallet et al., 2008aMallet et al., , 2008bNevado--Holgado et al., 2014;Sharott et al., 2005Sharott et al., , 2017. We have used DCM in this study since it allows for the quantification of effective connectivity changes underlying transient modulations in cortical beta power. Bhatt et al., 2016 have previously employed DCM to link interlaminar dynamics within the motor cortex to the modulation of beta activity, evoked by movement. Fitting MEG data from healthy subjects to a neural mass model of the motor cortex, Bhat and colleagues reported that the increase in beta power observed due to the transition from grip to rest was induced by an increase in the extrinsic input applied to deep and superficial layers of the cortex. Our study suggests that beta power enhancement in Parkinsonism can be attributed to an increase in excitatory inputs to SP and DP; specifically from MP and thalamic relay cells, respectivelyand a concomitant reduction of excitatory input to MP. In addition, our results highlight the importance of intrinsic interactions in the thalamus for beta power modulation as the excitatory projection from the thalamic relay cells to reticular cells also contributes shows the impact of changing the coupling strength of MP-SP and REL-DP on the beta spectral output, when the coupling between REL-RET is weak. Here, although an increase in synaptic strength from MP-SP is enough to generate relatively high levels of beta spectral output, the connection from REL-DP seems to have a modulatory effect, i.e., the weaker the extrinsic coupling between relay and deep pyramidal cells the higher the beta spectral output. Plot B. considers a (constantly) strong coupling between REL-RET with the same changes in coupling strength between MP-SP and REL-DP. This time, we observe that an increase in beta is achieved with a concurrent strengthening of both MP-SP and REL-DP connections. In both plots, axis x and y denote a reduced connectivity strength when values are between À1 and 0 and an increased connectivity strength when values are between 0 and 1 for connections from middle pyramidal cells to superficial pyramidal cells and from relay cells to deep pyramidal cells respectively. Axis Z and colormap depict the magnitude of beta power.
Similarly, focusing on cortical intrinsic dynamics, Sherman et al. (2016) used a computational model to generate transient high beta power events (i.e. beta bursts), which were temporally identical to those observed in the somatosensory and frontal cortices in the physiological state. Two circuit features have been proposed as crucial for the generation of beta bursts: 1) a drive from the lemniscal thalamus to the proximal dendrites of the pyramidal neurons and inhibitory interneurons in L2/3 and L5 (via the granular layer) and 2) a strong drive from the nonlemniscal thalamus to the distal dendrites of the pyramidal neurons and inhibitory interneurons found in supragranular and infragranular layers.
It should be noted that due to the nature of the neural mass models employed in this study, we were not able to model detailed dendritic dynamics that contribute to the generation and modulation of neural activity in the beta band. Instead, here we have assumed fixed conduction delays for all within-region connections (1 m s) and between-regions projections (8 m s) and did not account for variable propagation delays for inputs arriving to distal and proximal dendrites. This creates a distinction between the excitatory input received by the superficial pyramidal cells versus those received by the deep pyramidal cells and the common inhibitory population; since the latter are attributed to extrinsic projections from thalamus and hence are inherently modelled with longer conduction delays. Nonetheless, our results relate to the observations made in Sherman et al. (2016), given that comparable circuitry mechanisms yielded similar oscillatory effects. In other words, both studies propose that laminar specific glutamatergic inputs to the motor cortex must occur at two temporally separate instances in order to achieve high beta power oscillatory activity.
Furthermore, it is worth noting that alternative models of the TC circuit showing thalamic projections to both the superficial and deep layers of the motor cortex (i.e., all models with architecture number 6, shown in Fig. 3), had lower model evidences than the winning model. Possibly because, unlike the winning model, simultaneous projections from thalamic relay cells to superficial and deep layers do not allow for a differentiation in input delays. Our detailed analysis of the parameter space on gradual beta increase, opposed to a transition from extremely low to extremely high beta power, further corroborated that the intrinsic (shorter delay) input from middle to superficial cells should have high levels of synaptic strengthtogether with the extrinsic (longer delay) input from the thalamic relay to deep pyramidal cellsto explain an upregulation of beta synchrony when the coupling from thalamic relay to reticular is strong. (Fig. 7B). Nevertheless, taking both scenarios into accountstrong vs weak connectivity from relay cells to reticular cellsprojections to superficial pyramidal cells via middle pyramidal cells must assume a moderate coupling strength to avoid a ramping up of beta synchrony at the cortical level; highlighting a potential substrate that could be targeted in order to control and modulate cortical beta.
An important difference between Sherman et al. (2016) and our study stems from the assumptions made on thalamic activity patterns. Sherman et al. (2016) posit that thalamic activity should be in the alpha band to drive beta bursts in the somatosensory and frontal cortices. However, in our study, thalamic neurons exhibited activity in the beta band during both low and high beta power conditions (Fig. 8). Our results are supported by recent experimental work showing a substantial and coherent enhancement of beta activity (30-36 Hz) in the motor thalamus and motor cortex of behaving 6-OHDA-lesioned rats (Brazhnik et al., 2016). There is also evidence of aberrant beta synchrony in the thalamus of unmedicated PD patients (Kempf et al., 2009). Taken together, these results emphasise that thalamic neural activity in the beta band is likely to be a contributing circuit feature for the generation of aberrant beta synchronization in PD, and highlight a functional coupling between the thalamus and deep layers of the motor cortex.

Limitations and future directions
We would like to highlight that in DCM, "winning model" is a relative terminology. Winning model is effectively the most plausible model to generate the observed data among the set of models tested. As indeed there are plenty of configurations a neural mass model can assume (here 144 were tested), DCM is only a robust method if used as a hypothesisdriven approach (Stephan et al., 2010).
Finally, an interesting extension of this study would involve using a conductance based neural mass model of the thalamocortical circuit to investigate in detail the intrinsic neurotransmitter dynamics underlying beta enhancement in the parkinsonian thalamocortical circuit. While convolution-based models are useful in explaining the "macroscopic mechanisms" underlying a given data set, such as network architecture and effective connectivity; conductance-based models would capture the dynamics of intrinsic ion channel mediators, such as Glutamaergic and GABAergic receptors (refer to Moran et al., 2013). This could be extremely interesting given the modulatory role that GABA-A is thought to exert over beta oscillations (Brazhnik et al., 2016;Hall et al., 2010;Jensen et al., 2005;Yamawaki et al., 2008).

Conclusion
A broadly accepted postulate concerning healthy and Parkinsonian states of the CBGTC circuit is that they exhibit differential patterns of synchronization at beta frequencies (Brittain and Brown, 2014;Gatev et al., 2006;Hammond et al., 2007). While exaggerated beta synchronization has been associated with more frequent high power beta bursts in PD (Tinkhauser et al., 2017;Little et al., 2013), healthy states seem to manifest as an adequate proportion of high beta events and therefore a Fig. 8. Spectral output of subpopulations at Low Beta (LB) and High Beta (HB). All neural groups (in cortex and thalamus) have generated spectral responses within the beta band in both conditions as expected. Together with an increase in power a small increase in frequency peak of approximately 1-2 Hz is also apparentwhen comparing the output generated in condition 1 (LB) with condition 2 (HB). Full curves represent the median spectral output across the 36 animals and shaded regions refer to the 25th and 75th percentiles.
flexible motor behaviour (Feingold et al., 2015;Sherman et al., 2016). Following this reasoning, a recognition of the mechanisms adopted by the CBGTC network to regulate beta spectral undulations is vital to better understand healthy and diseased states; and consequently, inform novel therapeutic strategies. Here, using DCM, we highlight a set of synaptic alterations in the thalamocortical loop that elucidate how the transitions of beta synchrony from low to high levels might occur in Parkinson's disease. We provide a new perspective for the effective coupling of the Parkinsonian thalamocortical network, where a fine regulation of temporally different inputs to specific laminae of the motor cortex may underlie the spontaneous and transient variability in oscillatory neural activity in the beta band across the circuit.

Conflicts of interest
None.