The data used in this manuscript has been published previously in Bächinger et al. 2019; experiment 6. The focus of the analysis in the previous publication was on changes in fMRI BOLD activation associated with motor slowing. Here, we re-analysed the data with an emphasis on network modelling with DCM. The behavioural task, fMRI preprocessing methods, and the general linear model (GLM) analysis of BOLD activity (i.e., parametric and block design-based analysis) is identical to Bächinger et al. 2019. We reiterate the relevant methods for the reader’s convenience.

Participants

Of the 25 right-handed participants who took part in the experiment, 24 were included in the DCM analysis. The mean (SD) age was 23.8 (3.3) years and 50% were female. One participant was excluded because they did not show any motor slowing, but rather an increase in tapping speed, indicating that they did not perform the task as instructed. All participants were free of medication, had no history of neurological or psychiatric disease and were naïve to the purpose of the experiment. All experimental protocols were approved by the research ethics committee of the canton of Zurich (KEK-ZH 2015 − 0537) and participants gave written informed consent to the study.

Behavioural task and analysis

The experiment consisted of two different conditions: finger tapping for either 30 s (slowing condition) or 10 s (control condition), each followed by a 30 s break. Tapping was performed alternating between index and middle finger at maximum speed. Participants were informed about the condition prior to the start of tapping with a visual get-ready cue (randomly jittered between 2 and 3 s). The conditions were blocked within each fMRI run: One block was made up of four trials of the slowing condition, followed by four trials of the control condition, or vice versa. Each participant performed two fMRI runs consisting of two blocks each. This resulted in 16 trials per condition. The starting condition of the first run (slowing or control condition) was alternated across participants and the second run had a counterbalanced order in relation to the first run. Additionally, an implicit baseline of 20 s was measured after each block (Fig. 1A). Behavioural data was analysed as described previously (Bächinger et al. 2019). Tapping and break intervals were divided into 10 s bins and movement speed was normalised to the average speed of the control condition per participant. This normalised movement speed was subjected to a linear mixed effects model with the fixed factor time (i.e., time bins) and the random factor participant. Motor slowing was defined as a significant main effect of time (Fig. 1B).

fMRI acquisition and preprocessing

fMRI scans were acquired with a Philips Ingenia 3T whole body scanner. Prior to the functional runs, high resolution T1-weighted anatomical scans were acquired (voxel size = 1 mm3, 160 sagittal slices, matrix size = 240 × 240, TR/TE = 8.3/3.9 ms). These anatomical scans were used for functional image registration and normalisation. During the behavioural runs 360 volumes were acquired in each run (voxel size = 2.75 × 2.75 × 3.3 mm3, matrix size = 128 × 128, TR/TE = 2500/35 ms, flip angle = 82 degrees, 40 slices acquired in interleaved order for full brain coverage). Preprocessing was performed using SPM12 (Wellcome Trust) with default parameters.

First, functional images were realigned to the average functional image. Then, the anatomical image was segmented, of which a transformation to MNI space was obtained, and the structural image was skull stripped. The functional images were co-registered to the anatomical image using normalised mutual information, and to MNI space through the forward transformation. The normalised images (2 × 2 × 2 mm3) were spatially smoothed with an 8 mm isotropic Gaussian kernel at full-width-half maximum.

fMRI data analysis

fMRI analyses were also performed in SPM12. The first-level model of each participant consisted of a general linear model. The GLM design matrix included four regressors of interest: tapping, parametric modulation of tapping, recovery, and parametric modulation of recovery. The tapping regressor represented the time periods when the participant was tapping. The recovery reflects the 30 s rest condition after a tapping trial. The parametric modulation regressor consisted of a linear increase over the tapping periods (reflecting the increase in motor slowing) or a linear increase over the recovery period after a 30 s tapping trial (but not a 10 s tapping trial). The linear increase was the same across all participants and did not depend on the participant’s performance. Importantly, the parametric modulation regressor was orthogonalized with respect to the tapping regressor. Note that the 30 s slowing condition and the 10 s control condition were modelled together in each regressor. For the parametric modulator, the slowing condition consisted of a linear increase in six bins of 5 s, and the control condition was made up of a linear increase in two bins of 5 s. Regressors of no interest in the GLM consisted of get-ready periods and six head movement parameters (translation and rotation along the x, y, and z-axis). All regressors except the six head movement parameters were convolved with a canonical hemodynamic response function. The two regressors of interest were contrasted against the implicit baseline and were then subjected to a second-level random-effects analysis across participants. The second level analysis was a single one-sample t-test contrasting the regressors of interest against zero. P-values smaller than 0.05 family-wise error (FWE) corrected for multiple comparisons were considered statistically significant. Localisation of functional clusters was aided by the anatomy toolbox (Eickhoff et al. 2005).

Fig. 1

(A) Experimental Design of fMRI study. 24 participants were either tapping for 30 s (slowing condition) or 10 s (control condition) during fMRI scanning. (B) Behavioural Results. Behavioural results showing a significant decrease in movement speed over 30 s of tapping. (C) fMRI activations, associated with either tapping itself (motor network, blue) or motor slowing (increasing activation with decreasing movement speed, green). Regions of interest where defined based on the closest individual activations of PMd, SMA, and SM1. (D) Schematic overview of the model space for DCM

Dynamic causal modelling

To investigate the changes in effective (directed) connectivity with motor slowing, we performed DCM (Friston et al. 2003) using SPM12. In short, dynamic causal models are generative models that aim to capture directed interactions among brain regions or states based on à priori hypotheses. In DCM, changes in brain states over time are modelled in the form of a state-space equation:

$$\:\dot=\left(A+_^_^\right)x+Cu\:\:\:\#\left(1\right)\:$$

In this equation, \(\:x\) is the state vector representing the current neuronal state and \(\:\dot\) refers to the change in the neuronal state over time. The matrix \(\:A\) represents the underlying endogenous or intrinsic connectivity with fixed weights defined by the model, whereas \(\:^\) reflects the weights of task-dependent modulations of connectivity, driven by external modulatory inputs \(\:_\). \(\:C\) represents the weights for direct inputs, characterising how the extrinsic driving input \(\:u\) directly influences brain regions.

In the scope of this study, the matrix \(\:A\) represents the endogenous connectivity of the motor system (i.e., bidirectional connectivity between SM1, PMd, and SMA; see also regions of interests). The term \(\:_^\) characterises the strength of the modulatory changes that occur due to motor slowing, which are modelled as a linear increase reflecting motor slowing. The final term \(\:Cu\), describes the external driving input to the (pre-)motor system, modelled here as a constant input from prefrontal areas to either PMd, SMA, or both.

Regions of interest, endogenous connectivity, and its modulation

As a hypothesis driven method, DCM requires a neurobiologically-plausible model of connectivity to be defined à priori. We therefore selected regions of interest which were associated with motor slowing. Specifically, we previously found that activity (pFWE < 0.05) of the left SM1, left PMd, and bilateral SMA (Bächinger et al. 2019) were inversely correlated with motor slowing of the right hand: All these regions showed an activation increase with decreasing tapping speed. Based on this finding, we investigated here whether motor slowing is associated with changes in premotor-motor interactions. To that end we built several DCMs incorporating PMC, SMA, and SM1 (i.e., the three areas directly associated with motor slowing).

We extracted the BOLD signal time-series of our Effect of Interest from 4 mm radius spheres centred on the following three regions of interest: SM1, PMd, and SMA. The Effect of Interest consisted of our four regressors of interest (tapping, parametric modulation of tapping, recovery, parametric modulation of recovery). All regions of interest were defined by taking the coordinates from the group-level analysis (Supplementary Material 1) and then extracting the closest peak-level activation on the single subject level.

The endogenous connectivity matrix (matrix \(\:A\) in Equation \(\:\left(1\right)\)) was defined by previous anatomical studies: Specifically, we assumed that all regions are connected bidirectionally based on previous anatomical findings (Luppino et al. 1993; Rouiller et al. 1994; Michely et al. 2015). Also, all included regions were assumed to be self-modulatory. Self-modulations were chosen to represent two-state models, which model two neuronal populations, an excitatory and an inhibitory, per region. In these two-state models, connections between two areas A and B are assumed to be excitatory and connection within each area A or B are modelled as being self-inhibitory (Marreiros et al. 2008; see Supplementary Material 2 for more detailed information). The model thus aims to explain changes in BOLD signal via the interplay of excitatory and inhibitory dynamics within and across regions, thereby incorporating current knowledge of the cortical microcircuit in which functional neuronal dynamics arise from an interplay between excitatory pyramidal cells (which are mainly glutamatergic) and inhibitory interneurons (which are mainly GABAergic, McColgan et al. 2020; Douglas and Martin 2004). As such, the two-state model is biologically more plausible and closer to the underlying anatomy of the human cortical microcircuit than single-state models (McColgan et al. 2020).

The extrinsic regressor (term u in Equation \(\:\left(1\right)\)) that modulates connectivity of the network (term \(\:B\) in Equation \(\:\left(1\right)\)) reflected the effect of motor slowing. In the scope of this study, motor slowing was simplified as a linear change over time, as represented by the parametric modulation (see section fMRI data analysis) which served as input for this analysis. The driving input to the model was assumed to be a constant input from prefrontal areas (Michely et al. 2015).

Model space and model families

As our main interest was to test whether self-inhibition of SM1 is crucial to explaining modulation of connectivity during motor slowing, we split the model space into model families with and without self-inhibition of SM1. Further, we wanted to investigate whether the premotor areas (SMA, PMd or both) shape the self-inhibition in SM1. Therefore, we set up multiple model families: (1) The top-down model family (in accordance with our main hypothesis outlined in the introduction), in which connections from SMA to SM1 and from PMd to SM1 were modulated. (2) The bottom-up model family, to verify whether our hypothesis may be inversed, meaning self-inhibition of SM1 may modulate SMA and PMd in a bottom-up fashion. In this model family, connections from SM1 to SMA and from SM1 to PMd were modulated. (3) The selective premotor model family, to test whether one premotor area is much more strongly involved in modulations of motor slowing: only SM1-SMA or SM1-PMd connections were modulated in these models. (4) The null model family. In this model family, none of the connections between any of the premotor areas and SM1 were modulated.

As mentioned, our main interest was to determine the necessity of self-inhibition in SM1, which is why these 4 model families were further specified as either having self-inhibition of SM1 modulated or not. This resulted in 8 model families: Top-down models without self-inhibition of SM1 (36 models, Fig. 2A), bottom-up models without self-inhibition of SM1 (36 models, Fig. 2B), selective premotor models without self-inhibition of SM1 (12 models, Fig. 2C), null models without self-inhibition of SM1 (6 models, Fig. 2D), top-down models with self-inhibition of SM1 (36 models, Fig. 2E), bottom-up models with self-inhibition of SM1 (36 models, Fig. 2F), selective premotor models with self-inhibition of SM1 (12 models, Fig. 2G), null models with self-inhibition of SM1 (6 models, Fig. 2H).

In all model families, the models were set up (i) with and without modulation of premotor interactions between SMA and PMd (Fig. 1D, blue), and (ii) with the driving input set to either PMd, SMA, or both (Fig. 1D, red). The main modulations of the selective premotor model family were either a bidirectional modulation of SM1-SMA (Fig. 1D, purple) or SM1-PMd (Fig. 1D, green). In both cases, the involved premotor area was modelled with self-inhibition (Fig. 1D, orange). The top-down and bottom-up model family differed in whether the connections SM1-PMd, SM1-SMA, or both were modulated, with or without self-inhibition of the involved premotor area. All in all, the model space consisted of 180 models, split into 8 model families. A list of all the models can be found in Supplementary Material 3.

Fig. 2

Model Families. One model family consisted of several models differing in which connections were modulated. The connections that may have been modulated within one model family are marked in red. The selective premotor models (C and G) either had modulations of the SMA-SM1 (purple) or the PMd-SM1 (green) interactions, but not both in the same model. The premotor interaction (SMA to PMd and PMd to SMA) were always modelled together. The driving input was set either to PMd, SMA, or both. Self-inhibition of SMA and PMd were only modulated in combination with a modulated connection with SM1. However, for model families A and E, the premotor-motor connections were also modulated without self-inhibition of the corresponding premotor area. Self-inhibition of SM1 (orange) was assumed in all models of model families E-H

Model selection and statistical analysis

To identify the most likely model family given the data, we used random-effects family-level Bayesian model selection (Penny et al. 2004; Stephan et al. 2009). As the model selection did not reveal decisive evidence for a single winning model family (family exceedance probabilities < 0.95, Fig. 3), all models across all model families were averaged through Bayesian model averaging (BMA) with an Occam’s window of 0.05 to inspect the model parameters. These BMA parameter estimates were then subjected to two further analyses (Stephan et al. 2010). First, to identify the connections which were significantly modulated by motor slowing across participants, a group level post-hoc analysis on the maximum-à-posteriori (MAP) of the matrix \(\:B\) was performed using Bonferroni-corrected t-tests. Secondly, a stepwise linear regression was performed to identify which of these modulated connections were directly associated with individual differences in motor slowing as quantified by the behavioural data. The regression model tested whether behavioural changes in tapping speed can be explained by the MAPs of the modulated connections.