fMRI in Non-human Primate: A Review on Factors That Can Affect Interpretation and Dynamic Causal Modeling Application.

Dynamic causal modeling (DCM)-a framework for inferring hidden neuronal states from brain activity measurements (e. g., fMRI) and their context-dependent modulation-was developed for human neuroimaging, and has not been optimized for non-human primate (NHP) studies, which are usually done under anesthesia. Animal neuroimaging studies offer the potential to improve effective connectivity modeling using DCM through combining functional imaging with invasive procedures such as in vivo optogenetic or electrical stimulation. Employing a Bayesian approach, model parameters are estimated based on prior knowledge of conditions that might be related to neural and BOLD dynamics (e.g., requires empirical knowledge about the range of plausible parameter values). As such, we address the following questions in this review: What factors need to be considered when applying DCM to NHP data? What differences in functional networks, cerebrovascular architecture and physiology exist between human and NHPs that are relevant for DCM application? How do anesthetics affect vascular physiology, BOLD contrast, and neural dynamics-particularly, effective communication within, and between networks? Considering the factors that are relevant for DCM application to NHP neuroimaging, we propose a strategy for modeling effective connectivity under anesthesia using an integrated physiologic-stochastic DCM (IPS-DCM).

Introduction

Neuroimaging analyses in humans and non-human primates (NHP) have become increasingly sophisticated. One such innovative image analysis technique is dynamic causal modeling (DCM) which has been applied to human fMRI data (Friston et al., 2003; Rowe et al., 2010; Boly et al., 2012; Havlicek et al., 2017; Park et al., 2018; Tak et al., 2018). DCM is a method for inferring hidden neuronal states from brain activity measurements (e.g., fMRI) and their context-dependent modulation (Stephan et al., 2010). Using a Bayesian framework, DCM generates a predicted time serious using a set of differential equations to model neural dynamics. Then, one estimates model parameter by optimally fitting the predicted time series with the observed data. Using DCM, one can test mechanistic hypotheses about how the observed data was generated.

Utilizing a Bayesian approach, model parameters are estimated based on prior knowledge of conditions that might be related to neural and BOLD dynamics (e.g., requires empirical knowledge about the range of plausible parameter values). As such, we discuss: (1) what factors need to be considered when applying DCM to NHP (keeping in mind that it was designed for human fMRI data); and (2) considering said factors, what strategies can one implement when modeling effective connectivity to fMRI data recorded under anesthesia. The latter is an important consideration in NHP fMRI, since most imaging experiments are done under anesthesia and anesthetics have been demonstrated to impart changes in BOLD and neural dynamics, particularly the inhibitory drive (Martin et al., 2006; Masamoto et al., 2007; Moran et al., 2011; Aksenov et al., 2015; Paasonen et al., 2018).

Herein we discuss the following: (1) DCM overview and motivation for NHP data application; (2) comparison of cerebrovascular architecture and functional networks between humans and NHPs; (3) effects of anesthetics on vascular physiology, BOLD contrast, and neural dynamics—focusing on effective communication within and between networks. Lastly, we propose a strategy for modeling effective connectivity under anesthesia using an integrated physiologic-stochastic DCM (IPS-DCM). Taking all these together, we can avoid pitfalls in DCM application in NHP data and facilitate more accurate interpretation of the observed neural dynamics as measured through BOLD fMRI.

Dynamic Causal Modeling (DCM) in Non-human Primates

DCM Overview

This section discusses an overview of DCM, including its assumptions, and the motivations for NHP data application. DCM, developed by Friston et al. (2003), is a method that allows estimation of network dynamics and how the dynamic neuronal states give rise to the measured data (Heinzle and Stephan, 2018), as well as how functional coupling is affected by experimental factors. When applied to fMRI data, the BOLD signal is considered a measurable (observed) variable (y) of the underlying neural activity (z) that is not directly measured with fMRI; hence, neural activity is considered a “hidden state variable” (Kahan and Foltynie, 2013). As a measure of effective connectivity (EC), DCM quantifies the directed (causal) influence between regions (Friston, 2009) and considers the rate of change of neural activity over time (z) in response to incoming signals—afferents from other brain regions, experimental manipulation, or both. This is done by creating generative models which are plausible models of how the observed BOLD signal may be generated by those influences.

In DCM, one needs to specify which regions to include in the overall model. As such, DCM follows an analysis that can address which regions in the brain an experimental manipulation induces changes in BOLD response such as a general linear model (GLM) (Stephan et al., 2004). In GLM, time series from each voxel are fitted with an experimental design matrix wherein each condition is specified and the relationship between neural and BOLD response is modeled via a hemodynamic response function (Pernet, 2014). When there are no contrast differences between conditions, there is no motivation to do DCM. Once the relevant regions are identified with GLM, the time series from each region of interest is extracted. After which, one specifies the model architecture—the location of intrinsic connections, driving, and modulatory inputs. Then, one proceeds to estimate the parameters of the generative model.

DCM can offer a more accurate modeling of network dynamics from fMRI data as it considers both the distributed neuronal interactions, and the transformation of the neuronal dynamics to the measured (BOLD) signal (Havlicek et al., 2015). This is accomplished through three fundamental components of the DCM generative model: (1) neuronal (connectivity) model, (2) hemodynamic “Balloon” model, and (3) BOLD signal change equation (Stephan et al., 2007; Havlicek et al., 2015; Friston et al., 2017). A discussion on strategies for modifying the generative model for NHP fMRI under, particularly under anesthesia, will be discussed in a later section.

Each region or node is represented by (differential) neural state equations that provide an abstraction of the summed activity of a neuronal population. The neural state equation is comprised of three parameters that embody: (1) intrinsic connectivity among regions in the absence of input (A matrix); (2) direct influence of extrinsic input on neuronal activity (C matrix); and (3) context-dependent change in connectivity induced by the input (B matrix). Neural dynamics in DCM have been characterized by single-state, two-state, and adaptive two-state equations. The single-state neural equation (Friston et al., 2003) of the classical DCM models synaptic activity of a single excitatory neuronal population in a cortical region—with the rationale that most cortico-cortical connections are excitatory. Two-state (Marreiros et al., 2008) and adaptive two-state (Havlicek et al., 2015) neural equations model both excitatory (glutamatergic) and inhibitory (GABAergic) connections within each region, which makes it a more physiological representation of neural population dynamics. In contrast to the two-state model, the adaptive two-state model includes parameters which allow for adaptation and refractory effects of the neuronal response (Havlicek et al., 2015). These are described in detail in Friston et al. (2003); Marreiros et al. (2008); Havlicek et al. (2015).

The neural model is combined with the hemodynamic model (Friston et al., 2003)—which incorporates the dynamics of neurovascular coupling and the Balloon–Windkessel model (Buxton et al., 1998; Mandeville et al., 1999). This biophysical forward model provides region-specific estimates of the translation of neuronal activity into a predicted BOLD response—as such, hemodynamic states are a function of the neuronal state/s of each region (Friston et al., 2003). The hemodynamic model is comprised of four differential equations that characterize how for each region pre-synaptic activity drives hemodynamic responses—which are mediated by astrocytic Ca2+ signaling whose endfeet release vasodilatory metabolites (MacVicar and Newman, 2015). This vasodilatory signal is subject to auto-regulatory feedback (Attwell and Iadecola, 2002; Friston et al., 2003) and blood flow changes proportionately to the vasodilatory signal which leads to an increase in blood volume with concomitant decrease in deoxyhemoglobin (Stephan et al., 2004).

The output signal is defined by the BOLD signal change equation—which links blood volume and deoxyhemoglobin content to the BOLD signal change (Stephan et al., 2007). The resulting BOLD signal depends on the depends on the relative contributions of intra- and extravascular signals (Buxton et al., 2004; Friston et al., 2017) and follows the flow dynamics with a delay of ~1 s (Stephan et al., 2004). The BOLD signal change equation completes the DCM generative model.

In standard DCM, the parameters of the full forward model (combined neural and hemodynamic models) is estimated from the measured BOLD data using is estimated iteratively from the measured BOLD data using Variational Bayesian (VB) algorithm (Friston et al., 2007) to produce probabilistic estimates of the expected value of each parameter given the data. Details of the parameter estimation process are beyond the scope of this review and the reader is referred to Friston et al. (2002). The objective of the estimation process is to generate a predicted signal that closely matches the observed BOLD data (Kahan and Foltynie, 2013). From the DCM parameter estimates, one can then quantify and make inferences about connection strength and direction between regions based on prior assumptions about connectivity architecture.

Finally, Bayesian model selection (BMS) is done to compare a series of models that represent different prior hypotheses of connectivity architecture to examine which of the competing models most likely generated the observed data. This is done by comparing the evidence for each model. Model evidence—the normalization constant for the product of the likelihood of the data and prior probability of the parameters—is approximated using Laplace approximation in standard DCM (Stephan et al., 2007). This yields the log-evidence for each model, characterized as the conditional probability (posterior probability) of competing models given the observed data (Kahan and Foltynie, 2013). The optimal model is one that represents the best compromise between model fit (accuracy) and complexity (characterized by the number of free parameters in the model) (Stephan et al., 2007)—it is the best fitting, yet, most parsimonious model. As such, this model is more generalizable.

DCM Assumptions

There are three main assumptions in DCM: (1) deterministic assumption on the inputs that enter the system; (2) Gaussian assumption on the posterior density; and (3) assumption of equal detection of BOLD signal changes throughout the brain. Firstly, in classical DCM, all processes in the system are considered deterministic such that it is presumed that neural dynamics in a region is entirely due to incoming afferents from other regions and/ or experimental inputs. This assumption is not always tenable in NHP fMRI studies which are mostly done under anesthesia—this restricts task-based imaging to simpler visual/ auditory experiments and many imaging experiments are task-free. However, extensions of DCM in humans have allowed its application in resting state studies: (1) using classical DCM, one can stimulate nodes with fluctuations of specific frequencies (Di and Biswal, 2014); (2) modeling random fluctuations in neural dynamics explicitly using stochastic DCM (Li et al., 2011; Daunizeau et al., 2012); or (3) estimating the spectral density of neuronal fluctuations such that effective connectivity of hidden neuronal states is a function of observed functional connectivity from hemodynamic responses (Friston et al., 2014; Park et al., 2018). Secondly, parameter estimates of the generative model are assumed to be Gaussian—i.e., that the values are normally distributed. This may be a concern when using two-state DCM which use exponentiated scale parameters that introduce positivity constraints as these values likely do not have a normal distribution (Hillebrandt et al., 2014). Lastly, DCM rests on the assumption that BOLD signal detection is equally sensitive across brain regions. This was validated by Friston et al. (2003) by simulating region-specific dropout wherein they found that DCM does not accommodate substantial signal dropout (e.g., 50%). In the absence of profound dropout, DCM is robust to regional variations in sensitivity to BOLD signal changes.

Why Apply DCM to NHP FMRI?

The motivations for applying DCM to NHP data are founded on the prospective to make DCM models more accurate with animal experiments that are not yet feasible or are considered unethical in humans, as well as on being a more physiologically-informed characterization of network dynamics vs. previous effective connectivity models.

Animal imaging data and experiments possess the potential to make DCM models more accurate. Invasive procedures such as in vivo optogenetic and electrical stimulation can be combined with functional imaging to examine how disruptions at the microscale can affect whole brain network dynamics. Additionally, one can also investigate the temporal evolution of the effects of drugs (e.g., administered intravenously or through gas inhalation). Moreover, advances in high-field fMRI at the sub-millimeter scale, allows imaging and modeling at laminar resolutions (Heinzle et al., 2016; Friston et al., 2017).

In order to accurately disentangle neuronal- and hemo-dynamics, experimental manipulation can provide necessary constraints to estimate hemodynamic parameters, which may explain discrepancies in response shape between neuronal and BOLD response. This can be addressed by using multi-modal recordings (e.g., simultaneous measurements of CBF, CBV, and BOLD), which is often more feasible in animal studies, including NHP. For instance, Havlicek et al. (2017) showed that combined analysis of BOLD and CBF data yields more robust effective connectivity estimates.

Furthermore, DCM was originally developed specifically for fMRI data (Friston et al., 2003) which gives it an edge over other models such as Granger causality and structural equation modeling (SEM) that were initially applied in the fields of economics, psychology, and genetics (Wright, 1920; Granger, 1969). Compared with previous EC measures, DCM allows a more physiologically-informed characterization of network dynamics as it incorporates a hemodynamic model that has previously been experimentally validated (Buxton et al., 1998; Friston et al., 2003). On the other hand, the disadvantages of Granger causality in fMRI application is plentiful. Using four different algorithms, Smith et al. (2011) demonstrated that Granger causality exhibits poor (<20%) sensitivity in connection link detection, false positive identification and directionality estimation (Smith et al., 2011). Moreover, Witt and Meyerand (2009) found that Granger causality has poor sensitivity and specificity (close to chance levels) when modeling data including intrinsic variance from trimmed time series. While DCM is computationally costly, its ability to model non-linear and dynamic neuronal interactions (Bielczyk et al., 2019), as well as both unidirectional and bidirectional connections (Vaudano et al., 2013; Buijink et al., 2015) give a more accurate picture of underlying neuronal activity. While originally designed for task fMRI, recent developments such as stochastic or spectral DCM allow modeling of resting state fMRI (Li et al., 2011; Daunizeau et al., 2012; Friston et al., 2014; Park et al., 2018). Classical DCM (Friston et al., 2003) pose restrictions on network size as increasing the number of nodes considerably increases computational time (Bielczyk et al., 2019). However, ensuing extensions allow exploratory studies involving larger networks such as spectral DCM (Friston, 2011) for resting-state fMRI and DCM with sparsity constraints for task fMRI (Frässle et al., 2018). The many advantages of DCM over previous effective connectivity measures and the ability to combine invasive procedures in neuroimaging make DCM application to NHP data an exciting endeavor and may offer the possibility to improve the accuracy of causal models. A summary of the advantages and disadvantages of the commonly used effective connectivity models (structural equation modeling, Granger causality, transfer entropy, and dynamic causal modeling) is outlined in Table 1.

Table 1

Comparison of commonly used effective connectivity models.

	Pros	Cons
Structural equation modeling	* Can detect excitatory and inhibitory connections and connection strength (Bielczyk et al., 2019) * Sensitivity to small changes in path weight values due to large dynamic range (Witt and Meyerand, 2009)	* Difficulty in estimating reciprocal and cyclic connections (physiologically, reciprocal connections are ubiquitous in the brain) (Friston, 2011) * May not be as suitable to event-related design due to the assumption that random fluctuations change very slowly in relation to neuronal dynamics such that neuronal dynamics has already reached steady-state at the time of recording (Friston, 2011) * May be inappropriate in the context of disease or pharmacologic experiments that can affect hemodynamic response function (Rowe et al., 2010)
Granger causality	* Can detect excitatory and inhibitory connections and connection strength (Bielczyk et al., 2019) * Yields bidirectional connections (Bielczyk et al., 2019) * Results can be mapped onto the brain similar to fMRI (Goebel et al., 2003; Roebroeck et al., 2005; Witt and Meyerand, 2009)	* Poor sensitivity and specificty (close to chance levels) when modeling data including intrinsic variance from trimmed time series (Witt and Meyerand, 2009) * Assumption of signal stationarity (Seth et al., 2015) * Restriction on network size–the number of nodes divided by the number of shifts can never exceed the number of time points (Bielczyk et al., 2019) * Markovian assumption that random terms in the vector autoregression model are serially independent may not hold when the terms become temporally correlated upon converting from continuous to discrete time formulations (Friston, 2011) * The spatial distribution of GC has been associated with the Circle of Willis and identifies major arteries and veins as causal hubs (Webb et al., 2013) * Assumption of uniform hemodynamic response function across regions may elicit spurious causal relationship when one region has faster hemodynamic activity–the temporal precedence of the peak in one region may be mistaken for Granger causing the other (Bielczyk et al., 2019) * fMRI temporal resolution may be too slow for accurate depiction of neural dynamics using Granger causality (Witt and Meyerand, 2009) * Poor (<20%) sensitivity in connection link detection, false positive identification and directionality estimation (Smith et al., 2011)
Transfer entropy	* Can detect excitatory and inhibitory connections and connection strength (Bielczyk et al., 2019) * Captures linear and non-linear interactions between nodes (Bielczyk et al., 2019) * Computationally cost-efficient (Vicente et al., 2011)	* Restriction on network size–the number of nodes divided by the number of shifts can never exceed the number of time points (Bielczyk et al., 2019) * Imposes a time-lag in the inference procedure with similar disadvantages as Granger Causality in fMRI application (Schreiber, 2000)
Dynamic causal modeling	* Developed specifically for fMRI (Friston et al., 2003) data and incorporates a biologically-informed model of BOLD dynamics (Buxton et al., 1998), unlike other models Granger causality and SEM were originally applied in the fields of economics, psychology, and genetics (Wright, 1920; Granger, 1969) * Can detect excitatory and inhibitory connections and connection strength (Bielczyk et al., 2019) * Can model both unidirectional and bidirectional connections (Vaudano et al., 2013; Buijink et al., 2015) * Models nonlinear and dynamic neuronal interactions (Bielczyk et al., 2019) * Classical DCM is suitable for event-related designs (Rowe et al., 2010) * Stochastic or spectral DCM is suitable for resting state studies (Li et al., 2011; Daunizeau et al., 2012; Friston et al., 2014; Park et al., 2018) * For exploratory studies involving larger networks, spectral DCM (Friston, 2011) can be applied for resting-state fMRI while DCM with sparsity constraints can be applied for task fMRI (Frässle et al., 2018) * High reproducibility (Rowe et al., 2010; Schuyler et al., 2010; Bernal-Casas et al., 2013; Tak et al., 2018)	* Computationally-expensive (Bielczyk et al., 2019) * Restriction on network size (using classical DCM)—increasing the number of nodes considerably increases computational time (Bielczyk et al., 2019) * Depends on prior assumptions on connectivity architecture (Friston et al., 2003) * Assumes all models are equally likely (even implausible models) (Lohmann et al., 2012); hence, substantial knowledge is needed to define all plausible causal connections between nodes

Considerations Specific to DCM Application in NHP

Functional Networks and Cerebrovascular Architecture: Human vs. NHP

Cerebrovascular Architecture and Physiology

The cerebral vasculature of NHP is largely similar to humans in terms of the architecture of superficial pial vessels and intracortical vessels (Duvernoy et al., 1981; Weber et al., 2008). Similarly found in humans (Duvernoy et al., 1981), large vessels are found on the surface of the macaque brain (Weber et al., 2008). The density of superficial vessels have been found to vary across cortical regions—the occipital lobe surface is highly vascularized compared to the less dense vascularization at the top of the hemispheres near the interhemispheric fissure (Scharrer, 1960; Duvernoy et al., 1981). In terms of vessel diameter, pial veins generally have a larger diameter than arteries— central veins have an average diameter of 280–380 μm and peripheral veins average 130 μm, while central arteries have a diameter of 260–280 μm and peripheral arteries average 150–180 μm (Duvernoy et al., 1981; Guibert et al., 2010). As for intracortical vessels, the laminar distribution of vessels in NHP are similar to that found in humans. Studying the primate visual cortex, Bell and Ball (1985) found a high density of cortical vessels in layer IVC that ends at the boundary between primary and secondary visual cortices—akin to the vascular distribution in humans.

The fluctuations in deoxyhemoglobin concentration detected by the BOLD contrast depends on the combined changes in cerebral blood volume (CBV), cerebral blood flow (CBF), and cerebral metabolic rate of oxygen (CMRO2) (Buxton et al., 2014). In general, the values for each hemodynamic component is comparable between humans and NHP. Van Aken et al. (1986) found CBF at 48 ± 4 mL/100 g/min in baboons. In humans, Olsen et al. (1994) demonstrated 51 mL/100 g/min, while Ito et al. (2004) found 69.8 ± 15.4 mL/100 mL/min CBF values. CBV ranged from 3.5 to 4.7 ml/100 g of brain tissue in the macaque (Phelps et al., 1973; Grubb et al., 1974; Eichling et al., 1975) while CBV = 3.8 ± 0.7 ml/100 ml−1 (Ito et al., 2004). As for CMRO2, Van Aken et al. (1986) obtained 3.64 ± 1 (ml/100 g/min) in baboons and Olsen et al. (1994) found 3.5 (ml/100 g/min) in humans. The values for CBF and CMRO2 are modulated by anesthetics and this will be discussed in an ensuing section.

Functional Differences

Resting-state fMRI studies comparing human and NHP connectivity have revealed three main differences. First, there are differences in specific connectivity patterns that may potentially indicate cognitive specializations in humans. For example, comparing the organization of the dorsal frontal cortex between humans and macaques (Sallet et al., 2013), found a high degree of similarity in functional coupling patterns between the medial frontal cortex and other regions (i.e., frontal pole, medial prefrontal, and dorsal frontal convexity) in both. However, certain regions in the dorsolateral prefrontal cortex (areas 9/46) were coupled with the superior and medial parietal cortex in humans but not in macaques. Further, Mars et al. (2011) demonstrated resting-state functional connections between anterior prefrontal cortex and central inferior parietal lobule (IPL) in humans which are not found in macaques. Second, there may be species-specific differences in cortical hub distribution. Upon mapping putative hubs in humans, chimpanzees, and macaques, Li L. et al. (2013) demonstrated the ventrolateral prefrontal, medial parietal and retrosplenial cortices are hubs across three species. In contrast, medial prefrontal, inferior parietal, and V1 cortices were hubs in macaques and chimpanzees and not in humans. Additionally, superior parietal and medial premotor cortices were hubs in humans and not in the NHPs. Third, though it has been demonstrated that there are 11 functionally correspondent networks in both humans and macaques, three networks were found in humans that are missing in the latter (Mantini et al., 2013). While sensory-motor, attention, language, and default mode networks are evolutionarily conserved, two lateralized fronto-parietal networks are unique to humans (Van Essen and Dierker, 2007). These have been implicated in general intelligence (Duncan et al., 2000), abstract reasoning (Dehaene et al., 2003), and tool use, particularly retrieving and planning transitive actions for subsequent hand motor behavior (Frey, 2008). A third human-specific network includes the anterior insula and dorsal anterior cingulate cortex—both putatively involved in empathy (Singer and Lamm, 2009).

The Anesthetized Brain

Anesthetic Effect on Vascular Physiology and BOLD Contrast

Anesthetics have been found to impart changes in cerebrovascular function and the BOLD signal. As these are incorporated in the DCM generative model, it is important to examine how these are affected in the anesthetized brain, particularly since most NHP neuroimaging is done under anesthesia.

There are conflicting results in the regional distribution of anesthetic-induced changes in vascular physiology. Li C.-X. et al. (2013) examined the dose-dependent effect of isoflurane on regional CBF of cortical and subcortical structures in macaques. They found that high isoflurane concentrations (i.e., 1.5%) resulted in global CBF increase which was most evident in subcortical structures—specifically in the thalamus and cerebellum in macaques. Interestingly, under the 0.75–1.5% isoflurane maintenance doses, there were no observable CBF changes in cortical regions (i.e., anterior cingulated cortex, motor cortex, medial prefrontal cortex) and the caudate. These indicate that while CBF auto-regulation is intact in cortical regions and the caudate under isoflurane maintenance dose, it is impaired in the thalamus and cerebellum, and suggest that subcortical structures contribute the most to the increase in global CBF. On the other hand, Långsjö et al. (2005) demonstrated that ketamine increased whole brain CBF in humans—with the highest increase in the anterior cingulate. They found that CMRO2 increased only in the frontal cortex while glucose metabolism increased only the in the thalamus. Långsjö et al. (2005) assert that this indicates that majority of the increases in CBF most likely do not indicate neuronal activation. Another study by Van Aken et al. (1986) showed that the impact of isoflurane on CBF was biphasic—low levels (0.5 ± 0.35 vol%) resulted in vasoconstriction and decreased CBF while higher concentrations (0.95 ± 0.7 vol% and 1.4 ±1 vol%) caused vasodilation and increased CBF to baseline levels (no anesthesia). While they found that the effect of isoflurane on CBF was biphasic, CMRO2 continually decreased in a dose-dependent manner. The disparity in the findings of these studies may be attributed to methodological differences in anesthetic (isoflurane vs. ketamine), species (human vs. primate), and imaging modality (arterial spin labeling MRI vs. PET). Nonetheless, the results of these studies indicate anesthetic-related disruption of CBF-metabolism coupling characterized by more widespread CBF increases with minor changes in CMRO2 and glucose metabolism, with dose-related variabilities in response.

As for the impact of anesthetics on BOLD contrast, anesthetics have been demonstrated to cause changes in the spectral components and decrease the activated area and magnitude of the signal, as well as change the hemodynamic response temporal structure (Martin et al., 2006; Aksenov et al., 2015; Paasonen et al., 2018). Paasonen et al. (2018) found that awake mice exhibit higher spectral BOLD power at a wide frequency range. In contrast, they found that all six anesthetic conditions studied strongly suppressed power and BOLD fluctuations occurred at narrower frequency ranges, which potentially reflect more homogeneous activity. This was a predictable outcome as by definition, anesthetics exert an inhibitory effect on brain function. Aksenov et al. (2015) suggested that enhanced GABA-A receptor inhibition and diminished afferent input reduces the spread of stimulus-related activity, which results in more focal activation and decreased BOLD response area. Further, they surmised that decreased BOLD response magnitude reflects reduced thalamic input and intra-cortical processing associated with decreased neuronal excitation. A decrease in hemodynamic response magnitude under anesthesia has also been demonstrated in previous studies employing optical imaging spectroscopy (OIS) (Berwick et al., 2002). Aside from decreased duration of the BOLD signal reported by Aksenov et al. (2015), Martin et al. (2006) found that anesthesia increased hemodynamic response latency (~2 s awake vs. 4 s anesthetized).

Anesthetic Impact on Effective Communication Within and Between Networks

It is now widely acknowledged that anesthetics modulate communication within and between networks. Anesthetic-induced perturbations in effective connectivity have been observed across multiple neuroimaging techniques, species (humans and rodent models), and anesthetics. In particular, effective connectivity changes have been demonstrated in thalamo-cortical, cortico-thalamic, and cortico-cortical (both association and sensory cortices) connections. Under propofol, ketamine, or isoflurane, selective disruption of frontal to parietal feedback has been widely replicated in human subjects and rats (Imas et al., 2005; Lee et al., 2009; Boly et al., 2012)—particularly gamma frequencies (50 Hz) (Imas et al., 2005). In addition, Gómez et al. (2013) found that propofol decreased feedback from middle frontal gyrus to superior temporal gyrus. As for sensory-motor regions, multiple studies found impaired effective connectivity under anesthesia. White and Alkire (2003) demonstrated impaired effective drive from the supplementary motor area (SMA) to the primary motor (M1) cortex in humans, while Kang et al. (2016) showed abolished causal flow of 7–12 Hz activity from primary sensory (S1) to M1 and ventrobasal thalamic nucleus in mice. Decreased feedforward and feedback connections were found within auditory cortical regions (Heschl's gyrus and superior temporal gyrus) in humans under propofol (Gómez et al., 2013). Pertaining to anesthetic-related changes in thalamocortical effective connectivity, the results are more heterogeneous. Causal flow from the thalamus to the anterior cingulate and posterior parietal cortices was intact under propofol (Boly et al., 2012). On the other hand, effective connectivity from the thalamus to the SMA was impaired under halothane or isoflurane (White and Alkire, 2003). However, effective thalamic drive to M1 has shown variable results—remaining unchanged under ketamine in mice (Kang et al., 2016) compared to the awake state, while decreasing under halothane or isoflurane in humans (White and Alkire, 2003). This discrepancy may indicate species-specific differences in thalamic-motor connectivity in response to anesthetics, differences in regional distribution of different anesthetics (ketamine vs. halothane or isoflurane), or both. Altogether, these findings indicate that anesthetics induce: (1) disruption of higher-order information processing; and (2) reduced capacity for sensory-motor integration.

Relevance to DCM in NHP

For DCM to generate reliable neuronal signal estimates, having an accurate model of the hemodynamic response is crucial; thus, it is important to assess whether differences in cerebrovascular architecture and physiology in NHP—factors that contribute to the BOLD contrast—may affect priors on the biophysical parameters of the model.

The cerebrovascular architecture and laminar distribution of blood vessels in NHPs are similar to that of humans. The fMRI BOLD signal indicates changes in deoxyhemoglobin and the signal is sensitive to fluctuations in venous blood volumes. The cerebral vasculature of NHP is largely similar to humans in terms of the architecture of superficial pial vessels and intracortical vessels (Duvernoy et al., 1981; Weber et al., 2008; Adams et al., 2015). By extension, the BOLD signal-to-noise ratio in superficial to deep brain areas is expected to have similar distribution in both.

More importantly, the BOLD signal is a function of CBV, CBF, and CMRO2—which reflects fluctuations in deoxyhemoglobin content (Buxton et al., 2014), and DCM includes priors on these hemodynamic parameters. Since the CBV, CBF, and CMRO2 in NHP are comparable to human values (Phelps et al., 1973; Grubb et al., 1974; Eichling et al., 1975; Van Aken et al., 1986; Olsen et al., 1994), it is tempting to presume that one can proceed with applying DCM in NHP without adjustments on the biophysical parameters if one solely considers cerebrovascular architecture and physiology.

Notably, anesthetics such as isoflurane and ketamine have been found to increase blood flow (Van Aken and van Hemelrijck, 1991; Långsjö et al., 2005; Slupe and Kirsch, 2018), which is relevant in primate neuroimaging as it is usually done under anesthesia. The decision on how to proceed with DCM in the anesthetized brain depends on the research question. Upon comparing eight different hemodynamic models by modifying the coefficients of the BOLD signal change equation, Stephan et al. (2007) assert that variations in the hemodynamic model are relatively inconsequential when making inferences on underlying neuronal dynamics and their causal influences. This is due to the relative independence of the parameters of the neural state equation—characterizing intrinsic connectivity and their context-dependent modulation (A and B matrices)—from the amplitude of the hemodynamic response and the direct influence of extrinsic inputs on neuronal activity (matrix C). On the other hand, if the research question focuses on regional variations in hemodynamic parameters (e.g., in healthy vs. clinical populations), then preliminary tests need to be performed to determine the effect of variations on the priors of the biophysical model—related to anesthetic influence. For instance, if the experimenter can identify the baseline CBF and CBV using multi-modal recordings, the mean transit time (i.e., baseline CBV/baseline CBF)—which scales the CBV and deoxyhemoglobin changes in the hemodynamic model—could be directly calculated. Other strategies for modifying the hemodynamic model are discussed in the last paragraph of the next section on IPS-DCM parameters.

Moreover, caution is necessary when interpreting results if the model includes regions or networks in which effective connectivity may be modulated by anesthetics. Regions most affected by anesthetics are: (1) fronto-parietal (Ku et al., 2011; Boly et al., 2012; Kim et al., 2017); (2) sensory-motor (i.e., S1 to M1, SMA to M1, auditory cortical regions) (White and Alkire, 2003; Gómez et al., 2013); and (3) thalamocortical networks specifically involving somatosensory and motor function (White and Alkire, 2003; Kang et al., 2016). Since effective drive of these networks are reduced or abolished by anesthetics, it may result in underestimation of the effect of interventions (e.g., medications, brain stimulation) on connectivity within and between affected networks. It may also result in the overestimation of the modulatory effects of different interventions if aimed at reducing hyperconnectivity.

Additionally, anesthetics can modulate neural dynamics by altering inhibitory drive, neural refractory period, and cortical adaptation (Masamoto et al., 2007; Moran et al., 2011). These artifacts can be addressed by using an adaptive two-state neural connectivity model (Havlicek et al., 2015). This and other strategies for applying DCM in NHP under anesthesia are discussed in the succeeding section.

Integrated Physiologic-Stochastic DCM: Modeling Effective Connectivity Under Anesthesia

IPS-DCM Parameters

An important factor to consider upon implementing DCM in the anesthetized NHP is the type of DCM to use. Anesthetics can decrease BOLD signal-to-noise ratio and artificially lower spontaneous fluctuations and activity correlations (Hutchison et al., 2014). In addition to dose-dependent linear decrease in glutamatergic excitatory postsynaptic potentials (EPSPs) together with non-linear increase (saturating) in GABAergic inhibitory postsynaptic potentials (IPSPs) (increasing local inhibitory drive) (Moran et al., 2011), anesthetics modulate the neural refractory period and cortical adaptation with varying degrees depending on anesthetic type and dose (Masamoto et al., 2007). To overcome these limitations, we propose the combined use of two DCM extensions: stochastic DCM (Li et al., 2011) and physiologically-informed DCM (P-DCM) which has an adaptive two-state neural connectivity equation (Havlicek et al., 2015).

Stochastic DCM models fluctuations in hidden states (e.g., neuronal or hemodynamic states) that are due to endogenous (autonomous) dynamics, not exogenous experimental inputs (Li et al., 2011). By accounting for noise in the model using stochastic DCM, the risk of under-estimating effective connectivity parameters can be reduced (Gómez et al., 2013). Additionally, stochastic DCM is more useful when there are non-linear interactions among hidden states, such as the non-linear increase in GABAergic inhibitory drive (Moran et al., 2011; Daunizeau et al., 2012). Finally, Li et al. (2011) demonstrated that stochastic DCM using the generalized filtering (GF) inversion method showed higher sensitivity in detecting group differences compared to the variational Bayesian (VB) algorithm (Friston et al., 2007) of standard deterministic DCM.

Two-state (Marreiros et al., 2008) and adaptive two-state (Havlicek et al., 2015) neural equations model both excitatory (glutamatergic) and inhibitory (GABAergic) connections within each region, which makes it a more physiological representation of neural population dynamics. In contrast to the two-state model, the adaptive two-state model includes parameters which allow for adaptation and refractory effects of the neuronal response (Havlicek et al., 2015). This extension is relevant since anesthetics induce changes to the neural refractory period and cortical adaptation (Ogawa et al., 1992; Masamoto et al., 2007). Furthermore, explicitly modeling both excitatory and inhibitory signals may help capture the potentiation of GABAergic inhibition under anesthetics (Moran et al., 2011). This may be elucidated in the dynamics of the inhibitory subpopulation or the interaction between excitatory and inhibitory subpopulations (Gómez et al., 2013).

The adaptive two-state DCM is the neuronal connectivity component of P-DCM developed by Havlicek et al. (2015). The parameters that incorporate adaptation and refractory effects to neuronal response are: (1) the inhibitory gain factor λ, which modulates the amplitude and temporal smoothness of the inhibitory activity in relation to the excitatory drive, and (2) the inhibitory–excitatory connection μ which reflects the temporary imbalance in temporal smoothness between excitatory and inhibitory activity that can result in neuronal adaptation.

The second modification in P-DCM involves the hemodynamic model. The changes are 2-fold: (1) modeling feedforward neurovascular coupling (vs. feedback NVC in classical DCM); and (2) incorporating a viscoelastic effect in the Balloon model. Experiments by Lindauer et al. (2010) and Powers et al. (1996) show that manipulating oxygen and glucose levels in the blood do not regulate blood flow as per negative feedback hypotheses (Attwell et al., 2010; Havlicek et al., 2015). Additionally, hypercapnia experiments demonstrated that higher baseline CBF has minimal impact on the absolute stimulus-induced CBF change (while relative CBF became smaller) (Li et al., 2000; Brown et al., 2003; Zappe et al., 2008). This latter point is pertinent to DCM application in NHP since isoflurane—the most commonly used anesthetic in animal experiments—is a potent vasodilator and causes higher baseline CBF that may be uncoupled from cerebral energy metabolism (Van Aken and van Hemelrijck, 1991). Thus, feedforward neurovascular coupling is more relevant in NHP image analysis, particularly under anesthesia. On the other hand, the viscoelastic component in the hemodynamic model was added to account for transient responses (i.e., BOLD post-stimulus undershoot and overshoot) outside of the steady-state relationship between CBF and CBV, described by the power law where α = 0.38 (Grubb et al., 1974). The duration of the transient adjustment period is regulated by viscoelastic time constant τ –which allows for variations in outflow curve during balloon inflation and deflation, corresponding to BOLD response overshoot and undershoot (Buxton et al., 1998, 2004). The BOLD transients are then presumed to reflect both neuronal post-stimulus deactivation and vascular uncoupling related to slow recovery of venous CBV (Havlicek et al., 2015).

Lastly, the BOLD signal change equation was modified to accommodate different magnetic field strengths. k1, k2, and k3 are parameters that reflect baseline physiological properties of brain tissue and acquisition parameters and have been adjusted to depend on different magnetic field strengths. Moreover, they also suggest revised values for ε (ratio of intra—and extravascular signal) and r0 (regression slope of changes in intra-vascular signal relaxation rate with changes in oxygen saturation) depending on acquisition sequence (gradient echo vs. spin echo) (Uludag et al., 2009; Havlicek et al., 2015).

As stated above, we propose to combine P-DCM with stochastic DCM. Since stochastic DCM uses the same biophysical forward model as classical DCM, the changes in biophysical parameters in P-DCM can be integrated in the stochastic extension (Havlicek et al., 2015). Thus, in this application of DCM on NHP fMRI, we suggest estimating the full forward model using the biophysical parameters of P-DCM applying the generalized filtering (GF) inversion method of stochastic DCM (Li et al., 2011). A schematic of IPS-DCM is illustrated in Figure 1. A salient consideration here is that stochastic DCM tends to emphasize neuronal dynamics over hemodynamics (Reviewer 2 comment). Nonetheless, as aforementioned, stochastic DCM is a more accurate model when there are non-linear interactions among hidden states, such as the non-linear increase in GABAergic inhibitory drive (Moran et al., 2011; Daunizeau et al., 2012) which is relevant when imaging under anesthesia. For NHP neuroimaging under anesthesia, this is useful for task-based imaging with simple visual/ auditory experiments, and even more so for resting state studies. Therefore, in order to more accurately estimate both neuronal and hemodynamic changes, future extensions of P-DCM (which is so far formulated and applied as a deterministic model) to fully stochastic or partially stochastic (stochastic neuronal model and deterministic hemodynamic model) may be relevant.

Figure 1

Schematic of integrated physiologic-stochastic DCM (IPS-DCM). After identifying and extracting the time series from regions of interest, model parameters are then estimated using P-DCM equations (Havlicek et al., 2015). P-DCM incorporates: (1) an adaptive two-state neuronal model that allows adaptation and refractory effects to neuronal response; (2) a hemodynamic model that implements feedforward neurovascular coupling and a viscoelastic effect on the Balloon model; (3) a BOLD signal change equation that accounts for magnetic field differences. The model inversion is done using generalized filtering (stochastic DCM) (Li et al., 2011). Lastly, one proceeds to model comparison and selection of the winning model.

With this integrated physiologic-stochastic DCM, one can examine anesthetic impact on effective connectivity in NHP. This can be done by first building generative models of fMRI done pre and post anesthesia. Then, Bayesian model comparison allows one to examine which model best explains the changes in effective connectivity between awake and anesthetized states.

An important consideration here is that anesthetics such as isoflurane and ketamine increase baseline CBF due to increased vasodilation (Van Aken and van Hemelrijck, 1991; Långsjö et al., 2005; Slupe and Kirsch, 2018). It is also a common observation that increased baseline CBF results in smaller relative CBF change (Li et al., 2000; Cohen et al., 2002; Zappe et al., 2008). While current DCMs model relative CBF changes (Friston et al., 2003; Havlicek et al., 2015), the baseline CBF is not directly expressed. Thus, in order to account for the effect that anesthetics have on relative change in CBF, the equation modeling CBF change in DCM models could be further modified to reflect also baseline CBF values. It would be then also beneficial to consider multi-modal data that measure both CBF and BOLD signals (Havlicek et al., 2017). Even if CBF is not measured directly, it is helpful to have baseline CBF as a parameter. Additionally, as mentioned above, anesthesia can the increase latency of CBF response (also later reflected in BOLD response), e.g., from ~2 to ~4 s (Martin et al., 2006). Latency of CBF response can be controlled in DCM models. For example, Havlicek et al. (2017) accounted for differences between latency of positive and negative CBF responses measured in anesthesized NHP by optimizing parameter χ in the feedforward neurovascular coupling of the P-DCM, while other parameters (ϕ, φ) could remain fixed. These considerations could permit more physiologically accurate evaluation of the effect of pre- and post-anesthesia on effective connectivity, as modeled using IPS-DCM.

Comparison With Other DCM Applications Under Anesthesia

DCM has been applied in two other studies examining connectivity changes under propofol-induced loss of consciousness in fronto-parietal (Boly et al., 2012) and auditory (Gómez et al., 2013) networks. Boly et al. (2012) applied DCM for steady-state responses (SSR) and used neural mass models for each region of the fronto-parietal network (plus thalamic source)—three subpopulations in cortical regions (one excitatory and two inhibitor) and two for the thalamic source (excitatory relay cells and inhibitory reticular cells). On the other hand, Gómez et al. (2013) utilized combined stochastic (Daunizeau et al., 2009) and two-state (Marreiros et al., 2008) DCM to model one frontal and two temporal (auditory) cortical regions. In this section, we discuss the main caveats of both studies and how IPS-DCM may address said limitations.

The main weakness of both studies is the steady-state assumption of both DCMs which may not hold under anesthesia. Boly et al. (2012) found in initial increase in beta and gamma rhythms while delta to alpha frequencies continually increased as loss of consciousness emerged and acknowledge that this may reflect neuronal up and down states. Gómez et al. (2013) asserts that they maintained steady-state by keeping the anesthetic dose fixed after Ramsay evaluation during fMRI recordings. However, Yeom et al. (2017) demonstrated that even when patient-controlled anesthetic levels (propofol and midazolam) were unchanged once consciousness was lost, over time, there was gradually increasing power in frequencies <15 Hz together with decreasing power at >15 Hz. Increased delta and alpha power were most evident in frontal and parieto-occipital regions. Thus, the steady-state assumption of both DCMs may not hold in both studies with anesthetics.

The combined use of the adaptive two-state DCM by Havlicek et al. (2015) and stochastic DCM by Li et al. (2011) offers a number of advantages over the previous two DCM applications. Li et al. (2011) demonstrated that the generalized filtering inversion method in this DCM extension provides better effective connectivity estimates and higher sensitivity to detecting group differences than the expectation maximization (EM) or dynamic expectation maximization (DEM) of classical and variational Bayes stochastic DCM (i.e., it detected two additional connections exhibiting group differences).

As for the neuronal model, Havlicek et al. (2015) performed simulations comparing standard two-state DCM and P-DCM (adaptive two-state) and showed time courses of responses to either 1 s or 30 s stimulation in one region, as well as connectivity dynamics of a three-region network in response to a 30 s stimulation. They demonstrated that with standard two-state DCM, varying the neuronal post-stimulus deactivation does not translate to the BOLD undershoot, while P-DCM accurately models post-stimulus neuronal deactivation, both in response to 1 s or 30 s stimulation. This ability also holds even when CBV and CBF are uncoupled—the BOLD post-stimulus undershoot is stronger than the CBF response which resembles experimental results (Chen and Pike, 2009).

Pertaining to the connectivity dynamics of the simulated neuronal network, standard two-state DCM was unable to capture decreased neuronal activity below baseline—which suggests poor signal variance when both activations and deactivations are in the measured BOLD data. Meanwhile, P-DCM effectively expresses transients in neuronal and BOLD responses during faster and slower neuronal dynamics. This holds in both positive and negative responses (activation and deactivation).

Conclusion

In summary, the most important factors to consider when applying DCM in NHP under anesthesia are cerebrovascular physiology as well as anesthetic-induced changes in neural and BOLD response dynamics. Moreover, caution is necessary when interpreting results if the model includes regions or networks in which effective connectivity may be modulated by anesthetics such as: (1) fronto-parietal; (2) sensory-motor (i.e., S1–M1, SMA–M1, auditory cortical regions); and (3) thalamocortical networks specifically involving somatosensory and motor function.

The integrated P-DCM (Havlicek et al., 2015) and stochastic DCM (Li et al., 2011) (IPS-DCM) aims to address the modulatory effects of anesthetics on neural activity and the BOLD response such as changes in inhibitory drive, neural refractory period, and cortical adaptation (Masamoto et al., 2007; Moran et al., 2011). This is accomplished through: (1) an adaptive two-state neuronal model that incorporates adaptation and refractory effects to neuronal response; (2) a hemodynamic model that incorporates feedforward neurovascular coupling and a viscoelastic effect on the Balloon model; (3) a BOLD signal change equation that accounts for magnetic field differences; and (4) stochastic (generalized filtering) model inversion that addresses non-linear interactions among hidden states, such as the non-linear increase in GABAergic inhibitory drive (Moran et al., 2011; Daunizeau et al., 2012) under anesthesia.

This paper presents the motivations for applying DCM to NHP fMRI and potential strategies for addressing anesthetic effects on neuronal activity and BOLD response, which is pertinent in primate neuroimaging under anesthesia. Clearly, a number of challenges remain. For example, the validity of this physiologic-stochastic DCM integration needs to be established. Additionally, its applicability to resting state data is also an interesting consideration. Stochastic DCM has been validated in resting-state fMRI (Razi et al., 2015); however, P-DCM has only been applied to task-based recordings (Havlicek et al., 2017). Future investigations are needed to examine the face and construct validity of IPS-DCM, as well as its applicability to resting-state data.

Author Contributions

DJ and DD conceived the review focus, reviewed the literature, finalized the manuscript, and approved the final version of the manuscript. DJ summarized the literature review and wrote the first draft.

Conflict of Interest Statement

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.