Subjects with sub-acute ischemic stroke (n = 142) were recruited from the Foothills Medical Centre or Dr. Vernon Fanning Centre in Calgary Alberta, Canada. All subjects were 18 years or older at the time of assessment. Subjects were excluded for the following reasons: clinical diagnosis of stroke prior to the current one, hemorrhagic stroke, stroke affecting both sides of the brain, stroke in the brainstem, no identifiable lesion on MRI or CT, pre-existing neurological disorder (i.e., diagnosis of Parkinson's Disease, Multiple Sclerosis, etc.), orthopedic problems, neuropathy, or pain in either upper extremity, or inability to follow the instructions for the robotic assessment due to aphasia, language barriers, apraxia, or cognitive deficits. Cerebellar strokes were excluded due to software limitations that resulted in poor normalization to the Montreal Neurological Institute (MNI) template brain. Subjects were also excluded if they demonstrated motor impairments of the ipsilesional arm on clinical assessment. Subjects were required to be alert while performing the task and were excluded if fatigue was determined to be limiting participation by the therapist operating the robot. All individuals provided written informed consent prior to study participation in accordance with the Declaration of Helsinki. This study was approved by the University of Calgary Research Ethics Board.

Each subject underwent either MRI (n = 125) or noncontrast CT (n = 17) at the Foothills Medical Centre. T2-weighted fluid-attenuated inversion-recovery (FLAIR), diffusion weighted imaging (DWI), and apparent diffusion coefficient (ADC) sequences were performed according to the standard acute stroke MRI protocol at the Foothills Medical Centre. CT scans were also performed according to the standard acute stroke protocol at the Foothills Medical Centre. We specifically used CT images that were beyond the very acute stage post-stroke (mean = 2.1 days post-stroke), so we could better define the region of infarct, and only if an MRI was not done. If the region of infarct on CT was questionable or not well-defined according to a trained stroke neurologist, subjects were not included in our analyses. These radiological investigations were performed for clinical diagnostic purposes and obtained, with consent, for use in this study. Subjects with only CT scans were included in our analyses to increase our sample size and improve statistical power, as is commonly done in large VLSM studies (Karnath et al., 2009; Verdon et al., 2010; Winder et al., 2015). Those individuals where only a CT was performed typically had clear stroke symptoms and a well-defined area of infarct on CT. In our center, when there is a clear area of infarct on CT an MRI is often not done. Sperber and Karnath (2016) have argued that excluding subjects with only CT will likely inject bias into a sample in favor of smaller subcortical strokes. There is currently no consensus as to whether subjects with only CT imaging should be excluded from lesion analysis studies. Several studies have combined both CT and MR imaging (Bates et al., 2003; Karnath et al., 2004, 2011; Winder et al., 2015) while others use only MRI (Baier et al., 2010a,b; Meyer et al., 2016). We chose to include individuals with only CT. For MR imaging, either a Siemens or GE Medical Systems scanner at 1.5 or 3T was used, respectively. The in-plane resolution was consistent across all scan types at 1 mm² and slice thickness varied from 3 to 5 mm depending on the scanner used with 0 mm interslice gaps.

Lesion location was marked by J.K. and/or S.F. directly on each axial slice of the T2-weighted FLAIR or non-contrast CT images using MRIcron software (Rorden et al., 2007) (http://www.mccauslandcenter.sc.edu/mricro/mricron/) to obtain a binary volume of interest (VOI) indicating the region of damaged tissue. The DWI and ADC images were used to accurately determine the location of acute brain damage caused by ischemia. Only areas of acute ischemia (hyperintense on DWI and hypointense on ADC) were included in the VOI. FLAIR hyperintensities with no corresponding DWI hyperintensity were not included in the lesion drawing, as these likely represented small microangiopathic changes or age-related white matter changes. In instances where a DWI hyperintensity was observed with little FLAIR hyperintensity (i.e., very acute scans), we marked the area of DWI hyperintensity on the FLAIR image. An experienced stroke neurologist (J.D. or A.Y.) who was blinded to subjects' clinical and robotic scores confirmed each lesion delineation and corrected markings where appropriate. Each subject's scan and VOI were then normalized to the Montreal Neurological Institute (MNI) template brain using the clinical toolbox in SPM8 (http://www.nitrc.org/projects/clinicaltbx) (Rorden et al., 2012; Winkler and Kochunov, 2012). CT images were transformed to MNI space using the “CT Normalization” function in the Clinical Toolbox. Images were first converted from Hounsfield units to the image brightness range of the template and then normalized using SPM8's standard normalization function. These images were then converted back to Hounsfield units (Rorden et al., 2012). MR images were transformed to MNI space using the “MR Normalization” function, since high-resolution T1 images were not available. SPM8's standard normalization procedure was then used, which first determines the optimal 12-parameter (translations, rotations, zooms, and shears) affine transformation followed by estimation of non-linear deformations (Ashburner and Friston, 1999). Default SPM8 settings were used primarily (spatial smoothing of anatomical image with 8 mm full width at half maximum (FWHM), smoothing of lesion with 3 mm FWHM and 0.5 threshold, affine regularization using the International Consortium for Brain Mapping (ICBM) space template, non-linear frequency cutoff: 25 mm, number of non-linear iterations: 16, non-linear regularization: 1, 1 × 1 × 1 mm voxel size, and automatic setting of origin). If the normalized images appeared distorted, we performed the normalization steps again, setting the origin manually. FLAIR and CT specific template images (in MNI space) were used for normalization. To prevent warping of damaged tissue during the normalization process, cost function masks were used for lesioned brain areas (Brett et al., 2001). Final normalized VOIs (binarized, warped, and smoothed) were compared to the original imaging to ensure accuracy, and used in subsequent imaging analyses.

Mean robotic task scores for each parameter were compared to lesion location using voxel-based lesion symptom mapping (Bates et al., 2003). At each voxel, robotic scores between the group with damage and the group without damage were compared using a t-test. This was performed using the non-parametric mapping (NPM) software available in the MRIcron software package. For response latency and initial direction error the absolute magnitude of error was used, with higher numbers indicating increasing impairment. For peak speed ratio and path length ratio, we used the one-sided transformed values as this provided a continuous one-sided measure of impairment on these parameters. To maintain sufficient statistical power only those voxels where a minimum of five subjects had damage were analyzed, which is a common threshold in VLSM studies (Kalénine et al., 2010; Lo et al., 2010; Geva et al., 2011; Molenberghs and Sale, 2011; Herbet et al., 2014). This threshold also ensures that we don't limit our “search area” to only the most common areas of damage. To correct for multiple comparisons we used a voxel-wise false discovery rate (FDR) correction (q = 0.01). Statistical maps are presented on the T1–weighted MNI template brain available in MRIcron.

For the region of interest analysis, 150 regions were defined. Cortical and subcortical regions were defined based on the Automated Anatomical Labeling Atlas (Tzourio-Mazoyer et al., 2002) and white matter tract regions were defined based on Neuroanatomy and Tractography Laboratory atlases (http://www.natbrainlab.com) (Catani and de Schotten, 2008). Specially designed software (Niistat; http://www.nitrc.org/projects/niistat) was operated in MATLAB (MathWorks, Natick, MA). For each region where at least five subjects had damage, the proportion of damage resulting from stroke and kinesthetic performance for those subjects were entered into a general linear model. This model tested whether the proportion of damage to a given region of interest was significantly associated with impaired kinesthesia. Results were converted to z-scores for each region and family-wise error was controlled via 4000 permutations (p < 0.05). This method improves statistical power by reducing the number of statistical tests performed.

Standard clinical assessments were performed by a study physician or therapist with experience in stroke assessment. To assess proprioception of the stroke affected upper extremity the Thumb Localizing Test (TLT) was performed (Hirayama et al., 1999). During this test the subject's eyes were closed and their stroke affected arm was manipulated by the examiner and placed in space above eye level. The subject was then asked to use the thumb and forefinger of the ipsilesional hand to pinch the thumb of the hand that was being held in a static position by the therapist. Each subject was scored based on their ability to locate the thumb on a scale from 0 (no difficulty locating thumb) to 3 (unable to locate thumb).

A battery of other clinical assessments was performed within 1 day of the robotic assessment. An apraxia assessment (Zwinkels et al., 2004) was performed on individuals clinically suspected of having apraxia and they were excluded from the study. The Functional Independence Measure (FIM) (Keith et al., 1987) was used as an indicator of daily functional abilities. The Chedoke-McMaster Stroke Assessment (CMSA) Impairment Inventory for the upper extremities was completed to assess arm and hand motor function (Gowland et al., 1993). Finally, the Behavioral Inattention Test (BIT) provided an indication of visuospatial neglect in each subject (Wilson et al., 1987).

Statistical comparisons in clinical, demographic and robotic outcomes were made between left and right hemisphere lesion subjects using independent samples t-tests or chi-squared tests when appropriate. Normality of our robotic parameters was tested using Kolmogorov-Smirnov tests. Analysis of covariance was used to account for lesion volume when identifying which variables significantly affected kinesthesia post-stroke. Spearman correlations were used to determine the relationship between demographic, clinical and robotic measurements. Statistical analyses were performed using IBM SPSS Statistics version 20. Corrections for multiple comparisons were conducted using a Bonferroni correction for behavioral data, a voxel-wise false discovery rate correction for voxel-based analyses and permutation testing for region of interest analyses.

Subject demographics and clinical scores are presented in Table 1. The mean time between stroke onset and neuroimaging was 2.3 ± 2.8 days (MRI = 2.3 ± 2.8 days; CT = 2.1 ± 2.4 days) and the time between stroke and behavioral assessment was 8.9 ± 6.4 days. A total of 24 subjects were considered to have visuospatial neglect based on BIT scores of less than 130 (Halligan et al., 1991) (18 right hemisphere, 6 left hemisphere). Five subjects presented with involvement of more than one cerebral artery territory and thus are counted in all applicable categories in Table 1.

Initially, we compared various clinical and demographic parameters of individuals with left and right hemisphere lesions. No significant differences were observed between groups in terms of age, FIM score (unpaired t-tests, p > 0.05), CMSA, or thumb localizing tests for the affected limb (Mann-Whitney U-tests, p > 0.05). On average, right hemisphere stroke subjects had significantly larger lesion volume (37.4 ± 50.0 mL) compared to left hemisphere stroke subjects (20 ± 39.2 mL; unpaired t-test, p = 0.005). Right hemisphere stroke subjects also had lower BIT scores (right = 132 ± 18 vs. left = 139 ± 12; unpaired t-test, p = 0.014).

Failure on the kinesthesia task, defined as falling outside of the normal healthy control range on three or more of the eight measured parameters (Semrau et al., 2013), occurred in 56% (n = 79) of subjects. Among right hemisphere strokes, 72% failed the task whereas 37% of left hemisphere subjects failed (Fisher's Exact test: p < 0.001). Thus, subjects with right hemisphere strokes had 1.93 times the risk of failing the kinesthesia task compared to left hemisphere strokes (risk ratio = 1.93; 95% confidence interval = 1.37–2.71).

Subjects with right hemisphere stroke had poorer performance overall on the kinesthesia task (Table 3). After correcting for multiple comparisons (p < 0.006, Bonferroni corrected, n = 8, p = 0.05) lesion volume was positively correlated with the total number of kinesthesia parameters failed (r = 0.39, p < 0.001). We also found that most of the individual kinesthesia parameters were correlated with lesion volume in our sample (Table 2). Correlations between individual parameters and clinical scores are shown in Table 2. Age did not correlate with performance on any robotic parameters or clinical scores (p > 0.05). Using analysis of covariance, including lesion volume (mL) as a covariate, and comparing left hemisphere to right hemisphere stroke subjects on each robotic parameter, we found that right hemisphere stroke subjects still performed worse compared to left hemisphere stroke subjects on the kinesthesia task (Table 3).

We generated an overlap map of all individual lesion locations as an initial indicator of lesion distribution (Figure 2). There was greater overlap of lesions in the right hemisphere. Little overlap was observed in left inferior and posterior parietal areas. Since there were differences in lesion distribution between left and right hemisphere stroke subjects, subsequent analyses were performed on left and right hemisphere lesions separately. Results of the VLSM analysis for the mean values from the kinesthesia task are shown in Figure 3, Table 4. The areas of brain damage associated with increased response latency (Figure 3A) are distributed through the right inferior frontal gyrus, inferior post-central gyri, and bilateral insula and frontal parietal operculum. Areas of damage associated with poor performance on peak speed ratio (Figure 3B) involve a much smaller area centered on the right post-central and supramarginal gyri as well as the left frontal subcortical white matter. Initial direction error (Figure 3C) was associated with lesions to bilateral post-central gyri and right supramarginal gyrus, angular gyrus, superior temporal lobe and the posterior aspect of the insula. Errors in path length ratio (Figure 3D), a measure of movement amplitude perception, were associated with lesions to the right superior and middle temporal gyri, anterior insula and bilateral parietal operculum, pre-central gyrus and supramarginal gyrus. The peak speed ratio parameter did not survive the stringent false discovery rate correction and so is presented at a more liberal threshold of p < 0.01, uncorrected.

A comparison of three of the kinesthetic parameters (peak speed ratio, initial direction error, and path length ratio) are presented in Figure 5. The response latency parameter was not included in this figure since this parameter was simply the time taken to respond to the passive movement. The other three parameters measure the kinematics of the movement (speed, direction, amplitude). Impairments in all three parameters were associated with damage to the supramarginal gyrus in the right hemisphere.

The sROI analysis examined whether the proportion of damage to each region was associated with poor performance on the robotic task. Several regions were associated with poor response latency that were highly distributed throughout the MCA territory, including the right inferior frontal gyrus, insula, superior temporal gyrus, and arcuate fasciculus (Figure 4A, Table 5). Impaired peak speed ratio was not associated with proportion of damage in any region after correction for multiple comparisons (p < 0.05, 4000 permutations). These results are presented at p < 0.01, uncorrected (Figure 4B). Impaired initial direction error was significantly associated with damage to the right post-central, superior and inferior parietal, and supramarginal gyri (Figure 4C). Proportion of damage to the left calcarine sulcus was associated with initial direction errors. Interestingly, impaired path length ratio was associated with damage to the right pre-central gyrus (Figure 4D). The individual regions and z-scores are presented in Table 5. There was good agreement between the VLSM and sROI techniques that we employed.

Longer response latency was associated primarily with lesions to the right inferior frontal gyrus, bilateral insula and parietal operculum (Figures 3A, 4A). Conversely, lesions in frontal and parietal cortices were more associated with the kinematics of movement perception. Among the various functions of the insula, it has been proposed as a structure that integrates somatosensory information to create a sense of body and self-awareness (Karnath, 2005). This body awareness processing likely informs us of our current body state and any changes to that state, such as the kinesthetic stimulus in our task. The insula is likely an important structure for identifying that a change in body state has occurred (i.e., passive arm movement) while structures in frontal and parietal cortices are important for identifying the specific nature and properties of that change (i.e., speed, direction, and amplitude of movement).

Difficulty in matching the speed of movement was significantly associated with lesions in the post-central and supramarginal gyri (Figures 3B, 4B). Previous research involving single cell recordings in non-human primates has shown that neurons in Brodmann's area 2 and 5 of the parietal lobe contain cells that modulate their activity in relation to the velocity of limb movement during motor tracing tasks (Averbeck et al., 2005). The ability of subjects to perceive the speed of passively generated limb movements, in order to guide movements of the ipsilesional limb, was significantly affected by lesions to the post-central gyrus. Directional errors were also linked to lesions along the post-central gyrus, in addition to lesions in the posterior parietal cortex, superior temporal lobe and insula (Figures 3C, 4C). Along with the post-central gyrus, it has been observed that neurons in primate Broadmann's area 5 respond to specific directions of movement during reaching (Kalaska et al., 1983). Accurate perception of movement direction may be a more complicated and integrative process than perception of movement speed, and thus requires more processing from these highly integrative areas such as the posterior parietal cortex and parietal operculum.

Impairments in matching the amplitude of movement (PLR) were significantly associated with lesions to the right superior and middle temporal gyri, anterior insula, supramarginal gyrus, bilateral parietal opercula (Figure 3D) and pre-central gyrus (Figure 4D). It is interesting that the sROI analysis revealed that the amount of damage to the pre-central gyrus, and not the post-central gyrus, was a significant predictor of impairment in matching the amplitude of movement (Figure 4D). Previous functional MRI studies using passive movement stimuli in healthy controls have shown robust activations of both pre- and post-central gyri (Naito et al., 1999, 2007). Our results suggest that the amount of damage to the pre-central gyrus may be an important indicator for impaired perception of movement amplitude. With the VLSM analysis three of the kinematic parameters (peak speed Ratio, Initial Direction Error, Path Length Ratio) shared a common lesion location (right supramarginal gyrus) that was significantly associated with kinesthetic impairments (Figure 5, white area). The right supramarginal gyrus has recently been implicated in proprioceptive processing at the wrist using functional MRI in healthy controls and stroke subjects (Ben-Shabat et al., 2015).

It should be noted that our robotic task (mirror matching movements of the opposite limb without vision) is a complex behavior as it requires subjects to not only attend to the affected arm that is moved by the robot, but also execute a movement with their unaffected arm. Thus, it requires attentional, executive, and higher motor processing in addition to kinesthetic processing. These other processes may be a contributor to the distribution in lesion maps we observed. In particular, damage to the right superior and middle temporal gyri (as seen in our response latency and path length ratio parameters) and inferior parietal lobule have been associated with visuospatial neglect post-stroke (Karnath et al., 2011) and damage to left hemisphere frontal parietal areas have been associated with impaired executive function (Barbey et al., 2012).

Our analysis identified lesions in the right hemisphere that correlated with different behavioral parameters. In contrast, we were unable to find corresponding relationships for the left hemisphere. This difference between left and right hemispheres may be due to several factors. First, our behavioral results showed that right hemisphere stroke subjects failed more parameters overall on the kinesthesia task and thus provide a greater range of impairments for comparing to lesion location (Table 3). Second, there were significant differences in lesion locations and lesion volumes between left and right hemisphere subjects in our sample (Table 1, Figure 2). Critically, we limited recruitment of many left hemisphere cortical lesions due to aphasia (Bates et al., 2003), as it was important that they understood the instructions for the kinesthesia task. We suspect that unfortunately, we simply had an insufficient overlap of lesions to accurately map where kinesthetic processing occurs in the left hemisphere, or to determine whether any lateralization in kinesthetic processing exists in the present study.

Some evidence has suggested right hemisphere dominance in processing proprioceptive information in right-handers. Goble and Brown (2010) have demonstrated that healthy subjects performed better on a dynamic position matching task in the upper extremities when matching to their non-dominant arm, suggesting a right hemisphere specialization for proprioceptive feedback processing. Additionally, the right hemisphere has been proposed to play a specialized role in sensorimotor stabilization mechanisms for accurately reaching goal targets (Mutha et al., 2012). Further research is needed to determine the potential lateralization of kinesthetic processing, and its relationship to sensorimotor function.

SS is cofounder and chief scientific officer of BKIN Technologies, the company that commercializes the KINARM robotic device used in this study. The other authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.