For the computation of PAC we followed the parametric approach described in Penny et al. (2008) and van Wijk et al. (2015). Specifically, we first obtained low- and high-frequency components of the LFP signals by bandpass-filtering around centre frequencies between 5 and 35 Hz with 0.5 Hz steps (filter bandwidth ±1 Hz), and between 150 and 400 Hz in steps of 2 Hz (filter bandwidth ±35 Hz). Subsequently we extracted the instantaneous phase for each low-frequency component via θx=modanglehilbertx,2π, and amplitude of each high-frequency component via ay=abshilberty. To avoid filter ringing we removed the first and last 167 ms of each epoch. PAC was then estimated for each low-/high-frequency pair of frequency combinations using a general linear model of the form:ay=β1sin(θx)+β2cos(θx),where we took r=β12+β22 as the strength of PAC. We performed this analysis in two ways for each individual STN channel: (1) with time series concatenated across all epochs to assess the overall PAC; (2) for each epoch individually. The latter results in a set of β-coefficients, which we used to test for significance of the overall PAC via an
