A set of functionally and structurally connected brain areas that are activated in a resting state and deactivated by external stimulation or cognitively effortful tasks have been identified as the default mode network (DMN) in several early papers (Gusnard et al., 2001; Gusnard and Raichle, 2001; Raichle et al., 2001). A number of hypotheses on the functionality of the DMN have been formulated mostly relating it to higher-order aspects of consciousness and cognition (Gusnard and Raichle, 2001; Gusnard et al., 2001; Schilbach et al., 2008).

The main part of the DMN has been identified in the medial prefrontal cortex (MPFC), posterior cingulate cortex (PCC), and the inferior parietal cortex (IPC) of both hemispheres (Buckner et al., 2008). Key DMN regions are well-connected structurally (Vogt et al., 2006; van den Heuvel et al., 2008, 2009; Greicius et al., 2009) and functionally (Greicius et al., 2003; Biswal et al., 2010). The MPFC and PCC are linked by powerful cingulate paths, which are identified by diffusion tensor imaging (DTI) (van den Heuvel et al., 2008). IPC is a functionally heterogeneous region involved in visual-spatial orientation, attention, memory, and math knowledge (Uddin et al., 2010). Paths connecting IPC and PPC were identified using a probabilistic approach of DTI (Khalsa et al., 2014; Tao et al., 2015). With respect to hippocampal formation (HF), rank correlations of activity also reveal the basic pattern of activation/deactivation characteristic of the DMN (Greicius et al., 2004; Vincent et al., 2006). Functional connectivity (FC) of both hippocampi has been analyzed in studies (Cui et al., 2015; Robinson et al., 2015) e.g., in a recent meta-analytic study by Robinson et al. (2016). However, causal relations of HF to other brain structures in the resting state have been investigated far less systematically than interactions between the key DMN regions.

An appropriate method to study these interactions is dynamic causal modeling (DCM), a Bayesian approach typically used to explain effective connectivity changes underlying task-related brain responses (Friston et al., 2003). The DCM estimates effective connectivity, which is the measure that mediates the influence one neuronal system exerts on another. DCM treats the brain as a ‘black box’ which receives input and generates output (Razi and Friston, 2016). To model the resting-state, a new version of DCM was introduced based on a deterministic model that generates predicted cross spectra, the spectral DCM (Friston et al., 2014). The resting-state fMRI signals convey fluctuations in the low-frequency band typically within 0.01–0.08 Hz (Biswal et al., 1995). To identify active regions of the DMN and corresponding time-series, the resting state is modeled using a generalized linear model containing a discrete cosine basis set with frequencies in the corresponding range, in addition to the individual nuisance regressor (Fransson, 2005; Kahan et al., 2014). After this, DCM models are specified in terms of nodes and possible connections between them. For example, in the present study, we consider two model sets with LHIP/RHIP inclusion together with the basic four DMN nodes, and one set with LHIP inclusion and PCC exclusion from the basic nodes. The best DCM model at the group level is then determined using the Bayesian model selection (BMS) procedure (Penny et al., 2004; Stephan et al., 2009).

The first work on the dynamic modeling of the causal links of DMN revealed that MPFC exerts greater effect on PCC than vice-versa, while bilateral IPC transmit information in the PCC and the MPFC (Di and Biswal, 2014). In addition, the authors found that the endogenous exposure may be greater in the right hemisphere than in the left. The introduction of full connectivity, taking into account inter-hemispheric connections between IPC, partially destroys the symmetry of these components but retains the main trend – the predominance of MPFC impact on PPC (Bastos-Leite et al., 2015). In a recent DCM analysis of causal relations between the HF, MPFC, the anterior cingulate cortex (ACC) and dorsolateral prefrontal cortex (DLPFC) Cui et al. (2015) showed slight asymmetry with a predominance of effects from MPFC to the right and left parts of HF, at least in healthy subjects. At the same time, the left MPFC seems to exert a stronger influence which is more apparent on the same side.

The present study has two related objectives. Firstly, it aims at localizing the hippocampal regions within the known structure of DMN causal connections, i.e., to examine effective connectivity between the four core DMN nodes and left and right parahippocampal regions estimating the coupling parameters. The second objective is to evaluate the model stability when adding or excluding regions of interest. For example, does the established connectivity pattern between sources change by addition of a new source, such as the left hippocampus? If not, then it could be possible to increase model complexity beginning with the model for well-known DMN sources by adding other candidate structures to participate in coupling. When estimating a more complex model, one could be confident that the previous connectivity pattern does not change much and all changes in the extended model are associated with the new source (region of interest). In our previous work (Sharaev et al., 2016), we found a stable coupling pattern between the four basic DMN nodes. Here, we aimed at developing the model to assess its parameter stability when adding or removing active regions. As in the previous work, endogenous fluctuations were explicitly modeled using Discrete Cosine Set and spectral DCM was applied to find effective connectivity patterns. To test only biologically plausible hypotheses on effective connectivity and to reduce the model space, we used structural connectivity data described by Tao et al. (2015).

Firstly, it should be noted that after performing conventional SPM analysis we found that 4 of 30 subjects did not reveal significant activity in either the LHIP or RHIP region. Obviously, such cases need to be more deeply investigated, but for the current study, the data from these subjects were excluded from the modeling step and further consideration.

For extended models with five sources (including LHIP or RHIP), after inverting 96 DCMs for 26 subjects, we received 96 × 26 = 2496 F-values (the log-evidence approximation for each model for every subject) and for the reduced model (with LHIP but without PCC), after inverting 72 DCMs, we received 72 × 26 = 1872 F-values. With a large number of models (e.g., 96 or 72), a question arises: do these models behave alike across subjects? If they are stable, i.e., the same model behaves in a similar way when applied to different subject data, then one can expect that the model reflects some factual neural processes. Otherwise, when the model performs randomly across subjects, it probably does not describe the same underlying neural activity. To answer this question, we counted correlations between individual F-values for 96 (in the case of LHIP/RHIP) and 72 (in the case of the reduced model without PCC) models across all 26 subjects. This results in correlation matrices with 26 rows as shown in Figure 4A. The color encodes the pairwise correlation value. The posterior probabilities of model families are shown in Figure 4B, and the sums of the models’ F-values across subjects for the winning family a is shown in Figure 4C.

As can be seen from the matrices, for most subject pairs, the correlation is rather high (mean value about 0.7–0.8), except for a couple of subjects for whom correlation was somewhat less. This is true for all models sets. Thus, we can conclude that models are quite stable across the group, since the same model behaves in a similar way when applied to different subject’s data, producing highly correlated F-values. Because there are no negative values in correlation matrices, this means that no models perform in the opposite way across subjects.

The winning families are a and 6 for LHIP inclusion, a and 1 for RHIP inclusion (Figure 4B). Regarding family a, one may recall from Figure 2 it is the full connected base, which was the best model when analyzing four source models (Sharaev et al., 2016). This means that no matter how the LHIP/RHIP region is included, the best connection pattern between these four nodes remains the same. This is a significant finding, because it means that connectivity between four basic DMN nodes is not corrupted by adding the fifth node. Next, the 12 best performing models from family a are shown as 12 peaks in Figure 4C.

From Figure 4B (family a) and Figure 4C, it is clear that five models (a_5, a_9, a_6, a_10, a_1) are better than others, both for the LHIP and RHIP inclusion scheme. Though other models perform significantly worse and can be easily discarded, it becomes hard to distinguish between these five leading models. The same situation remains if we consider the number of wins, i.e., how often each model was the best one among 96 competing models in the group. The results are provided in Table 1 below:

In both groups, the model a_1 (full connected base and full connected LHIP/RHIP areas) wins by a narrow margin, though by the BMS results, this model is the best one only in the RHIP group; in the LHIP group, the best model is a_6. All five models from Table 1 imply that both hippocampal regions have causal connections to the majority of four key DMN regions. They also have a very similar architecture, which could explain why it becomes so difficult to find the best one among them. The BMA/BPA results are shown in Tables 2 and 3 for LHIP/RHIP inclusion, respectively. The BMA results from the previous four source analysis are shown in Table 4.

The tables above show that in the base nodes set, only the connection from mPFC to LIPC significantly changed in the case of LHIP inclusion. This means that adding a new source (LHIP or RHIP) practically did not change the basic connectivity pattern. The best models with their BMA results are shown in Figure 5 (note that BMA are plotted only for links from/to LHIP/RHIP).

Bayesian model selection on reduced models with LHIP but without PCC showed that the best model at the group level is a′_6′ – the same connectivity pattern as for the five source model including both LHIP and PCC. The number of wins by top performing models is: for a′_6′ and a′_1′ – both 5 wins; a′_5′ and d′_1′ – both 3 wins. There is a parity between a′_1′ and a′_6′ models, also a′_5′ and d′_1′ models are rather close to the winning ones. It is worth noting that a model with d′ base set outperforms many models with a full connected base set (a′). The most interesting part of the reduced models analysis is the BMA result in Table 5, comparing it to the five sources in LHIP BMA. Again, it can be seen that the majority of parameter estimates for the model with five sources and with four sources are very close, no new connections arose, and no connections disappeared. The significant difference (p < 0.05, insignificant after Bonferroni correction) between the parameters was found for only one connection: from LHIP to LIPC (0.28 versus 0.18, both positive – excitatory connections).

Overall, the important finding is that, with spectral DCM, there are stable patterns of effective connectivity even if a key region of the explored network is missed. The best reduced model at the group level is shown in Figure 6A.

A comparison of effective and functional connections is important as well. All significant functional connections are presented in Table 6; Figure 6B. The six regions of interest are interconnected with positive connections (correlations). The numbers are Fisher-transformed correlation coefficients. The model with four basic ROIs taken from (Sharaev et al., 2016) is shown in Figure 6C.

In order to assess the reproducibility of DCM parameter estimates, we examined the distributions of connection strengths together with their confidence intervals for the three strongest and most significant connections over subjects. For the model with RHIP inclusion, these are LIPC to PCC, LIPC to RIPC, and RIPC to RHIP. For the model with LHIP inclusion, these are LIPC to PCC, RIPC to PCC, and LHIP to PCC (Figures 7A,B, marked as “strongest”). We also examined connections which were strong and significant after the BMA and became non-significant after the BPA (Figures 7A,B, marked as “BPA^∗”). It can be seen that there is a high degree of consistency over subjects for the strongest connections. Usually not more than 2–3 connections had the opposite sign. For the insignificant correlations? after the BPA connections, the picture is not so clear: for example, the LIPC to mPFC connection (Figure 7A) in 6 subjects is negative (with a big negative outlier), and the LHIP to LIPC connection (Figure 7B) is negative in 5 subjects. Despite this fact, the behavior of RIPC to mPFC (Figure 7A), LHIP to RIPC and mPFC to LIPC connections seems quite stable and determined over subjects, assuming that the BPA analysis produces here counter-intuitive results. Possibly, this is due to strong outliers in the group, or there could be other reasons. One of the reasons of such counter-intuitive behavior is described in (Kasess et al., 2010): BPA takes into account the posterior covariance structure and at high signal-to-noise ratios (SNRs), these covariances can have a profound impact on BPA results. They also say that “for SNRs greater than 2 the group estimate of connection strength provided by BPA increasingly underestimates the mean of the parameter distribution even though the single-subject estimates were quite accurate. Moreover, for SNRs >5, the average of the modulatory connection lies outside the actual range of the individual parameter estimates.” Figure 7 shows that single-subject estimates are rather accurate. For all 26 subjects, in each voxel, we calculated the mean and the standard deviation of the corresponding time series to determine the SNR as in Welvaert and Rosseel (2013). The absolute voxel-wise minimum SNR value among all subjects was 2.4, ranging from this minimum to hundreds of units, which is consistent with resting-state time-series SNR (Welvaert and Rosseel, 2013). These findings, together with Figures 7A,B, lead us to the conclusion that BPA is not the best way to calculate parameter averages across subjects in our study. So, we preferred BMA (which is simply a weighted average) and performed further model analysis and discussion based on the BMA results.

The target regions in the current work were LHIP and RHIP. Our data on functional connectivity in part echoes that on causal connections showing that LHIP is more involved in the DMN than RHIP (see Figure 5). LHIP has connections with all key DMN regions. In contrast, RHIP has only two significant functional connections, whereby the stronger one is its connection to LHIP. Our functional connectome also includes links not reflected in the patterns of effective connectivity, e.g., a number of strong correlations of PCC with the activity of other key DMN nodes. Functional connections of PCC with LIPC and RIPC are stronger than LIPC and RIPS correlational links to hippocampal regions. Functional connections between RHIP and PCC, LIPC and RHIP, RHIP and mPFC are not significant.

Despite the fact that in the current work, the effective connectivity between two hippocampal regions was not directly assessed, one model combining the best LHIP and RHIP inclusion patterns (a_1 for RHIP and a_6 for LHIP with totally six nodes) was considered. Leaving the thorough analysis of both hippocampal regions inclusion for future research, we inverted only this particular model (without performing BMS), calculated parameter averages over subjects and directly compared LHIP and RHIP connections to the other four DMN nodes. It was found that: (a) LHIP/RHIP connections that were significant in the best 5-source model remain the same in 6-source model (except one RIPC to RHIP connection that fell slightly below the significance threshold), (b) connections from LHIP to mPFC, LIPC and RIPC are significantly (p < 0.01) stronger than the same connections from RHIP, (c) there is a very strong (0.44 Hz) and highly significant (p < 0.001) connection from LHIP to RHIP and no connection in the opposite direction. In our view, this preliminary analysis not only supports the overall picture of the discovered asymmetry, but also could be a good starting point to more exhaustive analysis of models with both hippocampal regions.

In this article, we extended our previous work on the causal relations within the DMN (Sharaev et al., 2016) by adding to its dynamic structure left and right parahippocampal regions. We also tested computational stability of our basic model by virtually removing one of the key DMN regions from the network. A discrete cosine basis set was used to model low frequency fluctuations in the brain, together with the spectral DCM to model the effective connectivity between the key DMN regions and the new target structures. The best models at the group level strongly suggested the same connectivity pattern between the four basic DMN regions – mPFC, PCC, LIPC, and RIPC, as revealed previously (Sharaev et al., 2016), and shed new light on the role of hippocampus in the resting state.

From the start of DMN studies, most researchers included HF in its composition (Greicius et al., 2004; Vincent et al., 2006). In addition, the DMN is often considered as a network implementing the basic state of consciousness (Gusnard and Raichle, 2001; Gusnard et al., 2001). In this role, it has to be closely related to HF which is crucial to basic cognitive functions of episodic memory (Tulving, 1985; Başar and Düzgün, 2016) and representation of surrounding space (Burgess et al., 2000; Moser and Moser, 2008). However, only a few reports addressed the lateralization of the DMN effective connections with respect to HF (Cui et al., 2015). By using spectral DCM as a tool, we for the first time found evidence for a strong asymmetric pattern in such relationships. LHIP demonstrated a high involvement in the DMN activity, with information outflow preponderant to all other DMN regions including RHIP, as shown by our preliminary analysis of 6-nodes interaction. Causal interactions of LHIP were bidirectional only in the case of LIPC. On the contrary, RHIP was mainly affected by inputs from LIPC and RIPC. In our view, this pattern of asymmetry in effective connections of the hippocampal regions may be related to lateralization phenomena in verbal and spatial domains documented in human neurophysiology, neuropsychology, and neurolinguistics (Howard and Templeton, 1966; Luria, 1966; Harrison, 2015). An obvious example of a drawback of such lateralized architecture is its potential vulnerability to unilateral destruction in functions of RHIP that could lead to the left-sided spatial hemi-neglect.

The authors of earlier studies often reported conflicting evidence on effective connectivity within the DMN (Jiao et al., 2011; Li et al., 2012; Di and Biswal, 2014). Models of effective connectivity from the present study are slightly different from the previous models and their power connections. For the first time, an inhibitory link was found between MPFC to LIPC (Figures 5 and 6), which may indicate the effort of the subjects to maintain resting state. This conclusion is supported by the data on dysfunction of the links between MFC and LIPC regions in patients with ADHD (Franzen et al., 2013). According to our results, bidirectional links between MPFC and PCC are relatively weak and statistically insignificant (Figure 5). Besides, our best reduced model remains stable even without the PCC node, when the link between mPFC and PCC is missing (Figure 6). The results are not consistent with the well-known strength of the structural links between the MPFC and PCC (Khalsa et al., 2014). We think that this is due to weakening effective connections of MPFC and PCC at rest, as were previously shown by Li et al. (2012) and Razi et al. (2015), although there are alternative data (Bajaj et al., 2013). Our conclusion can be supported by the data on the weakening of connections between functional networks in the resting state and their strengthening in the task performance (Di et al., 2013). On the other hand, our functional connectome provides a good average relatedness of these sites (Figure 6B), indicating a possible instability of the causal relationships in models which may indicate the effort with strong connections between MPFC and PCC.

Our extended models demonstrate the strong influence of both multimodal IPC zones on other DMN nodes (Figure 5), confirming previous reports (Jiao et al., 2011; Di and Biswal, 2014; Sharaev et al., 2016). Reduced models in their structure preserves these multimodal units and the majority of their bonds (Figure 6A). LIPC and RIPC influence MPFC and PCC, as well as interact with the hippocampal regions, doing this, as stressed above, in an asymmetric manner. The functional connectome also identifies asymmetric interactions with LIPC, RIPC and the hippocampal regions. LIPC interacts with LHIP more intensively than with RIPC with RHIP. At the same time, RIPC strongly interacts with LHIP than with RHIP (Figure 6B). The best reduced DCM model retains only the left hippocampal region, which affects all nodes of the reduced connectome and receives feedback from LIPC (Figure 6A). Functional connectivity confirms weak links of the right hippocampal region to other network nodes, except symmetrical structure of IPC. Some deviation in patterns of functional and effective connectivity is not surprising because of the large difference in temporal resolution of respective methods and a correlational character of the functional links.

It has to be noted that structural studies also demonstrate asymmetry in hippocampal regions such as differences in volumes (Pedraza et al., 2004; Shi et al., 2009). However, being dependent on a number of factors including handedness, age and anatomical position along the long axis (Lazić et al., 2015; Robinson et al., 2016), these differences seem to be less systematic than differences in functional activity. In addition, there is a kind of dissociation between the size of the hippocampus and its functional activity because LHIP, which functionally is more active, often has smaller volume than RHIP (Shi et al., 2009). One explanation for this discrepancy is that the increased activity of the left hippocampus in healthy subjects can lead to hippocampal volume reduction on the same side due to physiological pruning (Lazić et al., 2015; Okada et al., 2016). Functional lateralization seems to be a cause and not consequence of the anatomical difference. Our DCM data is consistent with the prominence of the LHIP functional activity. Of importance here is the known fact that the left parahippocampal region is primarily associated with verbal functions (Bernier et al., 2016), while the right one is involved in visual memory processes (Barnett et al., 2015). We are therefore inclined to regard the functional and effective lateralization of HF with the role of verbal memory and inner speech processes in the social meditation. The latter may be the very essence of the DMN activity as a whole (Li et al., 2014; Xie et al., 2016).

In general, our experimental DMN architecture can be interpreted as an optimally tuned system of operational rest (Bernhofer, 2016). It is characterized by a division of labor between left and right parahippocampal regions. RHIP plays a more receptive role by receiving streams of information from both IPCs and by building on the basis of this multimodal sensory input a bilateral representation of space. This spatial representation role seems to be restricted in the case of LHIP to representation of only the contralateral hemi-space. Having access to verbal episodic memory, LHIP takes over the function of the main driving force of activation within the DMN. In particular, it influences the decision-making center in MPFC (Lee and Seo, 2016). The latter is configured to behavioral rest and therefore inhibits sensory-motor structures in the leading hemisphere, i.e., LIPC (see green arrow in Figure 5 left). As our model summarizes the data of a relatively large group of healthy subjects, it does not include the individual and temporal characteristics that must be investigated further.

Our study design is not without limitations. Due to the number of models for the simultaneously considered regions, we were limited by the causal analysis to not more than five nodes at once. For this reason, the present study left some of the obvious candidate connections without an in-depth analysis. In addition, the posterior hub of DMN, PCC, happened to play a strictly receptive role, not being as active in influencing other regions, as expected initially. This could be partially reflected in the results of our extraction procedure (see Figure 6A). Finally, it would be interesting to find functional and causal relations between DMN and the task-positive network (i.e., elucidating negative DMN connectivity - see Di Perri et al., 2016). This would be another obvious objective of such type of analysis. Similarly, to a study of individual and situational DMN properties, however, research related to task performances and goal-directed behavior would need the application of a larger battery of methods combining slow and fast neuroimaging approaches. Another potentially questionable point is our choice of averaging method – BMA instead of BPA. Despite the fact that BPA takes into consideration parameters (co) variance structure and thus more information about the estimate precision and parameter interdependence, BPA seems to be less applicable when examining systems with high SNR (>5).

As mentioned above, this work has two objectives. Besides investigating the HF place within the effective connections of the key DMN regions, we tested the stability of our models when adding or removing regions of interest. In the end, we could emphasize the finding that adding additional sources (LHIP or RHIP), as well as removing one of the regions (PCC in this case), does not change dramatically the overall connectivity pattern. Even if one takes into account several limitations of the study, this computational stability is surprising. It shows the potential of the spectral DCM as a reliable tool for neurocognitive research in health and in clinical conditions. Our main conclusion is therefore that the method can be of importance for constructing complex causal hypotheses (models) from simple connectivity circuits. In a similar vein, possibly, sophisticated psychological effects could be decomposed and evaluated in a parametric way.

VU: Work concept, article text (introduction, results, discussion), final article approval, accountable for all aspects of the work. MS: Mathematical methods, programming, data analysis, modeling, article text (introduction, methods, results) and pictures preparation, accountable for all aspects of the work. VZ: Experimental design, data acquisition and preprocessing, final article approval, accountable for all aspects of the work. SK: Data acquisition, article text (introduction), final article approval, accountable for all aspects of the work. VV: Experimental design, article text (introduction, discussion), final article approval, accountable for all aspects of the work. BV: Overall guidance, project idea, article text (introduction, neuropsychological aspects of discussion), final article approval, accountable for all aspects of the work.

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.