How does the human brain represent the meaning of abstract symbols? Some theories postulate the existence of semantic spaces where concepts occupy positions that reflect their conceptual relationships. In the number domain, psychological evidence suggests that integers are represented along a mental number line which, with education, integrates higher-level number concepts such as fractions. To test this hypothesis, we recorded whole-brain 7T fMRI responses to integer and fraction symbols during an active magnitude comparison task designed to engage magnitude processing. Consistent with predictions, we found both behavioral and neural numerical distance effects. Activation vectors in intraparietal, inferior temporal, prefrontal, hippocampal, and parahippocampal cortices formed a two-dimensional semantic space organized by numerical magnitude and domain (fractions versus integers). Gaussian fits revealed a topographic map of numerical preferences in the anterior inferior parietal cortex, common to both domains. Our results suggest that, when educated adults engage in an explicit magnitude comparison task, their brain can construct an ordered representation that integrates fractions and integers, with a topographic organization in the anterior IPS/IPL.

Ethics statement

This 7T study was approved by the local ethics committee (Comité de Protection des Personnes Sud-Ouest et Outre-mer III; references: CEA 100 055, ID RCB 2020-A01713-36), and the study was conducted following the Declaration of Helsinki. Written informed consent was obtained from all participants, who were compensated for their participation.

Participants

Twenty-two healthy adult participants were recruited for this study (13 women and 9 men, age range: 19–39 years, mean age = 26.4, SD = 5.6). All participants were right-handed, native French speakers, and held a high-school degree with a scientific background. They had normal or corrected-to-normal vision and reported no history of psychiatric or neurological disorders.

One participant was excluded due to a cyst in the temporal lobe that could affect that region’s anatomy and functional connectivity. Another participant was excluded after experiencing a panic attack in the scanner and requesting to leave it. Additionally, for a different participant, one of the response buttons failed, preventing analysis of his behavioral responses. As a result, we included 20 participants in the fMRI analyses and 19 in the behavioral analyses.

Number comparison task

Participants were engaged in a serial continuous comparison task, inspired by Nelli and colleagues [30], in which each target had to be judged as larger or smaller than the previous one. To achieve that, we avoided consecutive numbers with the same magnitude. The targets consisted of 12 numbers: six integers (i.e., 0, 1, 2, 3, 4, and 5), and six fractions (, , , , , and ). The stimuli were presented as white Arabic numerals on a black screen for 1 s, with fractions shown using a horizontal line. The small set of familiar fractions was selected to avoid excessive task difficulty and discouragement among participants. Importantly, the fractions were constructed from the same integers, thus allowing us to test whether participants represented magnitude rather than perceptual features (e.g., whether is closer to smaller numbers than to 4). All possible transitions between the 12 target numbers occurred an equal number of times, with no number repeated consecutively, resulting in 11 presentations of each target and 132 unique transitions per run. This relatively small stimulus set facilitates interpretable results in MDS and RSA, while ensuring that each number was preceded exactly once by every other number.

Before starting the fMRI experiment, participants completed a few trials to familiarize themselves with the task and stimuli. Each participant was assigned a number based on their arrival order. Participants with an even number were instructed to press the “left button” when the current number was greater than the previous one, and the “right button” when the current number was lower. In the second run, the buttons were reversed. Participants with an odd number started with the opposite button assignments. Each participant completed two runs, each lasting 9 min and 42 s. After the number presentation, a pseudo-randomly jittered interval of 2–4 s followed, fixed by the experimenter to ensure a mean SOA of 4 s, during which the fixation cross turned green if the response was recorded. Subsequently, the fixation cross turned red 500 ms before the onset of the next stimulus to prepare participants (see Fig 1A). Participants received a 10 s rest time for every 45 stimuli, the fixation cross was on the screen the entire time, and the feedback with the percentage of correct answers appeared for 3 s. Before starting the next block, the fixation cross turned red to prepare participants. At the beginning of each block, the last number before the rest was repeated when participants returned to the task.

Math localizer

To localize math-responsive areas within each individual, we used a highly reproducible sentence-based paradigm previously employed multiple times in our lab [51–55], which has the advantage of not including symbolic numbers in Arabic format. Participants read 60 sentences and had to decide whether they were True or False. Half of the sentences were mathematical statements, either stating arithmetic calculations in verbal form or describing geometrical facts concerning basic shapes (e.g., “Number three multiplied by twelve is thirty-six” or “A rectangle has only one axis of symmetry”). The remaining sentences were nonmathematical, including general factual knowledge or sentences referring to a person’s perspective (e.g., “The camel is a large humpbacked mammal” or “For Einstein the atomic bomb is a harmless object”). In each category, half of the sentences were true (T), and the other half were false (F), with the response buttons (left and right) counterbalanced across participants. Participants completed a single run lasting ~11 min.

The visual stimuli were balanced to control the number of words, total number of characters, word length, and frequency of occurrence in French. They were presented using rapid serial visual presentation (RSVP) at 350 ms per word, with an average total duration of 3.2 s. Each sentence was followed by a fixation cross, jittered to maintain an average stimulus onset asynchrony (SOA) of 8 s. A red cross was presented 350 ms in advance of the upcoming stimulus. Participants were instructed to answer as quickly and as accurately as possible, and response times were measured relative to the onset of the last word.

Stimulus delivery and response collection of both Localizer and Numbers tasks were controlled by Expyriment [56].

fMRI acquisition

MRI data were acquired using a 7T MAGNETOM whole-body MR scanner (Siemens Healthineers, Erlangen, Germany) with a 32-channel head coil (Nova Medical, Wilmington, USA) at the NeuroSpin Center of the French Alternative Energies and Atomic Energy Commission. Structural MRI data were collected using T1-weighted rapid gradient echo (MPRAGE) sequence (repetition time (TR) = 5,000 ms, echo time (TE) = 2.51 ms, resolution = 0.65 mm isotropic, field of view (FoV) = 208 ⨯ 208 mm2, flip angle = 5/3, bandwidth (BW) = 250 Hz/px, echo spacing = 7 ms). Functional MRI (fMRI) data were acquired using a T2*- weighted gradient echo planar imaging (EPI) sequence (TR = 2,000 ms, TE = 21 ms, voxel size = 1.2 mm isotropic, multiband acceleration factor = 2, FoV = 192 ⨯ 192 mm2, flip angle = 75, BW = 1,488 Hz/px, partial Fourier = 6/8, echo spacing = 0.78 ms, number of slices = 70 (no gap)). To correct for EPI distortion, a five-volume functional run with the same parameters except for the opposite phase encoding direction (posterior to anterior) was acquired immediately before each task run. Manual interactive shimming of the B0 field was performed for all participants. The system voltage was set at 250 V for all sessions, and the fat suppression was decreased per run to ensure that the specific absorption rate did not surpass 62% for all functional runs.

Anatomical and functional preprocessing

Both anatomical and functional MRI data were preprocessed using fMRIPrep 20.0.2. This pipeline incorporates standard preprocessing steps within a robust and user-friendly workflow, ensuring reliable performance across datasets independently of the data idiosyncrasies and highly consistent results [57].

The T1-weighted (T1w) anatomical image was first corrected for intensity nonuniformity (INU) with N4BiasFieldCorrection [58], distributed with ANTs 2.3.3 [59], RRID: SCR_004757), and used as T1w-reference throughout the workflow. Next, the T1w-reference was then skull-stripped with a Nipype implementation of the ANTs antsBrainExtraction.sh workflow, with the OASIS30ANTs template as the target. Brain tissue segmentation into cerebrospinal fluid (CSF), white matter (WM), and gray matter (GM) was performed on the skull-stripped T1w image using fast from FSL [60]. Cortical surface reconstruction was carried out using recon-all from FreeSurfer 7.2.0 [61], and the resulting brain mask was subsequently refined with a custom variation of the Mindboggle method [62]. Finally, volume-based spatial normalization to the MNI152NLin2009cAsym standard space was performed through nonlinear registration with antsRegistration (ANTs 2.3.3), using brain-extracted versions of both the T1w reference and the T1w template.

A reference volume and its skull-stripped version were generated by aligning and averaging a single-band reference (SBRefs). Head-motion parameters relative to the BOLD reference image, including transformation matrices and the six corresponding rotation and translation parameters, were estimated before any spatiotemporal filtering using MCFLIRT (FSL 6.0.5; [63]). The estimated fieldmap was then aligned to the target echo-planar imaging (EPI) reference run. The resulting field coefficients were then mapped onto the reference EPI using the computed transformation. BOLD runs were corrected for slice-time differences using 3dTshift from AFNI [64], with all slices temporally realigned to 0.975 s, corresponding to the midpoint of the slice acquisition interval (0–1.95 s). The BOLD reference image was then co-registered to the T1w anatomical reference using bbregister (FreeSurfer), which implements boundary-based registration [65]. Co-registration was performed with nine degrees of freedom to account for residual distortions present in the BOLD reference image. A reference volume and its skull-stripped version were first generated using a custom methodology implemented in fMRIPrep. Several confounding time series were then calculated from the preprocessed BOLD data, including framewise displacement (FD), DVARS, and three region-wise global signals. FD was calculated using two formulations: the absolute sum of relative motions described by [66] and the relative root mean square displacement between affine transformations proposed by [63]. Both FD and DVARS were computed for each functional run using their implementations in Nipype, following the definitions by [66]. Global signal regressors were extracted from the CSF, the WM, and the whole-brain masks. In addition, physiological noise regressors were estimated using component-based noise correction method (CompCor, [67]). Principal components were derived from the preprocessed BOLD time series after high-pass filtering with a discrete cosine basis set using a 128 s cut-off period. Two CompCor variants were implemented: temporal (tCompCor) and anatomical (aCompCor). For tCompCor, principal components were calculated from the top 2% variable voxels within the brain mask. For aCompCor, probabilistic masks corresponding to CSF, WM, and combined CSF + WM compartments were generated in anatomical space. This mask was obtained by dilating a GM mask extracted from the FreeSurfer’s aseg segmentation, and it ensures components are not extracted from voxels containing a minimal fraction of GM. Finally, these masks were resampled into BOLD space and binarized using a threshold of 0.99. Components were also calculated separately within the WM and CSF masks. For each CompCor decomposition, the k components with the largest singular values were retained, such that the selected components’ time series explained at least 50% of the variance within the corresponding nuisance mask (CSF, WM, combined, or temporal mask). Components that did not contribute to this threshold were discarded. The head-motion estimates obtained during the correction step were included in the corresponding confounds file. To further model nuisance-related variance, the confound time series derived from head motion estimates and global signals were expanded by including temporal derivatives and quadratic terms [68]. Volumes with FD greater than 0.5 mm or standardized DVARS exceeding 1.5 were identified as motion outliers. Additional nuisance time series are calculated using principal components analysis of the signal found within a thin band (crown) of voxels around the edge of the brain, following the approach proposed by [69]. The BOLD time series were resampled into standard space, generating a preprocessed BOLD run in MNI152NLin2009cAsym space. First, a reference volume and its skull-stripped version were generated using a custom procedure implemented in fMRIPrep. All resampling can be performed with a single interpolation step by composing all the pertinent transformations, including head-motion transform matrices, susceptibility distortion correction (when available), and co-registration transforms between functional, anatomical, and output spaces. Volumetric resampling were performed using antsApplyTransforms (ANTs), employing Lanczos interpolation to reduce the smoothing effects commonly associated with other interpolation kernels [70]. Surface-based resampling was carried out using mri_vol2surf from FreeSurfer.

Data analysis

We analyzed behavioral data by integrating response times and error rates into a Balanced Integration Score (BIS; [71]), separately for each combination of the 12 numbers, and as a function of the numerical distance between the current and previous number. To ensure that all distance values were nonnegative, we subtracted, for each participant, the minimum BIS value from all cells in the matrix. To the BIS values, we fitted a linear mixed-effects model in which trial type (fraction versus integer domain for both the current and the previous stimuli) and numerical distance (the log of the absolute difference between current and previous numbers) were included as interacting fixed effects. To account for repeated measures and individual differences in baseline performance, participant was included as a random effect.

Statistical analyses of the fMRI data were performed using NiLearn (RRID: SCR_001362) and rsatoolbox [72]. In the general linear models (GLM), regressors were convolved with the canonical hemodynamic response function (HRF). GLM1 included 17 regressors of interest: one for each current target number (12 regressors); two for left and right button presses; one for RT difficulty, calculated as the current RT minus the participant’s average RT; one for feedback onset; and one for the first target presentation, where participants were not required to respond. Additionally, to regress out motion and physiological noise, we included six head motion parameters, the CSF and white matter signal, 10 anatomical CompCor (computed by fMRIPrep), and 14 cosine drifts. This GLM was used for the PCA, RSA, searchlight, MDS, and Gaussian fit analyses. The data remained unsmoothed for all analyses, but for visualization purposes, we smoothed the data at the group level with a 6 mm kernel. We chose to use unsmoothed data to preserve sensitivity to subtle patterns and fine-grained distinctions between numbers. An exception was made for the Gaussian fit analysis, where the data were smoothed at the individual level with a 4 mm kernel as well as the results of the group-averaged-data. Possibly due to a limited amount of data, smoothing was necessary, but the results remained when using either 2.4 mm or 4 mm smoothing.

GLM2 was very similar to GLM1, except that we excluded the RT difficulty regressor and included the logarithm of the absolute distance, calculated as log(|current number − previous number|). This model was used exclusively for Fig 2C, and the data were smoothed at the group level with a 6 mm kernel.

Principal component analysis (PCA)

We used PCA on the 12 target-related voxel-wise model weights to identify the most important orthogonal dimensions for the task. We averaged the t-maps for each target number across participants using unsmoothed data (similar results were obtained at the individual level in the majority of participants). Then, we created a matrix where each row corresponded to a target number, and each column corresponded to a voxel. We applied PCA to these weights. We determined that two components were optimal by examining the elbow of the Scree Plot. To project the results onto the cortical surface of a common space across participants, we selected the top 1% of voxels, and this data was subsequently smoothed with a 6 mm kernel. This threshold allowed us to isolate the voxels with the most robust and reliable response to each component, whereas more lenient thresholds preserved the results but extended the coverage to a larger portion of the brain.

Representational similarity analysis (RSA)

We identified brain areas showing greater activation for math-related sentences compared to nonmath-related sentences in the localizer task, at the whole-brain level for each participant (each thresholded at p < 0.001 uncorrected). From the contrast math > nonmath sentences, we also manually selected the parietal, frontal, and temporal ROIs at the individual level. BOLD representational dissimilarity matrices (RDMs) were constructed by calculating the correlation distance d = 1−r, where r is the Pearson correlation coefficient between the within-subject multivoxel patterns elicited by each number, resulting in a 12 x 12 RDM for each participant, which was then averaged across participants. To visualize neural geometry, we used a two-dimensional metric inter-individual difference scaling (INDSCAL) MDS model, which takes into account the individual participants’ RDMs, implemented using the R-package SMACOF [73,74].

To identify cortical sectors in which the neural RDM correlated with a specific RDM—particularly focusing on the distance across the domains of fractions and integers—while accounting for the influence of remaining RDMs, we performed a multiple-regression searchlight RSA. This analysis was performed within a mask of mathematically responsive areas (math > nonmath sentences, FDR q < 0.05) as defined by the localizer, to which we added Glasser-defined ROIs of hippocampi and parahippocampal cortex. The neural RDM for each searchlight was computed using correlation distances. For the searchlight analysis, we included all voxels within a radius of 4 mm of the central voxel. For each searchlight sphere, we computed neural RDMs from the condition-by-voxel matrix of estimated neural responses using Pearson correlation distance between pairs of conditions. This procedure was carried out at the individual participant level, followed by a second-level analysis across participants, using FDR correction q < 0.001 for multiple comparisons. The multiple regression model included five predictor RDMs: distance within fractions, category of fractions, distance within integers, category of integers, and distance across the domains of integers and fractions.

Voxel tuning

For this analysis, we averaged the t-maps for each target number across participants using smoothed data (4 mm FWHM). The same analysis was also performed at the individual level in each participant’s native space. We then created a matrix in which each column represented a target t-value and each row corresponded to a voxel within the mask of math-related regions defined by the localizer math > nonmath sentences, FDR q < 0.05, as well as the Glasser ROIs of IPS, pITG, MFG (see S4 Fig), parahippocampal, and hippocampus. Separate analyses were conducted for integers and fractions. Each voxel’s values were normalized between their minimum and maximum across conditions of integers or fractions, resulting in a range from 0 to 1. A Gaussian function was then fitted using the Levenberg–Marquardt method within a bounded parameter space, yielding for each voxel a preferred value on a continuous scale from 0 to 5 for integers and from 0.25 to 2.5 for fractions. The amplitude was constrained between 0 and 1, and the width between 0.1 and 2. At the single-subject level, maps were thresholded at R² > 0.3 because this threshold allowed us to reliably inspect maps across the 20 participants while ensuring consistency. Similar thresholds have been used in pRF analysis at the single-participant level, even after averaging data across multiple sessions (e.g.[7,75]). In contrast, group-averaged data provides a significantly higher signal-to-noise ratio, allowing us to use a stricter threshold of R² > 0.5. The same threshold was applied both to individual runs and to the average of the two runs.

Using data from one of the runs (1 or 2), we categorized voxels into 6 bins by rounding their preferred values for integers and, for fractions, by applying category boundaries based on the subtraction (next fraction - previous fraction)/2. Then, we plotted the normalized activation using data from the other run and averaged the two resulting plots (finding preferred numbers with run 1 and plotting the activation of run 2, and vice versa). To assess statistical significance, we organized the average data in a matrix (6 x 6 for integers, 6 x 6 for fractions) where columns indicated the preferred numbers and rows the average normalized activation to each target number. Then, we tested whether the diagonal elements of this matrix were significantly larger than the off-diagonal elements. For that, we shuffled the columns to break the correspondence between rows and columns and computed a permutation test (n = 720 permutations).

We correlated each voxel’s preferred value for fractions and for integers within math-responsive areas (pITG, MFG, and IPS) and in aIPS/IPL (Glasser ROIs: PF and IP2). To assess the significance of these correlations, we generated 10,000 random white-noise brain volumes, smoothed them with a 4 mm FWHM kernel to match the spatial smoothness of our real data, and extracted the corresponding voxels to each ROI. For each ROI, we computed the Pearson correlation between the fraction map and each smoothed noise map, generating a null distribution of correlation values.

Statistical significance was then calculated by comparing the observed Pearson correlation against this spatially preserved null distribution. This allowed us to test whether the observed correlation between fractions and integers was greater than chance while accounting for the spatial smoothness of the data. We then applied the same procedure to assess the correlation between the two runs for integers and fractions. To plot the tuning data across fractions and integers, we selected the voxels corresponding to numbers from 0 to 3, as there were too few voxels representing magnitudes of 4 and 5. Voxels were assigned to each category based on their rounded preferred magnitude (e.g., voxels categorized as 0 had preferred values between 0 and 0.49). We then extracted the normalized activation for fractions from these voxels (see S7 Fig).

To control for the cumulative risk of false positives across ROIs, each permutation analysis was corrected for multiple comparisons across ROIs using the Benjamini–Hochberg FDR test.

Behavioral evidence for a distance effect and a mental number line

Participants responded significantly more slowly to fractions (M_Frac-Frac = 1,192 ms; SD_Frac-Frac = 129) than to integers (M_Int-Int = 808 ms; SD_Int-Int = 64; t(28) = −17.12, p < 0.001), and they made more errors when both the previous and current numbers were fractions (M_Frac-Frac = 9.5%; SD_Frac-Frac = 6.90) than when both were integers (M_Int-Int = 3.7%, SD _Int-Int = 3.56; Mann-Whitney U = 368, z = 3.77, p < 0.001). In trials requiring comparisons between integers and fractions, participants were also slower when integers preceded fractions (M_Int-Frac = 1,053 ms; SD_Int-Frac = 191) than when fractions preceded integers (M_Frac-Int = 952 ms; SD_Frac-Int = 147, Mann−Whitney U = 469, z = 2.14, p = 0.03), with no corresponding differences in accuracy between the two conditions (M_Int-Frac = 8.11% of errors; SD_Int-Frac = 8.11; M_Frac-Int = 5.56% of errors; SD_Frac-Int = 5.40; Mann−Whitney U = 329, z = 1.23, p = 0.22). These results were stable across time, since reaction times were highly correlated between run 1 and run 2 (r = 0.95, p < 0.0001), as were the percentage of errors (rho = 0.49, p < 0.0001). Overall, these results indicate that participants found conditions involving fractions more difficult than those involving integers.

RTs and the percentage of errors are presented as a function of absolute distance in Fig 1B and 1C for each condition (i.e., Frac-Frac, Frac-Int, Int-Frac, and Int-Int). Regardless of condition, we found a statistically significant distance effect: participants were slower and less accurate when two numbers were closer on the number line. This effect was particularly pronounced when both the current and previous numbers were fractions (RTs: β_Frac-Frac = −220.26, t(28) = −6.54, p < 0.0001; Errors: β_Frac-Frac = −11.54, t(28) = −6.26, p < 0.0001; βs obtained with a linear regression model with log-transformed numerical distance as predictor), but a distance effect was also present when both numbers were integers (RTs: β_Int-Int = −146.68, t(28) = −3.60, p = 0.001; Errors: β_Int-Int = −6.65, t(28) = −2.75, p = 0.01) or when participants compared an integer to a previous fraction (RTs: β_Frac-Int = −252.01, t(34) = −5.12, p < 0.0001; Errors: β_Frac-Int = −9.69, t(34) = −5.54, p < 0.0001) and vice-versa (RTs: β_Int-Frac = −308.79, t(34) = −4.61, p < 0.0001; Errors: β_Int-Frac = −12.67, t(34) = −4.37, p < 0.0001).

To combine information from RTs and error rates, we calculated the BIS score for each pair of numbers [71]; see Fig 1D; see S1 Fig for an alternative measure, Efficiency Score). Values closer to 0 indicate more difficult trials. A participant level linear mixed-effects model revealed a logarithmic distance effect between the previous and current numbers (β = 0.79, t(2447) = 3.77, p = 0.0002). Additionally, the BIS score decreased when the previous item was a fraction compared to an integer (β = −0.38, t(2447) = −4, p < 0.0001), and also when the current item was a fraction compared to an integer (β = −0.68, t(2447) = −7.04, p < 0.0001). The only statistically significant interaction was between the previous and current number conditions, indicating further difficulty when both items were fractions (β = 0.27, t(2447) = 2.17, p = 0.03). Model comparison confirmed that numerical distance was best modeled by the logarithm of numerical distance rather than by a linear function (Δ AIC = 6), reflecting a classical tapering of the distance effect when distance increases.

MDS of the balanced integration score (Fig 1E) showed a one-dimensional curvilinear structure (number line) embedded in two dimensions, similar to that previously observed in similar magnitude comparison tasks [30,76]. Importantly, fractions and integers were integrated together, and both were ordered according to their magnitudes on the first dimension (x axis), with smaller numbers on the left side and larger numbers on the right side. When we correlated x-coordinates of the MDS solution with numerical magnitude using a Spearman correlation (to account for curved or otherwise nonlinear representations), we found a statistically significant correlation (rho(10) = 0.75, p = 0.007). The second dimension (y axis) seemed to encode difficulty, as it separated the intrinsically easier numbers 0 and 5 (because they lie at the extremes) from other numbers. To obtain a measure of difficulty, we calculated the average value of each column in the behavioral BIS matrix, where higher BIS values indicate easier comparisons. A spearman correlation between y-coordinates of the MDS solution and this difficulty measure revealed a significant relationship (rho(10) = 0.64, p = 0.03). Similar results, but perhaps noisier, were found with the Efficiency score (S1 Fig). When we recomputed the MDS excluding 0 and 5, their number line organization was largely unchanged (S2 Fig), further supporting the role of magnitude as a major organizing dimension of semantic space.

Univariate fMRI activations converge to a classic network of math-responsive areas

Three univariate contrasts converged to a largely common set of classic math-responsive regions (Fig 2). First, a contrast for greater activation to math-related sentences than to nonmath-related sentences in the localizer revealed an extensive bilateral network, including the full length of the IPS, posterior inferior temporal gyri (pITG), precentral gyri, middle frontal gyri (MFG), and precuneus. The opposite contrast (in blue), indicating a greater response to nonmath-related sentences, activated a classic and bilateral network of semantic regions, comprising the angular gyri, superior and middle temporal gyri, temporal poles, fusiform gyri, superior and inferior frontal gyri, and posterior cingulate (see Fig 2A). Second, in the numbers task, most of the above math-responsive network showed a greater activation for fractions compared to integers in bilateral regions, including the IPS, precentral gyri, inferior and middle frontal gyri, precuneus, pITG, now also extending to lateral fusiform and lateral occipital cortex (presumably due to the higher visual complexity of fractions; see Fig 2B). Third, similar regions, including bilateral IPS and pITG showed a distance effect during number comparison, i.e., a parametric modulation by the log of the absolute distance between the current and previous numbers, with greater activation to smaller distances, as previously described (e.g., [26], see Fig 2C).

Vector-based analyses of the neural representation of numbers

Using a GLM with one regressor per target, we obtained, within each participant, a reliable high-resolution estimate of the fMRI activation pattern or “vector” evoked by each of 12 individual number concepts. Our aim was to probe whether these vectors organized as a semantic map and what were its main axes of variations.

As a first step, we performed a whole-brain principal component analysis (PCA) of the average vectors across participants within a standard anatomical space (MNI). PCA revealed two main orthogonal components that explained most of the variance in our data: PC1 reflected a categorical distinction between fractions and integers, while PC2 reflected a near-perfect ordering by number magnitude, from small to large, distinguishing smaller and larger numbers while pooling together fractions and integers (Fig 3). In individual analyses, the first component was observed in 14 of our 20 participants, whereas the second was present in 13. Each component produced distinct spatial maps, although some parietal regions were involved in more than one component. Specifically, PC1 confirmed the univariate analysis of Fig 2B, indicated that fractions (probably because of their greater visual and/or semantic complexity), mobilized bilateral occipital, pITG, posterior IPS, and MFG. Crucially, PC2 revealed that a magnitude map mixing fractions and integers was present in much of the math-responsive network (Fig 3B, bottom), with a particularly strong small-to-large gradient in the anterior IPS (extending into the inferior parietal lobule for small numbers), with additional contributions of the MFG, pITG, posterior cingulate, parahippocampal cortex and hippocampus. Together, PC1 and PC2 explained 36.23% of the variance. To determine the unique contribution of domain (fractions versus integers) and numerical magnitude on the brain, we conducted separate multiple linear regressions for the first and second principal components. For the first principal component, the model fit was highly significant (F(2, 9) = 25.47, p = 0.0002, R² = 0.82), and domain was a significant predictor (b = 132.42, t(9) = 6, p = 0.0002), whereas magnitude was not (p = 0.74). For the second principal component, the model fit was also significant (F(2, 9) = 11.98, p = 0.003, R² = 0.67) but revealed the opposite effect: magnitude significantly predicted the second component (b = −28.17, t(9) = −4.82, p = 0.0009), whereas domain did not (p = 0.13). Thus, this analysis confirms that number magnitude is a major determinant of number representations, which affects a broad set of brain areas. Importantly, these findings remained unchanged when the analysis was repeated after excluding the larger integers 4 and 5, so that integers and fractions covered a similar numerical range (see S3 Fig).

To take advantage of the high-resolution individual activations available using 7T fMRI, we extracted single-participant vectors within the entire math-responsive network by first identifying subject-specific voxels showing greater activation to math than to nonmath sentences in the localizer (p < 0.001 uncorrected).

Next, we focused on the six main areas identified by Amalric and colleagues [53] (left and right IPS, pITG, and MFG). In those voxels, we extracted the activation vectors to each of the 12 target numbers in the main comparison task, then visualized those vector spaces using RDMs followed by MDS.

This analysis confirmed that neural vector spaces were organized by two dimensions: a main distinction between fractions and integers, and a gradient of continuous magnitude (Fig 4). The latter was seen in two ways: in the RDMs, neural vectors became increasingly similar as the numbers were closer, i.e., for cells closer to the diagonal of the matrix (for both integers and fractions); and in the 2-dimensional MDS solutions of those RDMs, integers and fractions, while separated, were similarly ordered by magnitude. Furthermore, information was represented similarly across the three math-responsive regions IPS, pITG, and MFG.

Importantly, these results remained stable even when the analysis was restricted to within-domain trials (i.e., integers preceded by integers, and fractions preceded by fractions; see S5 Fig). Moreover, this separation allowed us to assess the reliability of our brain-activity measures: there was a high correlation between the brain activity patterns measured on within-domain trials and on cross-domain trials (see S6 Fig).

Since the hippocampus and parahippocampal cortex were also identified in the PCA, and since the human hippocampus and mesial temporal lobe have been shown to contain number-tuned neurons [18,29,50], we also analyzed vector spaces in those regions. We used the bilateral ROIs defined by [77] and extracted the entire hippocampus and parahippocampal cortex. In these regions, an organization according to magnitude was even more prominent, with MDS dimension 1 ordering the targets according to their approximate magnitude regardless of whether they were fractions or integers (although a distinction between the two domains remained visible in MDS dimension 2). Note that, similar to the behavioral data, the integer number line seemed to bend after three in all MDS representations, which may reflect the absence of fractions in this range (see Fig 5). Importantly, this effect was not driven solely by the presence of “easy” extreme numbers, as the pattern remained consistent when we recomputed the MDS excluding the extremes 0 and 5 (S2 Fig).

To explicitly test for a shared representation of fractions and integers, we modeled the RDMs using multiple regression analysis focusing on the distance between fractions and integers, but also including 4 regressors controlling for categorical effects and for within-domain distance effects (Fig 6, top). In each region of the math-responsive network (i.e., IPS, ITG, and MFG) as well as medial temporal lobes, a significant effect of fraction-to-integer distance on neural dissimilarity was found (S1 Table). Thus, math-responsive regions were not only organized according to within-domain distances between two integers or two fractions (e.g., is more similar to than to ), but also by the between-domain distances (e.g., is more similar to 0 than to 4), indicating that in our educated adult participant, the magnitude of fractions was well integrated with that of integers. A searchlight RSA analysis within the math-responsive areas (Fig 6, bottom) further showed that the regions sensitive to cross-domain distances included the bilateral hippocampus and parahippocampal cortex, IPS, pITG, MFG as well as precentral cortex, posterior cingulate cortex, and precuneus.

Individual voxel tuning and numerotopic map

The similarities observed so far between neural activation vectors could, in principle, arise either from distributed representations or from sparse localized neural responses. Furthermore, in the latter case, cortical patches responding to different concepts could either be distributed randomly on the cortical surface or form an overall mesoscopic topographic map [3,5]. Indeed, in the number domain, previous research showed that the firing of individual neurons can be fitted by a Gaussian tuning curve, reflecting their preference for a specific quantity [18,29,50,78]. Furthermore, those preferences cluster on the cortical surface to form “numerotopic” maps that are detectable by fMRI [7,28,47]. Inspired by this prior work, we fitted the present single-voxel fMRI data with Gaussian tuning curves in order to characterize how numerical information is spatially represented and whether a spatial map exists in which voxels corresponding to similar magnitudes are located closer together, regardless of domain (integers or fractions).

We extracted the single-voxel group-average responses to all numbers within the math-responsive areas defined by the localizer. Each voxel’s response profile was fitted by a Gaussian, separately for integers and fractions and for each run, keeping only the voxels that showed greater than 50% variance explained. The presence of reliable voxel tunings was indicated by the fact that, across voxels, the numerical preferences computed in run 1 and run 2 were significantly and positively correlated (r = 0.28, p corr = 0.0001 for integers, and r = 0.53, p corr = 0.0001 for fractions). There was also a smaller but still reliable correlation between preferred fraction and preferred integer (r = 0.09, p corr = 0.004). For brain visualization, in Fig 7A and 7B, as well as S7 Fig, we used the two runs. The graphs of Fig 7A and 7B show the average cross-validated voxel response profiles, after pooling across all the aforementioned voxels (voxel preferences were determined by fitting the data from one run, and then the normalized activation was plotted from the other run). A permutation statistical test confirmed that the preferred number identified in one run indeed led to significantly stronger activation to that target number than to others in the other run (for integers p corr = 0.013, and for fractions p corr = 0.017).

For integers, a clear distance effect can be seen in Fig 7A, with activation decreasing as the numerical distance to the preferred number increased. This effect was prominent for voxels tuned to numbers 0, 1, 2, and to a lesser extent 3, but became undetectable for numbers 4 and 5, probably because very few voxels were tuned to the latter numbers. For fractions in Fig 7B, the distance effect was less evident, perhaps because magnitudes were closer, but a reliable distinction was observed between voxels tuned to fractions, those tuned to fraction , and those tuned to fractions larger than 1. No such reliable tuning was observed for integers in the hippocampus (p corr = 0.62) and parahippocampal cortex (p corr = 0.11), though a statistical effect was found for fractions in the hippocampus (r = 0.49, p corr = 0.0001) but not for the parahippocampal cortex (p corr = 0.51).

When analyzing separately each of the three main math-responsive regions through a Monte Carlo simulation approach, only the IPS showed a positive and statistically significant correlation between fraction and integer magnitudes (integers: r = 0.23, p corr = 0.0001; fractions: r = 0.80, p corr = 0.0001; across the two domains r = 0.28, p corr = 0.0003). Although within-domain correlations were significant in MFG and pITG (MFG: integers: r = 0.53, p corr = 0.0001; fractions: r = 0.49, p corr = 0.0001; pITG: integers: r = 0.31, p corr = 0.0001; fractions: r = 0.58, p corr = 0.0001), cross-domain correlations were not (MFG: r = −0.26; pITG: r = −0.11).

Within the parietal lobe, a continuous cortical map from larger to smaller numbers was observed, overlapping with areas PF and IP2 of the HCP atlas. As shown in Fig 7C and 7D, the spatial patterns were broadly similar for fractions and integers, with smaller numbers represented in the most anterior part of the IPS, extending into the inferior parietal lobule, and increasingly larger numbers represented increasingly posteriorly in the horizontal part of the intraparietal sulcus. In these areas, consistent tuning was observed both within domain across the two independent runs (integers: r = 0.44, p corr = 0.0001; fractions: r = 0.91, p corr = 0.0001), as well as across the two domains (r = 0.69, p corr = 0.0003). Furthermore, a topographic map was confirmed by correlations between voxel tuning preferences and the Z (inferior-superior) and X (left-right) coordinates of the MNI atlas, both for integers (Z: left r = 0.21, p corr = 0.06; right r = 0.41, p corr = 0.003; X: left r = 0.61, p corr = 0.001; right r = −0.31, p corr = 0.03) and for fractions (Z: left r = 0.83, p corr = 0.0001; right r = 0.70, p corr = 0.0001; X: left r = 0.56, p corr = 0.0001; right r = −0.58, p corr = 0.0001). In contrast, correlations with Y coordinate (anterior-posterior direction) were weaker (Y: left r = −0.34, p corr = 0.004; right r = 0.08, p corr = 0.7) and fractions (Y: left r = −0.36, p corr = 0.002; right r = −0.30, p corr = 0.01).

Finally, as a further proof of cross-domain consistency, when selecting voxels preferring numbers 0, 1, 2, or 3 in the integer map, we found that these voxels exhibited magnitude preferences consistent with fractional quantities. For example, voxels most responsive to 0 in the integer map also showed a higher activation for the fractions , , and (S7C Fig). Importantly, this topographic organization remained evident in the aIPS/IPL even when the analysis was restricted to within-domain number comparison trials (i.e., integers preceded by integers, and fractions preceded by fractions; see S8 Fig).

The above results were obtained at the group level. At the single-subject level, although there was considerable variability across participants, the same general topographic pattern in the aIPS/IPL was observed in 14 participants for integers and 15 participants for fractions (at R² > 0.3). Maps for fractions and integers from four representative participants are shown in S9 and S10 Figs shows the topographic progression for each hemisphere and map.

S1 Fig. Results of an alternative method of measuring task difficulty, efficiency score, calculated by combining information about RTs and the percentage of errors (Accuracy/ RTs) for each pair of numbers.

A) Efficiency score matrix, where smaller values (closer to purple) indicate greater difficulty in comparing numbers, while larger values (yellow) indicate easier comparisons. B) When we correlated the x-coordinates of the MDS solution with numerical magnitude using a Spearman correlation (to account for the curved trajectory of the representation), we found a statistically significant correlation. (rho(10) = 0.86, p = 0.0004). A Spearman correlation between y-coordinates of the MDS solution and the Efficiency score difficulty measure (defined as the average value of each column) also revealed a similar tendency that did not reach statistical significance (rho(10) = −0.48, p = 0.11). The data underlying this Figure can be found in https://doi.org/10.17605/OSF.IO/JQCNA.

https://doi.org/10.1371/journal.pbio.3003820.s001

(TIF)

S2 Fig. Two-dimensional MDS representation obtained without using the extreme numbers (0 and 5).

MDS solutions were obtained using either the behavioral balanced integration score, efficiency score matrix, math-responsive areas, as well as IPS, pITG, MFG, parahippocampal, and hippocampus.

https://doi.org/10.1371/journal.pbio.3003820.s002

(TIF)

S3 Fig. Whole-brain principal component analysis without using larger numbers (i.e., 4 and 5) to have similar ranges to integers and fractions.

The first principal component distinguishes fractions and integers, whereas the second component distinguishes smaller and larger numbers. To determine the unique contribution of domain (fractions versus integers) and numerical magnitude on the brain, we conducted separate multiple linear regressions for the first and second principal components. The model predicting the first component was highly significant (F(2, 7) = 17.15, p = 0.002, R² = 0.78), and the domain was a significant predictor (b = 138.96, t(7) = 5.84, p = 0.0006), whereas magnitude was not (p = 0.27). The model predicting the second component was also significant (F(2, 7) = 31.95, p = 0.003, R² = 0.87) but revealed the opposite effect. Magnitude significantly predicted the second component (b = −45.22, t(7) = −7.59, p = 0.0001), whereas domain not (p = 0.74). The data underlying this Figure can be found in https://doi.org/10.17605/OSF.IO/JQCNA. Abbreviations: SFG, superior frontal gyrus; MFG, middle frontal gyrus; IFG, inferior frontal gyrus; SMA, supplementary motor area; PCC, posterior cingulate cortex; ITG, inferior temporal gyrus; MTG, middle temporal gyrus.

https://doi.org/10.1371/journal.pbio.3003820.s003

(TIF)

S1 Table. Effect of fraction-to-integer distance on neural dissimilarity matrices within the regions of interest used in the RSA (S4 and S5 Figs).

For each region, a multiple regression of RSA matrices was performed as described in Fig 6. The table reports the mean beta coefficient for the crucial fraction-to-integer distance predictor (averaged across subjects, with standard error of the mean in brackets), the peak MNI coordinates of the region, and the corresponding Z-score (derived from a group-level one-sample t test of the subject betas against zero). P-values were corrected for multiple comparisons using Benjamini–Hochberg FDR.

https://doi.org/10.1371/journal.pbio.3003820.s005

(DOCX)

S5 Fig. Visualization of neural representational spaces for numbers within-domain (i.e., trials in which fractions were preceded by fractions and integers were preceded by integers).

Representational similarity matrices (RDMs) and their 2-dimensional multidimensional scale (MDS) solutions were used to visualize the 12-dimensional vector space spanned by number activations within individual voxels-of-interest identified by the math>nonmath contrast in the sentence localizer. A) On the left, the RDM shows the average dissimilarity across participants, computed using only within-domain trials. In the middle, voxel overlap across individuals is displayed with lighter colors indicating higher overlap (up to a maximum of 15 participants). On the right is the corresponding neural geometry representation derived from MDS. B) Results for the intraparietal sulcus (IPS), including average RDM (on the top); the bottom panel shows participant overlap and the resulting neural geometry. C) Analysis of the middle frontal gyrus (MFG). D) Analyses for the posterior inferior temporal gyrus (pITG). E, F) Corresponding analyses for the hippocampal and parahippocampal cortex ROIs, respectively. These regions were defined using Glasser ROIs.

https://doi.org/10.1371/journal.pbio.3003820.s006

(TIF)

S6 Fig. Reliability of representational similarity matrices (RDMs) estimated from independent subset of trials.

RDMs visualize the 24-dimensional vector space spanned by number-evoked activity within voxels-of-interest identified by the math>nonmath contrast in the sentence localizer (p < 0.001 unc, at the individual level). To assess reliability, the data were split according to the preceding stimulus, yielding independent estimates for fractions preceded by either a fraction (“f_”) or an integer (“i_”; top left, black-highlighted region), and for integers preceded by either a fraction or an integer (bottom right black-highlighted region). Despite this split, the resulting RDMs (blue dashed squares) preserved two key signatures: higher similarity for repeated presentations of the same target than for different targets (diagonal > nondiagonal) and a graded decrease in similarity with increasing numerical distance. These effects were observed for both fractions (diagonal > nondiagonal similarity: t(19) = −6.58, p < 0.0001; correlation with numerical distance for nondiagonal terms: t(19) = 4.06, p = 0.0003) and integers (diagonal > nondiagonal similarity: t(19) = −11.28, p < 0.0001; numerical distance: t(19) = 6.93, p < 0.0001). Cross-domain similarity patterns (blue dashed squares) were also highly reproducible across independent subset of trials. For integers, the cross-domain similarity matrices were strongly correlated whether the preceding stimulus was an integer or a fraction (Pearson’s r(13) = 0.84, p corr = 0.0002; Pearson’s r(13) = 0.75, p corr = 0.002, respectively). For fractions, the cross-domain similarity matrices (top left, blue dashed squares) were highly correlated whether the preceding stimulus was a fraction or an integer (Pearson’s r(13) = 0.66, p corr = 0.007; Pearson’s r(13) = 0.93, p corr < 0.0001, respectively).

https://doi.org/10.1371/journal.pbio.3003820.s007

(TIF)

S7 Fig. Further evidence for voxel tuning to specific numbers.

A) Whole-brain view of the voxel tuning preferences for integers and fractions, with all math-responsive regions in the sentence localizer. Overall, the correlation between the two maps is very low (r = 0.09), but it is high and significant in the aIPS/IPL (r = 0.69) (panel B). Furthermore, in the aIPS/IPL, voxels preferring 0, 1, 2, and 3, selected based on their preferences in the integer map, show systematic changes in their normalized activation in response to fractions (panel C). The data underlying this Figure can be found in https://doi.org/10.17605/OSF.IO/JQCNA.

https://doi.org/10.1371/journal.pbio.3003820.s008

(TIF)

S8 Fig. Surface view and coronal slices showing voxel-preferences for integers preceded by integers (left) and fractions preceded by fractions (right) across several voxels within anterior IPS and IPL (R² > 0.5).

This region comprises a congruent numerotopic map across integers and fractions, in which larger numbers are represented with bluish hues and smaller numbers with reddish hues. The map is shown using data averaged across participants. C) Correlation between the map of fractions and integers (r = 0.29). D) In the aIPS/IPL, voxels preferring 0, 1, 2, and 3, selected based on their preferences in the integer map (using only within-domain trials), show systematic changes in their normalized activation in response to fractions (within-domain). The data underlying this Figure can be found in https://doi.org/10.17605/OSF.IO/JQCNA.

https://doi.org/10.1371/journal.pbio.3003820.s009

(TIF)

S10 Fig. Topographic progression along the aIPS/IPL ROI for the participants 2, 8, 18, and 22 (see S9 Fig).

For each participant, we plotted the topographic progression for each hemisphere and each map (integers or fractions). We extracted the preferred numerical magnitude as a function of cortical distance (in mm), binning magnitudes every 5 mm, and the bars show the standard error mean. Cortical distance was computed geodesically starting from the vertex with the lowest preferred magnitude. Then, a spearman rank correlation was performed between distance and preferred magnitude, using permutation testing to assess statistical significance and applying Benjamini–Hochberg FDR correction for multiple comparisons across the 16 tests. After FDR correction, a statistically significant spatial gradient was confirmed in 14 of 16 individual trajectories (with the vast majority surviving at p corr < 0.01). LH, left hemisphere; RH, right hemisphere.

https://doi.org/10.1371/journal.pbio.3003820.s011

(TIF)