The relation between parietal GABA concentration and numerical skills

Several scientific, engineering, and medical advancements are based on breakthroughs made by people who excel in mathematics. Our current understanding of the underlying brain networks stems primarily from anatomical and functional investigations, but our knowledge of how neurotransmitters subserve numerical skills, the building block of mathematics, is scarce. Using 1H magnetic resonance spectroscopy (N = 54, 3T, semi-LASER sequence, TE = 32 ms, TR = 3.5 s), the study examined the relation between numerical skills and the brain’s major inhibitory (GABA) and excitatory (glutamate) neurotransmitters. A negative association was found between the performance in a number sequences task and the resting concentration of GABA within the left intraparietal sulcus (IPS), a key region supporting numeracy. The relation between GABA in the IPS and number sequences was specific to (1) parietal but not frontal regions and to (2) GABA but not glutamate. It was additionally found that the resting functional connectivity of the left IPS and the left superior frontal gyrus was positively associated with number sequences performance. However, resting GABA concentration within the IPS explained number sequences performance above and beyond the resting frontoparietal connectivity measure. Our findings further motivate the study of inhibition mechanisms in the human brain and significantly contribute to our current understanding of numerical cognition's biological bases.

www.nature.com/scientificreports/ frontal regions were demonstrated to subserve magnitude representations, thus suggesting that both frontal and parietal regions are crucial for number processing. A more recent investigation directly examined the IPS's role in representing and manipulating numerical information and showed that anterior IPS is involved in encoding and that posterior parietal regions such as parietal lobules construct a working memory representation or orient spatial attention during number comparison 17 . Despite these advancements, our knowledge of how the concentration of the brain's major inhibitory (GABA) and excitatory (glutamate) neurotransmitters relate to numerical cognition has not yet been directly investigated. Unlike other neurochemicals, GABA and glutamate are abundant in the brain and directly affect brain activity, including regional blood-oxygen-level-dependent (BOLD) signal and brain connectivity in several networks including the default mode and the motor [18][19][20][21][22] . Importantly, recent studies demonstrated the role of glutamate and GABA in cortical excitation and inhibition, including in a single case of a mathematical prodigy, and in perceptual training, visual perception, attention, and cognitive skill acquisition [23][24][25][26][27][28][29][30][31][32] . These prior findings provide the motivation to examine the role of GABA and glutamate concentrations in tracking individual variation in numerical cognition. Therefore, the primary aim of the present study was to examine the neurotransmitter bases of several numerical tests involving rule-identification, reasoning, fluency, and knowledge of numerical facts. This would help address the question of which of those aspects of numerical cognition are associated with neurotransmitters. To this end, single-voxel proton magnetic resonance spectroscopy ( 1 H-MRS) was employed to examine the capacity of GABA and glutamate concentration within key numerical parietal and frontal regions in explaining the variance in the performance of various numerical tasks (numerical operations, mathematical reasoning, tempo test, number sequences, computational estimation, numerical agility) in a group of young adults. Given that both left and right frontoparietal regions were shown to underpin numerical processing 10,33,34 , and since we did not aim to interrogate the effect of brain laterality in the present study, we merely focused on left frontoparietal regions to keep the duration of the study to an acceptable length. Given the well-established role of frontoparietal brain connectivity in numerical cognition and the well-documented impact of neurochemical concentration on resting brain connectivity, we additionally employed resting functional magnetic resonance imaging (fMRI) to examine the complementary role of frontoparietal connectivity 10,19,20,22 . Overall, our aims were to (1) examine the relation between numerical cognition and resting neurochemical concentration within frontal and parietal regions and to (2) investigate how these relations depend on resting frontoparietal brain connectivity.

Results
Descriptive statistics and correlations between the cognitive tests. As a first step, we examined the behavioural results. The descriptive statistics of our cognitive tests can be seen in Table 1. To assess whether performance in the aforementioned cognitive tests is significantly related to each other we performed bivariate correlations, and as can be seen in Table 2 there was a generally positive relation between them.  Table 3 for all effect sizes (correlation coefficients, r).

IPS GABA and number sequences association is cognitively specific.
To delve more into the underlying structure of our six mathematical measures, we ran principal components analyses where a single factor was extracted (see Methods section) which was not significantly associated to IPS GABA (β = − 0.    . Since frontoparietal connectivity is associated with number sequences (step 1), and IPS GABA (step 2), and since IPS GABA still predicts number sequences after controlling for frontoparietal connectivity (step 3), then the IPS GABA can be viewed as a mediator on the relation between frontoparietal connectivity and number sequences 36 , but a formal mediation test showed that this was not the case even though there was a trend (Fig. 3, CI of the indirect effect = [− 0.003 0.14]). Taken together our results suggest that resting IPS GABA tracks individual variation in number sequences above and beyond resting left frontoparietal connectivity.

Discussion
The present study investigated the relation between the brain's major inhibitory and excitatory neurotransmitters, GABA and glutamate, and numerical cognition as assessed using several key numerical measures. We found a negative association between a complex numerical measure, the number sequences, and the amount of GABA within the left IPS. This was followed by (1) a cognitive, (2) a neurotransmitter, and (3) a regional specificity supporting the explicitness of our finding. Subsequently, we demonstrated a positive association between left frontoparietal connectivity and number sequences performance. Lastly, we show that the resting GABA concentration within the IPS is a stronger and independent predictor of number sequences above and beyond the observed resting frontoparietal connectivity.
The main finding of the present study was the negative association between GABA concentration in the left IPS and number sequences performance. The number sequences task required the identification of a logical rule that regulates the relations between numbers in a sequence. The process of identifying such a rule involves a reiterative formulation and testing of hypotheses along with continuously performing calculations. As an example, testing the "increasing difference" hypothesis in the sequence 1 3 6 10 15 21 _, one would need to perform several subtractions (e.g., 3-1, 6-3, 10-6, etc.), while keeping active the partial results in the aforementioned frontoparietal networks, and verify whether the hypothesized rule is respected (i.e., differences increase by one unit) 10 .
Notably, the number sequences task contains underlying cognitive components also assessed in other administered tasks. First, numerical procedures are included in numerical operations, mathematical reasoning, and the tempo test. Second, numerical flexibility is included in other tests like the numerical agility test. Third, the seek for a logical rule is also measured outside the numerical domain using the matrix reasoning task. However, the logical rule in the number sequences is numeric unlike the one in the matrix reasoning task. Taken together, we conclude that the GABA concentration in the left IPS tracks a unique and specific set of abilities concerning the identification of a logical rule particularly in the numerical domain that additionally involves numerical flexibility.
An open question for future research is the delineation of how IPS GABA affects the numerical brain at the network level. The biological significance of GABA modulation during rest was established in several settings. First, it was shown that variation in baseline GABA level, as measured by MRS, correlates with blood oxygenation-dependent (BOLD) fMRI signal in the visual, somatosensory and anterior cingulate cortex regions [37][38][39][40] . Moreover, GABA concentration during rest was additionally associated to several sensory and cognitive functions, including monocular deprivation, perceptual training, visual perception, attention and cognitive skill www.nature.com/scientificreports/ acquisition 24,25,27,30,32 . Here, we extend the biological significance of regional resting GABA concentration in tracking individual variation in number sequences performance. Given the well-established role of frontoparietal brain connectivity in numerical cognition and the well-documented impact of neurochemical concentration on resting brain connectivity, we subsequently utilized resting fMRI to examine the role of frontoparietal connectivity. Our results suggest that frontoparietal connectivity does not uniquely contribute to number sequences once GABA concentration is taken into account, while GABA concentration can predict number sequences above and beyond frontoparietal connectivity. Therefore, it is currently unknown how IPS GABA affects the numerical brain network level. Future studies can delve more into the role of IPS GABA at the brain network level by linking IPS GABA concentration to frontal and parietal regional activations obtained in a task-based fMRI setting. This multi-modal approach has the potential to offer further insights concerning the neurofunctional relevance of IPS GABA in the context of numerical cognition at the brain network level. In particular, if the relation between IPS GABA and number sequences performance is mediated by the individual variation in task-based brain activity, this set of findings would shed some light on the neurodynamic effect of IPS GABA at the brain network level. By interrogating regional specificity, we found that parietal but not frontal GABA, and parietal GABA but not parietal glutamate, was related to number sequences. Why is parietal but not frontal GABA (and not glutamate) an important index of numerical cognition? Two competing hypotheses have been proposed concerning neurochemical mechanisms across the brain 41 . The first postulates that neurochemical concentration is uniquely determined in each region and each individual by genetic and environmental influences, thus neurochemicals in different regions may exert different influences on sensory and cognitive processing 30,[42][43][44] . The second theory postulates that neurochemical concentration across different regions may be very similar due to common embryonic origins or shared subcortical projections [45][46][47][48][49] . Such neurochemical selectivity was previously demonstrated in other domains of neuroscience including visual perception and attention 27,41 . Our results support both a regional and a neurochemical selective influence on number sequences, therefore, supporting the first theory. Overall, unlike the IPS, the neurochemical contribution of MFG did not track numerical performance. Meta-analyses suggested a hierarchical contribution of prefrontal cortex in numerical cognition 9 . Namely, the inferior frontal gyrus was typically engaged in relatively simple numerical tasks, while the MFG is involved in more complex tasks which require several procedural steps or increase storage load 9 . This role of MFG likely reflects shared links to working memory which is behaviourally related to numerical performance 50,51 and supported by the prefrontal cortex 52,53 . Indeed, a recent study employing functional MRS found elevated glutamate levels in dorsolateral prefrontal cortex during the execution of a 2-back task compared to passive visual fixation 54 . Given this contribution of MFG in demanding computations that mirror the computations underlying the present cognitive tasks, one potential reason we did not find an association between MFG neurotransmitter levels and performance may be accounted for the fact that neurotransmitter levels were measured at baseline rather than during the execution of the numerical tasks which is indeed a limitation of the present study. This consequently prevents us from drawing conclusions regarding the task-based contribution of neurotransmitter and frontoparietal connectivity measures which are likely to differ compared to resting contributions.
Another limitation of the present study is common to other MRS studies and relates to the difficulty of distinguishing between intracellular and extracellular neurotransmitter concentrations or even a portion of these based on the MRS signal alone 55 . Consequently making direct inferences of cortical excitability/inhibition based on the neurotransmitter concentrations alone should be done with caution. Similar limitations concerning the exact nature of the neural signal are not specific to MRS but exist with other modalities e.g., diffusion and structural MRI 56,57 . Therefore, our finding demonstrates the pressing need for a multi-modal imaging approach for examining GABAergic mechanisms involved in numeric cognition.
In sum, by investigating numerical cognition with MRS, we identified a specific neurotransmitter, GABA in the IPS, a key brain region involved in numeracy, that is associated to a number sequences task. Our findings shed light on the possible mechanisms that are involved in individual differences in numeracy. Given the involvement of GABA in neuroplasticity 58 , our findings raise the question of whether such mechanisms are involved in the www.nature.com/scientificreports/ learning and development of numerical skills. These findings also motivate future examination of the underlying causality mechanisms which could be assessed using non-invasive brain stimulation approaches 25,59-61 .

Methods
Participants. We recruited 54 university first-year undergraduate participants (21 females, Mean age: 18.89, SD: 0.62). The completion of the imaging session lasted (~ 60 min), and the completion of the cognitive-behavioural testing lasted (~ 60 min); these sessions were part of a more extensive battery that included several other cognitive and behavioural assessments. All imaging data were acquired on a single scanning session. Because the magnetic resonance spectroscopy sequence we utilized measures baseline-level neurochemical concentration, participants watched a movie during the acquisition rather than performing the numerical tasks. During the "eyes-open" resting fMRI scan, participants were asked to fixate on a white cross on a black background. All participants were predominantly right-handed, as measured by the Edinburgh Handedness Inventory 62 . We confirmed via a questionnaire that they were free of current or past neurological, psychiatric or learning disability or condition that might affect cognitive or brain functioning. As compensation, participants received £50. Informed written consent was obtained. All methods were carried out in accordance with relevant guidelines and regulations. MRS neurotransmitters were quantified with the LCmodel 65 , using a basis set of simulated spectra generated based on previously reported chemical shifts and coupling constants based on a VeSPA (versatile simulation, pulses, and analysis) simulation library 66 . Simulations were performed using the same RF pulses and sequence timings as in the 3T system described above. Absolute neurotransmitter concentrations were extracted from the spectra using a water signal as an internal concentration reference. The exclusion criteria for data were Cramer-Rao lower bounds (CRLB, the estimated error of the neurotransmitter quantification) > 50% were classified as not detectable. Additionally, we excluded cases with an SNR beyond 3 standard deviations and a concentration value (or a cognitive or a connectivity value, see sections below) that was beyond 3 standard deviations. Raw metabolite-to-water amplitude ratio neurotransmitter concentrations were then corrected using the structural properties of the selected regions. Namely, we segmented the images into different tissue classes, including grey matter (GM), white matter (WM), and cerebrospinal fluid (CSF) using the SPM12 segmentation facility. Next, we calculated the number of GM, WM, and CSF voxels within the two masks of interest separately around the left middle frontal gyrus (MFG) and the left intraparietal sulcus (IPS) in native space. Subsequently, we divided these six numbers (GM, WM, and CSF for IPS and MFG) by the total number of GM, WM, and CSF voxels creating the corresponding GM, WM, and CSF fraction values per participant and region. As a final computation step, we corrected the raw metabolite-to-water amplitude ratio neurotransmitter values to these structural fractions using the following LCmodel 65 computation, as can be seen in Eq. 1. The numbers in the equation (43,300, 55,556, and 35,880) represent the amount of water concentration in the white matter (35,880), grey matter (43,300) and CSF (55,556), and these values were selected based on the LCmodel manual recommendations.
To minimize the potential confounding effects T2-relaxation times we corrected the concentration from Eq. 1 based on the T2 of tissue water values as can be seen in Eq 2. Fully relaxed unsuppressed water signals were acquired at TEs ranging from 32 to 4040 ms (TR = 15 s) to water T2 values in each region (32 ms, 42 ms, 52 ms, 85 ms, 100 ms, 115 ms, 150 ms, 250 ms, 450 ms, 850 ms, 1650 ms, 3250 ms, 4040 ms). The transverse relaxation times (T2) of tissue water and percent CSF contribution to the region were obtained by fitting the integrals of the unsuppressed water spectra acquired in each region at different TE values with a biexponential fit 67 , with the T2 of CSF fixed at 740 ms and three free parameters: T2 of tissue water, the amplitude of tissue water, and amplitude of CSF water. These analyses were conducted on a subject-by-subject basis. Our analysis pipeline was also used previously in an overlapping cohort of participants 35 . We used the LCModel default output with the following parameters: ATTMET, ATTH2O, and WCONC, which in our control files were set as follows ATTMET = 32, For triangulation purposes, we additionally calculated the metabolite ratios with respect to total creatine, which was calculated as the raw metabolite-to-water amplitude ratio neurotransmitter concentration referenced the concentration of total creatine (i.e., tCr, where tCr = creatine and phosphocreatine concentration).
Resting fMRI. Functional images were acquired with a multi-band acquisition sequence (Multi-band accel. factor = 6, TR = 933 ms, TE = 33.40 ms, flip angle 64°, number of slices = 72, voxel dimension = 2 × 2 × 2 mm, number of volumes: 380). rsfMRI data were pre-processed and analyzed using the CONN toolbox (www. nitrc. org/ proje cts/ conn, RRID: SCR_009550) in SPM12 (Wellcome Department of Imaging Neuroscience, Institute of Neurology, London, UK) using the default pre-processing pipeline "MNI-space direct normalization" 68 . Functional volumes were motion-corrected, slice-timed corrected, segmented, normalized to a standardized (MNI) template, spatially smoothed with a Gaussian kernel (8-mm FWHM) and a pass filtered (0.01 Hz to Inf). The low-pass portion of the filter (0.01 Hz) was used to reduce the physiological and noise components contributing to the low-frequency segment of the BOLD signal, the high-pass portion of the filter was set to Inf to additionally allow higher frequencies to be included (above 0.1 Hz) given our fast acquisition. It also increased the degrees of freedom improving the quality of denoising and connectivity estimation. We also accounted for physiological noise in the time-series by regressing out the confounding effects of white matter, CSF, realignment, and scrubbing that were automatically calculated in CONN. All frontal regions of interest (ROI) were defined based on the Harvard-Oxford Atlas (superior frontal gyrus, (SFG), middle frontal gyrus (MFG), and inferior frontal gyrus (IFG)), apart from the IPS which was defined based on the Dorsal Attention network atlas, both atlases are featured within the CONN toolbox. The reason we used a different atlas for the IPS is that this region is not featured in the other atlas. We calculated the (frontoparietal) ROI-to-ROI connectivity which we then added into the simple regression and multiple regression models. Instead of merely using the MFG, the choice of also (2) T2 -corrected concentration = tissue -corrected concentration * exp(−TE/T2)   70 , the number sequences subtest of the numerical aptitude test 71 , the computational estimation 72 and the numerical agility 73 tests. The numerical operations subtest is composed of written numerical problems, which require the implementation of numerical procedures with minimal/absent time pressure.
The mathematical reasoning subtest is composed of mathematical problems, which require participants to create a mental model of the problem, extract relevant information, and then select and execute the appropriate operation 74 . We calculated the proportion of correct responses for the numerical operations and the mathematical reasoning subtests. We did not use the standardized scores as all participants were between 18 and 21 years of age.
The tempo test entails five columns (addition, subtraction, multiplication, division, and mixed), each composed of 40 numerical problems (e.g., 7 + 8 = ….). Each column is presented sequentially with the instruction to solve as many problems as possible within 60 s. The time constrains make the tempo test a widely used measure of numerical fluency. For the tempo test instead, we calculated the proportion of correct responses in the five columns, and then we divided it by the individual solving time divided by total time at disposal (i.e., 300 s). This calculation was particularly done to take into account the fact that several participants completed the columns within one minute.
The computational estimation task consists of ten numerical problems: five multiplications (e.g., 64.6 × 0.16) and five divisions (e.g., 0.76/0.89) which were presented simultaneously for one minute, and participants were asked to write down their best estimate for each numerical problem. Participants were explicitly told to estimate without calculating the exact answer. For each response, we calculated the proportion of its absolute deviation from the exact answer (the proportion of absolute deviation =|estimate-correct answer|/correct answer). When estimates differed from the correct answer by more than the value of the correct answer, such responses were deemed outliers and were replaced with '1' , as missing values. Therefore, '1' was the maximum error score. Outliers were replaced rather than removed to avoid unfairly disadvantaging participants who misunderstood the task (e.g., by failing to notice the change in numerical sign from multiplication to division). Missing values were also replaced with '1' rather than removed to avoid unfairly advantaging participants who decided not to complete the task and to spend more time on fewer calculations. Therefore, we decided to use this imputation method to avoid giving an advantage to slow responders. For instance, a given slow respondent may have provided three accurate estimates (i.e., a small deviation in proportion), while a given fast responder may have provided 10 less accurate estimates (i.e., large deviation in proportion). The calculation of average deviation would indicate a higher score for the slower than the fast responder. Therefore, by conversely giving a "1" to the missed responses of the slow participant somewhat balanced this discrepancy. To enhance the interpretability of the results we reversed the final computational estimation score by subtracting it from one, thus positive scores indicated higher performance.
In the numerical agility task, participants were asked to generate the number 24 from four presented numbers (for example with the numbers 7, 5, 5, 4 a solution could be (7-5/5) * 4). There were five problems in total, which were administered in order of difficulty. There was a two-minute time limit for each problem. The instructions made it clear that only the four basic operations were allowed, such that any individual having completed secondary school should be able to solve these problems. For each item, two points were awarded if the correct answer was achieved within one minute, one point was awarded if it was achieved within two minutes, and no points were awarded otherwise.
The number sequences subtest is composed of 26 complex number sequences reflecting a specific principle/ rule. Participants had to complete the sequences by adding the last number of the sequence from a pool of five numbers at disposal (e.g., 1 3 6 10 15 21 _; Choices: 26 27 28 29 30). Participants had 15 min at disposal to solve as many sequences as they could. For this test, we calculated the proportion of correct responses. To assess the cognitive specificity of our associations, we also assessed the participants' general cognitive ability using the matrix reasoning test 75 which is a 30-item (of increasing difficulty) assessment tool which requires identifying a logical pattern in a sequence of visuospatial stimuli. The subtest is interrupted when the participant provides 3 consecutive wrong responses. We calculated the proportion of correct responses.

Statistical analyses.
To assess the association between neurotransmitter concentration and numerical performance, we employed bivariate correlation analyses, as well as multiple regressions. As can be seen in Table 4, some cognitive variables were not normally distributed. In these cases and when applicable we performed Spearman's correlation (correlation analyses) and bootstrapping methods (multiple regression analyses), both of which are not sensitive to the assumption of normality. The bootstrapped P-values in the results section were obtained from 10,000 samples at 95% confidence intervals (CI) (SPSS v25) 76 . Our Bonferroni corrected P-value was 0.0021 derived from dividing our alpha value (0.05) by 24 [(4 brain measures (GABA MFG, GABA IPS, glutamate MFG, glutamate IPS)*6 numerical tasks)].
To delve more into the underlying structure of the six mathematical measures we performed principal component analyses. The observed correlation matrix was significantly different from an identity matrix as assessed with Bartlett's test of sphericity (Approx. Chi-Square(15) = 62.08, P < 0.001) and the Kaiser-Meyer-Olkin measure of sampling adequacy was 0.79 which indicates that this set of variables are appropriate for component or factor analyses [77][78][79] . One factor was extracted which accounted for 45.46% of the variation in our measures. In particular, this factor's association with each of the measures was as follows (tempo test, r S (50) = 0. 70 www.nature.com/scientificreports/ numerical operations, r S (50) = 0.57, P < 0.001; mathematical reasoning, r S (50) = 0.77, P < 0.001; number sequences, r(50) = 0.69, P < 0.001; numerical agility test, r S (50) = 0.56, P < 0.001; computational estimation, r(50) = 0.63, P < 0.001).