Neural correlates of cognitive variability in childhood autism and relation to heterogeneity in decision-making dynamics.

Heterogeneity in cognitive and academic abilities is a prominent feature of autism spectrum disorder (ASD), yet little is known about its underlying causes. Here we combine functional brain imaging during numerical problem-solving with hierarchical drift-diffusion models of behavior and standardized measures of numerical abilities to investigate neural mechanisms underlying cognitive variability in children with ASD, and their IQ-matched Typically Developing (TD) peers. Although the two groups showed similar levels of brain activation, the relation to individual abilities differed markedly in ventral temporal-occipital, parietal and prefrontal regions important for numerical cognition: children with ASD showed a positive correlation between functional brain activation and numerical abilities, whereas TD children showed the opposite pattern. Despite similar accuracy and response times, decision thresholds were significantly higher in the ASD group, suggesting greater evidence required for problem-solving. Critically, the relationship between individual abilities and engagement of prefrontal control systems anchored in the anterior insula was differentially moderated by decision threshold in subgroups of children with ASD. Our findings uncover novel cognitive and neural sources of variability in academically-relevant cognitive skills in ASD and suggest that multilevel measures and latent decision-making dynamics can aid in characterization of cognitive variability and heterogeneity in neurodevelopmental disorders.

Introduction

Autism Spectrum Disorder (ASD) is a complex, debilitating, and heterogeneous set of neurodevelopmental conditions that affects 1 in 68 children (Christensen et al., 2016). As a spectrum disorder, phenotypic variability in core symptomology, cognitive abilities, behavioral profiles, and brain function are a key feature of ASD, presenting unique clinical and educational challenges (Ameis, 2017; American Psychological Association, 2013; Leekam et al., 2011; Lenroot and Yeung, 2013; Wing, 1981). The characterization of behavioral and neural sources of such heterogeneity remains largely unexplored. Critically, the majority of studies to date have focused on group-averaged differences between ASD and neurotypical groups (see for example Iuculano et al., 2014) and findings from previous studies have been largely inconsistent (Goldberg et al., 2005; Hill, 2004; Kana et al., 2007; Koshino et al., 2005; Minshew et al., 1994; Molesworth et al., 2005; Solomon et al., 2009; Urbain et al., 2015). Here, we take a novel approach that combines multilevel behavioral and brain measures with computational modeling to investigate sources of heterogeneity in cognitive skills important for academic learning during childhood, in children with ASD and their well-matched neurotypical peers. We focus on numerical problem-solving abilities as they are critical for academic achievement and quantitative reasoning in everyday life. Critically, even more so than reading, poor abilities in this domain can have negative consequences for employability, wages, socio-economic wellbeing, and life expectancy (Parsons and Bynner, 2005).

Heterogeneity is a hallmark of ASD: the severity and range of restricted interests and repetitive behaviors vary considerably in ASD, ranging from stereotyped motor movements and preference to sameness, to extreme rigidity and cognitive inflexibility (Leekam et al., 2011). Social communication abilities also vary in ASD (Wing and Gould, 1979), including attention and responses to social cues and levels of nonverbal engagement (Tager-Flusberg and Kasari, 2013). Furthermore, while some individuals with ASD show a marked delay in spoken language acquisition (Tager-Flusberg and Kasari, 2013), others show exceptional abilities in language learning (Tammet, 2006). In addition to these phenotypic domains, there is evidence for significant heterogeneity in intellectual, cognitive and academic abilities (Chen et al., 2018; Goin-Kochel et al., 2008; Iuculano et al., 2014; Mayes and Calhoun, 2003, 2008; Mottron et al., 2014; Oswald et al., 2016; Treffert, 2009). While there has been some progress in mapping heterogeneity in core phenotypic features of ASD, especially with intrinsic functional connectivity analysis (Abrams et al., 2013; Keown et al., 2013; Uddin et al., 2013a), there is a dearth of cognitive task-based investigations of neurocognitive variability in childhood autism.

Children with ASD demonstrate significant heterogeneity in academically-relevant cognitive skills, including numerical problem-solving abilities (Assouline et al., 2012; Baron-Cohen et al., 2007; Cash, 1999; Chiang and Lin, 2007; Kim and Cameron, 2016; Oswald et al., 2016; Wei et al., 2015). Numerical cognition provides a powerful domain for investigating neurocognitive profiles of heterogeneity in ASD as both superior abilities (Baron-Cohen et al., 2007; Chiang and Lin, 2007; Iuculano et al., 2014; Minshew et al., 1994) and underachievement (Chiang and Lin, 2007; Estes et al., 2011; Mayes and Calhoun, 2003) have been reported in the extant literature. A recent quantitative analysis of heterogeneity, using gap statistics in a large sample of 114 children with ASD and 96 matched controls, revealed a unique profile of heterogeneity in ASD, consisting of a low-achieving subgroup with poor math skills compared to reading, and a high-achieving subgroup who showed superior math skills compared to reading (Chen et al., 2018). However, these studies do not provide a comprehensive characterization of neurobehavioral mechanisms that lie at the core of heterogeneity of skills in ASD. Multilevel characterization of heterogeneity in numerical problem-solving skills using functional brain imaging and cognitive-behavioral modeling has the potential to improve our understanding of individual differences in cognition and inform early remediation of academically-relevant cognitive deficits in affected children.

Here we investigate sources of heterogeneity in numerical problem-solving abilities in ASD by integrating multiple levels of analyses (Fig. 1). Our general approach was to probe whether sources of heterogeneity in cognitive abilities in ASD are similar to the neurotypical population, or whether they manifest in a way that is unique to the disorder. Specifically, we determined whether differences in activation of task-related brain systems and decision-making processes underlie heterogeneity of problem-solving abilities in ASD. An important aspect of our investigation relates to the role of latent decision-making processes supporting numerical problem-solving, and their neural correlates. Analysis of latent variables underlying decision-making using drift diffusion modeling has the potential to provide a more comprehensive framework for understanding cognitive processes beyond observed behavioral metrics such as accuracy and reaction time (Ratcliff and McKoon, 2008; Ratcliff and Smith, 2004). These models assess how a decision-choice depends on the accumulation of evidence for various response alternatives (Froehlich et al., 2016; Oganian et al., 2016; Ratcliff and McKoon, 2008; Ratcliff and Smith, 2004; Spaniol et al., 2006; Thapar et al., 2003; Vandekerckhove et al., 2010). A previous study of visual perception demonstrated that individuals with ASD show differences in latent decision-making processes despite similarities in overall accuracy and reaction time to neurotypical individuals (Pirrone et al., 2017).

Fig. 1

Schematic illustration of study design and analyses. We combined clinical, cognitive, behavioral, and neural levels of investigation to assess their contributions to individual differences in numerical problem-solving in children with ASD. Clinical measures (ADOS and ADI-R scores) were used to characterize the ASD group. Academically-relevant cognitive problem-solving skills were assessed using the Numerical Operations test of the WIAT-II. This test was administered to both the ASD and TD groups. Behavioral measures used a two-alternative forced-choice task of numerical problem-solving in which every child – in the ASD and TD groups – was asked to assess the validity of single-digit addition equations (e.g. 5 + 4 = 9). Neural. Event-related functional Magnetic Resonance Imaging (fMRI) task.Hierarchical drift diffusion modeling (HDDM) of reaction time data during the forced-choice fMRI task was used to assess decision-making dynamics in both groups. All analyses included average group-level comparisons between the ASD and TD groups at each level of investigation. Individual variability was assessed using general linear models, as a function of WIAT-II Numerical Operations scores. Schematic graphs illustrate possible relations between Numerical Operations and clinical, behavioral, and neural features. Moderation analyses assessed the influence of behavioral decision-making parameters to the relation between functional brain activation and cognitive abilities in ASD and TD children. Abbreviations: ASD = Autism Spectrum Disorder; TD = Typically Developing; ADOS = Autism Diagnostic Observation Schedule; ADI-R = Autism Diagnostic Interview - Revised; RIRB = Restricted Interests and Repetitive Behaviors; WIAT-II = Wechsler Individual Achievement Test - Second Edition; RTs = Reaction Times; HDDM = Hierarchical Drift Diffusion Model.

At the neural level, we investigated the relation between cognitive abilities and functional activation during a numerical problem-solving task. In TD children, numerical problem-solving is shown to be supported by a widely distributed network of brain regions, including those involved in (i) visual perception and symbol recognition, anchored in the fusiform gyrus (FG) within the ventral temporal-occipital cortex (VTOC); (ii) semantic processing of quantity, anchored in the intraparietal sulcus (IPS) region of the posterior parietal cortex (PPC); (iii) attention and cognitive control processes supported by prefrontal cortical regions including the anterior insula and the ventrolateral and dorsolateral prefrontal cortices (VLPFC, DLPFC); and (iv) mnemonic processing supported by the medial temporal lobe (MTL) (Iuculano and Menon, 2018; Menon, 2014). We previously demonstrated that children with ASD show different patterns of multivariate activity across these brain systems during numerical problem-solving (Iuculano et al., 2014). Moreover, there is evidence of dysfunction of prefrontal control systems, and the anterior insular cortex in particular, in ASD (Uddin et al., 2013a).

We tested a group of children (7−12 years old) with ASD and an age-, sex-, IQ-, and in-scanner motion parameters-matched Typically Developing (TD) control group. We hypothesized that latent decision-making variables would differ in children with ASD from TD children during problem-solving. We further hypothesized that individual differences in problem-solving abilities in children with ASD would be associated with aberrant profiles of brain response in VTOC, PPC, MTL and PFC regions important for numerical cognition (Iuculano et al., 2014). Finally, we tested the hypothesis that latent decision-making variables would contribute to individual differences in numerical problem-solving abilities and brain responses in ASD.

Methods

We combined clinical, cognitive, behavioral, and neural levels of investigation to assess their contributions to individual differences in numerical problem-solving in children with ASD (Fig. 1).

Participants

For all participants, informed written consent was obtained from the legal guardian of the child and all study protocols were approved by the Stanford University Review Board. Twenty-four 7−12 year-old children with Autism Spectrum Disorders (ASD) were recruited from the San Francisco Bay area, along with a large group of Typically Developing (TD) children without known genetic, psychiatric, or neurological disorders (N = 138). The TD group was recruited as part of a neurodevelopmental study of mathematical skill development (Chen et al., 2018; Qin et al., 2014). Eight children with ASD were excluded because they did not meet inclusion criteria for: (i) Full-Scale IQ (FSIQ > 70) (N = 1); (ii) quality of un-normalized functional images (N = 2); (iii) in-scanner motion parameters (total frames interpolated < 15 %) (N = 3); (iv) comorbidity with other neuropsychiatric disorders (Tourette syndrome) (N = 1); and (v) scanning-session compliance (N = 1). The final sample consisted of sixteen children with ASD (fifteen boys; mean age = 9.46; SD = 1.80), who were all verbal and within a normal IQ-range (FSIQ range: 93–142). Additional data on the ASD sample is included in (Table S1). FSIQ was based on the Wechsler Abbreviated Scale of Intelligence (WASI) (Wechsler, 1999). Diagnosis of ASD was based on scores from the Autism Diagnostic Interview-Revised (ADI-R) (Le Couteur et al., 1989; Lord et al., 1994; Rutter et al., 2003) and/or the Autism Diagnostic Observation Schedule (ADOS) (Gotham et al., 2007; Lord et al., 2000). A group of TD controls who were matched on age-, sex-, IQ-scores (Table 1) and in-scanner motion parameters (Supplementary Information – Methods – fMRI Preprocessing) to the ASD group, was selected from the larger sample of 7−12 year old TD children using a matching algorithm (Uddin et al., 2013a). Specifically, participants were matched on FSIQ scores (WASI-scores), age, sex and 6 fMRI motion parameters (x, y, z, roll, pitch, yaw) using the genetic algorithm (Uddin et al., 2013a). The final matched TD group consisted of sixteen children (fifteen boys; mean age = 9.83; SD = 1.75) (Table 1) within a normal IQ range (FSIQ range: 85–137) .

Table 1

Demographic, cognitive and diagnostic measures. Demographic, mean IQ scores, and Working Memory (WM) scores are shown for the ASD and TD groups. Mean ADI-R and ADOS scores are shown for the ASD group only. Standard deviations (SD) are shown in parentheses.

Measure	ASD (N = 16)	TD (N = 16)	p-value
Male to Female ratio	15:1	15:1
Age (years)	9.46 (±1.80)	9.83 (±1.75)	.56
WASI scale
Verbal IQ	117.31 (±16.39)	121.88 (±11.19)	.37
Performance IQ	119.00 (±19.83)	114.37 (±15.97)	.47
Full IQ	120.25 (±15.25)	120.31 (±11.72)	.99
WIAT-II scale
Numerical Operations	116.88 (±21.77)	106.37 (±17.22)	.14
Math Reasoning	123.25 (±17.51)	115.31 (±14.37)	.17
Word Reading	117.75 (±8.05)	113.44 (±10.65)	.21
Reading Comprehension	111.75 (±12.07)	111.75 (±11.83)	1.00
WMTB - C
Digit Recall	114.88 (±19.38)	109.05 (±16.45)	.42
Block Recall	93.56 (±21.50)	93.60 (±13.20)	.99
Count Recall	96.33 (±20.70)	94.75 (±15.40)	.81
Backwards Digit Recall	98.69 (±21.34)	92.28 (±17.85)	.54
ADI-R
Social (a)	18.87 (±5.55)
Verbal (b)	14.94 (±6.01)
Repetitive Behavior (c)	5.81 (±3.31)
Development (d)	2.62 (±1.59)
ADOS
Social/Affect	10.60 (±3.02)
Restricted and Repetitive Behavior	2.33 (±1.54)
Severity Scores	7.47 (±1.50)
Total	12.93 (±3.41)

Abbreviations: WASI = Wechsler Abbreviated Scale of Intelligence; WIAT – II = Wechsler Individual Achievement Test – Second Edition; WMTB – C = Working Memory Test Battery for Children; ADI-R = Autism Diagnostic Interview – Revised (diagnostic scores) (Rutter et al., 2003); ADOS = Autism Diagnostic Observation Schedule – new algorithm (Gotham et al., 2007); df = (1,30) for all statistics.

Children in the ASD group did not differ from their TD peers on Verbal and Performance IQ measures (all p > .37), or the visuospatial, phonological, and central executive components of working memory (all p > .42). They also did not differ on math reasoning, word reading and reading comprehension (Table 1) (all p > .17).

Neuropsychological assessments

Each participant was administered additional neuropsychological assessments including the Working Memory Test Battery for Children WMTB-C (Pickering and Gathercole, 2001), and the Wechsler Individual Achievement Test - Second Edition (WIAT-II) (Wechsler, 2001) for mathematical abilities. Standardized scores from the Numerical Operations subtest of the WIAT-II was used as our primary metric of numerical abilities.

Event-related numerical problem-solving functional Magnetic Resonance Imaging (fMRI) task

The event-related numerical problem-solving fMRI task consisted of two arithmetic verification conditions: (i) addition problems, and (ii) ‘plus 1′ problems, which served as a high-level control. Previous research has shown that simple ‘N+1′ addition is solved by incremental counting (Campbell and Metcalfe, 2007) with higher accuracy and faster reaction times relative to more complex addition problems (Cho et al., 2011). Moreover, because stimuli in ‘plus 1′ control problems have the same format as the addition problems, they provide a high-level control for sensory and number processing, as well as decision-making and response selection.

Participants were presented with an equation involving two addends and asked to indicate, via a button-press, whether the answer shown was valid (e.g. 3 + 4 = 7) or invalid (e.g. 3 + 4 = 8). One operand ranged from 2 to 9 and the other ranged from 2 to 5 (tie problems, such as 5 + 5 = 10 were excluded). Control problems were identical in format to the addition problems except that one of the addends was always ‘1′. There were a total of 26 trials per condition (addition, control). Each equation was presented for five seconds followed by a jittered fixation from 2.5 to 3.5 seconds and the total length of the experimental run was 6 min 30 s (Supplementary Information – Methods – fMRI stimuli and task). Performance (accuracy and reaction times) during this task was analyzed using a two-way analysis of variance with within-subject factor problem-type (‘addition’, ‘control’) and between-subject factor group (ASD, TD). Based on recent recommendations to use Bayesian statistics for evaluating evidence against the null hypothesis (i.e. no group differences) (Wagenmakers et al., 2018a, b), we calculated Bayes Factors (BFs) (Kass and Raftery, 1995; Wetzels and Wagenmakers, 2012) on significant results (JASP Team, 2018) (“Bayes Factor” R Package, R-Version 3.4.1; 2017). We also computed Cohen’s f2 effect sizes (Cohen, 1988) (R-Version 3.4.1; 2017) to further assess robustness of our findings.

Hierarchical Drift-Diffusion model of decision-making dynamics

A Hierarchical Drift Diffusion Model (HDDM) was applied to the distribution of reaction times from the numerical problem-solving fMRI task using the python toolbox ‘HDDM’ (Wiecki et al., 2013). Drift diffusion models account for the decision process as a continuous sampling of information that accumulates over time until a decision threshold is reached, and have been widely used to analyze latent decision-making variables during forced-choice tasks (Froehlich et al., 2016; Oganian et al., 2016; Ratcliff et al., 2004; Ratcliff and McKoon, 2008; Thapar et al., 2003), also in ASD (Pirrone et al., 2017; South et al., 2014). Output model parameters of interest include: decision threshold ‘a’ characterizing the amount of information that needs to be accumulated for the decision; drift-rate ‘v’ representing the speed of evidence-accumulation over time; and non-decision time ‘t’ representing the aspects of response time not included in active deliberation (e.g. stimulus encoding, motor response). The present study did not model response-bias ‘z’ as we had no reason to believe that an a-priori bias existed here.

HDDM is better suited to handle data with smaller number of trials, such as those from fMRI tasks. A Bayesian approach is used to estimate model parameters by assigning prior probability distributions to each parameter. Markov Chain Monte-Carlo (MCMC) sampling is then used to approximate the posterior distribution for each parameter at both the individual- and group- levels. We initialized HDDM with five chains (5000 samples, first 200 discarded as burn-in each, for a total of 24,000 samples and 1000 burn-in) and estimated the parameters ‘a’, ‘v’, and ‘t’ for each task-condition (addition and control) in each group (ASD, TD). Model convergence was determined using the Gelman-Rubin statistic (R-hat) (Gelman and Rubin, 1992) (Supplementary Information – Methods – HDDM Model Convergence). To compare parameter estimates between groups on the addition task-condition, we used Bayesian hypothesis-testing and calculated the proportion of overlap between the group-level posterior distributions of the parameters ‘a’, ‘v’, and ‘t’. Finally, we examined the relation between each latent decision variable and Numerical Operations scores.

Here also, based on recent recommendations to use Bayesian statistics for evaluating evidence against the null hypothesis (Wagenmakers et al., 2018a, b) (i.e. no group differences), we calculated Bayes Factors (BFs) (Kass and Raftery, 1995; Wetzels and Wagenmakers, 2012) on significant results (JASP Team, 2018) (“Bayes Factor” R Package, R-Version 3.4.1; 2017). We also computed Cohen’s d effect sizes (Cohen, 1988) (R-Version 3.4.1; 2017), and power to further assess robustness of our findings.

Functional MRI data analysis

Functional MRI preprocessing

fMRI data were preprocessed using SPM8 (http://www.fil.ion.ucl.ac.uk/spm/) (Supplementary Information – Methods – fMRI Preprocessing).

Brain activation during numerical problem-solving

Numerical problem-solving task-related brain activation was identified using the General Linear Model (GLM) implemented in SPM8. Voxel-wise contrast and t-statistic images were generated for each participant by contrasting correct addition trials versus correct control trials (Supplementary Information – Methods – fMRI - General Linear Model). Differences in brain activation between the ASD and TD groups were examined using a t-test on these contrast images (correct addition trials vs. correct control trials). Significant clusters of activation were identified using a height threshold of p < .001, with family-wise error (FWE) correction for multiple spatial comparisons at the cluster level (p < .05, spatial extent: 389 voxels), as implemented in SPM8.

Relation between brain activation and numerical problem-solving abilities

To assess brain activation as a function of numerical abilities and group (ASD or TD), we conducted a whole brain one-sample t-test on the contrast images (correct ‘addition’ versus correct ‘control’ trials) with Numerical Operations scores of the WIAT-II and Group (ASD and TD) as covariates of interest in the model. Significant clusters of activation were identified using the same activation threshold (p < .001; FWE p < .05, spatial extent: 389 voxels).

To assess the robustness of our findings we computed Cohen’s f2 (Cohen, 1988), Cohen’s d (Cohen, 1988) and Bayes Factor (BF) values (Kass and Raftery, 1995; Wetzels and Wagenmakers, 2012) (R-Version 3.4.1; 2017) for reported effects in each brain region in which task-related functional brain activity showed a significant interaction between Group and Numerical Operations scores.

Influence of decision-making dynamics on the relationship between brain activation and numerical problem-solving abilities

Moderation analysis as implemented in the R-package (Version 3.4.1; 2017) was used to assess how decision-making dynamics influence the relationship between brain activation and numerical problem-solving abilities. We first demeaned individual Numerical Operations scores and decision parameters using the grand mean of the combined ASD and TD samples. This analysis was conducted in the combined group to avoid circularity and under the assumption that some proportion of children with ASD may have profiles similar to TD children, while others may differ significantly. This allowed us to more fully characterize sources of heterogeneity in the ASD group. Numerical Operations scores were entered as the dependent variable; brain activation for ‘addition’ vs. ‘control’ problems was entered as the independent variable, and HDDM scores that distinguished the ASD and TD groups were entered as the moderator variable. One participant was excluded based on outlier detection procedures (i.e. > 3SD based on Cook’s Distance). A total of 31 subjects were therefore included in the moderation analysis. Robustness of moderation effects was assessed using Cohen’s f2 (Cohen, 1988) (R-Version 3.4.1; 2017), and evidence against the null hypothesis (i.e. no moderation) using BF (Kass and Raftery, 1995; Wetzels and Wagenmakers, 2012) implemented in the “Bayes Factor” R Package (R-Version 3.4.1; 2017).

Results

Numerical problem-solving abilities did not differ between the ASD and TD groups

ASD and TD groups did not significantly differ on any standardized measures (Table 1), including our main measure of numerical ability, Numerical Operations from the WIAT-II (p = .14). Range of problem-solving abilities in the ASD and TD groups showed a similar profile: 82–158 in the ASD group, and 83–140 in the TD group. Critically, variance on the Numerical Operations scores did not differ between groups: [F(15,15) = 1.5978, p = .37].

Individual differences in numerical abilities are not related to standardized measures of IQ, Working Memory or symptom severity

IQ, visuospatial, phonological, and central executive components of working memory were not correlated with Numerical Operations scores in either the ASD or TD groups (all p > .1) (Table S2). None of the clinical measures from the ADI-R and ADOS were significantly correlated with Numerical Operations scores in children with ASD (all p > .26) (Table S2).

Performance during fMRI task did not differ between the ASD and TD groups

Accuracy

There was a main effect of problem-type (‘addition’, ‘control’) [F(1,30) = 20.55, p < .0001; Cohen’s f2 =0.676; 95 % Confidence Limits (C.L.) = .156 ∼ 1.434; BF = 3.09 × 102], with reduced accuracy on addition problems relative to control problems for both ASD and TD groups, consistent with previous studies showing that addition problems are significantly more difficult than control problems (Campbell and Metcalfe, 2007; Cho et al., 2011). There was no main effect of group [F(1,30) = 1.691, p = .203; Cohen’s f2 =0.056; 95 % C.L. = 0 ∼ 0.336], and no group by problem-type interaction [F(1,30) = 0.407, p = .528; Cohen’s f2 =0.008; 95 % C.L. = 0 ∼ 0.209] (Fig. 2A).

Fig. 2

ASD and TD groups differ on latent but not overt behavioral measures. A. No significant differences were found on accuracy or reaction times (RTs) between ASD and TD groups during the numerical problem-solving task. B.Left: Schematic representation of the Hierarchical Drift Diffusion Model (HDDM) parameters for RTs: non-decision time ‘t’, drift rate ‘v’, and threshold ‘a’. Right: HDMM parameters for current ASD and TD groups. C & D. A significant difference was found between ASD and TD groups for the decision threshold parameter ‘a’ (p = .00013). No difference was found for non-decision time (t) and drift rate (v) parameters. (C) posterior probability density, and (D) mean parameter estimates of the three HDDM parameters.

Reaction Times

There was a main effect of problem-type [F(1,30) = 80.609, p < .0001; Cohen’s f2 =2.682; 95 % C.L. = 1.107 ∼ 4.453; BF = 1.52 × 107]; all participants responded faster to control than addition problems, consistent with previous studies showing that addition problems are significantly more difficult than control problems (Campbell and Metcalfe, 2007; Cho et al., 2011). There was no main effect of group [F(1,30) = .043, p = .836; Cohen’s f2 =0.001; 95 % C.L. = 0 ∼ 0.119], nor an interaction between group and problem-type [F(1,30) = .057, p = .812; Cohen’s f2 =0.0005; 95 % C.L. = 0 ∼ 0.129] (Fig. 2A).

Drift-diffusion modeling of reaction times revealed differences in decision-making dynamics between the ASD and TD groups

Convergence statistics (R-hat) for all drift diffusion parameters were between 0.98 and 1.02 indicating good model fit (Gelman and Rubin, 1992). We found that the threshold parameter ‘a’, was significantly higher in children with ASD compared to TD children (t30) = -4.596, p = .00013; Cohen’s d = 1.625; 95 % C.L. = 0.810 ∼ 2.420; BF = 2.77 × 102; power = 0.9935177) (Fig. 2B–D). These values reflect strong evidence for the reported findings (Kass and Raftery, 1995). There were no group differences in non-decision time ‘t’ (t(30) = 1.959, p = .060) or drift rate ‘v’ (t(30)= -1.323, p = .1966) (Fig. 2B–D). None of the three latent decision variables were correlated with individual differences in Numerical Operations scores in either ASD (all p > .24) or TD groups (all p > .08) (Table S3).

Individual differences in numerical problem-solving abilities are supported by dissociable profiles of functional brain activation in the ASD and TD groups

No brain regions showed significant differences in functional activation between the ASD and TD groups during the numerical problem-solving task (i.e. did not meet the threshold of p < .001 height, FWE p < .05 cluster extent). However, analysis of interaction between Group (ASD, TD) and Numerical Operations scores revealed significant and opposite patterns of brain activation in the ASD and TD groups in multiple cortical and subcortical regions (Fig. 3A). Specifically, in the ASD group, successful problem-solving was associated with increased activation of right anterior insular cortex, and bilateral dorsolateral and dorsomedial prefrontal cortices (DLPFC, DMFC) in the prefrontal cortex; left lateral occipital cortex (LOC) and bilateral fusiform gyrus (FG) in the ventral temporal-occipital cortex (VTOC); and bilateral intraparietal sulci (IPS, hIP1 and hIP3) in the posterior parietal cortex (Fig. 3B, Table 2). In contrast, the TD group showed exactly the opposite pattern, with a negative correlation between Numerical Operations score and brain activation in prefrontal, VTOC, and posterior parietal regions (Fig. 3B, Table 2). Differences between groups as a function of Numerical Operations scores were also evident in subcortical regions, including the bilateral thalamus (Fig. 3B, Table 2). BFs, quantifying the strength of the evidence for significant interaction between Group and Numerical Operations revealed values ranging from 2.46 to 3.28 × 102, with the highest values in the right anterior insular cortex (Table 2), reflecting strong evidence (Kass and Raftery, 1995) for differential engagement of this region in the two groups as a function of numerical abilities. Analysis of simple slopes within each group further clarified these results (Tables S4 and S5).

Fig. 3

Relation between task-related brain activation and numerical problem-solving abilities differs between ASD and TD groups. A. Surface rendering of brain regions showing a significant Group x Numerical Operations interaction in multiple prefrontal, parietal, ventral temporal-occipital, and subcortical regions; B. Children with ASD showed increased brain activation with numerical abilities (red line) while TD children (blue line) exhibited decreased activation with numerical abilities. The Numerical Operations subtest of the standardized WIAT-II assessments was used to assess numerical abilities in each child. Abbreviations: L = Left; R = Right; DLPFC = Dorsolateral Prefrontal Cortex; DMPFC = Dorsomedial Prefrontal Cortex; LOC = Lateral Occipital Cortex; FG = Fusiform Gyrus; IPL = Inferior Parietal Lobule; hIP1 = horizontal segment of the Intraparietal sulcus, subdivision 1; IPS = Intraparietal Sulcus; hIP3 = horizontal segment of the Intraparietal sulcus, subdivision 3.

Table 2

Brain areas that showed a significant interaction between problem-solving ability – assessed using . Effect sizes, Confidence Limits (C.L) and Bayes Factors (BFs10) are reported for each Region Of Interest (ROI). Note: BFs10 > 5 (in bold) provide moderate to strong support against the null hypothesis (Kass and Raftery, 1995).

Brain Region	MNI coordinates			Max	Cluster size	Effect size	C.L. of Cohen’s f2	BF10
	x	y	z	t		Cohen’s f2
Prefrontal Cortex
R Anterior Insular Cortex	36	18	−6	5.32	1743†	0.9151	0.173 ∼ 1.637	3.28 × 102
R IFG, pars triangularis/VLPFC	44	18	8	4.55		0.7297	0.104 ∼ 1.363	7.8 × 10
R SFG, posterior segment/DMPFC	22	−6	56	5.27	7230††	1.0715	0.234 ∼ 1.865	2.71 × 102
L Precentral Gyrus/Premotor Cortex	−20	−14	68	4.85		0.7156	0.099 ∼ 1.342	7.3 × 10
L SFG, posterior segment/DMPFC	−18	−6	64	4.80		0.7380	0.107 ∼ 1.375	8.4 × 10
R DLPFC	40	6	46	4.66		0.4827	0.025 ∼ 0.987	11.55
L DLPFC	−38	14	54	4.45		0.2825	NA ∼ 0.665	2.46
L SMA	−6	6	58	4.42		0.5281	0.038 ∼ 1.058	18.77
R SMA	6	6	58	4.33		0.7964	0.128 ∼ 1.462	1.35 × 102
L Paracingulate Cortex	−8	38	24	4.55	1350	0.6446	0.075 ∼ 1.235	38.14
R Paracingulate Cortex	4	46	12	4.32		0.5779	0.053 ∼ 1.134	18.79
R Anterior Cingulate Cortex	12	38	18	3.86		0.5399	0.041 ∼ 1.076	14.96
L Anterior Cingulate Cortex	−8	18	32	3.47		0.3427	NA ∼ 0.765	3.54
Parietal Cortex
L SMG	−56	−38	48	5.11	††	0.568	0.050 ∼ 1.119	24.51
L IPL/hIP1	−34	−64	44	4.86	††	0.498	0.029 ∼ 1.011	15.74
L Postcentral Gyrus	−56	−28	48	4.75	††	0.602	0.061 ∼ 1.171	34.83
L Precuneous	−12	−42	66	4.55	††	0.348	NA ∼ 0.774	4.49
L IPS	−24	−64	46	4.39	††	0.412	0.006 ∼ 0.876	7.65
R IPL/PGa	44	−60	50	4.75	471	0.494	0.029 ∼ 1.006	14.09
R IPS	30	−60	58	3.70		0.348	NA ∼ 0.773	4.44
Lateral/Medial Temporal Cortex
R STG/Planum Temporale	58	−22	10	4.07	†	0.407	0.004 ∼ 0.868	6.97
R Amygdala	30	−2	−16	3.98	†	0.495	0.028 ∼ 1.007	9.36
R STG/Heschl’s Gyrus	44	−20	12	3.64	†	0.297	NA ∼ 0.689	2.66
Ventral Temporal-Occipital Cortex
L LOC inferior division/FG	−48	−66	−14	5.10	4882*	0.558	0.047 ∼ 1.104	23.85
R Intracalcarine Cortex	12	−84	6	4.93		0.796	0.128 ∼ 1.462	10.9 × 10
L LOC inferior division	−32	−80	4	4.91		0.581	0.054 ∼ 1.138	24.03
L Occipital Fusiform Gyrus	−30	−74	−12	4.56		0.631	0.070 ∼ 1.214	29.42
R Lingual Gyrus	24	−72	0	4.33		0.532	0.040 ∼ 1.064	17.26
L Intracalcarine Cortex	−20	−64	10	4.30		0.571	0.051 ∼ 1.123	16.90
L LOC superior division	−30	−76	30	4.29		0.461	0.019 ∼ 0.954	11.51
L Temporal Occipital Fusiform Cortex	−54	−48	−22	4.22		0.460	0.018 ∼ 0.953	10.21
R Occipital Fusiform Gyrus	16	−74	−8	4.19		0.428	0.010 ∼ 0.902	7.93
Subcortical Structures
L Thalamus	−12	−20	−2	4.77	1260	0.490	0.027 ∼ 0.999	14.71
R Cerebellum (lobule IV)	24	−40	−24	4.26		0.510	0.032 ∼ 1.029	13.98
R Cerebellum (lobule VI)	34	−68	−22	4.88	*	0.598	0.060 ∼ 1.165	28.04
L Cerebellum (lobule VI/Crus I)	−4	−72	−20	4.55	*	0.567	0.050 ∼ 1.116	23.84
L Cerebellum (Vermis)	−2	−56	−14	4.10	*	0.379	NA ∼ 0.824	5.31
R Putamen	30	2	2	4.53	†	0.561	0.048 ∼ 1.107	25.30
R Thalamus	14	−18	4	4.40	389	0.457	0.018 ∼ 0.947	11.11
R Midbrain	10	−22	−10	3.73		0.554	0.046 ∼ 1.097	14.04
R Fornix	22	−30	10	3.73		0.321	NA ∼ 0.730	3.45

Abbreviations: BF10 = Bayes Factor in favor of H1 over H0; IFG = Inferior Frontal Gyrus; VLPFC = Ventrolateral Prefrontal Cortex; SFG = Superior Frontal Gyrus, DMPFC = Dorsomedial Prefrontal Cortex; DLPFC = Dorsolateral Prefrontal Cortex; SMA = Supplementary Motor Area; SMG = Supramarginal Gyrus; IPL = Inferior Parietal Lobule; hIP1= horizontal segment of the Intraparietal sulcus, subdivision 1; IPS = Intraparietal Sulcus; STG = Superior Temporal Gyrus; LOC = Lateral Occipital Cortex; FG = Fusiform Gyrus; †,††,* These are anatomically distinct sub-peaks of the same activated cluster. Sub-peaks are italicized. Note: p < .001 height, p < .05 extent. Moderate to strong BFs are bolded.

Distinct and heterogenous influence of decision-making dynamics on the relationship between brain activation and numerical problem-solving abilities in ASD

Next, we used decision-making parameters from the HDDM model to further investigate mechanisms underlying the relation between brain activation and problem-solving abilities, and to elucidate additional sources of heterogeneity in the ASD group. We performed a moderation analysis to investigate how the decision threshold ‘a’, the decision-making parameter that differed between ASD and TD groups (Fig. 2B–D), influenced the relationship between Numerical Operations scores and brain activation during the numerical problem-solving task. This analysis focused on the right anterior insula for two reasons: (i) it has been previously reported to be aberrant in autism (Uddin et al., 2013b, 2015), and (ii) it has been demonstrated to have a critical role in high-level cognitive control and attentional processes (Cai et al., 2014; Menon and Uddin, 2010), which are critical for problem-solving (Supekar and Menon, 2012).

In the combined-group sample, we found that the decision threshold ‘a’ significantly moderated the relationship between Numerical Operations and anterior insula activation (B = 11.337, t = 2.95, p = .000645; Cohen’s f2 = 0.486; 95 % C.L. = 0.021 ∼ .999; BF = 8.61; power = 0.871) (Fig. 4A). To further explore heterogeneity related to the decision threshold, we next split the entire sample of participants into two groups based on the decision-threshold parameter. Specifically, we divided the combined group based on the median split of the decision threshold parameter and performed post-hoc analyses on the Low and High decision threshold groups. In both the Low and High decision threshold groups, brain activation in the right anterior insular cortex was associated with Numerical Operations scores, albeit in opposite directions (Low: B = -7.848, t = -2.524, p = .0243; Cohen’s f2 = 0.455; 95 % C.L. = 0.001 ∼ 1.354; BF = 2.816) and High (B = 8.056, t = 2.456, p = 0.0289, Cohen’s f2 = 0.464; 95 % C.L. = 0.000 ∼ 1.404; BF = 2.523) (Fig. 4B). Crucially, the high decision threshold group was over-represented by children with ASD (filled dots) while the low decision threshold group was over-represented by TD children (χ² = 7.2419; p = .0071). This result reveals the heterogenous influence of decision-making dynamics on the relationship between brain activation and numerical problem-solving abilities in children with ASD.

Fig. 4

Decision threshold moderates the relation between task-related brain activation and problem-solving. A. The decision threshold parameter (‘a’) from the hierarchical drift diffusion model moderated the relation between right insula activation during problem-solving and individual differences in numerical abilities assessed using the Numerical Operations subtest of the WIAT-II. ** p < .01. B. Significant correlation between task-related activation in the right anterior insula (axial slice) differed between the high (red; p = .0289) and low (blue; p = .0243) decision-threshold groups. The high decision threshold group was over-represented by children with ASD (filled dots) while the low decision threshold group was over-represented by TD children (p < .01). Abbreviations: L = Left; R = Right.

Discussion

Heterogeneity in clinical, cognitive, and behavioral profiles is a prominent feature of ASD, yet our understanding of the nature of individual variation, especially in the context of cognitive problem-solving and academic achievement, is severely limited. Here, using a well-characterized group of children ages 7–12 with ASD, and an age-, sex-, IQ-, and in-scanner motion parameters-matched group of TD children, we probed behavioral and neural sources of heterogeneity in ASD using multiple levels of investigation (Fig. 1). Compared to TD children, children with ASD exhibited a higher decision threshold during numerical problem-solving, despite a lack of differences in accuracy and reaction times, and overall numerical abilities (Fig. 2; Table1). Compared to a group of well-matched TD children, children with ASD also showed dramatically different patterns of individual differences in brain activation as a function of numerical problem-solving abilities (Fig. 3). Critically, effect sizes and Bayes Factor values, particularly in the anterior insular cortex, reflected strong evidence for these findings, despite a smaller sample. Lastly, the relationship between numerical abilities and activation of the right anterior insula, a brain region important for salience detection and cognitive control (Cai et al., 2014; Menon and Uddin, 2010), and a dysfunctional hub in ASD (Uddin et al., 2013a, b, 2015) was moderated by the decision threshold during problem-solving with heterogeneous profiles in children with ASD (Fig. 4). Our findings highlight unique neural and behavioral sources of heterogeneity in ASD and emphasize how latent behavioral processes can uncover different profiles of individual differences. Our study highlights quantitative analysis of heterogeneity as a critical step in autism research and provides a novel framework for examining individual differences in ASD that can be extended to other cognitive domains and clinical populations.

Children with ASD show aberrant decision-making dynamics during numerical problem-solving

Our first goal was to investigate differences in numerical problem-solving abilities between the ASD and TD groups. The two groups did not differ on standardized measures of numerical ability assessed using the WIAT-II Numerical Operations subtest (Table 1). We also did not find any group differences in performance on the forced-choice arithmetic verification task that was completed during fMRI scanning (Fig. 2A). However, children with ASD demonstrated significant differences from TD children in decision-making processes during problem-solving, with strong effect sizes (Cohen’s d = 1.625) and high Bayes factors (BF = 2.77 × 102) (Fig. 2C,D).

Decision-making processes were assessed using HDDM, which provided quantitative information about three latent variables: ‘t’, the amount of time needed to execute a response, ‘v’, the speed of evidence accumulation in order to reach a decision, and ‘a’, the amount of information required to make a decision (Fig. 2B). Although reaction times did not differ between the groups (Fig. 2A), the ASD group demonstrated a significantly higher decision threshold than the TD group (Fig. 2C,D). Higher levels of decision threshold likely reflect the need to accumulate greater amounts of information in the face of competing alternatives and internal mental processes (Voss et al., 2004). Although the increased decision threshold in children with ASD did not affect their ability to accurately solve arithmetic problems, it might impair their decision-making process in other cognitive contexts when finding an efficient speed-accuracy trade-off is critical. More generally, our findings point to the effectiveness of decision-making models in disentangling important components of information processing during cognitive problem-solving in children with ASD. As we demonstrate below, individual differences in decision-making are an important source of variability in brain-behavior relations in ASD.

Individual differences in numerical problem-solving abilities are supported by unique patterns of functional brain activation in children with ASD

Individual differences in numerical abilities were supported by a distributed network of brain regions in both ASD and TD groups (Fig 3A, Table 2). These included the left LOC and bilateral FG in the VTOC that support high-level perceptual representations of visual number-form (Ansari, 2008; Cantlon et al., 2011); bilateral posterior IPS subdivisions hIP1, hIP3 involved in representations and manipulations of quantity (Ansari, 2008; Arsalidou and Taylor, 2011; Dehaene et al., 2003), together with fronto-parietal working memory systems anchored in the IPL/SMG and bilateral DLPFC, and right anterior insula regions important for attention and cognitive control (Chang et al., 2016; Iuculano & Menon, 2018; Menon, 2014; Supekar and Menon, 2012).

Critically, the ASD and TD groups showed completely opposite patterns-relations (Fig. 3B), with high BF10 values (Table 2), reflecting moderate to strong evidence for our findings (Kass and Raftery, 1995). In the ASD group, brain activity in the aforementioned brain regions increased with numerical ability (Fig. 3B), whereas in the TD group, brain activity decreased with numerical ability (Fig. 3B). The negative association between brain activation during numerical problem- solving and numerical abilities in TD children is consistent with previous reports that children with lower math abilities show hyper-activation of posterior parietal, prefrontal, and ventral temporal-occipital areas during numerical problem-solving (Iuculano et al., 2015; Rosenberg-Lee et al., 2015). The positive association between brain activation during numerical problem-solving and numerical abilities in children with ASD suggest that, while children with ASD engage similar brain regions during numerical problem-solving as TD children, the nature of variability underlying individual skills differs markedly in these children. This unique pattern of neural processing points to fundamental differences in how children with ASD process and manipulate numerical information, and may reflect the need to accumulate greater amount of information for successful problem-solving (Fig. 2C,D). More generally, our results suggest that factors that may contribute to vulnerability in cognitive deficits are different in affected children and need to be taken into account while evaluating individual strategies for remediation in children at the low-end of the distribution of problem-solving skills.

Decision-making dynamics moderate the relationship between numerical problem-solving abilities and anterior insula engagement in ASD

The decision threshold ‘a’ significantly moderated the relationship between brain activation in the right anterior insula and numerical ability (Fig. 4A). Exploratory analysis revealed that the high threshold group consisted mostly of children with ASD, while the low threshold group consisted mostly of TD children (Fig. 4B). Notably, children with the highest decision thresholds demonstrated a significant positive association between numerical ability and right anterior insula activation (Fig. 4B), while children in the low threshold group showed a significant negative association (Fig. 4B). This analysis, albeit exploratory, reveals yet another source of individual variability in ASD and furthermore identifies a subgroup of children with ASD with similar cognitive profiles to TD children. Critically, our findings highlight a key role for the right anterior insula in accumulating information needed for making a decision and further demonstrates that this process contributes a significant source of variability in numerical problem-solving in children with ASD.

The right anterior insula is a major dysfunctional hub in children with ASD (Uddin et al., 2013a,b) and has been linked to their reduced cognitive flexibility (Menon, 2011). Specifically, the right anterior insula is a major causal hub for initiating control signals, switching from low to high-load cognitive states, and its functional and structural interactions with frontal and parietal systems have been linked to efficient numerical processing in neurotypical children and adults (Supekar and Menon, 2012). Taken together, these results point to the right anterior insula as a particular locus of vulnerability in a subtype of children with ASD and identifies an additional source of heterogeneity arising from aberrant decision-making dynamics.

Implications for remediation

Our study advocates for an individual differences approach to uncover the complex phenotype of ASD, rather than treating heterogeneity of data as a source of ‘noise’ (Kanai and Rees, 2011). More generally, our results suggest that unique behavioral and neural sources of variability may underlie the striking individual variation in cognitive skills and academic achievement in ASD (Assouline et al., 2012; Baron-Cohen et al., 2007; Cash, 1999; Chen et al., 2018; Chiang and Lin, 2007; Estes et al., 2011; Kim and Cameron, 2016; Mayes and Calhoun, 2003, 2008; Minshew et al., 1994; Oswald et al., 2016). Greater focus on these sources of variability may explain why individuals with ASD show accelerated learning in some tasks (Plaisted et al., 1998; Sears et al., 1994), but degraded learning in others (Froehlich et al., 2012; Klinger and Dawson, 2001), and why a higher proportion of individuals with ASD tend to achieve lower levels of post-secondary education, employment, and independent living, even relative to individuals with intellectual and learning disabilities (Newman et al., 2011; Troyb et al., 2014). Our findings also point to the important role of the anterior insula in modulating cognitive problem-solving skills and suggest that interventions focused on insular-mediated functions may be an important avenue for improving cognitive dysfunction in children with ASD.

Limitations

First, although effect sizes in our study fell within the medium to large range (0.2825 < Cohen’s f2 < 1.0715) (Cohen, 1988; Selya et al., 2012), with Bayes factors ranging from 2.46 to ∼ 3.28 × 102, future work with larger sample sizes is required to replicate and assess stability and robustness of our findings, particularly on the moderation effects of latent-decision variables on numerical skills and insular activity. Larger samples are also needed to better characterize subgroups of ASD children based on aberrant decision-making dynamics. Second, our sample was limited to children with normal IQ. Future studies that address hetereogeneity of ASD should also include children with low IQ and those who are minimally verbal. However, it is important to note that neuroimaging studies with low-functioning individuals with ASD are still difficult to conduct, especially with task-based fMRI. Third, our sample was primarily composed of male participants, which is consistent with prevalence rates in the disorder. However, it is important that future studies examine sex differences to assess whether heterogeneous profiles of cognitive abilities are moderated by sex. Finally, further studies are needed to determine whether behavioral, cognitive and neural sources of heterogeneity in numerical problem-solving in ASD are generalizable to other cognitive domains.

Conclusions

Our investigation of individual differences using a multi-level brain-behavior-cognitive approach reveals patterns of differences and unique profiles that are not detectable when averaging data across participants and comparing overall group differences. A particular emphasis on domain-specific neural sources of heterogeneity is warranted because none of the general cognitive (IQ, working memory), clinical, or even decision-making measures directly predicted individual differences in cognitive problem-solving in ASD. Our findings highlight the unique potential of combining brain activation and latent decision-making variables as more sensitive metrics for assessing sources of individual differences in ASD. This is profoundly important not only for defining and better characterizing the nature of ASD (Happé, 1999), but also for identifying areas of relative strength and weaknesses, which could in turn help developing targeted educational interventions (Barnett and Cleary, 2015) and facilitate long-term academic and professional success in affected individuals. Finally, our multilevel approach provides a novel methodology for characterization of heterogeneity in neurodevelopmental and learning disabilities.

Availability of data and material

Behavioral and fMRI data will be made available on the NIMH Data Archive upon publication. The genetic algorithm used in this study for matching the groups is available at: https://med.stanford.edu/content/dam/sm/scsnl/documents/groupmatch-scripts.tar.zip.

Funding

This work was supported by grants MH084164, HD059205, and MH121069 to VM, and a Postdoctoral Award from the Stanford Child Health Research Institute and the Stanford NIH-National Center for Advancing Training Sciences-Clinical and Translational Science Awards program (TR001085) to LC. TI is currently supported by grants from the (ANR-18-CE28-0010) and from the Université de Paris (IdEx).

Authors’ contributions

V.M., and T.I. designed the research; T.I. performed the research; T.I. L.C., S.M., C.A. and J.N. analyzed the data; and T.I., A.P., and V.M. wrote the paper.

Authors’ information

Teresa Iuculano, Institutional address: Centre National de la Recherche Scientifique & Université de Paris – UMR CNRS 8240, La Sorbonne, 46 rue Saint-Jacques, 75005 Paris, France. Past address: Stanford Cognitive and Systems Neuroscience Laboratory, Department of Psychiatry and Behavioral Sciences, Stanford University School of Medicine, 401 Quarry Rd, Stanford, CA, 94305, USA.

Declaration of Competing Interest

The authors report no biomedical financial interests or potential conflicts of interest. All authors of the manuscript have read and agreed to its content.