Hemodynamic variability in soldiers with trauma: Implications for functional MRI connectivity studies.

Functional MRI (fMRI) is an indirect measure of neural activity as a result of the convolution of the hemodynamic response function (HRF) and latent (unmeasured) neural activity. Recent studies have shown variability of HRF across brain regions (intra-subject spatial variability) and between subjects (inter-subject variability). Ignoring this HRF variability during data analysis could impair the reliability of such fMRI results. Using whole-brain resting-state fMRI (rs-fMRI), we employed hemodynamic deconvolution to estimate voxel-wise HRF. Studying the impact of mental disorders on HRF variability, we identified HRF aberrations in soldiers (N = 87) with posttraumatic stress disorder (PTSD) and mild-traumatic brain injury (mTBI) compared to combat controls. Certain subcortical and default-mode regions were found to have significant HRF aberrations in the clinical groups. These brain regions have been previously associated with neurochemical alterations in PTSD, which are known to impact the shape of the HRF. We followed-up these findings with seed-based functional connectivity (FC) analysis using regions-of-interest (ROIs) whose HRFs differed between the groups. We found that part of the connectivity group differences reported from traditional FC analysis (no deconvolution) were attributable to HRF variability. These findings raise the question of the degree of reliability of findings from conventional rs-fMRI studies (especially in psychiatric populations like PTSD and mTBI), which are corrupted by HRF variability. We also report and discus, for the first time, voxel-level HRF alterations in PTSD and mTBI. To the best of our knowledge, this is the first study to report evidence for the impact of HRF variability on connectivity group differences. Our work has implications for rs-fMRI connectivity studies. We encourage researchers to incorporate hemodynamic deconvolution during pre-processing to minimize the impact of HRF variability.

Introduction

Functional MRI (fMRI) is used extensively for studying the neural correlates of brain functioning. FMRI is an indirect measure of neural activity as it measures changes in blood oxygenation level. Blood oxygenation is impacted by neural activity, the neurochemical signals which couple neural activity with blood flow, the properties of blood flow, and the biochemistry of blood's response to oxygen demand from the neuron (Handwerker et al., 2004). The non-neural components of the hemodynamic response vary across brain regions, which are in turn variable across individuals (Handwerker et al., 2004, Aguirre et al., 1998). With neural activity being the subject of interest in fMRI studies, interpretation of results is often less reliable due to the aforementioned non-neural sources of variability in fMRI.

The mathematical transfer function between local neural activity and corresponding blood oxygenation level dependent (BOLD) fMRI signal is called the hemodynamic response function (HRF). Most fMRI studies assume a standard canonical HRF (usually made up of two gamma functions) during analysis. However, prior research shows variability in HRF for different brain regions across subjects (Handwerker et al., 2004, Aguirre et al., 1998). This challenges the interpretation of fMRI results since it is unclear if observed changes are due to neural activity or HRF variability. There are three main dimensions of variability in HRF: (i) intra-subject spatial variability (the HRF being different in different brain regions within the same individual), (ii) inter-subject intra-group variability (for a given location in the brain, the HRF being different across different healthy individuals), and (iii) inter-group variability (for a given location in the brain, the HRF being different between a healthy group and a pathological group, arising at least in part from neurochemical disturbances due to pathophysiology). Each of these dimensions can lead to misleading results during fMRI data analysis.

Intra-subject variability, for example, could lead to detection of false activations or mistaken strong connectivities, as well as missed true activations or mistaken weak connectivities. Inter-subject variability can lead to noisy variations in variables of interest, thus reducing the statistical significance of true effects, while attributing higher significance to false effects. Inter-group variability causes detection of wrong group differences in activations or connectivities, as well as missed detection of true effects. The effect of HRF variability on activation analysis can be alleviated, in part, by using time and dispersion derivatives in the general linear model (Poldrack et al., 2011). Much attention has also been received on the effect of HRF variability on lag-based connectivity models (Deshpande et al., 2010). However, its effects on zero-lag connectivity models based on correlation measures have not been explored. In this work, we address this issue by investigating the effect of inter-group HRF variability on functional connectivity differences between the groups. In order to do so, we considered the case of soldiers with posttraumatic stress disorder (PTSD) and post-concussion syndrome (PCS) (Cicerone and Kalmar, 1995), a chronic outcome associated with mild traumatic brain injury (mTBI).

PTSD and mTBI arising from combat exposure are highly relevant to the society, with a significant percentage of soldiers acquiring them during warfare. In the U.S. alone, > 2.7 million soldiers served in Iraq and Afghanistan, with about 20% acquiring PTSD, 19% acquiring mTBI and 7% acquiring both (http://www.veteransandptsd.com/PTSD-statistics.html). PTSD has high comorbidity with mTBI (Hoge et al., 2008), added to the fact that they have similar symptomatology (Eierud et al., 2014). Recent evidences using Doppler ultrasound and infrared spectroscopy suggested alterations in cerebrovascular reactivity in mTBI (Len and Neary, 2011). Neurochemical alterations in PTSD are well established (Southwick et al., 1999), though it is important to explore if these changes affect cerebrovascular reactivity. We hypothesized that the HRF, which depends on cerebrovascular reactivity and neurovascular coupling, may be altered in soldiers with PTSD and PCS. We tested this primary hypothesis by obtaining significant group differences in voxel-specific HRF parameters, which were estimated by performing hemodynamic deconvolution of resting-state fMRI (rs-fMRI) data obtained from these populations. As a corollary, we also tested the hypothesis (= secondary hypothesis) that functional connectivity differences between groups are, at least, partially driven by HRF differences, if HRF variability is not removed through deconvolution.

The HRF is characterized by three main parameters (Handwerker et al., 2004, Aguirre et al., 1998) (see Fig. 1): (i) response height (RH), (ii) time-to-peak (TTP), and (iii) full-width at half-max (FWHM). Recent works have shown that reduced TTP and FWHM as well as increased RH are attributable to disruption in metabolism and microvasculature associated with brain pathology (Mayer et al., 2014). Additionally, aging causes reduction in TTP/FWHM and increase in RH (Arichi et al., 2012). Taken together, this profile of HRF alterations could be indicative of degraded neurochemical metabolism in the brain. With our primary hypothesis, we predicted that the HRF in the PTSD and PCS groups would be taller (RH), quicker (TTP), and narrower (FWHM) as compared to healthy combat controls in certain regions.

Fig. 1

Typical HRF with its three characteristic parameters. FWHM = full-width at half max.

Upon identifying HRF differences between the groups, it would be necessary to elucidate its impact on rs-fMRI data analysis and subsequent inferences. This would impact rs-fMRI analysis at large, as well as studies on PTSD and mTBI. In this work, we consider the impact of HRF variability on functional connectivity (FC) analysis. Specifically, we investigated the negative effects of HRF variability on group-level connectivity differences.

We illustrate two possible negative effects using example time series from experimental rs-fMRI data (see Fig. 2): (i) there exists no correlation between latent neural signals, but BOLD fMRI time series show high correlation, and (ii) there is true high correlation between latent neural signals, but BOLD fMRI time series show no correlation. The former leads to possible false positives while the latter leads to possible false negatives in traditional rs-fMRI FC analysis that does not remove HRF variability. Since the ground-truth HRF is not known, and the HRF is being estimated from a blind deconvolution approach, we do not take the leap of calling them false positives or negatives. Instead they are likely false positives and false negatives, or what we term as pseudo-positives and pseudo-negatives respectively.

Fig. 2

Illustration of the effect of HRF variability on connectivity analysis. Using two timeseries picked from real fMRI data, we demonstrate that: (a) the BOLD fMRI timeseries are highly correlated while the underlying neural signals are not (giving false high correlation when the true correlation is low), leading to likely false positives; and (b) the underlying neural signals are highly correlated while the BOLD fMRI timeseries are not (giving false low correlation when the true correlation is high), leading to likely false negatives.

In this work, we identified the cluster(s) which showed significant differences in HRF between the groups (our primary hypothesis), and performed seed-based connectivity using the cluster(s) as seed(s). Similar to the example, we hoped to identify pseudo-positives and pseudo-negatives arising from traditional connectivity analysis, which ignores HRF variability, given the fact that those seeds have different HRF profiles across the groups. Significant group differences in connectivity were obtained. This procedure was performed for two separate pipelines: (i) data pre-processed without hemodynamic deconvolution, (ii) data pre-processed with hemodynamic deconvolution. Deconvolution minimizes HRF variability in BOLD fMRI, giving latent neural variables. We then compared the group differences in connectivity for the two pipelines. We hypothesized that (= secondary hypothesis), owing to HRF variability, data without deconvolution would show misleading connectivity differences (both pseudo-positives and pseudo-negatives) as compared to the data with deconvolution.

Methods

Participants

Active-duty male U.S. Army soldiers between the ages of 18 and 50 years were recruited (N = 87) from Fort Rucker, AL, USA and Fort Benning, GA, USA to voluntarily participate in the study. Participants were grouped into 17 with PTSD, 42 with both PCS and PTSD (PCS+PTSD), and 28 being combat controls (groups matched in race, education and age), all having combat experience in Iraq (Operation Iraqi Freedom, OIF) and/or Afghanistan (Operation Enduring Freedom, OEF). The study was conducted in accordance with the Declaration of Helsinki (latest version). The procedures and protocol were approved by the Headquarters U.S. Army Medical Research and Materiel Command, IRB (HQ USAMRMC IRB) and Auburn University's Institutional Review Board.

Subject grouping was done based on PTSD symptom severity using the PTSD Checklist-5 (PCL5) score, clinician referral, post-concussive symptoms using the Neurobehavioral Symptom Inventory (NSI) score and medical history. Each participant was clinically assessed in the clinics at Fort Benning, GA and Fort Rucker, AL (where the participants were recruited from), and were either identified as healthy soldiers or were diagnosed with either PTSD, PCS, or both by a licensed medical practitioner. The PCL-5 and NSI self-report questionnaires were administered to the participants when they arrived at the Auburn University MRI Research Center for scanning. (i) Subjects without any history of mTBI in prior five years, clinician referral, a PCL5 score ≥ 38 and NSI score < 26 were classified as PTSD. (ii) Subjects with history of mTBI, post-concussive symptoms, clinician referral and PCL5 score ≥ 38 and NSI score ≥ 26 were classified as PCS+PTSD (comorbid group). (iii) Subjects with a PCL5 score < 38 and NSI score < 26, no mTBI within the previous 5 years, no DSM-V or DSM-IV-TR diagnosis of a psychiatric disorder (based on medical records), and no history of a moderate to severe TBI were classified as combat controls upon clinical evaluation. Participants with mood, psychotic or substance dependency disorders were excluded. All participants reported being deployed in a combat environment. PCL5 scores were statistically significant (F(1, 172) = 20.6443, p = 3.64 × 10− 44) between control group and the PTSD and PCS+PTSD groups combined. We performed such a comparison because the common factor between PCS+PTSD and PTSD groups is PTSD, while only PTSD symptom severity is reflected in the PCL5 score. Similarly, NSI scores were found to be statistically significant (F(1, 172) = 32.6878, p = 1.32 × 10− 29) between PCS+PTSD group, and the PTSD and control groups combined. Fig. 3 provides a flowchart illustrating the classification of subjects into different groups.

Fig. 3

Flowchart illustrating the selection criteria of subjects and their classification into the three groups: PTSD, comorbid PCS and PTSD (PCS+PTSD) and healthy combat controls.

Measures

PTSD Checklist-5 (PCL5 (Dickstein et al., 2014))

PCL5 is a self-report measure which assesses DSM-5 symptoms of PTSD. It screens individuals for PTSD, helps make PTSD diagnoses and aids in monitoring change in symptoms during and after treatment. Items are rated on a 5-point Likert scale, with 1 being “Not at all“ to 5 being “Extremely“. With 20 items, a total score is obtained in the range 20–100 by summing the scores of each of the 20 items, with a cut-score of 38 for the diagnosis of PTSD (Weathers et al., n.d.).

Neurobehavioral Symptom Inventory (NSI (Cicerone and Kalmar, 1995))

NSI is a self-report questionnaire which assesses post-concussive symptoms in persons who have experienced a TBI. For each symptom, participants rate its severity within the past month on a 5-point Likert scale, which ranges from 0 (none) to 4 (very severe). With 22 items, a total score is obtained in the range 0–88 by summing the scores of the 22 items.

Procedures

Upon arriving at the research lab, participants went through re-screening for eligibility, screening for MRI contraindications and re-consenting to guarantee full comprehension of the study's benefits, procedures and their rights.

fMRI

A 3 T MAGNETOM Verio scanner was used (Siemens Healthcare, Erlangen, Germany). Participants were scanned using T2* weighted multiband echo-planar imaging (EPI) sequence in resting-state. They were requested to have their eyes open and fixated on a white cross displayed on a dark background on the display, using an Avotec projection system. They were asked to not think of anything specific. Scanning parameters were as follows: TR = 600 ms, TE = 30 ms, FA = 55°, multiband factor = 2, voxel size = 3 × 3 × 4 mm3 and 1000 volumes. Brain coverage was restricted to the cerebral cortex, midbrain, pons and subcortical structures (cerebellum excluded). Two separate scans were done for every subject (same day), hence providing 174 sessions of rs-fMRI data for 87 subjects. This mathematically boosted the statistical power of our analysis.

FMRI data pre-processing

Standard rs-fMRI pre-processing sequence was carried out, which included realignment, normalization of rs-fMRI volumes to MNI space, detrending, and regressing out nuisance covariables (six head motion parameters, white matter and cerebrospinal fluid signals). Spatial smoothing was not performed. Pre-processing was carried out in Data Processing Assistant for Resting-State fMRI (DPARSF v1.7) (Chao-Gan and Yu-Feng, 2010), which is based on Statistical Parametric Mapping (SPM8) (Friston et al., 2007) and Resting State fMRI Data Analysis (REST) Toolkit (Song et al., 2011).

The 3D+time data were deconvolved across time at every voxel to get latent neural variables using a popular method proposed by Wu et al. (Wu et al., 2013). This method has gained increasing popularity and acceptance, thanks to its interpretability, robustness, simplicity of implementation, validity, and a rising awareness in the research community on the necessity for deconvolution. Several recent works have employed it (for example, see (Boly et al., 2015, Lamichhane et al., 2014, Amico et al., 2014)). This deconvolution is blind since we have access to only one variable (rs-fMRI timeseries), using which it estimates both the HRF and the latent neural timeseries. In short, the method models rs-fMRI data as event-related timeseries with randomly occurring events employing point processes (Saad et al., 2012, Power et al., 2015), and then evaluating voxel-wise HRFs using Weiner deconvolution. All data analysis was performed on the Matlab® platform (version R2014a).

HRF analysis

Deconvolution provided the estimated HRF at each voxel in each subject, which was characterized by three parameters – response height (RH), time-to-peak (TTP), and full-width at half-max (FWHM), as illustrated earlier. The voxel-wise HRF parameters for all subjects have been made publicly available (Rangaprakash et al., 2017a). Tests for statistical significance were performed separately on each of the three parameters to obtain group-wise voxel-specific differences in HRF parameters (p < 0.05, cluster-level thresholded, controlled for age, race, education and head-motion). This was done separately for the three pairwise comparisons between groups, that is, Control vs PTSD, Control vs PCS+PTSD and PTSD vs PCS+PTSD. Thus, we obtained nine maps: three pairwise comparisons, with three HRF parameters for each of these comparisons.

These maps were then used to obtain the following final overlapped (intersection) maps of interest. (i) O-1: Overlap between Control vs PTSD and Control vs PCS+PTSD maps (we call this control vs disease comparison), for each of the three HRF parameters separately. This would elucidate HRF differences in disease compared to healthy controls, which would directly validate or invalidate our primary hypothesis. (ii) O-2: Overlap between all the three HRF parameters' maps obtained from both Control vs PTSD and Control vs PCS+PTSD comparisons (i.e. control vs disease). These differences would obviously be a subset of the first case, but would identify those regions with alteration of all HRF parameters in PTSD and PCS. (iii) O-3: Overlap between all three group-wise comparisons for each of the three HRF parameters separately. This would illustrate those HRF differences which were altered in all the groups. These differences would naturally be a subset of the first case. These differences represent important information regarding PTSD and PCS since identifying commonalities and differences between them are of deep interest, given that the two disorders have high comorbidity and similar symptomatology (Eierud et al., 2014). In all cases, the overlapped maps were obtained by finding the common regions in the maps being considered (i.e. intersection), which had a cluster size of at least 50 mm3 in the overlapped map (to eliminate false positives).

Seed-based functional connectivity analysis

The overlapped maps obtained from the second and third cases mentioned above were used in further seed-based functional connectivity (FC) analysis, with the identified cluster(s) being chosen as the regions of interest (ROIs). The time series from all the voxels in each ROI were averaged to obtain a single time series per ROI. ROI time series were orthogonalized with respect to each other (Di Martino et al., 2008, Margulies et al., 2007). Seed-based connectivity was performed by evaluating Pearson's correlation coefficient between the ROI time series and rest of the voxels in the brain. Significant group differences in whole-brain FC maps were obtained between all groups using multivariate-ANOVA (p < 0.05, FDR corrected, controlled for age, race, education and mean head-motion). This pipeline was implemented separately for two cases: (i) NDC (no deconvolution): data pre-processed without deconvolution, and (ii) DC: data pre-processed with deconvolution. The NDC case, which is the traditional approach in most studies, contains data contaminated by HRF variability. As hypothesized by us (secondary hypothesis), we expected to see pseudo-positives and pseudo-negatives in the connectivity map obtained from NDC data as compared to the DC data. The pseudo-positives were obtained from final NDC and DC maps as “NDC > DC”. Similarly, pseudo-negatives were obtained as “NDC < DC”. We tested our secondary hypothesis through these two maps.

Results

Demographics

Group-wise demographics are presented in Table 1. There were no significant differences between the groups in age, p = 0.699, or education, p = 0.152. The results indicated that there was a difference in the frequency of reported psychotropic use between the groups, τb = 0.24, p = 0.011, with the comorbid group having the highest percentage of medicated subjects. There was a significant difference between the groups in the number of reported lifetime mTBIs, F(2) = 5.81, p = 0.004, specifically between the control group and PCS+PTSD group, but not the PTSD and PCS+PTSD groups or the control and PTSD groups, p > 0.05.

Table 1

Basic demographics.

Variable		Controls	PTSD	PCS+PTSD
Age, years	Mean	32.6	32.2	33.7
	Median	31	32	33
	SD	6.7	7.6	6.8
	Range	24	24	30
Race	White	18 (66.7%)	11 (64.7%)	26 (66.7%)
	Black	2 (7.4%)	3 (17.6%)	9 (22.0%)
	Hispanic	3 (11.1%)	3 (17.6%)	2 (4.9%)
	Asian	2 (7.4%)	0	1 (2.4%)
	Other	0	0	1 (2.4%)
Education, years	Mean	15.1	14.5	14.1
	Median	16	14	14
	SD	1.9	2.2	1.9
	Range	8	9	8
Lifetime mTBIs	Mean (range)	0.3 (2)	1.1 (6)	2.5 (15)⁎
Medication		2 (7.4%)	4 (23.5%)	13 (31.7%)⁎

Statistically significant (p < 0.05, Bonferroni corrected).

Inter-group HRF differences

As mentioned earlier, group differences in voxel-wise values of each of the three HRF parameters were obtained, and three categories of overlapped maps were derived from them. In all the regions with altered HRF, we found that RH increased in the disease groups compared to controls, while TTP and FWHM decreased in the disease groups. We first elucidate the differences for Control vs Disease comparison, which refers to an overlap of Control vs PTSD and Control vs PCS+PTSD comparisons (overlap O-1 mentioned earlier). Differences in RH were found in (see Fig. 4a) the thalamus, midbrain, precuneus, posterior cingulate cortex (PCC), secondary visual areas and parts of insula (anterior and posterior). Further, differences in FWHM (Fig. 4b) and TTP (Fig. 4c) were largely similar, with key default-mode network (DMN) regions being disrupted (PCC and precuneus) along with secondary visual areas.

Fig. 4

(a) Regions with significantly altered response height (RH) of the hemodynamic response function, HRF. They were significant for PTSD > Control and PCS+PTSD > Control comparisons. Thalamus, midbrain, insula, visual and default-mode network regions were altered. PCC = posterior-cingulate cortex. Please refer to Supplemental Table S1 for further details. (b) Regions with significantly altered FWHM in HRF. They were significant for Control > PTSD and Control > PCS+PTSD comparisons. Visual and default-mode network regions were altered. PCC = Posterior-cingulate cortex; IPL = inferior parietal lobule (angular gyrus). Please refer to Supplemental Table S2 for further details. (c) Regions with significantly altered time-to-peak in HRF. They were significant for Control > PTSD and Control > PCS+PTSD comparisons. Visual and default-mode network regions were altered. PCC = Posterior-cingulate cortex. Please refer to Supplemental Table S3 for further details. (d) Regions which had significant alterations in all three HRF parameters. They were significant for Control > PTSD and Control > PCS+PTSD comparisons. Posterior-cingulate cortex (PCC) and precuneus were identified. Please refer to Supplemental Table S4 for further details. (e) Regions which were significantly different between all three groups, implying that both PTSD and mTBI caused alterations in them. This difference was observed only with FWHM. Posterior-cingulate cortex (PCC) and precuneus were identified. Please refer to Supplemental Table S5 for further details.

Next, we present the results for overlap between RH, TTP and FWHM for Control vs Disease (overlap O-2 mentioned earlier). We found the common regions to be the left PCC and right precuneus (see Fig. 4d). Then, identifying those differences which were significantly different between all three groups (overlap O-3), we found FWHM to be significantly different between all three groups in the left PCC and right precuneus again (see Fig. 4e, note: RH and TTP did not show differences). Please refer to Supplemental Tables S1 through S5 for further details about region centroids, volumes and statistical values.

In accordance with our secondary hypothesis, we performed seed-based FC analysis using the left PCC (as well as the right precuneus, independently) as the seeds, thus obtaining two separate sets of connectivity maps, one for each seed. Group differences in FC were then obtained for each seed. We performed this procedure for two separate pipelines: (i) NDC (no deconvolution): data pre-processed without deconvolution, and (ii) DC: data pre-processed with deconvolution. Results from the two pipelines were then compared, with the hope of identifying pseudo-positives (connectivities which were not significantly different [or considerably less statistically significant] after removal of HRF variability, but were significantly different in NDC case), and pseudo-negatives (connectivities significantly different in DC case but not (or less) in NDC case). Pseudo-positives were obtained as: NDC map > DC map, and pseudo-negatives as: NDC map < DC map. Here we present the results for the PCC ROI, while results for the precuneus ROI can be found in Supplemental Information since the latter also leads to the same conclusion as the former (Figs. S1 through S4 and Tables S10 through S13).

Group differences for the NDC (Fig. 5a) and DC (Fig. 5b) pipelines clearly show that the identified significant functional connectivity alterations in PTSD and PCS groups differ appreciably between the two pipelines. This translates to providing us the evidence that HRF variability drives a sizeable portion of group differences reported in rs-fMRI studies ignoring such HRF variability. We also notice that there are more number of pseudo-positives (Fig. 5c) than pseudo-negatives (Fig. 5d), which is especially undesirable since false positives are more detrimental than missed detections. These observations are formally presented in Table 2.

Fig. 5

(a) Brain regions whose functional connectivity with the left posterior cingulate (L_PCC) seed ROI (marked blue region) was significantly different between the groups for data without hemodynamic deconvolution. Please refer to Supplemental Table S6 for further details. (b) Brain regions whose functional connectivity with the left posterior cingulate (L_PCC) seed ROI (marked blue region) was significantly different between the groups for data with hemodynamic deconvolution. Please refer to Supplemental Table S7 for further details. (c) Pseudo-positives, that is, functional connectivity group differences which were greater (higher T-value) in data without deconvolution as compared to that with deconvolution performed (for the left posterior cingulate [L_PCC] seed ROI). Please refer to Supplemental Table S8 for further details. (d) Pseudo-negatives, that is, functional connectivity group differences which were smaller (lower T-value) in data without deconvolution as compared to that with deconvolution performed (for the left posterior cingulate [L_PCC] seed ROI). Please refer to Supplemental Table S9 for further details. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 2

Number of non-zero voxels in the thresholded T-maps and its derivatives (please refer to Fig. 4 for L_PCC ROI and Figs. S1 through S4 for R_Prec ROI). There are notable differences in the significance maps (as seen in the figures) as well as corresponding notable differences in the number of voxels between non-deconvolved and deconvolved data, and their comparisons. Further details are available in Supplemental Tables S6 through S13.

Seed ROI	Number of significant voxels		Number of voxels
	Data without deconvolution	Data with deconvolution	Without deconv > with deconv	With deconv > without deconv
L_PCC	11599	5141	9968	2526
R_Prec	6407	5828	4753	3814

Discussion

In this work, we tested two hypotheses. First, we hypothesized that the HRF, which depends on cerebrovascular reactivity and neurovascular coupling, may be altered in soldiers with PTSD and PCS. We then tested this primary hypothesis by obtaining significant group differences in voxel-specific HRF parameters which were estimated by performing blind hemodynamic deconvolution of rs-fMRI data obtained from these populations. Second, we also tested the hypothesis that functional connectivity differences between groups are at least partially driven by HRF differences, if HRF variability is not removed through deconvolution.

We found substantial evidence to support our hypotheses. First, we found altered HRF parameters in the disease groups when combined compared to controls in subcortical and DMN regions. We also found a subset of these alterations to be significantly different between all three groups, which characterized HRF differences between PTSD and PCS+PTSD. In all the regions with altered HRF, we found that RH increased in the diseased groups compared to controls, while TTP and FWHM decreased in the diseased groups. Throughout, “control vs diseased” refers to an overlap of control vs PTSD and control vs PCS+PTSD comparisons. This finding conforms to our prediction made earlier that HRF in PTSD and PCS would be taller, quicker and narrower compared to controls. This profile of HRF alteration has been attributed to disrupted metabolism and microvasculature associated with brain disorders, for example in mTBI (Mayer et al., 2014). Additionally, in accordance with our findings, Thompson et al. (Thompson et al., 2014) showed that there is a negative relationship between RH and TTP/FWHM, that is, whenever the height of HRF increases, it is highly likely that its ascent and descent are quicker.

In the current study, we found RH differences between controls versus diseased mainly in the thalamus, midbrain, default-mode regions (PCC, precuneus) and secondary visual areas. Alterations in TTP and FWHM, were found in default-mode regions (PCC, precuneus) and secondary visual areas. Prior work has shown abnormal GABAergic and glutamatergic neurotransmitter systems implicated in anxiety disorders such as PTSD (Mifflin et al., 2015). Neuromodulators released by glutamatergic and GABAergic interneurons are known to directly modulate local cerebral blood flow (Buzsáki et al., 2007), and ultimately the HRF (Brown et al., 2003). Glutamate acts on N-methyl-d-aspartate (NMDA) receptors, which causes dilation of blood vessels associated with activated brain regions (Busija et al., 2007), ultimately impacting the HRF. Importantly, lower gamma-Amino butyric acid (GABA) concentration has been shown to result in taller, quicker and narrower HRFs (Muthukumaraswamy et al., 2012), which was also observed by us in a population of soldiers with PTSD and PCS. As such, at least part of the HRF alterations observed in the current study could be attributed to lower GABA concentration (as also suggested in (Muthukumaraswamy et al., 2012)).

The thalamus, which is necessary to produce anxiety, showed altered RH in our diseased groups. In fact, the thalamus is anatomically well situated for mediating pre-frontal activity with subcortical and midbrain structures to produce the experience of anxiety (Cohen, 2009). Additionally, serotonin (5-hydroxytryptamine) in the midbrain is known to play a key role in anxiety disorders (Nikolaus et al., 2010). Serotonin is a vasoconstrictor which provides blood-brain barrier permeability for modulating neurovascular coupling, and thus the HRF via the neuronal-astrocytic-vascular tripartite functional unit (Cohen et al., 1996).

Earlier studies have reported neurochemical alterations in key areas of the DMN and insula to be associated with PTSD. Using PET imaging, Ramage et al. (Ramage et al., 2015) found heightened glucose metabolism in precuneus in PTSD compared to combat controls, which correlated with PTSD symptom severity. Using magnetic resonance spectroscopy, Rosso et al. (Rosso et al., 2014) found decreased GABA in the insula in individuals with PTSD, which correlated with anxiety levels. In addition to our findings, this evidence provides substantiation for neurochemical and vascular alterations, and thus HRF differences, in brain regions found to be significantly different between healthy controls and our PTSD groups.

Interestingly, findings of both higher brain activation (Eierud et al., 2014) and hyper-connectivity (Simmons and Matthews, 2012, Cisler et al., 2014, Rangaprakash et al., 2015, Rangaprakash et al., 2017b) in PTSD and mTBI, reflecting reduced inhibition, also corroborates with the regions identified with altered HRF in this work. In light of the prior literature discussed above, reduced GABA likely causes both reduced neural inhibition and altered HRF; and the HRF alterations identified in this work largely overlapped with neural alterations identified in prior works (PCC, precuneus, secondary visual and thalamus). This observation conforms to prior studies that attribute part of HRF variability to neural activity differences (Handwerker et al., 2004), which as noted earlier might be an indirect relationship mediated by neurotransmitters like GABA (Muthukumaraswamy et al., 2012). These findings raise important considerations - if group-differences in HRF parameters and neural activity are indeed largely similar, owing to underlying neurochemistry, then interpretation of findings from rs-fMRI studies which ignore HRF variability would not be straightforward.

Next, for the overlap between RH, TTP and FWHM for Control vs Disease (overlap O-2 mentioned earlier), we identified the common regions to be the left PCC and right precuneus. Additionally, significant differences between all the three groups (overlap O-3) were also found in the left PCC and right precuneus (with FWHM). This shows that the hemodynamic response in the PCC and precuneus were possibly affected by both PTSD and mTBI. These regions were naturally a subset of the regions identified in overlap O-1. Previous studies have reported neurochemical alterations in these key areas in soldiers affected by trauma (Ramage et al., 2015, Rosso et al., 2014). These regions have also been largely implicated in PTSD (Simmons and Matthews, 2012) and mTBI (Eierud et al., 2014) studies. This is a substantial finding given that neural underpinnings of comorbid PTSD and mTBI are poorly understood (Simmons and Matthews, 2012). The PCC and precuneus showed altered HRF between all three groups with the three HRF parameters.

We also found strong evidence in support of our secondary hypothesis. Seed-based functional connectivity analysis with the PCC and precuneus ROIs revealed perceptible distinction in group differences obtained from data without deconvolution as compared to data with deconvolution. This implies that part of the functional connectivity group differences reported in rs-fMRI studies (which do not perform deconvolution) could potentially be attributable to non-neural HRF variability. Several previous works have speculated on this aspect (Handwerker et al., 2004, Aguirre et al., 1998, Deshpande et al., 2010), since BOLD fMRI is not a direct measure of neural activity. We provide quantitative formal evidence for the impact of hemodynamic variability on group differences in rs-fMRI-derived measures like functional connectivity. In view of this, we urge the community to employ deconvolution in their pre-processing pipeline to remove/reduce HRF variability. Researchers must exercise caution in interpreting their results when the effect of HRF variability has not been accounted for in their resting-state functional connectivity analysis, more so when dealing with PTSD and mTBI populations. Future studies could also investigate the effect of HRF variability on functional connectivity measures in other clinical populations.

To summarize the conclusions of this paper in one simple sentence: HRF differences were found to exist between healthy and our clinical groups in certain brain regions, and such differences, if unaccounted for, could potentially drive false connectivity findings. Evidently, our study has both clinical and methodological conclusions which are related to each other. This work achieves two distinct but connected objectives: (i) identifying HRF alterations in PTSD and comorbid PCS+PTSD compared to healthy soldiers, which is a consequence of the altered neurochemistry, vasculature and neurovascular coupling in the identified brain regions, and (ii) the impact of such HRF differences between the groups on functional connectivity modeling. The former and latter objectives are connected in that the former finds evidence for HRF alterations in mental disorders like PTSD and PCS+PTSD, and the latter finds evidence for its consequence on connectivity modeling. As such, this paper contributes both clinically and methodologically, and the conclusions drawn are thus also twofold: (i) certain brain regions identified in this study exhibit HRF differences between PTSD, PCS+PTSD and healthy soldiers, which has implications for understanding these disorders, and (ii) if functional connectivity modeling were to be performed without accounting for such HRF variability then a portion of the connectivity findings (especially those associated with the regions identified in this study, but not limited to it) could be potentially false findings in the sense that they may not be entirely neural in origin, which has implications for rs-fMRI data analysis that employs functional connectivity modeling, specifically to study these disorders as well as to study mental disorders in general. Large number of fMRI studies often report a new method, apply it on clinical data and draw clinical conclusions; while in this paper we reported clinical findings, and followed it up with a methodological conclusion that might impact the way in which connectivity modeling is performed in future studies.

In conclusion, we showed that PTSD and PCS are associated with overlapping and distinct HRF alterations in subcortical structures and the DMN. The HRF is known to vary with the balance between excitatory and inhibitory neurotransmitters, with differences in size of neighboring vasculature and with neurovascular coupling (Handwerker et al., 2004); and given that these cerebral characteristics are known to be altered in mental conditions, we found evidence for HRF alterations in PTSD and PCS in this work. Our findings also corroborate with prior findings, in addition to providing new insights and directions. Since our findings were obtained from an overlap/intersection of results for the PTSD and the PCS+PTSD groups, the observations and conclusions are equally applicable to the study of PTSD alone.

Given these findings, future studies on PTSD and PCS, and rs-fMRI studies in general, must practice prudence in reporting and interpreting results obtained from resting-state functional connectivity analysis of non-deconvolved BOLD fMRI data; especially if they assume a fixed canonical HRF. Although only a handful of regions showed HRF alterations in this work (since we used a conservative statistical threshold which might have ignored smaller effects), one cannot rule out that HRF variability between groups does not exist elsewhere. Based on our findings, we encourage researchers to employ hemodynamic deconvolution to better understand the neural basis of pathophysiological differences between various disease states.