Resting state functional MRI in Parkinson’s disease: the impact of deep brain stimulation on ‘effective’ connectivity

Abstract

Depleted of dopamine, the dynamics of the parkinsonian brain impact on both ‘action’ and ‘resting’ motor behaviour. Deep brain stimulation has become an established means of managing these symptoms, although its mechanisms of action remain unclear. Non-invasive characterizations of induced brain responses, and the effective connectivity underlying them, generally appeals to dynamic causal modelling of neuroimaging data. When the brain is at rest, however, this sort of characterization has been limited to correlations (functional connectivity). In this work, we model the ‘effective’ connectivity underlying low frequency blood oxygen level-dependent fluctuations in the resting Parkinsonian motor network—disclosing the distributed effects of deep brain stimulation on cortico-subcortical connections. Specifically, we show that subthalamic nucleus deep brain stimulation modulates all the major components of the motor cortico-striato-thalamo-cortical loop, including the cortico-striatal, thalamo-cortical, direct and indirect basal ganglia pathways, and the hyperdirect subthalamic nucleus projections. The strength of effective subthalamic nucleus afferents and efferents were reduced by stimulation, whereas cortico-striatal, thalamo-cortical and direct pathways were strengthened. Remarkably, regression analysis revealed that the hyperdirect, direct, and basal ganglia afferents to the subthalamic nucleus predicted clinical status and therapeutic response to deep brain stimulation; however, suppression of the sensitivity of the subthalamic nucleus to its hyperdirect afferents by deep brain stimulation may subvert the clinical efficacy of deep brain stimulation. Our findings highlight the distributed effects of stimulation on the resting motor network and provide a framework for analysing effective connectivity in resting state functional MRI with strong a priori hypotheses.

Introduction

Progressive asymmetric degeneration of nigrostriatal dopaminergic innervation is a primary hallmark of Parkinson’s disease. Clinically, this typically produces asymmetric motor symptoms that impact on behaviours both ‘during action’ and ‘at rest’. Of this latter group, rigidity and resting tremor are most common. Chronic high frequency deep brain stimulation (DBS) of the subthalamic nucleus (STN) has become an established therapy for managing these symptoms, when dopaminergic medications alone are no longer sufficient (Limousin et al., 1995; Krack et al., 2003), primarily improving ‘OFF’ periods. Although most conventional therapeutics aim to restore dopamine concentrations to physiological levels, the mechanism of action of STN-DBS is less clear. As the effect of DBS mimics that of ablative lesions, it was suggested that DBS ‘inhibits activity’ in the target, which accorded with rate-based models of basal ganglia circuits (Albin et al., 1989; DeLong, 1990; Beurrier et al., 2001; Meissner et al., 2005). However, the literature—ranging from animal to computational models—suggests that stimulation has a myriad of effects on various neural elements in and around the STN, culminating in clinical improvement (Perlmutter and Mink, 2006; Deniau et al., 2010; McIntyre and Hahn, 2010; Vedam-Mai et al., 2012).

Blood oxygen level-dependent (BOLD) functional MRI signal recorded while the subject lies at rest with eyes closed could represent an important tool for clinical diagnosis and understanding brain disorders (Biswal et al., 1995; Fox and Greicius, 2010; Deco et al., 2011). Brain regions of similar functional specialization have been shown to display functional connectivity (i.e. correlated BOLD signal) during rest (e.g. motor and visual networks, respectively; Biswal et al., 1995; Lowe et al., 1998). Parkinsonian patients OFF medication have been shown to display reduced BOLD functional connectivity amongst the pre/motor cortex and putamen (Wu et al., 2009; Esposito et al., 2013), increased connectivity between the motor cortex (M1) and cerebellum (Wu et al., 2009), reduced striato-thalamic connectivity (Hacker et al., 2012) as well as increased M1–STN connectivity (Baudrexel et al., 2011). Electrophysiological studies report similar changes in coherence between the STN and motor cortical areas in resting parkinsonian patients (Shimamoto et al., 2013) that can be restored by dopaminergic medication (Litvak et al., 2011). Thus, there is a wealth of evidence implicating dopamine’s neuromodulatory role in resting cortico-subcortical circuits.

When investigating neuromodulatory effects of DBS at rest, previous human imaging studies have largely used either PET or single photon emission tomography. These studies have demonstrated STN-DBS-induced changes in both blood flow and glucose metabolism at rest in key constituents of the motor network including M1, putamen, and thalamus (Limousin et al., 1997; Ceballos-Baumann et al., 1999; Jech et al., 2001; Hershey et al., 2003; Stefurak et al., 2003; Payoux et al., 2004; Haslinger et al., 2005; Asanuma et al., 2006; Grafton et al., 2006; Hilker et al., 2008; Karimi et al., 2008; Cilia et al., 2009; Geday et al., 2009; Boertien et al., 2011). Imaging in animal models has yielded similar results (Min et al., 2012; Lai et al., 2013). One short report described increased functional connectivity among premotor regions in response to STN-DBS (Mueller et al., 2013); however, functional MRI has been limited in these patients because of safety concerns. We have previously demonstrated that any risk to the patient can be minimized under a strict acquisition protocol (Carmichael et al., 2007; Kahan et al., 2012).

Functional MRI has clear advantages over PET and single photon emission tomography in view of the post processing analytical methods that can be used. Dynamic causal modelling (DCM) is a Bayesian scheme typically used to explain connectivity changes underlying task-related brain responses (Friston et al., 2003). Conceptually, DCM usually treats the brain as a ‘deterministic system’, meaning the response in a region is determined entirely by inputs to that region (for a brief clinician-friendly introduction to DCM; Kahan and Foltynie, 2013). Deterministic DCM, therefore, depends on experimental stimuli driving the system. The investigator typically identifies regions of interest that respond to these stimuli, extracts time-series that best summarize their activity, and then constructs a series of models representing competing hypotheses about the underlying functional architecture. The connectivity parameters are then estimated based on the model structure. Models are compared using Bayesian model selection and those optimized parameters reported (Penny et al., 2004). In contrast to functional connectivity studies, DCM estimates ‘effective connectivity’, that is, how the activity in one region causes activity changes in another, thereby modelling how information passes through networks. This permits a mechanistic understanding of observed data; in other words, what is causing the change we witness in the data. For example, from neurophysiological recordings, the spectral shifts in Parkinson’s disease-related β-band responses within the basal ganglia motor loop can be explained by increased connectivity of STN afferents from M1 with accompanying changes in STN-pallidal activity (Moran et al., 2011; Marreiros et al., 2012).

Resting state functional MRI data do not possess any exogenous (driving) inputs as in task-based studies. Thus the evolution of activity in a given region must be driven by endogenous fluctuations. Stochastic DCM accounts for this by including stochastic fluctuations in the differential state equations and uses Bayesian filtering to estimate the hidden neuronal states, coupling parameters and the precision of observation noise (Friston et al., 2011; Li et al., 2011; Daunizeau et al., 2012a, b).

In this work, we deconstruct the basal ganglia motor loop using stochastic DCM of functional MRI data acquired while human patients lay at rest, and examine the effect of therapeutic STN-DBS on the underlying effective connectivity. We finesse methodological issues pertaining to region selection in resting state studies, permitting analysis of data that are both functionally verified and anatomically localized. Given the clear clinical response of these patients to DBS at rest, we hypothesized that clinical improvement could be explained by changes in the way neural populations within the motor loop impact upon one another; in other words, DBS has modulatory effects on ‘extrinsic’ (between region) effective connectivity. Specifically, we sought to answer whether the effect was precisely related to particular pathways of the cortico-basal ganglia-cortical loop, or whether there was a diffuse effect on all connections. Furthermore, we explored whether clinical measures of Parkinson’s disease impairment could be explained by differences in extrinsic effective connectivity, within the basal ganglia motor loop.

Materials and methods

This study was approved by the National Hospital and Institute of Neurology Joint Ethics committee (09/H0716/51). All participants provided written informed consent.

Patients

Twelve patients who met UK brain bank criteria for idiopathic Parkinson’s disease were studied (Table 1). Patients had received chronic bilateral STN-DBS for >6 months. Stereotactic MRI, for both preoperative targeting and immediate postoperative verification was used during surgery (Foltynie and Hariz, 2010), and ensured electrode contacts were well-sited within the STN. All patients received bilateral electrodes (Model 3389, Medtronic) and a dual channel pacemaker (‘implanted pulse generator’ Kinetra^TM, Medtronic) implanted. Stimulation parameters were set to produce optimal clinical responses. Medication was withdrawn for 10–12 h (overnight) before scanning. Inclusion was limited to those patients who could tolerate lying flat with minimal head tremor while being both OFF medication and off stimulation.

Unified Parkinson’s Disease Rating Scale part III (UPDRS-III) scores (clinical measure of Parkinson’s disease motor impairment; higher score confers greater impairment) were recorded both on and off stimulation before scanning. Additionally, stimulation parameters and system impedance were noted, and implanted pulse generator counters were reset. The pacemaker monitors how many times it has been switched on and off. We reset the counters before scanning, and checked the counter after scanning to ensure the pacemaker was not turned on or off during scanning.

Magnetic resonance imaging data acquisition

Following on-site tissue-equivalent test-object thermometry experiments confirming that (under strict protocol) sequences used in functional MRI studies posed no risk to the patient (Carmichael et al., 2007), scanning was performed in a Siemens Avanto 1.5 T MRI scanner using a transmit-receive (Tx/Rx) head coil. The specific absorption ratio in the head was limited to <0.1 W/kg.

Subjects received three functional MRI scans during each stimulation condition: (i) resting state with eyes closed (repetition time = 2420 ms; echo time = 40 ms; flip angle = 90°; field of view = 192 × 192 mm²; matrix size = 64 × 64; 32 axial slices 3.5 mm thick, gap between slices of 0.7 mm; spatial resolution = 3 × 3 × 4.2 mm³; duration = 8 min; 200 scans); (ii) motor task (right hand); and (iii) motor task (left hand). A vacuum moulded cushion was used to securely support the head, and limit head movement. Patients had an alarm if they experienced discomfort. The motor task has previously been presented in detail (Kahan et al., 2012). In brief, subjects lay in the scanner with a joystick in one hand and were instructed to move the joystick in a random direction in response to auditory tones during ‘go’ blocks, and to ignore the tones and rest their hand on the joystick during ‘rest’ blocks. The task was repeated for both hands, serving as functional localizers (one for each hemisphere) for subsequent resting state analysis.

Patients were scanned both during active therapeutic (ON) and inactivated stimulation (OFF). The scanning order was counterbalanced across subjects. Thus, two resting state sessions were collected per patient, one for each stimulation condition, in addition to two field map scans, an anatomical T₁-weighted MP-RAGE structural scan, and the aforementioned functional (motor) localizers.

After scanning, active stimulation was restored, normal medication was administered, and a UPDRS-III examination was repeated to confirm patients had returned to their clinical baseline. The settings, counters and impedance of the DBS system were recorded to confirm there were no additional activations induced by the scanner.

Dynamic causal modelling of resting state functional MRI data

Dynamic causal modelling is a framework for fitting differential equation (state space) models of neural states to neuroimaging data using Bayesian inference (Friston et al., 2003). We used DCM to infer how observed resting state BOLD data from pre-specified regions of the brain were generated from an underlying network of interdependent neural populations (nodes). Importantly, DCMs are fully generative models that are quantitative descriptions of how the data were generated from neural activity (‘hidden states’) that cannot be recorded directly. Directed interactions between and within regions, i.e. effective connectivity, occurs at the neural level (Marreiros et al., 2010; Kahan and Foltynie, 2013). Notably, connections between regions are not necessarily monosynaptic, thus not all regions of a network are strictly necessary to specify a complete model. The hidden neural states enter a haemodynamic forward model (with a number of haemodynamic parameters), resulting in a predicted BOLD signal that would be observed if the model parameters were true. Parameters are estimated iteratively using a standard Bayesian (generalized filtering) scheme, maximizing the fit between observed and predicted data under model complexity constraints. Models are scored in terms of ‘model evidence’ (more accurately, an approximation of the model evidence called variational free energy), and is a compromise between model accuracy and model complexity, avoiding bias towards over-fitting models that are less generalizable. Models of the same data can thus be compared on the basis of their model evidence using Bayesian model selection (Penny et al., 2004; Stephan et al., 2009), allowing the investigator to compare different hypothetical functional architectures.

Traditionally, DCM models the brain as a ‘deterministic system’, meaning a change in activity in one region is determined completely by its inputs, forgoing ‘endogenous’ activity that may characterize certain brain regions, or inputs from sources not included in the model. A recent extension (stochastic DCM for functional MRI) incorporates spontaneous noisy fluctuations into the state equations, allowing one to model dynamics without any experimental inputs, for example, in the ‘resting state’ (Li et al., 2011; Daunizeau et al., 2012a, b).

In this work, we constructed a series of models of the basal ganglia motor loop consisting of M1, putamen, motor thalamus, and STN (hidden node, see below) in a single hemisphere. We hypothesized that STN-DBS has modulatory effects on extrinsic connections within the network, and thus we constructed a model space (for each hemisphere of each subject), where models differed in terms of which connections were modulated by DBS. The artefact caused by the electrode precluded the acquisition of precise BOLD data from the STN region, and this loss of precision was included in our model. This means the STN node can be regarded as ‘hidden’ from measurement and is referred to as a ‘hidden node’ (David et al., 2011).

The processing stream is summarized in Fig. 1. Processing and analysis was performed using SPM12 (Wellcome Trust Centre for Neuroimaging, London, UK; http://www.fil.ion.ucl.ac.uk) and DCM12. Anatomical masks of the precentral gyrus from each hemisphere were created using the Harvard-Oxford cortical structural atlas, as available in the FSL suite (http://fsl.fmrib.ox.ac.uk/). Masks from each hemisphere of the motor putamen and motor thalamus were created using probabilistic white matter connectivity atlases thresholded at 50% probability (Behrens et al., 2003; Tziortzi et al., 2013), constraining our analysis to regions that exhibit strong structural connectivity at a population level. For the purpose of this analysis, we considered both hemispheres independently, given (i) Parkinson’s disease causes asymmetric degeneration of the substantia nigra pars compacta; (ii) severity of motor symptoms are usually asymmetric; and (iii) DBS settings are often asymmetric.

Functional localization of subject-specific M1

Using the task data, a contrast defining the main effect of movement for each hand was specified, and the peak voxel of this contrast (constrained by an anatomical mask of the precentral gyrus contralateral to hand movement) was selected as the hemisphere’s M1 coordinate. Analysis of the task data is discussed elsewhere (Kahan et al., 2012). Given DCM is ultimately a single subject analysis repeated on many subjects, problems of overlapping loss-of-signal artefact produced by the DBS hardware were avoided (Kahan et al., 2012).

Preprocessing

The first five scans of each resting state session were removed and data were corrected for field inhomogeneity using the field maps. Data were then realigned, coregistered, anatomical images were normalized to MNI space, the resultant normalization matrix was then used to normalize the functional data. Finally, the data were visually inspected and spatially smoothed using an 8 mm Gaussian kernel. This specific preprocessing routine was used to ensure that any artefacts produced by the DBS hardware did not impact normalization to MNI space. The ON and OFF sessions were concatenated to produce a single 390 scan time-series. Ultra-low frequency fluctuations were removed in the usual way using a high-pass filter (1/128 s, ≈0.0078 Hz). Confound time-series were extracted from predefined coordinates of extra-cerebral compartments (the lateral ventricle and eye globe).

General linear model of resting state dynamics

We extracted data exhibiting physiologically-relevant resting state (i.e. low frequency) dynamics from functionally specialized motor regions. The resting state was thus modelled using a General Linear Model with a discrete cosine basis set (GLM-DCT) consisting of 189 functions with frequencies characteristic of resting state dynamics (0.0078–0.1 Hz; Biswal et al., 1995; Fransson, 2005; Fox and Raichle, 2007; Deco et al., 2011), a regressor encoding the effect of DBS, six nuisance regressors from each session capturing head motion, and the confound time-series from the extra-cerebral compartments. The regional BOLD signal was summarized with the principal eigenvariate (adjusted for confounds: head movements and extra-cerebral compartments) of voxels within 4 mm of the subject’s M1 coordinate, as identified using statistical parametric mapping of motor sessions. For those familiar with the process of extracting volume(s) of interest, this was achieved by using an F-contrast including the effect of stimulation regressor, as well as the discrete cosine set modelling the resting state. This procedure allowed us to extract physiologically relevant resting state data from the functionally defined M1 for each hemisphere.

To identify voxels within the subcortical nuclei that were sensitive to DBS, we used a psychophysiological interaction analysis (Friston et al., 1997). The corresponding general linear model included the main effects of DBS regressor, the BOLD activity from M1 and their interaction. Statistical parametric mapping testing for the psychophysiological interaction was masked first with the putamen mask, and then the thalamus mask (of the same hemisphere). The BOLD signal (corrected for the same confounds as above) was extracted from a sphere (radius 4 mm) centred on the peak T-value within each mask, producing three volumes of interest per hemisphere and subject (M1, putamen, thalamus). Here, the putamen and thalamus volumes of interest exhibit a psychophysiological interaction with M1 and DBS. BOLD data from the STN could not be considered because of its small size and loss-of-signal artefact caused by the DBS electrode.

Model space and comparison

Volumes of interest from each hemisphere were used to construct a series of 32 DCMs representing different hypothetical architectures. Two-state (Marreiros et al., 2008), stochastic DCM for functional MRI was used, endowing each node with excitatory and inhibitory subpopulations in receipt of noisy fluctuating inputs or endogenous activity. The STN was included as a hidden node, whose noise precision (given the electrode artefact) was effectively zero, permitting estimation of its hidden states and coupling parameters in the normal way. This reflects a strength of DCM; the Bayesian inversion of these models solves the complex problem of estimating hidden or latent variables. In DCM all parameters of interest are hidden and their expression in data serves to estimate these latent variables, regardless of whether their effects on the data are local or distributed. This has been capitalized on in electrophysiological DCMs (David et al., 2011; Moran et al., 2011; Marreiros et al., 2012) and is exploited here in DCM for functional MRI.

The globus pallidus pars internus and globus pallidus pars externus were not included in our models; rather it was assumed that thalamic afferents arrived through GABAergic projections from the globus pallidus pars internus. In other words, striatal GABAergic medium spiny neurons projecting directly to the globus pallidus pars internus (constituting the direct pathway) were modelled with a net excitatory effect on the thalamus, whereas the glutamatergic STN-thalamic inputs had a net inhibitory effect (Fig. 2A). Striatal afferents to the STN summarize the polysynaptic putamen-globus pallidus pars externus-STN connection that produces a net excitatory effect on the STN in accordance with electrophysiological findings (Kitai and Deniau, 1981; Rouzaire-Dubois and Scarnati, 1985; Albin et al., 1989; DeLong, 1990; Nambu et al., 1996, 2000; Gradinaru et al., 2009; Kravitz et al., 2010). The direct pathway was thus defined as the striato-thalamic excitatory connection, and the indirect pathway thus consisted of both the striato-STN and STN-thalamic connections. Finally, the hyperdirect pathway was defined as the M1-STN connection. The effect of DBS entered the models by modulating a subset of connections. The model space compared DCMs that included modulatory effects on the direct, indirect, hyperdirect, cortico-striatal, or thalamo-cortical pathways, or combinations of those five pathways (comprising six connections), resulting in (2⁵) 32 models per hemisphere (Fig. 2).

Models were inverted using generalized filtering (Friston et al., 2010; Li et al., 2011), providing an estimate of the coupling parameters and model evidence. Generalized filtering is a Bayesian filtering scheme for non-linear state–space models in continuous time, i.e. dynamic causal models (Friston et al., 2010). The 32 models from each of the 24 hemispheres entered a Bayesian model selection procedure (fixed effects assumptions; see Supplementary material) that computed the posterior probabilities over competing models (Penny et al., 2004; Stephan et al., 2009). Models were subsequently grouped into families depending on whether they expressed modulatory effects on the five pathways discussed above. A post hoc Bayesian model selection family analysis was used to evaluate the posterior probabilities of a modulatory effect on each (set) of the pathways.

Relationship between winning dynamic causal modelling and clinical response to deep brain stimulation

To quantify and validate the effective connectivity estimates generated, classical regression models were then used to test for the effect of DBS at the between subject level and as predictors of clinical status and therapeutic response to DBS. The model with the greatest log-evidence at a group level (i.e. across hemispheres) was considered the winning model. Extrinsic connectivity values (and their DBS dependent modulation) from the six connections comprising the five pathways of interest were taken from the winning model of each hemisphere. The on coupling for each parameter was calculated using the off coupling (DCM A-matrix values), and the DBS modulatory effects (DCM B-matrix); see Supplementary material for the parameterization. Paired t-tests were then used to compare coupling parameters in the on versus off conditions.

The coupling parameters were subsequently entered into a multiple linear regression model, as independent variables, with the contralateral UPDRS-III score (excluding axial score) as the dependent variable. This was performed separately for the on and off conditions. The direction of the coupling (i.e. excitatory or inhibitory) was not considered in the model, we simply asked whether the strength of the coupling predicted impairment. Additionally, the DBS-induced modulation was calculated for each connection (on coupling strength / off coupling strength), thus modulations >1 denoted DBS increased the coupling, regardless of whether the connection targeted excitatory or inhibitory subpopulations (see Supplementary material for details about parameterization). These DBS effects entered a final stepwise regression model (employing backwards elimination) as independent variables, predicting the percentage improvement in clinical phenotype (i.e. the larger the percentage, the greater the clinical efficacy).

Results

Patients

Scanning proceeded with no adverse effects or change in post-scan UPDRS-III scores, implanted pulse generator function or unexpected activations, or change in circuit impedance. Reintroduction of medication led to restoration of baseline motor function. Detailed patient and DBS parameter information can be found in Table 1. The mean clinical improvement to DBS OFF medication was 23.8 points [95% confidence interval (CI): 16.9–37.8, P < 0.0001]. When considered by hemisphere, and ignoring axial subscores, DBS improved contralateral signs in the limbs by 8.9 points (95% CI: 6.8–11.0, P < 0.000001).

Model fit

Models from two hemispheres (from two separate subjects) failed to fully converge during DCM’s model inversion procedure. The time-series from these problematic analyses were scaled down by a factor of 2 to suppress high amplitude ‘spikes’. Following this, all 32 models across the 24 individual hemispheres were inverted successfully, furnishing predicted BOLD time-series for each of the volumes of interest. These were similar to the observed signal. A representative plot showing the observed and predicted data from one model is provided in Fig. 3.

Bayesian model selection

Fixed effects Bayesian model selection across 24 hemispheres revealed that Model 32 had the greatest relative log-evidence, beating its closest competitor by >17. This means the posterior probability (i.e. the likelihood of the model generating the data, given the model space) of the winning model was >99% in relation to its nearest competitor (Fig. 4A). Model 32 included modulatory effects on all five of the pathways explored. Subsequent family analysis—where models were grouped by the presence of modulatory effects on the five pathways—confirmed that all five pathways were >99% likely to be modulated by DBS (Fig. 4B).

Direction of neuromodulatory effect

Estimated coupling parameters from the six extrinsic connections were extracted from Model 32 of each hemisphere. Paired t-tests were used to compare the means of the connection strengths on and off DBS (Fig. 5). Stimulation increased the strength of the cortico-striatal (95% CI: 0.01–0.06 Hz, P < 0.05), direct (95% CI: 0.03–0.07 Hz, P < 0.001), and thalamo-cortical pathways (95% CI: 0.03–0.07 Hz, P < 0.001). In contradistinction, STN DBS reduced the strength of all STN afferents and efferents; the hyperdirect (95% CI: −0.0011 to −0.0008 Hz, P < 0.001), striatal afferents (95% CI: −0.0010 to −0.0006 Hz, P < 0.001), and STN-thalamic (95% CI: −0.0008 to −0.0005, P < 0.001) connections. The magnitudes of changes to connections involving the STN (i.e. connections constituting the hyperdirect and indirect pathways), although highly significant, were markedly smaller than the other modulatory effects discovered. This was because of the fact that the STN was estimated as a hidden region (see ‘Materials and methods’ section). All P-values are corrected for multiple comparisons using the Bonferroni procedure.

Connection strengths predict clinical impairment

Separate on and off multiple linear regression analyses revealed consistent results with regards to the direction of the regression coefficient for each connection (Fig. 6). Three connections showed significant predictive capability in both stimulation conditions; the hyperdirect, direct, and striatal-STN connections. As the strength of the direct pathway increased, clinical impairment was reduced. Surprisingly, the same was true of the hyperdirect pathway, despite STN-DBS reducing the strength of this connection. The reverse was true of the striato-STN pathway; the stronger this pathway, the more disabled the patient.

Change in connection strengths predict clinical efficacy

The three connections described above were the only connections of the original six that were included in the most parsimonious (after backwards elimination) regression model (Fig. 7). DBS scaling of coupling parameters were coded such that scaling >1 conferred increased coupling, whereas <1 conferred reduced coupling; thus as scaling increased numerically, the coupling got stronger. Consistent with the previous regression, increased scaling of the direct pathway predicted increased efficacy (P < 0.05). Increased scaling of the hyperdirect pathway also predicted increased efficacy (P < 0.05), whereas increased scaling of the basal ganglia afferents predicted reduced efficacy (although this was only trend significant, P = 0.067).

Discussion

Widespread neuromodulation

In this work, we used a well-established (if simplified) model of the functional architecture of the motor cortico-striato-thalamic loop and estimated which connections, or combinations of connections, are modulated by STN-DBS to produce the observed BOLD signal recorded at rest. Using Bayesian model selection, we were able to demonstrate that a model with DBS-related modulatory effects on the extrinsic direct, indirect, hyperdirect, thalamo-cortical, as well as cortico-striatal pathways, consistently out-performs other (plausible) models. Analysis of the connection strength parameters revealed significant differences between active and inactive STN-DBS in the six extrinsic connections investigated. The relevance of extrinsic coupling parameters to clinical phenotype was subsequently demonstrated using two orthogonal regression models; extrinsic coupling predicted clinical phenotype, and changes in extrinsic coupling predicted DBS-related clinical improvements (summarized in Fig. 8).

Coupling parameters in DCM reflect the sensitivity of a target to its afferent. Modulatory effects on coupling can thus be conceptualized as an afferent-specific gain modulation of the target; in other words, the target responds to the afferent in question more/less during stimulation. Our finding that STN-DBS increases the sensitivity of the putamen to cortical afferents is interesting given that Parkinson’s disease is known to reduce striatal medium spiny neuron dendritic spine density and impair cortico-striatal glutamatergic transmission (Garcia et al., 2010), and computational approaches have implied dopamine potentiates cortico-striatal synaptic strength (Leblois et al., 2006). Neuroimaging studies have demonstrated increased metabolism in the striatum under STN-DBS (Hilker et al., 2002; Geday et al., 2009). Animal models of STN-DBS have suggested that stimulation normalizes cortico-striatal glutamatergic synaptic densities and distributions in 6-hydroxydopamine-lesioned rats (Walker et al., 2012), and regularizes striatal discharge patterns in the resting 1-methyl-4-phenyl-1,2,3,6- tetrahydropyridine (MPTP)-lesioned macaque (Santaniello et al., 2012). It remains possible that such enhancements could contribute to clinical improvements; however, the modulation of this connection was excluded from our parsimonious model of stimulation efficacy.

We found that cortical sensitivity to thalamic afferents showed even greater, highly significant enhancements. PET and single photon emission tomography functional imaging in STN-DBS patients have previously documented reduced perfusion and/or metabolism in the resting motor cortex during stimulation, and several authors have suggested this is caused either by restoration of pallidal inhibition of thalamic outflow or reduced synchronized oscillatory activity between motor cortex and basal ganglia (Limousin et al., 1997; Ceballos-Baumann et al., 1999; Hershey et al., 2003; Payoux et al., 2004; Haslinger et al., 2005; Grafton et al., 2006; Cilia et al., 2009; Geday et al., 2009). Our results are not necessarily in disagreement: increased effective connectivity between the thalamus and cortex does not necessarily confer greater metabolic activity at the target, rather it suggests that during stimulation, while at rest, the cortex may respond more efficiently to afferents from thalamus. Nevertheless, changes in this connection similarly failed to predict clinical efficacy, thus this effect may not be key to the immediate therapeutic mechanism of STN-DBS.

Enhanced thalamic sensitivity to the direct pathway

STN-DBS increased the effective connectivity of the direct pathway, summarized here as a net excitatory connection from the putamen to the thalamus. Additionally, the strength of this connection predicted clinical severity in both on and off DBS states, and its scaling predicted clinical improvement; stronger direct pathway effective connectivity diminished Parkinson’s disease impairments. This finding is highly concordant with early rate-based models of basal ganglia function that propose dopamine depletion results in an underactive direct pathway, and an overactive indirect pathway, culminating in thalamic inhibition (Albin et al., 1989; DeLong, 1990). Optogenetic investigation in 6-hydroxydopamine-lesioned rats has lent further evidence to the beneficial effects of enhancing the direct pathway (Kravitz et al., 2010). The mechanism underlying this proposed gain tuning cannot be discerned here; our models do not include the globus pallidus pars internus, thus, we are unable to confirm whether modulation takes place further upstream of the thalamus. However, MPTP-lesioned monkeys receiving STN-DBS display a more regular pattern of neuronal firing in the pallidal-receiving thalamic nuclei (Xu et al., 2008), and reduced neuronal entropy (Dorval et al., 2008), suggesting that STN-DBS has downstream effects as far as the thalamus.

Subthalamic nucleus afferents and deep brain stimulation

The role of the hyperdirect pathway in health and Parkinson’s disease remains unclear (Nambu et al., 2002; Nambu, 2005; Leblois et al., 2006; Baudrexel et al., 2011). Parkinson’s disease is electrophysiologically characterized by enhanced β-band oscillations in the STN, basal ganglia nuclei and motor cortex (Brown et al., 2001). The hyperdirect pathway has been implicated in the generation of this biomarker: First, simultaneous STN local field potential and M1 electrocorticography suggest an exaggerated phase-amplitude coupling between cortical β-phase and broadband γ-amplitudes in Parkinson’s disease. Importantly, cortical γ was shown to precede STN β, potentially suggesting hyperdirect drive maintains STN β. This altered coupling was reduced by STN-DBS (de Hemptinne et al., 2013). Furthermore, single-unit recordings in 6-hydroxydopamine-lesioned rats receiving STN-DBS have shown antidromic hyperdirect activation of M1 layer V neurons, modifying their firing probability, and suppressing β-synchrony at a population level (Li et al., 2012). STN local field potential and cortical MEG have similarly shown that cortico-STN coherence is predominantly cortex leading (Litvak et al., 2011), and optogenetic stimulation of the hyperdirect pathway in 6-hydroxydopamine-lesioned rats has been shown to ameliorate parkinsonian symptoms, supporting the idea that interrupting hyperdirect cortico-STN transmission is therapeutic (Gradinaru et al., 2009).

Secondly, DCMs of local field potential data in 6-hydroxydopamine-lesioned rats (Moran et al., 2011) and human patients (Marreiros et al., 2012) suggest that parkinsonian spectral patterns at β frequencies can be explained by increased effective connectivity of the hyperdirect pathway. Both these results support the hypothesis that STN β oscillations are caused (in part) by increased effective drive from the cortex, something that we found DBS to decrease. Notably, in these studies, DCM was used to explain observed increased β-band power at the nodes interrogated, whereas BOLD data (such as ours) are not known to possess a reliable correlate of this electrophysiological biomarker. Concurrent results from stochastic DCM for functional MRI data are encouraging (and rather remarkable), and suggest DCM is able to furnish similar underlying functional architectures from data sets with markedly different temporal structures. This is potentially valuable, especially when considering the non-invasive nature of functional MRI.

The relationship between hyperdirect coupling strength and clinical severity is particularly interesting. Comparing coupling on and off stimulation reveals results that are in line with previous studies of effective connectivity, i.e. stronger coupling is associated with parkinsonian phenotypes (Moran et al., 2011; Marreiros et al., 2012). However, when all six connections are considered in a regression model predicting impairment including tremor, rigidity and bradykinesia, as rated by the UPDRS-III in humans, stronger effective hyperdirect coupling is beneficial to patient symptoms. Further studies are needed to confirm our findings; however, considering hemibody scores—as we have done—may provide a more comprehensive evaluation of clinical effects that have not previously been captured by animal models or simple on versus off comparisons. Taken together, our results suggest that the most beneficial effects of stimulation on the resting motor system seem to be explained by strengthening the effective coupling along the direct pathway, and not reducing coupling along the hyperdirect pathway. Of course, we are not able to address whether the amelioration of the hyperdirect pathway is predictive of unwanted effects (or improvements) that are not indexed by the UPDRS-III score, and this issue could be the subject of further investigation. Disruption of the hyperdirect pathway has previously been used to explain increased impulsivity observed in STN-DBS patients during high-conflict decision making tasks (Frank et al., 2007); however, this has been studied with medial prefrontal projections as opposed to M1 efferents (Cavanagh et al., 2011). Our work suggests that disruption to hyperdirect effective connectivity may not be limited to prefrontal projections. Further work would be required to establish whether disruption to the hyperdirect pathway depends on the cortical source. While DBS may modulate the effective connectivity of the hyperdirect pathway in a way that may not be optimal, we acknowledge that the hyperdirect pathway may still be important in mediating some of the physiologically beneficial changes seen in the circuit.

The indirect pathway displayed similar significant modulation. Specifically, DBS reduced the effective connectivity of STN afferents arising from the striatum, and efferents to the thalamus. These pathways summarize the known polysynaptic connections between these structures (through the globus pallidus pars externus for STN afferents, and globus pallidus pars internus for STN efferents), although they do not preclude a direct influence. Two-state DCM for functional MRI enabled us to place prior constraints on the direction (excitatory or inhibitory) of the coupling (Marreiros et al., 2008), in agreement with a wealth of electrophysiological and anatomical tracing literature. It is difficult to compare this finding to the existing literature on the STN afferents from the globus pallidus pars externus, especially as these two nuclei possess dense reciprocal connections (Moran et al., 2011; Marreiros et al., 2012), although it does suggest that indirect pathway afferents are an important determinant of clinical severity. Unlike the hyperdirect pathway though, this scaling was in agreement with the direction of neuromodulation induced by stimulation.

Most of the STN-DBS neuroimaging literature has tried to explain the benefits of DBS in terms of changes in regional activity. Regarding the resting motor circuit, although results appear heterogeneous, most conclude STN-DBS increases thalamic activity (Ceballos-Baumann et al., 1999; Jech et al., 2001; Hershey et al., 2003; Hesselmann et al., 2004; Phillips et al., 2006; Karimi et al., 2008), while simultaneously decreasing motor cortical activity (Limousin et al., 1997; Ceballos-Baumann et al., 1999; Hershey et al., 2003; Payoux et al., 2004; Haslinger et al., 2005; Grafton et al., 2006; Cilia et al., 2009; Geday et al., 2009), suggesting respective increases and decreases in synaptic activity. In the context of this literature, we propose that STN-DBS causes these regional changes through afferent specific gain modulation at key neuronal populations composing the basal ganglia motor loop. Computational accounts have proposed that dopamine encodes the precision of sensorimotor prediction errors by modulating the post synaptic gain of these afferents (Friston et al., 2012). In the absence of dopamine, STN-DBS may serve to partially restore this computational function.

Given the wealth of evidence for more global changes in brain activity under STN-DBS (e.g. dorsolateral prefrontal cortex, Cilia et al., 2009; supplementary motor area and parietal cortex, Hershey et al., 2003; anterior cingulate cortex, Limousin et al., 1997), it is possible that our findings may generalize to other loops of the basal ganglia (interestingly, mostly reducing their cortical metabolic signatures at rest in PET/single photon emission tomography studies). Additionally, higher order afferents to the motor cortex (e.g. premotor and prefrontal input) could also be subject to modulation, as well as other subcortical nuclei (Asanuma et al., 2006). However, finding any behavioural correlates of changes within these circuits at rest may be more difficult.

Limitations and model assumptions

Our model makes a number of simplifying assumptions, most notably the independence of the right and left hemisphere basal ganglia motor loops. As discussed above, this was motivated by the inherent clinical asymmetry observed in patients. However, there is evidence to suggest that components of the loops are functionally connected across hemispheres (De Solages et al., 2010). Furthermore, this analysis assumes that DBS of one hemisphere is no more effective than DBS of another (i.e. there is no STN-dominant hemisphere).

It is important to note that we model DBS as a modulatory effect on extrinsic coupling, not as a driving input to individual nodes. Thus, this work specifically addresses the effects of DBS on extrinsic coupling; it does not address how these changes are mediated or delivered to the therapeutic targets. For example, as discussed above, STN-DBS has been found to induce antidromic effects on the cortex (Li et al., 2012). Although this finding may explain how DBS reaches its target, it does not mean there is a DBS-dependent effective connection from the STN to cortex—the ‘consequences’ of antidromic stimulation would be expressed orthodromically in cortical efferents (and afferents) to/from basal ganglia or other regions. In summary, we believe that DBS modulates functional integration within the motor loop, and we present evidence to support this hypothesis; however, these effects may be mediated by orthodromic effects on the thalamus and/or antidromic effects on the cortex or basal ganglia.

Another simplifying assumption is the lack of pallidal nodes in our network. As discussed, DCM does not necessarily quantify monosynaptic coupling, thus not all nodes are required to estimate effective connectivity. Furthermore, stochastic DCM for functional MRI uses computationally demanding estimation routines, where processing time scales with the number of nodes modelled.

As a result of signal drop-out around the electrode, it is not possible to record BOLD data from the STN itself. Therefore, we modelled the STN as a hidden node, enabling inference on its afferents and efferents based on the influence they exert on nodes from which precise recordings were available. This is standard practice in DCM of EEG data, where some sources are hidden or silent because they cannot be ‘seen’ by scalp electrodes (David et al., 2011; Marreiros et al., 2012). In principle, our hidden node could be any brain region with the connectivity fingerprint specified by the model (i.e. any brain region excited by both M1 and the putamen, and that exerts inhibition on the thalamus). Given the anatomical and electrophysiological literature on the functional anatomy of the basal ganglia, our hidden node was attributed to the STN. Including hidden nodes can reduce associated effect sizes (because the parameters of hidden nodes are not informed by empirical data and shrink to their prior expectations of zero).

The range of clinical severity was also fairly small. This was because of the inherent limitation of using DBS patients and placing them in MRI scanners; mildly affected patients with Parkinson’s disease do not receive DBS, and those who have severe symptoms during OFF (both medication and stimulation) periods are less likely to engage in MRI research. Furthermore, we did not consider axial symptoms in our regression models as they are intrinsically difficult to lateralize.

Conclusion

Our findings highlight the distributed effects of DBS on the parkinsonian resting motor network in human patients in a non-invasive manner, and are largely in agreement with invasive animal experiments. Our integration of clinical data, distinguishing our work from the animal literature, suggests that the hyperdirect, direct, and STN afferents arising from the striatum are the most important predictors of clinical improvement. Sensitizing the thalamus to direct pathway afferents, while simultaneously desensitizing the STN to its afferents appear to increase stimulation efficacy. Intriguingly, STN-DBS appears to achieve these effects on the direct and striato-STN pathways, but actually has the reverse effect on the hyperdirect pathway, potentially subverting its therapeutic potential. It is tempting to hypothesize that sparing this pathway of modulatory effects would improve the efficacy of STN-DBS; potentially ameliorating unwanted effects of stimulation.

Although previous work has focused on statistical dependencies (functional connectivity) between regions of the brain, and how these change with onset of disease/therapy (Deco et al., 2011), this study uses validated modelling to derive effective connectivity between regions; that is, how neural populations impact on one another (David et al., 2008; Daunizeau et al., 2012a). From a methodological standpoint, this work is an example of how researchers can use DCM to provide mechanistic insights into neurological/psychiatric diseases and their therapies in a non-invasive manner using functional MRI data acquired in the resting state.

Acknowledgements

The authors acknowledge the use of the UCL Legion High Performance Computing Facility (Legion@UCL), and associated support services, in the completion of this work.

Funding

This study was funded by the Brain Research trust (www.brt.org.uk). J.K. receives funding from The Astor Foundation, The Rosetrees Trust, and the MHMS General Charitable Trust. T.F., L.Z., M.H. and the Unit of Functional Neurosurgery are funded by the Parkinson’s Appeal and the Monument trust. T.F. also receives funding from Cure Parkinson’s Trust, Parkinson’s UK and the Michael J Fox Foundation. K.F. receives funding from the Wellcome Trust. The work was undertaken by UCL/UCLH, who receives a proportion of funding from the UK Department of Health’s NIHR Biomedical Research Centres funding scheme. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Supplementary material

Supplementary material is available at Brain online.

Abbreviations

Abbreviations
 
• BOLD

  blood oxygen level-dependent

 
• DBS

  deep brain stimulation

 
• DCM

  dynamic causal modelling

 
• M1

  primary motor cortex

 
• MPTP

  1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine

 
• STN

  subthalamic nucleus

 
• UPDRS

  Unified Parkinson’s Disease Rating Scale

References