Counteracting Effects of Glutathione on the Glutamate-Driven Excitation/Inhibition Imbalance in First-Episode Schizophrenia: A 7T MRS and Dynamic Causal Modeling Study

Oxidative stress plays a key role in the pathophysiology of schizophrenia. While free radicals produced by glutamatergic excess and oxidative metabolism have damaging effects on brain tissue, antioxidants such as glutathione (GSH) counteract these effects. The interaction between glutamate (GLU) and GSH is centered on N-Methyl-D-aspartate (NMDA) receptors. GSH levels increase during glutamate-mediated excitatory neuronal activity, which serves as a checkpoint to protect neurons from oxidative damage and reduce excitatory overdrive. We studied the possible influence of GSH on the glutamate-mediated dysconnectivity in 19 first-episode schizophrenia (FES) patients and 20 healthy control (HC) subjects. Using ultra-high field (7 Tesla) magnetic resonance spectroscopy (MRS) and resting state functional magnetic resonance imaging (fMRI), we measured GSH and GLU levels in the dorsal anterior cingulate cortex (dACC) and blood-oxygenation level-dependent activity in both the dACC and the anterior insula (AI). Using spectral dynamic causal modeling, we found that when compared to HCs, in FES patients inhibitory activity within the dACC decreased with GLU levels whereas inhibitory activity in both the dACC and AI increased with GSH levels. Our model explains how higher levels of GSH can reverse the downstream pathophysiological effects of a hyperglutamatergic state in FES. This provides an initial insight into the possible mechanistic effect of antioxidant system on the excitatory overdrive in the salience network (dACC-AI).


Introduction
The role of oxidative stress in molecular mechanisms of neurodegenerative diseases has been investigated thoroughly [1]; however, the details of its involvement in pathophysiology of schizophrenia are not completely understood [2]. While free radicals produced by glutamatergic excess and oxidative metabolism have damaging effects on brain tissue, antioxidants such as glutathione (GSH) counteract the "toxic effect" of oxidative stress [3]. This antioxidant protection is reportedly aberrant in schizophrenia, with some patients having a notable reduction of GSH in the brain [4][5][6] and others with relatively better outcomes having higher than expected levels [7]. Several pharmacological approaches to improve glutathione-mediated antioxidant capacity are currently being studied in schizophrenia [8].
The interaction between glutamate (GLU) and glutathione is centered in the N-Methyl-D-aspartate (NMDA) receptors as well as the neuro-glial metabolic shuttling. GSH levels increase in response to glutamate-mediated excitatory neuronal activity. This increase is mediated by the NMDA receptor system [9] as well as by the conversion of glutamate to GSH in glial cells [10]. GSH has also been reported to have a direct signaling effect, by facilitating NMDA function [11,12] and increasing the inhibitory tone of microcircuits [13]. In humans, correlations between GSH and GLU levels in the dorsal anterior cingulate cortex (dACC) and the anterior insula (AI) have been reported [7]. Taken together, a glutamate-mediated GSH increase serves as a checkpoint to protect neurons from oxidative damage and reduce excitatory overdrive.
GSH depletion in the developing brain affects the interneurons that provide the NMDA-mediated inhibitory checkpoint for glutamatergic activity. This early NMDA hypofunction is considered to disrupt the normal excitation-inhibition balance and prime cortical networks for glutamatergic excess and excitatory overdrive in schizophrenia [7]. At present, it is unclear if higher levels of GSH can overcome this excitatory overdrive by restoring the inhibitory tone in patients with schizophrenia.
We have recently provided the first imaging evidence for the NMDA hypofunction model by demonstrating that glutamate levels are indeed related to a reduced inhibitory tone in schizophrenia [14]. We studied glutamate levels from the dorsal anterior cingulate cortex (dACC) and studied the salience network that connects dACC with the anterior insula (AI), using a biologically realistic neural model of resting-state functional magnetic resonance imaging (fMRI) data. In the current work, we extended this observation to study if GSH influences the glutamate-mediated dysconnectivity in first-episode schizophrenia (FES). Such an influence, if demonstrated, will add credibility to the notion that GSH can physiologically counteract the glutamate-mediated excitation-inhibition imbalance in schizophrenia.

Magnetic Resonance Spectroscopy (MRS) Acquisition and Analysis
All data was acquired using a 680-mm neuro-optimized 7 T MRI scanner (Siemens MAGNETOM Plus, Erlangen, Germany) equipped with an AC84 II head gradient coil and an 8-channel Tx, 32-channel Rx radiofrequency coil. We defined a 2.0 × 2.0 × 2.0 cm (8 cm 3 ) 1 H-MRS voxel on the bilateral dACC ( Figure 1). To this aim, we used a two-dimensional sagittal anatomical image (37 slices, TR = 8000 ms, TE = 70 ms, flip-angle (α) = 120 • , thickness = 3.5 mm, field of view = 240 × 191 mm) as reference. We defined the voxel position both by setting the posterior face of the voxel in coincidence with the precentral gyrus and by setting the position of the inferior face of the voxel to the most caudal point not part of the corpus callosum. We set the voxel angle tangentially to the corpus callosum. A semi-LASER 1 H-MRS sequence (TR = 7500 ms, TE = 100 ms, bandwidth = 6000 Hz, N avg = 2048) was used to acquire 32 channel-combined, [17] VAPOR water-suppressed spectra as well as a water-unsuppressed spectrum (N avg = 1) to be used for spectral editing and quantification. We asked all participants to fix their gaze on a white cross (50% gray background) during MRS acquisition. All scanning took place at the Centre for Functional and Metabolic Mapping of Western University, London, Ontario. LASER 1 H-MRS sequence (TR = 7500 ms, TE = 100 ms, bandwidth = 6000 Hz, Navg = 2048) was used to acquire 32 channel-combined, [17] VAPOR water-suppressed spectra as well as a water-unsuppressed spectrum (Navg = 1) to be used for spectral editing and quantification. We asked all participants to fix their gaze on a white cross (50% gray background) during MRS acquisition. All scanning took place at the Centre for Functional and Metabolic Mapping of Western University, London, Ontario. Based on Near and colleagues [18], we phase-and frequency-corrected the 32 spectra. Following, we computed a single average spectrum which was used in all subsequent analyses. Spectrum's line shape deconvolution and removal of a residual water signal was Based on Near and colleagues [18], we phase-and frequency-corrected the 32 spectra. Following, we computed a single average spectrum which was used in all subsequent analyses. Spectrum's line shape deconvolution and removal of a residual water signal was performed via QUECC [19] (a combination of quantification improvement by converting lineshapes to the Lorentzian type, QUALITY, and eddy current correction, ECC) and Han-kel singular value decomposition (HSVD) [20], respectively. Spectral fitting was done via fitMAN [21] (a time-domain fitting algorithm that uses a non-linear, iterative Levenberg-Marquardt minimization algorithm and echo-time, field strength, and pulse sequence specific prior knowledge templates). The metabolite-fitting template included 17 brain metabolites including glutamate and glutathione reported here. The other metabolites were N-acetyl aspartate, N-acetyl aspartyl glutamate, alanine, aspartate, choline, creatine, γ-aminobutyric acid (GABA), glucose, glutamine, glycine, lactate, myo-inositol, phosphorylethanolamine, scyllo-inositol, and taurine. Due to the long echo time used, no significant macromolecular contribution was expected. Metabolite quantification was then performed using Barstool [22] with corrections made for tissue-specific (gray matter, white matter, CSF) T 1 and T 2 relaxation through partial volume segmentation calculations of voxels mapped onto T 1 -weighted images acquired using a 0.75-mm isotropic MP2RAGE sequence (TR = 6000 ms, TI 1 = 800 ms, TI 2 = 2700 ms, flip-angle 1 (α 1 ) = 4 • , flip-angle 2 (α 2 ) = 5 • , FOV = 350 × 263 × 350 mm, T acq = 9 min 38 s, iPAT PE = 3 and 6/8 partial k-space). All spectral fits underwent visual quality inspection as well as Cramer-Rao lower bounds (CRLB) assessment for each metabolite.
The quality of metabolite quantification was measured using CRLB percentages for both groups using a CRLB threshold < 30% for glutathione to determine inclusion toward further analyses, in line with our prior study [23]. There was no significant difference in CRLB between the FES patients and HC subjects for both metabolites being reported in this study. A sample of fitted spectrum for a single participant is presented in Figure 1.

Bayesian Analysis
We estimated the posterior distribution of the (estimated) between-group differences in GLU and GSH by means of the generalized linear model within the context of hierarchical the Bayesian parameter estimation as follows: where the data conformed to a normal distribution around the predicted value (metabolite concentration) with a (wide) data-scaled uniform prior distribution for the standard deviation (σ i ). The baseline parameter (β 0 ) had a data-scaled normal prior distribution with mean equal to the data mean and (wide) standard deviation relative to the standard deviation (SD data ) of the data (1/(SD data × 5) 2 ). Group deflection parameters (β group ) had normal prior distributions with mean zero and a Gamma prior distribution for the standard deviation σ β with data-scaled shape and rate parameters (SD data /2 and 2 × SD data respectively). This meant that σ β provided informed priors on each group's (deflection) parameter. In other words, groups would act as priors between each other. In total, we estimated posterior distributions of five free parameters (σ i , β 0 , β HC , β FES , and σ β ). Posteriors were estimated in the R-software equivalent of "just another Gibbs sampler" (RJAGS) [24] using Markov chain Monte Carlo methods, drawing 11,000 samples (thinning = 10). We reported the proportion of the posterior distribution (i.e., posterior proportion, PP) of the between-groups difference in GSH and GLU levels along with the 95% highest density interval (HDI) of the posterior proportion. The posterior distributions and HDIs of the relevant effect sizes were also reported.

Spectral Dynamic Causal Modeling of Network Connectivity
We fit spectral dynamic causal models to the fMRI time series data to quantitatively infer how the fMRI timeseries were generated by (unobserved) neural activity of coupled neuronal populations between the dACC and the AI during resting state [25]. At present, dynamic causal modeling is considered the most physiologically grounded technique to infer the effective connectivity between brain regions [26]. Specifically, a spectral dynamic causal model is a special case of "generative models" in which the neural causes are hidden states and the blood-oxygenation level-dependent (BOLD) signals are observed measurements. Therefore, the (generative) dynamic causal model comprised one evolution function (2) where x'(t) is the rate of change of the neuronal states x(t), θ represents the unknown parameters of the effective connectivity, and v(t) represents the states noise. The output of the evolution function, x(t), was mapped onto an observed function (3) where y(t) is the measured BOLD signal, ϕ represents the unknown parameters, and e(t) is the observation noise. Crucially, the diffusion (or noise) terms in (2) and (3) could be parameterized. Therefore, the evolution function became a random differential equation. For a thorough mathematical description of the generative model we refer the interested reader to both Friston, et al. [27] and Razi, Kahan, Rees and Friston [25].
Our dynamic causal model of the dACC-AI network represented both intrinsic or within-region (GABAergic) connections and extrinsic (between-region) glutamatergic neuronal populations within each region [28,29]. Each population comprised self-inhibition connections (which are fixed parameters). Two free parameters were fit to the fMRI data: Interregional excitatory-to-excitatory connections and within-region inhibitory-to-excitatory connections ( Figure 2). Each of these parameters was the log of a scaling factor, which was multiplied by the default connection strength: 1/8Hz for between-region connections and −1/8Hz for within-region connections. This formulation enforces positivity or negativity constraints on the connections, and gave the parameters a simple interpretation, as follows. Between-region connections were excitatory, so more positive values corresponded to greater excitation and more negative values corresponded to less excitation. Conversely, positive values of inhibitory connections indicated greater inhibition and less positive values indicated less inhibition. At a subject level, the analysis was the same as we have previously reported [15]. Specifically, we estimated the resting-state effective connectivity within the dACC-AI network by fitting a fully connected model [25,29]. To this aim, realignment, normalization At a subject level, the analysis was the same as we have previously reported [15]. Specifically, we estimated the resting-state effective connectivity within the dACC-AI network by fitting a fully connected model [25,29]. To this aim, realignment, normalization (to MNI space), and spatial smoothing (4-mm full width at half maximum with a Gaussian Kernel) were performed on the functional images. A general linear model (including six head movement parameters and time series corresponding to the white matter and cerebrospinal fluid as regressors) was fit to the images. A cosine basis set with frequencies ranging from 0.0078 to 0.1 Hz was also included in the general linear model [30]. Images were high-pass filtered to remove slow frequency drifts (<0.0078 Hz). By using an F contrast, we identified regions with blood oxygen level fluctuations within frequencies ranging from 0.0078 to 0.1 Hz [30]. Time series that summarized the activity within spheres (8-mm radius) in the right AI (MNI coordinates X = 38, Y = 20, Z = −4) and in the right DACC (MNI coordinates x = 1, y = 16, z = 38) were extracted and used to specify the dynamic causal models.
At a group level, we relied on parametric empirical Bayes (PEB) [31][32][33] to estimate the effect of GLU and GSH on connectivity parameters. We estimated a "two-metabolite" model with which we aimed to evaluate the evidence in support of the hypothesis that both GLU and GSH best explained the effective connectivity within the two-node network. The design matrix of the two-metabolite model comprised one column coding for group membership, one column comprising the mean-centered GLU levels, one column comprising the GLU × group interaction, one column comprising the mean-centered GSH levels, and one column comprising the GSH × group interaction (in this order). The design matrix also comprised a constant (column of ones). We compared the evidence in support of this model against the evidence in support of an "only-group" model, a reduced model comprising only the effect of the group on connectivity parameters.
We adjudicated between the two-metabolite and the only-group models by means of Bayesian model selection [34]. Specifically, we evaluated the evidence of each model (as estimated by the negative variational free energy, F). In principle, the strongest evidence is ascribed to the model with the least negative free energy. However, it is useful to assess the evidence of a given model relative to the evidence ascribed to other models. This is achieved by means of Bayesian model comparison in which the evidence of a given model (F1) is compared to the evidence of the model with the most negative free energy (F2), yielding the log of the Bayes factor (lnBF 1 = F 1 −F). In terms of posterior probability (PP), a BF > 20 is equivalent to a PP > 0.95 [34] which indexes very strong evidence. Therefore, we relied on a PP > 0.95 as a decision rule (i.e., threshold) for model selection [35]. Finally, the sum of the posteriors of all models' posteriors equaled to 1.

Between-Group Comparison in Metabolite Levels
The Bayesian linear model revealed higher GSH levels in the FES group than in the HC group (mode of the between-groups difference = 0.25, PP = 0.98; mode of the effect size = 0.71, PP = 0.98). The Bayesian analysis did not reveal an effect of the group on GLU levels (mode of the between-groups difference = 0.17, PP = 0.84; mode of the effect size = 0.1, PP = 0.84). Summary statistics of the posterior distributions of the model's parameters are reported in Table 2, and Figure 3 shows the posterior distributions of the estimated between-group difference in GSH and GLU levels.

Spectral Dynamic Causal Models of Effective Connectivity
The two-metabolite model, comprising the effect of GSH and GLU, performed better than the group-only model (PP > 0.99). As shown in Figure 4 and in line with our previous work [15], the activity of the inhibitory neurons in the dACC decreased as a function of GLU in the FES group (PP > 0.95). Crucially, in the current model, this effect was reversed by GSH, which was associated with increased inhibitory activity. This effect of GSH on IE connections (see also Figure 2 for reference) was observed not only in the dACC (PP >

Spectral Dynamic Causal Models of Effective Connectivity
The two-metabolite model, comprising the effect of GSH and GLU, performed better than the group-only model (PP > 0.99). As shown in Figure 4 and in line with our previous work [15], the activity of the inhibitory neurons in the dACC decreased as a function of GLU in the FES group (PP > 0.95). Crucially, in the current model, this effect was reversed by GSH, which was associated with increased inhibitory activity. This effect of GSH on IE connections (see also Figure 2 for reference) was observed not only in the dACC (PP > 0.95) but also in the inhibitory neural population of the AI (PP > 0.95).
Antioxidants 2021, 10, x 10 of 14 Figure 4. Counteracting effect of the GSH on the hyperglutamatergic state in the dACC-AI network. As previously reported [15], in the dACC the effect of GLU on IE connections (see also Figure 2 for reference) was weaker in FES patients than in HC subjects (indexed by the negative parameter estimate). This indicates stronger disinhibition of excitatory neuronal population with higher levels of GLU in patients than in controls, leading to a hyperglutamatergic state. Crucially, positive parameter estimates in blue indicate the effect of GSH on inhibitory (i.e., GABAergic) activity within the dACC and AI, this effect being stronger in FES patients than in HC subjects. Since the net activity of a given neuronal population is a linear function of the relevant parameter estimates, the magnitudes of these estimates indicate that the glutamatergic influence on intrinsic connectivity in the dACC is compensated by the "antioxidative" state.

Discussion
In drug-naïve patients with first episode psychosis, GSH levels were higher than in HC subjects. Higher GSH levels were related to stronger intrinsic inhibition within dACC and AI nodes of the salience network, a large-scale network known to play a cardinal role in schizophrenia symptoms [36]. This effect was in direct contrast to the relationship between higher GLU levels and putative disinhibition (i.e., reduced intrinsic inhibition) within the dACC. Our model provided an explanation for how higher levels of GSH can reverse the downstream pathophysiological effects of a putative hyperglutamatergic state in FES.
The presence of higher levels of GSH in patients compared to controls was in contrast to our meta-analytic observation of a small reduction of GSH in schizophrenia [37]. Nevertheless, as we reported in the same meta-analysis, patients with bipolar disorder had a small increase in GSH levels compared to healthy controls, leading to the speculation that GSH levels may mark the outcomes of psychotic disorders rather than the diagnosis per se. In fact, our prior observation from an overlapping sample with longitudinal clinical Figure 4. Counteracting effect of the GSH on the hyperglutamatergic state in the dACC-AI network. As previously reported [15], in the dACC the effect of GLU on IE connections (see also Figure 2 for reference) was weaker in FES patients than in HC subjects (indexed by the negative parameter estimate). This indicates stronger disinhibition of excitatory neuronal population with higher levels of GLU in patients than in controls, leading to a hyperglutamatergic state. Crucially, positive parameter estimates in blue indicate the effect of GSH on inhibitory (i.e., GABAergic) activity within the dACC and AI, this effect being stronger in FES patients than in HC subjects. Since the net activity of a given neuronal population is a linear function of the relevant parameter estimates, the magnitudes of these estimates indicate that the glutamatergic influence on intrinsic connectivity in the dACC is compensated by the "antioxidative" state.

Discussion
In drug-naïve patients with first episode psychosis, GSH levels were higher than in HC subjects. Higher GSH levels were related to stronger intrinsic inhibition within dACC and AI nodes of the salience network, a large-scale network known to play a cardinal role in schizophrenia symptoms [36]. This effect was in direct contrast to the relationship between higher GLU levels and putative disinhibition (i.e., reduced intrinsic inhibition) within the dACC. Our model provided an explanation for how higher levels of GSH can reverse the downstream pathophysiological effects of a putative hyperglutamatergic state in FES.
The presence of higher levels of GSH in patients compared to controls was in contrast to our meta-analytic observation of a small reduction of GSH in schizophrenia [37]. Nevertheless, as we reported in the same meta-analysis, patients with bipolar disorder had a small increase in GSH levels compared to healthy controls, leading to the speculation that GSH levels may mark the outcomes of psychotic disorders rather than the diagnosis per se. In fact, our prior observation from an overlapping sample with longitudinal clinical data supports this idea [23]. Higher levels of GSH are likely to indicate a more favorable prognosis, with a quicker response to antipsychotics in first-episode psychosis [23]. As such, the current sample of first-episode patients likely comprised subjects with more favorable outcomes than the chronic schizophrenia samples studied in our prior meta-analysis.
Our results also indicated that in early phases of psychosis, GSH may operate to reverse glutamate-mediated dysconnectivity. Specifically, stronger disinhibition of GABA neurons with higher levels of GLU reflected a hyperglutamatergic state in FES subjects-as indexed by the negative value of the parameter estimate representing inhibitory connections within the dACC. However, positive values of the effect of GSH indicated a direct relationship between the GSH level and a much stronger inhibitory (i.e., GABAergic) activity within both the dACC and AI. Since, as per the model's assumptions, the net activity of a given neuronal population is a linear function of the relevant parameter estimates, the magnitudes of these estimates suggest that the hyperglutamatergic state is compensated (or restrained) by an "antioxidative" state. This is important especially because our FES subjects were untreated when these data were collected. On this basis, we speculate that a targeted increase in dACC GSH levels via antioxidant supplementation or targeting the Nrf2 pathway could improve patients' response to antipsychotics [38]. More speculatively, this may assist in achieving an adequate response at lower-than-usual doses and cut down the total duration of higher dose exposure, both of which are now argued by some as key strategies to improve functional recovery in psychosis [39]. In the context of antioxidant trials, this also speaks to a stratification strategy based on baseline levels of GSH-as suggested in previous works [40].
The robustness of our results rests on several methodological strengths. First, we used 7T-MRS sequence with improved specificity to detect GSH resonance with reduce macromolecular interference [41]. This level of specificity contrasts with the specificity achieved using 3T-MRS [42,43]. Furthermore, the effective connectivity model is biologically grounded despite the fact that it does not consider the variability of inhibitory neuronal populations that have been recently reported [44]. Regardless of this limitation, the two-neuronal-population model was enough to evaluate our hypothesis. Finally, from our cross-sectional data, we cannot infer if the GSH increase is secondary to increased intrinsic inhibition or vice versa. Longitudinal fMRI and MRS data on GLU and GSH could potentially address this limitation in the future.

Conclusions
In summary, our data and computational model provide initial clues to understand the mechanistic effect of GSH on the previously reported hyperglutamatergic state within the dACC-AI network. As summarized in Figure 5, redox imbalance in early life may prime the brain for excitatory overdrive in schizophrenia; but if an appropriate increase in GSH accompanies glutamatergic excess, the inhibitory tone may be strengthened in compensation. Figure 5. Summary of the counteracting effect of GSH on redox imbalance in the salience network. In a healthy state, the excitation-inhibition balance is achieved via the interplay of excitatory and feedback inhibitory neuronal populations. Glutamate-mediated dysconnectivity in psychosis (likely caused by GSH depletion in the developing brain) would manifest in terms of interneuron deficit, leading to excitatory overdrive (i.e., hyperglutamatergic state). The resulting excitatory overdrive can be reduced by increasing inhibition via a compensatory increase in the levels of GSH, at least in a subset of patients. This illustration was produced using biorender.com

Data Availability Statement:
The data presented in this study are available on request from Dr. Lena Palaniyappan (lpalaniy@uwo.ca) Acknowledgments: We thank Trevor Szekeres, Scott Charlton, and Joseph Gati for their assistance in data acquisition and archiving. We thank Rob Bartha and Dickson Wong for consultation provided on MRS analysis. We thank all research team members of the NIMI lab and all the staff members of the PEPP London team for their assistance in patient recruitment and supporting clinical care. We gratefully acknowledge the participants and their family members for their contributions. Requests for data should be addressed to Lena Palaniyappan lpalaniy@uwo.ca. Figure 5. Summary of the counteracting effect of GSH on redox imbalance in the salience network. In a healthy state, the excitation-inhibition balance is achieved via the interplay of excitatory and feedback inhibitory neuronal populations. Glutamate-mediated dysconnectivity in psychosis (likely caused by GSH depletion in the developing brain) would manifest in terms of interneuron deficit, leading to excitatory overdrive (i.e., hyperglutamatergic state). The resulting excitatory overdrive can be reduced by increasing inhibition via a compensatory increase in the levels of GSH, at least in a subset of patients. This illustration was produced using biorender.com Author Contributions: Conceptualization, R.L. and L.P.; methodology, R.L., P.J. and J.T.; formal analysis, R.L.; writing-original draft preparation, R.L. and L.P. All authors have read and agreed to the published version of the manuscript.