Ethics statement

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

Patients

Ten PD patients took part in this study (Table 1). All patients had PD meeting UK brain bank criteria, and had received bilateral STN DBS for at least 6 months. Surgery had been performed using stereotactic MRI for both preoperative targeting and immediate postoperative verification of lead location prior to implantation of the extension cables and the implanted pulse generator (IPG) [1], [24].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Clinical Response to STN DBS: Details of 10 PD participants and their motor response to STN DBS, assessed immediately before the fMRI scan while off medication.

https://doi.org/10.1371/journal.pone.0050270.t001

Stimulation equipment

All patients had bilateral STN electrodes (model 3389, Medtronic, Minneapolis) and a dual channel IPG (Kinetra™ or ActivaPC™, Medtronic, Minneapolis) implanted. Stimulation parameters had been previously optimised according to clinical response. Inclusion in this study was restricted to those patients who (1) could tolerate lying flat while being both off medication and off stimulation, (2) exhibited minimal head tremor, and (3) demonstrated an immediate >35% improvement in UPDRS part 3 (UPDRS-III) off-medication score when stimulation was switched ON compared with OFF. Medication was withdrawn for 10–12 hours (overnight) before the scanning session.

Before scanning, (1) UPDRS-III motor scores were documented both ON and OFF stimulation (OFF was scored approximately 10 minutes after stimulation was stopped), (2) stimulation parameters and system impedance were recorded, and (3) IPG counters were reset.

Participants wore MRI compatible isolating headphones and held an MRI compatible joystick in one hand (Cambridge Research Systems, Kent, England: model No: HH-JOY-4. Angular range: 30 degrees (+/−15 degrees), Grip: 11.5×3 cm). The position of the joystick in time and space was recorded at a sampling rate of 20 Hz. During the task, participants were instructed to move the joystick consistently in response to auditory signals and to avoid excessively fast or large movements. Their heads were securely supported using a vacuum moulded cushion to dampen any head movement. Patients held an alarm in their non-moving hand to alert the clinical team if they experienced any discomfort during the scan. Patients were asked to keep their eyes closed throughout scanning.

The task was performed both with therapeutic stimulation active (ON), and again when their stimulation inactivated (OFF). During each stimulation condition, the task was performed twice, once with each hand. In other words, every patient performed a right and left hand movement task while stimulation was ON and OFF. The order of stimulation (ON versus OFF) and the movement (right versus left) were randomised over subjects. DBS was switched ON or OFF using the patients' own Access™ controller, which we ensured functioned normally within the MRI environment.

MRI data acquisition

All scans were performed with a Siemens Avanto 1.5T MRI scanner (Siemens, Erlangen, Germany) using a Siemens-supplied transmit-receive (T/R) head coil, similar to the one that detailed tissue-equivalent test-object thermometry experiments had been performed with in our previous safety study [23]. The specific absorption ratio (SAR) in the head was limited to under 0.1 W/Kg.

Each participant completed the session with four movement task time-series (sessions), one for each hand during each stimulation condition. This corresponds to a factorial design with three factors; task (movement versus no movement); laterality (right versus left) and stimulation (ON versus OFF). The whole session took approximately 90 minutes. The connection between the electrode lead and the extension cable, commonly sited above the left parietal bone caused a loss of signal artefact resulting in data not being acquired in left hemispheric sensorimotor areas. Given these regions were a priori regions of interest, particularly when examining right hand movements, we elected only to analyse the left hand movement data. Acquisition parameters were as follows:

1. T1 weighted magnetization-prepared rapid gradient-echo (MPRAGE) structural scan (repetition time TR = 1590 ms, echo time TE = 3.3 ms, inversion time TI = 1100 ms, flip angle = 15°, field of view FOV = 250×250 mm², matrix size = 192×192, 144 sagittal slices 1.3 mm thick, for a spatial resolution of 1.3 mm isotropic) lasting approximately 10 minutes. This scan, and an additional 8 minute resting scan (to be reported elsewhere), allowed a constant period of equilibrium to follow each patient's stimulation adjustment.
2. GE-EPI Movement session, Hand 1: (TR = 3695 ms, TE = 40 ms, flip angle = 90°, FOV = 192×192 mm², matrix size = 64×64, 49 axial slices 2.5 mm thick, gap between slices of 0.5 mm, for a spatial resolution of 3×3×3 mm³, 96 volumes, acquisition time = 6 minutes). The fMRI task paradigm consisted of 12 blocks lasting ∼30 seconds each. During each block, a series of 15 audio stimuli (beeps) were sounded through the headphones. The time between beeps was randomised to between 1–3 seconds. The blocks alternated between a “rest” and a “go”. At the beginning of each “rest” block the participant heard the word “rest” and was instructed to rest their hand on the joystick, ignore the beeps and keep still. At the beginning of each “go” block the participant heard the word “go” and was instructed to move the joystick in one of four random directions of their choice, and then return the manipulandum back to the central resting position. A single movement was defined as moving the joystick from its position of equilibrium and then returning the joystick back to this position. The exact timings of the beeps were also recorded.
3. GE-EPI Movement session, Hand 2: The joystick was then moved to the opposite hand, and acquisition 2 was repeated.

Additionally, field maps were acquired to correct for field inhomogeneity. Patients then had their stimulation switched to the opposite condition. The joystick was returned to the hand that had first performed the task and the aforementioned acquisitions were repeated.

At the end of the session, DBS was switched back ON if OFF during the second session, and the patient was examined (including a repeat evaluation of UPDRS-III). The DBS system was interrogated to check the settings and impedance, and to check for additional activations. The patients were given their regular PD medication and had a final clinical assessment after their medication had started to take effect to confirm they had returned to their baseline level of Parkinsonian disability before leaving the department.

Movement durations and reaction times were extracted from the joystick dataset. Paired t tests were used to judge significant changes in joystick movements comparing ON and OFF stimulation periods.

Image processing and regional BOLD signal analysis

Data was pre-processed and analysed using the SPM8 (Wellcome Trust Centre for Neuroimaging, London, UK; http://www.fil.ion.ucl.ac.uk); for effective connectivity analyses, DCM12 was used. The SPM Anatomy toolbox [25] was used to translate peak MNI coordinates into anatomical and functional regions based on probabilistic cytoarchitectonic maps.

Data were first corrected using the acquired field maps; accounting for field inhomogeneity caused both by the skull's air sinuses, and the implanted DBS equipment. When we examined the field maps, we noted that the extent and amplitude of the expected distortions caused by the DBS hardware on the skull are approximately the same as the extent and amplitude of the distortions caused by the presence of the sinuses. This suggests that the DBS hardware causes no more distortion than the sinuses do, and the field maps are sufficient to correct for them. Data were then realigned, correcting for motion. Each subject's T1 image was then normalised to MNI space. The normalisation matrix was then used to normalise the functional GE-EPI data. Each of the images was then visually inspected to confirm they had been correctly normalised. This order of normalisation ensured that the functional and anatomical scans were correctly translated into MNI space. The data were smoothed using an 8 mm Gaussian kernel, accounting for variation across subjects in structural and functional anatomy. Low frequency fluctuations were modelled using a high-pass filter set to the standard threshold (128-s).

Standard SPM (whole brain) statistical analyses were then performed using an epoch-related design, where each activation epoch (block) was defined as the time period from the beginning of the first movement in a “go” block, to the end of the last movement in that block. Each movement session thus consisted of six motor epochs, corresponding to the six “go” blocks. The resulting boxcar task function was then convolved with a canonical haemodynamic response function to form expansionary variables or regressors that constitute the design matrix. Both (ON and OFF) movement sessions, for each participant, were analysed in one design matrix. Six nuisance regressors were included for each session modelling the confounding effects of head motion in the design matrix.

We performed a standard random effects analysis by first computing contrasts of effects at the “first level” (within subject) and then analysing these summary statistics at the “second level” (between subjects) using one sample t-tests. Intrinsic masking was used to exclude voxels affected by DBS hardware-related artefact.

We examined for the contrast corresponding to the main effect of movement to (1) ensure that this could be detected in the DBS setup with a suitable degree of sensitivity and anatomical precision, and (2) to define a network of brain regions engaged by the motor task. We then explored the interaction between task and stimulation. This resulted in two contrasts (Main Effect of Movement – Left hand, Movement×Stimulation interaction – Left hand), and ensuing statistical parametric maps (SPMs).

All 10 subjects' normalised structural T1 scans (taken during ON) were combined to create a group structural T1 normalised to MNI space. One-sample t-tests were performed on group data separately for each of the contrasts to produce SPMs that were then superimposed on the group structural image. Second level tests on the main effect of movement contrast were adjusted for handedness and UPDRS-III ON score by including mean-centred cofounds in the second level design matrix (multiple linear regression model). This accounts for confounding effects due to inter-subject variability in Parkinsonian disability. In the same manner, second level tests on interaction contrasts were corrected for handedness and percentage improvement in UPDRS-III score when going from OFF to ON.

The main effect of movement contrast served to define the network of brain regions related to voluntary movement. The fifteen peak voxels of clusters larger than 5 voxels with the highest z statistics (range: 4.87–3.81) were defined as ‘nodes’ of the motor network. Any clusters consisting of purely white matter voxels were omitted.

A restricted volume analysis was then performed to assess the interaction between movement and stimulation within the nodes of the aforementioned network. Restricted volumes (8 mm radius spheres) were centred on the peak voxel of each node. A statistical threshold of p<0.05 (FWE corrected, with a cluster threshold of 5 voxels) was used to assess significance. Regions surviving this threshold were considered to show robust interactions between movement and stimulation.

Given that extensive PET imaging studies have previously proposed a network of areas involved in the therapeutic response to STN DBS and that this is the first report of fMRI activations and their modulation by therapeutic DBS, we also performed a whole-brain statistical search. From this analysis, regions that survived a criterion of p<0.0005 (uncorrected, cluster threshold of 5 voxels) are reported as discovered areas of interest that will be examined in future studies.

Effective connectivity analysis

Dynamic causal modelling (DCM) is a Bayesian framework that aims to explain how observed BOLD responses are generated by estimating the effective connectivity between specified regions of interest [26], [27]. DCM models hidden neuronal dynamics using an explicit forward model based on the balloon model [26], [28]. In brief, realistic models of the functional architecture are constructed, the BOLD signal from these regions is extracted and the neuronal (hidden) states are inferred. The effective connectivity between the regions is then expressed in the form of differential equations using three parameter sets; “average connection” parameters (values of the DCM A-matrix), represent latent or average coupling strengths in the absence of experimental manipulation (in our case, average connectivity represents the coupling during voluntary movement), “modulatory” or “bilinear” parameters (values of the DCM B-matrix) denote changes to the average connectivity associated with experimental manipulations (i.e. the additive effect of DBS on coupling strength), and thirdly, “input” parameters (values of the DCM C-matrix) control the effect of driving stimuli by external perturbations (in our case, movements). These parameters are then estimated using Bayesian estimators and are given in Hertz [29]. The coupling parameters represent changes in the sensitivity of one region to afferents from other regions, conceptually comparable to electronic gain; i.e. how much its output changes in response to a given input. DCM has become the method of choice for modelling effective connectivity in neuroimaging data and has been used widely across the literature [30]–[34]. Thus we modelled the effective connectivity between nodes of the motor network that demonstrated regional movement-related increases in BOLD signal that were sensitive to DBS.

The design matrix was finessed (rotated) for the DCM analysis. The left-hand movement ON and OFF scans were concatenated into a single (movement) regressor. Parametric modulators were used to model the movement×stimulation interaction. The main effect of DBS was modelled as a boxcar, with values of one during stimulation ON and zero otherwise.

Subject-specific peak coordinates of the regional interactions were used to identify nodes or regions in the DCM. The inclusion criteria for the DCM analysis required each subject to show a non-trivial interaction in both nodes (n = 7); within subject peaks (p<0.05, uncorrected) were within 16 mm of the second-level (between subject) peaks. Regional activity was summarised as the principal eigenvariate – adjusted for slow fluctuations and other nuisance variables – based on voxels within 4 mm of the subject-specific peaks.

Figure 1 summarises the different dynamic causal models we evaluated with Bayesian model comparison. Our movement effect entered all models as a driving input to the cortical node. All the areas had intrinsic (within region) and reciprocal extrinsic (between region) connectivity. The modulatory input was set to modulate a subset of connections in each model. These were thus; the forward and backward extrinsic connections between the two regions (identified in the prior restricted volume analysis at FWE corrected p<0.05), and the intrinsic self-connectivity within each region. This resulted in (2⁴ = ) 16 different DCMs per subject, and thus (16 DCMs×7 subjects  = ) 112 DCMs in total. Models were inverted and scored – in terms of their model evidence – using ‘DCM12’.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 1. The dynamic causal models compared using Bayesian model selection (model space).

Model 15 – the winning model – is shown enlarged. The blue node represents the right insula cortex, and red node, the right thalamus. Green arrows indicate the connection/s that DBS modulates. The ‘movement input’ is likely made up of both motor inputs arriving from M1, PM and SMA, as well as sensory inputs elicited by on-going movements. Thalamic ‘intrinsic subcortical projections’ refer to thalamic afferents from BG nuclei. Cortical ‘intrinsic insula cortical projections’ refer to cortical afferents from within the cortex. Average DCM parameters are included on the enlarged model 15, units are in Hertz (Hz). Positive A-matrix parameters represent an excitatory effect on the target, whereas negative values indicate an inhibition of the target area. Positive B-matrix parameters (value in green) represent an increased target response to input (i.e. an increased gain), whereas negative values indicate a decreased target response to input (i.e. a reduced gain). The coupling during movements with DBS is equal the sum of the A and B value on that connection, e.g. during movements with DBS, the cortico-thalamic drive switches from −0.21 Hz to (−0.21+0.30) = 0.09 Hz, i.e. it switches from an inhibitory to an excitatory drive.

https://doi.org/10.1371/journal.pone.0050270.g001

Bayesian model selection (BMS, fixed-effects assumptions) was then employed to select which of the 16 models had the greatest evidence, given the data collected [27], [35]. BMS computes a Bayes factor for each of the models, allowing us to make inferences about which of several biologically plausible models is optimal given the data [35]. Fixed-effects model comparison was chosen because we selected our subjects under the assumption they have the same functional architecture and that DBS had consistent effects within this anatomy. The free-energy of each model (F) corresponds to the log of the model evidence and indicates the accuracy of the model corrected for its complexity. To assess for the effects of outliers on the BMS, we plotted the relative F values for each model for each subject. Relative values were generated by subtracting the F value of the model with the least evidence, from each model in each subject. In addition, we performed a random-effects BMS for verification to allow for the possibility that different subjects had different connectivity architectures. For quantitative interpretation, the coupling parameters of the DCMs were averaged using Bayesian Model Averaging (BMA), in which parameter estimates are weighted by the model evidence [27].

Finally, we examined the relationship between connectivity parameters and the clinical UPDRS scores by performing correlation analyses between connectivity parameters during DBS-OFF and clinical score during OFF, connectivity parameters during DBS-ON and clinical score during ON, as well as percentage change in connectivity parameters and percentage change in clinical score.

Clinical response and motor task data

Clinical responses as measured by the UPDRS-III scores are shown in Table 1. The mean improvement was 27.5 points (56.7% improvement, p<10⁻⁶). Similar improvements were also seen when the task data was analysed. Left hand movement durations and reaction times were decreased in the ON condition by an average of 28.87% (p = 0.002), and 20.33% (p = 0.025) respectively. The mean movement duration during ON and OFF were 0.82 s and 1.27 s respectively. The mean reaction times during ON and OFF were 0.63 s and 0.83 s respectively. Post-operative MRI – employing fine cuts through the STN – confirmed that each electrode contact lay within or overlapped the anatomical border of the STN in both axial and coronal views.

Scanning proceeded with no adverse events or change in post-scan UPDRS-III scores. Re-introduction of PD medications led to restoration of baseline motor function. Post-scan inspection of the IPG revealed DBS stimulation parameters and circuit impedance were unchanged.

Hardware-related Artefact

All GE-EPI scans suffered dropout artefact thought to be caused by the subgaleal connectors between the leads and extension cables sited over the left parietal bone (see Figure 2a–c).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. Imaging results.

A typical drop-out artefact in a single subject's GE-EPI acquisition viewed from (a) axial, (b) coronal, and (c) sagittal sections; cross-hair position = −34.8, −21.5, 53.3 mm (MNI coordinates). SPMs in (d) summarize the movement network on a rendered MNI brain (p<0.001 uncorrected). Clusters representing BOLD signal increases in the insula cortex (e, green arrow), and thalamus (f, green arrow).

https://doi.org/10.1371/journal.pone.0050270.g002

Regional interactions between movement and stimulation.

The main effect of movement was in accordance with previously published accounts (Limousin et al., 1997). The purpose of this contrast was to establish a network of functionally specialised nodes associated with task performance in our cohort. A summary rendered SPM of the network is displayed in Figure 2d.

Adjustment for clinical response to stimulation did not affect the regions that demonstrated highest levels of peak level BOLD response, although it did increase the significance levels of most of the areas described.

Our initial restricted volume analysis of the interaction (movement×stimulation) revealed robust increases in BOLD responses in the right insula cortex, and right thalamus (p<0.05, FWE corrected) when DBS was active. See Figure 2d,e.

Subsequent whole brain analysis of the interaction (movement×stimulation), using uncorrected thresholds (p<0.0005) revealed additional increases in the left superior frontal gyrus (Premotor area, PM, BA 6) and middle frontal gyrus (BA 10/46, DLPFC), right intra-parietal sulcus (hIP1), and inferior frontal gyrus pars triangularis (BA 45). See Table 2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. The interaction between movement and stimulation identified using a whole-brain analysis.

https://doi.org/10.1371/journal.pone.0050270.t002

Dynamic Causal Modelling: Bayesian Model Selection

In order to explore the effective connectivity between the insula cortex and thalamus, we constructed 16 models of connectivity and used Bayesian Model Selection to determine the most likely model to produce our data. The relative log-evidences across all models for all participants are shown in Figure 3a,b.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. BMS results.

FFX = Fixed Effects Assumptions, RFX = Random Effects Assumptions. (a) The relative log-evidences across all 16 specified models with model 3 showing the highest log-evidence. (b) Given the observed data and the models specified, one can be >95% certain that model 15 is the data generator. (c, d) The relative log-evidences between the 4 most likely models – again highlighting model 15 as the most likely model. This is repeated using RFX BMS, confirming the FFX findings (e, f).

https://doi.org/10.1371/journal.pone.0050270.g003

Across all models, model 15 emerges as the most likely, followed by 3 and 1. The difference in relative log-evidences, ΔF, are 3.77 and 8.25 respectively, indicating that there is very strong evidence in favour of model 15 [37]. Restricting the BMS to the four most likely models illustrates this more clearly (3e,f). Relative F values for each model – for each subject – confirmed models 15, 3 and 1 consistently scored highly. Random effects BMS produced similar results, confirming model 15 as the most likely of all the models to be the generator of the data (3g,h).

Dynamic Causal Modelling: Parameter Estimates

Effective connectivity estimates (corresponding to the DCM A-matrix) reveal that during voluntary movements activity in the insula cortex had an inhibitory (cortico-thalamic) drive on the thalamus, whereas thalamic activity had an excitatory (thalamo-cortical) effect on the insula cortex. Note, in DCM, regions are often used to ‘stand in’ for multiple systems, and effective connectivity is polysynaptic. In our context, the insula region is taken as representative of a cortical system, while the thalamus summarises subcortical responses. The thalamo-cortical-cortico-thalamic loops most probably comprise several synaptic relays.

Crucially, DBS effectively reversed this extrinsic coupling, changing cortico-thalamic connectivity from inhibitory to excitatory, whereas the thalamus assumed an inhibitory effect on the activity in the insula cortex under DBS. In other words, DBS appears to sensitise subcortical responses to cortical afferents, with a reciprocal desensitisation of cortical responses to subcortical projections. This is a remarkable reversal of effective connectivity that is the quantitatively largest (to our knowledge) that has been reported in the DCM literature. There was no significant correlation between cortico-thalamic connectivity parameters and the clinical scores (P<0.05) to suggest a direct linear relationship between these indices.

Limitations of this study

The extension cables sited over left parietal bone created a drop-out susceptibility artefact on GE-EPI acquisitions, partially obscuring left sensorimotor regions. Given our careful preprocessing, including the use of field map correction and a bespoke normalisation procedure, there is no evidence to suggest that BOLD signal from the remaining brain regions were affected by this artefact, especially given the ‘main effect of movement’ contrast produced a network of well-described motor regions. Previous studies involving implanted electrodes have confirmed that similar artefacts do not significantly impair the functional data [58].

Only patients who had a significant therapeutic response to stimulation (minimum UPDRS-III improvement in our sample was 38.5%) were included in our study. While this permitted novel investigation of the modulatory effects of confirmed therapeutic stimulation, our conclusions should only be applied to patients who have shown such improvements.

We included a covariate for hand dominance to minimize variability vis-à-vis laterality of function. Other studies investigating motor control generally use right hand movements, complicating comparisons with other studies, but this was unavoidable given our standard surgical practice of placing the connector to electrode extension cables subcutaneously over the left parietal bone.

The T/R head coil was used to minimize RF exposure to the DBS instrumentation circuit in situ. Such coils are not regularly used in fMRI studies as they forgo the signal to noise ratio advantages of conventional multi-channel receive-only head coils. This may explain why regions previously reported in the literature and identified in our whole-brain search at uncorrected thresholds (including the premotor cortex, and DLPFC) did not survive FWE correction. While this was unavoidable given the safety concerns, we nevertheless identified two regions that are irrefutably associated with therapeutic STN stimulation, which formed the subsequent focus of the more sophisticated network modelling possible with fMRI data.

Our results cast the therapeutic effects of DBS in a fundamentally new light, emphasising a role in changing distributed cortico-subcortical interactions in a way that has not been previously explored. Investigating the effective connectivity changes induced by DBS in vivo represents a new avenue of study that may shed light on its underlying mechanisms of action. Given our modest sample size and small network, we stress that further work is required to verify and validate our findings. Selection of the target nucleus for DBS is pivotal to producing the desired therapeutic effect. Historically this has relied upon stimulating nuclei that had previously been targeted for ablative procedures. Understanding the impact that stimulation is having on the networks that course through the target may allow for improvements in current targeting, as well as rational selection of novel targets to extend its use to patients with other disabling conditions.