Dynamic Effective Connectivity using Physiologically informed Dynamic Causal Model with Recurrent Units: A functional Magnetic Resonance Imaging simulation study

), this set can be represented as: gi = (see Figures 1B, C and Equation 8 below). The values of these parameters serve as the starting values or initializations for the next, i.e., (i + 1)th window. In other words, the parameters for window (i + 1) are initialized with the predicted values for its previous window i, preserving continuity between adjacent windows. It is to be noted that before performing model inversion, we do zero-padding with half-window length on each end of the data so that we can cover the entire extent of the actual signal. Therefore, the computed connectivity at every window is centered at that window.

The output of each dP-DCM recurrent unit is the fMRI BOLD response. A recurrence is being followed because the same unit with different parameter values repeatedly performs the same task or operation (on input sequences) till the final window.

Dynamic effective connectivity is estimated using overlapping windows. That is, partitioning the time-series (including the inputs and the observed/measured BOLD fMRI responses) into N number of overlapping windows such that ∀ i ∈ , we get:

d⁢P-D⁢C⁢Mi=(8)

{Xi(t+Δt)= WXiXi(t)+WUiUdi(t)Hi(t+Δt)= WHiθHi(Hi(t))+WH⁢XixEi(t)+ωHiyi(t+Δt)= WB⁢O⁢L⁢DiθB⁢O⁢L⁢Di([qi⁢(t+Δ⁢t)vi⁢(t+Δ⁢t)])+ ωB⁢O⁢L⁢Di

A schematic illustrating the updates of the neuronal, neurovascular, hemodynamic, and BOLD variables (following the above Equation 8) in each dPDCM recurrent unit has been shown in the Figure 1C. Our proposed workflow has also been demonstrated in Figure 2 in a nutshell.

Figure 2. Workflow of discretized Physiologically informed Dynamic Causal Model (dP-DCM)-recurrent unit (RU).

3. Simulations set-up and results

To check the face validity of the approach, first, we simulated connectivity profiles for 3- and 10-region graphical models (representing cognitive hypotheses), to show the ability of dP-DCM-RU to estimate dynamic effective connectivity and to distinguish different causal functional graphs using model evidence for a simple and a complex case, respectively.

3.1. Three-region model 3.1.1. Case a: Time-varying connectivity

This model comprises three regions (R1, R2, and R3) (as shown in Figure 3A). A sinusoidal input u is applied to R1 which then activates R2 and R3. Area-specific time-varying fMRI BOLD responses, given as a percentage signal change are shown in the Figure 3B. The colored boxes correspond to the recurrent window size for each of these time courses.

Figure 3. Three-region model and corresponding area-specific simulated functional magnetic resonance imaging (fMRI) blood oxygenation level-dependent (BOLD) responses. (A) Three-region model for which connections exist from Region 1 (R1) to Region 2 (R2) and from Region 1 (R1) to Region 3 (R3). The sinusoidal input u is applied to R1 and then activity propagates to both R2 and R3. The model inside the white circle (with black border) shows the internal sample model representation of R2. (B) Corresponding area-specific simulated fluctuating fMRI BOLD time courses with windows centered at the 30th second. The gray horizontal line represents the duration of the input stimulus.

3.1.1.1. Forward simulation

In this example, we have considered a fast-varying connection from R1 to R2 and a slow-varying connection from R1 to R3. The assumed connectivity time courses between R1 and R2 and between R1 and R3 are shown in Figure 4A (colored plots). Supplementary Table 1A shows the piece-wise continuous functions used to simulate the connectivity values. The effective connectivity as a function of time (t) is denoted as eff_conn(t). Using the simulated connectivity pattern, we get the corresponding area-specific fMRI BOLD responses as shown in the Figure 4B (colored plots). The value of Δt is 1/32s for all simulations (Wang et al., 2018). For all other parameters of the model, we please see Supplementary Table 1B (also please refer to Supplementary Table 1A of Havlicek et al., 2015).

Figure 4. Predicted estimates along with the ground truth values for Case a. The black dashed lines represent the predicted responses and the colored time courses are the ground truth values. (A) Connectivity time courses, (B) area-specific functional magnetic resonance imaging (fMRI) blood oxygenation level-dependent (BOLD) time courses each expressed as a percentage change in response. The gray horizontal bar represents the stimulus duration.

3.1.1.2. Model inversion

For this model inversion, we assume the same connectivity graph as was used for the forward simulation. In Figure 4, the black dashed lines represent the predicted estimates from the model on top of the colored ground truth (i.e., forward simulated) values. It can be noticed that for both connectivity time courses (Figure 4A) and fMRI BOLD responses (Figure 4B) the fitting is highly accurate, with slight inaccuracies in the effective connectivity estimates during the initial rise and during return to baseline, indicating that sharp increases or decreases are smoothed out in the BOLD signal and model inversion.

The Normalized Root Mean Squared Error [NRMSE (Shcherbakov et al., 2013)] value averaged over the two connections (R1 and R2, R1 and R3) is 2.14% and the NRMSE value averaged over the fMRI BOLD responses from the three regions (R1, R2, and R3) is 0.91%. We have repeated the above simulation set-up for a higher frequency input and have done corresponding model inversion, whose results have been provided in the Supplementary Section 1.1.1.

3.1.1.3. Model comparison

To show a noticeable difference in performance between two models during model comparison, we have considered an additional scenario. In this scenario, we have considered 2 competing hypotheses models m1 and m2 as shown in Figures 5A, B. In m1, driving input u1 (in red) is applied to R1, which is connected to both R2, and R3. Feedback connection exists from R2 to R1, which is influenced by modulatory input u2 (in purple). In m2, driving input u1 is applied to R1, which is connected to R2, which in turn is connected to R3. Feedback connection exists from R3 to R2. Connection from R2 to R3 is modulated by input u2. Hypothesis model m1 has been used for generating the ground truth BOLD data.

Figure 5. Competing hypotheses models (A) m1 (ground truth model) and (B) m2 (randomly chosen model) along with the corresponding area-specific functional magnetic resonance imaging (fMRI) blood oxygenation level-dependent (BOLD) responses. The black dashed lines represent the predicted responses and the colored time courses are the ground truth values. The red and purple arrows represent direct and modulatory inputs, respectively.

3.1.1.3.1. Predictions using m1 and m2

In Figures 5A, B, the black dashed lines represent the predicted estimates from the model on top of the colored ground truth (simulated) values.

For m1, it can be noticed that the fitting of the BOLD responses is accurate for all three regions (Figure 5A). The NRMSE value for this reconstructed fMRI BOLD time-series with respect to the ground truth time-series is 1.45%. For m2, the errors of the predicted BOLD responses are higher compared to those of m1 (Figure 5B). This can be attributed to the absence of feedback connection from R2 to R1 in model m2. The NRMSE value thus has increased to 10.31%. Hence, in terms of accuracy, m1 performed better than m2, ΔNRMSE (= NRMSEm2–NRMSEm1 = 8.86%) is high. This demonstrates that our method is able to adequately distinguish between two cognitive hypotheses.

In addition to the above m2, we have the also evaluated performances of 10 more three-region models with randomly chosen configurations (Supplementary Section 1.1.2). ΔNRMSE values of these models with respect to m1 (Supplementary Figure 2) clearly show that m1 is superior in terms of performance and m2 is one of the less inferior models in the group (m2–m12).

3.2. 10-region model

In a typical fMRI experiment, more than three brain areas are active. Thus, to demonstrate scalability, we evaluated a 10-region model as shown in the Figure 6 (center). The connectivity graph for the forward simulation is illustrated in the Supplementary Figure 4 (colored plots). Two time-varying inputs u1 and u2 (see Supplementary Figure 3) are applied to R1 and R2, respectively (see Figure 6 for the time courses). There is a delay of 20 s between these two inputs (see Supplementary Figure 3). Activity then propagates from R1 and R2 to the remaining regions.

Figure 6. Functional magnetic resonance imaging (fMRI) blood oxygenation level-dependent (BOLD) time courses for 10-region model where connections exist between R1 to R3, R3 to R8, R1 to R5, R1 to R6, R1 to R7, R2 to R6, R2 to R4, R4 to R9, and R4 to R10. Two time-varying inputs u1 and u2 are applied to R1 and R2, respectively. The black dashed lines represent the predicted responses, and the colored lines represent the simulated time courses.

3.2.1. Forward simulation

The assumed connectivity time courses are shown in the Supplementary Figure 4 (colored plots). Using the model and the connectivity time courses, we get the corresponding area-specific fMRI BOLD responses.

3.2.2. Model inversion

In both Figure 6 and Supplementary Figure 4, the black dashed lines represent the predicted estimates from the model on top of the colored ground truth (simulated) values. The prediction has a low NRMSE value (averaged over all the 10 regions) of 1.17%. The predicted connectivity time courses are also shown in Supplementary Figure 4. As can be seen, the predictions follow the ground truth time courses very closely for all brain areas.

3.2.3. Model comparison

For model comparison purposes, we additionally performed model inversion for randomly selected model m2 as shown in Figure 7 (center).

Figure 7. Predicted area-specific functional magnetic resonance imaging (fMRI) blood oxygenation level-dependent (BOLD) time courses along with the ground truth values for m2 which is a randomly chosen model. The black dashed lines represent the predicted responses and the colored time courses are the ground truth values. There are fitting inconsistencies for R3, R4, R5, R6, R7, R8, R9, and R10 fMRI BOLD time courses.

3.2.3.1. Prediction using m2

In Figure 7, the black dashed lines represent the predicted estimates from the model on top of the colored ground truth (simulated) values. Please notice that for fMRI BOLD responses the reconstruction and hence the fitting is good in the earlier brain areas, such as R1 and R2, but the discrepancies become more pronounced and errors larger for the later brain areas. The predicted connectivity time-courses are shown in Supplementary Figure 5. The NRMSE value (averaged over all the regions) for this reconstructed fMRI BOLD time-series with respect to the ground truth time-series is 25.75%. Hence, in terms of accuracy, m1 performed better than m2, ΔNRMSE (= NRMSEm2–NRMSEm1 = 24.58%) is high showing that our method can sufficiently differentiate between two cognitive hypotheses. Additionally, we have also considered another 10-region model with reciprocal connections and have performed corresponding model inversion with a randomly chosen configuration as shown in the Supplementary Section 3, confirming that m1 can be distinguished from a model with erroneous connectivity graphs.

4. fMRI BOLD time courses with measurement noise

In this case, we have used two models- a Three-region model (from “section 3.1.1 Case a: Time-varying connectivity”) as a simple model and a 10-region model (from “section 3.2.1 Forward simulation”) as a complex model. Here, the fMRI BOLD time courses are noisy with measurement or physical noise being present. Measurement or physical noise is an external noise which gets added to the signal while it is being acquired. For simulation purposes, this measurement noise used by us is a zero mean gaussian random noise with varying levels of standard deviation. Contrast-to-Noise Ratio (CNR) values were computed with respect to the ground truth values for fMRI BOLD responses. It is worth noting that we have defined CNR as the ratio between the standard deviation of the signal and the standard deviation of the noise (Definition four of CNR, Welvaert and Rosseel, 2013) and this definition is consistent with the previous DCM works (Friston et al., 2003; Frässle et al., 2017).

The extracted fMRI time series stem from Principal Component Analysis (PCA) over numerous voxels in local volumes of interest which suppresses noise (Frässle et al., 2017). Therefore, as outlined in Frässle et al. (2017), the typical CNR levels (SNR in DCM terminology, see footnote 3 for further clarification) of fMRI time series used for DCM are three or more. It is to be noted that this value is highly dependent on the tasks and brain areas investigated. That is, for some tasks and/or brain areas, DCM and its variants and the derived connectivity values are unreliable. This is already true for DCM estimating only one connectivity profile for a run and even more so for estimation of time-varying connectivity. Please note that this is also true for estimation of time-varying functional connectivity using resting-state. That is, it is not our claim that our approach (and no other statistical approach) will work for any task and/or brain area or any MRI acquisition parameters (such as field strength, sequence, etc.) but only if certain conditions, such as CNR levels, are met that time-varying effective connectivity can be estimated using our approach. Therefore, we have conducted the simulations for 3- and 10-region models with different CNR values [CNR = (1, 3, 5, 10, 20)] and repetition times [TR = (4, 2, 1, 0.25, 0.1 s)]. In each case, we resampled (via linear interpolation) the timeseries to sampling frequency of 32 Hz resulting in an increase in the number of samples. Using our method, we did model inversions for each of these settings for both three-region and 10-region models. We compared the estimated effective connectivity time courses to those of the simulated ground truth using Normalized Root Mean Squared Error (NRMSE), expressed as a percentage (%) as shown in Figure 8. We repeated the above setup for 10 runs with newly sampled noise and reported the mean and standard deviations (s.d.) of the NRMSE values over these 10 runs in Figure 8.

Figure 8. Normalized Root Mean Squared Error (NRMSE) (mean ± s.d.) expressed as a % vs. CNR for. (A) Three-region model and (B) 10-region model for five different repetition time (TR) settings.

We make the following two major observations for both three-region (Figure 8A) and 10-region (Figure 8B) models:

(i) The respective NRMSE values decrease with increase in the CNR levels for each TR setting. Moreover, at any TR value we observe that the standard deviations decrease with an increase in the CNR level, as expected. This indicates that the reliability of predictions increases toward high CNR values (i.e., when the data is less noisy).

(ii) The corresponding NRMSE values decrease with a decrease in the TR value for each CNR level. Notably, at any CNR level, the standard deviations also typically decrease with a decrease in the TR value (i.e., when the sampling rate is high) indicating more reliable predictions.

5. Discussion

The fMRI signal is an indirect reflection of neuronal activity, mediated by neurovascular coupling and hemodynamics. Generative models describe the biophysical basis underlying fMRI and present a framework to interpret empirical observations. Through model inversion, generative models enable investigations of underlying neuronal dynamics and functional integration in the brain. One such state-of-the-art generative model is the Physiologically informed Dynamic Causal Modeling (P-DCM, Havlicek et al., 2015). Most existing DCM studies (Friston et al., 2003; Stephan et al., 2008; Moran et al., 2009; Havlicek et al., 2015) typically consider the effective connectivity to be static for a cognitive task within an experimental run. However, experimental conditions can vary with time, especially in cases of complex stimuli, e.g., movie, music, etc. Consequently, the connectivity strengths between disparate brain regions involved in processing these complex stimuli may fluctuate with time. Please note that in the conventional DCM framework (Friston et al., 2003; Stephan et al., 2008; Moran et al., 2009), it may also be possible to model dynamic connectivity by utilizing dynamic B or C matrices (please refer to Equations 2a–f for definition of B and C). However, such an approach requires prior knowledge of the time-varying connectivity profiles and just estimates the strength of these dynamic connectivity predictors. In contrast, our method does not require prior knowledge of connectivity profiles.

In the recent years, there has been an increasing number of studies to elucidate the dynamic (functional) connectivity in fMRI by investigating the temporal correlations of resting-state BOLD fluctuations in distributed brain areas (Cribben et al., 2012; Handwerker et al., 2012; Calhoun et al., 2014; Monti et al., 2014). In such Dynamic Functional Connectivity (DFC) studies, one of the predominant methods is to employ a sliding window-based approach to find the time-varying correlations (Chang and Glover, 2010; Kiviniemi et al., 2011; Jones et al., 2012). However, lacking a generative model, the correlations between the areas are determined on the level of observations but not on the level of the underlying causes (Stephan et al., 2010). In contrast, DCM accounts for the indirect nature of the BOLD signal and fits BOLD signals in the different ROIs using a system of differential equations (Friston et al., 2003; Stephan et al., 2008; Havlicek et al., 2015, Havlicek et al., 2017b).

In this paper, we have introduced discretized Physiologically Informed Dynamic Causal Model with Recurrent Units (dP-DCM-RU) to characterize dynamic effective connectivity of various brain regions during tasks. This method is a combination of two approaches, namely, a Euler based discretization technique and a recurrent sliding window approach for dynamically modeling fMRI BOLD responses and for exploring the causal interactions between different neuronal populations. To validate, we have carried out simulations with 3- and 10-region models. To that end, we have decomposed effective connectivity into static and dynamic components. The static component acts as a baseline component and the dynamic component varies with time sinusoidally. However, please note that our recurrent window-based parameter estimation method can predict any connectivity profiles.

For the Three-region model, we have considered two different connectivity graphs. Using the first example, we have simulated noiseless time-varying effective connectivity between the regions. The results show that the fits of fMRI BOLD responses and the effective connectivity have low error values, compared to the ground truth (i.e., forward simulated) BOLD responses. This is not a trivial result as even assuming the correct connectivity graphs does not guarantee model invertibility due to potential ill-posed problem.

To show a noticeable difference in performance between models during model comparison, we have selected another example scenario (Three-region model) with feedback from R2 to R1 modulated by an input (i.e., the ground truth model, m1) and 11 randomly chosen models for comparison. Results show that the randomly chosen models did not perform well in terms of fitting accuracy values (given by NRMSE) relative to the ground truth model (m1). The results indicate a clear distinction between the hypothesis models and demonstrates that the dP-DCM-RU approach does not necessarily guarantee a good fit to any BOLD signal and is therefore capable of distinguishing models.

Typically, more than three brain regions are active during a complex cognitive task. Therefore, to illustrate a more complex case, we have performed simulations for a 10-region model (also see 10-region model example with reciprocal connections in the Supplementary Section 3). Results clearly suggested that the method was able to predict effective connectivity with low error values and to fit fMRI BOLD responses with high accuracy. We did model a competing 10-region connectivity graph (hypothesis model) for comparison. Model inversion results demonstrated that this randomly chosen hypothesis model was unable to reconstruct fMRI BOLD signals accurately leading to large deviation values from the ground truth fMRI BOLD responses, in particular, in areas most distant from the input areas. In this 10-region example (“section 3.2 10-region model”), we did not provide any delay between the input and connectivity (e.g., R1 to other regions). However, the window/duration of connectivity and input do not necessarily have to match. That is, connectivity from a certain brain region to another brain region may change from baseline value later than the input depending on the experiment (for such an example, please refer to “section 3.1.1 Case a: Time-varying connectivity”).

In addition, we have considered noisy scenarios by adding measurement noise (see “section 4 fMRI bold time courses with measurement noise”). We have conducted simulations for 3- and 10-region models with different Contrast-to-Noise Ratio (CNR) values [CNR = (1, 3, 5, 10, 20)] and repetition times [TR = (4, 2, 1, 0.25, 0.1 s)]. We observed that the connectivity prediction error decreases and the reliability of predictions increases with an increase in CNR and a decrease in TR values, respectively (see Figure 8). Depending on the threshold for accuracy used, estimation of dynamic connectivity using our method may be prohibited (the same argument also applies for dynamic functional connectivity calculation). In the cases of low CNR and high TR values, standard denoising procedures may be used, which may lead to less erroneous and more reliable parameter estimates.

One of the foremost limitations of DCM in general is that they are restricted to a fixed number of regions because of their computational demand (Frässle et al., 2017). In particular, the number of coupling parameters (i.e., connectivities) grows quadratically with the number of nodes, and therefore the computational time required to invert these models grows exponentially with the number of free parameters (Seghier and Friston, 2013; Frässle et al., 2017). Since our method is based on the P-DCM framework, the above problem also persists in our case. Additionally, the windowing scheme makes the approach even more computationally intensive. Furthermore, higher-order integrators are potentially slow in practice and computational (memory) requirements become even larger because of explicit Jacobian-based update schemes, which are evaluated numerically at each time point (Friston et al., 2003). Therefore, to increase computational speed and reduce memory, we utilize the lower-order Euler method. However, a disadvantage of the Euler method compared to higher order ODE solvers is lower numerical accuracy. Nevertheless, the numerical errors can be kept low for small step size Δt values. Following Wang et al. (2018), we have assessed the impact of different step sizes for the three-region model. As illustrated in the Supplementary Figure 10, NRMSE is the lowest (with a value of 0.82%) for Δt = 1/32 s and it increases (more than linearly) with step-size (e.g., for Δt = 1/4 s, NRMSE = 19.97%). Hence, we selected Δt = 1/32 s for all our simulations.

Another challenge for dP-DCM-RU is the selection of optimal window size. If the window sizes are too large, then the transients may not be captured, whereas too small window sizes may lead to overfitting of the model (Sakoğlu et al., 2010; Shakil et al., 2016). We investigated the effect of window sizes on the (three-region) model performance given in terms of NRMSE (%). Supplementary Figure 11 shows that for a window size of 5 s, NRMSE is low (1.45%), whereas model performance substantially degrades for a window size of 11 s (NRMSE increases to 22.68%). However, the window size can be easily adjusted to the complexity of the experimental design. A 5 s window size worked well for varying connectivity chosen in our simulations (Supplementary Figure 11). Please note that—as the BOLD signal can be represented as a smoothing kernel—very fast neuronal dynamics cannot be recovered by fMRI for any window size, independent of acquisition speed and analysis method (Polimeni and Lewis, 2021).

Park et al. (2018), Zarghami and Friston (2020) have proposed dynamic functional connectivity estimation frameworks for resting state fMRI (rs-fMRI). Park et al. (2018) utilized discrete cosine transform eigenvariates and Hierarchical Parametric Empirical Bayes (PEB) approach (Friston et al., 2016) to model dynamic functional connectivity at two levels. In the first level, they have inverted a spectral DCM (spDCM) separately for each window to obtain (within-window) connectivity parameters. Subsequently, in the second level they have applied PEB to estimate between-window effects on these connectivity parameters. Motivated by the nonlinear dynamical systems theory, Zarghami and Friston (2020) proposed a hybrid generative framework consisting of Hidden Markov Model (HMM), PEB and spDCM (discrete and continuous hierarchical models) to explain metastable dynamics in the brain via modeling the temporal evolution of different connectivity states. Their paradigm utilized the variational message passing technique, for which the HMM provided Bayesian model averages for the intermediate PEB level, which successively supplied priors to each spDCM. In this manner, by assigning an itinerant prior to the state-transition matrix, they estimated the transition and state-dependent effective connectivity parameters. On the contrary, our approach (“section 2.2 Time varying calculations of the model parameters”) neither requires “two-level” connectivity estimation nor involves a “hybrid” generative approach. Furthermore, these previous studies have typically considered larger window sizes for their estimations, which they have shown to work well for rs-fMRI. It is important to note that in contrast to our approach, these works do not capture the entire temporal extent of connectivity dynamics (i.e., connectivity values do not exist for all time-points) and are only suitable for tracking slow dynamics (Park et al., 2018, Zarghami and Friston, 2020). Furthermore, our method considers the stride of the overlapping windows to be 1 s and, therefore, covers the entire extent of the signal (in the order of seconds), and is able to track relatively faster dynamics (as demonstrated in Figure 4, where we have considered slow and fast varying connectivities). Finally, our work deals with task-based fMRI and not rs-fMRI.

Another recent work (Wang et al., 2018) utilizing Recurrent Neural Networks (RNNs) (Rumelhart et al., 1985; Hochreiter and Schmidhuber, 1997) proposed a biophysically interpretable DCM-RNN. Although both dP-DCM-RU and DCM-RNN take inspiration from the recurrence concept as in RNNs, DCM-RNN differs from dP-DCM-RU in many aspects: In their method, they have used Truncated Backpropagation Through Time (TBPTT) for computing parameter updates, whereas we have simply used gradient descent as done in standard DCM implementations (Havlicek et al., 2015). Their definition of recurrence is the same as in standard RNNs utilizing segmented batches. Therefore, they have used multiple batches in parallel (as in deep learning) and updated model parameters via TBPTT. To ensure that the gradients obtained by TBPTT are reliable, each of these segments has to be sufficiently long and the sampling time has to be sufficiently small. Since they use parallel batches for parameter update, the gradients may often not be accurate (Wang et al., 2018). This is because each of those batches does not represent the characteristic of the whole signal. Due to this batch parallel processing, they can only estimate static connectivity. In our case, the recurrence lies between the successive P-DCM units. We do sequential processing of each such unit estimating time-varying connectivity parameters. Furthermore, unlike DCM-RNN, we have used the state-of-the-art DCM model, i.e., P-DCM (Havlicek et al., 2015) in our framework instead of Single-State DCM (S-DCM) (Friston et al., 2003). Finally, the authors of DCM-RNN claimed that their model can be extended for complex paradigms such as movie watching using representations of the complex stimuli. However, at this stage their model cannot estimate dynamic connectivity without further modifications.

One of the major DCM steps is to conduct a Model Selection (using group Bayes factor) between several alternative competing models to establish which model accounts best for the experimental observations. After selection of the optimal model, making further inferences about its parameter estimates (e.g., connectivity) is typically not done in DCM studies and usually some statistical values (e.g., mean) of the parameter estimates for the group is reported (Stephan et al., 2009). However, inferences about model parameters can still be made with either a fixed effects or a random effects approach. For fixed effects parameter inference, a common way is to use Bayesian average (for example “DCM average” function in SPM) (Friston et al., 1994). For random effects parameter inference, subject-specific parameter estimates can be used with a classical frequentist test, such as a paired t-test (between model parameters) or repeated measures ANOVA in case of multiple sessions per subject (Stephan et al., 2009).

It is important to note that our proposed method is not necessarily restricted to block designs. Naturalistic stimulus typically comprises block (e.g., fluctuation of light intensity) and event-related (e.g., presence of face for a limited period of time) components. However, that due to sluggishness of the hemodynamic response, very fast events may not be detected, similarly as in standard fMRI data and analysis. Please note that our approach can be used with arbitrary time-varying inputs and the choice of sinusoidal inputs (in our simulations) is for illustrative purposes and can easily be modified.

In our approach, the estimates from the previous window serve as the initial values (and not constraints) for the next window. Therefore, the model gets a good starting point which makes it easier to optimize. However, for any window, model inversion is performed independently, therefore, an estimation error in the previous window will unlikely be reflected in the next window. Furthermore, for optimization we have resorted to using the conventional gradient based scheme, i.e., gradient descent. Although gradient descent can potentially be slow, it has a better generalization performance (Zhou et al., 2020). Alternatively, one can use momentum based (Sutskever et al., 2013) or adaptive algorithms (Kingma and Ba, 2015) to further speed up performance while maintaining accuracy (Ruder, 2016). Nonetheless, we have not explored different optimization algorithms and strongly feel that it is beyond the scope of the current paper since this is a general topic for all DCM and model inversion approaches and not specific to our paper. To summarize, when a subject is exposed to a complex stimulus (e.g., watching a movie), human brains show dynamic effective connectivity between remote areas on the neuronal level, which can be indirectly measured using fMRI and which can be effectively recovered using the d-PDCM-RU approach. In the future, we will demonstrate the validity of our method in clinical and cognitive neuroscience studies.

Data availability statement

The original contributions presented in this study are included in the article/Supplementary material, further inquiries can be directed to the corresponding authors.

Author contributions

SN and KU made substantial contributions to the conceptualization, methodology, design of the experiments, data analysis, visualization, and drafting the manuscript. Both authors provided final approval of the submitted version of the manuscript.

Funding

This study was supported by the Institute for Basic Science, Suwon, Republic of Korea (IBS-R015-D1) to KU.

Acknowledgments

We thank Labeeb Talukder and Shawn Carere for their constructive feedbacks on the manuscript.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnhum.2023.1001848/full#supplementary-material

Footnotes References

Attwell, D., Buchan, A. M., Charpak, S., Lauritzen, M., MacVicar, B. A., and Newman, E. A. (2010). Glial and neuronal control of brain blood flow. Nature 468, 232–243.

Google Scholar

Buxton, R. B. (2012). Dynamic models of BOLD contrast. Neuroimage 62, 953–961.

Google Scholar

Buxton, R. B., and Frank, L. R. (1997). A model for the coupling between cerebral blood flow and oxygen metabolism during neural stimulation. J. Cereb. Blood Flow Metab. 17, 64–72. doi: 10.1097/00004647-199701000-00009

PubMed Abstract | CrossRef Full Text | Google Scholar

Buxton, R. B., Uludağ, K., Dubowitz, D. J., and Liu, T. T. (2004). Modeling the hemodynamic response to brain activation. Neuroimage 23, S220–S233.

Google Scholar

Calhoun, V. D., Miller, R., Pearlson, G., and Adalı, T. (2014). The chronnectome: Time-varying connectivity networks as the next frontier in fMRI data discovery. Neuron 84, 262–274. doi: 10.1016/j.neuron.2014.10.015