Abstract
Aphasia is a prevalent cognitive syndrome caused by stroke. The rarity of premorbid imaging and heterogeneity of lesion obscures the links between the local effects of the lesion, global anatomic network organization, and aphasia symptoms. We applied a simulated attack approach in humans to examine the effects of 39 stroke lesions (16 females) on anatomic network topology by simulating their effects in a control sample of 36 healthy (15 females) brain networks. We focused on measures of global network organization thought to support overall brain function and resilience in the whole brain and within the left hemisphere. After removing lesion volume from the network topology measures and behavioral scores [the Western Aphasia Battery Aphasia Quotient (WAB-AQ), four behavioral factor scores obtained from a neuropsychological battery, and a factor sum], we compared the behavioral variance accounted for by simulated poststroke connectomes to that observed in the randomly permuted data. Global measures of anatomic network topology in the whole brain and left hemisphere accounted for 10% variance or more of the WAB-AQ and the lexical factor score beyond lesion volume and null permutations. Streamline networks provided more reliable point estimates than FA networks. Edge weights and network efficiency were weighted most highly in predicting the WAB-AQ for FA networks. Overall, our results suggest that global network measures provide modest statistical value beyond lesion volume when predicting overall aphasia severity, but less value in predicting specific behaviors. Variability in estimates could be induced by premorbid ability, deafferentation and diaschisis, and neuroplasticity following stroke.
SIGNIFICANCE STATEMENT Poststroke, the remaining neuroanatomy maintains cognition and supports recovery. However, studies often use small, cross-sectional samples that cannot fully model the interactions between lesions and other variables that affect networks in stroke. Alternate methods are required to account for these effects. “Simulated attack” models are computational approaches that apply virtual damage to the brain and measure their putative consequences. Using a simulated attack model, we estimated how simulated damage to anatomic networks could account for language performance. Overall, our results reveal that global network measures can provide modest statistical value predicting overall aphasia severity, but less value in predicting specific behaviors. These findings suggest that more theoretically precise network models could be necessary to robustly predict individual outcomes in aphasia.
Introduction
Aphasia is one of the primary cognitive symptoms following left hemispheric strokes, affecting 180,000 new individuals a year in the United States (Ovbiagele and Nguyen-Huynh, 2011). Despite decades of research, the brain basis of aphasia outcomes and recovery remain only partially understood. The majority of stroke research has focused on the relationship between the regional anatomic influences of stroke on cognitive symptoms and outcomes (Bates et al., 2003; Rorden et al., 2007; Meyer et al., 2016; Mirman and Thye, 2018). More recently, investigators have studied the relationships between individual anatomic tracts, the topology of complex brain networks, the connectome (Hagmann, 2005; Sporns et al., 2005; Bassett and Sporns, 2017; van den Heuvel and Sporns, 2019), and behavior (Yourganov et al., 2016; Gleichgerrcht et al., 2017; Del Gaizo et al., 2017; Fox, 2018; Gollo et al., 2018).
Poststroke, the remaining neuroanatomy maintains cognition and supports recovery. Anatomical network connectivity in the lost and residual (spared) connectome after stroke is related to behavior (Kunimatsu et al., 2003; Heller et al., 2005; Mukherjee, 2005; Crofts et al., 2011; Kim and Jang, 2013; Forkel et al., 2014; Yourganov et al., 2016; Gleichgerrcht et al., 2017). In particular, single-connection analyses have demonstrated that regions with links to classical hub regions such as the temporoparietal junction are crucial for overall language function assessed with clinical measures (Yourganov et al., 2016). Strokes that directly impact network hubs disproportionately lead to global cognitive deficits poststroke on tasks that place significant semantic or language-production demands on patients (Warren et al., 2014). In addition, cognitive outcomes are associated with the preservation of the brain's modular configuration, the tendency for brain regions to group into well-connected clusters (Marebwa et al., 2017; Schlemm et al., 2020). Overall, these findings suggest that the role of single regions and their connections in network topology, as well as overall network topology, are related to stroke symptomatology.
A primary difficulty in assessing stroke-induced effects on network topology is that researchers often lack premorbid data within-subjects, leading them to rely on cross-sectional analyses. This is problematic when comparing behaviors in stroke to that observed in healthy subjects. Lesions interact with other factors about the individual, such as their development, demographics, and brain organization. As a complement to observing the consequences of stroke and other types of brain injury, “simulated attack” models are computational approaches that apply virtual damage to the brain and measure their putative consequences (Kaiser et al., 2007; Joyce et al., 2013). These models can be used to systematically quantify the influences of damage to regions and connections on brain network organization. Examining brain networks of individual participants with brain lesions conflates the direct effects of lesions with idiosyncratic individual differences that may have predated the lesions with the effects of changes in the brain after the lesion. In contrast, simulated attack models can isolate the direct effect of lesions on brain networks, quantifying how a given lesion effects the network properties of a range of otherwise typical brains. In the context of aphasia because of stroke where sample sizes comprise a few tens of individuals, simulated attacks also allow us to estimate a range of expected brain-behavior relationships that can be directly attributed to the direct effects of the lesion on the brain network, leaving aside individual differences that may have predated the stroke or changes in the brain after the stroke. After simulating damage, hypotheses about network robustness, cognitive resilience, and recovery can be tested in the residual connectomes (Aerts et al., 2016). Measures characterizing the disconnectivity of circuits and networks (Gleichgerrcht et al., 2017), the overall efficiency of the network (Van Den Heuvel et al., 2009; Bullmore and Sporns, 2012; Ajilore et al., 2014), and the balance between local and distributed processing (small-worldness; Bassett and Bullmore, 2017) could relate to behavioral performance. In addition, the deviation in these properties from that expected in a comparison model of healthy subjects might also characterize variation in resilience to cognitive decline.
To examine these possibilities in aphasia severity, we used probabilistic diffusion tractography to create anatomic connectomes in 39 subjects with left-hemispheric strokes. Then, we computed measures that quantify five network properties of anatomic connectivity poststroke thought to be related to the integrity of observed topology. Using a simulated attack model, we computed the effects of each stroke's specific pattern of connection losses to quantify its effects on the whole brain and intra-left hemisphere connections in a sample of healthy subjects. Then, we computed models estimating the behavioral variance measured with clinical language measures accounted for by simulated anatomic network measures. This technique allowed us to obtain confidence intervals for the strength of brain-behavior relationships between lesioned network topology and behavior. Above and beyond lesion volume, we hypothesized that total edge weights, network modularity, global shortest path length, higher local clustering, and small-worldness would be related to better language performance. We further hypothesized that the global network measures would be more related to global measures of the severity of language deficits than factor scores representing specific lexical, auditory comprehension, phonology, and cognitive/semantic deficits.
Materials and Methods
Subjects
The Western aphasia battery, revised (Kertesz, 2007) was obtained for each individual with stroke. In addition, participants with stroke performed a broader battery of tasks previously described in detail (Lacey et al., 2017). To reduce the scores from the battery, a principal components factor analysis was performed in SPSS 25 using the individual test scores from the Western Aphasia Battery Revised (WAB-R) and the other battery tasks on the 59 participants with stroke from a larger ongoing study who were able to provide complete behavioral data. Factor analysis was performed on the correlation matrix, factors were extracted based on the standard cutoff of eigenvalue >1, and Varimax rotation with Kaiser normalization was applied to achieve orthogonal factors. Consistent with a previously reported factor analysis on a subset of these participants (Lacey et al., 2017), the factor analysis revealed four factors cumulatively accounting for 83.7% of variance in the scores that we interpreted to reflect lexical production, auditory comprehension, phonology, and cognitive and semantic aspects of behavior (see Table 1). Factor scores for each participant were calculated using the regression method.
Factor loadings for the behavioral data
All procedures were approved in a convened review by Georgetown University's Institutional Review Board and were conducted in accordance with the guidelines of the Institutional Review Board/Human Subjects Committee, Georgetown University. All participants volunteered and provided informed consent in writing before data collection.
Behavioral data
The WAB-R (Kertesz, 2007) was obtained for each individual with stroke. In addition, participants with stroke performed a broader battery of tasks previously described in detail (Lacey et al., 2017). To reduce the scores from the battery, a principal components factor analysis was performed in SPSS 25 using the individual test scores from the WAB-R and the other battery tasks on the 59 participants with stroke from a larger ongoing study who were able to provide complete behavioral data. Factor analysis was performed on the correlation matrix, factors were extracted based on the standard cutoff of eigenvalue >1, and Varimax rotation with Kaiser normalization was applied to achieve orthogonal factors. Consistent with a previously reported factor analysis on a subset of these participants (Lacey et al., 2017), the factor analysis revealed four factors cumulatively accounting for 83.7% of variance in the scores that we interpreted to reflect lexical production, auditory comprehension, phonology, and cognitive and semantic aspects of behavior (see Table 1). Factor scores for each participant were calculated using the regression method.
Neuroimaging
Diffusion images were acquired on a Siemens 3.0T Magnetom Trio for all subjects along with a T1-weighted 1-mm resolution MPRAGE anatomic scan at each scanning session as part of a larger imaging protocol. We used a high-angular resolution diffusion imaging (HARDI) acquisition scheme with a maximum b-value of 1100 (80 dirs, 10 b = 0; 10 b = 300; 60 b = 1100) and a 2.5-mm isotropic voxel size. We used a transversal acquisition of 55 axial slices with the following parameters: repetition time (TR) = 7.5 s; echo time (TE) = 87 ms; field of view (FoV) = 240 × 240, 138 mm, matrix = 96, total acquisition time of 10:00. MPRAGE scans were collected with TR = 1900 ms, TE = 2.52 ms, 176 sagittal slices with 0.9 mm slice thickness, FoV = 240 × 240, matrix = 256, inversion time (TI) = 900 ms and flip angle = 9°, total acquisition time of 5:34.
Anatomical image imputation
We conducted full-brain tractography with techniques that reduce tractography artifacts. To achieve this, tractography must be constrained anatomically to seed or terminate streamlines at the gray matter/white matter border (Smith et al., 2012; Song et al., 2014; Jeurissen et al., 2019). However, identifying tissue types in stroke cases is problematic because of the abnormal signal intensity at the gray-white matter border because of encephalomalacia and other pathologic effects. These intensity differences and the tissue lost because of stroke can be problematic for surface-based registration and parcellations such as those used in the current study. We resolved these issues by imputing estimates of healthy tissue in the anatomic T1 images from subjects in the stroke sample.
We imputed anatomic images in two steps. In step 1, lesioned voxels were identified using lesion tracings provided by an experienced cognitive neurologist (coauthor PET). We then flipped the lesioned brain along the left-right plane and registered the flipped brain onto the nonflipped brain. Next, we filled the lesioned area with the healthy tissue of the homotopic contralesional hemisphere. This procedure can leave visible marks of the filled area because of the sudden change in signal, which may cause artifacts when identifying tissue types. Thus, we imputed a new brain with highly similar morphologic features as the original subject's brain. The imputation procedure used the morphologic structure of the reference image (the filled brain from step 1) and the voxel values of a set of healthy control images to produce a new image. Each voxel value in the new image was determined by combining the values from all the healthy images using ANTs' joint image fusion procedure, where images more similar around the voxel of interest received more weight (similar to multiatlas label fusion; H. Wang et al., 2013). We conducted a search of the optimal number of healthy brains and the optimal radius of similarity around each voxel to obtain the best result. We obtained an optimal outcome with 22 healthy brains and a radius of 1 (i.e., a single layer of voxels around each voxel is used to check the similarity between images and assign weights to healthy images). We inspected the resulting imputed image to make sure there were no artifacts; none were found. Importantly, the nonlesioned gyri and sulci in the original image followed the gyri and sulci of the imputed image without any visible deviation.
We performed all the imputation procedures in ANTs (v. 2.2.0). Before any processing, all images were skull-stripped (antsBrainExtraction.sh), corrected for magnetic field inhomogneity (N4BiasFieldCorrection), and denoised with an edge preserving algorithm (PeronaMalik, denoising amount: 0.7, iterations: 10). We added back the lesion mask to the brain mask after skull-stripping to ensure that the lesion area was included in the imputation. Each imputation required one registration of the flipped image and 22 registrations of the healthy brains onto the filled image. We conducted all registrations using the SyN nonlinear algorithm (Avants et al., 2011) with cost function masking to remove the lesion mask from consideration during the registration computations (Brett et al., 2001).
Diffusion tractography
We used MRtrix3 to process the diffusion data (Tournier et al., 2019). to denoise the diffusion images (function: dwidenoise -extent 9,9,9), correct for motion and eddy currents (function: dwipreproc), and correct for field inhomogeneity (function: dwibiascorrect). We then computed response functions for multiple tissues using the tissue information available in the DWI data (function: dwi2response dhollander). Finally, we computed fiber orientation distributions (FOD) via a multishell multitissue constrained spherical deconvolution (function: dwi2fod msmt_csd; Jeurissen et al., 2014).
To find the GM/WM tissue, we applied tissue classification to the imputed anatomic image (function: 5ttgen fsl) and brought the tissue information into DWI space after registering the original (lesioned and imputed) T1w image of the subject onto the mean b = 0 image (function: antsRegistration, order: translation, rigid, SyN) and applying the transformations to the tissue types. We performed white matter tractography by seeding 15 million streamlines probabilistically from the white matter based on estimated fiber densities (tckgen algorithm: iFOD2, step: 1 mm, min length: 10 mm, max length: 300 mm, angle: 45°, seeding: dynamic, backtracking allowed, streamlines cropped at GM/WM border).
To decrease tractography artifacts, spherical deconvolution informed filtering of tractograms (SIFT2) was conducted to determine the relative apparent fiber density associated with each streamline (Smith et al., 2015a, b). Each subject's anatomic connectivity was then quantified through FA connectomes that accounted for the apparent fiber densities associated with each streamline. The contribution of each streamline's FA value to an edge's mean FA was weighted by the streamline's SIFT2 cross-sectional multiplier, which represents the streamline's relative apparent fiber density. Specifically, the mean FA along the path of each streamline was calculated using the subject's tractogram and FA image derived from the DWI.
Network construction
Anatomical scans (imputed in the case of patients) were segmented using FreeSurfer (Fischl, 2012) and parcellated using the connectome mapping toolkit (Cammoun et al., 2012). A parcellation scheme including N = 234 regions was registered to a single b = 0 volume from each subject's native-space diffusion spectrum imaging (DSI) data. The b = 0 to MNI voxel mapping produced via Q-Space Diffeomorphic Reconstruction (QSDR) was used to map region labels from native space to MNI coordinates so that individual subject data could be combined and analyzed in a shared standard space. To extend region labels through the gray-white matter interface, the atlas was dilated by 4 mm (Cieslak and Grafton, 2014). Dilation was accomplished by filling nonlabeled voxels with the statistical mode of their neighbors' labels. In the event of a tie, one of the modes was randomly selected. Each streamline was labeled according to its terminal region pair. From these data, we constructed an anatomic connectivity adjacency matrix, A whose element Aij represented either (1) the weighted mean fractional anisotropy (FA) of the streamlines connecting that pair of regions based on the underlying fiber densities or (2) the streamline connectome edge values representing the sum of weights associated with each streamline computed by SIFT2 (Hagmann et al., 2008; Smith et al., 2015b).
To enable group analyses, intersubject apparent fiber density and connection density normalization was conducted. Specifically, each subject's connectome was multiplied by the geometric mean of the ratio of the individual's response function size at each b value to the group average response function size at each b value. Individual differences in white matter b0 intensity were accounted for by multiplying each connectome by the ratio of the mean median b0 value within the subject's white matter mask to the grand mean median b0 value for the whole group. Intersubject connection density normalization was then achieved through scalar multiplication of each connectome by the subject's “proportionality coefficient” derived by SIFT2, denoted by μ, which represents the estimated fiber volume per unit length contributed by each streamline.
To visualize the effects of lesions on parcels in the Lausanne anatomic atlas, we registered the lesion masks to each individual's T1 image (the same space as the Lausanne parcel registration). We computed whether the lesion intersected > 0 voxels in that parcel, and counted the number of subjects at which that parcel was intersected by the lesion. See Figure 1 for a visualization of the distribution of lesions across subjects and Figure 2 for a summary of the tractography pipeline including the imputation.
Distribution of lesioned parcels across subjects with stroke. Lesions mapped prominently to parcels in left perisylvian regions with decreasing frequency in the superior, inferior, anterior, and posterior directions. The degree of red is proportional to the number of subjects with lesions at that location. N = the number of subjects with a lesion at that parcel label.
Schematic of stroke imputation and diffusion tractography. A, Processing scheme for healthy subjects. Diffusion tractography was computed in subjects' native space, and the Lausanne multiscale parcellation was fit to subjects' anatomic T1 images. Connectomes were defined based on the FA or streamline counts of the edges connecting each region pair and advanced to analyses. B, The processing scheme for stroke subjects was the same as the healthy subjects with an additional preprocessing step. Specifically, the anatomic T1 image was imputed using the stroke subject's right hemisphere and healthy subjects' data to estimate the prelesion T1 anatomic image. The parcellation was computed on this imputed anatomic image to guide connectome extraction through the same regions as the controls.
Connectome edge inclusion mask
Given well-described false positives issues in diffusion tractography, it is difficult to ensure that every individual streamline is valid in the absence of ground-truth data (Maier-Hein et al., 2017). We only permitted edges to participate in our analysis if they were present in 100% of the healthy control sample. This procedure ensured that any changes in topology observed in the stroke sample were likely to be driven by lesion-related effects rather than spurious patterns attributable to unreliable tractography findings between each pair of parcels. All triangle matrices produced via MRTrix3 were symmetrized across the diagonal before network analyses.
Simulated attack
Our goal was to simulate connectome attacks using estimates of the consequences of real strokes on the connectome relative to the connectivity observed in control subjects. To identify a set of potentially lesioned edges for each subject, we compared the edge values from each stroke subject's observed connectome to those observed in healthy subjects (Fig. 3). Because the definition of lesioned tissue depended on a binary threshold, we computed the simulated attacks at thresholds of 2, 3, and 4 SDs below the mean FA or the log of the streamlines (to account for lognormal edge distributions, see Smith et al., 2015a; Zhang et al., 2018; Yeh et al., 2019; Rosen and Halgren, 2021) relative to the control sample. For each threshold, a mask was created for all lesioned edges for each stroke subject. Next, we applied the edge lesion mask to each connectome in the control sample in addition to the same stroke subject. By applying the mask to the stroke subject, we ensured that the number and configuration of edges included in the analysis was equal between each stroke subject and the controls. Finally, each connectome measure was averaged across the thresholds to obtain a representative value for each subject. We computed the pairwise similarity between lesion masks across thresholds using intrasubject Pearson's correlations for each threshold pair. Across the thresholds, we found very high consistency for the streamline networks [R(mean, SD) = 0.88 (0.14), p(mean, SD) ≪0.001 (≪0.001)] and FA networks (R(mean, SD) = 0.82 (0.13), p(mean, SD) ≪0.001 (≪0.001)] across the intrasubject correlations. After simulating the attacks, an empirical distribution of the expected effects of lesions on network measures was obtained by computing the network measures (described below) for each possible control subject-to-lesion pairing.
Schematic of network lesion masking. A, Top, Each element Ai j from each subject with aphasia (Strokei) was compared with (bottom) the observed values in all control subjects (Control1 to Controln). B, The elements with FA or streamlines 2, 3, or 4 standard deviations (SDs) less than controls were labeled as lesioned edges. C, Then, the lesion mask was applied to the stroke subject and all control subjects, and the resulting networks were advanced to connectomic analyses.
Connectome measures
We focused on five global network measures that are thought to characterize overall network communication (Fig. 4). We examined (1) the network edge weights: the sum of observed edge values, which could be related to overall network intactness and can drive other global network measures. The sum of edge weights served as a network proxy that estimates the overall integrity of edges in the residual connectome. Intuitively, edge weight should be positively related to function in the context of stroke and aphasia.
Schematic of network measures. Moving Left to Right from the Top, We began with the (1) sum of edge weights in each network as an overall metric capturing the density of the networks, including any edges lost because of stroke. Four other measures were of interest. (2) Modularity: measures the extent to which nodes in the network are grouped into modules (sometimes, “communities”) as a function of highly-connected nodes. (3) Global efficiency: one long path is represented by the set of consecutive edges highlighted in green. (4) Transitivity: one possible triplet's edges are represented in green. (5) SWP: involves a high degree of local clustering (represented by the set of nodes connected by purple edges) and short path lengths (e.g., higher weights along the green path represents a shorter network path between the prefrontal and occipital nodes).
There are additional measures that can exhibit complex interactions with one another and behavior, and our goal was to benchmark these relationships with our simulated attack approach. Thus, we examined other commonly examined characteristics of the overall network topology that have been linked to variability in global cognitive performance. The (2) modularity of a network is thought to support local computations within tightly-connected subgroups of nodes within a network (Meunier et al., 2009; Alexander-Bloch et al., 2010; Marebwa et al., 2017; Siegel et al., 2018) and can predict intervention-related cognitive plasticity (Gallen and D'Esposito, 2019). It is computed to estimate the relative within-module connectivity of a network relative to between-module connectivity with the modularity value Q (Reichardt and Bornholdt, 2006; Leicht and Newman, 2008).
In addition, to examine overall network processing efficiency, we examined (3) network global efficiency (Iturria-Medina et al., 2008; Griffa et al., 2013; Lawrence et al., 2014; Berlot et al., 2016; Beare et al., 2017). Here, global efficiency is a weighted measure that quantifies average inverse shortest path length across the connectome for all pairs of nodes (Latora and Marchiori, 2001). It is inversely related to the network measure path length. A pair of nodes with a short path length are connected by sequences of stronger edges (higher FA or streamline values). Intuitively, stronger connections between nodes can theoretically represent the strength of information flow between regions. Thus, the average path length of a network represents the extent to which all pairs of nodes are associated via short hops through the network. Accordingly, networks with high global efficiency are thought to have increased long-distance information processing capacity across all nodes mediated by short paths.
In brain networks, local clustering among u-fibers constitute most of the brain's white matter (Schüz and Braitenberg, 2002). Accordingly, short-distance, local clustering is another important aspect of information processing in brain networks. Therefore, we also examined network (4) transitivity (Lo et al., 2011; Prasad et al., 2013; Llufriu et al., 2017). Transitivity is the ratios of triangles, which are groups of three nodes connected by three non-zero edges, to all possible triplets (sets of three nodes). Networks typically have many more triplets than they have triangles. Therefore, greater transitivity means that there are more local clusters in a network. Networks with high transitivity are thought to have increased local communication efficiency (Rubinov and Sporns, 2010).
Finally, healthy brain networks are characterized by an optimal use of available anatomic connections to support short path lengths and high clustering, which is often referred to as small-worldness (Watts and Strogatz, 1998; Stam, 2014). Small-world brain networks are thought to confer many of the processing advantages that support diverse and dynamic cognitive functions (Bassett and Bullmore, 2017). To investigate this property, we used a robust measure of small-worldness, (5) small-world propensity (SWP; Muldoon et al., 2016). SWP is a weighted metric for small-worldness that accounts for networks of different densities, standardizing the measure against individualized network null models. This technique makes SWP appropriate for measuring small-worldness in weighted networks by mitigating the network density-dependence of other measures.
Statistical analyses
The effects of lesions on connectome measures and behavior
First, we examined whether observed and simulated strokes had significant effects on each network measure in the whole brain and within the left hemisphere. We computed Welch's t tests assuming unequal variances using Satterthwaite's approximation for degrees of freedom for each measure against those observed in the control subjects corrected for multiple comparisons at an α level of 0.05. Then, we used bootstrapping to estimate the proportion of network measure sample means from the simulated attacks that fell within the range of the observed lesion for each measure. We used this technique because we intended the simulation to sample from all lesion-control subject pairs to yield a distribution of possible lesion profiles in a much larger simulated sample. Specifically, we performed 10,000 resamplings with replacement of 39 subjects and quantified the proportion of simulated attack sample means for each network measure (FA and streamline) and each size network (whole brain or left hemisphere). In these and similar permutation analyses, there is always a trade-off between establishing a null distribution and computing time and complexity, and inferential significance. In this case of these analyses, there were 39 subjects with stroke and 36 healthy controls, totaling 1404 unique lesion-subject pairings for the simulated nulls. We chose a permutation set of 10,000 re-samplings per test at an order of magnitude higher than the possible unique pairs, allowing each pair to contribute to multiple re-samplings. This sample size allowed us to obtain effect size and variance estimates for univariate models with a maximum p-value of 0.0001 per test. Finally, to estimate the behavioral variance accounted for by lesion volumes, we fit separate linear regression models using lesion volume as an independent variable and either WAB-AQ, the factor sum score, or each behavioral factor score corrected for multiple comparisons at an α level of 0.05.
Preparing network measures to identify behavioral variance beyond lesion volume
Our objective was to obtain and present an empirical estimate for the full range of possible lesion-behavior relationships observed in the real, simulated, and null analyses.
For all connectome analyses, we were interested in the total variance accounted for using all five network measures for the observed and simulated data. Before analyses, we tested the network measures for violations of normality with the Kolmogorov–Smirnov test. To correct for skewed distributions in the observed statistics, we used a log-transformation for global efficiency and SWP. In the simulated attack statistics, we observed negative values and skew for each statistic; thus, we added a constant value of 1 to each measure before a log-transformation. Finally, network measures were standardized using z scores before all analyses.
To test the hypothesis that real and simulated measures of network topology were related to behavioral scores beyond lesion volume, we first computed separate linear regressions using stroke subject's lesion volume as the independent variable and each behavioral scale as the dependent variable. Then, we used the residualized behavioral scores as the dependent variable for all connectome-behavior analyses. We performed the same procedure for each network measure to mitigate any remaining influences of lesion volume effects.
Network-behavior relationships in observed, simulated, and null regression models
Next, all analyses associating network measures with behavioral scores were performed using linear regression in R statistical software (R Core Team, 2013). Our simulated attacks broke the relationship between stroke subjects and behavior by randomly sampling residual anatomic connectomes after simulating strokes in control subjects, but preserved the relationship between each behavioral score and simulated lesion. Specifically, in the simulated attacks, the stroke connectome lesion mask was defined from a subject with a real stroke. This lesion mask was applied to each healthy control. Each healthy control was not the individual who suffered a real stroke. Thus, the sample of simulated attacks represented a range of possible relationships between that specific lesion profile and the behavior but did not include the original (real) connectome associated with the stroke. In each case, the distribution was compared with the observed relationship between the real anatomic connectomes with a stroke and the behaviors of interest. The null model completely randomized the relationships between simulated lesions and behaviors.
Specifically, we computed linear regression models for (1) the observed stroke network topology, (2) the simulated stroke network topology (10,000 permutations per lesion edge threshold), and (3) a randomized shuffling of all simulated network measures against the behavior (10,000 permutations per lesion edge threshold). We examined the effects of observed and simulated lesions on each of the topological measures in the whole brain and within intra-left hemisphere connectomes (i.e., only the whole-brain connectomes included interhemispheric and right-hemispheric fibers). Then, we computed the relationships between network measures and the behavioral scores for the whole brain and intra-left hemisphere connectomes.
In analyses of the observed data (i.e., data from subjects with real strokes), we used each of the five network measures as independent variables (z-scored across subjects) and each behavioral score as a dependent variable (raw WAB scores, the factor sum scores, or one of the four behavioral factor scores) in separate linear regression models. Because the network measures and behavioral scores were the residuals obtained after regressing out the influence of lesion volume, we obtained specific parameter estimates for each network measure and the total variance accounted for in the models (R2 value) beyond lesion volume.
To obtain estimates for network-behavior relationships in the simulated attack, we computed the linear regression models with the same dependent behavioral variables, but with independent network variables sampled from the control-lesion pairings for the simulated attack (z-scored across subjects). Specifically, in 10,000 permutations, we randomly assigned each stroke lesion to a healthy brain from the control sample while preserving the link between that lesion and the behavioral outcome. We then computed each network measure on that sample of simulated attacks, and fit a regression model. Across the 10,000 models, this approach provided a full representation of the absolute minimum, maximum, and of the predictive value (R2) of the anatomic connectomes beyond lesion volume. In addition, we obtained the range of β weights for each network measure in the simulations to reveal their relative contributions to the prediction.
In high-dimensional data analyses such as these, it is often helpful to have an empirical null distribution to contextualize the models of interest. By fully randomizing the relationships between the simulated and network measures and behavior (i.e., breaking the lesion-behavior pairing in the simulated strokes), we were able to create a distribution of the expected variance accounted for (R2) if the data were to be completely randomized. In this kind of null permutation analysis, we would expect that models could trivially account for more variance than 0 by chance. Further, the null could include permutations equaling or similar to real lesion-behavior network pairings, potentially meeting or exceeding the variances obtained from the observed or simulated attack. If the observed or simulated models exhibited effects that were higher than the null's central tendency, it would increase our confidence that the network topology across the range of simulated attack outcomes is nontrivially related to behavior above and beyond lesion volumes.
Our primary goal was to test whether the simulated attacks based on observed lesions differed from the null distribution. After computing the R2 value for each null and simulated model in each permutation, we used a 2 (simulated vs null) × 6 (behavioral variables) ANOVA to test the effects of (1) the simulated attack relative to the null permutations and (2) behavioral domain on the estimated R2 values.
Results
The effects of lesions on network measures
The observed strokes influenced network topology for each measure. In FA networks in the whole brain, we observed a reduction in mean edge weight and network efficiency and increased SWP. In the left hemisphere, we observed all of these effects in addition to reduced modularity (Fig. 5). In the streamline networks in the whole brain, we observed reduced mean edge weight, network efficiency, and transitivity as well as increased modularity. In the left hemisphere, we observed reduced edge weights, modularity and increased transitivity (Fig. 5; Table 2). The network measures from the simulated strokes were similar to those observed after real strokes (see Table 3), suggesting that they were reasonable approximations of stroke effects on the connectomes.
Two-samples t tests assuming unequal variances comparing network measures between the control and stroke samples for the whole brain FA and streamline networks (df = Satterthwaite's approximation for samples with unequal variances)
The proportion of bootstrap sample means of simulated attacks compared the range of observed lesion values
The effects of stroke on network measures in observed and simulated attack connectomes. The leftmost column of each plot facet shows the network statistic observed in controls, followed by that observed in strokes, then the simulated attacks. Network measures are presented in their raw (untransformed) values before inclusion in network-behavior analyses. Asterisks indicate a significant Welch's two-sample t test between the control and stroke network measures at p < 0.001 (a stringent threshold after Bonferroni correction for 40 total tests in FA and streamline data). The top and bottom edges of the boxes represent the 25th and 75th percentiles, respectively. SWP = small world propensity.
Relationships between observed and simulated connectome topology and aphasia-related behaviors
Lesion volume accounted for ∼44% of WAB-AQ, 53% of the factor sum, and 10–16% of the variance in factor scores (Table 4). In addition, lesion volume was negatively correlated with FA network edge weights, modularity, and efficiency, and positively associated with small world propensity. Lesion volume was negatively correlated with streamline network edge weights, efficiency, and positively associated with transitivity and small world propensity (Table 5).
The relationships between lesion volume and behavioral scores
The relationships between lesion volume and network measures
R-values represent Pearson's correlation coefficients between lesion volume and the network measure.
We regressed lesion volume from the behavioral and network data to obtain the additional unique variance between the network measures and behavior. The full model results are presented in Figure 6 For each of the FA and streamline whole brain and intra-left hemisphere networks, the variance accounted for (R2) by the simulated attack networks was greater than the nulls in the omnibus ANOVA [see Table 6 for condition-wise marginal means, which quantify the difference in the observed vs the null (R2) across behaviors].
Marginal means for the simulated attack relative to the null for FA and streamline networks in the whole brain and left hemisphere
Network measures and behavioral variance in FA and streamline networks. Each plot facet illustrates histograms of the simulated null, histograms of the simulated attack distributions, and the observed R2 with solid vertical lines. Asterisks indicate significant post hoc Welch's one-tailed t tests assuming unequal variances comparing the R2 values in the simulated attacks to the null distribution at p < 0.001. Daggers indicate cases where the observed R2 value was outside the range obtained in the simulated attack models. See Extended Data Figures 6-1, 6-2, 6-3, and 6-4 for specific effects of behavior on variance explained in the whole brain and left hemisphere.
Figure 6-1
Supplementary Extended Data Figures. Download Figure 6-1, DOCX file.
Figure 6-2
Supplementary Extended Data Figures. Download Figure 6-2, DOCX file.
Figure 6-3
Supplementary Extended Data Figures. Download Figure 6-3, DOCX file.
Figure 6-4
Supplementary Extended Data Figures. Download Figure 6-4, DOCX file.
In addition, there were main effects of behavior on the model R2 values. The central tendency of the null models revealed that network measures would be expected to account for nearly 10% of the behavioral variance on global or specific language performance at random (i.e., when the link between the lesion being simulated and behavioral score was broken). Outperforming the null, the simulated attack models generally suggested that ∼20% of global aphasia outcomes measured with the WAB-AQ could be accounted for by residual anatomic network topology. Among the language behavioral factor scores, lexical processing exhibited the strongest relationships with network topology, at 20% or more variance accounted for by the whole brain or left hemisphere network measures. In most cases, the observed R2 estimate was within the range estimated in the simulated attacks. Exceptions were observed in several cases, and more frequently in FA networks, where observed stroke estimates were outside the simulated estimates for the whole-brain lexical, phonology, and cognitive/semantic factors and the left-hemisphere cognitive/semantic factors.
β Weights for specific network measures
We additionally obtained the β values for the simulated attack models to observe which measures contributed the most weight to model R2. β Weights for the whole brain and intra-left hemisphere models are illustrated in Figures 7. In the whole-brain models, edge weights and efficiency were most consistently associated with higher βs across behaviors, with some variation across the individual factor scores. Within the left hemisphere, edge weights and efficiency remained relatively stronger contributors to the global WAB and factor score sum behavioral measures. Among the four factor subscores, network measure β weights exhibited more variation across specific factors.
Simulated attack estimated beta weights for each network measure in FA and streamline networks. Each plot represents the range of betas obtained from the simulated attack models for the network measure. The top and bottom edges of the boxes represent the 25th and 75th percentiles, respectively. W. = edge weights; Mod. = modularity; Eff. = efficiency; Trans. = transitivity; SWP = small world propensity. See Extended Data Figure 7-1 for observed and simulated model betas.
Figure 7-1
Supplementary Extended Data Figures. Download Figure 7-1, DOCX file.
Discussion
In subjects with left-hemispheric strokes, we used an anatomic network simulated attack analysis to examine the relationships between topological network measures, a widely used clinical aphasia outcome measure (the WAB-AQ), and dimensional factor scores of language performance. We found that (1) simulating lesions in healthy connectomes can provide a good approximation of the effects of real lesions on anatomic networks, (2) the network properties of the simulated lesions can explain variance in behavior above and beyond lesion size, (3) in most cases observed lesions do not explain more variance in behavior than simulated lesions, and (4) the relationships between simulated and observed networks and behavior varied across language functions.
In general, lesions in the left hemisphere disrupted several measures of global topological organization. In whole brain and left hemispheric FA networks, stroke reduced the overall edge weight and efficiency relative to controls. In contrast, SWP tended to increase poststroke, consistent with reduced network efficiency relative to relatively unaffected global clustering (transitivity; Muldoon et al., 2016). Interestingly, modularity was decreased in the intrahemispheric connectomes. The relative increases in modularity in the whole brain were likely driven by the absorption of residual left-hemispheric networks into right hemispheric homotopic communities via interhemispheric fibers mediated through the corpus callosum. These patterns were similar in the simulated attacks, suggesting that the simulations were reasonable approximations of real stroke effects. A minority of sampled simulated cases in the left hemisphere had SWP values averaging outside the range of stroke subjects, potentially reflecting neuroplastic or premorbid differences in the hemispheric balance between long distance and local communication. Across the analyses, as lesion volume increased, efficiency decreased (meaning that path lengths increased) along with increases in transitivity. SWP quantifies clustering relative to path length in null models. Thus, our findings suggested that the positive association between lesion volume and SWP is because of the increase in transitivity with increasing lesion size outpacing the increase in path lengths.
As the focus of several previous studies, modularity changes in persons with stroke could vary based on the location of lesions. For instance, left hemisphere anatomic modularity has been found to increase in subjects with upper limb motor deficits (Cheng et al., 2019), and increased left hemisphere anatomic modularity have been associated with more severe chronic aphasia (Marebwa et al., 2017). In contrast, reduced modularity in functional connectomes observed in multiple stroke phenotypes (Corbetta et al., 2018) has been shown to partially recover in the transition from the acute to chronic phase (Lim and Kang, 2015; Siegel et al., 2018; Schlemm et al., 2020). Anatomical connections and network topology predict region-to-region functional connectivity (Griffis et al., 2020), and it will be important to clarify how specific lesion distributions interact with anatomy, and joint anatomy-function relationships. For instance, sensorimotor cortices are highly interconnected within each hemisphere, and precentral regions are often revealed to participate in the brain's anatomic hub system (Van Den Heuvel and Sporns, 2011; McColgan et al., 2015). Thus, disrupting sensorimotor regions and their connections is likely to enhance the modularity of the remaining intrahemispheric network. Reduced modularity was not strongly related to the language behaviors examined here, suggesting that this measure is not uniquely a feature of language function. Rather, individual patterns of region or connection damage could drive overall changes in topology that influence modularity. Overall, the topology of the left hemisphere was sparser and more tightly clustered, which is often thought to limit the general ability for a network to transmit information and reduce the interference between competing demands on the network because of the loss of specialized processing modules (Stam, 2014). However, this change in topology is not guaranteed to increase modularity, which is computed with respect to the configuration of clusters rather than overall network clustering.
In the simulated attack regression models associating brain network measures with behavior, we found that anatomic network measures accounted for behavioral variance above and beyond the effects of lesion volume. Unsurprisingly, lesion volume accounted for a moderate amount of variance in the global aphasia measures (WAB-AQ and factor sum), and less variance in the specific factor scores. Among the network-behavior analyses, in most cases, the observed model value was within the simulated attack distribution, placing the observed R2 value within a few percent of the central tendency of the simulated distribution. We did not obtain evidence that the observed brain-behavior relationships always significantly overperform or underperform estimates from simulated attacks for the measures we examined. These findings suggest that the majority of brain-behavior relationships in the simulated permutations are driven by the direct effects of the lesion on the connectome, i.e., that the simulated attacks were an informative basis to obtain confidence intervals for the effects of prototypical lesions on the connectome and behavior. Cases where specific deviations between the observed stroke and simulated attack models were found (e.g., in the FA lexical, phonological, and semantic factors) could reveal the influences of sampling effects, premorbid network organization and behavior, deafferentation, diaschisis, adaptive neuroplasticity, or related neurologic effects and secondary damage cascades (Turkeltaub, 2019). Complementary physiological data sources would be necessary to test these possibilities.
Across the simulations, there was an intuitive relationship between measures of global network topology and overall aphasia severity, accounting for 20% of the variance on the WAB-AQ. Interestingly, this network-behavior relationship was stronger than that observed with the total factor sum score. Within the left hemisphere, this pattern remained, and lexical processing and auditory comprehension tended to have strong relationships with network topology, potentially reflecting the anatomically distributed demands of these tasks in left hemispheric perisylvian circuits (Hickok and Poeppel, 2007), association regions (Dronkers et al., 2004), and sensory-perceptual pathways (Desai et al., 2010; Friederici, 2012). Perhaps because of the relatively circumscribed circuits thought to mediate phonological processing (Poldrack et al., 1999; Pollack and Ashby, 2018) and relatively preserved prefrontal circuitry that might mediate the functions in our cognitive factor (Krieger-Redwood et al., 2016; Ralph et al., 2017), we observed weaker relationships between these behaviors and the topology in the simulated attacks.
Overall, streamline networks estimated using apparent fiber density appeared to offer substantially more reliable point estimates for global and dimensional behavioral outcomes despite similar central tendencies to FA networks. This could be because our use of a reliable healthy connectome ensured consistent sets of streamline edges in healthy controls, whereas FA values are derived from estimated streamlines and offer an additional source of variance in the simulations. In addition, streamline distributions tend to be heavy-tailed with few highly connected pairs of regions (Smith et al., 2015a; Zhang et al., 2018; Yeh et al., 2019; Rosen and Halgren, 2021) with node degree (number of connections) and strength (the total weight of the connections) distributions that follow exponentially-truncated power laws (Clauset et al., 2009; Crossley et al., 2014). Lesions induce a significant loss in the number of estimated fibers, and the measured topology of these losses across subjects within the reliable healthy connectome will be strongly influenced by the heavy-tailed distribution of edges and presence of a subset of high-connection-strength nodes. In contrast, FA is computed over the estimated streamlines connecting region pairs regardless of their number, exhibiting significantly less skew and consequently relatively fewer high-connection-strength nodes. Qualitatively, FA can represent the integrity of axonal pathways (Seehaus et al., 2015; Kantarci et al., 2017), offering a distinct interpretive value relative to streamlines. However, our results suggest that increased caution when evaluating the effect sizes of brain-behavior relationships could be advised for FA relative to streamline connectomes.
Across the FA networks, it was not clear that any one of the investigated network measures uniquely corresponds to a single dimension of language function poststroke. When examined using the whole brain FA networks, edge weights were one measure that contributed to most of the behavioral measures. This finding could represent the possibility that individual variation in the myriad nearby, short-distance connections such as u-fibers and direct connections are relevant to mediating recovery (Mukherjee, 2005; Pallast et al., 2020). U-fibers dominate the brain's white matter but their links to cognition are conspicuously understudied (Y. Wang et al., 2020). Intuitively, a higher degree of intact local bypasses could facilitate adaptations to lost functions in general (Lizarazu et al., 2020). It is likely that as the behavioral measures increase in specificity, the unique edges contributing to losses in each function vary, driving differences in topology-behavior relationships (Yourganov et al., 2016). In lexical processing (where the models accounted for the most variance among the individual factor scores) edge weights were most prominently related to behavior for both streamlines and FA, suggesting that overall loss of connections independent of their relationship with lesion volume is a key mediator of deficits relative to other topological measures. More generally, damage to u-fibers could influence network local and global clustering, potentially explaining differences in observed network transitivity and corresponding changes in SWP and modularity.
Several limitations to our work can motivate future studies. We focused on a narrow set of commonly used network topology measures that characterize some of the aspects of global network organization as an initial benchmark for the connectome bases of language performance. Numerous other measures are available, but the link between specific network measures and cognitive functions remains an active area of inquiry and debate (Sporns, 2014; Medaglia et al., 2015; Bassett and Sporns, 2017; Bassett et al., 2018). In aphasia research, it remains to be seen whether valid network-based models that predict variance and inform us about the mechanisms of language function and recovery will emerge. More specific hypotheses that allow researchers to rule out spurious or nonspecific network-behavior effects, ideally informed by theoretical models, should be a focus of applied network studies. In particular, a priori theoretical language networks and data-driven networks could be used to focus analyses on one or more set of brain regions and connections that drive topological changes and could be more sensitive and specific to discrete language functions. In addition, other measures that represent connectome edges (mean or radial diffusivity, etc.; Smith et al., 2012, 2015b; Donahue et al., 2016) could be investigated. We recommend that these efforts will be best supported by collaborative efforts to pool patient samples and test the robustness of brain-behavior relationships. It would also be necessary to include larger permutation tests to distinguish the influences of all pooled covariance in the data, which would be encouraged in larger datasets if specific and sensitive differences between models are a priority. Such tests could additionally investigate fine distinctions between lesion threshold effects. Finally, we used a well-established anatomic atlas to guide our parcellation with an imputation procedure to compare healthy and stroke subjects, but numerous atlases are now available. Given that there is no consensus that a particular atlas is ideal for any specific purpose (Salehi et al., 2020), we further encourage researchers to collaboratively pool data to examine the reliability and validity of different processing decisions in anatomic and functional studies (Caiazzo et al., 2018; Sinke et al., 2018; Botvinik-Nezer et al., 2019; Boukadi et al., 2019; Elliott et al., 2020; Reid et al., 2020).
In conclusion, our simulations revealed that several anatomic connectome measures thought to be related to global network processing can be expected to account for 10–20% of the variance in language performance on clinical measures above and beyond lesion volume. Importantly, measures of whole brain and left hemisphere anatomic connectomes have stronger relationships with global language function than specific language dimensions, reflecting an intuitive relationship between network-wide integrity and overall functioning. More specific measures of anatomic circuits could be necessary to gain more sensitivity to distinct language processes. Simulated attacks are useful in leveraging matched comparison samples to obtain confidence estimates for observed effects. Differences between observed and simulated values could identify the influences of premorbid status, deafferentation, diaschisis, and neuroplasticity following stroke.
Footnotes
J.D.M. was supported by National Instituts of Health (NIH) Grants DP5-OD-021352, R01-DC-16800, and R01-AG-059763 and by the Department of the Army Grant PRMRP 12902164. P.E.T. and J.D.M. were supported by the NIH Grant R01-DC014960.
The authors declare no competing financial interests.
- Correspondence should be addressed to John D. Medaglia at johnmedaglia{at}gmail.com
