## Abstract

Neural representations of approximate numerical value, or numerosity, have been observed in the intraparietal sulcus (IPS) in monkeys and humans, including children. Using functional magnetic resonance imaging, we show that children as young as 3–4 years old exhibit neural tuning to cardinal numerosities in the IPS and that their neural responses are accounted for by a model of numerosity coding that has been used to explain neural responses in the adult IPS. We also found that the sensitivity of children's neural tuning to number in the right IPS was comparable to their numerical discrimination sensitivity observed behaviorally, outside of the scanner. Children's neural tuning curves in the right IPS were significantly sharper than in the left IPS, indicating that numerical representations are more precise and mature more rapidly in the right hemisphere than in the left. Further, we show that children's perceptual sensitivity to numerosity can be predicted by the development of their neural sensitivity to numerosity. This research provides novel evidence of developmental continuity in the neural code underlying numerical representation and demonstrates that children's neural sensitivity to numerosity is related to their cognitive development.

**SIGNIFICANCE STATEMENT** Here we test for the existence of neural tuning to numerosity in the developing brain in the youngest sample of children tested with fMRI to date. Although previous research shows evidence of numerical distance effects in the intraparietal sulcus of the developing brain, those effects could be explained by patterns of neural activity that do not represent neural tuning to numerosity. These data provide the first robust evidence that from as early as 3–4 years of age there is developmental continuity in how the intraparietal sulcus represents the values of numerosities. Moreover, the study goes beyond previous research by examining the relation between neural tuning and perceptual tuning in children.

## Introduction

Numerosity is a fundamental concept that enables humans to estimate, without counting, the number of items in a set. Within the first year of life, infants can represent the numerosity of a set of objects (Jordan and Brannon, 2006a; Izard et al., 2009). However, infants' abilities to discriminate numerosities are limited by the ratio of numerical change between sets (Weber's law; Xu and Spelke, 2000; Lipton and Spelke, 2003). This ratio-dependent numerical discrimination ability is maintained throughout childhood into adulthood and gradually improves with age (Halberda and Feigenson, 2008).

The human ability to discriminate numerosities is thought to be evolutionarily primitive. Nonhuman animals discriminate numerosities with a ratio-dependent precision comparable to that of humans (Cantlon and Brannon, 2006, 2007; Jordan and Brannon, 2006b). In monkeys, single neurons exhibit numerically modulated tuning functions with a Gaussian shape and a width proportional to the value they represent (Sawamura et al., 2002; Nieder and Miller, 2003, 2004). In the intraparietal sulcus (IPS), these neurons are tuned to cardinal numerosities such that their responses peak for their preferred numerosity and decrease as the numerical quantity becomes more disparate from the preferred numerosity. Directly comparing neural and behavioral discrimination functions from monkeys showed similar sensitivity to numerical differences, suggesting that perceptual and neural representations of number are directly related (Nieder and Miller, 2003).

The data from monkeys are consistent with the predictions of a summation coding neural network model proposed to explain neural representations of numerosities in humans (Dehaene and Changeux, 1993). Dehaene and colleagues extended this model to predict the neural representation of numbers in adult humans during functional magnetic resonance imaging (fMRI) adaptation (Piazza et al., 2004). They predicted that the collective neural responses to numerosity changes in an array should depend on the ratio of numerical change following a peaked tuning function. Importantly, because this model was tested with fMRI adaptation, the amount of recovery to numerically deviant stimuli was predicted to follow an inverted Gaussian shape across ratio. In fact, the pattern of recovery activation to numerical change was fit by the Gaussian shape predicted by the numerosity detector model. These tuning curves were found in the IPS, which is consistent with single neuron data from monkeys (Sawamura et al., 2002; Nieder and Miller, 2004). Several studies have reported evidence of numerosity tuning in the adult IPS not only for representations of sets (Piazza et al., 2004) but also for numerals and count words (Notebaert et al., 2011; Roggeman et al., 2011; Holloway et al., 2013; Demeyere et al., 2014) and even fractions (Jacob and Nieder, 2009). Together these studies show that intraparietal tuning curves related to numerical representation are widespread across tasks in adults.

Neural tuning curves of numerosity representations have never been shown in children. The only evidence of numerical distance effects comes from older children (6+ years; Ansari and Dhital, 2006; Cantlon et al., 2009; Holloway and Ansari, 2010; Bugden et al., 2012). However, numerical distance effects could arise from neural patterns that do not conform to a tuning curve. Although previous neuroimaging work finds sensitivity to numerical over shape processing in posterior regions of the infant brain (Izard et al., 2008; Libertus et al., 2009; Hyde et al., 2010; Edwards et al., 2016) and in the IPS of 4-year-old children (Cantlon et al., 2006), those regions cannot be said to represent cardinal values of numerosities if they do not show numerosity tuning. Understanding whether IPS neurons show Gaussian tuning to numerosity in young children is important for determining whether cardinal representations of numerosities are represented similarly throughout development (Nieder and Miller, 2003; Diester and Nieder, 2007; Nieder, 2009; Cantlon, 2012). If neural tuning to numerosity is a fundamental process in the human IPS then (1) this process should emerge in early childhood and (2) neural numerosity tuning curves should parallel perceptual numerosity tuning curves. The precision of children's neural tuning to numerosity may also predict their perceptual sensitivity to numerosity.

## Materials and Methods

##### Participants.

Forty-three typically developing children (3.6–6.99 years old; mean age, 5.36 years; 27 females) participated in this study. All participants had normal or corrected to normal vision and no history of neurological impairments. The parents of all children gave informed written consent in accordance with the University of Rochester's Research Subjects Review Board. Children were excluded from the analyses due to failure to complete at least one run of the adaptation paradigm (*n* = 4), excessive button pressing during the adaptation scans (*n* = 2; responding >12% of the trials or >4 times more than expected), a speech impediment (*n* = 1), and chance performance at all levels of the numerical discrimination task (*n* = 1; see below, Behavioral assessment). This resulted in a final sample size of 35 children (3.6–6.99 years old; mean age, 5.45 years; 3 years: *n* = 3; 4 years: *n* = 10; 5 years: *n* = 11; 6 years: *n* = 11; 22 females). Among those 35 children, the maximum amount of movement averaged across translation and rotation for an individual child was <4 mm (3.86 mm; range, 0.84–3.86 mm; for motion correction methods, see below, Preprocessing).

##### fMRI session.

Before scanning, children were familiarized with the MRI environment and practiced lying still in a mock scanner. Following the 30 min practice session, children proceeded to the actual MR scanner where their heads were secured with headphones, foam padding, and medical tape.

During the MR scanning session, neural activity (BOLD signal) was measured during one (*n* = 9) or two runs (8.4–9 min each) of an adaptation paradigm. One high-resolution anatomical scan was also collected from each participant.

##### fMRI adaptation paradigm.

Dot arrays were presented on a black background in the top center of the viewing screen. Standard dot arrays consisted of 16 purple dots (RGB values: 128, 0, 128) each with an area of 500 pixels^{2}. Location and density of the dots varied continuously. Standard densities randomly varied from 0.00017 to 0.0003 dots per area [in pixels squared (px^{2})] of the screen in which the dots were placed. Standard arrays were presented for 6–8 trials to allow for sufficient neural adaptation before the deviant dot array. Deviant arrays varied by one of three dimensions: dot area, number of dots, or dot color. Dot color was used as a control condition to disentangle recovery to magnitude of the dots (number and area) from recovery to any sort of perceptual change. Only one dimension of interest varied at a time to ensure that the resulting neural changes could only be attributed to the dimension of interest. Importantly, this means that although cumulative surface area of the dot array varied for both numerical and area deviants, the cumulative surface area changes were matched across the numerosity and area conditions, allowing us to attribute the differences in response to a particular dimension. Dot area was represented using the area of the dot in pixels squared. Color was represented as a binary system in which each of the three RGB values was either 0 or 128. Area and number of dots varied using large (1:2) and small (2:3) ratios and their inverses, while color varied based on the number of changed RGB values, resulting in four alternatives for each deviant dimension (Table 1). Density for the deviant trials was calculated by using the average density from the preceding adaptation trials. Each deviant alternative was presented four times throughout a run, resulting in 48 deviant trials per run. See Figure 1 for an example adaptation period and examples of deviant arrays.

Because the adaptation trials were consistent in all aspects except for density and location, regions that process density and location would not adapt to the standard arrays, and the resulting neural changes in response to the deviants must be attributed to the deviant dimension (surface area of dots, number of dots, or color of dots). All deviant stimuli also had the average density of the adaptation stimuli from the preceding trials and thus there should be no increase in activity to deviants based on density. Although the cumulative area of the dots changed for both the area and number deviants, the number of dots in the area deviant stimuli was the same as in the adaptation stimuli, which enabled us to separate the effects of changes in cumulative area from the effects of changes in number.

All dot arrays were presented for 400 ms with an interstimulus interval of 800 ms. Dot arrays were generated using Psychtoolbox functions in Matlab (Psychtoolbox RRID:SCR_002881; Matlab RRID:SCR_001622). A small smiley-face icon (1.26 cm diameter) was placed in the center of the dot array to serve as a fixation point. On rare occasion (11 trials; 3% of total trials), the smiley face winked. To ensure fixation, participants were instructed to press a button with their right index finger when the smiley face winked. Two to three adaptation trials were added as a buffer following each winking face (≥6 adaptation trials total).

##### MR parameters.

Whole-brain BOLD imaging was conducted on a 3 tesla Siemens Magnetom Trio scanner with a 12-channel head coil at the Rochester Center for Brain Imaging. High-resolution structural T1 contrast images were acquired using a magnetization-prepared rapid gradient echo pulse sequence at the start of each session [TR = 2530 ms; TE = 3.44 ms; flip angle, 7°; FOV, 256 mm; matrix, 256 × 256, 192, or 176 (depending on head size); 1 × 1 × 1 mm sagittal left-to-right slices].

An echo-planar imaging pulse sequence was used for T2* contrast (TR = 2000 ms; TE = 30 ms; flip angle, 90°; FOV, 256 mm; matrix, 64 × 64; 30 axial oblique slices, parallel to the anterior commissure–posterior commissure plane; voxel size, 4 × 4 × 4 mm).

##### Preprocessing.

fMRI data were analyzed using BrainVoyager (Goebel et al., 2006; RRID:SCR_013057) and in-house scripts drawing on the BVQX toolbox (RRID:SCR_009532; http://support.brainvoyager.com/available-tools/52). The first two volumes of functional data in each run were discarded before analysis. Preprocessing consisted of slice scan time correction (cubic spline interpolation), motion correction with respect to the first volume in the first run, and linear trend removal in the temporal domain (cutoff: two cycles within the run). Functional data were registered to high-resolution anatomy on a participant-by-participant basis in native space. Echo-planar and anatomical volumes were then transformed into Talairach space (Talairach and Tournoux, 1988). Data were normalized into Talairach space by first aligning images with the stereotactic axes and then transforming them to the Talairach grid using a piecewise affine transformation based on manual identification of anatomical landmarks. Analyses were performed on preprocessed data in Talairach space. A Gaussian spatial filter was applied to each volume of functional data at 1.5 voxels (6 mm) FWHM. Average framewise displacement (Grill-Spector et al., 2008; Power et al., 2012) was regressed for each subject using their volume-by-volume realignment parameters. This controls for any increases in signal intensity due to volume-to-volume changes in motion.

Functional data were then analyzed using the general linear model. Experimental events were convolved with a standard dual gamma hemodynamic response function. Two models were used to identify the effects of the adaptation period and the recovery to the deviant conditions. The model for adaptation consisted of two regressors of interest corresponding to the main effect of adaptation and a parametric effect of adaptation. The parametric effect of adaptation assigned weights to each trial of adaptation such that the first three trials of a block of adaptation were weighted the highest and all remaining trials received a weight one unit less than the previous trial. Any trial of adaptation during which the smiley face changed to a wink was not included in the adaptation predictor. In addition, the three trials of adaptation following a wink were also excluded. Excluding these trials prevented the neural response to the wink from interfering with the modeling of a monotonically decreasing effect of adaptation. One regressor of no interest was included for the button press and an additional six regressors of no interested corresponded to the motion parameters obtained during preprocessing. In the recovery model, there were 12 regressors of interest that corresponded to the four deviant alternatives for each of the three deviant conditions. This model also included seven regressors of no interest: one for the button press and six for the motion parameters obtained during preprocessing. Any deviant trial that followed a wink was excluded from its corresponding predictor to ensure that the neural response to the wink did not interfere with the deviant response.

##### Analyses of fMRI data.

Random effects analyses were used to analyze the group data. One set of analyses is a whole-brain analysis. The critical analysis for our hypothesis is a conjunction analyses: an analysis that tests the conjunction of adaptation and numerosity recovery effects in children. A second critical analysis for our hypothesis, also a conjunction analysis, tests for independence in patterns of numerosity activation compared with other deviant conditions. For all whole-brain conjunction analyses, the voxel-level threshold is *p* < 0.0025, cluster corrected for multiple comparisons to *p* < 0.05. In some cases, the individual map data are displayed at a fixed cluster threshold of 25 mm^{2} and voxel level threshold of *p* < 0.05 to report the full extent of activation in the individual maps that contributed to the conjunction and affirming that any lack of an effect is not due to a type II error.

A second set of analyses tests whether the model of numerosity encoding previously used by Piazza et al. (2004) to predict neural responses to numerosity changes in adults also accounts for patterns in neural data from young children. These critical tests of our neural tuning hypotheses are region-of-interest (ROI) analyses. The ROIs were defined in two ways: (1) by extracting individual data from the IPS activation in whole-brain group conjunction activation to the adaptation and recovery periods described above and (2) by using a leave-one-out analysis to measure each child's numerosity recovery activation within the adaptation region with an unbiased, independent ROI analysis. Conjunction analyses were cluster corrected by implementing a Monte Carlo simulation of the activation patterns. This procedure was modeled after the BrainVoyager Cluster-Level Statistical Threshold Estimator. However, because the BrainVoyager thresholding procedure cannot simulate the likelihood of a cluster size after taking the overlap of two maps, we implemented our own, similar version, in Matlab using BVQX tools. This enabled us to accurately calculate a cluster threshold for our analyses. First, two images of the same size as the contrast maps of interest (58 × 40 × 46 voxels) were filled with random values along a normal distribution. These images were then spatially smoothed (FWHM, two voxels), and the values in the voxels of each image were scaled using the mean and SD of the contrast values for one of the two original maps (the map of the adaptation effect for one image and the map of the number-preferring effect for the other image). This ensured that the randomly generated images had the same distribution of contrast values as the original maps. Then we took the overlap of the two randomly generated maps and recorded the sizes of clusters significant at a voxelwise threshold of *p* < 0.05. This process was iterated 10,000 times. After 10,000 iterations, we calculated the probability that a cluster of a particular size appeared in a randomly generated overlap map. The cluster threshold was determined by finding the largest cluster size with a probability greater than 0.05 and adding one voxel. This threshold represents the smallest possible cluster size that is unlikely to be due to chance.

Finally, we conducted analyses relating individual differences in neural activation to individual differences in behavior. These analyses were conducted on the ROI data. To ensure an unbiased ROI analysis, we used a leave-one-out procedure described in the Results for all these analyses.

##### Numerosity-encoding model.

The model of numerosity coding that we used to test for evidence of neural tuning to number in children is the model of neural responses to changes in numerosity that predicts adult neural activity (Piazza et al., 2004). This model for adult neural data assumes that the neural signal obtained from fMRI is linearly related to the firing rate of the neurons in a voxel and fits three free parameters: λ, μ, and *w*. In this model (Eq. 1), *I* represents the neural activation as a ratio of the deviant array with respect to the standard array, λ represents the product of the baseline neural signal measured by using fMRI (brain imaging signal) and the underlying baseline amplitude of cell firing, μ represents the baseline brain imaging signal in relation to the rate of habituation and the rate of recovery from habituation, and *w* represents the neural Weber fraction. The λ and μ parameters scale the predicted tuning curve model to the range of the observed data. The neural Weber fraction (*w*) represents the width of the neural tuning curve. The model takes the deviant numerosity and habituated (standard) numerosity as input and then generates the predicted amplitude of recovery activation at each deviant ratio from the starting model parameters. The model calculates the average of the observed neural amplitudes at each numerical pair at each ratio across participants before the model is fit. This step reduces the idiosyncratic noise in the neural (and behavioral) Weber estimates. The best fit for those data are obtained by using the nonlinear squares function (nls, “port” algorithm) in R (version 3.1.1; via R-Studio version 0.98.1073; RRID:SCR_001905) to minimize the nonlinear least-squares estimates of the three parameters with respect to the observed data. Fits between predicted and observed neural data were tested using a traditional goodness-of-fit test. High fits would be expected if the model well explains the data, particularly if the observed data conform to the predicted Gaussian shape and width. The model (Eq. 1) is expressed as follows:
The model and statistics are fully parallel to prior research with adults (Piazza et al., 2004), which makes our results comparable to the adult results, to test for developmental continuity in neural mechanisms.

##### Behavioral assessment.

Following the scan, children completed a numerical discrimination task on a touchscreen computer to measure the acuity of their approximate number system. During the numerical discrimination task, children were shown two ladybugs side by side and were instructed to touch the ladybug with the most dots. The ladybugs' backs were constant in size (125,600 px^{2}) across all trials. The dot arrays presented on the ladybugs' backs varied in location and number of dots from 1 to 30. On half of the trials, the dots on the two ladybugs were the same size (510 px^{2}), and on the other half of the trials the dots were varied in size such that the cummulative surface area of the dots was the same (10,300 px^{2}). Comparisons were classified as having a small, medium, or large number of dots and were made across five different ratios: 0.25, 0.33, 0.5, 0.7, and 0.9. Dot densities (average interitem distance) varied between a large and small density such that half of the time the ladybug with the most dots had the larger density and half the time it had the smaller density. Correct answers were equally presented on the left and right sides of the screen. Numerical discrimination blocks consisted of 30 trials. Four blocks were presented for a total of 120 trials. One child, who successfully completed the adaptation scans, was excluded for lack of motivation during the behavioral session, resulting in a final sample size of 35 children for both the fMRI and behavioral analyses. Individual Weber fractions were derived by fitting model curves to children's behavioral data. Individual curves were fit using the “optim” function (“L-BFGS-B” method; *w* constrained from 0.01 to 2) in R to minimize the sum of squared error.

## Results

### Behavioral results: perceptual tuning to numerosity

To assess the acuity of children's numerical tuning during an active numerical discrimination task, we measured children's performance. Average accuracy across children was 79% (SD, 8%). We derived a group perceptual tuning curve for children's performance on the numerical discrimination task. Children's mean accuracies were calculated at each of the ratios in the numerical discrimination task with respect to the array on the left side of the screen. To fit these averages to a tuning curve and calculate a Weber fraction for the group, we followed the method used by Piazza et al. (2004) for adults. This revealed a numerical tuning curve with a Weber fraction of 0.22 at the group level [*R*^{2} (coefficient of determination) = 0.99, *p* < 0.00001, *n* = 10]. Group perceptual curves were similar across conditions (Fig. 2) with a Weber fraction of 0.20 when the cumulative surface area was equal and 0.25 when the dot size was equal across the two arrays being compared.

We compared the children's behavioral numerosity curve to the behavioral curve obtained for adults who completed a similar task (Piazza et al., 2004). The adult data were rescaled to match the range of our data from children. Data from adults and children were then plotted against each other as shown in Figure 3*A*. Children's performance on the numerical discrimination task follows a similar pattern to that of adults, but exhibits a wider deviation than the adult curve, which is consistent with the conclusion that numerical discrimination improves between childhood and adulthood.

### fMRI results

#### Whole-brain analyses

##### Adaptation effect.

First, we identified regions that showed a monotonically decreasing effect of adaptation (parametric decrease in BOLD signal). This revealed large regions of the parietal cortex, the ventral temporal cortex, and the inferior frontal cortex in both hemispheres (Fig. 4, blue). Thus, a substantial swath of cortex shows a monotonic decrease in neural responses on each trial to the constant number, surface area, and color arrays.

##### Recovery from adaptation.

To determine where in the brain children exhibited neural recovery from adaptation during the presentation of deviant stimuli, we examined neural recovery to the deviant stimuli that presented large changes in number, surface area, and color (ratio changes of 0.5 and 2). We found regions in the bilateral parietal and visual cortex and regions in the left ventral temporal cortex that significantly responded to changes in number (Figs. 4, 5, red), visual and parietal regions in both hemispheres that significantly responded to changes in color (Fig. 5, yellow), and regions in the left parietal cortex that significantly responded to changes in surface area (Fig. 5, green). Furthermore, subregions in both the left and right visual and parietal cortices responded to both color and number changes (Fig. 5, orange) and one region in the left postcentral cortex responded to changes in all three dimensions, though this region did not survive cluster correction (Fig. 5, navy). These patterns resemble those previously reported in studies with adults as “distributed and overlapping” networks for processing information about number, color, and size (Pinel et al., 2004).

##### Number preferences.

To identify regions that responded significantly more strongly to changes in number than to changes in other dimensions, the neural response to numerical changes was compared with the neural response to changes in area and color (number deviants > area and color deviants, balanced). This revealed the number-preferring cortex in the bilateral intraparietal sulcus (Fig. 4, pink).

#### ROI analyses

##### Neural tuning in the IPS.

To determine whether neural tuning to number occurs in children, we defined ROIs in the right and left IPS that showed (1) a monotonically decreasing effect of adaptation and (2) a preference for numerical changes (Fig. 4, green with black outlines; Table 2).

We extracted beta values from the left and the right IPS to calculate children's neural response amplitudes at each of the numerical deviants (0.5, 0.67, 1, 1.5, and 2; the 1.0 ratio represents no change in number and was calculated from the average activation to color and area deviants, which did not vary in number from the adaptation stimuli). As with the behavioral data, these values were tested for their fit to a neural tuning curve following the method used by Piazza et al. (2004; see Materials and Methods). Neural Weber fractions were obtained from the model of children's tuning curves in the right and left IPS. Children's neural Weber fraction in the right IPS was 0.35 (*R*^{2} = 0.97, *p* = 0.003, *n* = 5 ratios), which was similar to children's perceptual Weber fractions from their behavior (behavioral *w* = 0.22; correlation between neural and perceptual tuning curves: *r* = 0.98, *t*_{(3)} = 8.01, *p* = 0.004). In the left IPS, the neural tuning model failed to converge, suggesting that unlike the right IPS, the left IPS is only weakly tuned to numerosity in children (*R*^{2} = 0.87; *p* = 0.02; *n* = 5 ratios). The similarity between the neural tuning curve in the left IPS and the perceptual tuning curve was still high (*r* = 0.94, *t*_{(3)} = 4.85, *p* = 0.017), suggesting some connection to children's numerical perception. The model fit and the relation between children's behavioral and neural tuning, particularly in the right IPS, shows that both functions are explained by a common model of numerosity tuning that predicts Gaussian shape response profiles with similar widths.

Separating the child data by age group (3-year-olds and 4-year-olds vs 5-year-olds and 6-year-olds) revealed that a strongly tuned Gaussian-shaped curve is observed in even the youngest children in the right IPS (3–4-year-olds: *R*^{2} = 0.85, *p* = 0.03, *n* = 5 ratios; 5–6-year-olds: *R*^{2} = 0.97, *p* = 0.003, *n* = 5 ratios) and neural responses were only weakly tuned to numerosity in the left IPS (3–4-year-olds: *R*^{2} = 0.77, *p* = 0.052, *n* = 5 ratios; 5–6-year-olds: *R*^{2} = 0.97, *p* = 0.003, *n* = 5 ratios).

The neural tuning curves for the left and right IPS from children are plotted against the adult data redrawn from Piazza et al. (2004) in Figure 3*B*. The adult data from Piazza et al. (2004) were rescaled to match the amplitude range of the children's data by first standardizing the adult data to a mean of 0 and an SD of 1 and then adding the mean of the children's data and dividing by the SD of the children's data. Children's neural tuning curves and Weber fractions were more similar to the adults' in the right IPS than in the left IPS. These data suggest that neural tuning to number matures faster in the right IPS than in the left IPS.

Next, we examined whether children's tuning curves in the IPS were better explained by a logarithmic scale, as predicted by either approximate linear–scalar or logarithmic encoding of numerical values, versus a precisely linear scale. We tested whether the neural responses to changes in number could be explained by a precise linear version of the model of numerosity encoding (Eq. 2). When the neural data were fit to the precise linear model, the model was found to converge in both the left and the right IPS. However, goodness-of-fit tests for the linear model revealed that it only explained ∼70% of the variance in neural data (right IPS: *R*^{2} = 0.76, *p* = 0.055; left IPS: *R*^{2} = 0.64, *p* = 0.09) compared with ∼90% of the variance explained by the logarithmic model (Fig. 6). Equation 2 is expressed as follows:
Finally, we used a “leave-one-out” approach to independently define the ROIs from which to extract each child's data in order avoid introducing bias between voxel selection and the extracted beta values (Kriegeskorte et al., 2009). We individually defined the regions of the left and right IPS that showed an overlap between adaptation and a preference for number using 34 (*n* − 1) subjects. We then extracted the data from the 35th (left out) subject. This procedure was iterated 35 times, each time leaving a different participant out of the analysis.

The leave-one-out approach resulted in independently defined estimates of neural amplitudes for each subject, at each ratio. Fitting these data to tuning curves based on the response to number validated the result found in the regions of the IPS defined using a group contrast. Specifically, we reproduced the same pattern of Weber fractions in the individually defined leave-one-out ROIs as found in the group-defined ROIs: the right IPS exhibited sharp tuning to number (*w* = 0.23; *R*^{2} = 0.95, *p* = 0.006, *n* = 5 ratios) and the left IPS exhibited immature tuning and did not converge at the group level.

##### Control Analyses.

To ensure that neural tuning was unique to numerosity, we tested whether the IPS regions that showed numerosity tuning curves in children would also show tuning curves across the other deviant dimensions (color and area; Fig. 7). We found neither neural tuning to brightness nor to area in the numerosity-tuned IPS regions as children's neural responses to those dimensions were not modulated by color or area ratio in the numerosity regions (right IPS brightness: *R*^{2} = 0.08, *p* = 0.65, *n* = 5 ratios; left IPS brightness: *R*^{2} = 0.10, *p* = 0.6, *n* = 5 ratios; right IPS area: *R*^{2} = 0.02, *p* = 0.82, *n* = 5 ratios; left IPS area: *R*^{2} = 0, *p* = 1, *n* = 5 ratios; models for brightness and area in the right IPS failed to converge and the model for area in the left IPS resulted in singular convergence). Children's neural tuning to numerosity is thus dimension-specific, rather than a general property of cognitive or neural processing in the IPS. These results are important to show because they provide information about the profile of activation in this specific region. However, they do not represent an independent test of non-numerical tuning because the region was defined by its activation to numerosity. Next, we describe an independent test of the specificity of neural tuning by conducting a whole-brain control analysis of color-preferring cortex. We focus on color-preferring regions because there were substantial color-preferring cortical regions that overlapped number-preferring regions and virtually no area-preferring regions.

As a positive control, we identified regions of the bilateral ventral temporal cortex and the bilateral posterior parietal cortex that showed both an effect of adaptation and a preference for changes in color over number and area. We tested whether those color-preferring regions would show neural tuning to brightness. Fitting curves to the neural response to color deviants revealed strong neural tuning to brightness in all four ROIs (Fig. 8). In contrast, the neural response to numerosity in those color-preferring regions did not converge on a tuning curve function and resulted in very large Weber fractions (numerosity Weber fractions in color-preferring regions: left fusiform: *w* = 2.78; right fusiform: *w* = 2.94; left posterior parietal cortex: *w* = 2.23; right posterior parietal cortex: *w* = 2.15), suggesting poor tuning to numerosity in these regions. These results indicate that dimension-specific neural tuning curves are observed in regions of cortex that show a preference for processing that dimension in very young children. Regions that show these patterns for numerosity tuning are distinct from those of brightness tuning.

### Individual variability

To date, numerosity tuning curves have not been fit to neural data from individual subjects, even in adults, nor have attempts been made to directly relate measures of perceptual and neural tuning at the individual level. Most likely this is because the behavioral output of perceptual tuning is the result of a combination of neural responses throughout the entire brain and is unlikely to be predicted by neural responses from just one small brain region. Additionally, neural data at the individual subject level are noisy. Therefore, it is possible that there is no relation, or only a weak relation, between neural tuning in the IPS and perceptual tuning. Here we present the first evidence of neural tuning to numerosity at the individual level in humans and conduct exploratory analyses to examine the relation between perceptual tuning and neural tuning in individual subjects.

First, we calculated perceptual and neural tuning curves for each subject. Neural tuning curves were fit using the optim algorithm in R (L-BFGS-B method to constrain the parameter space to plausible *w*, λ, and μ based on the group data) and were calculated over the independently defined leave-one-out right and left IPS ROIs. Figure 9 shows each subject's neural and behavioral curves. More children showed evidence of neural tuning to numerosity in the right IPS than in the left IPS, which exhibited more flat curves (for μ> 0; right IPS: mean observed/predicted *r* = 0.48; range, 0.09–0.93; *t* test vs 0: *t*_{(24)} = 10.07, *p* < 0.000001; binomial test *z*-ratio = 2.37, *p* = 0.009; left IPS: mean observed/predicted *r* = 0.55; range, 0.55–0.99; *t* test vs 0: *t*_{(20)} = 10.24, *p* < 0.000001; binomial test *z*-ratio, 1.01; *p* = 0.156). This suggests that neural tuning curves are evident at the individual level as measured using fMRI adaptation, but are more consistent across children in the right IPS than in the left IPS. Across all children, the neural Weber fraction was significantly higher in the left IPS than in the right IPS (*t*_{(34)} = 2.52, *p* = 0.017; mean difference, 0.297), which is in line with the higher neural *W* in the left IPS observed at the group level.

We derived behavioral Weber fractions at the individual level by fitting the model of numerosity encoding to each child's accuracy across numerical ratios. Individual behavioral curves converged for all children (mean observed/predicted *r* = 0.84; range, 0.27–0.99; *t* test vs 0: *t*_{(34)} = 34, *p* < 0.000001). As expected, we found that numerosity sensitivity increases with age: as age increases, children's Weber fractions decrease (*F*_{(1,33)} = 11.94, *R*^{2} = 0.27, *p* = 0.001; age estimate, −0.25, *t* = −3.46, *p* = 0.001). Individual Weber fractions did not differ across trial type (equal cumulative surface area vs equal dot size, *t*_{(34)} = 0.25, *p* = 0.8). Subsequent analyses were performed collapsing across trial type (mean behavioral *w* = 0.25; *w* range, 0.07–0.65) with data from children showing evidence of perceptual tuning at the individual level (μ > 0, *w* < 2).

To evaluate the degree of similarity between neural and behavioral curves at the individual level, we correlated each child's neural tuning curve with his/her behavioral curve. *T* tests versus 0 revealed that these correlations were significant (right IPS: mean *r* = 0.93; *R* range, 0.58–0.99; *t*_{(22)} = 44, *p* < 0.00001; left IPS: mean *r* = 0.88; *R* range, 0.56–0.99; *t*_{(18)} = 26, *p* < 0.00001). This suggests a close relation between children's neural tuning during numerosity adaptation and behavioral discrimination of numerosity.

Next, we conducted exploratory analyses to test whether children's neural responses to numerical deviants predicted individual performance on a numerical task. We correlated children's neural Weber fractions with their perceptual Weber fractions. Perceptual Weber fractions were significantly correlated with neural Weber fractions from the right IPS (*r* = 0.40, *t*_{(21)} = 2.05, *p* = 0.056), but not from the left IPS (*r* = 0.35, *t*_{(17)} = 1.54, *p* = 0.142; Fig. 10). This relation was unique to Weber fractions as age was not related to neural Weber fractions (right IPS: *r* = −0.08, *t*_{(23)} = −0.393, *p* = 0.697; left IPS: *r* = −0.08, *t*_{(19)} = −0.334, *p* = 0.742). Children's average neural Weber fractions from the right and left IPS also related to their perceptual Weber fractions (perceptual *W*: *r* = 0.46, *t*_{(16)} = 2.088, *p* = 0.053; age: *r* = −0.40, *t*_{(16)} = −1.763, *p* = 0.097). Thus, both the average *w* and the right IPS *w* alone predicted behavioral *w* in children. These results reveal a relation between individual differences in neural Weber fractions and perceptual Weber fractions.

## Discussion

To determine whether nonsymbolic numerosity processing in early childhood relies on neural computations similar to those seen in adults, we applied a model of numerosity coding based on adult neural processes to neural data from 3–6-year-old children. We found that the adult model explained children's neural responses to numerosity in the same neural regions observed in adults, suggesting that analog representations of number are a foundational aspect of nonsymbolic numerosity processing. This model of neural tuning was tested separately on the youngest children in our sample, 3–4-year-olds, and accounted for neural responses even in those children alone. Our data provide the most robust evidence to date that the neural mechanisms underlying numerosity perception are stable across human development. Moreover, the results show that individual variability in perceptual tuning is related to individual variability in neural tuning in children. Thus, the IPS can be said to represent the approximate cardinal values of numerosities in young children and is a likely source of individual variability in numerosity perception. Here we discuss the continuities and discontinuities in neural coding between children and adults and their relation to numerical development.

A source of neural continuity in numerosity representations over development is the numerosity code of the right IPS. The right IPS exhibits adaptation to constant numerosity, recovery in neural responding for numerical changes, and a neural preference for numerosity changes compared with changes along other dimensions. The right IPS also exhibits a robust tuning curve to variation in numerical value in children as young as 3–4 years. Although previous studies have demonstrated evidence of neural distance effects in older children (Ansari and Dhital, 2006; Cantlon et al., 2009; Holloway et al., 2010; Bugden et al., 2012; Vogel et al., 2015), neural distance effects could arise from a variety of patterns that do not conform to a tuning curve. The current study is the first to fit a mathematical model of cardinal numerosity tuning to neural data from children. Neural tuning to numerosity in the right IPS in children was comparable to adult-level sensitivity though slightly weaker (*w* = 0.35 in children vs 0.25 in adults), and paralleled the children's current degree of numerical sensitivity in perception at the group level (*w* = 0.22), which was comparable to perceptual Weber fractions previously reported from this age group (Halberda and Feigenson, 2008). Prior research suggests that numerosity is initially represented on a logarithmic scale (approximate) and shifts to a linear (precise) representation with education and formal training with a symbolic system of numbers (Siegler and Opfer, 2003; Opfer and Siegler, 2007; Dehaene et al., 2008; Siegler et al., 2009). Our data show that young children's pattern of neural tuning to numerosity is better explained by a logarithmic representation of numerosities (90% of variance) compared with a precise linear representation (70% of variance). Importantly, the same logarithmic model that accounted for children's perceptual tuning to numerosity also accounted for neural tuning to numerosity in the right IPS, implicating a relation between the development of numerosity perception and the numerosity code of the right IPS.

This study finds an important difference between numerical representation in the right and the left IPS. Unlike adults, whose neural responses exhibited finer sensitivity in the left IPS compared with the right (Piazza et al., 2004), children's neural sensitivity in the left IPS was far cruder than their neural sensitivity in the right IPS. This contrasts with findings from adults where the neural tuning model of numerosity representation accounted for neural responses in both the right and left IPS (Piazza et al., 2004). Hemispheric differences in numerical processing during early childhood could be an important source of developmental change in numerical cognition (Cantlon et al., 2006; Ansari, 2008; Hyde et al., 2010). The observation that the right IPS matures earlier than the left IPS is consistent with previous developmental work showing more adult-like patterns of number-related activation in the right IPS compared with the left (Cantlon et al., 2006; Ansari, 2008; Izard et al., 2008; Libertus et al., 2009; Holloway and Ansari, 2010; Hyde et al., 2010; Cantlon and Li, 2013; Emerson and Cantlon, 2015). Previous work also indicates that activation in the left IPS is more related to symbolic mathematics development than the right IPS (Rivera et al., 2005; Ansari and Dhital, 2006; Rosenberg-Lee et al., 2011; Bugden et al., 2012; Cantlon and Li, 2013; Emerson and Cantlon, 2015; Vogel et al., 2015). The current data from younger children suggest that numerosity representation follows separate developmental trajectories in the right and left hemispheres. Ultimately, in adulthood, the left IPS exhibits sharper tuning than the right IPS (Piazza et al., 2004, 2007). One implication of our data, in light of previous findings, is that numerosity tuning in the left IPS is refined later in development, during exposure to symbolic numerical concepts.

Although the behavioral and right IPS neural data share approximately the same inverted Gaussian shape and width and are explained by a common model, we observed a finer sensitivity in children's behavioral responses than in their neural responses, a pattern that was also seen in previous work with adults (Piazza et al., 2004). In addition, our data show only a moderate relation between behavioral and neural Weber fractions at the individual subject level (*r* = 0.40). These results suggest that neural Weber fractions in the IPS from an adaptation task are not a direct, high-fidelity “read out” of individual perceptual sensitivity. This could be due to the active versus passive nature of the behavioral and fMRI tasks. That is, making explicit numerical judgments in a behavioral task could heighten the sensitivity of the neural mechanisms that underlie numerosity coding beyond what we recorded at the neural level from the passive viewing fMRI task. Alternatively, behavior could be the output of only a high-quality subset of the total neural population measured with fMRI. Both possibilities could result in the slightly greater behavioral sensitivity we observed in children compared with their neural sensitivity and the weak relation between individual subject Weber fractions in behavior and neural activity.

Here, for the first time, we derived individual curves for neural tuning in the right and the left IPS. Previously this had not been done, even in adults, because when Weber fractions are calculated for each individual, those estimates compound idiosyncratic differences in performance between the behavioral and neural tasks, exaggerate the influence of data points that might stray from a Gaussian shape, and generally add noise. Nonetheless, we observed both group-level and individual-level relations between perceptual and neural tuning curves in the youngest sample of children tested to date and a significant relation between behavioral and neural sensitivity to numerosity. Specifically, children's right IPS neural Weber fractions predicted their behavioral Weber fractions. This provides some of the first evidence of a direct relation between neural and perceptual tuning and suggests that even though there may be differences between the perceptual Weber fraction and neural Weber fractions (Dehaene, 2007), they are related, albeit weakly, in early childhood.

Finally, we observed developmental continuity in the “distributed and overlapping” pattern of activation to number and brightness previously observed in the adult parietal cortex (Pinel et al., 2004). As in adults, children recruited bilateral regions of parietal and occipital cortices for number and brightness detection. Regions of the anterior IPS exhibited a bias for numerosity processing whereas the posterior IPS exhibited a bias for color and brightness processing. These findings from children parallel the distribution of neural responses from number and brightness in the study of adults (Pinel et al., 2004). We further showed that regions of the anterior IPS that exhibit neural tuning to numerosity do not exhibit tuning to brightness in children, and our positive control analysis showed that color-preferring regions of the posterior IPS are tuned to brightness but not numerosity. Thus, like adults, children's representations of number and brightness can be segregated in the cortex. These findings indicate that young children possess distinct mechanisms for representing numerosity.

Together, our results show that (1) approximate cardinal tuning to numerosity emerges in the right IPS by 3–4 years of age, (2) the laws of neural tuning to numerosity parallel the laws of perceptual tuning in the right IPS, and (3) the precision of children's neural tuning predicts their perceptual sensitivity to numerosity. These data support the conclusion that there is developmental continuity in the neural representation of cardinal numerosities over development, and that the right IPS underlies early development of numerical perception.

## Footnotes

This work was supported by National Science Foundation (NSF) Graduate Research Fellowship Program Grant DGE-1419118 to A.J.K., NSF Grant DRL1459625 to J.F.C., National Institutes of Health Grant R01 HD064636 to J.F.C., Alfred P. Sloan Foundation Fellowship #220020300 to J.F.C., and the James S. McDonnell Foundation. The authors thank Pat Weber and the Concepts, Actions, and Objects Lab, especially Santiago Alonso-Diaz, Elon Gaffin-Cahn, and Julia Yurkovic.

The authors declare no competing financial interests.

- Correspondence should be addressed to Jessica F. Cantlon, Department of Brain and Cognitive Sciences, University of Rochester, Rochester, NY 14627. jcantlon{at}rcbi.rochester.edu