Data-driven brain network models differentiate variability across language tasks

The relationship between brain structure and function has been probed using a variety of approaches, but how the underlying structural connectivity of the human brain drives behavior is far from understood. To investigate the effect of anatomical brain organization on human task performance, we use a data-driven computational modeling approach and explore the functional effects of naturally occurring structural differences in brain networks. We construct personalized brain network models by combining anatomical connectivity estimated from diffusion spectrum imaging of individual subjects with a nonlinear model of brain dynamics. By performing computational experiments in which we measure the excitability of the global brain network and spread of synchronization following a targeted computational stimulation, we quantify how individual variation in the underlying connectivity impacts both local and global brain dynamics. We further relate the computational results to individual variability in the subjects’ performance of three language-demanding tasks both before and after transcranial magnetic stimulation to the left-inferior frontal gyrus. Our results show that task performance correlates with either local or global measures of functional activity, depending on the complexity of the task. By emphasizing differences in the underlying structural connectivity, our model serves as a powerful tool to assess individual differences in task performances, to dissociate the effect of targeted stimulation in tasks that differ in cognitive demand, and to pave the way for the development of personalized therapeutics.


Introduction
Cognitive responses and human behavior have been hypothesized to be the outcome of complex interactions between regional populations of neurons [1,2] and show significant variability across individuals. While certain patterns of brain activity are robust [3], many patterns change with learning and aging [4][5][6][7][8], and an underlying inter-subject variability in neural activity has been observed [9][10][11][12][13]. Importantly, the underlying anatomical connectivity of the brain provides a crucial backbone that drives neuronal dynamics and thus behavior [4,14,15]. Given recent and ongoing advancements in imaging techniques, such as diffusion spectrum imaging (DSI), which estimates the presence of white matter tracts connecting brain regions, mesoscale maps of anatomical brain connectivity can now be obtained [14,16]. While differences in brain connectivity have long been known to exist between diseased and healthy populations [17][18][19], recent findings indicate measurable differences in patterns of white matter connectivity across healthy individuals [20][21][22][23]. Although work is beginning to link individual variability in white matter structure, functional activity, and task performance [24][25][26][27][28][29][30], there is currently no standard methodology for evaluating the interplay between the brain's structural topology, dynamics, and function, and many open questions remain about how these features are coupled.
The new field of network neuroscience [31] provides a coherent framework in which to model and investigate this coupling. In the context of modeling human brain networks, network nodes can be chosen to represent brain regions and network edges can represent either physical connections (anatomical networks) or statistical relationships between nodal dynamics (functional networks) [17,31,32]. Analytical tools from network science have been successfully utilized to quantify both structural and functional brain network organization and to gain insight into topics as diverse as brain development [33], disease states [18, 19,34], learning [4], and intelligence [35].
Combining a network representation of the brain with computational modeling of brain dynamics allows one to further investigate the links between brain structure, dynamics, and performance by providing a controlled environment in which to perform in silico experiments and make predictions about real-world brain function. Using experimentally obtained structural brain network data combined with biologically motivated computational models of brain dynamics [36], one can build personalized brain network models of human brain activity [37][38][39]. This modeling approach has been used to gain insight into structure-function relationships in disease populations [40,41], perform virtual lesioning or resection experiments [42][43][44], and assess the differential impact of stimulation to different brain regions [45,46].
Here, we use this data-driven modeling approach to investigate the interplay between structural variation, brain activity, and task performance. Our computational model is built upon subject-specific connectomes that are combined with biologically informed Wilson-Cowan nonlinear oscillators (WCOs) [47,48] to produce simulated patterns of personalized brain activity. Additionally, this controlled computational environment allows us to perform targeted stimulation experiments-motivated by actual laboratory experiments-and quantify the emerging neural activity patterns that represent both global brain network activation and task-specific local subnetwork activation. Using this model, we construct computational measures in order to relate structure and individual performance variability across a cohort of ten healthy individuals who performed three language-demanding cognitive tasks before and after transcranial magnetic stimulation (TMS) [49] targeted at the left inferior frontal gyrus (L-IFG) of the brain. The performed tasks involved verb generation, sentence completion, and number reading, and are known to vary in their cognitive complexity [50,51].
Based upon the patterns of simulated brain activity driven by the structural networks in our model, we find that task performance can be correlated with either local or global circuitry depending on the complexity of the task. Further, we observe that, post experimentally applied TMS, the correlation between model output and task performance is weakened. Finally, we show that the eigenspectrum of the observed structural brain networks plays a key role in global brain dynamics which can additionally provide predictive insight into performance of some, but not all, tasks. Taken together, our results reveal that using personalized brain network models to simulate brain dynamics provides an important tool for studying and understanding human performance in cognitively demanding tasks and represents an important step towards the development of personalized therapeutics.

Results
In order to assess the link between variability in brain connectivity, activity patterns, and behavior, we build data-driven personalized brain network models of ten individuals who performed three language-demanding cognitive tasks before and after TMS targeted at the L-IFG (see Materials and methods). Subject-specific anatomical connectivity was derived from DSI imaging data, combined with a brain parcellation scheme based on the Lausanne atlas [52] (Fig 1a). We studied two different resolutions scales of the Lausanne atlas that parcellate the brain into either 83 or 234 regions (see Material and methods). In the main text, we present our findings using the 234-region parcellation and include the 83-region results in the Supporting Information, as the findings were largely preserved across scales. In Fig 1b, we present the structural (anatomical) network of an individual as a weighted connectivity matrix whose entries represent the strength of the connection between two brain regions. The dynamics of each brain region are modeled using nonlinear Wilson-Cowan oscillators (WCOs), coupled through a subject-specific anatomical brain network (Materials and methods). A WCO is a biologically motivated model of local brain activity, developed to describe the mean behavior of small neuronal populations [47]. The model therefore simulates a specific individual's spatiotemporal macro-level brain activity.
By holding the mathematical modeling of regional brain dynamics constant across individuals but allowing the underlying (subject-specific) structural connectivity matrix, A, to vary, the model provides a causal link between differences in the structural organization of the brain and differences in the resulting simulated brain dynamics. The global strength of the coupling between different brain regions via A is controlled by a global coupling parameter (c 5 ; Materials and methods). The dynamical state of the brain can be tuned by varying this parameter as shown in Fig 1c, which depicts the average excitatory dynamics as a function of global coupling. When the global coupling parameter exceeds a threshold value (c T 5 ), the brain dynamics abruptly transition to an active state that is characterized by the oscillatory activity of Fig 1d. Mathematically, this equates to oscillators switching from hovering near a fixed point to jumping to a limit cycle. (See [45] for more details about this transition.) Previous work has shown that the computational model is particularly sensitive to the point at which model dynamics undergo this transition: the value of c T 5 at which the transition takes places is subject-specific, and the inter-subject variability in c T 5 is greater than that seen using anatomical networks derived from different scans of the same individual [45]. The transition value c T 5 can be thought of as a proxy for the global excitability of the brain. Brains with a lower transition value require less external input (e.g. stimulation) to the system to enter the active state, whereas brains with Data-driven brain network models. (a) The brain connectivity used as a basis for the computational model is obtained by combining tractography estimates from diffusion spectrum imaging data of a specific individual's brain and a parcellation of the brain into 234 regions. (b) The resulting anatomical connectivity matrix where entries indicate the density of connections between two brain regions. (c) The dynamic state of brain activity can be tuned by the global coupling parameter, c 5 . As this parameter crosses a threshold, a sudden transition to an excited state is observed as marked by the orange arrow. (d) Representative excitatory dynamics of selected brain regions in the excited state, demonstrating a rich spectrum of temporal activity that is driven by the underlying anatomical connectivity. https://doi.org/10.1371/journal.pcbi.1006487.g001 Data-driven brain network models differentiate variability across language tasks a higher transition value require greater external input to enter the active state. Here, we use this as a parameter to measure differences in global network dynamics between individuals.

Individual variability: Model and experiments
Because we are interested in linking variability in brain structure and behavior, we first measured the extent of individual variability in anatomical connectivity, simulated brain activity, and task performance in our data set. Across our cohort of subjects, we observed measurable variability in anatomical network structure as seen in Fig 2a, which shows the standard deviation in edge weights between two brain regions, normalized by the mean edge weight. This structural variability is also manifested in the simulated brain activity depicted in Fig 2b and 2c. We observe variability across individuals in the specific patterns of brain activity in the active state (Fig 2b) as well as in the transition values (c T 5 ) (Fig 2c). Since the nonlinear WCOs are all identical in the model, these observed differences in simulated brain activity are a direct result of variation in the underlying anatomical connectivity. In Fig 2d we show the spatial map of the average regional variability in structural connectivity and functional activity across Assessing variability in structure, dynamics, and performance. (a) Variability in the connection strengths between pairs of brain regions across ten subjects, measured as the standard deviation of edge weights between brain regions across individuals, normalized by the mean of edge weights. Although there are few regions of low variability (blue), there are a significant number of regions with moderate to high variability (green to yellow). (b) Spatial patterns of excitatory activity vary across individuals due to regional variations in the structural connectivity between subjects. Each column represents the temporal average of the excitatory activity across brain regions for a given individual in the excited state and is unique in terms of the overall activity pattern. (c) The global coupling transition value (c T 5 ) also varies across our cohort of ten individuals. (d) Spatial mapping of the average variability in structural connectivity and functional activity across individuals. (e) Cognitive task performance for the ten subjects was assessed by experimentally measuring the response times for three language-demanding tasks: verb generation (VG), sentence completion (SC), and number reading (NR), both before (filled symbols) and after (open symbols) a targeted transcranial magnetic stimulation (TMS) to the left-inferior frontal gyrus (L-IFG) of the brain. individuals. Interestingly, we observe a diverse range of regional variation in functional dynamics that does not necessarily match the pattern or level of structural variation.
Finally, we assessed the extent of variability in the cognitive performances of individuals across three language-demanding tasks: (i) verb-generation (VG); (ii) sentence-completion (SC); and (iii) number-reading (NR). Performance was measured as the median response time across multiple (~50) trials (see Materials and methods) and shows variability across both subjects and tasks (Fig 2e). Moreover, performances were altered after transcranial magnetic stimulation (TMS) was applied to the left inferior frontal gyrus (L-IFG). The L-IFG, also known as Broca's area, is traditionally believed to play an important role in language comprehension and syntactic processing, and specifically in the selection and retrieval of words [53]. It is therefore expected that external stimulation to this region should affect task performance. However, it is important to note that while subject performance on a task did often change after TMS, we did not consistently observe an improvement (or degradation) in task performance across individuals, nor was the effect consistent between tasks within a given subject. This suggests that although the L-IFG plays an important role in the context of language comprehension, the actual cognitive response reflects contributions from a larger part of the brain network.
Given the observed variability in the structural, dynamical, and behavioral aspects of our data, we next focused on assessing how this variability was related across these three domains. We therefore examined network features measured at both the global network level (using the entire brain network) and within task-specfic subnetworks that were selected to represent the specific circuitry involved in task completion and asked how these measures related to task performance.

Global network features correlate with certain task performances
We first examined the relationship between global network properties and task performance by estimating the correlation between global brain activation and task performance. For each subject, we measured the threshold value of the global coupling parameter (c T 5 ) at which the individual's brain transitions to the excited state, and we calculated the correlation between this value and their performance on each task (before the application of TMS). As seen in Fig 3 and Table 1, we observed a significant positive correlation (r = 0.86, 90% CI [0.68 0.95], p = 0.001) between model transition values and task performance in the sentence completion (SC) task. Thus for the SC task, individuals with a lower value of c T 5 (more easily excitable brain) are likely to perform better (as measured through a short response time). However, we did not observe a significant correlation in the verb generation (VG) or number reading (NR) tasks, indicating that the performance of these tasks cannot be predicted by a global network property. To ensure that these results were dependent upon the organization of the subjectspecific anatomical connectivity used as a basis for the computational model, for each individual, we created randomized brain networks by preserving the distribution of edge weights but randomly reassigning connection strengths between brain regions (Materials and methods). We recalculated the c T 5 values for simulations using these randomized connectivity matrices, but did not observe any significant correlations between transition values and task performance in this case.
To further explore the link between global brain dynamics and behavior, we additionally assessed the relationship between specific patterns of brain activity and task performance. Since the L-IFG is involved in controlled language processing [54][55][56], one can argue that the pattern of brain synchronization as a result of targeted stimulation to the region might also be predictive of task performance. We therefore computationally stimulated the brain regions comprising the L-IFG (Fig 4a), and quantified how the stimulation spread throughout the global brain network (see Materials and methods for details). As shown in Fig 4a, computational stimulation pushes the dynamics of the region into oscillatory activity which then drives the functional dynamics of other brain regions through the underlying structural connections. We measure the resulting pattern of brain activity by calculating the pairwise functional connectivity using the maximum normalized correlation between brain regions [45,57,58] (Materials and methods). In Fig 4b, we show the spatial mapping of the variability in average functional connectivity across subjects (measured as the standard deviation divided by the mean) resulting from computational stimulation of the L-IFG, and in Fig 4c, we show the subject-specific patterns of functional connectivity for three subjects. Due to the variability in the underlying anatomical connectivity matrices, the resulting patterns of functional connectivity differ between individuals. We measure the extent of the global spread of synchronization by calculating the functional effect [45] which measures the average value of synchronization across the entire brain (Materials and methods). Interestingly, unlike our observation with the The variables r and p denote the Pearson correlation coefficient and associated p-value, respectively. The 90% confidence interval for r is reported below each correlation. Here a Ã denotes that the observed correlation is significant under FDR correction for multiple comparisons across tasks (for p < 0.05) and a • denotes a significant correlation across the two scales of the brain parcellation studied in this paper (S1 Table). VG = verb generation, SC = sentence completion, and NR = number reading. global coupling parameter, the global functional effect does not show a significant correlation with task performance for any of the cognitive tasks (Table 1).

Localized task-specific circuit features correlate with other task performance
While we did observe a significant correlation between the global threshold value and task performance for the SC task, we saw no correlations between global brain dynamics and task performance in the remaining two tasks, and the global functional effect was not correlated with performance in any of the tasks. However, the three language tasks performed in this study differ in semantic demands, and the absence of a significant correlation for the VG and NR tasks could be due to either a drastically different cognitive mechanism for performing these tasks, or the dependence of these tasks on a more localized brain circuit. To investigate the latter possibility, we assessed the role of task-specific subnetworks in task performance.
To construct a task dependent, spatially localized measure, we follow the work of Roux et al. [59] to identify the possible brain regions involved in reading alphabets and numbers. Roux et al. used intracranial direct cortical stimulation paired with behavioral measurements to spatially map the brain regions that were differentially involved in reading both alphabets and numbers. This approach provides a much more direct and causal form of evidence that implicates the role of specific brain regions in the behaviors of interest as compared to anatomical or functional networks based on task-related activity measured using fMRI, which are Computational stimulation of a single brain region pushes its dynamics into a limit cycle (oscillations). (b) Spatial mapping of the variability in functional connectivity across subjects after computational stimulation to the L-IFG. Variability is measured as the standard deviation in regional functional connectivity across subjects divided by mean. (c) Individual functional connectivity matrices for three subjects resulting from computational stimulation of the L-IFG. Note the variation in the observed connectivity patterns across subjects. https://doi.org/10.1371/journal.pcbi.1006487.g004 Data-driven brain network models differentiate variability across language tasks indirect measures of cognitive circuits. We mapped these regions to the Lausanne atlas and constructed two task circuits: one involved in VG and SC (alphabets-related, Fig 5a), and one involved in NR (number-related, Fig 5b). We found that these circuits are also consistent with other studies mapping brain regions involved in language processing [53,[60][61][62]. Note that both of these sub-networks are contained entirely in the left hemispheric language network.
We then calculated the functional effect within these task circuits (averaging the functional connectivity values only within the subnetwork as opposed to the entire brain network as done previously) and correlated this local measure with task performances (Fig 5c-5e and Table 1). We observed a significant positive correlation between the task-specific functional effect and task performance for the NR task (r = 0.74, 90% CI [0.20 0.94], p = 0.016), indicating that performance on the NR task depends on the localized spread of activation throughout the task sub-network and is less dependent on the global brain network structure. Our results indicate that individuals with a lower functional effect (less synchronization within the task circuit) also have a lower response time (better performance). When synchronization within the task circuit is increased (a high functional effect), task performance degrades, suggesting that high levels of synchronized activity within the task circuit could potentially impede the ability of the circuit to perform localized computations necessary for task completion.  Table 1. There is a significant correlation between the functional effect and task performance for the NR task (r = 0.74, 90% CI [0.20 0.94], p = 0.016) that is preserved across the different scales of brain parcellation (S1 Fig and S1 Table). https://doi.org/10.1371/journal.pcbi.1006487.g005 Data-driven brain network models differentiate variability across language tasks If we compare the functional effect measured only within brain regions outside of the task circuit, the significant correlation with the NR task is lost (Table 1), indicating the specificity of the task circuit in the model. Although we also observed a significant correlation between the task-specific functional effect and task performance for the SC task, this result is driven by a single subject and does not hold if this subject is removed from the data set. Performing the analysis on a larger data set would therefore be necessary to confirm this finding for the SC task. No significant correlations were observed for the VG task.
To validate the specificity of our selected task circuits, we constructed 10,000 random subnetworks by randomly selecting the same number of brain regions as in each task circuit and then calculated the functional effect within these random circuits after stimulating the L-IFG (Materials and methods). The randomized circuits had a variable degree of overlap with the actual task circuit (0 to 35%). We observed that only 1.8% of the randomly selected circuits gave significant correlations (r > 0.5 and p < 0.05) between the circuit-specific functional effect and task performance, which we estimate to be the false positive rate in our computational predictions. This low error rate signifies that the observed significant correlation in the NR task is due to the selection of brain regions in the task-specific circuit.
We also verified that the observed effects were related to our choice of stimulating the L-IFG as opposed to some other brain structure. We chose different sets of brain regions (equal in size to the number of brain regions that compose the L-IFG) that were randomly distributed within the task circuit and applied targeted computational stimulation to these randomly selected regions. When randomly selected brain regions were stimulated, we did not observe a significant correlation between the task performance and functional effect, confirming the importance of specifically targeting the L-IFG in our in silico experiments. Additionally, consistent with previous findings [45], we observed that the patterning of activation as a result of stimulation of randomly selected brain regions varied with the selection of brain regions and across subjects. These findings support the possibility that in the future our modeling approach can be used to help design and optimize therapeutic strategies that use external activation such as TMS to treat neurological disorders [63].

Correlating task performance post-TMS
We also asked if our computational model correlated with individual performance after the application of experimentally applied TMS targeted at the L-IFG. The underlying mechanisms of how TMS affects the brain are not well understood, but it is believed that TMS locally influences neuronal firing which can then propagate within the brain through inter-regional neuroanatomical pathways [45,64]. We therefore examined the correlation between model features and behavioral performance during the post-TMS task. As seen in Fig 6 and Table 2, while we still observe a positive correlation between the transition value and performance in the SC task (r = 0.68, 90% CI [0.32 0.87], p = 0.03) and between the functional effect within the task circuit and performance in the NR task (r = 0.63, 90% CI [-0.07 0.88], p = 0.05), in both cases, the strength of the correlation is decreased when compared to correlations with task performance before TMS. It should also be noted that these correlations exist across the two spatial scales of parcellation but their significance does not survive multiple comparison corrections. Speculatively, the fact that the correlation strength is weakened post-TMS potentially reflects contributions of noise to the system mediated through inhibitory stimulation to the L-IFG.

Correlations with graph theoretical measures of network connectivity
Finally, to be able to draw a direct connection between anatomical brain connectivity and behavior, we asked if the correlations observed in our computational model of brain dynamics  Table 2. There is a significant correlation between the global transition value and task performance for the SC task  The variables r and p denote the Pearson correlation coefficient and associated p-value, respectively. The 90% confidence interval for r is reported below each correlation. Here a • denotes a significant correlation across the two scales of the brain parcellation studied in this paper (S2 Table). VG = verb generation, SC = sentence completion, and NR = number reading. https://doi.org/10.1371/journal.pcbi.1006487.t002 Data-driven brain network models differentiate variability across language tasks could be revealed using only graph theoretical measures of network structure applied directly to the anatomical connectivity matrices alone (in the absence of model dynamics). We therefore calculated three measures of network structure that have been shown to be related to the spread of synchronization in networks: the average degree, inverse spectral radius, and synchronizability (see Materials and methods). The average degree measures the average strength of connections per brain region, while the spectral radius and synchronizability are related to the eigenspectrum of the adjacency matrix or graph Laplacian respectively, and relate to the ease and extent to which the network is expected to synchronize under the assumption that regions behaved as oscillators.
In Table 3, we report the Pearson correlation coefficients for these measures of network structure with task performance both before and after TMS was applied. Both the average degree and the inverse spectral radius show significant correlations with task performance for the SC task, and as observed earlier, show a decrease in their correlation magnitude after the application of TMS. Both of these measures are known to relate to the ease of global synchronization in a network [65,66], and therefore this finding is somewhat expected given our prior observations about c T 5 . Interestingly, we find that the synchronizability of the network does not correlate with performance of the SC task. Instead, this measure, which has also been shown to be more sensitive to changes in local network structure [67,68], shows a correlation with performance of the NR task that increases after the application of TMS. However, this correlation is weaker than that observed with the local measurements of the functional effect within the task circuit reported earlier in Table 1.

Discussion
We presented a computational model of individualized brain dynamics that allows us to make predictions about cognitive task performance based on the variability of anatomical brain connectivity between people. By analyzing global model excitability and localized patterns in the spread of a targeted stimulation in our data-driven model, we explain individual variability in The variables r and p denote the Pearson correlation coefficient and associated p-value, respectively. The 90% confidence interval for r is reported below each correlation. Here a Ã denotes that the observed correlation is significant under FDR correction for multiple comparisons across tasks (for p < 0.05) and a • denotes a significant correlation across the two scales of the brain parcellation studied in this paper (S3 Table). VG = verb generation, SC = sentence completion, and NR = number reading. https://doi.org/10.1371/journal.pcbi.1006487.t003 Data-driven brain network models differentiate variability across language tasks the performance of cognitively demanding tasks and establish the significant role of brain network organization in driving individual behavior. We examined three cognitive tasks and found that different measures of model dynamics gave rise to correlations with task performance. Our results indicate that as the complexity of the cognitive task increases, a larger portion of the brain network, including interhemispheric connections, contribute in determining the overall response. The SC task, which has additional language processing demands compared to the other two tasks, requires understanding different words, constructing the overall meaning of the sentence, and determining the appropriate word to fit in the sentence. For this cognitively demanding task, a global parameter of the model, c T 5 , which could be related to the overall excitability of the network, correlates with individual task performance. We also observed significant correlations between global network properties of the anatomical connectivity matrices (specifically the average degree and inverse spectral radius) and performance of the SC task. The fact that these global measures are correlated with performance in the SC task but not the other two tasks possibly reflects the global nature and complexity of the network of brain regions required to complete a sentence [69].
On the other hand, performances on the NR task, which has been established to be more localized in terms of the brain regions involved (i.e., involves the activation of only a sub-network of the brain as opposed to the entire brain network) [59], can be predicted by tracking the spread of a targeted stimulation within a task-specific circuit in the brain. Further, we observed that a higher functional effect of the stimulation correlates with a higher response time (poorer performance) on the NR task. In our model, a high functional effect indicates a high degree of synchronization within the task circuit. This heightened level of synchronization could effectively constrain the degrees of freedom for brain circuit activity and computation, resulting in poorer task performance. Better task performance might be achieved through an integration and segregation mechanism which allows brain regions to mutually communicate but does not require synchrony over a longer time scale. Such a modular functional organization has been associated with efficient learning of motor skills [5].
Interestingly, no tested measures of model dynamics or structural network features were able to explain the performance of the VG task. It is likely that VG requires a different coupling mechanism of global and local dynamics which is not captured in the model features constructed here. For example, while verb generation and sentence completion are both "open ended" language tasks, SC requires our ability to accumulate the probability of a response over an entire sentence, whereas VG is cued by a single word. Thus, even if both involve the L-IFG, SC might recruit quite a few more distributed resources, while VG might require more focal resources. It is possible that our selection of regions comprising the language task circuit was not sensitive enough to these focal areas. We propose that this problem could potentially be resolved by a more extensive research effort that combines experimental and theoretical approaches from cognitive and network neuroscience in order to better understand the specific circuitry involved in task completion.
A somewhat surprising finding was that the application of TMS did not have a directional effect; task performance both improved and degraded post-TMS. One might ask if this could be interpreted as a regression to the mean. However, here, our primary focus for behavioral analyses was on the correlation between inter-subject differences in response times and the simulated dynamic measures. We found evidence that our measures of simulated network excitability independently correlated with behavior before and after stimulation. That is, the relationship between behavior and our excitability measures cannot be explained by regression to the mean because they do not involve a correlation with changes in behavior over time, which would indeed be potentially susceptible to regression to the mean. However, it is true that predicting variability in change scores where a brain stimulation effect is strong across subjects is an important future direction.
A strength of our modeling approach is that the construction of localized task circuits is not limited to those involved in the tasks studied here, and could easily be extended to represent different brain circuitry. Ultimately, our model magnifies the effects of differences in the observed anatomical brain structure, and it is not always feasible to measure differences in network structure in terms of general network statistics. While some network statistics also provided predictive power of task performance, the strength of the correlation was only comparable for measures of global network structure, and our computational model can generate measures of structure and dynamics specific to the task that is being performed. We also note that while in the SC task we saw correlations between global measures of network structure such as the average degree and inverse spectral radius, we did not observe consistently significant correlations between measures of network structure and the NR task; correlations between synchronizability and NR task performance existed only at the 234-region scale and did not survive multiple comparison corrections. This suggests that for tasks that invoke local computations involving specific circuitry, generic measures of network structure might not be sensitive to anatomical variations in the local circuitry. However, within the modeling framework presented here, we are able to define task-specific subnetworks and examine the patterns of brain activity within these local circuits, making the use of computational modeling preferable to assess potential relationships between variability in individual brain structure and task performance.
It should also be noted that the functional measures examined here (for both global and task circuit specific calculations) involved measuring the spread of synchronization after computational stimulation of the L-IFG. The choice to stimulate the L-IFG represents essential a priori knowledge from cognitive neuroscience studies, and stimulating randomly chosen sets of brain regions within the task circuits did not produce significant correlations with task performances. In other paradigms, one might prefer to study the effects of stimulation to a different brain region, or to compare the effects of applying stimulation to multiple different brain regions. In each scenario, the exact design of the appropriate computational experiment will depend on the cognitive question at hand and the type of experimental data with which one wishes to compare model outputs. The ability to integrate knowledge from cognitive neuroscience into computational experiments represents a strength of our data-driven modeling approach, as it makes the model flexible and opens the door to studies across a range of cognitive paradigms. As stimulation strategies are increasingly being used to treat neurological disorders, we hope that our modeling framework can guide the design, personalization, and understanding of such treatment strategies.
Although the current study was constrained by its small subject population size, it is encouraging that despite this limitation, we were still able to observe significant correlations between model features and certain task performance. These findings should therefore serve as a proof of principle study to promote the use of similar data-driven modeling approaches in larger data sets with more subjects and a higher diversity of task conditions. The use of personalized brain network models will serve as an increasingly valuable tool to establish explicit links between brain connectivity, dynamics, and behavior, and to develop personalized therapeutic strategies.

Ethics statement
All procedures were approved in a convened review by the University of Pennsylvania's Institutional Review Board and were carried out in accordance with the guidelines of the Institutional Review Board/Human Subjects Committee, University of Pennsylvania. All participants volunteered with informed consent in writing prior to data collection.

Subjects and cognitive tasks
Ten healthy individuals (mean age = 25.4, St.D. = 4.5, 6 female) from a larger neuroimaging study [70] returned to participate in the present study. Participants performed two openended language tasks and one closed-ended number naming task. The language tasks included a verb generation task [50] and a sentence completion task [51]. For the verb generation task, subjects were instructed to generate the first verb that came to mind when presented with a noun stimulus (e.g., 'cat'). The verb could be either something that the noun does or something that can be done with the noun. Response times (RTs) were collected from the onset of the noun cue to the onset of the verb response. For the sentence completion task, subjects were presented with a sentence, for example "They left the dirty dishes in the --?", and were instructed to generate a single word that appropriately completes the sentence, such as 'sink'. Words in the sentences were presented serially in 1 s segments consisting of one or two words. RTs were computed as the latency between the onset of the last segment, which always contained two words (i.e., a word and an underline), and the onset of the participant's response. For all items in the sentence completion task, items in the high vs. low selection demand conditions were matched on retrieval demands (association strength) [51]. For both language tasks, each trial began with the presentation of a fixation point (+) for 500 ms, followed by the presentation of the target stimulus, which remained on the screen for 10 s until the subject made a response. Subjects were given an example and five practice trials in the first administration of each language task (i.e., before TMS), and were reminded of the instructions before performing the task a second time (i.e., after TMS). In each of the before and after TMS conditions, subjects completed 50 trials for a total of 100 trials.
In the number reading task, participants produced the English names for strings of Arabic numerals presented on the screen. On each trial, a randomized number (from tens of thousands to millions; e.g., 56395, 614592, 7246856) was presented in black text on a white background. The numbers were uniformly distributed over three lengths (17 per length for each task administration). The position of items on the screen was randomized between the center, left, and right of the screen to reduce the availability of visual cues to number length and syntax [50]. RTs were collected from the onset of the stimulus presentations to the onset of the subject's response. The number appeared in gray following the detection of a response (i.e., voice key trigger), and remained on the screen thereafter to reduce the working memory demands required for remembering the digit string. At the start of the experiment, subjects performed 50 trials of the number naming task to account for initial learning effects [50]. Prior to performing the task for the first time, subjects were given an example and five practice trials, and were later reminded of the instructions before performing the task a second (i.e., before TMS) and a third (i.e., after TMS) time. In each of the before and after TMS conditions, subjects completed 51 trials for a total of 102 experimental trials. The items for the verb generation task were identical to those used in [71] and the items for the sentence completion task were those from [72]. The difficulty of items was sampled to cover a distribution of values computed via latent semantic analysis (LSA) applied to corpus data.
Verbal responses for all tasks were collected from a computer headset microphone. The microphone was calibrated to reduce sensitivity to environment background noise prior to the collection of data for each session such that the recording software was not triggered without clear verbalizations. List order (before or after TMS) was counterbalanced across participants.
Item presentation order within each task was fully randomized across participants. Task performance was assessed based on the subject's median response time across all the trials.

Experimental design
A schematic of the experimental design can be seen in Fig 7. Subjects performed the verb generation, sentence completion, and number reading tasks in a randomized order across subjects. Then, the subjects received 40s of continuous theta burst stimulation. Finally, the subjects performed alternate difficulty-matched forms of each task (i.e., matched distributions of association strengths and entropy values across items for the sentence completion and verb generation task, and matched distributions of number lengths in the number reading task) after stimulation in the same order as before stimulation.

Transcranial magnetic stimulation
The Brainsight system (Rogue Research, Montreal) was used to co-register MRI data with the location of the subject and the TMS coil. The stimulation site was defined as the posterior extent of the pars triangularis in each individual subject's registered T1 image. A Magstim Super Rapid2 Plus1 stimulator (Magstim; Whitland, UK) was used to deliver continuous theta burst stimulation (cTBS) via a 70 mm diameter figure-eight coil. cTBS was delivered at 80% of each participant's active motor threshold [73]. Each participant's threshold was determined prior to the start of the experimental session using a standard up-down staircase procedure with stimulation to the motor cortex [M1]. The measured thresholds (given in units measured as the percentage of machine output from the Magstim device) for the 10 participants were 68,60,46,57,37,36,44,39,46, and 65, respectively.

Human DSI data acquisition and preprocessing
Diffusion spectrum images (DSI) were acquired for a total of 10 subjects along with a T1-weighted anatomical scan at each scanning session, in line with previous work [64]. DSI scans sampled 257 directions using a Q5 half-shell acquisition scheme with a maximum bvalue of 5,000 and an isotropic voxel size of 2.4 mm. We utilized an axial acquisition with the following parameters: repetition time (TR) = 5 s, echo time (TE) = 138 ms, 52 slices, field of view (FoV) (231, 231, 125 mm). DSI data were reconstructed in DTI Studio (www.dsi-studio. labsolver.org) using q-space diffeomorphic reconstruction (QSDR) [74]. QSDR first reconstructs diffusion-weighted images in native space and computes the quantitative anisotropy (QA) in each voxel. These QA values are used to warp the brain to a template QA volume in Montreal Neurological Institute (MNI) space using the statistical parametric mapping (SPM) nonlinear registration algorithm. Once in MNI space, spin density functions were again reconstructed with a mean diffusion distance of 1.25 mm using three fiber orientations per voxel. Fiber tracking was performed in DSI Studio with an angular cutoff of 35˚, step size of 1.0 mm, minimum length of 10 mm, spin density function smoothing of 0.0, maximum length of 400 mm and a QA threshold determined by DWI signal in the colony-stimulating factor. Data-driven brain network models differentiate variability across language tasks Deterministic fiber tracking using a modified FACT algorithm was performed until 1,000,000 streamlines were reconstructed for each individual.
Anatomical (T1) scans were segmented using FreeSurfer [75] and parcellated using the connectome mapping toolkit [76]. We studied two different scales of regional parcellation (n = 83 or n = 234 brain regions) which were registered to the B0 volume from each subject's DSI data. The B0 to MNI voxel mapping produced via QSDR was used to map region labels from native space to MNI coordinates. To extend region labels through the grey-white matter interface, the atlas was dilated by 4 mm [77]. Dilation was accomplished by filling non-labelled voxels with the statistical mode of their neighbors' labels. In the event of a tie, one of the modes was arbitrarily selected. Each streamline was labeled according to its terminal region pair.

Construction of anatomical brain networks
To construct the subject-specific anatomical connectivity networks each brain region was represented as a network node. A list of brain regions included in our study is given in S4 Table  for both the 83-region and 234-region parcellations. As in prior studies, we define pairwise connection weights between nodes based on the number of streamlines connecting brain regions and normalized by the sum of the volumes of the nodes [45,52]. This procedure results in a sparse, weighted, undirected structural brain network for each subject (N = 10), where network connections represent the density of white matter tracts between brain regions (Fig 1b).

Computational model of brain dynamics
In our data-driven network model, regional brain dynamics are given by Wilson-Cowan oscillators [45,47]. In this biologically motivated model of neuronal populations, the fraction of excitatory and inhibitory neurons active at time t in the i th brain region are denoted by E i (t) and I i (t) respectively, and their temporal dynamics are given by: We note that A ij is an element of the subject-specific coupling matrix, A, whose value is the connection strength between brain regions i and j as determined from DSI data (see above). The global strength of coupling between brain regions is tuned by excitatory and inhibitory coupling parameters c 5 and c 6 respectively. We fix c 6 = c 5 /4, representing the approximate ratio of excitatory to inhibitory coupling. P i (t) represents the external inputs to excitatory state activity and is used to perform computational stimulation experiments (see below). The parameter t ij d represents the communication delay between regions i and j. If the spatial distance between regions i and j is d ij , t ij d ¼ d ij =t d , where t d = 10m/s is the signal transmission velocity. Additive noise is input to the system through the parameters w i (t) and v i (t) which are derived from a normal distribution and σ = 10 −5 . Other constants in the model are biologically derived: c 1 = 16, c 2 = 12, c 3 = 15, c 4 = 3, a E = 1.3, a I = 2, θ E = 4, θ I = 3.7, τ = 8 as described in references [45,47].
To numerically simulate the dynamics of the system we use a second order Runge Kutta method with step size 0.1 and initial conditions E i (0) = I i (0) = 0.1. All analysis is performed after allowing the system to stabilize for 1s.

Assessing individual variability
In order to calculate the model transition values for each individual, we ran 1s simulations (after allowing the system to stabilize) for a range of c 5 parameters (0.05 c 5 0.25, with a step-size of 0.001) in which no external input was applied (P = 0). The average excitatory activity was recorded for each region as a function of c 5 as shown in Fig 1c. The value of c 5 at which we observed a sudden increase in the average activity (marked by an arrow in Fig 1c), was identified as the transition value, c T 5 . Structural variability was measured by calculating the standard deviation of a given connection strength A ij between brain regions i and j for all the subjects and then normalizing by the average connection strength, <A ij >. In Fig 2d we show the regional average value of the structural variability across individuals. Similarly, the average functional variability, as also shown in Fig 2d, was estimated for brain regions by calculating the standard deviation of the excitatory activity of a given brain region E i across individuals and then normalizing by the average excitatory activity of that region across individuals.
In Fig 4b, the variability in regional functional connectivity is calculated as the standard deviation of the average regional functional connectivity across individuals divided by the regional mean of functional connectivity across individuals.

Targeted computational activation
To activate (stimulate) a particular brain region, we applied a constant external input P i = 1.15 which drives the regional activity into a limit cycle as shown in Fig 4b. Before targeted activation, the global coupling parameter c 5 was fixed just below its transition value c T 5 (which differs between subjects) to perturb the system in a state with maximum sensitivity for perturbations. When stimulating the L-IFG, we simultaneously activated four brain regions (pars orbitalis (POB), pars triangularis (PTA) and pars opercularis (POC, two regions), Fig 4a) by applying an external input of P = 1.15.

Functional connectivity and functional effect of stimulation
To quantify the spread of computational stimulation, functional connectivity was determined by calculating the pairwise maximum normalized cross-correlation [57,58] between the excitatory activity E i (t) and E j (t), for brain regions i and j. We used a window size of 1s and calculated correlations over a maximum lag of 250ms.
In order to quantify the functional effect of stimulation, we considered the regional dynamics within a window of 2s once the system is stabilized after the initial transient period. The system was first allowed to evolve without any external input (P i = 0) for 1s and then an external input of P i = 1.15 was applied for 1s to the set of brain regions selected for activation or stimulation. Functional connectivity was calculated separately for the stimulation-free period and the period of stimulation. We then calculated the difference between the pairwise values of functional connectivity measured during and before stimulation, resulting in a matrix that describes the pairwise changes in functional connectivity resulting from the stimulation. The functional effect of stimulation [45] was then calculated using this matrix either globally by averaging over the entire matrix, or within a task circuit by averaging within the task circuit, or outside of a task circuit by averaging over regions outside of the task circuit.

Defining task circuits
To construct a task dependent localized measure, we followed the work of Roux et al. [59] to identify the possible brain regions involved in reading alphabets and numbers. We extended their findings to propose that the brain regions involved in reading alphabets contribute towards performing VG and SC tasks, and the brain regions involved in reading numbers contribute towards performing the NR task. We mapped these regions to the Lausanne atlas and constructed possible task circuits involved in VG and SC (alphabets-related , Fig 5a), and NR (number-related, Fig 5b). These circuits are also consistent with other studies mapping brain regions involved in language processing [53,[60][61][62].

Randomizing network structure and stimulation effect
To assess the effect of global network organization, we repeated our entire computational analysis with randomized anatomical connectivity matrices. Anatomical connectivity matrices were randomized separately for each individual in order to alter the original brain network organization. Here, we preserved the overall connectivity distribution but randomly reassigned connection strength values to each pair of the brain regions from this distribution.
To test the specificity of our results to the choice of regions comprising the task circuit, we constructed 10, 000 sub-networks by randomly picking the same number of brain regions as in the NR circuit (22) from all of the 234 brain regions. We then stimulated the L-IFG and calculated the functional effect within these randomly constructed sub-networks for each individual. These sub-networks had varying degree of overlap with the actual task circuit, ranging from 0 to 35%. We observed that only 1.8% of the random sub-networks produced a significant correlation between the localized functional effect and task performance with r > 0.5 and p < 0.05, indicating that our results are indeed related to the specific construction of the task circuit.
In order to assess the importance of stimulating the L-IFG for the cognitive tasks chosen in this particular work, we compared our findings with those when stimulation was applied to a group of regions randomly chosen (under spatial constraints) within the actual task circuits with a regional brain volume equivalent to L-IFG (4 brain regions). The 4 brain regions were chosen such that they formed a continuous spatial volume within the brain. We could identify 7 such volumes within the task circuits that we constructed, and used these 7 alternate stimulation sites to assess the specificity of our findings to the choice of stimulation site. The model did not produce significant correlations when any of these 7 volumes were stimulated, indicating that the model is indeed sensitive to the choice of stimulation site.

Confidence intervals and multiple comparison testing
In order to calculate 90% confidence intervals of the Pearson's correlation coefficient, r, we bootstrapped each sample 5000 times (with replacements). We used the resulting distribution of r values to obtain the 0.05 and 0.95 quantiles.
We tested the significance of correlations under multiple comparison corrections by applying the false discovery rate (FDR) correction to p-values across tasks. Correlations were deemed to be significant under multiple comparisons if the corrected p-value was less than 0.05. In all tables throughout the text, we report original (uncorrected) p-values, and indicate whether or not the findings remained significant as a result of applying the FDR correction.

Network statistics
We calculated network statistics for each subject using the structural matrices derived from their DSI data. The weighted degree centrality, k i , for a given region i is defined as [78]. The average across the degree centralities of all network nodes was used to obtain the average degree of a given subject. The spectral radius, S r , is given by the maximum eigenvalue of the connectivity matrix A, S r = λ| Max , and is a global measure of network structure that is related to the spread of synchronization in a network [65,66]. Synchronizability, S, is defined as, S ¼  Fig 3. (d-e) Correlations between the functional effect within the task-specific circuits and task performance for comparison with the 234-region scale shown in Fig 5. The red lines represent a linear fit to the data for visual guidance. The corresponding Pearson correlation coefficients are given in S1 Table. There is a significant correlation between the transition value and task performance for the SC task (panel b), and functional effect and task performance for the NR task (panel e). Model features are calculated for the 83-region scale of parcellation. (a-c) Correlations between the global transition value and task performance, and (d-e) correlations between the functional effect within the task-specific circuits and task performance for comparison with the 234-region scale shown in Fig 6. The red lines represent a linear fit to the data for visual guidance. The corresponding Pearson correlation coefficients are given in S2 Table. We observe somewhat significant correlations between the transition value and task performance for the SC task, and functional effect and task performance for the NR task. (EPS) S1 Table. Correlations between model features and task performance for the 83-region parcellation. The variables r and p denote the Pearson correlation coefficient and associated p-value, respectively. The 90% confidence interval for r is reported below each correlation. Here a Ã denotes that the observed correlation is significant under FDR correction for multiple comparisons across tasks (for p < 0.05) and a • denotes a significant correlation across the two scales of the brain parcellation studied in this paper. VG = verb generation, SC = sentence completion, and NR = number reading. (DOCX) S2 Table. Post-TMS correlations between model features and task performance for the 83-region parcellation. The variables r and p denote the Pearson correlation coefficient and associated p-value, respectively. The 90% confidence interval for r is reported below each correlation. Here a Ã denotes that the observed correlation is significant under FDR correction for multiple comparisons across tasks (for p < 0.05) and a • denotes significant correlation across the two scales of the brain parcellation studied in this paper. VG = verb generation, SC = sentence completion, and NR = number reading. (DOCX) S3 Table. Correlations between anatomical network features and task performance for the 83-region parcellation. The variables r and p denote the Pearson correlation coefficient and associated p-value, respectively. The 90% confidence interval for r is reported below each correlation. Here a Ã denotes that the observed correlation is significant under FDR correction for multiple comparisons across tasks (for p < 0.05) and a • denotes a significant correlation across the two scales of the brain parcellation studied in this paper. VG = verb generation, SC = sentence completion, and NR = number reading.