Acute stroke model
This study involved a retrospective analysis of data acquired previously as part of a study for which all procedures were approved by the Subcommittee for Research Animal Care (SRAC), the Institutional Animal Care and Use Committee (IACUC) of our institution that follows the National Institutes of Health Guide for the Care and Use of Laboratory Animals. Unilateral stroke was induced in seven adult male macaques (Macacca fascicularis, 7.7 ± 1.2 kg, 6–12 years old) by obstruction of the M1 branch of the middle cerebral artery (MCA) either by injection of a small volume of cyanoacrylate creating a permanent MCAo (n = 2), or by transient insertion of a micro infusion catheter. The catheter was subsequently removed after 3 h to effect reperfusion for a transient MCAo model (n = 5). Stroke induction followed procedures previously described [31]. Briefly, animals were sedated with diazepam (1 ml), then anesthesia was induced with atropine (0.04 mg/kg, i.m.) and ketamine (10 mg/kg, i.m.) and maintained either with isoflurane (2–3 %) in a 80/20 air/oxygen mixture or with propofol (300 µg/kg per hour, i.v.) in combination with remifentanil (0.1 µg/kg per hour, i.v.). Animals were mechanically ventilated with 20 % oxygen in air mixture to maintain end-tidal CO2 between 30 and 40 mm Hg, and physiological signs were monitored continuously.
Three 3-h transient MCAo group animals were excluded from analysis because of incomplete occlusion (N = 1), extensive imaging artifacts (N = 1), and MRI scanner problems in the hyperacute stage after stroke onset (N = 1).
MRI
Serial MRI data were acquired on a 1.5 T MRI scanner (GE Signa) at 1 h intervals up to 6 h post-MCAo, and at 1, 3, 7, 10, 17, and 30 days. DTI (monopolar single-shot echo planar imaging, repetition time (TR)/echo time (TE) 8400/65.9 ms, 128 × 128 matrix, number of scans (NEX) 1, b-values of 0 (b0) and 1000 s/mm2 (6 directions), field-of-view (FOV) 200 mm, 3 mm slice thickness), and dual echo T2-weighted MRI (fast spin echo, TR/TE1/TE2 4200/10.7/95.9 ms, echo train length 16, 512 × 512, NEX 2, FOV 200 mm, 2 mm slice thickness, 0.5 mm slice gap) were obtained. The apparent diffusion coefficient (ADC) and fractional anisotropy (FA) maps were calculated from eddy-current corrected DTI datasets (MRVision, Winchester, MA, USA); T2 maps were calculated from the T2-weighted images.
Data analysis
Pre-processing
Serial T2-maps and DTI were spatially aligned to the 1-h post MCAo b0 image using a two-step co-registration procedure. In the first step global alignment was achieved using a full affine, mutual information based registration procedure (Montreal Neurological Institute (MNI), Autoreg) [32]. In the second step, ventricular distortions, mid-line shifts, and sequence-induced distortions were largely compensated for using a mutual information based B-splines approach with a specialized cost function that enforced rigid transformations for particular volumes of interest (VOI) [33]. This prevented local lesion volume expansion or compression, effectively reducing potential co-registration induced artifacts.
Before cluster analysis, surrounding muscle and skull tissue was removed from the MRI (FMRIB Software Library (FSL) Brain Extraction Tool) [34]. Tissue with ADC values greater than 1.2 × 10−3 mm2/s at the first time-point, mostly arising from CSF and large vessels, were excluded [2]. The resulting segmentation served as a mask for the cluster analysis. Relative ADC (rADC), FA (rFA) and T2 (rT2) maps were calculated by dividing each map by its mean contralateral hemispheric value.
ISODATA cluster analysis
Temporal lesion evolution was assessed using a spatially and temporally adjusted ISODATA (ST-ISODATA) approach [6, 7]. ST-ISODATA determined the number of clusters (K) based on the underlying data rather than by prior specification of the expected number of clusters. A detailed description of the original algorithm can be found elsewhere [2, 9]. In brief, ISODATA is a K-means related unsupervised segmentation algorithm that automatically determines the number of clusters used for segmentation by iteratively merging and splitting clusters based on intra- and inter-cluster dispersion terms. These terms—that include the minimum number of voxels within a cluster (ΦN), the minimum inter-cluster distance (ΦC), and the maximum allowed intra-cluster dispersion (ΦS)—can be set prior to and may allow for certain control over the segmentation process. For this study values for ΦC were derived by calculating the Mahalanobis distance between contra-lesional white matter (WM) and grey matter (GM) [2, 14]. ΦS was obtained by calculating the standard deviation of WM values over all datasets within each brain. An initial guess of the number of clusters (k) was necessary to initialize the algorithm. The initial k cluster means were calculated using a semi-randomized approach where the new cluster mean (µ
j
) was chosen to be proportional to a randomly selected fraction of the total distance to all the previously selected cluster means [35]. K defined the cluster range in which splitting and merging was allowed and was defined as 0.5*K < K < (2*K + 1). At each iterative step, voxels were assigned to the clusters where the Mahalonobis distance to the centroid (i.e. cluster mean) was smallest. Clusters with less than ΦN voxels were discarded and their voxels were redistributed over the remaining clusters based on their smallest relative distance. Subsequently, the intra-cluster distances (Dintra: the average Mahalanobis distance between voxel vectors and the cluster centroid) and inter-cluster distances (Dinter: the Mahalanobis distance between two cluster centroids (µ
i
, µ
j
)) were calculated. When Dintra was more than ΦS or Dinter was less than ΦC, clusters were split or merged respectively. Subsequently, the cluster means were recalculated and a new iteration was initiated provided the stopping criteria had not been met. Clustering was completed when the algorithm converged or a maximum number of iterations (I) was reached. The stopping and merging criteria were based on maximum number of iterations (I), convergence error threshold (ɛr), and number of merging steps per iteration (L). For this study ΦN,k,K I, ɛr, and L were empirically derived and set to: 100, 6, 8, 100, 0.0001, and 1, respectively, to ensure convergence of the algorithm.
Additional spatial contiguity constraints were previously suggested to reduce the influence of noise and outliers on the ISODATA clustering results [2]. In this study, ISODATA was extended with spatial contiguity constraints that for each voxel calculated the local neighborhood intensity homogeneity weighted by the distance of the neighboring voxels to the voxel in the center prior to ISODATA clustering [36]. For each animal’s MRI dataset, a t x f feature matrix was created for each voxel (v)—for which t represents the number of time points and f the number of features per time point—and evaluated using ISODATA. To further reduce temporal artifacts particularly at tissue boundaries, the resultant ISODATA maps were pruned according to the amount of temporal dispersion measured using coefficient of variance (CoV) analysis. CoV was employed to reduce the influence of small noisy fluctuations in MRI indices over time. Temporal cluster dispersion was measured as the ratio between its temporal standard deviation and its mean. A cluster with dispersion above a pre-defined threshold was considered abnormal; clusters with less dispersion were merged into the normal tissue cluster.
Because the number of clusters identified may depend on the data used for clustering and the order in which the cluster values are assigned by ISODATA may vary per segmentation, the pruned maps were normalized using previously described techniques [16]. Pruned maps were scaled by assigning the cluster matching the cerebrospinal fluid (CSF) region to 100 and the cluster matching the contralateral WM region to 1, the remaining clusters identified by the algorithm were scaled to values between 1 and 100 depending on the difference of each cluster’s mean centroid and those of the CSF or contralateral WM regions [16].
For group comparisons, global group evolution profiles or signatures were obtained by binning in groups of 10 the normalized signatures within the ipsilateral hemisphere. Values below signature 5 were classified as a signature consistent with “normal”. Remaining signatures were binned in steps of 10 and represented regions of abnormal tissue signatures.
Statistical analysis
Region-of-interests (ROIs) were manually outlined by two experienced researchers blinded to the ISODATA results. These manually delineated lesion volumes were enlarged with three 3-dimensional dilation steps. Four tissue classes were operationally defined as: ‘Core’, ‘Growth’, ‘Recovery’, or ‘Edema’. Voxels within the 1-h ADC lesion that overlapped with the “chronic” lesion outlined on the 17-day T2 map were used to define ‘Core’ voxels. Voxels that did not overlap were subdivided in voxels that were within the “chronic” lesion, but not within the ‘Core’. These were defined as lesion Growth voxels. Voxels included in the ‘Core’, but not in the “chronic” lesion, were marked as lesion Recovery regions. Areas of Edema were defined as voxels that appeared abnormal at the time-point for which the maximum T2 lesion volume occurred (‘Maximal Lesion’), but were normal at both the 1-h and 17-day time points (Fig. 7). All remaining voxels were classified as normal tissue.
Five combinations of input MRI-parameters were considered: ADC, FA, T2, ADC + FA, ADC + T2, and ADC + FA + T2. Additionally, CoV thresholds were varied from 0.1 to 0.01 with steps of 0.01 to determine the threshold that resulted in highest overlap between calculated abnormal ST-ISODATA regions and the ‘Maximal Lesion’ ROI. Overlap was calculated by subdividing the clustered voxels in voxels that were correctly classified as abnormal (TP); incorrectly classified as abnormal (FP); incorrectly as classified normal (FN); or correctly classified as normal (TN) and expressed as Dice’s similarity index (DSI = (2*TP)/(2*TP + FP + FN)) [37]. DSI of 1 represented excellent overlap, DSI of 0 no overlap.
Repeated measures one-way analysis of variance (ANOVA) with post hoc Tukey correction was used to compare DSI of the ST-ISODATA models using the different MRI input parameters or CoV thresholds. Tissue class volumes were compared using Chi Square test. Repeated-measures ANOVA with post hoc Tukey correction was used to evaluate differences between temporal signatures. Values are reported as mean ± SD. P < 0.05 was considered statistically significant.
Histological brain tissue preparation and analysis
Animals were euthanized by an intravenous overdose of sodium pentobarbital at 30 days post-stroke induction with the exception of one animal that was sacrificed earlier (17 days) due to poor prognosis. Brains were removed and fixed with 10 % formalin followed by gross sectioning of the brains into 2.5 mm blocks, to match the MRI slice thickness and orientation. Coronal blocks were embedded in paraffin and consecutively sectioned into 6 µm thick slices from the cut face throughout the entire lesion area. Successive slices were stained with hematoxylin and eosin (H&E), myelin basic protein (MBP), and luxol fast blue (LFB), scanned (Epson® Perfection 3170 Photo Scanner Epson America, Miami FL, USA), and pictures digitally stored. A template of the brain showing the boundaries of the affected regions and the outline of the brain was manually traced on the pictured sections. The stained sections were examined and rated by an experienced primate neuropathologist. Histopathological tissue features within abnormal MRI tissue signatures were assessed.