Automated Ischemic Lesion Segmentation in MRI Mouse Brain Data after Transient Middle Cerebral Artery Occlusion

Introduction

In pre-clinical research on ischemic stroke, histological evaluation of brain infarct volume is accepted as the gold standard. However, errors introduced due to changes in brain morphology during processing of brain sections (swelling/shrinkage of tissue) combined with its manual-labor-intensive nature make it a suboptimal evaluation method. Moreover, animal sacrifice makes longitudinal studies and multiple readout times impossible, or, alternatively, will considerably increase the number of animals used in case parallel groups are investigated. Hence, magnetic resonance imaging (MRI) has become increasingly important to assess infarct volume development in animal experiments. With MRI, ischemic stroke can be analyzed in the first few days after induction using spin-spin relaxation time contrast T2 (Jacobs et al., 2000), which is sensitive to vasogenic edema (Dijkhuizen and Nicolay, 2003). Using MRI, infarct characteristics such as infarct volume and progression can be examined in a longitudinal manner (Hoehn-Berlage et al., 1995a,b; Weber et al., 2006). Efficient use of computational image processing methods allows for automated reproducible quantitative analysis of such data.

In clinical research, several algorithms for detection, segmentation and classification of different brain abnormalities (brain tumors, trauma lesions, hematoma, edemas, Alzheimer's disease) in both MRI and CT have been developed (Ghosh et al., 2014). However, segmentation of brain lesion in MRI animal (especially mouse) imaging data still heavily relies on manual tracing and semi-automated thresholding techniques (Rekik et al., 2012; Leithner et al., 2015; Yao et al., 2015; Yin et al., 2015; Zheng et al., 2015; Donath et al., 2016; Moraga et al., 2016; Yip et al., 2016), making the analysis time-consuming, with low inter- and intra-observer reproducibility (Niimi et al., 2007; Ghosh et al., 2011; Rekik et al., 2012). Also, there is considerable variability between research centers with respect to the methods used to calculate lesion volume using MRI (Milidonis et al., 2015), which can result in unjustified conclusions.

When it comes to ischemic brain tissue segmentation of pre-clinical data, only a few (semi-) automated algorithms have been developed. These algorithms typically use data from multiparametric MR imaging, more specifically by combining the apparent diffusion coefficient (ADC) maps with T1, T2, T1- and T2-weighted images to distinguish and characterize healthy and ischemic brain tissues (Hoehn-Berlage et al., 1995a; Li et al., 2002; Sotak, 2002; van Dorsten et al., 2002; Wang et al., 2007). Ghosh et al. (2011, 2012, 2014) developed an automated method—Hierarchical Region Splitting (HRS)—where adaptive thresholds are automatically selected to detect, quantify and distinguish between core and penumbral tissue regions in T2-weighted MRI data of HII Sprague-Dawley rats. Jacobs et al. (1999a,b, 2000, 2001a,b) developed an unsupervised segmentation algorithm based on K-means clustering—iterative self-organizing data analysis technique (ISODATA)—for analysis of multiparametric MRI data of focal cerebral ischemia in Wistar rats that was validated with histology data. Ding et al. (2004, 2006) applied a modified version of the same method to analyze embolic stroke rat data. A comprehensive overview of all available image analysis methods for ischemic stroke lesion and a discussion of their pros and cons was provided by Rekik et al. (2012).

Most of the aforementioned methods were developed for rats, are age- and disease-specific and require manual intervention. The only two automated segmentation methods developed for pre-clinical data—HRS (Ghosh et al., 2011, 2012, 2014) and ISODATA (Jacobs et al., 1999a,b, 2000, 2001a,b)—obtained promising results but specifically focused on neonatal and adult rat and were validated on relatively small data sets.

Moreover, current pre-clinical stroke studies are compromised by reproducibility issues (Dirnagl, 2016; Llovera and Liesz, 2016). Having a reproducible automated method available might be of particular value when considering published guidelines for reporting animal research studies (Kilkenny et al., 2010). To the best of our knowledge, there are no automated methods that have been developed to segment ischemic lesions in MRI mouse brain data.

Thus, the goal of this study was to develop an automated approach to quantify ischemic lesion volumes in MRI data of mouse brains and to make it publicly available. Our algorithm is built upon existing segmentation paradigm (level sets) with a minimal set of assumptions and requires only one MRI sequence (T2) as input. We validated our algorithm on heterogeneous data consisting of 121 mouse brain scans, that covered various age groups, time points after infarct induction and were obtained with different MRI hardware and acquisition parameters. Performance of our automated approach was compared against the ground truth obtained in a semi-automated fashion by two observers, evaluating the volumetric accuracy and regional overlap of segmented ischemic lesions. Our approach showed good agreement with results obtained by the human observers and high accuracy on heterogeneous data, whilst having much shorter execution times. Thus, it has considerable potential in replacing the biased manual labor in quantifying ischemic lesion volume in mouse brain MRI data and might be of value for MRI stroke volume analysis in a multicenter setting.

Materials and Methods

Animals

Our computational approach was tested on three main groups of male mice that were labeled as “Leiden-Set,” “Cologne-Set-1” and “Cologne-Set-2” depending on the origin of the data and data acquisition protocol.

Leiden-Set: C57BL/6J mice (n = 65) were further subdivided into three age groups (3- to 5-, 12- to 14- and 20- to 24-month-old). Mice were repeatedly scanned at different time points: 4 h, 24 h, 48 h and 8 d after infarct induction (see Section “Experimental Infarct Model”). Histological specimens were obtained on a small subset of animals (n = 6) and used for visual validation of the infarct area.

Cologne-Set-1: C57BL/6J mice (n = 6; 3-month-old) were scanned at 18 h and 4 d after infarct induction.

Cologne-Set-2: Transgenic mice expressing luciferase under doublecortin control (DCX-Luc, Couillard-Despres et al., 2008) (n = 10; 2- and 12-month-old) were scanned at 48 h after infarct induction.

The group of 3- to 5-month-old mice from the Leiden-Set, 24 h after stroke induction, was used as a primary validation cohort as it was the largest group (Table 1) and because the 24-h time point is often chosen in cerebral infarction experiments. For better assessment of the performance of our method on different infarct shapes and volumes, we subdivided this cohort into striatal (Str) and corticostriatal (Ctx+Str) lesions.

TABLE 1

Table 1. Overview of all the used mice and accompanying MRI data sets: origin/acquisition protocol, number, age, time points after experimental infarct induction.

Animals were housed with littermates, in a temperature-controlled environment, with food and water ad libitum. All animal experiments performed at the Leiden University Medical Center (LUMC) were approved by the local committee for animal health, ethics and research of LUMC. All animal experiments conducted at the Max Plank Institute for Metabolism Research in Cologne were performed in accordance with the German Animal Welfare Act and approved by the local authorities (Landesamt für Naturschutz, Umwelt und Verbraucherschutz NRW).

Experimental Infarct Model

Infarcts were induced using a modified transient middle cerebral artery occlusion (tMCAo) model first described by Longa et al. (1989). Mice were anesthetized using isoflurane (3% induction, 1.5% maintenance) in 70% pressurized air and 30% O₂. Painkiller carprofen (5 mg/kg, s.c.; Carporal, 50 mg/mL, AST Farma BV, Oudewater, the Netherlands) was given before surgery. During surgery, the mouse body temperature was maintained at 37°C using a rectal probe and feedback system. During the surgical procedure, a silicone-coated nylon monofilament (7017PK5Re; Doccol Company, Redlands, CA, USA) was inserted into the right common carotid artery and advanced via the internal carotid artery and circle of Willis to eventually block the middle cerebral artery (MCA) at its origin (decreasing blood flow substantially in the MCA territory, in the right hemisphere) and the skin was sutured. During the occlusion period, the mouse was allowed to wake up in a temperature-controlled incubator (V1200; Peco Services Ltd, Brough, UK). After 30 min of occlusion, the mouse was re-anesthetized in order to remove the suture and withdraw the monofilament to allow reperfusion. After surgery, the animal was allowed to recover for 2 h in the incubator to maintain body temperature at 37°C, with easy access to food and water. At the end of the experiment, brains were fixated by transcardially perfusing the mice using fresh cold 4% buffered PFA (Paraformaldehyde P6148; Sigma-Aldrich Co. LLC, Saint Louis, MO, USA). Brains were collected, processed and sectioned. Sections were stained with Nissl (Cresyl Violet, Merck Millipore, Billerica, MA, USA) using standard protocol.

Magnetic Resonance Imaging

Scans were acquired with small-animal Bruker MRI systems using a Multi-Slice Multi-Echo sequence protocol at different time points after infarct induction. Animals from the Leiden-Set were scanned at 7 T (Pharmascan, Bruker BioSpin, Ettlingen, Germany), whilst animals from the Cologne-Sets were scanned at 11.7 T (Biospec 11.7 T/16, Bruker BioSpin). Quantitative T2 maps were calculated from the multi-echo trains using Paravision 5.1 software (Bruker Pharmascan) for the Leiden-Set and IDL software for the Cologne-Sets. Table 1 shows a complete overview of all 121 scans, together with a summary of the main imaging acquisition parameters.

Ischemic Lesion Segmentation Challenges

The ischemic lesion is characterized by elevated T2 values with respect to the healthy brain tissue. However, it is not possible to segment the ischemic lesion with a simple threshold as several other objects, like ventricles, structures surrounding the brain and other regions in the periventricular zone of the brain, share similar T2 values to that of the infarct lesion. Also, both the ventricles and the stroke regions vary considerably in shape and size between different subjects and different time points after infarct induction, and in many cases the stroke region engulfs the ventricles making the segmentation task even more difficult. Manual delineation of the ischemic lesion by experts takes into account all this information, together with possible stroke density differences and contiguity of the lesion area throughout the MRI image stack.

Semi-Automated Ischemic Lesion Segmentation

Semi-automated segmentation was performed in a slice-by-slice manner by two trained observers (IM and SdJ) using freely available ImageJ software (https://imagej.nih.gov/ij/); see Figure 1.

M1. Threshold determination. The threshold, the same for both observers, was determined as the mean plus two standard deviations of a vector containing average T2 values within a ROI in the contralateral hemisphere for every group of animals (concerning time point, age and data set origin) (Figure 1M1).

M2. Threshold mask. A mask was compiled of all high-intensity pixels above the threshold value (Figure 1M2).

M3. Infarct ROI segmentation. Confounding regions (ventricles and other non-infarct-related high-intensity areas) were manually removed by both observers individually (thereby highly relying on anatomical knowledge, left-right symmetry, density and experience), resulting in the delineated infarct lesion (Figure 1M3).

M4. Infarct ROI volume quantification. For each observer, the number of voxels located inside the ROIs representing the infarcted area of each MRI slice was multiplied by the voxel size to obtain the total infarct lesion volume (Figure 1M4).

FIGURE 1

Figure 1. The pipeline for the animal model, data acquisition and lesion segmentation. The transient middle cerebral artery occlusion (tMCAo) model was used for infarct induction in mice of different age groups. A Multi-Slice Multi-Echo MRI scan was acquired for determination of quantitative T2 maps at multiple time points (T = 4 h, 18 h, 24 h, 48 h, 4 d and 8 d) after tMCAo. A subset of infarcted brains was used for histological validation. All MRI scans were analyzed in a semi-automated fashion by two observers (IM and SdJ) using threshold determination (M1), threshold mask (M2), manual exclusion of non-infarct-related areas (M3) and lesion volume determination (M4). Automated lesion segmentation was performed using image registration (A1), brain segmentation (A2.a), contralateral ventricle segmentation (A2.b), ventricles segmentation (A2.c), infarct ROI segmentation (A3) and infarct ROI volume quantification (A4).

Automated Ischemic Lesion Segmentation

The infarct lesion was segmented on the T2 map based on its elevated T2 values with respect to the healthy brain tissue. However, as mentioned above, several other objects appearing on the T2 map share similar T2 values to that of the infarct lesion. To discriminate the infarcted area from other brain structures and to quantify the infarct lesion volume, we developed a fully automated approach that integrates a series of image processing steps depicted in Figure 1.

Our approach is based on two general assumptions that: (i) T2 value distribution of the infarct lesion differs from that of the healthy tissue; and (ii) only the ipsilateral hemisphere is affected by the infarct, whereas the contralateral hemisphere remains unaffected. Thus, we use the latter to learn the T2 value distributions of the ventricles and of the healthy brain tissue and apply this information to segment the infarcted area. All segmentations were performed on 3D image volumes using a modification of the region-based level sets method of Chan and Vese (2001). The parameters of our method (listed in Table 2) were optimized on the Leiden-Set and fixed for all three data sets during validation. The remainder of this section provides a detailed description of each particular step.

TABLE 2

Table 2. Parameters of the level-set-based segmentation algorithm for different segmentation steps.

Image Registration

In this step, each brain scan was registered to a template brain consisting of a number of manually drawn labels: whole-brain (M_WB), ipsilateral hemisphere (M_IBH), contralateral hemisphere (M_CBH), ipsilateral ventricle (M_IV), contralateral ventricle (M_CV) and periventricular zone (M_PVZ); see Figure 1A1. The template labels were propagated to each subject and used to initialize the subsequent steps. In the following, by referring to a region being occupied by a certain label we mean the result of the label propagation.

For each subject, the sum of all its echo images was used to register the scan of that particular subject to the reference brain scan. Consequently, the template labels were propagated to the individual data sets using the information provided by the deformation field for each subject-to-reference registration. The labels were used to initialize the segmentation of the whole brain and the ventricles as described in the next section. The quality and success of the registration was visually inspected by three independent observers (AK, OD and IM).

Registration was performed in a coarse-to-fine fashion. Initially, rigid registration was performed to compensate for translation and rotation. Afterwards, affine registration was conducted to compensate for differences in brain size. Because large deformations occur in stroke brains, a non-rigid B-spline registration was necessary to compensate for the large local changes (especially in the ipsilateral hemisphere and the ventricles region). A Gaussian image pyramid was employed in all registration steps, applying four resolutions for the rigid and two for the affine and B-spline registrations each. Normalized Correlation Coefficient was used as a similarity metric.

Segmentation

For each particular segmentation, we used a level set function ϕ(x) to partition the image space Ω into two classes, further referred to as “object” (Ω_O = {x ∈ Ω : ϕ(x) ≥ 0}) and “background” (Ω_B = {x ∈ Ω : ϕ(x) < 0}), respectively. Here x ∈ Ω ⊂ ℝ³ are the Cartesian coordinates. The level set function ϕ(x) was evolved from its initial state ϕ₀(x) for the predefined number of iterations n_iter or till the stopping criterion was reached. The energy functional E(ϕ) contains both image-based and regularization-based terms and, in its most general form, is given by the following equation (Rousson and Deriche, 2002):

$\begin{array}{l}E\left(\varphi \right)=\alpha ·g(\text{x})·\mathrm{\text{Length }}\left(\partial {\Omega }_O\right)+\mu ·\text{Area }\left({\Omega }_O\right)\\ -logP\left(I(\text{x});{\Omega }_O\right)-logP\left(I(\text{x});{\Omega }_B\right).\end{array}$

Here g(x) is the gradient map (Dufour et al., 2005) that can be optionally used to drive the segmentation toward the boundaries of the structures of interest, and α and μ are weights. Parameter values for different segmentation steps described in the remainder of this section are summarized in Table 2.

Data preprocessing

After the label propagation was completed, each scan was cropped to the rectangle containing the volume of interest. The cropping rectangle was defined as the smallest rectangle that contains the brain mask M_WB obtained as result of the label propagation. The volume intensities were consequently scaled to the [0;1] range. For the echo images we performed an additional per-slice intensity normalization by mapping the cumulative intensity distribution of each slice to that of the chosen reference slice (the one with the maximum entropy) and rescaling the image intensity accordingly.

Whole brain segmentation

In this step, the segmentation was performed by minimizing the energy functional E(ϕ) on the fourth echo image; see Figure 1A2.a. The level set function in this case was initialized by the whole-brain mask M_WB obtained as result of label propagation. The final result, denoted as R_WB, was obtained by: (1) Applying the morphological hole filling and morphological opening with a disk of radius of two pixels as the structure element on the result of the level set propagation; and (2) Selecting the largest connected component. The parameter values for this step are provided in Table 2.

Contralateral ventricle segmentation

Contralateral ventricle R_CV was segmented from the T2 map by running the level sets starting from the result of the label propagation; see Figure 1A2.b. In this case, the level set evolution was restricted to the contralateral hemisphere R_CBH = R_WB\M_IBH, the rest of the parameter values are listed in Table 2.

Ventricles segmentation

Ventricles R_V were segmented on the T2 map by evolving the level sets inside R_WB; see Figure 1A2.c. The object (ventricle) energy was calculated from the segmented contralateral ventricle and kept fixed. To prevent the segmentation from leaking into the stroke area touching the ventricles, we used the gradient map that was defined as g$(\text{x})=1-e_v^2(\text{x})$, where e_v(x) = −log P_Gaussian(I(x); R_CV) is the energy (assuming the Gaussian intensity distribution) scaled to the [0;1] interval, and parameters of the Gaussian distribution were calculated from R_CV (the region occupied by the contralateral ventricle).

Stroke segmentation

Finally, the stroke area was segmented from the T2 map; see Figure 1A3. To initialize the level set function, we initially calculated the area where the histogram-based energy of the brain region R_WB with the ventricles and the periventricular area excluded (R_WB\(R_V ∪ M_PVZ)) was larger than that of the contralateral hemisphere (R_CBH\(R_V ∪ M_PVZ)):

$\begin{array}{l}{R}_s^{init}=P\left(I\left(\text{x}\right);R_B\backslash (R_V\cup M_{PVZ})\right)\\ >P\left(I(\text{x});R_{CBH}\backslash (R_V\cup M_{PVZ})\right).\end{array}$

Consequently, connected components that did not intersect with dense areas on $R_s^{init}$ were filtered out. The dense areas were defined as binary masks composed of all voxels in $R_s^{init}$ for which at least 75% of their neighbors (the neighborhood was in this case defined as the circle of radius 4 around the voxel of interest) belong to $R_s^{init}$. In this case, evolution of the level set function was restricted to R_IBH\(R_V ∪ M_PVZ) (the ipsilateral hemisphere without the ventricles and the periventricular area). The energy of the background was assumed fixed and equal to that calculated from the contralateral hemisphere R_CBH. The rest of the parameters are provided in Table 2.

Parameter selection

The parameters (α, μ, n_iter) for the level-set-based segmentation were optimized on the Leiden-Set. Suitable parameter values for segmenting the whole brain region, the contralateral ventricle and both ventricles were determined by trial and error. The parameters for the stroke segmentation were determined by maximizing the Dice index (Dice, 1945) via exhaustive search within an empirically selected range of feasible values for each parameter. The final values used for validation of our method, the same for all three data sets, are reported in Table 2.

Implementation

All the data are publically available and published as a data report (Mulder et al., in review). The template labels were manually drawn based on the Allen Brain Atlas (Sunkin et al., 2013; http://www.brain-map.org/) using AMIRA (v5, FEI Software, Hillsboro, OR, USA). Both registration and segmentation were implemented in MATLAB R2012b (The MathWorks, Inc., Natick, MA, USA) and compiled into a toolbox that can be downloaded from the following webpage (www.lkeb.nl). The registration scheme was implemented using the open-source image registration toolbox elastix (v4.700; Klein et al., 2010). Information on the used registration parameters can be found on the elastix website (http://elastix.bigr.nl/wiki/index.php/Par0038). Average computational time for executing the entire segmentation routine on one data set on 3.60 GHz Intel(R) Xeon(R) computer with 32 GB RAM was: 4 s for segmentation and 309 s for prior registration and label propagation.

Performance Measures and Statistical Analysis

Ischemic lesion regions segmented by the observers and the proposed automated approach were compared using two primary measures: total volume of the ischemic lesion and Dice index (Dice, 1945) that measures overlap between two regions R₁ and R₂:

$\begin{array}{l}Diceindex(R_1,R_2)=2\frac{|R_1\cap R_2|}{\left|R_1\right|+\left|R_2\right|}.\end{array}$

Intraclass Correlation Coefficient (ICC) and the accompanying trend line were also calculated for each group of interest. ICC (2-way mixed with absolute agreement) was determined with SPSS (SPSS Statistics 23; IBM Corp., Armonk, NY, USA) to investigate the correlation between the mean of the two observers and the calculated automated volumes.

Results

Sample segmentation results on the data sets of different origin and acquired at different time points after infarction induction are shown in Figure 2. Supplementary Material 1 provides corresponding segmentation results on all 121 data sets. Figure 3 illustrates segmented striatal (Str) and corticostriatal (Ctx+Str) lesions, as well as those with fragmented infarct areas and with large edema and morphed ventricles. The results show good agreement between both segmentation methods as well as with the accompanying histological sections.

FIGURE 2

Figure 2. Examples of segmentation results. Data sets of different origin [Leiden-Set (A–D), Cologne-Set (E)] were acquired at different time points after infarction induction [4 h (A), 24 h (B), 48 h (C,E) and 8 d (D)]. For each of the image panels, the first row illustrates reformatted raw image stack and the rest of the rows provide overlaid segmentation results by: automated (green) vs. Observer 1 (red), automated (green) vs. Observer 2 (red), Observer 1 (green) vs. Observer 2 (red), respectively. The regions where two segmentations overlap are colored in yellow. Image intensity was enhanced for visualization purposes.

FIGURE 3

Figure 3. Example traces of different infarct anatomies in a single slice and complete 3D brain. Lesions (striatal, corticostriatal, fragmented and with large edema compressing the ventricle) were segmented using the semi-automated (red delineated area) and automated (green delineated area) approaches. The fifth column shows the overlap between both methods, in yellow. The last column shows the overlay of both approaches on the corresponding Nissl stained histological section. “A” stands for anterior. V_Obs1 corresponds to the volume obtained by performing semi-automated stroke segmentation by Observer 1 (IM). V_Aut corresponds to the volume obtained by the automated method. The Dice index between Observer 1 and Observer 2 for the same data sets was: 0.82 for striatal infarct, 0.91 for corticostriatal infarct, 0.92 for fragmented infarct and 0.88 for the large edema case.

The results of experiments for all age groups at the various time points after stroke induction from Leiden-Set and Cologne-Sets in terms of the absolute volume difference, ICC and Dice index are reported in Table 3. The mean per-group absolute volume difference in the Leiden-Set ranged from 5.9 to 18.4 mm³, for the Cologne-Set-1 it was 15.3 mm³ (for 18 h) and 14.2 mm³ (for 4 d) and for the Cologne-Set-2 it was 22.0 mm³. Infarct volumes of all groups, from both observers and by the automated method are presented in Figure 4.

TABLE 3

Table 3. Comparison of performance for infarct lesions for different groups of interest in terms of absolute volume difference, ICC and Dice index.

FIGURE 4

Figure 4. Comparison of total stroke volume between our automated method and human observers, for all test data sets, different times after stoke induction and different mouse ages.

The mean per-group Dice index for Observer 1 vs. Observer 2 was between 0.79 and 0.98, whereas for automated vs. observers it was ranging from 0.63 to 0.88 (Observer 1) and from 0.63 to 0.89 (Observer 2).

Figure 5A illustrates the mean infarct volume obtained by the two observers vs. the corresponding automated volume for our primary group, for both Str and Ctx+Str lesions. This plot shows very high degree of correlation with an ICC of 0.957 (ICC = 0.966 for Str lesions and ICC = 0.944 for Ctx+Str lesions). Figures 5B–E show the corresponding information for all the other groups. At 24 h and 48 h, the data match the primary validation group. However, at 4 h and 8 d after stroke induction the regression line is less perfect: ICC = 0.806 at 4 h and ICC = 0.491 at 8 d for 3- to 5-month-old mice, ICC = 0.384 at 4 h for 12- to 14-month-old mice and ICC = 0.445 at 4 h for 20- to 24-month-old mice. Figure 5E shows that our method is also robust with respect to data obtained from different imaging hardware using different acquisition parameters, time after stroke induction and age: ICC = 0.793 at 18 h, ICC = 0.964 at 4 days and ICC = 0.552 for 1- to 2-year-old mice scanned at 48 h after stroke induction. Figure 6 shows distribution of Dice index for different time points. Here, again, the results of the 24 h and 48 h groups are superior to those of the 4 h and 8 d groups: mean Dice index of 0.68 ÷ 0.87 at 24 h and 48 h and of 0.63 ÷ 0.80 for 4 h and 8 d for all groups, respectively. Finally, Figure 7 provides segmentation results on the data sets that exhibit large disagreement between two observers. These data sets are outliers of the corresponding box-whisker plot on Figure 6.

FIGURE 5

Figure 5. Correlation of lesion volumes calculated in a semi-automated fashion by observers vs. the automated method. (A) 3- to 5-month-old mice, at 24 h after tMCAo for striatal (Str) and corticostriatal (Ctx+Str) lesions. (B) 3- to 5-month-old mice, MRI at 4 h, 48 h and 8 d after tMCAo. (C) 12- to 14-month-old mice, MRI at 4 h, 24 h and 48 h after tMCAo. (D) 20- to 24-month-old mice, MRI at 4 h, 24 h and 48 h after tMCAo. (E) MRI using different hardware and software at 18 h and 4 d (Cologne-Set-1) and 48 h (Cologne-Set-2) after tMCAo. The solid lines provide the trend of the corresponding data group and the dashed gray line corresponds to equal values on both axes and is shown as a reference.

FIGURE 6

Figure 6. Distribution of Dice index for each of the analyzed time points after stroke induction and for all three data sets: Leiden-Set, Cologne-Set-1 (18 h, 4 d) and Cologne-Set-2 (48 h).

FIGURE 7

Figure 7. Examples of segmentation results that exhibit large disagreement between two observers. Data sets were acquired at different time points after infarct induction [4 h (A,B), 24 h (C), 48 h (D,E) and 8 d (F)]. For each of the image panels, the first row illustrates reformatted raw image stack and the rest of the rows provide overlaid segmentation results by: automated (green) vs. Observer 1 (red), automated (green) vs. Observer 2 (red), Observer 1 (green) vs. Observer 2 (red), respectively. The regions where two segmentations overlap are colored in yellow. Image intensity was enhanced for visualization purposes.

We have also performed a parameter sensitivity study to assess stability of our algorithm with respect to the parameters. For this, we changed, in turn, the value of each of the parameters listed in Table 2 by ± 50% and recalculated the results. Our approach turned out to be highly robust with respect to both validation metrics: maximal change for the mean absolute volume difference was less than 6%, and for the mean Dice index it was less than 1.2%.

Discussion

The major challenges in stroke segmentation in pre-clinical MRI are (i) irregular lesion volume and shape, (ii) low resolution compared to clinical imaging data, (iii) fluctuating contrast and (iv) considerable brain deformation and asymmetry because of the induced stroke. The presented approach addresses these challenges by estimating intensity distributions per scan from the contralateral, unaffected, hemisphere and comparing these with the affected hemisphere. With our approach, an ischemic lesion in a mouse brain can be successfully segmented from T2 MR images (Milidonis et al., 2015) in a fully automated manner.

Our approach has a number of advantages over existing methods for stroke segmentation in small animal data:

1) A main asset of our automated analysis approach is that it takes only few minutes to analyze an MRI volume (containing 16 slices) and is therefore at least an order of magnitude faster than the semi-automated method (~ 60 min). It is also more robust because of its objective character excluding potential observer bias. Thus, although histology remains the gold standard for assessment of lesion formation in pre-clinical stroke research, its downsides (that one needs to sacrifice the animal for a given time point and the fact that results are confounded by morphological changes due to the processing of fragile tissue) make MRI more desirable when one investigates the development of lesions over time. It is not unexpected that MRI is establishing itself as the main diagnostic modality for such studies, as it allows automated data processing and was shown to correspond well with histology (Milidonis et al., 2015).

2) Our method is robust with respect to different infarct shapes and volumes, even for small fragmented lesions or extremely large lesions where ventricles become deformed due to edema formation in and around the lesion. It performed well for a heterogeneous group of 121 brain scans, including different MRI scanners (noise levels, resolution), different lesion anatomies (cortical, corticostriatal, with large edema, fragmented), over different time points after infarct induction (from 4 h up to 8 d). In particular, it exhibited robust performance on data acquired from different MRI hardware—with differences in slice number and thickness—without the need for modification of the model or registration/segmentation parameter settings. We expect that our method is also suitable for other MRI lesion detection sequences, such as DWI, and perhaps even for different rodent species, with only limited model and parameter adjustments required. It is important to point out that, as we have mentioned in Section “Parameter Selection,” all parameters of our algorithm were optimized with respect to the entire Leiden-Set, which means that our method was completely blinded to the properties of each particular data set. However, in our analysis we subdivided the Leiden-Set into smaller groups, based on age and time after infarct induction, and also report results on two unseen sets of images (Cologne-Sets). This naturally results in suboptimal performance on each particular group, confirmed, e.g., by lower ICCs reported in Figure 5 and Table 3. This effect is especially pronounced in relatively small groups as these were underrepresented at the parameter tuning stage. Optimizing the parameters for each particular group will help improving the performance, which, however, will inevitably be limited by complexity of the data.

3) Unlike numerous published methods in the field (Jacobs et al., 1999a,b, 2000, 2001a,b; Bernarding et al., 2000; Li et al., 2002; Sotak, 2002; van Dorsten et al., 2002; Soltanian-Zadeh et al., 2003; Wang et al., 2007), our algorithm operates on a single contrast (T2) and does not require a pre-scan.

4) Unlike other methods developed for small animal data (HRS by Ghosh et al., 2011, 2012, 2014 and ISODATA by Jacobs et al., 1999a,b, 2000, 2001a,b), the ischemic lesion segmentation approach presented here only makes use of the quantitative T2 maps. Contrary to T2-weighted MRI data, true quantitative MRI maps (T2 or others) are comparable across research centers, whereas parameter-weighted images (T2-weighted or others) are subjective estimates derived for better discrimination of the object of interest from the background and are not comparable across research centers due to arbitrary operator choice of TE and TR values. A combination of complementary quantitative parameters, for instance, T2 and ADC values, on the other hand, can in some cases improve discrimination not only between the different objects present in an image, but also between the different subcategories of the object (Hoehn-Berlage et al., 1995a, 1997). This improvement, however, comes at the expense of slightly increased scan times and additional data analysis.

5) Our approach is generic as we do not make any specific assumptions or create a general model that would describe all of our data. Instead, the stroke area is segmented on each mouse brain MR volume by using the intensity distributions of the background and the ventricles calculated from the very same image. In all validation experiments, our method was completely blinded to the properties (origin, age, time point) of each particular brain volume, meaning that the same parameter settings were used for all data sets (including those acquired with different hardware and software compared to the training set).

The main pitfall of our approach is its multi-step nature. This type of complex algorithm, and our method in particular, is sensitive to propagation of errors made at earlier stages. Our experiments show that success or failure of the registration with label propagation, the first step in our algorithm, has significant impact on subsequent segmentation steps. More precisely, the largest segmentation errors were caused by inability to accurately segment the ventricles, e.g., when they were virtually invisible due to a large stroke area. In particular, inability to achieve high-quality registration due to large difference between the atlas scan and the rest of the data, also explains suboptimal performance of our method on Cologne-Sets.

In the absence of histology, semi-automated segmentation by experts has always been used as the reference standard in stroke quantification, which was also the case for this work. Thus, during development of any automated stroke segmentation method the main goal is to achieve performance comparable to that by experts. This would allow bypassing the manual observer bias and, hence, lead to more objective quantification of the stroke region. In this work, we evaluated the performance of our method by comparing it to two observers (separately and combined) and by analyzing the inter-observer variability.

Our results indicate that quantification of the infarct in its acute phase (4 h after induction) remains challenging, both for our automated approach and for the human observers. During this early phase of lesion development, T2 enhancement is still very weak as edema is only slowly evolving, so that in some cases it may barely reach above the normal contralateral value. It should be stated that also histological lesion demarcation at this very early time point is unreliable, which admittedly complicates infarct detection during the first hours after stroke induction in experimental animal models.

Author Contributions

Conceived and designed the study: IM, AK, OD, MH, BL, AvdM. Performed animal experiments, tMCAo surgery and acquired MRI data: IM. Performed histopathology: IM, NR. Performed semi-automated ischemic lesion segmentation: IM, SdJ. Developed the automated ischemic lesion segmentation in MRI mouse brain data after tMCAo occlusion algorithm: AK, OD. Analyzed the data: IM, AK, OD, MH, BL, AvdM. Wrote the manuscript: IM, AK, OD. Discussed the results and commented on the manuscript: IM, AK, OD, MW, MH, BL, AvdM.

Funding

The authors acknowledge funding from the: Dutch Heart Foundation (2011T055; MW), ZonMW Veni grant (MW), Dutch Brain Foundation (F2014(1)-22; MW), Centre for Medical Systems Biology (CMSB) in the framework of the Netherlands Genomics Initiative (NGI) (AvdM), FP7 EUROHEADPAIN (no. 602633; AvdM), Marie Curie IAPP Program BRAINPATH (no. 612360; AK, AvdM, MH), FP7/2007-2013 under grant agreement no. 604102—Human Brain Project (AK, AvdM, BL), H2020-Marie Skłodowska-Curie Action Research and Innovation Staff Exchange (RISE) Grant 644373-PRISAR (BL), and Dutch Technology Foundation STW (as part of the STW project 12721: “Genes in Space” under the IMAGENE perspective program; OD).

Conflict of Interest Statement

The 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.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/article/10.3389/fninf.2017.00003/full#supplementary-material

Supplementary Material 1. Segmentation results on the entire validation set.

References

Bernarding, J., Braun, J., Hohmann, J., Mansmann, U., Hoehn-Berlage, M., Stapf, C., et al. (2000). Histogram-based characterization of healthy and ischemic brain tissues using multiparametric MR imaging including apparent diffusion coefficient maps and relaxometry. Magn. Reson. Med. 43, 52–61. doi: 10.1002/(SICI)1522-2594(200001)43:1<52::AID-MRM7>3.0.CO;2-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Chan, T. F, and Vese, L. A. (2001). Active contours without edges. IEEE Trans. Image Process. 10, 266–277. doi: 10.1109/83.902291

PubMed Abstract | CrossRef Full Text | Google Scholar

Couillard-Despres, S., Finkl, R, Winner, B., Plotz, S., Wiedermann, D., Aigner, R., et al. (2008). In vivo optical imaging of neurogenesis: watching new neurons in the intact brain. Mol. Imaging 7, 28–34. doi: 10.2310/7290.2008.0004

PubMed Abstract | CrossRef Full Text | Google Scholar

Dice, L. R. (1945). Measures of the amount of ecologic association between species. Ecology 26, 297–302. doi: 10.2307/1932409

CrossRef Full Text | Google Scholar

Dijkhuizen, R. M., and Nicolay, K. (2003). Magnetic resonance imaging in experimental models of brain disorders. J. Cereb. Blood Flow Metab. 23, 1383–1402. doi: 10.1097/01.WCB.0000100341.78607.EB

PubMed Abstract | CrossRef Full Text | Google Scholar

Ding, G., Jiang, Q., Li, L., Zhang, L., Zhang, Z. G., Soltanian-Zadeh, H., et al. (2006). Characterization of cerebral tissue by MRI map ISODATA in embolic stroke in rat. Brain Res. 1084, 202–209. doi: 10.1016/j.brainres.2006.02.054

PubMed Abstract | CrossRef Full Text | Google Scholar

Ding, G., Jiang, Q., Zhang, L., Zhang, Z., Knight, R. A., Soltanian-Zadeh, H., et al. (2004). Multiparametric ISODATA analysis of embolic stroke and rt-PA intervention in rat. J. Neurol. Sci. 223, 135–143. doi: 10.1016/j.jns.2004.05.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Dirnagl, U. (2016). Thomas Willis lecture: is translational stroke research broken, and if so, how can we fix it? Stroke 47, 2148–2153. doi: 10.1161/strokeaha.116.013244

PubMed Abstract | CrossRef Full Text | Google Scholar

Donath, S., An, J., Lee, S. L., Gertz, K., Datwyler, A. L., Harms, U., et al. (2016). Interaction of ARC and Daxx: a novel endogenous target to preserve motor function and cell loss after focal brain ischemia in mice. J. Neurosci. 36, 8132–8148. doi: 10.1523/jneurosci.4428-15.2016

PubMed Abstract | CrossRef Full Text | Google Scholar

Dufour, A., Shinin, V., Tajbakhsh, S., Guillén-Aghion, N., Olivo-Marin, J. C., and Zimmer, C. (2005). Segmenting and tracking fluorescent cells in dynamic 3-D microscopy with coupled active surfaces. IEEE Trans. Image Process. 14, 1396–1410. doi: 10.1109/tip.2005.852790

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghosh, N., Recker, R., Shah, A., Bhanu, B., Ashwal, S., and Obenaus, A. (2011). Automated ischemic lesion detection in a neonatal model of hypoxic ischemic injury. J. Magn. Reson. Imaging 33, 772–781. doi: 10.1002/jmri.22488

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghosh, N., Sun, Y., Bhanu, B., Ashwal, S., and Obenaus, A. (2014). Automated detection of brain abnormalities in neonatal hypoxia ischemic injury from MR images. Med. Image Anal. 18, 1059–1069. doi: 10.1016/j.media.2014.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghosh, N., Yuan, X., Turenius, C. I., Tone, B., Ambadipudi, K., Snyder, E. Y., et al. (2012). Automated core–penumbra quantification in neonatal ischemic brain injury. J. Cereb. Blood Flow Metab. 32, 2161–2170. doi: 10.1038/jcbfm.2012.121

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoehn-Berlage, M., Eis, M., Back, T., Kohno, K., and Yamashita, K. (1995a). Changes of relaxation times (T1, T2) and apparent diffusion coefficient after permanent middle cerebral artery occlusion in the rat: temporal evolution, regional extent, and comparison with histology. Magn. Reson. Med. 34, 824–834. doi: 10.1002/mrm.1910340607

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoehn-Berlage, M., Hossmann, K. A., Busch, E., Eis, M., Schmitz, B., and Gyngell, M. L. (1997). Inhibition of nonselective cation channels reduces focal ischemic injury of rat brain. J. Cereb. Blood Flow Metab. 17, 534–542. doi: 10.1097/00004647-199705000-00007

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoehn-Berlage, M., Norris, D. G., Kohno, K., Mies, G., Leibfritz, D., and Hossmann, K. A. (1995b). Evolution of regional changes in apparent diffusion coefficient during focal ischemia of rat brain: the relationship of quantitative diffusion NMR imaging to reduction in cerebral blood flow and metabolic disturbances. J. Cereb. Blood Flow Metab. 15, 1002–1011. doi: 10.1038/jcbfm.1995.126

PubMed Abstract | CrossRef Full Text | Google Scholar

Jacobs, M. A., Knight, R. A., Soltanian-Zadeh, H., Zheng, Z. G., Goussev, A. V., Peck, D. J., et al. (2000). Unsupervised segmentation of multiparameter MRI in experimental cerebral ischemia with comparison to T2, diffusion, and ADC MRI parameters and histopathological validation. J. Magn. Reson. Imaging 11, 425–437. doi: 10.1002/(SICI)1522-2586(200004)11:4<425::AID-JMRI11>3.0.CO;2-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Jacobs, M. A., Knight, R. A., Windham, J. P., Zhang, Z. G., Soltanian-Zadeh, H., Goussev, A. V., et al. (1999b). Identification of cerebral ischemic lesions in rat using Eigenimage filtered magnetic resonance imaging. Brain Res. 837, 83–94. doi: 10.1016/S0006-8993(99)01582-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Jacobs, M. A., Mitsias, P., Soltanian-Zadeh, H., Santhakumar, S., Ghanei, A., Hammond, R., et al. (2001b). Multiparametric MRI tissue characterization in clinical stroke with correlation to clinical outcome: part 2. Stroke 32, 950–957. doi: 10.1161/01.STR.32.4.950

PubMed Abstract | CrossRef Full Text | Google Scholar

Jacobs, M. A., Windham, J. P., Soltanian-Zadeh, H., Peck, D. J., and Knight, R. A. (1999a). Registration and warping of magnetic resonance images to histological sections. Med. Phys. 26, 1568–1578. doi: 10.1118/1.598671

PubMed Abstract | CrossRef Full Text | Google Scholar

Jacobs, M. A., Zhang, Z. G., Knight, R. A., Soltanian-Zadeh, H., Goussev, A. V., Peck, D. J., et al. (2001a). A model for multiparametric MRI tissue characterization in experimental cerebral ischemia with histological validation in rat: part 1. Stroke 32, 943–949. doi: 10.1161/01.STR.32.4.943

PubMed Abstract | CrossRef Full Text | Google Scholar

Kilkenny, C., Browne, W. J., Cuthill, I. C., Emerson, M., and Altman, D. G. (2010). Improving bioscience research reporting: the ARRIVE guidelines for reporting animal research. PLoS Biol. 8:e1000412. doi: 10.1371/journal.pbio.1000412

PubMed Abstract | CrossRef Full Text | Google Scholar

Klein, S., Staring, M., Murphy, K., Viergever, M. A., and Pluim, J. (2010). elastix: a toolbox for intensity-based medical image registration. IEEE Trans. Med. Imaging 29, 196–205. doi: 10.1109/TMI.2009.2035616

PubMed Abstract | CrossRef Full Text | Google Scholar

Leithner, C., Füchtemeier, M., Jorks, D., Mueller, S., Dirnagl, U., and Royl, G. (2015). Infarct volume prediction by early magnetic resonance imaging in a murine stroke model depends on ischemia duration and time of imaging. Stroke 46, 3249–3259. doi: 10.1161/strokeaha.114.007832

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, F., Liu, K. F., Silva, M. D., Meng, X., Gerriets, T., Helmer, K. G., et al. (2002). Acute postischemic renormalization of the apparent diffusion coefficient of water is not associated with reversal of astrocytic swelling and neuronal shrinkage in rats. ANJR Am. J. Neuroradiol. 23, 180–188.

PubMed Abstract | Google Scholar

Llovera, G., and Liesz, A. (2016). The next step in translational research: lessons learned from the first preclinical randomized controlled trial. J. Neurochem. doi: 10.1111/jnc.13516

PubMed Abstract | CrossRef Full Text | Google Scholar

Longa, E. Z., Weinstein, P. R., Carlson, S., and Cummins, R. (1989). Reversible middle cerebral artery occlusion without craniectomy in rats. Stroke 20, 84–91. doi: 10.1161/01.STR.20.1.84

PubMed Abstract | CrossRef Full Text | Google Scholar

Milidonis, X., Marshall, I., Macleod, M. R., and Sena, E. S. (2015). Magnetic resonance imaging in experimental stroke and comparison with histology: systematic review and meta-analysis. Stroke 46, 843–851. doi: 10.1161/strokeaha.114.007560

PubMed Abstract | CrossRef Full Text | Google Scholar

Moraga, A., Gómez-Vallejo, V., Cuartero, M. I., Szczupak, B., San Sebastián, E., Markuerkiaga, I., et al. (2016). Imaging the role of toll-like receptor 4 on cell proliferation and inflammation after cerebral ischemia by positron emission tomography. J. Cereb. Blood Flow Metab. 36, 702–708. doi: 10.1177/0271678X15627657

PubMed Abstract | CrossRef Full Text | Google Scholar

Niimi, T., Imai, K., Maeda, H., and Ikeda, M. (2007). Information loss in visual assessments of medical images. Eur. J. Radiol. 61, 362–366. doi: 10.1016/j.ejrad.2006.09.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Rekik, I., Allassonnière, S., Carpenter, T. K., and Wardlaw, J. M. (2012). Medical image analysis methods in MR/CT-imaged acute-subacute ischemic stroke lesion: Segmentation, prediction and insights into dynamic evolution simulation models. A critical appraisal. NeuroImage Clin. 1, 164–178. doi: 10.1016/j.nicl.2012.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Rousson, M., and Deriche, R. (2002). “A variational framework for active and adaptative segmentation of vector valued images,” in MOTION'02: Proceedings Workshop on Motion and Video Computing (Washington, DC: IEEE Computer Society), 56.

Soltanian-Zadeh, H., Pasnoor, M., Hammoud, R., Jacobs, M. A., Patel, S. C., Mitsias, P. D., et al. (2003). MRI tissue characterization of experimental cerebral ischemia in rat. J. Magn. Reson. Imaging 17, 398–409. doi: 10.1002/jmri.10256

PubMed Abstract | CrossRef Full Text | Google Scholar

Sotak, C. H. (2002). The role of diffusion tensor imaging in the evaluation of ischemic brain injury – a review. NMR Biomed. 15, 561–569. doi: 10.1002/nbm.786

PubMed Abstract | CrossRef Full Text | Google Scholar

Sunkin, S. M., Ng, L., Lau, C., Dolbeare, T., Gilbert, T. L., Thompson, C. L., et al. (2013). Allen brain atlas: an integrated spatio-temporal portal for exploring the central nervous system. Nucl. Acids Res. 41, D996–D1008. doi: 10.1093/nar/gks1042

PubMed Abstract | CrossRef Full Text | Google Scholar

van Dorsten, F. A., Olàh, L., Schwindt, W., Grüne, M., Uhlenküken, U., Pillekamp, F., et al. (2002). Dynamic changes of ADC, perfusion, and NMR relaxation parameters in transient focal ischemia of rat brain. Magn Reson. Med. 47, 97–104. doi: 10.1002/mrm.10021

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Cheung, P. T., Shen, G. X., Bhatia, I., Wu, E. X., Qiu, D., et al. (2007). Comparing diffusion-weighted and T2-weighted MR imaging for the quantification of infarct size in a neonatal rat hypoxic-ischemic model at 24 h post-injury. Int. J. Dev. Neurosci. 25, 1–5. doi: 10.1016/j.ijdevneu.2006.12.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Weber, R., Ramos-Cabrer, P., and Hoehn, M. (2006). Present status of magnetic resonance imaging and spectroscopy in animal stroke models. J. Cereb. Blood Flow Metab. 26, 591–604. doi: 10.1038/sj.jcbfm.9600241

PubMed Abstract | CrossRef Full Text | Google Scholar

Yao, X., Derugin, N., Manley, G. T., and Verkman, A. S. (2015). Reduced brain edema and infarct volume in aquaporin-4 deficient mice after transient focal cerebral ischemia. Neurosci. Let. 584, 368–372. doi: 10.1016/j.neulet.2014.10.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Yin, J., Han, P., Tang, Z., Liu, Q., and Shi, J. (2015). Sirtuin 3 mediates neuroprotection of ketones against ischemic stroke. J. Cereb. Blood Flow Metab. 35, 1783–1789. doi: 10.1038/jcbfm.2015.123

PubMed Abstract | CrossRef Full Text | Google Scholar

Yip, H. K., Yuen, C. M., Chen, K. H., Chai, H. T., Chung, S. Y., Tong, M. S., et al. (2016). Tissue plasminogen activator deficiency preserves neurological function and protects against murine acute ischemic stroke. Int. J. Cardiol. 205, 133–141. doi: 10.1016/j.ijcard.2015.11.168

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, S., Bai, Y. Y., Liu, Y., Gao, X., Li, Y., Changyi, Y., et al. (2015). Salvaging brain ischemia by increasing neuroprotectant uptake via nanoagonist mediated blood brain barrier permeability enhancement. Biomaterials 66, 9–20. doi: 10.1016/j.biomaterials.2015.07.006

PubMed Abstract | CrossRef Full Text | Google Scholar