Functional Mechanisms of Recovery after Chronic Stroke: Modeling with the Virtual Brain123

Abstract We have seen important strides in our understanding of mechanisms underlying stroke recovery, yet effective translational links between basic and applied sciences, as well as from big data to individualized therapies, are needed to truly develop a cure for stroke. We present such an approach using The Virtual Brain (TVB), a neuroinformatics platform that uses empirical neuroimaging data to create dynamic models of an individual’s human brain; specifically, we simulate fMRI signals by modeling parameters associated with brain dynamics after stroke. In 20 individuals with stroke and 11 controls, we obtained rest fMRI, T1w, and diffusion tensor imaging (DTI) data. Motor performance was assessed pre-therapy, post-therapy, and 6–12 months post-therapy. Based on individual structural connectomes derived from DTI, the following steps were performed in the TVB platform: (1) optimization of local and global parameters (conduction velocity, global coupling); (2) simulation of BOLD signal using optimized parameter values; (3) validation of simulated time series by comparing frequency, amplitude, and phase of the simulated signal with empirical time series; and (4) multivariate linear regression of model parameters with clinical phenotype. Compared with controls, individuals with stroke demonstrated a consistent reduction in conduction velocity, increased local dynamics, and reduced local inhibitory coupling. A negative relationship between local excitation and motor recovery, and a positive correlation between local dynamics and motor recovery were seen. TVB reveals a disrupted post-stroke system favoring excitation-over-inhibition and local-over-global dynamics, consistent with existing mammal literature on stroke mechanisms. Our results point to the potential of TVB to determine individualized biomarkers of stroke recovery.


Introduction
Previous research has provided key insights into the disease process in stroke. Studies in mammals have uncovered basic mechanisms of ischemic injury, inflammatory responses, and cellular recovery (Carmichael, 2012;Nudo, 2013). In humans, researchers have suggested predictive imaging biomarkers for disease progression and recovery, mapped associated changes in brain networks, and developed new rehabilitative therapies (Reiss et al., 2012). Despite this, stroke remains a major source of disability in the United States, with ϳ6.5 million people living with stroke, with some level of hemiparesis present in ϳ50% (Go et al., 2014). This is neither the fault of mammal nor human studies, as both are constrained by their respective study populations. Studies in mammals are well controlled yet homogeneous, limiting their translational abilities. Human studies reflect the population at hand, yet often rely on indirect measures, obscuring the full picture. Although both share a common goal of curing stroke via the repair and reorganization of the injured brain, what is missing is a translational bridge to effectively span the divide between basic mechanisms and dynamic human brain systems.
At the same time, the neuroscience community is immersed in collecting large datasets to provide greater understanding of brain function and dysfunction. Such initiatives span normal function (Human Connectome Project), development (NIH Pediatric Database), and brain disorders, such as Alzheimer's disease (ADNI) and mental illness (Research Domain Criteria Project). Although these initiatives provide the necessary empirical foundation, quantitative tools are missing to integrate these multiple datasets to "reconstruct" the brain, and provide the link between these data and those from a single person.
Over the last 6 years, a neuroinformatics platform has been developed: The Virtual Brain (TVB) (Sanz Leon et al., 2013). The defining feature of TVB is that it generates personalized functional neuroimaging data based on individual structural connectome data to create personalized virtual brains. These models are specific to each individual person, and contain the connectivity between parts of the brain and the dynamics of local neural populations. TVB uses structural MRI data to create the custom brain surface, diffusion-weighted MRI data to infer the anatomical connections between brain areas, and then functional MRI data as the target to modify the parameters of the model to reproduce the observed functional data. The neuroinformatics architecture of TVB houses a library of models, which catalogues the biophysical parameters that produce different empirical brain states (Ritter et al., 2013). Global biophysical parameters represent biological mechanisms governing dynamics between brain regions, whereas the local biophysical parameters describe the properties of small populations of neurons integrating dynamics at the local mesoscopic level. That is, modeling in TVB comprises multiple scales of brain dynamics that are invisible to brain imaging devices, and therefore TVB acts as a "computational microscope", allowing the inference of internal states and processes of the system. TVB thus offers a novel platform to formulate biologically interpretable hypotheses on the effects of stroke and its recovery based on biophysical mechanisms governing brain dynamics. Beyond the direct clinical implications of network dysfunction in stroke, these insights can contribute a first step to the understanding of fundamental mechanisms of the brain's structure-function relationship. TVB has been established and applied to normative datasets (Deco et al., 2012) and for learning and plasticity (Roy et al., 2014), yet a proof of concept needs to be established based on pathological states.
The objective of the present study using the TVB platform was to determine changes in local and global biophysical parameters to better understand individualized brain dynamics after stroke. In this approach, the model parameters act as a means to assess brain health, analogous to blood samples assessing physical health, and hence, parameter changes could ideally be used as potential biomarkers of stroke and/or stroke recovery. So far, such biomarkers have mostly focused on stable architectures, from behavior to fine anatomical and functional levels (Burke and Cramer, 2013). In contrast, our aim is to create a synergistic amalgamation of mathematical models with neuroimaging, where the biomarker derives from the dynamical model itself.

Subjects
Twenty volunteers with chronic stroke (ages 23-74, 8 females) in the middle cerebral artery (MCA) territory and 11 age-matched controls were included in the study. This study and all procedures for recruitment and consent were approved by the Institutional Review Board of the University of Chicago and the University of California Irvine Medical School. Demographic details and stroke characteristics of our cohort can be found in Table 1.
Motor performance was assessed with the following: the Functional Ability Scale of the Wolf Motor Function Test (WMFT), Nine-hole peg test, the Fugl-Meyer upper arm test, and the Motor Activity Log (MAL-14). These assessments were collected at baseline (pre-therapy), af-ter 1 month of intensive hand therapy (post-therapy) and 6 -12 months after therapy (maintenance).

Virtual brain transplantation
Because of mechanical deformation consequent to large cortical strokes, the anatomical parcellation on T1w images using semiautomated methods is very difficult to achieve. Hence, a "virtual brain transplant" process was performed in accordance with a previous approach (Solodkin et al., 2010). This method replaces the cortical lesion with the homologous image from the contralesional hemisphere from the same subject. With this, brain parcellation is possible using semiautomatized software. The process consisted of the following steps: (1) Lesion segmentation by hand.
(2) Using the AFNI 3dcalc function (Cox, 1996), the homologous region in the nonlesioned hemisphere was dissected and transplanted into the stroke region, effectively filling in the missing portions of the brain. (3) Manual corrections were then done in the interface between the native and transplanted T1-w images by visually examining each voxel and making voxel intensities uniform using AFNI's 3dLocalStat and 3dcalc commands. (4) The brain was then parcellated into 96 cortical and subcortical regions. The original parcellation based on a macaque template (Van Essen, 2004) was transformed to the human MNI template via PALS (Van Essen, 2005). To increase accuracy, the deformation process was carried out using landmarks (based on CARET) and functional activation patterns considered homologous between the two species (Van Essen and Dierker, 2007).

Diffusion tensor imaging
Preprocessing of DTI data consisted of the following: (1) motion correction using the FSL eddy current correction (Leemans and Jones, 2009), (2) generation of a binary brain mask from the b0 image and application of the mask to all diffusion images using the Brain Extraction Tool from FSL (Smith, 2002), (3) fitting of a diffusion tensor at each voxel using FSL's dtifit function, (4) nonlinear coregistration of T1 data to the MNI brain and coregistration of T1 images to their respective DTI images producing an MNI to DTI transformation using ANTS (Avants et al., 2011), (5) white and gray matter segmentation performed on the MNI-to-T1 atlas using FAST (Zhang et al., 2001), and (6) parcellation of the gray matter into 96 regions as described above and registration of these regions to the DTI using the T1-to-DTI transformation with a nearest neighbor interpolation.

Tractography and structural connectivity matrix generation
Probabilistic tractography was performed to trace the fiber bundles associated with pairs of cortical regions in the MNI space, which were defined as edges in the network (Ritter et al., 2013;Zalesky and Fornito, 2009).Two connectivity measures were extracted: (1) capacities, depicting the maximum rate of transmission of information through edges, were calculated using the number of streamlines at the minimum cross-sectional area of an edge (Zalesky and Fornito, 2009); and (2) distances, defined by the lengths of each edge, were calculated by averaging the lengths of all streamlines in an edge. These measures were used to generate two 96 ϫ 96 structural connectivity matrices. Quality assurance to reduce falsepositives was performed on each structural connectivity matrix by a trained neuroanatomist (A.S.).

Resting state fMRI preprocessing
Preprocessing was done in AFNI (Cox, 1996) and included the following steps: motion correction of functional and anatomical datasets (Cox and Jesmanowicz, 1999), 3D spatial registration to a reference acquisition from the rsfMRI run, registration of functional images to the T1-w volume, de-spiking and mean normalization of the time series, motion correction (Ͼ1 mm; Johnstone et al., 2006) and regression of CSF and white matter signals to remove slow-wave components (eg, physiological noise; Lund et al., 2006).

Resting state fMRI postprocessing
Average time series were extracted for each of 96 MNI regions. For each subject, a 96 ϫ 96 functional connectivity matrix was generated by calculating the pairwise correlation of the time series for each region (Ritter et al., 2013) using the "corr" function in MATLAB.

Modeling in TVB
The Virtual Brain (TVB v1.08; Fig. 1) was used for all simulations (Sanz Leon et al., 2013) where the principal empirical input to the platform is the structural connectivity matrix derived from each individual subject's tractography. Based on this input, TVB simulates field potentials by integrating global dynamics with a local (mesoscopic) model that determines the dynamics within brain regions. Following, BOLD signals are derived from the generated field potentials. In this work, we used the Stefanescu-Jirsa 3D (SJ3D; Fig. 2) local model, as the resulting mean field model does not rely heavily on synaptic delays (Stefanescu and Jirsa, 2008;Jirsa and Stefanescu, 2011;Sanz-Leon et al., 2015), making it compatible with the poor time resolution associated with BOLD signals. Specifically, the SJ3D model is derived from populations of bursting neurons and includes six states describing excitatory and inhibitory dynamics via the inclusion of a variety of biophysical parameters defining the local mean fields (for a list of the parameter values used in the present study see Table 2; Hindmarsh and Rose, 1984;Stefanescu and Jirsa, 2008). The following sequential steps were performed for each individual subject: (1) Importing of a subject-specific connectivity matrix into the TVB platform.
(2) Selection of the SJ3D local model.
(3) Parameter space estimation (exploration and fitting). We sequentially performed systematic parameter space explorations and fitting to determine the optimal values for global and local parameters in all subjects. (1) Parameter space exploration: we used heat maps of global variance (mean variance of time series across all brain regions) to constrain the range of values for each model parameter (Fig. 3). The range of values considered is assessed based on those values with high global variance flanked by bifurcation points (Breakspear and Jirsa, 2007). An additional advantage of this approach is that it is not only pragmatic but it can also provide information on the degree of variability and sensitivity that parameter values Figure 2. Equations of the Stefanescu-Jirsa 3D model. A, Evolution equation implemented in TVB to simulate brain activity. The mean field potential x i (t) of a region i at time t is dependent on the local dynamics f(x i (t)) provided by the Stefanescu-Jirsa-3D model, the long-range structural connectivity w, which links regions i and j and is provided by the input of individual structural connectivity matrices (weights), and noise ͑t͒. Time delays (⌬t) are distance dependent and are provided by the structural connectivity matrices (lengths). All mathematical details of the model and its numerical implementation are provided by Sanz-Leon et al. (2015). B, Equations comprising Stefanescu-Jirsa 3D. The first three equations (, , ) represent the excitatory subpopulation of neurons within a local region, whereas the last three equations (␣, ␤, ␥) represent the inhibitory subpopulation of neurons in that region. IE and II denote the input current to the excitatory and inhibitory populations of each node, respectively. The first of each of the two sets of equations accounts for neuron potentials. The second and third equations account for the transport of ions across the membrane through ion channels. Note that the dynamics of these populations are dependent on the interactions between inhibitory and excitatory influences (K 12 , K 21 , K 11 ). have onto the simulated signals.
(2) Parameter fitting: the final optimal value was subsequently obtained by assessing the specific value for the parameters that resulted in the best fit between the empirical and simulated signals based on three metrics described below (Step 6). The global parameters explored included conduction velocity and global coupling and the local parameters included K 12 (excitatory on inhibitory coupling), K 21 (inhibitory on excitatory coupling), and K 11 (excitatory on excitatory coupling). The local parameters were chosen as they have the strongest impact on the dynamics of the SJ3D model (Stefanescu and Jirsa, 2008).
(4) Stochastic network simulation: based on the values obtained in the parameter space exploration, we generated field potentials with the same duration (4 min) and sampling rate (TRϭ2 s) as the empirical rsfMRI acquisition. The length of the simulated data was kept equal to the length of the empirical data to minimize the influence of variability over the course of the time series, as it is becoming increasingly patent that values of functional connectivity are not stable over time (Hutchison et al., 2013). White noise with Gaussian amplitude (mean ϭ 0, SD ϭ 1) was added to each node. Numerical integration of the system was performed using stochastic Heun's method (Mannella 2002), with an integration step size of 0.0122 ms.
(5) The BOLD signals were derived from the field potentials using a hemodynamic response function implemented with a gamma kernel (Boynton et al., 1996;Sanz-Leon et al., 2015).
(6) Assessing reliability of the simulated time series: comparison between the empirical and simulated BOLD time series was done in terms of amplitude, frequency, and phase. Amplitude: we calculated the range of amplitude by identifying the highest and lowest peaks present in the time series across all regions. Frequency: fast Fourier transforms of the raw and simulated time series were obtained using MATLAB's "fft" function with a sampling frequency of 0.5 Hz, to determine the range, profile, and peak frequencies (Ritter et al., 2013). Phase: this was assessed by comparing the functional connectivity matrices of the simulated and empirical time series. We averaged all matrices from healthy controls to obtain a group control matrix, and calculated the pairwise linear correlation coefficient between the simulated functional connectivity matrix for each individual to the group.
(7) Differences in parameter values between healthy controls and stroke cases were evaluated with Wilcoxon sum rank test corrected for multiple comparisons (Bonferroni).
(8) Relationship with clinical phenotype. To determine whether there was any relationship between TVB parameters and the clinical phenotype, multiple linear regression was performed between model parameters (dependent variables) and the following independent variables: motor outcome measures (Fugl-Meyer, WMFT, 9-hole peg, and MAL-14), patient demographics (age, sex, presence of depression), and lesion characteristics (size, location, time after stroke, side of stroke).

BOLD simulations generated with TVB correlated with the empirical BOLD responses
The frequency spectrums of the simulated and the empirical BOLD responses had similar ranges (0 -0.25 Hz) and mean peak (empirical ϭ 0.05 ϩ 0.035 Hz; simulated ϭ 0.03 ϩ 0.023 Hz; Fig. 5). Although the mean amplitudes were similar (empirical ϭ 8.15; simulated ϭ 9.49), the range of values was wider in the empirical signals (0.17-87.43) than those found in the simulated BOLD (3.79 -22.64). The relative phases of the regions within simulated and empirical time series were similar as assessed by the mean correlation coefficient between their respective functional connectivity matrices (mean ϭ 0.27Ϯ0.02; p b ϭ 0. 9e-12 Fisher Z-transformation). These validated simulations provided us with specific parameter values at both the global and the local levels associated with healthy control subjects and after stroke.

Stroke was associated with reliable changes in global and local parameters
Although qualitative in nature, the color-coded graphic representation of the variance distribution done as part of the parameter space exploration (Table 3; Fig. 3) provides a glimpse into differences of combined values for the two global parameters: global coupling (x-axis) and conduc-tion velocity (y-axis) in healthy controls and in stroke subjects, with warm colors representing higher variance. These explorations demonstrated at this early stage of analysis that the range of optimal parameter values (hot colors) in controls had similar topology of the distribution of variance, as well as concrete values. In contrast, stroke cases displayed high variations in both topology and values, where although some had similar distribution patterns as the healthy controls, others had scattered, fragmented patterns. Similar observations were found with respect to local parameters (Table 4).
Numerically, differences in parameter values between healthy controls and the stroke cohort are as follows:

Global parameters
(1) Conduction velocity: the range of modeled conduction velocities obtained via TVB in healthy controls ranged from 45 to 90 m/s with a mean of 62 Ϯ 10 m/s. In contrast, the conduction velocities in stroke subjects had a range between 12 and 80 m/s with a mean of 46 Ϯ 21 m/s. Comparison between the two groups with Wilcoxon rank sum test (p c ϭ 0.05) was marginally significant after correction for multiple comparisons.(2) Global coupling (rescale factor of incoming activity linking global with local dynamics): in healthy controls, the mean was 0.053 Ϯ 0.009 (range, 0.044 -0.047) and in cases with stroke the mean was 0.061 Ϯ 0.016 (range, 0.04 -0.09). Wilcoxon sum rank test showed this difference was significant after correction for multiple comparisons (p c ϭ 0.013).
In addition, it is important to note that the trend in all stroke cases where the values were different from those in controls was consistent: that is, it presented always as a decrease in conduction velocities (N ϭ 12) and an increase in global coupling (N ϭ 14). The rest of the stroke cases did not show differences with healthy controls.
(3) K 11 (influence between excitatory populations): the values of K 11 in controls had a mean of 0.833 Ϯ 0.142 (range, 0.6 -0.95) and in stroke cases had a mean of 0.613 Ϯ 0.301 (range, 0.1-0.99). Comparison between the two groups with Wilcoxon sum rank test gave p c ϭ 0.1.
In summary, compared to values in healthy controls, there was a higher global coupling and a decrease of local inhibitory dynamics represented by the local parameter K 21 along with a trend toward a reduction of conduction velocity.

Global and local parameters were correlated with clinical phenotype
Multiple linear regression analysis to establish a relationship between modeling parameters and some clinical metrics did not show a correlation. The following clinical elements were considered in this preliminary assessment: stroke phenotype (size, location, time after stroke, side of stroke), depression, patient demographics (age, sex), and severity of impairment.
Next, we assessed the relationship between parameter values with recovery from stroke immediately after therapy and after 1 year (maintenance) using a multiple linear regression. This analysis showed a negative relationship between K 12 and Fugl-Meyer scores both post-therapy (t ϭϫ2.386; p d ϭ0.038) and at maintenance 1 year later (t ϭϪ3.824; p d ϭ0.005). In addition, global coupling had a positive relationship with the Wolf Motor Function Test (tϭ 2.461; p d ϭ0.039) at maintenance. Thus, these two parameters derived from modeling based on pre-therapy conditions were related to long-term motor gains rather than the physical features of the stroke or the patient's demographics (Fig. 6).

Discussion
The main result of the study showed that the simulation of BOLD signals using TVB in stroke enables the identification of key changes associated with large-scale neural dynamics in individual patients. Overall, our results showed that, compared to healthy controls, individuals with stroke have a consistent reduction in conduction velocity and a relative increase in local-over-global brain dynamics. Further, the identified parameters were related to functional outcomes such that these parameters predicted long-term recovery after therapy. Together, these results not only back TVB as an effective tool in identifying dynamic brain changes in stroke spanning multiple scales, but also specifically identify potential predictors of recovery in stroke at the individual level. This study suggests that TVB may be a powerful platform for the application of large-scale modeling in understanding brain mechanisms at an individual subject level.

Stroke is related to consistent global and local parameter changes
The successful simulation of empirical rfMRI data in this study facilitated a particularly salient finding; the dynamic model derived from stroke subjects had a significant decrease in the local parameter K 21 and a consistent global coupling increase, accompanied by a trend in decreased conduction velocity. Two aspects of these results are of special interest: the first relates to the nature of the statistical outcomes and the second to the biological interpretation of these changes.
(1) Imaging-derived metrics in humans in general have high variance (Mueller et al., 2013); consequently, analytical measures have been developed to minimize it (Fischl et al., 1999). Further, this variance is amplified by stroke (Rehme et al., 2012), and has compelled researchers to stratify patients with precise criteria (Cramer, 2010), resulting in low sample sizes and high inter-study variability. In contrast, even when we used minimal exclusion criteria when selecting participants, changes seen after stroke were highly consistent, where all the cases that had a parameter change with respect to controls had the same directionality and relatively low variance. Given the high level of subject variability (as expected for a cohort including a large range of clinical phenotypes), we find this consistency somewhat surprising. However, we are not suggesting high reliability of our modeling, as the definitive answer will result from expanding the assessment to a larger population where the predictive value of the parameter changes can be formally assessed.
(2) Stroke survivors exhibited a significant decrease in K 21 , a parameter at the mesoscopic level that represents the influence of inhibitory on excitatory neuronal populations. A decrease in K 21 thus indicates local disinhibition. These results are highly consistent with existing data on the basic mechanisms of stroke at the cellular level. For example, rodent models of MCA stroke show an imbalance in the density of excitatory and inhibitory receptors in tissue surrounding the lesion (Schiene et al., 1996). Specifically, they suggest a decrease in gammaaminobutyric acid receptor expression in widespread ipsilesional cortical areas and a concomitant increase of N-methyl-D-aspartate receptor expression in the contralesional hemisphere.
In the context of stroke in humans, hyperexcitability has been described in two experimental paradigms: (1) Studies using TMS to test cortical excitability after stroke have shown a decrease in the current needed to elicit motor evoked potentials and an increase in their amplitude (Hallett, 2007) along with an expansion in the area producing them (Liepert et al., 2000) suggesting disinhibition in motor cortices (Shimizu, 2002). Furthermore, decreasing the hyperexcitability via repetitive lowfrequency stimulation (Takeuchi et al., 2005) along with a reduction of the TMS stimulation area (Liepert et al., 2000) has been related to motor recovery (Hallett, 2007).
(2) Increased activity in motor and non-motor regions has been reported in fMRI studies after stroke (Rehme and Grefkes, 2013). Specifically, increased contralesional activity has been observed (Weiller et al., 1992;Ward, 2003;Grefkes et al., 2008). Although this has been explained as a recruitment of supplementary areas to assist movement (Rehme and Grefkes, 2013), others have related it to widespread cortical hyperexcitability (Buchkremer-Ratzmann et al., 1996), suggesting long-range corticocortical inputs (Logothetis et al., 2001) with increased activation via decreased inhibition (Liepert, 2003;Blicher et al., 2009). Functional recovery has in turn been associated with the degree of recovery of activity in the affected cortical areas (Cramer, 2008). Complementing the above, our results show a correspondence between local and global levels. Indeed, the reduction in local inhibitory influence over excitatory populations was accompanied by an increase in global coupling, reflecting an imbalance after stroke between global and local brain dynamics, favoring the latter. That is, local dynamics exert a stronger influence than global dynamics following stroke. In this case, the imbalance could be exacerbated by the decrease in conduction velocity. Interestingly, this imbalance has also recently been modeled in other brain diseases. For example, early stages of schizophrenia have been associated with a breakdown of local dynamics occurring prior to the disruption of global dynamics occurring later on in disease progression (Rubinov et al., 2009;van den Berg et al., 2012).
A particularly interesting finding was the trend associated with a decrease in conduction velocity in individuals with stroke, as it has previously been described through measurements of central motor conduction times (CM-CTs) via transcranial magnetic stimulation (TMS) in the primary motor cortex. Immediately following stroke, CMCT decreases and correlates with functional measures (Abbruzzese et al., 1991;Pennisi, 2002) tending toward an incomplete normalization over the long-term (Heald et al., 1993). That being said, there is a paucity of information on decreased conduction velocity on corticocortical connections. The bulk of knowledge derives from studies in rodents showing structural changes to axons and oligodendrocytes in the primary lesion and the ischemic penumbra (Rosenzweig and Carmichael, 2015). In addition, although some degree of remyelination occurs in the recovery phase, the process is often arrested before completion (Syed et al., 2008). In human autopsy samples, there is an increase in nodal and paranodal lengths adjacent to lacunar lesions (Hinman et al., 2015), which may lead to decreased conduction velocities (Rasband, 2011). Our results thus provide direction for future animal studies, exemplifying the translational nature of TVB findings.
TVB thus appears to be effective at modeling brain activity in healthy brains and those impacted by disease processes, and has the novel capability of studying brain dynamics at multiple scales, including at a level that has thus far only been available via animal models or surrogate neuroimaging markers in humans. Applying this method of modeling, which is tied directly to biological mechanisms, to existing large datasets opens up the possibility to experiment with expanded models of brain states, including a myriad of diseases and their potential treatments.

Potential predictors of motor recovery after stroke
Our results demonstrated that local (K 12 ) and global (global coupling) parameters, derived from pre-therapy conditions, were significantly correlated with motor gains post-therapy and at maintenance. Furthermore, both parameters point in the same direction, as poor recovery was associated with an increase in local excitatory influences and with an emphasis on local dynamics, whereas values closer to controls correlated with better recovery.
Interestingly, TVB parameters in stroke did not correlate with severity of disease at the pre-stroke time point, even though the structural connectivity matrix used in the modeling coincided with this time point. In addition, other physical features of the stroke (size, location) or patient demographics (sex, age) did not correlate with the modeled parameters. Finally, neither lesion characteristics nor patient demographics correlated with recovery, highlighting the unique predictive potential of these parameters.
The question then becomes to what extent these parameter estimates can be used as predictors of recovery at the individual patient level. Although a cross-validation approach using the current dataset could serve to answer this question, a new and larger stroke cohort is ideal in obtaining estimates of the sensitivity and the specificity of our markers, due to high variance in stroke. However, there is clear value of our observations even with this limitation. At present, biomarkers for stroke recovery have been limited by the use of "substitute or surrogate" measures derived from brain imaging or electrophysiology, mainly due to the inability to measure in vivo more ideal basic elements, ie, at molecular or cellular levels (Burke and Cramer, 2013). Indeed, such elements may be observed more closely in animal models, but are difficult to translate to humans due to the limited homology between species. Specifically, the Stefanescu-Jirsa 3D model used in this study evolved from the mesoscopic level Hindmarsh-Rose model. The Hindmarsh-Rose model itself is rooted in the principles of the Hodgkin-Huxley neuron equations, in addition to dynamics based on bursting neurons found in invertebrate circuitry (Hindmarsh and Rose, 1984). Further, the neural behaviors described by the Hindmarsh-Rose model have been biologically verified in other animal models (Selverston and Ayers, 2006;Gu, 2013). Therefore, although any model of the mesoscale does not encompass the complexity of brain processes at the cellular level, there is likely emergence of behavior from the cellular level to the mesoscopic level, exhibiting deterministic behavior that can be modeled and also observed in vivo.
That is, the transition between the macroscopic and microscopic level is represented by population dynamics at the mesoscopic level (Mitra, 2014). From this, one could conclude that the path toward basic biomarkers should include the intermediate mesoscopic level. Indeed, TVB allows one not only to estimate parameters at that level but also to link it to the macroscopic global wholebrain level. TVB is not unique in considering biophysical parameters as exemplified by inference models based on DCM (Moran et al., 2011). Basically, there are no conceptual differences in the inferential goals between TVB and DCM but they do differ in the detailed mechanics. For example, whereas TVB develops the model at the level of large-scale networks, DCM focuses on portions of these networks. Second, and perhaps the key contrast is that whereas DCM fits the parameter of the model but does not generate data, TVB uses the model to generate data, making these two approaches highly complementary.
An interesting and unique aspect of TVB is its highly individualized approach, as parameter estimates are derived from individualized structural connectivity matrices obtained from each subject, and hence, it can provide the first step to customize individual therapeutic interventions. For example, our ongoing work is beginning to test potential "virtual interventions" by modifying specific parameters changed after stroke and determining the degree of restoration of brain dynamics on each stroke patient.
A second ability of this modeling approach is to use the model of an individual patient's brain connectivity that can be objectively measured and evaluated as an indicator of normal biological processes (such as, resting state activity, rsfMRI), pathogenic processes, or pharmacologic responses to therapeutic intervention (Biomarkers Definitions Working Group, 2001). Dynamics of rsfMRI are highly nonstationary (Allen et al., 2014) and existing metrics, including the direct correlation between functional and structural connectivity, are so far incapable of addressing this issue satisfactorily (Goni, 2014). A number of studies have therefore used generative modeling to parse the relationship between structural and functional connectivity. A recent study (Andersen et al., 2014) demonstrated that the fusion of TVB-like network modeling with structural neuroimaging explains the nonstationary dynamics observed in rsfMRI. Thus, we propose a conceptual paradigm shift, in which the dynamic model shifts the nonstationary functional data from imaging at the mesoscopic scale to a more deterministic set of coefficients in a brain model. In other words, complex dynamics cannot be captured by stationary imaging analyses, but can be generated by a data-constrained mechanistic model of brain-circuit dynamics, as seen in the generative modeling approach detailed in stroke (Brodersen et al., 2011). Thus, the mathematical model could be seen as a compact generator of dynamics-based biomarkers, or even as the biomarker itself. The primary benefit, as we demonstrated here, is that it becomes easier to understand disease mechanisms by evaluating the coefficients of the model.
Of note, the approach used in this study to validate the simulated time series was to compare frequency, amplitude, and phase of the simulated and empirical signals. After the refinement of the TVB models, future studies will incorporate a larger variety of multidimensional analyses, particularly with respect to temporal variability in resting state signals. Furthermore, the current study determined optimal values of local parameters applied to all brain regions. Future studies will focus on local parameters for subsets of brain regions, eg, changing parameters of nodes within and/or around a stroke lesion to determine how this impacts the resultant simulated brain activity. We also note that the translational power of our findings depends upon the reproducibility of parameters for a given brain state, the answer for which will emerge with expanded application of TVB to other cohorts. The results from this study thus confirm that TVB allows the assessment of biophysical variables previously unattainable in human studies. This method provides a potentially important and novel application of large-scale modeling, in which we can probe brain dynamics and biomarkers on an individual level. Therefore, The Virtual Brain has the potential to become an important step toward the development of individualized medicine in stroke.