Quantifying motor recovery after stroke using independent vector analysis and graph-theoretical analysis.

The assessment of neuroplasticity after stroke through functional magnetic resonance imaging (fMRI) analysis is a developing field where the objective is to better understand the neural process of recovery and to better target rehabilitation interventions. The challenge in this population stems from the large amount of individual spatial variability and the need to summarize entire brain maps by generating simple, yet discriminating features to highlight differences in functional connectivity. Independent vector analysis (IVA) has been shown to provide superior performance in preserving subject variability when compared with widely used methods such as group independent component analysis. Hence, in this paper, graph-theoretical (GT) analysis is applied to IVA-generated components to effectively exploit the individual subjects' connectivity to produce discriminative features. The analysis is performed on fMRI data collected from individuals with chronic stroke both before and after a 6-week arm and hand rehabilitation intervention. Resulting GT features are shown to capture connectivity changes that are not evident through direct comparison of the group t-maps. The GT features revealed increased small worldness across components and greater centrality in key motor networks as a result of the intervention, suggesting improved efficiency in neural communication. Clinically, these results bring forth new possibilities as a means to observe the neural processes underlying improvements in motor function.

Introduction

The loss of hand function is often cited as the most devastating consequence following a stroke (Barker and Brauer, 2005). Over the past several decades, many novel rehabilitation approaches for this functional loss have developed in parallel with advances in neuroscience (Cauraugh et al., 2010; French et al., 2007; Langhorne et al., 2011; Thieme et al., 2013; Veerbeek et al., 2014; Waldner et al., 2009). The most effective treatments generally involve high doses of functionally meaningful tasks beginning in the early weeks post stroke. However, while practice guidelines for effective treatment are emerging, the lack of strong evidence regarding the content of these treatment paradigms remains a critical concern. Findings that certain treatments produce beneficial outcomes in some subjects, while are ineffective or even detrimental in others (Cauraugh et al., 2010; Coupar et al., 2010; Langhorne et al., 2009; Winter et al., 2011), underscores a crucial need to characterize individual differences in treatment responsiveness.

Given the extent of brain reorganization due to a stroke lesion (Westlake and Nagarajan, 2011; Westlake et al., 2012), it is essential to consider the marked variability in brain networks as a primary contributor to these differences. However, frequently used neuroimaging analysis tools are unable to adequately capture and efficiently use the statistical variability — within or across subjects with various lesion locations and reorganization patterns due to stroke. Hence, it is very important that the analysis approach is multivariate rather than mostly univariate as in the frequently used general linear model (GLM) and seed-based connectivity approaches. To this end, both data-driven and multivariate independent components analysis (ICA) have surfaced as effective alternatives. The commonly used group ICA (Calhoun et al., 2001) has been successful for the analysis of data from multiple subjects, but its ability to capture individual subject differences — which is especially critical in group comparisons of subjects with stroke — is limited. IVA of multiple datasets is a recent generalization of ICA that has been shown to be superior in capturing subject variability and making full use of the statistical information within and across multiple datasets (Adali et al., 2014; Ma et al., 2013; Michael et al., 2014). Moreover, the addition of graph theoretical (GT) approaches to IVA components creates a unique and powerful metric for interpreting brain states after stroke.

The aim of this study was to identify changes in intrinsic brain networks as a result of a task-oriented unilateral intervention of the affected upper extremity in subjects with stroke. We applied the GT framework to generate neural connectivity metrics within canonical cognitive and sensorimotor neural networks, as identified using IVA. We hypothesized that following the intervention, subjects would demonstrate increased sensorimotor network communication efficiency (i.e., increased clustering and decreased path lengths) that would correlate with clinical tests of affected upper extremity function. Data related to the behavioral effects as a result of the unilateral intervention in this study are part of a larger study comparing two different interventions, but these data have not yet been published or presented.

Methods

Study population

Ten volunteer subjects with a first ever monohemispheric stroke in the distribution of the middle cerebral artery participated. Additional inclusion criteria were as follows:

Score of 25–35/66 on the Upper Extremity Fugl–Meyer Assessment Scale, indicating moderate to moderately severe impairment in arm function.

50–80 years of age.

Completion of all conventional rehabilitation.

Able to actively perform a range grasp (finger flexion) with at least two fingers and the wrist maintained in a neutral position.

Subjects were excluded from the study for the following reasons:

Participation in an arm exercise program more than 20 min 3 times a week (to avoid bias from training).

Medical history of (a) recent hospitalization for severe disease or surgery, or (b) significant orthopedic or chronic pain conditions limiting upper extremity exercise.

Neurological history of (a) untreated major post-stroke depression, (b) dementia based on Folstein Mini-Mental Status Score or (c) severe receptive or global aphasia that confounded testing and training (unable to follow 2 point commands).

A non-stroke neuromuscular disorder restricting exercise (e.g., Parkinson's syndrome).

Independent ability to grasp, lift, transport and release a 6 inch diameter foam ball.

Indwelling metals, medical implants, or claustrophobia incompatible with MRI.

The study protocol was approved by the Institutional Review Board and was conducted in accordance with the Declaration of Helsinki. All participants gave written informed consent prior to study participation.

Rehabilitation intervention

The intervention consisted of two 6-week training blocks of three, 1-h training sessions/week for a total of 36 sessions. Because the subjects included in this study demonstrated moderate to moderately severe arm function, the aim of the first 6 weeks was to prime the affected upper extremity to improve voluntary motor control and active range of motion such that functional tasks could be more effectively practiced during the second 6 week session. Therefore, during the first block, training was focused on repetitive, active movement of the affected arm into shoulder flexion and elbow extension. The Tailwind, bilateral arm training with rhythmic auditory cueing (BACTRAC) device was used for all sessions, which was previously shown to improve arm function in subjects with chronic (> 6 months) stroke (McCombe Waller et al., 2008). Training at each visit included four 5-minute training bouts with a 10-minute rest period between bouts. Rate was initially set at a preferred speed and then increased as tolerated. Adjusting the height of the BATRAC device added resistance to the training. The device was flat during the first week and the height was increased as tolerated for the remaining weeks of treatment. During the second 6-week session block, the focus was on task-specific training that included use of the affected arm and hand in grasp, reach, and release activities. This study focused on the group before and after the first 6 weeks of training.

Behavioral measures

A testing battery of clinical assessments was conducted on each subject before and after the intervention and included the following measures:

The Upper Extremity Fugl–Meyer Assessment Scale was used to assess motor impairments.

The Wolf Motor Function Test was used to assess motor function in terms of movement time and quality.

The University of Maryland Arm Questionnaire for Stroke (UMAQS) was selected to measure daily use of the paretic arm.

Isometric grip strength was assessed with a hand dynamometer and pinch strength (lateral pinch and 3-jaw pinch) was assessed using the Jamar Pinch Gauge.

The Box and Block Test was used to assess hand function.

To determine the combined clinical effects of affected upper extremity training on measures of brain connectivity, a composite measure of the difference in all clinical outcomes was calculated and used as the single primary outcome measure. Intersubject baseline variability was normalized by use of a change score in which the difference between baseline and post-intervention was divided by the baseline score of each outcome measure.

Neuroimaging acquisition

A 3 T Phillips Archieva MR scanner was used for all magnetic resonance (MR) data acquisition. fMRI data were acquired using an echo-planar imaging (EPI) sequence with the following scanning parameters: TR = 3000, TE = 30, flip angle = 75, voxel size is 2.5/2.5/3, 50 slices/3 mm slice thickness with no gap. Data were collected while participants performed an auditory-cued affected hand grasp and release task. The hand motor task was standardized based on available range of finger flexion achieved with minimal effort. An MRI compatible SAEBO-FLEX exoskeleton (Saebo, Inc., Charlotte, NC) was worn by each participant to assist with finger extension at the completion of each finger flexion task. The task was cued once every 3 s in a block design of 24 s of motor task followed by 30 s of rest. Subjects were instructed to squeeze lightly at each beep, then to relax and wait for the next beep. Eight blocks were completed during each session. Practice blocks outside of the scanner ensured consistency in task performance and the absence of mirror movements of the unaffected hand. In addition to the fMRI data, a high-resolution 3-dimensional T1-weighted, image was acquired for anatomical reference.

fRMI preprocessing

The SPM8 software package (Friston, 2013) was employed to perform fMRI preprocessing. Slice timing was performed with a shift relative to the acquisition time of the middle slice using sinc-interpolation. All images were spatially realigned to the 1st volume to correct for inter-scan movement. To remove movement-related variance the realigned images were processed using the unwarp SPM8 function. Given the observed variance (after realignment) and the realignment parameters, estimates of how deformations changed with subject movement were made, which subsequently were used to minimize movement-related variance. Data were then spatially normalized to a standard echo-planar imaging template based on the standard Montreal Neurological Institute space (Friston et al., 1995) with an affine transformation followed by a non-linear approach with a 4 × 5 × 4 basis functions. Images, originally collected at 2.5 mm × 2.5 mm × 2.5 mm, were resampled at 2 mm × 2 mm × 2 mm voxels. Finally, data were spatially smoothed with a Gaussian kernel of full-width half maximum (FWHM) of 5 mm × 5 mm × 5 mm.

fMRI analysis

IVA

IVA is a recent extension of ICA to multiple datasets. A key feature of IVA is that it takes dependence across datasets into account while achieving joint blind source separation. IVA has also been shown to better capture subject variability in ICA comparison studies (Dea et al., 2011; Ma et al., 2013; Michael et al., 2014). For K datasets, IVA models each as linear mixtures of N statistically independent sources. For spatial IVA and ICA, which is the most commonly used model, a source can also be referred to as a spatial component or spatial map. Each dataset is represented as a random vector using superscript notation as x[ = [x1[, …, x[, …, x[], k = 1, …, K. The IVA generative mixture model for x[ iswhere each A[ is an N by N mixing matrix and s[ = [s1[, …, s[, …, s[] is the random vector representing the original sources for the kth dataset. When K = 1, IVA is equivalent to ICA:

For K > 1, IVA simultaneously estimates a demixing matrix W[ for each datasetwhere y[ is the random vector representation for the estimated sources for the kth dataset. Dependence across datasets is exploited in the IVA model by allowing each source to have statistical dependence with one source from each other dataset. Sources across datasets are placed in a vector called a source component vector (SCV) and the mutual information within each SCV is maximized as part of the IVA cost function while maximizing the independence between each SCV. In this study, sources refer to spatial components.

Order selection

Order selection refers to selecting the number of spatial components to be estimated during analysis; the number of components is referred to as the order. Order selection is important due to the high dimensionality and high noise level in fMRI data. However, the problem of estimating the number of clinically meaningful components is not straightforward due to noise, such as head movement artifacts and ambient noise, and the dependence among samples, i.e., voxels in the spatial domain. For data obtained from subjects with stroke, high subject variability due to lesions and the process of brain plasticity makes the order estimation more difficult. One approach, as in Li et al. (2007), has been to subsample to a set of independent and identically distributed (i.i.d.) voxel samples in order to take advantage of information-theoretic criteria that require this assumption. More recent work jointly estimates order and downsampling depth to produce i.i.d. samples using information-theoretic criterion leading to slightly higher orders (Li et al., 2011). The higher order used in Li et al. (2011) produced components of interest in the frontal, parietal, and temporal regions that were not observed at lower order ICA. Using this approach (Li et al., 2011), the average order for the current study was found to be 61 with a standard deviation of 35. Due to the large standard deviation, the effect that order has on component estimates was investigated. At lower orders with low t-values, few meaningful components were estimated. As the order was increased, t-values increased and components began to reflect canonical networks. Therefore, components at a range of orders from 20 to 100 were estimated in steps of 10. At order 80, the component estimates produced more meaningful components with higher t-values than those at lower orders. At orders higher than 80, the components become fragmented into separate spatial maps. The low to high frequency power ratio for each of the 80 components was then analyzed using the Group ICA of fMRI Toolbox (GIFT). Meaningful components were retained by removing artifacts found outside of gray matter (in ventricles, skull and surrounding tissue, white matter, and eye sockets), high frequency time courses due to scanner noise, and rings of signal representing motion. Other work has also made use of higher orders to provide a more detailed view of functional activity (Allen et al., 2011; Erhardt et al., 2011; Kiviniemi et al., 2009; Ma et al., 2012).

Validation

Calculated group features are validated by subsampling each group at one-third the group size to form a subgroup. A third of the subjects are randomly selected from the pre-intervention dataset and the corresponding subjects are selected at post-intervention. Repeating this random sampling with replacement, twenty subgroups are formed. The given feature is recalculated for both the pre- and post-subgroups. A two-sample t-test is then used to validate whether random subgroups demonstrate the same difference as when the entire population is considered. A similar approach to functional network connectivity (FNC) validation has been used by Jafri et al. (2008). FNC is computed by calculating the temporal dependence among spatial components usually through covariance and is thus a measure of connectivity across components. Spatial FNC (sFNC) measures connectivity in terms of the dependence among spatial components calculated using mutual information (Ma et al., 2014). Each spatial component can be thought of as a maximally independent brain state. Understanding how each state shares mutual information with each other leads to an understanding of sFNC. Thus, changes in FNC and sFNC give us a way to quantify differences in functional networks between pre- and post-intervention.

Graph-theoretical analysis

GT analysis simultaneously takes all clinically meaningful components into account to form group features such as the clustering coefficient, path length, small worldness, and node centrality. For each subject, an undirected graph is formed where each node corresponds to an estimated temporal or z-scored spatial component for the analysis of FNC and sFNC, respectively. For sFNC, each edge is defined by the mutual information between nodes and for FNC, edge values are defined by the correlation coefficient between nodes. Graph densities are then formed by removing edges with values below a chosen threshold. When the threshold is low, the graph retains many connections and is comparable to a random network with an equal number of edges and nodes. As the threshold increases, graphs with different connection densities are formed resulting in trends in the graph-theoretical features.

Random subject sampling is used for pre-intervention and post-intervention with replacement as described in Section 3.3. At each threshold, the GT features are calculated for each subject then the corresponding subgroup feature is the average taken over the subjects' features. Differences between groups were assessed using a two-sample t-test across subgroups with the FDR corrected level of significance set at alpha 0.05.

Clustering coefficient

The clustering coefficient is calculated as in Liu et al. (2008)

A cluster exists if two neighbors of the nth node are connected. The number of edges connected to the nth node is N. Therefore, N(N − 1)/2 is the maximum number of possible clusters and E represents the actual number clusters at the nth node. This clustering coefficient is then normalized by the corresponding coefficient from the comparable random graph with the same number of nodes and edges as the observed graph. These values are then averaged over subjects for each threshold.

Path length

At each threshold the average shortest path length for the ith node is calculated as defined in (Ma et al., 2013)where L is the number of edges along the shortest path from node i to node j. Then the graph's average shortest path length L is computed by averaging over L for i = 1, …, N. Smaller values of the average path length indicate greater efficiency in the connectivity structure (Hyvärinen et al., 2001). The average path length is then normalized by the corresponding average path length from the comparable random graph with the same number of nodes and edges as the observed graph.

Small worldness

Small worldness is calculated by dividing each subject's clustering coefficient by the average path length. Small worldness values are also averaged across subjects. Small worldness is a quality found in the FNC of healthy subjects and is a feature that takes all meaningful components into account simultaneously. As noted in Wang et al. (2010), the quality of small worldness in neural networks is indicative of efficient communication. Small values suggest functional network shifts toward random networks and have previously been demonstrated in stroke and other brain pathologies, such as brain tumors (Bartolomei et al., 2006), Alzheimer's disease (de Haan et al., 2009), and epilepsy (van Dellen et al., 2009). Larger values of small worldness indicate local cliquishness among groups of nodes with long-range paths connecting node clusters (Bassett et al., 2012).

Centrality

Investigation of how spatial brain components share mutual information with each other leads to an understanding of sFNC. Spatial nodes with high centrality are important in terms of the functional organization of the brain as these nodes are found along the shortest path that connect other components to each other and hence connect component clusters throughout the brain. The centrality of a node is determined by its role in the shortest path between all pairwise nodes. The shortest path from one node to another is calculated by first converting edge values into distance values. The edge values are the mutual information values between nodes which represent the z-scored spatial components. These edge values are mapped to distance values between zero and one. Let e and d represent the MI and distance values of the ith edge, respectively. The distance value for the ith edge is calculated as

The shortest path from node p to node q is the path that has the smallest sum of distances along the path. The centrality of the ith node is representative of the number of times that node i is found on the shortest path between all other pairwise nodes. High centrality suggests that the component is important in terms of efficiency of the brain's functional network connectivity. The centrality for the ith node is calculated as in Sporns et al. (2007),where E is the number of shortest paths between pairwise nodes that include node i and E is the number of shortest paths from node p to node q which is equal to one when using distance edges.

Task correlation

Two evenly sized groups were formed from the 10 subjects. Half of these subjects exhibited motor improvement with composite clinical scores ranging from 2.9 to 11.5 and will be referred to herein as responders. The other five subjects had scores ranging from − 0.4 to 1.8 and will be referred to as non-responders. In order to determine whether the fMRI data exhibit a similar trend, and hence correlate with clinical results, the two groups were separated and a time course correlation analysis was performed. The block design of 30 s rest followed by 24 s task for 8 cycles forms the hand-task time course where resting times are represented by zeros and task times are represented by ones. The hand-task time course is convolved with the hemodynamic-response function generating the model time course for task activation. The sensorimotor and cerebellum components are task-activated components that correlate with the task paradigm. The correlation coefficient is calculated between these components' time courses and the model time course for each subject. Then the average value is calculated for the responders and non-responders at pre- and post-intervention.

Results

t-Maps

Example t-maps are shown in Fig. 1 for pre- and post-intervention. Improvements can be seen in terms of larger t-values and more activated voxels for the visual, frontal, and DMN components. It is important to note that not all component t-maps revealed differences between pre- and post-intervention. However, differences are evident in the GT features despite the lack of difference in t-maps thus making GT analysis an attractive solution for quantifying post-intervention improvements. Each t-map was thresholded at p < 0.05 and the color bars beneath the t-maps in Fig. 1 indicate the t-value range.

Fig. 1

Example t-maps for some components generated by IVA. Post-intervention shows higher t-values and more activated voxels in several components such as the default mode network (DMN). The t-maps are thresholded at p < 0.05.

The average values of small worldness were calculated at a range of thresholds, i.e., edge values between nodes, for the components generated by IVA. See Fig. 2a and b for temporal and spatial small worldness, respectively. Blue circles indicate pre-intervention values and red asterisks indicate post-intervention values. Triangles on the x-axis indicate that the difference was found to be statistically significant at p < 0.01. An increase in small worldness was identified at post-intervention compared to pre-intervention, suggesting improved local efficiency in neural communication. As the threshold increased, trends in small-worldness values suggested an increase in local cliquishness among node clusters, for both temporal and spatial components. Validation methods as described in Section 3.3 revealed a false discovery rate (FDR) corrected p-value of 0.0094 for Fig. 2a and 0.0059 for Fig. 2b. The difference in small-worldness values is considered significant if the p-value at the respective threshold is at or below the FDR p-value. Note that IVA generates spatial components which are maximally statistically independent, which is why the thresholds for mutual information take on lower values than for the temporal correlations used in FNC.

Fig. 2

Small worldness vs threshold plots, where graph edges with values below the threshold are removed. A red triangle on the x-axis indicates a statistically significant increase in small worldness. FDR corrected p-values are 0.0094 for (a) and 0.0059 for (b).

Node centrality was calculated for 20 randomly selected subgroups selected from the pre- and post-group for sFNC validation and a two-sample t-test was then performed across these values. A red triangle pointing right along the x-axis of the plots in Fig. 3 indicates significantly higher centrality for the post-group whereas a blue triangle pointing left indicates that the pre-group exhibited more centrality at the given threshold. In particular, an increase in centrality was noted for the sensorimotor component across the entire range of thresholds. Increased centrality was also found for the frontal, left frontal–parietal, visual, and cerebellum components. Interestingly, the task-deactivated DMN component and the right frontal–parietal component both showed decreased centrality at post-intervention.

Fig. 3

Centrality plots for components of interest. The FDR adjusted p-values for each component are below p = 0.001. Increased centrality for the sensorimotor, cerebellum, and left frontal–parietal components. A decrease in centrality is noted for the DMN and the right frontal–parietal components.

Responders and non-responders were compared using differences in cross correlation values between the time courses of the sensorimotor and cerebellar components and the time course of the hand motor task. The average correlation value improved for the sensorimotor component of the responders by 0.1 and by 0.05 for non-responders. The standard deviation decreased for both groups at post-intervention. For the cerebellar component, the correlation value increased by 0.06 for the responders, but decreased by 0.07 for non-responders, reflecting negative improvement in the clinical results of some subjects. The standard deviation also decreased for both groups for this component, providing further confidence in the results of our correlation analysis.

Discussion

Results of this study demonstrate, for the first time, the usefulness of a novel fMRI analysis in the identification of neuroplastic changes in network communication induced by rehabilitation after stroke. Data driven IVA was used to capture between and within subject variability of distinct neural networks to which graph theoretical analysis was then applied. Several of these networks demonstrated significantly greater, albeit subtle, within network activation differences as a result of the upper extremity intervention. Despite these small changes of within network activation, marked differences in neural communication efficiency were determined. First, an increased in small worldness across networks was identified, suggesting an increased overall efficiency in the transfer of neural information. Analysis of individual spatial components then revealed increased centrality within the sensorimotor and cerebellum networks, providing further evidence of improved within-network communication. Conversely, decreased centrality was observed in the DMN and right frontal-parietal components, which likely reflected a reduced reliance on cognition as motor task performance improved. Lastly, responders were separated from non-responders based on clinical outcomes. Results revealed greater correlation between both the sensorimotor and cerebellar network and hand motor task time course in responders compared to non-responders, again suggesting a more efficient use of these two networks, which mirrored the behavioral results.

The increased values of small worldness and centrality in this study indicate improved local cliquishness among groups of neural nodes as a result of the rehabilitation intervention (Bassett et al., 2012).

In light of previous findings of functional network shifts toward random networks in stroke and other brain pathologies (Bartolomei et al., 2006; de Haan et al., 2009; van Dellen et al., 2009), our results suggest a return toward healthy brain function during motor task performance. Improvements were most pronounced within sensorimotor and cerebellar network efficiency, which is consistent with results from previous intervention studies (Grefkes and Fink, 2011; Ovadia-Caro et al., 2013; Westlake and Nagarajan, 2011; Westlake et al., 2012). However, most intriguing is that these changes would not have been as evident if the IVA approach were not paired with graph theoretical analysis. This point is especially important given the subtle changes that result from rehabilitation in the chronic post stroke stage or even among higher functioning individuals immediately after stroke. Until this point, an understanding of neural changes after rehabilitation has been confounded by group analysis approaches that collapse group data, thereby limiting their ability to capture individual variability. By preserving this information within canonical neural networks, important changes in neural communication within these networks could then be observed. Unique to this study is the decrease in communication efficiency noted within the DMN and frontal–parietal cognitive networks during the hand motor task. To the best of our knowledge, the involvement of cognitive neural resources and the changes induced by a rehabilitation intervention has not been previously examined. Following stroke, individuals often demonstrate impaired cognitive function such as decreased working memory, attention deficit, and executive function disorder. Improvements in such cognitive function were recently found to parallel increases in DMN connectivity in the resting brain up to 3 months post stroke (Park et al., 2014). However, unlike the commonly observed increase in DMN connectivity found when the brain is in a resting state, a relative decrease (i.e., suppression) in connectivity is a well-established observation during cognitively engaging tasks. Although our data-driven component analysis did not model differences between rest and task-activated brain states, the pre–post intervention increase in DMN suppression is intriguing and may suggest a return toward a more normative brain state. Future studies including age-matched healthy controls and an analysis of the spatiotemporal dynamics during the transition from rest to motor task will help us better interpret this result. A decrease in network communication within the right fronto-parietal network was also identified and reflects a reduction in the attentional requirements for the hand motor task. This network is generally engaged during tasks that require shifts of spatial attention and target detection (Shulman et al., 2010). Since our training included a motor task with auditory cueing, it was not surprising that the fMRI auditory cued hand task became more automated with a reduced need for attentional resources. Evidence of reduced attentional requirements as a result of motor learning has previously been demonstrated from a behavioral perspective, but our observation provides the first evidence of the changes in neurocognitive communication that underlies motor rehabilitation in the stroke-affected brain.

Future work should consist of much larger groups; our sample size was a limitation of this study. These groups should contain control subjects and patients with differing lesion locations, clinical improvements, and treatment interventions. The primary advantage of our group fMRI analysis is that the inherent neural variability after stroke can be modeled and maintained. GT analysis can then be applied to further analyze how GT features within and between key neural networks differ between groups as a result of different interventions. Investigations that demonstrate changes in these features with respect to clinical outcomes can ultimately lead to a better understanding of neural mechanisms and rehabilitation targets to improve motor function in individuals with stroke.