Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Combining unsupervised and supervised learning for predicting the final stroke lesion

Combining unsupervised and supervised learning for predicting the final stroke lesion Predicting the final ischaemic stroke lesion provides crucial information regarding the volume of salvageable hypoperfused tissue, which helps physicians in the dicult decision-making process of treatment planning and intervention. Treatment selection is in- fluenced by clinical diagnosis, which requires delineating the stroke lesion, as well as characterising cerebral blood flow dynamics using neuroimaging acquisitions. Nonetheless, predicting the final stroke lesion is an intricate task, due to the variability in lesion size, shape, location and the underlying cerebral haemodynamic processes that occur after the ischaemic stroke takes place. More- over, since elapsed time between stroke and treatment is related to the loss of brain tissue, assessing and predicting the final stroke lesion needs to be performed in a short period of time, which makes the task even more complex. Therefore, there is a need for automatic methods that predict the final stroke lesion and support physicians in the treatment decision process. We propose a fully automatic deep learning method based on unsupervised and supervised learning to predict the final stroke lesion after 90 days. Our aim is to predict the final stroke lesion location and extent, taking into account the underlying cerebral blood flow dynamics that can influence the prediction. To achieve this, we propose a two-branch Restricted Boltzmann Machine, which provides specialized data-driven features from di erent sets of standard parametric Magnetic Resonance Imaging maps. These data-driven feature maps are then combined with the parametric Magnetic Resonance Imaging maps, and fed to a Convolutional and Recurrent Neural Net- work architecture. We evaluated our proposal on the publicly available ISLES 2017 testing dataset, reaching a Dice score of 0.38, Hausdor Distance of 29.21 mm, and Average Symmetric Surface Distance of 5.52 mm. Keywords: Deep Learning, Image Prediction, Magnetic Resonance Imaging, Stroke 1. Introduction caused by thrombolysis, haemodynamic factors, or em- bolic causes (Grysiewicz et al., 2008). Due to vessel oc- Stroke is the second leading cause of death worldwide clusion, the insucient supply of oxygenated blood to (World Health Organization et al., 2014), being classified brain cells leads to hypoperfused brain tissue, trigger- in two types: ischaemic and haemorrhagic (Grysiewicz ing cellular mechanisms to preserve the integrity of the et al., 2008). Ischaemic stroke is the most common type, cell. The hypoperfused area consists of tissue at risk that resulting from an occlusion of a vessel, which can be can be salvaged, being designated penumbra. As time passes, in the absence of flow restoration or sucient col- lateral blood flow supply, the hypoperfused tissue even- Corresponding author: Department of Industrial Electronics, Cam- tually reaches a non-salvageable state designated core or pus Azurem, ´ Guimaraes, ˜ Portugal. infarct tissue (Memezawa et al., 1992). Email addresses: id6376@alunos.uminho.pt (Adriano Pinto ), csilva@dei.uminho.pt (Carlos A. Silva) Diagnosis and treatment of ischaemic stroke relies on Preprint submitted to Elsevier January 5, 2021 arXiv:2101.00489v1 [eess.IV] 2 Jan 2021 neuroimaging acquisitions, where Computed Tomogra- two-pathway architecture, trained with two subsets of phy (CT) and Magnetic Resonance Imaging (MRI) are the standard parametric MRI maps. One subset encompasses preferred imaging modalities (Gonzalez et al., 2007). CT the Time-To-Peak (TTP), Mean Transit Time (MTT), imaging remains the most used acquisition due to its ra- Time-to-Maximum (Tmax), and Apparent Di usion Co- pidity and availability (Gonzalez et al., 2007). However, ecient (ADC). The second set contains the ADC, the multi-parametric MRI provides a higher sensitivity in de- relative Cerebral Blood Volume (rCBV), and the relative tecting early ischaemic stroke and assessing the penum- Cerebral Blood Flow (rCBF). In a second stage, the fea- bra region (Gonzalez et al., 2007). Treatment consists in ture maps computed by the RBMs are combined with the restoring tissue perfusion levels, also known as reperfu- standard parametric MRI maps to form the input of a su- sion, by performing mechanical thrombectomy or throm- pervised deep learning architecture composed by Convo- bolysis. Since ischaemic stroke is a dynamic process that lutional Neural Networks (CNNs) and Recurrent Neural evolves over time, the treatment is only possible up to 24 Networks (RNNs). The proposed architecture was evalu- hours, where viable neurones still persist (El Tawil and ated using the publicly-available ISLES 2017 dataset. Muir, 2017; Zivelonghi and Tamburin, 2018). So, expert 1.1. On the complexity of predicting the final infarct physicians must evaluate the benefits and risks of mechan- stroke lesion ical thrombectomy before an intervention, since it may cause haemorrhage, vascular injury, and other complica- In acute ischaemic stroke, the clinical evaluation of the tions (Powers et al., 2018). If performed, the success of standard parametric maps (e.g. ADC and Tmax) can iden- the intervention is assessed radiologically via angiogra- tify infarct tissue and tissue that will infarct in the ab- phy imaging and scored by a qualitative expert-generated sence of therapeutic intervention. In this analysis, the in- scale designated the standardized Thrombolysis in Cere- farct tissue, is identified by the hypointense regions of the bral Infarction (TICI) scale (Higashida et al., 2003). Dur- ADC map, which characterize tissue with limited di u- ing the decision-making process, the physician needs to sion (Butcher and Emery, 2010a). Hypoperfused tissue, assess the nature and location of the lesion alongside i.e. tissue that will infarct, is identified by hyperintense pathophysiological factors such as age, presence of co- regions of the Tmax map, indicating an increased arrival morbidities, and collateral circulation (Liebeskind, 2003). time of contrast agent (Butcher and Emery, 2010b). How- The latter is of utmost importance in ischaemic stroke. ever, to correctly predict the final ischaemic stroke lesion, The presence of collateral circulation, where a secondary besides considering the complex time-evolving transfor- network of vessels is responsible for granting cerebral mation of hypoperfused tissue to infarcted tissue, it is blood flow to the lesioned tissue, increases the chances also necessary to appraise the impact of the clinical inter- of a successful reperfusion (Liebeskind, 2003). Asserting vention, thrombectomy, on the underlying brain perfusion the potential ecacy of treatment can be time-consuming and di usion. and prone to inter- and intra-variability among physicians, A successful thrombectomy should restore the perfu- which is further potentiated when performed in a clinical sion levels, recovering the hypoperfused tissue. How- emergency environment (Coutts et al., 2003). Moreover, ever, several factors may a ect the reperfusion, limiting since time is critical, MRI acquisitions are optimized for the degree of success of the intervention. To better un- speed, which is accomplished by reducing the resolution derstand the nuances of the clinical intervention, consider (Gonzalez ´ et al., 2011), making the prediction of the final the cases presented in Figure 1. In the first case, Figure stroke lesion an intricate task. Thus, automatic prediction 1a, the ADC does not present any hypointense region, so of a stroke lesion at a given time since stroke has a great no infarct tissue may be identified, and we should expect potential to guide physicians in this time-critical decision- a complete recovery of the hypoperfused tissue indicated making process. by Tmax; however, the follow-up delineation obtained af- We propose a novel automatic method based on unsu- ter thrombectomy presents a large final lesion, which is pervised and supervised deep learning. We utilize Re- explained by an unsuccessful intervention. In the second stricted Boltzmann Machines (RBMs) to jointly charac- case, Figure 1b, we observe a final infarct lesion that is terise the lesion and blood flow information through a smaller than the hypointense region present in the ADC 2 (a) (b) Figure 1: ADC and Tmax parametric maps of two patient cases from ISLES 2017 training set, and the final lesion delineated at a 90-day follow- up, overlapped with the onset ADC: patient 0036 (Figure 1a) with an unsuccessful reperfusion, and patient 0006 (Figure 1b), where the clinical intervention was successful. (Figure 1b arrow). This indicates reversible di usion re- Challenge in 2016 and 2017, new methods have been pro- striction, which is a rare case (Labeyrie et al., 2012) and posed. These aim to predict at a 90-day time-window. was only possible to identify by a follow-up T2-weighted Rose et al. (2001) proposed a two-stage method based acquisition. So, an automatic method for predicting the fi- on parametric perfusion and di usion MRI maps. On the nal stroke lesion has not only to capture the time-evolving first stage, the method defines a region of interest (ROI) process of di usion and perfusion, but also to consider di- based on the intensity signal of the standard parametric rectly or indirectly the degree of success of the thrombec- maps, the MTT, Cerebral Blood Flow (CBF), Cerebral tomy, which may condition the final lesion either to be Blood Volume (CBV), and Di usion-Weighted Imaging confined to the hypointense region of the ADC map, or (DWI). The second stage performs stroke tissue predic- to grow to brain tissue areas that are hyperintense in the tion using Gaussian mixture models trained in di erent Tmax. Due to the time-evolving process of di usion and sets of parametric maps. Bauer et al. (2014) used Random perfusion, the complexity of predicting the lesion will ag- Forests to segment or predict the final stroke lesion de- gravate as we move from a target window of some days to pending on whether acute stroke imaging or three-month several months. follow-up imaging was available, respectively. McKin- The complexity of the evaluation process may be also ley et al. (2017) also used a two-stage classification ap- observed in the inter-rater agreement of expert radiolo- proach as in Rose et al. (2001) for lesion characterisa- gists in ISLES 2017 dataset, which obtained a Dice score tion and lesion prediction, where each stage consists of of 0:58  0:20 on delineating the lesion using a 90-day two sets of Random Forests (RFs) classifiers. The first follow-up T2-weighted acquisition (Winzeck et al., 2018). stage aims to define a ROI that encompasses the hypop- erfused region. In the first set, each classifier is trained with features extracted from di erent sets of MRI para- 1.2. Previous Work metric maps. Having defined the location and extension Contrary to stroke lesion segmentation, where several of the lesion, a second set of two RFs performs stroke methods have already been proposed (Rekik et al., 2012; tissue prediction. Such classifiers were trained on dif- Maier et al., 2017), the complexity of predicting the fi- ferent sets of patients, stratified by the TICI score. One nal stroke lesion has only recently attracted attention in classifier is trained in patients with unsuccessful reperfu- the medical imaging community. For predicting the final sion interventions, whereas a second classifier is trained stroke lesion several methods have been already proposed in patients with successful reperfusion. The final predic- based on multivariate linear regression models (Scalzo tion is obtained by combining the results of both clas- et al., 2012; Rose et al., 2001; Kemmling et al., 2015), sifiers, using a logistic regression model. Scalzo et al. decision trees (McKinley et al., 2017; Bauer et al., 2014), (2012) proposed a framework for stroke tissue prediction, and CNNs (Choi et al., 2016). Furthermore, with the re- which characterises the state of the lesion four days af- lease of Ischaemic Stroke LEsion Segmentation (ISLES) ter clinical intervention (thrombectomy). From the Fluid 3 Attenuation Inversion Recovery (FLAIR) MRI sequence, olution, while in the second branch the input resolution ADC and Tmax MRI maps, the method applies a regres- was lowered by a factor of 3. The output of each branch sion model that learns the behaviour of neighbouring vox- is transformed to the same scale and merged by two fully els within a cuboid. Kemmling et al. (2015) proposed connected layers. The network is trained with four dif- a multi-modality approach based on CT and MRI maps ferent sets of hyper-parameters. These four networks are with non-imaging clinical meta-data, namely the TICI used as an ensemble, whose prediction is obtained by av- score and the time to treatment of each patient, to perform eraging the output of each one. Similarly, Niu et al. (2018) stroke tissue prediction. used multiple scales of overlapping 3D patches to capture In another line of research, authors have investigated local and global spatial information. In the review pa- the use of deep learning (Choi et al., 2016; Nielsen et al., per of Winzeck et al. (2018), Rivera et al. also built on 2018; Robben et al., 2020) for stroke tissue prediction. the work of Kamnitsas et al. (2017) and Milletari et al. Choi et al. (2016), the winner approach at ISLES 2016 (2016), by proposing a scheme to extract di erent patch Challenge, proposed an ensemble of twelve CNN archi- resolutions, independent of each other, that are fed into tectures, grouped into two sets of networks. The first four di erent paths. Afterwards, a fully connected layer group comprehends four 3D U-Nets (Ronneberger et al., combines all the outputs to perform stroke tissue predic- 2015) performing voxel-wise tissue prediction. The sec- tion. Pisov et al. (2017) employed an ensemble strategy ond group of networks uses two-pathway Fully Con- by combining di erent CNN-based architectures to over- nected Networks (FCNs) performing two types of patch- come the strong anisotropy of the data. As summarized by wise classification. One set of FCNs classifies a patch Winzeck et al. (2018), the work of Yoon et al. proposed a as lesion if it includes any lesion voxel. The other set two-stage gated CNN. In a first stage, the authors perform of FCNs classifies a patch as lesion if the central voxel lesion detection and delineation. Afterwards, based on the is a lesion. After merging the two pathway FCN, the probability maps of the first stage, a second CNN archi- method incorporates meta-data by adding a dense layer tecture processes the regions where the probability maps of clinical predictors merged with the imaging output of of healthy tissue and lesion are close to each other. Pinto each network. The final stroke lesion prediction results et al. (2018b) made use of temporal perfusion imaging, from a weighted merging of all models. Mok and Chung the Dynamic Susceptibility Contrast-MRI, in a U-Net ar- (2017) applied deep adversarial training for stroke tis- chitecture. This architecture aims to temporally process sue prediction in an ensemble of U-Nets. Monteiro and and extract deep features, which are then combined with Oliveira (2017) proposed a method based on the V-Net a second feature step of another U-Net network, which architecture (Milletari et al., 2016). The training was con- was trained on the standard parametric maps. Using a ducted with a custom loss function that applies a weighted large CT dataset, Robben et al. (2020) predicted the fi- sum between Dice score and cross entropy. Lucas and nal infarct stroke lesion with a temporal window ranging Heinrich (2017) proposed the use of a U-Net architec- from 24 hours to 5 days. The authors considered spatio- ture, which combines patches from the MRI maps in the temporal CT perfusion as input to a deep neural network same slice, with patches from 3 neighbouring slices and inspired in the architecture proposed by Kamnitsas et al. 2 hemispheric flips. In the expanding path of the U-Net, (2017). Additionally, the model combines CT neuroimag- each level computes a Dice loss for the healthy tissue and ing with clinical meta-data. Nielsen et al. (2018) pro- for the final lesion, after the softmax activation. After- posed a method based on the SegNet architecture (Badri- wards, all losses are summed up, having the loss of the narayanan et al., 2015), predicting on a 30-day follow-up lesion and healthy tissue weighted according to a prior acquisition based on a private dataset. probability (Winzeck et al., 2018). Robben and Suetens Principal and collateral blood flow has been consid- (2017) employed a CNN-based architecture inspired by ered either directly by modelling the temporal perfusion Kamnitsas et al. (2017). The authors proposed to com- imaging (Pinto et al., 2018b), or indirectly by perfusion bine the MRI inputs with clinical meta-data, before feed- and di usion parametric maps (Choi et al., 2016; Maier ing them to each branch of a two-pathway 3D network. et al., 2017; Scalzo et al., 2012), or through clinical in- In the first branch the input is kept with the original res- formation that characterises the success of the revascu- 4 larization (McKinley et al., 2017). We hypothesize that 2. Methods modelling the haemodynamics of the brain when artery In this work, predicting the final infarct stroke lesion occlusion occurs can be beneficial for predicting the final consists of delineating the lesion’s spatial extent at a 90- stroke lesion. So, in this work, we investigate the rep- day follow-up time-point, using multi-parametric MRI resentation of the haemodynamics through an unsuper- imaging, namely the ADC, MTT, TTP, Tmax, rCBF, and vised learning model. Contrary to previous approaches, rCBV, which are acquired at the onset time-point. The we propose grouping the input maps according to their architecture of the proposed system and its main compo- subjacent physical meaning and encoding each group sep- nents are described in the following subsections. arately with an RBM. As groups, we investigated the time-resolved perfusion maps (Tmax, TTP, MTT), and the 2.1. Architecture blood-flow-dynamic related maps (rCBF, rCBV) (Butcher and Emery, 2010a,b). Our proposal of combining features The overall architecture of the proposed method can be obtained unsupervisedly and supervisedly was motivated divided into two functional blocks, as shown in Figure 2. by the knowledge that unsupervised models learn struc- The first functional block performs unsupervised rep- tural features of the original image, while the supervised resentation learning using two unsupervised models, models learn features conditioned on the label, so there namely RBMs. This unsupervised block provides new is potential for obtaining richer and more discriminative features that represent structural information that comple- features by joining both types of models. ments the standard parametric MRI maps, enhancing the capacity of our model to predict the final infarct lesion. In our approach, we aim to model the clinical procedure, 1.3. Contributions which first locates and delineates the lesion at current time, and then considers the blood flow haemodynamic This work presents an automatic approach for predict- that might influence the final stroke lesion prediction. ing the final stroke lesion, using onset neuroimaging data. This procedure is encoded in our two-path RBM. The The main contributions are: first RBM is responsible for capturing information on le- sion location and extension, referred to as the RBM . - The use of unsupervised models for extracting Lesion The second RBM, RBM , aims to capture blood flow structural features of time-resolved perfusion and Haemo haemodynamics information (e.g. collateral circulation), blood-flow-dynamic related MRI maps for predict- which has been identified as a key factor by physicians ing stroke lesion. when assessing stroke final infarct lesion in clinical re- - The use of local and long spatial context provided ports (Berkhemer et al., 2016; Menon et al., 2015). On by gated recurrent neural networks for relating struc- one hand, to locate the onset ischaemic stroke lesion, the tural features and image information when learn- RBM considers standard parametric maps that char- Lesion ing features conditioned on the label in a supervised acterise the arrival times and mean transit times of the model. contrast agent. In the presence of an ischaemic lesion, the occluded vessel can decrease or interrupt the normal brain - The proposal of a competitive system which outper- perfusion, translating into hyperintense regions on time- forms state-of-the-art methods to predict the final in- related parametric maps (Butcher and Emery, 2010b). On farct stroke lesion, in ISLES 2017 Challenge dataset. the other hand, the RBM considers standard para- Haemo metric maps that characterise the amount of blood being The remainder of the paper is organized as follows. delivered in unit of time, which correlates to the cerebral Section 2 describes the fundamental components of the blood flow haemodynamics (Butcher and Emery, 2010b). proposed method. Section 3 describes the dataset, the Thus, the RBM considers the MTT, TTP and Tmax Lesion evaluation procedure and the setup. The results and the perfusion maps, while the RBM the rCBV and rCBF Haemo discussion are addressed in Section 4. Finally, in Section perfusion maps. Regarding the ADC standard di usion 5 we present the main conclusions. map, it is present in both RBM and RBM , since Lesion Haemo 5 Figure 2: Overview of the proposed method for predicting the final stroke lesion. In the supervised learning block, the input data dimensions are defined for each operation. it provides higher brain structural information and allows data into a feature vector is performed through the interac- the identification of tissue that is already infarcted. This tion of states between the visible and hidden units, which separation of the input imaging allows the RBM to learn is learned by minimizing an energy function. specific feature sets, which may enable the method to The complete pipeline of the unsupervised block is analyse dicult cases where information concerning the shown in Figure 3 and detailed in Section 3.4. The blood flow can have a favourable impact on the lesion pre- RBM and RBM function as feature generators Lesion Haemo diction. that output two complementary sets of feature maps N The second functional block consists of a deep learning and N . These features characterise the structure of the architecture that comprehends 2D convolutional blocks in images; however, we are interested only on the most dis- a U-Net-based structure, alongside recurrent blocks. As tinctive details. So, after training the RBMs, we perform imaging input data, we combine the standard parametric feature selection to reduce the generated feature space, maps with feature maps from each RBM, totalling 18 in- obtaining smaller but representative feature sets M and put feature maps. M , such thatjM j  jN j, for i 2 [1; 2], where the opera- 2 i i tor j:j denotes the cardinality of a set. In the feature selec- 2.2. Restricted Boltzmann Machines tion, we would like to select the features from the RBM that encodes the MRI maps, but also that correlates with The RBM is an undirected graphical model consti- the stroke prediction. Since the RBM is an unsupervised tuted by two layers of nodes: a visible layer and a hid- method, we compute the Normalized Mutual Information den layer (Rumelhart and McClelland, 1986). Each node to quantify the statistical dependence between each gener- has a weighted connection to all nodes in the other layer ated feature and the respective input MRI map, as defined (Rumelhart and McClelland, 1986). However, there are by Equation 1 (Vinh et al., 2010): no connections among nodes of the same layer. Orig- inally, Rumelhart and McClelland (1986) proposed the RBM to learn from binary data on both layers. How- M I(MRI ; Feat ) x y N M I (MRI ; Feat ) = 2 ; (1) sum x y ever, this does not represent well continuous real-valued H(MRI ) + H(Feat ) x y input data, which is the case of MRI data. Therefore, we model the visible nodes as linear units with independent where M I(:) is the mutual information between an MRI Gaussian noise. The hidden nodes are modelled as Noisy parametric map, MRI , and an output feature, Feat ; H(:) x y Rectifier Linear Units (NReLU), since they have been re- defines the entropy of a map, namely, MRI and Feat . x y ported to be suitable for feature extraction (Hinton, 2012). To relate the features of the RBM with the class label, This kind of RBM was previously used in segmentation we could use a classifier supervisedly trained. Since the tasks, such as in Pereira et al. (2019). Mapping the input neural network is trained iteratively, we use a RF clas- 6 Figure 3: Overview of the proposed unsupervised learning block. For each RBM of the unsupervised learning block, the selected features were M = M = 6. 1 2 sifer trained with the Mean Decrease Impurity (MDI) as decoder deep CNNs provide high levels of abstraction a surrogate to make the feature selection tractable. After, from the input data, increasing the global notion of con- we compute the MI and and MDI , normalize the text as the network grows deeper. However, it comes RBM RF MI by the maximum value, add both ranks and sort at a cost of a high receptive field (Zeiler and Fergus, RBM decreasingly. The best set will be the first M features. 2014). Thus, we used a 2D architecture in the plane Our selection method was inspired on the work of Pereira with the highest resolution, since the acquisition resolu- et al. (2018); however, their method cannot be directly tion is anisotropic in the dataset. Also, in the end of applied, since it would generate too many features for our the decoding path we expanded our learning block with problem. Gated RNNs. Due to their nature, Gated RNNs can cap- ture short- and long-term spatial relations, by retaining in- formation from previous nodes encoded in the time-steps. 2.3. Convolutional and Recurrent Neural Networks Hence, Gated-RNNs consider information from all previ- Our supervised functional block is based on the U-Net ous nodes when analysing the current one. This property, architecture as proposed by Ronneberger et al. (2015). when applied to imaging data, allows considering intra- The input of the U-Net considers the concatenation of slice contextual dependencies. In our work, we used a par- standard parametric maps with the sets of feature maps ticular Gated-RNN, namely the Long-Short Term Mem- extracted from the unsupervised block. In the first level ory (LSTM) (Hochreiter and Schmidhuber, 1997). How- of our encoder architecture we use four 2D convolutional ever, the LSTM was intrinsically developed to process blocks with kernel size of 3  3 and 32 channels. Af- 1D data (Hochreiter and Schmidhuber, 1997) (e.g. time- terwards, the output of the final convolutional block is series). To be applicable to 2D data, we developed an down-sampled by a factor of 2, starting the second encod- online 2D Partition layer that transforms a grid-structure ing level formed by two convolutional blocks with equal input (e.g. an image) into a one-dimensional sequence. In- kernel size but doubling the number of feature maps. The spired by Visin et al. (2016), the 2D Partition layer was third level of encoding follows the same pattern. The de- predefined with a neighbourhood of 2  2, where each coder level mimics the encoder counterpart. As in Ron- time-step is characterised by a feature space of four vox- neberger et al. (2015) we only used long skip connec- els. After, two Bidirectional LSTM layers are employed tions among encoder and decoder levels. These encoder- 7 along the left-right and frontal-dorsal directions followed 3.2. Evaluation Metrics by an up-sampling layer. These four layers, referred as We evaluated our proposal with five metrics, which are the Gated Recurrent block, are shown in Figure 2. In our the same ones computed by the online ISLES 2017 bench- supervised functional block, two Gated Recurrent blocks mark platform: Dice Similarity Score, Hausdor Distance were used, where the Bidirectional LSTMs have 64 and (HD), Average Symmetric Surface Distance (ASSD), Pre- 32 hidden layers, respectively. The impact of the main cision, and Recall (Kistler et al., 2013). components is evaluated in an ablation study in the exper- Dice score measures the spatial overlap between two iments. volumes. HD corresponds to the highest distance between surface points of di erent volumes, which characterise spatial outliers in the prediction. ASSD quantifies the 3. Experimental Setup average distances between the volumes’ surface. Preci- sion quantifies the proportions of correctly classified cases We evaluated the proposed approach on the publicly within a class, while Recall corresponds to the proportion available ISLES 2017 dataset and on a private dataset. of positive cases correctly identified as such. ISLES has an online benchmark platform (Kistler et al., 2013) that performs automatic evaluation (SMIR On- 3.3. Image pre- and post-processing line Platform, 2017). In this section we describe the Since MRI acquisitions were acquired from di erent dataset, the training and evaluation, and the main hyper- centers and configurations (Winzeck et al., 2018), for parameters of our method. each patient we resized all maps to a common volume of dimension of 256  256  32. Afterwards, the ADC 3.1. Data 6 2 maps were clipped between [0; 2600]  10 mm =s and ISLES 2017 dataset encompasses 75 ischaemic stroke the Tmax maps were clipped to [0; 20s], since values be- patients, which are separated into two sets: training (n = yond these ranges are known to be biologically meaning- 43) and testing (n = 32). Both sets have patients that un- less (McKinley et al., 2017). Finally, a linear scaling was derwent mechanical thrombectomy. Each patient is char- applied across all maps, to the range [0; 255]. The im- acterised by six 3D parametric MRI maps: di usion ADC ages are resized to its original size, after we perform the map, perfusion rCBF, rCBV, TTP, MTT and Tmax maps. prediction. In addition to the standard parametric maps, each case is We applied a morphological filtering as post- also characterised by a manual delineation of the lesion. processing, but since the final stroke lesion presents This refers to the 90-day stroke lesion delineated with ac- a wide range of lesion volumes (Winzeck et al., 2018), cess to the follow-up T2-weighted acquisition. However, we removed only small connected components with less the manual delineation is only available for the training than 250 voxels. This step was kept fixed for all the set, while the follow-up T2-weighted imaging is not dis- evaluated models. closed for any set. All parametric MRI maps are already 3.3.1. Data Augmentation co-registered and skull-stripped (Winzeck et al., 2018). Figure 4 (top row) shows an example of MRI maps, Data augmentation can be used to increase the number alongside the manual lesion delineation, the Ground Truth of training samples and reduce over-fitting (Krizhevsky (GT), of a patient. et al., 2012). Due to the relatively small size of the train- The private dataset considers 23 acute ischaemic stroke ing dataset, we employed artificial data augmentation in patients that underwent clinical therapy, acquired at Bern the supervised portion of our proposal. For each sample, University Hospital in Switzerland. As in ISLES 2017, we applied rotations of 90 , 180 , 270 . each patient is characterized by the same six parametric 3.4. Settings and model training maps, being the final lesion manually delineated at 90-day follow-up T2. The parametric maps were co-registered Unsupervised functional block. The unsupervised func- followed by skull-stripping with FSL BET2 on the co- tional block was trained by optimizing the negative log- registered follow-up T2 image (Jenkinson et al., 2012). likelihood of the data. However, since computing the 8 gradient is generally intractable, we performed the train- jVj p g ing by approximating the gradient with Contrastive Diver- i i Soft Dice loss = : (2) P P jVj jVj gence with one step of alternating Gibbs sampling (Hin- 2 2 p + g i i i i ton, 2012). The training process of an RBM can be dif- In the soft dice loss, the sum occurs over the set V ficult if one tries to learn the parameter  of the energy of voxels belonging to the predicted output patch, where function, which corresponds to the standard deviation of p 2 P denotes the probability of a voxel i in the output the Gaussian noise of a visible node i (Hinton, 2012). Ac- patch and g 2 G corresponds to the respective ground- cording to Hinton (2012), we normalize each component truth label voxel. of the data with zero mean and unit variance, and define The method was implemented using Keras with Ten- = 1. In Table 1, we present the settings used for train- sorflow backend, in a workstation equipped with a GTX ing the unsupervised model. 1080 Ti 11 GB. Prediction time takes around 20s per pa- For training each RBM, we randomly extract 3D tient. patches of shape 7  7  3 from the respective input set of MRI maps, C. Then, the 3D patches are reshaped into a 1D vector and fed into the visible layer of the RBM, hav- 4. Results and Discussion ing an input of size m = 7 7 3jCj, as shown in Figure In this section, we discuss the impact of the main 3. After training, we extract features from the NReLU contributions, namely, the incorporation of unsupervised units noise-free activations. These units exhibit intensity learning with supervised learning and the Gated Recur- equivariance when the bias has zero value, and they are rent blocks. Then, we compare our method with the state noise free units (Nair and Hinton, 2010). Due to the large of the art in ISLES 2017 Challenge. Finally, we delve on number of extracted feature maps (jN j = jN j = 600), 1 2 we perform a feature selection step, as described in Sec- the diculty of predicting the final infarct stroke lesion. tion 2.2, where M = M = 6. The most appropriate 1 2 cardinality of M is discussed in Section 4.1.1. 4.1. Ablation Study The ablation study aims to gradually measure the im- Table 1: Model training parameters for the unsupervised and supervised portance of the main components and consequently assert functional blocks. on the contribution of each component to the overall per- formance. Thus, we start by evaluating the importance of the unsupervised feature generator and the proposed in- Functional Block Parameter Description put grouping. After, we focus on the use of the Gated Optimizer SGD with momentum (lr = 1 10 ) Unsupervised Patch shape 7 7 3 Recurrent block and the choice of the dimensionality of Batch size 32 the spatial context. Optimizer ADAM (lr = 1 10 ) Supervised Patch shape 84 84 4.1.1. Unsupervised feature generation Batch size 4 We hypothesize that grouping the parametric MRI maps according to their physical meaning and encoding each group with an RBM has potential to extract better Supervised functional block. As for the supervised func- features to characterise the stroke lesion and the blood tional block, the complete settings of the training are haemodynamics. We perform several experiments to cor- given in Table 1. For each subject, 350 patches were ran- roborate this working hypothesis. In all experiments, the domly sampled. The training comprehended 36 subjects, while the remaining 7 subjects were used for validation. The settings were optimized through cross-validation in a Additional details of setting and model training are provided in the previous work (Pinto et al., 2018a). For training, we used supplementary material. Also, the source code for reproducing the seg- soft Dice loss function (Milletari et al., 2016). It is defined mentations, the models’ weights and segmentations can be found at: as: https://github.com/apinto92/stroke_prediction.git. 9 parametric MRI maps are also used as input to the super- Grouping parametric MRI maps according to the subja- vised block. The results are presented in Table 2. Figure cent physical meaning. In this experiment, we grouped 4 presents feature maps encoded by the RBMs and the the parametric maps according to their underlining physi- respective MRI maps. cal meaning together with ADC map in each group. Each group was encoded with an RBM. Comparing isolatedly Grouping all parametric MRI maps in a single group. the use of each group of features, we verify that RBM Lesion We considered, first, the e ect of encoding all paramet- had a higher average Dice score compared to using only ric maps using a single RBM. We varied the number of the parametric maps as input to the supervised block. selected features from the RBM, observing that in all The increase in the average Dice score was obtained by cases, the average Dice score is equal or lower than us- a higher average Recall. Also, we observe an improve- ing only the parametric maps as input to the supervised ment in all distance metrics. The experiment of using block. Also, using 12 features presented the lowest av- RBM presented the lowest average Dice and Re- Haemo erage Dice score. The use of 3 or 6 obtained the same call, as well as higher average distance metrics. However, average Dice score, having the second a lower average RBM presented higher average Precision, contrary Haemo Hausdor distance. So, based on the metrics, we may to RBM , which motivated the study on the combina- Lesion conclude that there is no clear gain in using the features tion of features from RBM with RBM besides Haemo Haemo generated by the RBM, at least, when we encode all the the parametric maps. The results of this experiment are parametric maps with a single RBM. presented in Table 2. We may observe that this combi- Since, the selection of 6 features also includes the pre- nation obtained the highest average Dice and Precision, vious top 3 features, we compared the normalized mutual as well as the lowest average distance metrics. How- information between them. As shown in Figure 4, the top ever, this improvement could have been originated from 3 features have low values of normalized mutual informa- the combination of maps according to a specific com- tion in relation to the additional 3 features, which indi- mon property, subjacent physical meaning of the para- cates that there is additional information. For this reason, metric maps, in each group, or because we reduced the we chose 6 as the number of features in the subsequent number of maps from 6 to 3 in each group. And this re- experiments. Table 2: Results obtained with di erent configurations of the unsupervised feature generator block in ISLES 2017 testing set. Each metric represents the mean  standard deviation. Underlined values correspond to the highest mean. Supervised Block Unsupervised Block Dice HD ASSD Precision Recall FCN G-RNN – U-Net LSTM 0.30  0.21 36.58  16.62 6.96  5.08 0.30  0.26 0.55  0.31 RBM (3 Feat.) U-Net LSTM 0.30  0.21 38.93  18.80 6.55  4.22 0.29  0.24 0.61  0.31 Single RBM (6 Feat.) U-Net LSTM 0.30  0.21 36.94  19.19 6.72  4.43 0.29  0.24 0.59  0.31 Single RBM (12 Feat.) U-Net LSTM 0.28  0.20 41.07  18.67 6.81  3.88 0.24  0.21 0.65  0.30 Single RBM U-Net LSTM 0.28  0.24 38.50  22.78 11.09  14.79 0.35  0.30 0.44  0.34 Haemo RBM U-Net LSTM 0.31  0.21 35.38  15.75 6.44  4.43 0.30  0.24 0.59  0.30 Lesion RBM + RBM U-Net LSTM 0.38  0.22 29.21  15.04 5.52  5.06 0.41  0.26 0.53  0.29 Lesion Haemo Two-RBMs U-Net LSTM 0.27  0.21 40.89  14.63 6.92  3.64 0.25  0.23 0.68  0.28 Mixed Three-RBMs U-Net LSTM 0.35  0.23 29.32  14.33 5.27  3.54 0.34  0.27 0.59  0.30 RBM + RBM + RBM Haemo/Less Lesion/Less ADC 10 Figure 4: Onset parametric maps of patient case 0011 in ISLES 2017 training set, alongside the final stroke lesion, at a 90-day follow-up, over the onset ADC map. The subsequent rows show the RBM features selected from the RBM , RBM and RBM , respectively. The last Lesion Haemo Single column shows the normalized mutual information, across whole dataset, among features of the same RBM. duction could have allowed a better training of the RBM. ble 2, the first experiment presented the lowest average So, we performed two complementary experiments. In Dice score and higher average distance metrics, while the the first experiment, we formed two groups with similar second experiment attained the second highest average size, but we randomly chose the parametric maps to in- Dice score, thus showing the importance of splitting the clude in each group. In the second experiment (Three- parametric perfusion maps and including the ADC map in RBMs), we changed the groups of MRI maps encoded both the RBM and RBM . Haemo Lesion in RBM and RBM by removing the ADC map Considering these experiments together, we may draw Lesion Haemo from each one. These two new groups were encoded in some conclusions. First, although CNNs are very e ec- RBM and RBM , respectively. The ADC tive in generating features from raw data, they can gen- Lesion/Less Haemo/Less was separately encoded in RBM . As presented in Ta- erate even better features if rich and complementary in- ADC 11 formation is provided. A similar conclusion was inferred the average distance metrics. by Oliveira et al. (2018) that observed improvement when Based on these experiments, we may conclude that the the coecients of the Wavelet were added as input in the CNN layers were able to extract additional information problem of retinal vessel segmentation. Here, we observe from the RBM features; however, at least to the problem a similar e ect, but using the encoding provided by an of inferring the extension of the lesion months ahead, long RBM trained unsupervisedly for the problem of stroke le- and local distance spatial relations among input voxels in- sion prediction. Second, at least to the problem of stroke troduced by Gated RNN was critical to reduce the detec- lesion prediction, when we have data with di erent latent tion of false positives, increasing substantially the average factors and we are able to group it, according to those Dice score by 6%. factors, then there is potential to extract complementary information from each group, but to mix them all together 4.1.3. Spatial context: 2D or 3D? can be detrimental. MRI images are 3D by nature, so the use of 3D fil- ters would allow capturing more context, which has the 4.1.2. Context aggregation based on gated recurrent potential to provide better prediction. Since 2D filters blocks are confined to a plane, unnatural discontinuous contour In medical imaging segmentation, which is similar to may occur in the perpendicular axis. However, as pre- our problem of inferring the extension of the lesion 90 sented previously, the resolution of MRI images in ISLES days ahead, the use of a cascade of convolutional layers dataset is not equal in all axis, being coarser along the to elaborate the features is the prevalent practise. How- axial axis. So, we studied the e ect of the spatial con- ever, as discussed previously, Gated-RNN layers are able text in our architecture. As we have two blocks, unsuper- to capture long distance spatial relations among input vox- vised and supervised blocks, the e ect on each one was els, so we performed some experiments to evaluate its evaluated separately. The results are presented in Table 4. contribution. The results are presented in Table 3. Considering the results, we observe that using 2D patches Analysing Table 3, we verify that when we just had in both blocks has lower average Dice score, than using parametric maps as input to the supervised block, adding only the parametric maps as input (baseline), because the a LSTM layer increased the average Precision, but the increase in the average Precision was not enough to com- average Recall decreased, resulting in the same average pensate the drop in the average Recall. Using 3D patches Dice score. But, a di erent behaviour is observed when for both blocks had the same performance as our base- we added the features computed from the RBMs. In this line. However, when we used 3D patches for the RBM scenario, we verify that using only CNN layers improved but 2D blocks for the U-Net block, we improved over our over having just parametric maps, which came by a higher baseline. This is the model with the highest average Dice average Precision. However, when we add the LSTM, we score without LSTM. So, we may conclude that for our have an even higher improvement, which is observed in a architecture, larger context using 3D patches was more larger increase in the average Precision, and a decrease in e ective for encoding features in the unsupervised block, Table 3: Results obtained when considering the Gated Recurrent block with and without the unsupervised learning block with ISLES 2017 testing set. Each metric represents the mean  standard deviation. Underlined values correspond to the highest mean. Supervised Block Unsupervised Block Dice HD ASSD Precision Recall FCN G-RNN U-Net – 0.30  0.21 38.83  21.10 7.08  5.15 0.26  0.23 0.64  0.30 U-Net LSTM 0.30  0.21 36.58  16.62 6.96  5.08 0.30  0.26 0.55  0.31 U-Net – 0.32  0.23 34.09  16.51 7.60  7.14 0.35  0.27 0.48  0.32 RBM + RBM [3D] Lesion Haemo U-Net LSTM 0.38  0.22 29.21  15.04 5.52  5.06 0.41  0.26 0.53  0.29 12 Table 4: Evaluation metrics obtained with di erent spatial context configurations in the unsupervised and supervised learning blocks in ISLES 2017 testing set. Each metric represents the mean  standard deviation. Underlined values correspond to the highest mean. Supervised Block Unsupervised Block Dice HD ASSD Precision Recall FCN G-RNN RBM + RBM [2D] U-Net [2D] – 0.27  0.23 36.35  14.89 9.14  12.35 0.31  0.28 0.53  0.34 Lesion Haemo U-Net [2D] – 0.32  0.23 34.09  16.51 7.60  7.14 0.35  0.27 0.48  0.32 RBM + RBM [3D] Lesion Haemo U-Net [3D] – 0.30  0.21 34.17  14.86 6.16  3.82 0.32  0.27 0.54  0.30 while 2D patches were better suited for encoding features which was manually delineated based on a follow-up T2 in the supervised U-Net-based block. MRI acquisitions, are not disclosed for public access. Considering the results, we observe that our baseline is competitive with an average Dice, being among the 4.2. Private dataset top 3 methods, and surpassing the ensemble methods of To further evaluate the generalization capacity of our Pisov et al. (2017) and Robben and Suetens (2017). Our proposal, we tested it on a private dataset and compare proposed method presented the lowest distance metrics it with the baseline method. Table 5 presents the results among all methods, especially for the Hausdor distance. obtained by the two methods. It obtained the second-best average Precision score, be- On the overall, our proposal was capable of surpassing ing surpassed by Robben and Suetens (2017) The authors the baseline model, attaining an higher average Dice, Pre- proposed the integration of meta-data information, using cision and distance metrics, which were statistically sig- a two-pathway 3D network in an ensemble; however, our nificant (Wilcoxon Signed Ranked test with p value < experiments did not indicate any improvement using 3D 0:05). Comparing to ISLES 2017 testing set, there was a patches for the U-Net, at least for our architecture. So, slight decrease in performance. This could be explained this improvement could have come from a combination by the shift on the intensity distribution of the MRI maps, of the e ect of the ensemble and the meta-data. But we due to di erent acquisition protocols or the di erences in note that their method presented a lower average Recall, the preprocessing step. which explains their lower average Dice score. Regard- ing the average Recall score, our method was fourth, but 4.3. State-of-the-art: ISLES 2017 Challenge when we consider the top 3 methods, specially Pinto et al. The results of published methods for final infarct stroke (2018a), we conclude that it was obtained with a much lesion prediction using ISLES 2017 testing set (Winzeck lower average Precision, which means that to increase the et al., 2018), together with our baseline and proposal true positive detections, they had to increase substantially methods are presented in Table 6. The metrics were com- the false positives. So, comparing with the state of the puted by the online platform, so the ground-truth data, art, our method presented a better balance between Pre- Table 5: Results obtained by our proposal and baseline method in the private dataset. Each metric is represented by the mean  standard deviation. Underlined values correspond to the highest mean, while bold values represent statistically significant values (p-value < 0:05). Supervised Block Unsupervised Block Dice HD ASSD Precision Recall FCN G-RNN – U-Net [2D] – 0.32  0.18 32.58  20.09 5.17  3.34 0.31  0.26 0.68  0.28 RBM + RBM [3D] U-Net [2D] LSTM 0.36  0.18 26.68  15.60 3.88  2.17 0.38  0.30 0.68  0.27 Lesion Haemo 13 Table 6: Published methods in ISLES 2017 testing dataset and our proposal. Each metric is represented by the mean  standard deviation. Underlined values correspond to the highest mean. Dice HD ASSD Precision Recall Mok et al. * 0.32  0.23 40.74  27.23 8.97  9.52 0.34  0.27 0.39  0.27 Kwon et al. * 0.31  0.23 45.26  21.04 7.91  7.31 0.36  0.27 0.45  0.30 Robben et al. * 0.27  0.22 37.84  17.75 6.72  4.10 0.44  0.32 0.39  0.31 Pisov et al. * 0.27  0.20 49.24  32.15 9.49  10.56 0.31  0.27 0.39  029 Monteiro et al. * 0.30  0.22 46.60  17.50 6.31  4.05 0.34  0.27 0.51  0.30 Pinto et al. (2018a) 0.29  0.21 41.58  22.04 7.69  5.71 0.21  0.21 0.66  0.29 Lucas et al. * 0.29  0.21 33.85  16.82 6.81  7.18 0.34  0.26 0.51  0.32 Choi et al. * 0.28  0.22 43.89  20.70 8.88 8.19 0.36  0.31 0.41  0.31 Niu et al. * 0.26  0.20 48.88  11.20 6.26  3.02 0.28  0.25 0.56  0.26 Sedlar et al. * 0.20  0.19 58.30  20.02 11.19  9.10 0.23  0.24 0.40  0.29 Rivera et al. * 0.19  0.16 63.58  18.58 11.13  7.89 0.27  0.25 0.21  0.17 Islam et al. * 0.19  0.18 64.15  28.51 14.17  15.80 0.29  0.28 0.25  0.25 Chengwei et al. * 0.18  0.17 65.95  25.94 9.22  6.99 0.37  0.30 0.21  0.23 Yoon et al. * 0.17  0.16 45.23  19.14 12.43  11.01 0.23  0.27 0.36  0.32 Baseline 0.30  0.21 36.58  16.62 6.96  5.08 0.30  0.26 0.55  0.31 Proposal 0.38  0.22 29.21  15.04 5.52  5.06 0.41  0.26 0.53  0.29 Methods presented in Winzeck et al. (2018), whose results were retrieved from SMIR Online Platform (2017). cision and Recall, which reflected into a higher average Dice score. Based on the results, we may conclude that the use of complementary features provided by the RBMs and the use of LSTM for a larger context allowed our baseline to surpass current state-of-the-art methods. Results from ChallengeR Benchmark. The SMIR plat- Figure 5: Boxplot of the top-10 ranking methods ordered by average form of ISLES 2017 provides a weekly benchmark report Dice score in ISLES 2017 testing set. of the current top-10 methods in the testing set, accord- ing to the average Dice score. So, some of the methods may not be published, lacking a description on their im- In Figure 6 we have the significance maps of the pair- plementation, and, for this reason, were not included in wise significant test with one-side Wilcoxon signed rank the previous discussion. test (p-value = 0:05), showing that our method reached Figure 5 presents the boxplots of each method consid- higher Dice score statistically significant when compared ered in the report. with other five ranked methods of the top-10. We observe that the top-10 methods failed to predict the lesion of one or more cases (lowest outliers), which Figure 7 shows the podium plot of each method for each may indicate the degree of complexity of predicting in- case in the testing set, and its ranking. We observe that farct stroke lesion 90 days ahead in ISLES 2017 Chal- our proposal is the method, which ranked first most of the lenge dataset. But we verify that our method is the only times, as well as second and third. Also, when we con- one to have the first quartile above 0.20 in the Dice score. sider the methods ranked bellow fourth, our method is in Single Model Ensemble art, presenting the highest average Dice score and lowest average distance score. Considering the ablation study, this performance was attained due to the combination of adding extra features obtained by encoding the parametric maps with RBMs, according to the underlining physical meaning, and the elaboration provided by the long context of the LSTM layers. 5. Conclusions Figure 6: One-side Wilcoxon signed rank test in ISLES 2017 testing set. Statistically significant tests are marked with the blue colour, while the In this work, we present a deep learning approach for red colour designates statistically non-significant tests. predicting the final stroke lesion, based on unsupervised and supervised learning. We proposed to group the input maps according to the underlying physical principle be- general among those with the lowest counts. Analysing hind their creation, namely, the time-resolved perfusion the cases individually, we note two trends, for some cases maps (Tmax, TTP, MTT), and the blood-flow-dynamic re- all methods presented similar performance, while for oth- lated maps (rCBF, rCBV). Each group was encoded using ers, we find a large variation from the first to the other an unsupervised model to obtain structural features spe- methods. The first trend may be found in the most di- cific to its underlying physical principle. These structural cult case, where all methods had zero or a close value for features together with the standard parametric maps were the Dice score. In the second trend, we observe that our fed to a supervised model to learn features conditioned method is ranked as first most of the cases. on the label, which in our problem, means to condition Based on the results of the benchmark, we may infer on the results of the medical intervention — lesion at 90- that our method is competitive among current state of the days follow-up. We also investigated the use of Gated Figure 7: Podium plot of each testing case in ISLES 2017. For each ISLES 2017 testing subject, defined by a coloured line with circles, the Podium plot orders decreasingly the Dice score obtained by each of the top 10 methods that are represented by coloured circles. 15 Recurrent Neural Networks to provide long spatial con- Sprengers, M.E., Jenniskens, S.F., Lycklama a ` Nije- text, which were critical in relating the structural features holt, G.J., et al., 2016. Collateral status on baseline to the information on input parametric maps. Our results computed tomographic angiography and intra-arterial showed that either the encoding or the long spatial context treatment e ect in patients with proximal anterior cir- improved over our baseline. Also, these two together in- culation stroke. Stroke 47, 768–776. teracted positively increasing the performance when con- Butcher, K., Emery, D., 2010a. Acute stroke imaging part sidering separately each one. i: Fundamentals. Canadian Journal of Neurological When evaluating our proposal on ISLES 2017 testing Sciences 37, 4–16. dataset, we observe a prediction improvement over cur- rent state-of-the-art methods. The proposed method ob- Butcher, K., Emery, D., 2010b. Acute stroke imaging part tained the first place in Dice and also in HD and ASSD. ii: the ischemic penumbra. Canadian Journal of Neu- Recent works (Pinto et al., 2018a; Robben et al., 2020) rological Sciences 37, 17–27. have shown the importance of clinical meta-data to pre- Choi, Y., Kwon, Y., Lee, H., Kim, B.J., Paik, M.C., dict the final stroke lesion in di erent revascularization Won, J.H., 2016. Ensemble of deep convolutional neu- scenarios. So as future work, we aim to study how such ral networks for prognosis of Ischemic Stroke, in: In- meta-data (i.e. TICI score) could be incorporated in our ternational Workshop on Brainlesion: Glioma, Mul- architecture, to consolidate the impact of the clinical in- tiple Sclerosis, Stroke and Traumatic Brain Injuries, tervention and to further improve the 90-day lesion pre- Springer. pp. 231–243. diction. Coutts, S.B., Simon, J.E., Tomanek, A.I., Barber, P.A., Chan, J., Hudon, M.E., Mitchell, J.R., Frayne, R., Acknowledgement Eliasziw, M., Buchan, A.M., et al., 2003. Reliability of assessing percentage of di usion-perfusion mismatch. Adriano Pinto was supported by a scholarship from the Stroke 34, 1681–1683. Fundac ¸ao ˜ para a Ciencia ˆ e Tecnologia (FCT), Portugal (scholarship number PD/BD/113968/2015). This work El Tawil, S., Muir, K.W., 2017. Thrombolysis and was supported by FCT national funds, under the national thrombectomy for acute ischaemic stroke. Clinical support to R&D units grant, through the reference project Medicine 17, 161–165. UIDB/04436/2020 and UIDP/04436/2020. Gonzalez, R., Hirsch, J., Koroshetz, W., Lev, M., Schae- fer, P., 2007. Acute ischemic stroke: imaging and in- tervention. American Journal of Neuroradiology 28, References Badrinarayanan, V., Kendall, A., Cipolla, R., 2015. Gonzalez, ´ R.G., Hirsch, J.A., Koroshetz, W., Lev, M.H., Segnet: A deep convolutional encoder-decoder ar- Schaefer, P.W., 2011. Acute ischemic stroke. Springer. chitecture for image segmentation. arXiv preprint arXiv:1511.00561 . Grysiewicz, R.A., Thomas, K., Pandey, D.K., 2008. Epi- demiology of ischemic and hemorrhagic stroke: inci- Bauer, S., Gratz, P.P., Gralla, J., Reyes, M., Wiest, R., dence, prevalence, mortality, and risk factors. Neuro- 2014. Towards automatic MRI volumetry for treat- logic clinics 26, 871–895. ment selection in acute ischemic stroke patients, in: En- Higashida, R.T., Furlan, A.J., Roberts, H., Tomsick, T., gineering in Medicine and Biology Society (EMBC), Connors, B., Barr, J., Dillon, W., Warach, S., Broder- 2014 36th Annual International Conference of the ick, J., Tilley, B., et al., 2003. Trial design and report- IEEE, IEEE. pp. 1521–1524. ing standards for intraarterial cerebral thrombolysis for Berkhemer, O.A., Jansen, I.G., Beumer, D., Fransen, acute ischemic stroke. Journal of Vascular and Inter- P.S., Van Den Berg, L.A., Yoo, A.J., Lingsma, H.F., ventional Radiology 14, E1–E31. 16 Hinton, G.E., 2012. A practical guide to training restricted Bentley, P., Chen, L., et al., 2017. ISLES 2015-a Boltzmann machines, in: Neural networks: Tricks of public evaluation benchmark for ischemic stroke lesion the trade. Springer, pp. 599–619. segmentation from multispectral MRI. Medical image analysis 35, 250–269. Hochreiter, S., Schmidhuber, J., 1997. Long short-term memory. Neural computation 9, 1735–1780. McKinley, R., Hani, ¨ L., Gralla, J., El-Koussy, M., Bauer, S., Arnold, M., Fischer, U., Jung, S., Mattmann, K., Jenkinson, M., Beckmann, C.F., Behrens, T.E., Woolrich, Reyes, M., et al., 2017. Fully automated stroke tis- M.W., Smith, S.M., 2012. Fsl. Neuroimage 62, 782– sue estimation using random forest classifiers (faster). Journal of Cerebral Blood Flow & Metabolism 37, 2728–2741. Kamnitsas, K., Ledig, C., Newcombe, V.F., Simpson, J.P., Kane, A.D., Menon, D.K., Rueckert, D., Glocker, B., Memezawa, H., Smith, M.L., Siesjo, ¨ B.K., 1992. Penum- 2017. Ecient multi-scale 3d CNN with fully con- bral tissues salvaged by reperfusion following middle nected CRF for accurate brain lesion segmentation. cerebral artery occlusion in rats. Stroke 23, 552–559. Medical image analysis 36, 61–78. Menon, B.K., Qazi, E., Nambiar, V., Foster, L.D., Yeatts, Kemmling, A., Flottmann, F., Forkert, N.D., Minnerup, S.D., Liebeskind, D., Jovin, T.G., Goyal, M., Hill, J., Heindel, W., Thomalla, G., Eckert, B., Knauth, M.D., Tomsick, T.A., et al., 2015. Di erential e ect M., Psychogios, M., Langner, S., Fiehler, J., 2015. of baseline computed tomographic angiography collat- Multivariate dynamic prediction of ischemic infarction erals on clinical outcome in patients enrolled in the in- and tissue salvage as a function of time and degree terventional management of stroke iii trial. Stroke 46, of recanalization. Journal of Cerebral Blood Flow & 1239–1244. Metabolism 35, 1397–1405. Milletari, F., Navab, N., Ahmadi, S.A., 2016. V-net: Fully Kistler, M., Bonaretti, S., Pfahrer, M., Niklaus, R., convolutional neural networks for volumetric medical Buchler ¨ , P., 2013. The virtual skeleton database: an image segmentation, in: 3D Vision (3DV), 2016 Fourth open access repository for biomedical research and col- International Conference on, IEEE. pp. 565–571. laboration. Journal of medical Internet research 15. Mok, T.C., Chung, A.C., 2017. Deep adversarial net- Krizhevsky, A., Sutskever, I., Hinton, G.E., 2012. Ima- works for stroke lesion segmentation . genet classification with deep convolutional neural net- Monteiro, M., Oliveira, A.L., 2017. Fully convolutional works, in: Advances in neural information processing neural network for 3d stroke lesion segmentation . systems, pp. 1097–1105. Nair, V., Hinton, G.E., 2010. Rectified linear units im- Labeyrie, M.A., Turc, G., Hess, A., Hervo, P., Mas, J.L., prove restricted Boltzmann machines, in: Proceedings Meder, J.F., Baron, J.C., Touze, E., Oppenheim, C., of the 27th international conference on machine learn- 2012. Di usion lesion reversal after thrombolysis: a ing (ICML-10), pp. 807–814. mr correlate of early neurological improvement. Stroke 43, 2986–2991. Nielsen, A., Hansen, M.B., Tietze, A., Mouridsen, K., 2018. Prediction of tissue outcome and assessment Liebeskind, D.S., 2003. Collateral circulation. Stroke 34, of treatment e ect in acute ischemic stroke using deep 2279–2284. learning. Stroke 49, 1394–1401. Lucas, C., Heinrich, M.P., 2017. 2d multi-scale res-net Niu, Y., Gong, E., Xu, J., Pauly, J., Zaharchuk, G., for stroke segmentation. 2018. Improved prediction of the final infarct from Maier, O., Menze, B.H., von der Gablentz, J., Hani, ¨ L., acute stroke neuroimaging using deep learning, in: Heinrich, M.P., Liebrand, M., Winzeck, S., Basit, A., STROKE. 17 Oliveira, A., Pereira, S., Silva, C.A., 2018. Retinal vessel Robben, D., Boers, A.M., Marquering, H.A., Langezaal, segmentation based on fully convolutional neural net- L.L., Roos, Y.B., van Oostenbrugge, R.J., van Zwam, works. Expert Systems with Applications 112, 229– W.H., Dippel, D.W., Majoie, C.B., van der Lugt, A., 242. et al., 2020. Prediction of final infarct volume from native ct perfusion and treatment parameters using deep Pereira, S., Meier, R., McKinley, R., Wiest, R., Alves, learning. Medical image analysis 59, 101589. V., Silva, C.A., Reyes, M., 2018. Enhancing inter- Robben, D., Suetens, P., 2017. Dual-scale fully convo- pretability of automatically extracted machine learning lutional neural network for final infarct prediction, in: features: application to a RBM-Random Forest system Ischemic stroke lesion segmentation-ISLES challenge on brain lesion segmentation. Medical image analysis 2017, held in conjunction with MICCAI 2017, Date: 44, 228–244. 2017/09/10-2017/09/10, Location: Quebec City, Que- Pereira, S., Pinto, A., Amorim, J., Ribeiro, A., Alves, bec, Canada. V., Silva, C.A., 2019. Adaptive feature recombination Ronneberger, O., Fischer, P., Brox, T., 2015. U-net: and recalibration for semantic segmentation with fully Convolutional networks for biomedical image seg- convolutional networks. IEEE Transactions on Medi- mentation, in: International Conference on Medical cal Imaging 38, 2914–2925. doi:10.1109/TMI.2019. image computing and computer-assisted intervention, Springer. pp. 234–241. Pinto, A., McKinley, R., Alves, V., Wiest, R., Silva, C.A., Rose, S.E., Chalk, J.B., Grin, M.P., Janke, A.L., Chen, Reyes, M., et al., 2018a. Stroke lesion outcome pre- F., McLachan, G.J., Peel, D., Zelaya, F.O., Markus, diction based on MRI imaging combined with clinical H.S., Jones, D.K., et al., 2001. MRI based di usion information. Frontiers in Neurology 9, 1060. and perfusion predictive model to estimate stroke evo- lution. Magnetic resonance imaging 19, 1043–1053. Pinto, A., Pereira, S., Meier, R., Alves, V., Wiest, R., Silva, C.A., Reyes, M., 2018b. Enhancing clinical MRI Rumelhart, D.E., McClelland, J.L., 1986. Parallel dis- perfusion maps with data-driven maps of complemen- tributed processing: explorations in the microstructure tary nature for lesion outcome prediction, in: Medical of cognition. volume 1. foundations . Image Computing and Computer Assisted Intervention Scalzo, F., Hao, Q., Alger, J.R., Hu, X., Liebeskind, D.S., – MICCAI 2018, pp. 107–115. 2012. Regional prediction of tissue fate in acute is- chemic stroke. Annals of Biomedical Engineering 40, Pisov, M., Belyaev, M., Krivov, E., 2017. Neural networks 2177–2187. ensembles for ischemic stroke lesion segmentation . SMIR Online Platform, 2017. Ischemic stroke lesion Powers, W.J., Rabinstein, A.A., Ackerson, T., Adeoye, segmentation 2017. https://www.smir.ch/ISLES/ O.M., Bambakidis, N.C., Becker, K., Biller, J., Brown, Start2017. [Acessed: 2019-11-18]. M., Demaerschalk, B.M., Hoh, B., et al., 2018. 2018 guidelines for the early management of patients Vinh, N.X., Epps, J., Bailey, J., 2010. Information the- with acute ischemic stroke: a guideline for health- oretic measures for clusterings comparison: Variants, care professionals from the American Heart Associa- properties, normalization and correction for chance. tion/American Stroke Association. Stroke 49, e46–e99. Journal of Machine Learning Research 11, 2837–2854. Rekik, I., Allassonniere, ` S., Carpenter, T.K., Ward- Visin, F., Ciccone, M., Romero, A., Kastner, K., Cho, K., law, J.M., 2012. Medical image analysis methods Bengio, Y., Matteucci, M., Courville, A., 2016. Re- in MR/CT-imaged acute-subacute ischemic stroke le- seg: A recurrent neural network-based model for se- sion: Segmentation, prediction and insights into dy- mantic segmentation, in: Proceedings of the IEEE Con- namic evolution simulation models. a critical appraisal. ference on Computer Vision and Pattern Recognition NeuroImage: Clinical 1, 164–178. Workshops, pp. 41–48. 18 Winzeck, S., Hakim, A., McKinley, R., Pinto, J.A.A.D.S., Alves, V., Silva, C., Pisov, M., Krivov, E., Belyaev, M., Monteiro, M., et al., 2018. ISLES 2016 & 2017- benchmarking ischemic stroke lesion outcome predic- tion based on multispectral MRI. Frontiers in Neurol- ogy 9, 679. World Health Organization, et al., 2014. Global status re- port on noncommunicable diseases 2014. World Health Organization. Zeiler, M.D., Fergus, R., 2014. Visualizing and under- standing convolutional networks, in: European confer- ence on computer vision, Springer. pp. 818–833. Zivelonghi, C., Tamburin, S., 2018. Mechanical thrombectomy for acute ischemic stroke: the therapeu- tic window is larger but still “time is brain”. Functional neurology 33, 5. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Computing Research Repository arXiv (Cornell University)

Combining unsupervised and supervised learning for predicting the final stroke lesion

Download PDF
19 pages

Loading next page...
 
/lp/arxiv-cornell-university/combining-unsupervised-and-supervised-learning-for-predicting-the-ZOk42aCxsC

References (68)

ISSN
1361-8415
eISSN
ARCH-3344
DOI
10.1016/j.media.2020.101888
Publisher site
See Article on Publisher Site

Abstract

Predicting the final ischaemic stroke lesion provides crucial information regarding the volume of salvageable hypoperfused tissue, which helps physicians in the dicult decision-making process of treatment planning and intervention. Treatment selection is in- fluenced by clinical diagnosis, which requires delineating the stroke lesion, as well as characterising cerebral blood flow dynamics using neuroimaging acquisitions. Nonetheless, predicting the final stroke lesion is an intricate task, due to the variability in lesion size, shape, location and the underlying cerebral haemodynamic processes that occur after the ischaemic stroke takes place. More- over, since elapsed time between stroke and treatment is related to the loss of brain tissue, assessing and predicting the final stroke lesion needs to be performed in a short period of time, which makes the task even more complex. Therefore, there is a need for automatic methods that predict the final stroke lesion and support physicians in the treatment decision process. We propose a fully automatic deep learning method based on unsupervised and supervised learning to predict the final stroke lesion after 90 days. Our aim is to predict the final stroke lesion location and extent, taking into account the underlying cerebral blood flow dynamics that can influence the prediction. To achieve this, we propose a two-branch Restricted Boltzmann Machine, which provides specialized data-driven features from di erent sets of standard parametric Magnetic Resonance Imaging maps. These data-driven feature maps are then combined with the parametric Magnetic Resonance Imaging maps, and fed to a Convolutional and Recurrent Neural Net- work architecture. We evaluated our proposal on the publicly available ISLES 2017 testing dataset, reaching a Dice score of 0.38, Hausdor Distance of 29.21 mm, and Average Symmetric Surface Distance of 5.52 mm. Keywords: Deep Learning, Image Prediction, Magnetic Resonance Imaging, Stroke 1. Introduction caused by thrombolysis, haemodynamic factors, or em- bolic causes (Grysiewicz et al., 2008). Due to vessel oc- Stroke is the second leading cause of death worldwide clusion, the insucient supply of oxygenated blood to (World Health Organization et al., 2014), being classified brain cells leads to hypoperfused brain tissue, trigger- in two types: ischaemic and haemorrhagic (Grysiewicz ing cellular mechanisms to preserve the integrity of the et al., 2008). Ischaemic stroke is the most common type, cell. The hypoperfused area consists of tissue at risk that resulting from an occlusion of a vessel, which can be can be salvaged, being designated penumbra. As time passes, in the absence of flow restoration or sucient col- lateral blood flow supply, the hypoperfused tissue even- Corresponding author: Department of Industrial Electronics, Cam- tually reaches a non-salvageable state designated core or pus Azurem, ´ Guimaraes, ˜ Portugal. infarct tissue (Memezawa et al., 1992). Email addresses: id6376@alunos.uminho.pt (Adriano Pinto ), csilva@dei.uminho.pt (Carlos A. Silva) Diagnosis and treatment of ischaemic stroke relies on Preprint submitted to Elsevier January 5, 2021 arXiv:2101.00489v1 [eess.IV] 2 Jan 2021 neuroimaging acquisitions, where Computed Tomogra- two-pathway architecture, trained with two subsets of phy (CT) and Magnetic Resonance Imaging (MRI) are the standard parametric MRI maps. One subset encompasses preferred imaging modalities (Gonzalez et al., 2007). CT the Time-To-Peak (TTP), Mean Transit Time (MTT), imaging remains the most used acquisition due to its ra- Time-to-Maximum (Tmax), and Apparent Di usion Co- pidity and availability (Gonzalez et al., 2007). However, ecient (ADC). The second set contains the ADC, the multi-parametric MRI provides a higher sensitivity in de- relative Cerebral Blood Volume (rCBV), and the relative tecting early ischaemic stroke and assessing the penum- Cerebral Blood Flow (rCBF). In a second stage, the fea- bra region (Gonzalez et al., 2007). Treatment consists in ture maps computed by the RBMs are combined with the restoring tissue perfusion levels, also known as reperfu- standard parametric MRI maps to form the input of a su- sion, by performing mechanical thrombectomy or throm- pervised deep learning architecture composed by Convo- bolysis. Since ischaemic stroke is a dynamic process that lutional Neural Networks (CNNs) and Recurrent Neural evolves over time, the treatment is only possible up to 24 Networks (RNNs). The proposed architecture was evalu- hours, where viable neurones still persist (El Tawil and ated using the publicly-available ISLES 2017 dataset. Muir, 2017; Zivelonghi and Tamburin, 2018). So, expert 1.1. On the complexity of predicting the final infarct physicians must evaluate the benefits and risks of mechan- stroke lesion ical thrombectomy before an intervention, since it may cause haemorrhage, vascular injury, and other complica- In acute ischaemic stroke, the clinical evaluation of the tions (Powers et al., 2018). If performed, the success of standard parametric maps (e.g. ADC and Tmax) can iden- the intervention is assessed radiologically via angiogra- tify infarct tissue and tissue that will infarct in the ab- phy imaging and scored by a qualitative expert-generated sence of therapeutic intervention. In this analysis, the in- scale designated the standardized Thrombolysis in Cere- farct tissue, is identified by the hypointense regions of the bral Infarction (TICI) scale (Higashida et al., 2003). Dur- ADC map, which characterize tissue with limited di u- ing the decision-making process, the physician needs to sion (Butcher and Emery, 2010a). Hypoperfused tissue, assess the nature and location of the lesion alongside i.e. tissue that will infarct, is identified by hyperintense pathophysiological factors such as age, presence of co- regions of the Tmax map, indicating an increased arrival morbidities, and collateral circulation (Liebeskind, 2003). time of contrast agent (Butcher and Emery, 2010b). How- The latter is of utmost importance in ischaemic stroke. ever, to correctly predict the final ischaemic stroke lesion, The presence of collateral circulation, where a secondary besides considering the complex time-evolving transfor- network of vessels is responsible for granting cerebral mation of hypoperfused tissue to infarcted tissue, it is blood flow to the lesioned tissue, increases the chances also necessary to appraise the impact of the clinical inter- of a successful reperfusion (Liebeskind, 2003). Asserting vention, thrombectomy, on the underlying brain perfusion the potential ecacy of treatment can be time-consuming and di usion. and prone to inter- and intra-variability among physicians, A successful thrombectomy should restore the perfu- which is further potentiated when performed in a clinical sion levels, recovering the hypoperfused tissue. How- emergency environment (Coutts et al., 2003). Moreover, ever, several factors may a ect the reperfusion, limiting since time is critical, MRI acquisitions are optimized for the degree of success of the intervention. To better un- speed, which is accomplished by reducing the resolution derstand the nuances of the clinical intervention, consider (Gonzalez ´ et al., 2011), making the prediction of the final the cases presented in Figure 1. In the first case, Figure stroke lesion an intricate task. Thus, automatic prediction 1a, the ADC does not present any hypointense region, so of a stroke lesion at a given time since stroke has a great no infarct tissue may be identified, and we should expect potential to guide physicians in this time-critical decision- a complete recovery of the hypoperfused tissue indicated making process. by Tmax; however, the follow-up delineation obtained af- We propose a novel automatic method based on unsu- ter thrombectomy presents a large final lesion, which is pervised and supervised deep learning. We utilize Re- explained by an unsuccessful intervention. In the second stricted Boltzmann Machines (RBMs) to jointly charac- case, Figure 1b, we observe a final infarct lesion that is terise the lesion and blood flow information through a smaller than the hypointense region present in the ADC 2 (a) (b) Figure 1: ADC and Tmax parametric maps of two patient cases from ISLES 2017 training set, and the final lesion delineated at a 90-day follow- up, overlapped with the onset ADC: patient 0036 (Figure 1a) with an unsuccessful reperfusion, and patient 0006 (Figure 1b), where the clinical intervention was successful. (Figure 1b arrow). This indicates reversible di usion re- Challenge in 2016 and 2017, new methods have been pro- striction, which is a rare case (Labeyrie et al., 2012) and posed. These aim to predict at a 90-day time-window. was only possible to identify by a follow-up T2-weighted Rose et al. (2001) proposed a two-stage method based acquisition. So, an automatic method for predicting the fi- on parametric perfusion and di usion MRI maps. On the nal stroke lesion has not only to capture the time-evolving first stage, the method defines a region of interest (ROI) process of di usion and perfusion, but also to consider di- based on the intensity signal of the standard parametric rectly or indirectly the degree of success of the thrombec- maps, the MTT, Cerebral Blood Flow (CBF), Cerebral tomy, which may condition the final lesion either to be Blood Volume (CBV), and Di usion-Weighted Imaging confined to the hypointense region of the ADC map, or (DWI). The second stage performs stroke tissue predic- to grow to brain tissue areas that are hyperintense in the tion using Gaussian mixture models trained in di erent Tmax. Due to the time-evolving process of di usion and sets of parametric maps. Bauer et al. (2014) used Random perfusion, the complexity of predicting the lesion will ag- Forests to segment or predict the final stroke lesion de- gravate as we move from a target window of some days to pending on whether acute stroke imaging or three-month several months. follow-up imaging was available, respectively. McKin- The complexity of the evaluation process may be also ley et al. (2017) also used a two-stage classification ap- observed in the inter-rater agreement of expert radiolo- proach as in Rose et al. (2001) for lesion characterisa- gists in ISLES 2017 dataset, which obtained a Dice score tion and lesion prediction, where each stage consists of of 0:58  0:20 on delineating the lesion using a 90-day two sets of Random Forests (RFs) classifiers. The first follow-up T2-weighted acquisition (Winzeck et al., 2018). stage aims to define a ROI that encompasses the hypop- erfused region. In the first set, each classifier is trained with features extracted from di erent sets of MRI para- 1.2. Previous Work metric maps. Having defined the location and extension Contrary to stroke lesion segmentation, where several of the lesion, a second set of two RFs performs stroke methods have already been proposed (Rekik et al., 2012; tissue prediction. Such classifiers were trained on dif- Maier et al., 2017), the complexity of predicting the fi- ferent sets of patients, stratified by the TICI score. One nal stroke lesion has only recently attracted attention in classifier is trained in patients with unsuccessful reperfu- the medical imaging community. For predicting the final sion interventions, whereas a second classifier is trained stroke lesion several methods have been already proposed in patients with successful reperfusion. The final predic- based on multivariate linear regression models (Scalzo tion is obtained by combining the results of both clas- et al., 2012; Rose et al., 2001; Kemmling et al., 2015), sifiers, using a logistic regression model. Scalzo et al. decision trees (McKinley et al., 2017; Bauer et al., 2014), (2012) proposed a framework for stroke tissue prediction, and CNNs (Choi et al., 2016). Furthermore, with the re- which characterises the state of the lesion four days af- lease of Ischaemic Stroke LEsion Segmentation (ISLES) ter clinical intervention (thrombectomy). From the Fluid 3 Attenuation Inversion Recovery (FLAIR) MRI sequence, olution, while in the second branch the input resolution ADC and Tmax MRI maps, the method applies a regres- was lowered by a factor of 3. The output of each branch sion model that learns the behaviour of neighbouring vox- is transformed to the same scale and merged by two fully els within a cuboid. Kemmling et al. (2015) proposed connected layers. The network is trained with four dif- a multi-modality approach based on CT and MRI maps ferent sets of hyper-parameters. These four networks are with non-imaging clinical meta-data, namely the TICI used as an ensemble, whose prediction is obtained by av- score and the time to treatment of each patient, to perform eraging the output of each one. Similarly, Niu et al. (2018) stroke tissue prediction. used multiple scales of overlapping 3D patches to capture In another line of research, authors have investigated local and global spatial information. In the review pa- the use of deep learning (Choi et al., 2016; Nielsen et al., per of Winzeck et al. (2018), Rivera et al. also built on 2018; Robben et al., 2020) for stroke tissue prediction. the work of Kamnitsas et al. (2017) and Milletari et al. Choi et al. (2016), the winner approach at ISLES 2016 (2016), by proposing a scheme to extract di erent patch Challenge, proposed an ensemble of twelve CNN archi- resolutions, independent of each other, that are fed into tectures, grouped into two sets of networks. The first four di erent paths. Afterwards, a fully connected layer group comprehends four 3D U-Nets (Ronneberger et al., combines all the outputs to perform stroke tissue predic- 2015) performing voxel-wise tissue prediction. The sec- tion. Pisov et al. (2017) employed an ensemble strategy ond group of networks uses two-pathway Fully Con- by combining di erent CNN-based architectures to over- nected Networks (FCNs) performing two types of patch- come the strong anisotropy of the data. As summarized by wise classification. One set of FCNs classifies a patch Winzeck et al. (2018), the work of Yoon et al. proposed a as lesion if it includes any lesion voxel. The other set two-stage gated CNN. In a first stage, the authors perform of FCNs classifies a patch as lesion if the central voxel lesion detection and delineation. Afterwards, based on the is a lesion. After merging the two pathway FCN, the probability maps of the first stage, a second CNN archi- method incorporates meta-data by adding a dense layer tecture processes the regions where the probability maps of clinical predictors merged with the imaging output of of healthy tissue and lesion are close to each other. Pinto each network. The final stroke lesion prediction results et al. (2018b) made use of temporal perfusion imaging, from a weighted merging of all models. Mok and Chung the Dynamic Susceptibility Contrast-MRI, in a U-Net ar- (2017) applied deep adversarial training for stroke tis- chitecture. This architecture aims to temporally process sue prediction in an ensemble of U-Nets. Monteiro and and extract deep features, which are then combined with Oliveira (2017) proposed a method based on the V-Net a second feature step of another U-Net network, which architecture (Milletari et al., 2016). The training was con- was trained on the standard parametric maps. Using a ducted with a custom loss function that applies a weighted large CT dataset, Robben et al. (2020) predicted the fi- sum between Dice score and cross entropy. Lucas and nal infarct stroke lesion with a temporal window ranging Heinrich (2017) proposed the use of a U-Net architec- from 24 hours to 5 days. The authors considered spatio- ture, which combines patches from the MRI maps in the temporal CT perfusion as input to a deep neural network same slice, with patches from 3 neighbouring slices and inspired in the architecture proposed by Kamnitsas et al. 2 hemispheric flips. In the expanding path of the U-Net, (2017). Additionally, the model combines CT neuroimag- each level computes a Dice loss for the healthy tissue and ing with clinical meta-data. Nielsen et al. (2018) pro- for the final lesion, after the softmax activation. After- posed a method based on the SegNet architecture (Badri- wards, all losses are summed up, having the loss of the narayanan et al., 2015), predicting on a 30-day follow-up lesion and healthy tissue weighted according to a prior acquisition based on a private dataset. probability (Winzeck et al., 2018). Robben and Suetens Principal and collateral blood flow has been consid- (2017) employed a CNN-based architecture inspired by ered either directly by modelling the temporal perfusion Kamnitsas et al. (2017). The authors proposed to com- imaging (Pinto et al., 2018b), or indirectly by perfusion bine the MRI inputs with clinical meta-data, before feed- and di usion parametric maps (Choi et al., 2016; Maier ing them to each branch of a two-pathway 3D network. et al., 2017; Scalzo et al., 2012), or through clinical in- In the first branch the input is kept with the original res- formation that characterises the success of the revascu- 4 larization (McKinley et al., 2017). We hypothesize that 2. Methods modelling the haemodynamics of the brain when artery In this work, predicting the final infarct stroke lesion occlusion occurs can be beneficial for predicting the final consists of delineating the lesion’s spatial extent at a 90- stroke lesion. So, in this work, we investigate the rep- day follow-up time-point, using multi-parametric MRI resentation of the haemodynamics through an unsuper- imaging, namely the ADC, MTT, TTP, Tmax, rCBF, and vised learning model. Contrary to previous approaches, rCBV, which are acquired at the onset time-point. The we propose grouping the input maps according to their architecture of the proposed system and its main compo- subjacent physical meaning and encoding each group sep- nents are described in the following subsections. arately with an RBM. As groups, we investigated the time-resolved perfusion maps (Tmax, TTP, MTT), and the 2.1. Architecture blood-flow-dynamic related maps (rCBF, rCBV) (Butcher and Emery, 2010a,b). Our proposal of combining features The overall architecture of the proposed method can be obtained unsupervisedly and supervisedly was motivated divided into two functional blocks, as shown in Figure 2. by the knowledge that unsupervised models learn struc- The first functional block performs unsupervised rep- tural features of the original image, while the supervised resentation learning using two unsupervised models, models learn features conditioned on the label, so there namely RBMs. This unsupervised block provides new is potential for obtaining richer and more discriminative features that represent structural information that comple- features by joining both types of models. ments the standard parametric MRI maps, enhancing the capacity of our model to predict the final infarct lesion. In our approach, we aim to model the clinical procedure, 1.3. Contributions which first locates and delineates the lesion at current time, and then considers the blood flow haemodynamic This work presents an automatic approach for predict- that might influence the final stroke lesion prediction. ing the final stroke lesion, using onset neuroimaging data. This procedure is encoded in our two-path RBM. The The main contributions are: first RBM is responsible for capturing information on le- sion location and extension, referred to as the RBM . - The use of unsupervised models for extracting Lesion The second RBM, RBM , aims to capture blood flow structural features of time-resolved perfusion and Haemo haemodynamics information (e.g. collateral circulation), blood-flow-dynamic related MRI maps for predict- which has been identified as a key factor by physicians ing stroke lesion. when assessing stroke final infarct lesion in clinical re- - The use of local and long spatial context provided ports (Berkhemer et al., 2016; Menon et al., 2015). On by gated recurrent neural networks for relating struc- one hand, to locate the onset ischaemic stroke lesion, the tural features and image information when learn- RBM considers standard parametric maps that char- Lesion ing features conditioned on the label in a supervised acterise the arrival times and mean transit times of the model. contrast agent. In the presence of an ischaemic lesion, the occluded vessel can decrease or interrupt the normal brain - The proposal of a competitive system which outper- perfusion, translating into hyperintense regions on time- forms state-of-the-art methods to predict the final in- related parametric maps (Butcher and Emery, 2010b). On farct stroke lesion, in ISLES 2017 Challenge dataset. the other hand, the RBM considers standard para- Haemo metric maps that characterise the amount of blood being The remainder of the paper is organized as follows. delivered in unit of time, which correlates to the cerebral Section 2 describes the fundamental components of the blood flow haemodynamics (Butcher and Emery, 2010b). proposed method. Section 3 describes the dataset, the Thus, the RBM considers the MTT, TTP and Tmax Lesion evaluation procedure and the setup. The results and the perfusion maps, while the RBM the rCBV and rCBF Haemo discussion are addressed in Section 4. Finally, in Section perfusion maps. Regarding the ADC standard di usion 5 we present the main conclusions. map, it is present in both RBM and RBM , since Lesion Haemo 5 Figure 2: Overview of the proposed method for predicting the final stroke lesion. In the supervised learning block, the input data dimensions are defined for each operation. it provides higher brain structural information and allows data into a feature vector is performed through the interac- the identification of tissue that is already infarcted. This tion of states between the visible and hidden units, which separation of the input imaging allows the RBM to learn is learned by minimizing an energy function. specific feature sets, which may enable the method to The complete pipeline of the unsupervised block is analyse dicult cases where information concerning the shown in Figure 3 and detailed in Section 3.4. The blood flow can have a favourable impact on the lesion pre- RBM and RBM function as feature generators Lesion Haemo diction. that output two complementary sets of feature maps N The second functional block consists of a deep learning and N . These features characterise the structure of the architecture that comprehends 2D convolutional blocks in images; however, we are interested only on the most dis- a U-Net-based structure, alongside recurrent blocks. As tinctive details. So, after training the RBMs, we perform imaging input data, we combine the standard parametric feature selection to reduce the generated feature space, maps with feature maps from each RBM, totalling 18 in- obtaining smaller but representative feature sets M and put feature maps. M , such thatjM j  jN j, for i 2 [1; 2], where the opera- 2 i i tor j:j denotes the cardinality of a set. In the feature selec- 2.2. Restricted Boltzmann Machines tion, we would like to select the features from the RBM that encodes the MRI maps, but also that correlates with The RBM is an undirected graphical model consti- the stroke prediction. Since the RBM is an unsupervised tuted by two layers of nodes: a visible layer and a hid- method, we compute the Normalized Mutual Information den layer (Rumelhart and McClelland, 1986). Each node to quantify the statistical dependence between each gener- has a weighted connection to all nodes in the other layer ated feature and the respective input MRI map, as defined (Rumelhart and McClelland, 1986). However, there are by Equation 1 (Vinh et al., 2010): no connections among nodes of the same layer. Orig- inally, Rumelhart and McClelland (1986) proposed the RBM to learn from binary data on both layers. How- M I(MRI ; Feat ) x y N M I (MRI ; Feat ) = 2 ; (1) sum x y ever, this does not represent well continuous real-valued H(MRI ) + H(Feat ) x y input data, which is the case of MRI data. Therefore, we model the visible nodes as linear units with independent where M I(:) is the mutual information between an MRI Gaussian noise. The hidden nodes are modelled as Noisy parametric map, MRI , and an output feature, Feat ; H(:) x y Rectifier Linear Units (NReLU), since they have been re- defines the entropy of a map, namely, MRI and Feat . x y ported to be suitable for feature extraction (Hinton, 2012). To relate the features of the RBM with the class label, This kind of RBM was previously used in segmentation we could use a classifier supervisedly trained. Since the tasks, such as in Pereira et al. (2019). Mapping the input neural network is trained iteratively, we use a RF clas- 6 Figure 3: Overview of the proposed unsupervised learning block. For each RBM of the unsupervised learning block, the selected features were M = M = 6. 1 2 sifer trained with the Mean Decrease Impurity (MDI) as decoder deep CNNs provide high levels of abstraction a surrogate to make the feature selection tractable. After, from the input data, increasing the global notion of con- we compute the MI and and MDI , normalize the text as the network grows deeper. However, it comes RBM RF MI by the maximum value, add both ranks and sort at a cost of a high receptive field (Zeiler and Fergus, RBM decreasingly. The best set will be the first M features. 2014). Thus, we used a 2D architecture in the plane Our selection method was inspired on the work of Pereira with the highest resolution, since the acquisition resolu- et al. (2018); however, their method cannot be directly tion is anisotropic in the dataset. Also, in the end of applied, since it would generate too many features for our the decoding path we expanded our learning block with problem. Gated RNNs. Due to their nature, Gated RNNs can cap- ture short- and long-term spatial relations, by retaining in- formation from previous nodes encoded in the time-steps. 2.3. Convolutional and Recurrent Neural Networks Hence, Gated-RNNs consider information from all previ- Our supervised functional block is based on the U-Net ous nodes when analysing the current one. This property, architecture as proposed by Ronneberger et al. (2015). when applied to imaging data, allows considering intra- The input of the U-Net considers the concatenation of slice contextual dependencies. In our work, we used a par- standard parametric maps with the sets of feature maps ticular Gated-RNN, namely the Long-Short Term Mem- extracted from the unsupervised block. In the first level ory (LSTM) (Hochreiter and Schmidhuber, 1997). How- of our encoder architecture we use four 2D convolutional ever, the LSTM was intrinsically developed to process blocks with kernel size of 3  3 and 32 channels. Af- 1D data (Hochreiter and Schmidhuber, 1997) (e.g. time- terwards, the output of the final convolutional block is series). To be applicable to 2D data, we developed an down-sampled by a factor of 2, starting the second encod- online 2D Partition layer that transforms a grid-structure ing level formed by two convolutional blocks with equal input (e.g. an image) into a one-dimensional sequence. In- kernel size but doubling the number of feature maps. The spired by Visin et al. (2016), the 2D Partition layer was third level of encoding follows the same pattern. The de- predefined with a neighbourhood of 2  2, where each coder level mimics the encoder counterpart. As in Ron- time-step is characterised by a feature space of four vox- neberger et al. (2015) we only used long skip connec- els. After, two Bidirectional LSTM layers are employed tions among encoder and decoder levels. These encoder- 7 along the left-right and frontal-dorsal directions followed 3.2. Evaluation Metrics by an up-sampling layer. These four layers, referred as We evaluated our proposal with five metrics, which are the Gated Recurrent block, are shown in Figure 2. In our the same ones computed by the online ISLES 2017 bench- supervised functional block, two Gated Recurrent blocks mark platform: Dice Similarity Score, Hausdor Distance were used, where the Bidirectional LSTMs have 64 and (HD), Average Symmetric Surface Distance (ASSD), Pre- 32 hidden layers, respectively. The impact of the main cision, and Recall (Kistler et al., 2013). components is evaluated in an ablation study in the exper- Dice score measures the spatial overlap between two iments. volumes. HD corresponds to the highest distance between surface points of di erent volumes, which characterise spatial outliers in the prediction. ASSD quantifies the 3. Experimental Setup average distances between the volumes’ surface. Preci- sion quantifies the proportions of correctly classified cases We evaluated the proposed approach on the publicly within a class, while Recall corresponds to the proportion available ISLES 2017 dataset and on a private dataset. of positive cases correctly identified as such. ISLES has an online benchmark platform (Kistler et al., 2013) that performs automatic evaluation (SMIR On- 3.3. Image pre- and post-processing line Platform, 2017). In this section we describe the Since MRI acquisitions were acquired from di erent dataset, the training and evaluation, and the main hyper- centers and configurations (Winzeck et al., 2018), for parameters of our method. each patient we resized all maps to a common volume of dimension of 256  256  32. Afterwards, the ADC 3.1. Data 6 2 maps were clipped between [0; 2600]  10 mm =s and ISLES 2017 dataset encompasses 75 ischaemic stroke the Tmax maps were clipped to [0; 20s], since values be- patients, which are separated into two sets: training (n = yond these ranges are known to be biologically meaning- 43) and testing (n = 32). Both sets have patients that un- less (McKinley et al., 2017). Finally, a linear scaling was derwent mechanical thrombectomy. Each patient is char- applied across all maps, to the range [0; 255]. The im- acterised by six 3D parametric MRI maps: di usion ADC ages are resized to its original size, after we perform the map, perfusion rCBF, rCBV, TTP, MTT and Tmax maps. prediction. In addition to the standard parametric maps, each case is We applied a morphological filtering as post- also characterised by a manual delineation of the lesion. processing, but since the final stroke lesion presents This refers to the 90-day stroke lesion delineated with ac- a wide range of lesion volumes (Winzeck et al., 2018), cess to the follow-up T2-weighted acquisition. However, we removed only small connected components with less the manual delineation is only available for the training than 250 voxels. This step was kept fixed for all the set, while the follow-up T2-weighted imaging is not dis- evaluated models. closed for any set. All parametric MRI maps are already 3.3.1. Data Augmentation co-registered and skull-stripped (Winzeck et al., 2018). Figure 4 (top row) shows an example of MRI maps, Data augmentation can be used to increase the number alongside the manual lesion delineation, the Ground Truth of training samples and reduce over-fitting (Krizhevsky (GT), of a patient. et al., 2012). Due to the relatively small size of the train- The private dataset considers 23 acute ischaemic stroke ing dataset, we employed artificial data augmentation in patients that underwent clinical therapy, acquired at Bern the supervised portion of our proposal. For each sample, University Hospital in Switzerland. As in ISLES 2017, we applied rotations of 90 , 180 , 270 . each patient is characterized by the same six parametric 3.4. Settings and model training maps, being the final lesion manually delineated at 90-day follow-up T2. The parametric maps were co-registered Unsupervised functional block. The unsupervised func- followed by skull-stripping with FSL BET2 on the co- tional block was trained by optimizing the negative log- registered follow-up T2 image (Jenkinson et al., 2012). likelihood of the data. However, since computing the 8 gradient is generally intractable, we performed the train- jVj p g ing by approximating the gradient with Contrastive Diver- i i Soft Dice loss = : (2) P P jVj jVj gence with one step of alternating Gibbs sampling (Hin- 2 2 p + g i i i i ton, 2012). The training process of an RBM can be dif- In the soft dice loss, the sum occurs over the set V ficult if one tries to learn the parameter  of the energy of voxels belonging to the predicted output patch, where function, which corresponds to the standard deviation of p 2 P denotes the probability of a voxel i in the output the Gaussian noise of a visible node i (Hinton, 2012). Ac- patch and g 2 G corresponds to the respective ground- cording to Hinton (2012), we normalize each component truth label voxel. of the data with zero mean and unit variance, and define The method was implemented using Keras with Ten- = 1. In Table 1, we present the settings used for train- sorflow backend, in a workstation equipped with a GTX ing the unsupervised model. 1080 Ti 11 GB. Prediction time takes around 20s per pa- For training each RBM, we randomly extract 3D tient. patches of shape 7  7  3 from the respective input set of MRI maps, C. Then, the 3D patches are reshaped into a 1D vector and fed into the visible layer of the RBM, hav- 4. Results and Discussion ing an input of size m = 7 7 3jCj, as shown in Figure In this section, we discuss the impact of the main 3. After training, we extract features from the NReLU contributions, namely, the incorporation of unsupervised units noise-free activations. These units exhibit intensity learning with supervised learning and the Gated Recur- equivariance when the bias has zero value, and they are rent blocks. Then, we compare our method with the state noise free units (Nair and Hinton, 2010). Due to the large of the art in ISLES 2017 Challenge. Finally, we delve on number of extracted feature maps (jN j = jN j = 600), 1 2 we perform a feature selection step, as described in Sec- the diculty of predicting the final infarct stroke lesion. tion 2.2, where M = M = 6. The most appropriate 1 2 cardinality of M is discussed in Section 4.1.1. 4.1. Ablation Study The ablation study aims to gradually measure the im- Table 1: Model training parameters for the unsupervised and supervised portance of the main components and consequently assert functional blocks. on the contribution of each component to the overall per- formance. Thus, we start by evaluating the importance of the unsupervised feature generator and the proposed in- Functional Block Parameter Description put grouping. After, we focus on the use of the Gated Optimizer SGD with momentum (lr = 1 10 ) Unsupervised Patch shape 7 7 3 Recurrent block and the choice of the dimensionality of Batch size 32 the spatial context. Optimizer ADAM (lr = 1 10 ) Supervised Patch shape 84 84 4.1.1. Unsupervised feature generation Batch size 4 We hypothesize that grouping the parametric MRI maps according to their physical meaning and encoding each group with an RBM has potential to extract better Supervised functional block. As for the supervised func- features to characterise the stroke lesion and the blood tional block, the complete settings of the training are haemodynamics. We perform several experiments to cor- given in Table 1. For each subject, 350 patches were ran- roborate this working hypothesis. In all experiments, the domly sampled. The training comprehended 36 subjects, while the remaining 7 subjects were used for validation. The settings were optimized through cross-validation in a Additional details of setting and model training are provided in the previous work (Pinto et al., 2018a). For training, we used supplementary material. Also, the source code for reproducing the seg- soft Dice loss function (Milletari et al., 2016). It is defined mentations, the models’ weights and segmentations can be found at: as: https://github.com/apinto92/stroke_prediction.git. 9 parametric MRI maps are also used as input to the super- Grouping parametric MRI maps according to the subja- vised block. The results are presented in Table 2. Figure cent physical meaning. In this experiment, we grouped 4 presents feature maps encoded by the RBMs and the the parametric maps according to their underlining physi- respective MRI maps. cal meaning together with ADC map in each group. Each group was encoded with an RBM. Comparing isolatedly Grouping all parametric MRI maps in a single group. the use of each group of features, we verify that RBM Lesion We considered, first, the e ect of encoding all paramet- had a higher average Dice score compared to using only ric maps using a single RBM. We varied the number of the parametric maps as input to the supervised block. selected features from the RBM, observing that in all The increase in the average Dice score was obtained by cases, the average Dice score is equal or lower than us- a higher average Recall. Also, we observe an improve- ing only the parametric maps as input to the supervised ment in all distance metrics. The experiment of using block. Also, using 12 features presented the lowest av- RBM presented the lowest average Dice and Re- Haemo erage Dice score. The use of 3 or 6 obtained the same call, as well as higher average distance metrics. However, average Dice score, having the second a lower average RBM presented higher average Precision, contrary Haemo Hausdor distance. So, based on the metrics, we may to RBM , which motivated the study on the combina- Lesion conclude that there is no clear gain in using the features tion of features from RBM with RBM besides Haemo Haemo generated by the RBM, at least, when we encode all the the parametric maps. The results of this experiment are parametric maps with a single RBM. presented in Table 2. We may observe that this combi- Since, the selection of 6 features also includes the pre- nation obtained the highest average Dice and Precision, vious top 3 features, we compared the normalized mutual as well as the lowest average distance metrics. How- information between them. As shown in Figure 4, the top ever, this improvement could have been originated from 3 features have low values of normalized mutual informa- the combination of maps according to a specific com- tion in relation to the additional 3 features, which indi- mon property, subjacent physical meaning of the para- cates that there is additional information. For this reason, metric maps, in each group, or because we reduced the we chose 6 as the number of features in the subsequent number of maps from 6 to 3 in each group. And this re- experiments. Table 2: Results obtained with di erent configurations of the unsupervised feature generator block in ISLES 2017 testing set. Each metric represents the mean  standard deviation. Underlined values correspond to the highest mean. Supervised Block Unsupervised Block Dice HD ASSD Precision Recall FCN G-RNN – U-Net LSTM 0.30  0.21 36.58  16.62 6.96  5.08 0.30  0.26 0.55  0.31 RBM (3 Feat.) U-Net LSTM 0.30  0.21 38.93  18.80 6.55  4.22 0.29  0.24 0.61  0.31 Single RBM (6 Feat.) U-Net LSTM 0.30  0.21 36.94  19.19 6.72  4.43 0.29  0.24 0.59  0.31 Single RBM (12 Feat.) U-Net LSTM 0.28  0.20 41.07  18.67 6.81  3.88 0.24  0.21 0.65  0.30 Single RBM U-Net LSTM 0.28  0.24 38.50  22.78 11.09  14.79 0.35  0.30 0.44  0.34 Haemo RBM U-Net LSTM 0.31  0.21 35.38  15.75 6.44  4.43 0.30  0.24 0.59  0.30 Lesion RBM + RBM U-Net LSTM 0.38  0.22 29.21  15.04 5.52  5.06 0.41  0.26 0.53  0.29 Lesion Haemo Two-RBMs U-Net LSTM 0.27  0.21 40.89  14.63 6.92  3.64 0.25  0.23 0.68  0.28 Mixed Three-RBMs U-Net LSTM 0.35  0.23 29.32  14.33 5.27  3.54 0.34  0.27 0.59  0.30 RBM + RBM + RBM Haemo/Less Lesion/Less ADC 10 Figure 4: Onset parametric maps of patient case 0011 in ISLES 2017 training set, alongside the final stroke lesion, at a 90-day follow-up, over the onset ADC map. The subsequent rows show the RBM features selected from the RBM , RBM and RBM , respectively. The last Lesion Haemo Single column shows the normalized mutual information, across whole dataset, among features of the same RBM. duction could have allowed a better training of the RBM. ble 2, the first experiment presented the lowest average So, we performed two complementary experiments. In Dice score and higher average distance metrics, while the the first experiment, we formed two groups with similar second experiment attained the second highest average size, but we randomly chose the parametric maps to in- Dice score, thus showing the importance of splitting the clude in each group. In the second experiment (Three- parametric perfusion maps and including the ADC map in RBMs), we changed the groups of MRI maps encoded both the RBM and RBM . Haemo Lesion in RBM and RBM by removing the ADC map Considering these experiments together, we may draw Lesion Haemo from each one. These two new groups were encoded in some conclusions. First, although CNNs are very e ec- RBM and RBM , respectively. The ADC tive in generating features from raw data, they can gen- Lesion/Less Haemo/Less was separately encoded in RBM . As presented in Ta- erate even better features if rich and complementary in- ADC 11 formation is provided. A similar conclusion was inferred the average distance metrics. by Oliveira et al. (2018) that observed improvement when Based on these experiments, we may conclude that the the coecients of the Wavelet were added as input in the CNN layers were able to extract additional information problem of retinal vessel segmentation. Here, we observe from the RBM features; however, at least to the problem a similar e ect, but using the encoding provided by an of inferring the extension of the lesion months ahead, long RBM trained unsupervisedly for the problem of stroke le- and local distance spatial relations among input voxels in- sion prediction. Second, at least to the problem of stroke troduced by Gated RNN was critical to reduce the detec- lesion prediction, when we have data with di erent latent tion of false positives, increasing substantially the average factors and we are able to group it, according to those Dice score by 6%. factors, then there is potential to extract complementary information from each group, but to mix them all together 4.1.3. Spatial context: 2D or 3D? can be detrimental. MRI images are 3D by nature, so the use of 3D fil- ters would allow capturing more context, which has the 4.1.2. Context aggregation based on gated recurrent potential to provide better prediction. Since 2D filters blocks are confined to a plane, unnatural discontinuous contour In medical imaging segmentation, which is similar to may occur in the perpendicular axis. However, as pre- our problem of inferring the extension of the lesion 90 sented previously, the resolution of MRI images in ISLES days ahead, the use of a cascade of convolutional layers dataset is not equal in all axis, being coarser along the to elaborate the features is the prevalent practise. How- axial axis. So, we studied the e ect of the spatial con- ever, as discussed previously, Gated-RNN layers are able text in our architecture. As we have two blocks, unsuper- to capture long distance spatial relations among input vox- vised and supervised blocks, the e ect on each one was els, so we performed some experiments to evaluate its evaluated separately. The results are presented in Table 4. contribution. The results are presented in Table 3. Considering the results, we observe that using 2D patches Analysing Table 3, we verify that when we just had in both blocks has lower average Dice score, than using parametric maps as input to the supervised block, adding only the parametric maps as input (baseline), because the a LSTM layer increased the average Precision, but the increase in the average Precision was not enough to com- average Recall decreased, resulting in the same average pensate the drop in the average Recall. Using 3D patches Dice score. But, a di erent behaviour is observed when for both blocks had the same performance as our base- we added the features computed from the RBMs. In this line. However, when we used 3D patches for the RBM scenario, we verify that using only CNN layers improved but 2D blocks for the U-Net block, we improved over our over having just parametric maps, which came by a higher baseline. This is the model with the highest average Dice average Precision. However, when we add the LSTM, we score without LSTM. So, we may conclude that for our have an even higher improvement, which is observed in a architecture, larger context using 3D patches was more larger increase in the average Precision, and a decrease in e ective for encoding features in the unsupervised block, Table 3: Results obtained when considering the Gated Recurrent block with and without the unsupervised learning block with ISLES 2017 testing set. Each metric represents the mean  standard deviation. Underlined values correspond to the highest mean. Supervised Block Unsupervised Block Dice HD ASSD Precision Recall FCN G-RNN U-Net – 0.30  0.21 38.83  21.10 7.08  5.15 0.26  0.23 0.64  0.30 U-Net LSTM 0.30  0.21 36.58  16.62 6.96  5.08 0.30  0.26 0.55  0.31 U-Net – 0.32  0.23 34.09  16.51 7.60  7.14 0.35  0.27 0.48  0.32 RBM + RBM [3D] Lesion Haemo U-Net LSTM 0.38  0.22 29.21  15.04 5.52  5.06 0.41  0.26 0.53  0.29 12 Table 4: Evaluation metrics obtained with di erent spatial context configurations in the unsupervised and supervised learning blocks in ISLES 2017 testing set. Each metric represents the mean  standard deviation. Underlined values correspond to the highest mean. Supervised Block Unsupervised Block Dice HD ASSD Precision Recall FCN G-RNN RBM + RBM [2D] U-Net [2D] – 0.27  0.23 36.35  14.89 9.14  12.35 0.31  0.28 0.53  0.34 Lesion Haemo U-Net [2D] – 0.32  0.23 34.09  16.51 7.60  7.14 0.35  0.27 0.48  0.32 RBM + RBM [3D] Lesion Haemo U-Net [3D] – 0.30  0.21 34.17  14.86 6.16  3.82 0.32  0.27 0.54  0.30 while 2D patches were better suited for encoding features which was manually delineated based on a follow-up T2 in the supervised U-Net-based block. MRI acquisitions, are not disclosed for public access. Considering the results, we observe that our baseline is competitive with an average Dice, being among the 4.2. Private dataset top 3 methods, and surpassing the ensemble methods of To further evaluate the generalization capacity of our Pisov et al. (2017) and Robben and Suetens (2017). Our proposal, we tested it on a private dataset and compare proposed method presented the lowest distance metrics it with the baseline method. Table 5 presents the results among all methods, especially for the Hausdor distance. obtained by the two methods. It obtained the second-best average Precision score, be- On the overall, our proposal was capable of surpassing ing surpassed by Robben and Suetens (2017) The authors the baseline model, attaining an higher average Dice, Pre- proposed the integration of meta-data information, using cision and distance metrics, which were statistically sig- a two-pathway 3D network in an ensemble; however, our nificant (Wilcoxon Signed Ranked test with p value < experiments did not indicate any improvement using 3D 0:05). Comparing to ISLES 2017 testing set, there was a patches for the U-Net, at least for our architecture. So, slight decrease in performance. This could be explained this improvement could have come from a combination by the shift on the intensity distribution of the MRI maps, of the e ect of the ensemble and the meta-data. But we due to di erent acquisition protocols or the di erences in note that their method presented a lower average Recall, the preprocessing step. which explains their lower average Dice score. Regard- ing the average Recall score, our method was fourth, but 4.3. State-of-the-art: ISLES 2017 Challenge when we consider the top 3 methods, specially Pinto et al. The results of published methods for final infarct stroke (2018a), we conclude that it was obtained with a much lesion prediction using ISLES 2017 testing set (Winzeck lower average Precision, which means that to increase the et al., 2018), together with our baseline and proposal true positive detections, they had to increase substantially methods are presented in Table 6. The metrics were com- the false positives. So, comparing with the state of the puted by the online platform, so the ground-truth data, art, our method presented a better balance between Pre- Table 5: Results obtained by our proposal and baseline method in the private dataset. Each metric is represented by the mean  standard deviation. Underlined values correspond to the highest mean, while bold values represent statistically significant values (p-value < 0:05). Supervised Block Unsupervised Block Dice HD ASSD Precision Recall FCN G-RNN – U-Net [2D] – 0.32  0.18 32.58  20.09 5.17  3.34 0.31  0.26 0.68  0.28 RBM + RBM [3D] U-Net [2D] LSTM 0.36  0.18 26.68  15.60 3.88  2.17 0.38  0.30 0.68  0.27 Lesion Haemo 13 Table 6: Published methods in ISLES 2017 testing dataset and our proposal. Each metric is represented by the mean  standard deviation. Underlined values correspond to the highest mean. Dice HD ASSD Precision Recall Mok et al. * 0.32  0.23 40.74  27.23 8.97  9.52 0.34  0.27 0.39  0.27 Kwon et al. * 0.31  0.23 45.26  21.04 7.91  7.31 0.36  0.27 0.45  0.30 Robben et al. * 0.27  0.22 37.84  17.75 6.72  4.10 0.44  0.32 0.39  0.31 Pisov et al. * 0.27  0.20 49.24  32.15 9.49  10.56 0.31  0.27 0.39  029 Monteiro et al. * 0.30  0.22 46.60  17.50 6.31  4.05 0.34  0.27 0.51  0.30 Pinto et al. (2018a) 0.29  0.21 41.58  22.04 7.69  5.71 0.21  0.21 0.66  0.29 Lucas et al. * 0.29  0.21 33.85  16.82 6.81  7.18 0.34  0.26 0.51  0.32 Choi et al. * 0.28  0.22 43.89  20.70 8.88 8.19 0.36  0.31 0.41  0.31 Niu et al. * 0.26  0.20 48.88  11.20 6.26  3.02 0.28  0.25 0.56  0.26 Sedlar et al. * 0.20  0.19 58.30  20.02 11.19  9.10 0.23  0.24 0.40  0.29 Rivera et al. * 0.19  0.16 63.58  18.58 11.13  7.89 0.27  0.25 0.21  0.17 Islam et al. * 0.19  0.18 64.15  28.51 14.17  15.80 0.29  0.28 0.25  0.25 Chengwei et al. * 0.18  0.17 65.95  25.94 9.22  6.99 0.37  0.30 0.21  0.23 Yoon et al. * 0.17  0.16 45.23  19.14 12.43  11.01 0.23  0.27 0.36  0.32 Baseline 0.30  0.21 36.58  16.62 6.96  5.08 0.30  0.26 0.55  0.31 Proposal 0.38  0.22 29.21  15.04 5.52  5.06 0.41  0.26 0.53  0.29 Methods presented in Winzeck et al. (2018), whose results were retrieved from SMIR Online Platform (2017). cision and Recall, which reflected into a higher average Dice score. Based on the results, we may conclude that the use of complementary features provided by the RBMs and the use of LSTM for a larger context allowed our baseline to surpass current state-of-the-art methods. Results from ChallengeR Benchmark. The SMIR plat- Figure 5: Boxplot of the top-10 ranking methods ordered by average form of ISLES 2017 provides a weekly benchmark report Dice score in ISLES 2017 testing set. of the current top-10 methods in the testing set, accord- ing to the average Dice score. So, some of the methods may not be published, lacking a description on their im- In Figure 6 we have the significance maps of the pair- plementation, and, for this reason, were not included in wise significant test with one-side Wilcoxon signed rank the previous discussion. test (p-value = 0:05), showing that our method reached Figure 5 presents the boxplots of each method consid- higher Dice score statistically significant when compared ered in the report. with other five ranked methods of the top-10. We observe that the top-10 methods failed to predict the lesion of one or more cases (lowest outliers), which Figure 7 shows the podium plot of each method for each may indicate the degree of complexity of predicting in- case in the testing set, and its ranking. We observe that farct stroke lesion 90 days ahead in ISLES 2017 Chal- our proposal is the method, which ranked first most of the lenge dataset. But we verify that our method is the only times, as well as second and third. Also, when we con- one to have the first quartile above 0.20 in the Dice score. sider the methods ranked bellow fourth, our method is in Single Model Ensemble art, presenting the highest average Dice score and lowest average distance score. Considering the ablation study, this performance was attained due to the combination of adding extra features obtained by encoding the parametric maps with RBMs, according to the underlining physical meaning, and the elaboration provided by the long context of the LSTM layers. 5. Conclusions Figure 6: One-side Wilcoxon signed rank test in ISLES 2017 testing set. Statistically significant tests are marked with the blue colour, while the In this work, we present a deep learning approach for red colour designates statistically non-significant tests. predicting the final stroke lesion, based on unsupervised and supervised learning. We proposed to group the input maps according to the underlying physical principle be- general among those with the lowest counts. Analysing hind their creation, namely, the time-resolved perfusion the cases individually, we note two trends, for some cases maps (Tmax, TTP, MTT), and the blood-flow-dynamic re- all methods presented similar performance, while for oth- lated maps (rCBF, rCBV). Each group was encoded using ers, we find a large variation from the first to the other an unsupervised model to obtain structural features spe- methods. The first trend may be found in the most di- cific to its underlying physical principle. These structural cult case, where all methods had zero or a close value for features together with the standard parametric maps were the Dice score. In the second trend, we observe that our fed to a supervised model to learn features conditioned method is ranked as first most of the cases. on the label, which in our problem, means to condition Based on the results of the benchmark, we may infer on the results of the medical intervention — lesion at 90- that our method is competitive among current state of the days follow-up. We also investigated the use of Gated Figure 7: Podium plot of each testing case in ISLES 2017. For each ISLES 2017 testing subject, defined by a coloured line with circles, the Podium plot orders decreasingly the Dice score obtained by each of the top 10 methods that are represented by coloured circles. 15 Recurrent Neural Networks to provide long spatial con- Sprengers, M.E., Jenniskens, S.F., Lycklama a ` Nije- text, which were critical in relating the structural features holt, G.J., et al., 2016. Collateral status on baseline to the information on input parametric maps. Our results computed tomographic angiography and intra-arterial showed that either the encoding or the long spatial context treatment e ect in patients with proximal anterior cir- improved over our baseline. Also, these two together in- culation stroke. Stroke 47, 768–776. teracted positively increasing the performance when con- Butcher, K., Emery, D., 2010a. Acute stroke imaging part sidering separately each one. i: Fundamentals. Canadian Journal of Neurological When evaluating our proposal on ISLES 2017 testing Sciences 37, 4–16. dataset, we observe a prediction improvement over cur- rent state-of-the-art methods. The proposed method ob- Butcher, K., Emery, D., 2010b. Acute stroke imaging part tained the first place in Dice and also in HD and ASSD. ii: the ischemic penumbra. Canadian Journal of Neu- Recent works (Pinto et al., 2018a; Robben et al., 2020) rological Sciences 37, 17–27. have shown the importance of clinical meta-data to pre- Choi, Y., Kwon, Y., Lee, H., Kim, B.J., Paik, M.C., dict the final stroke lesion in di erent revascularization Won, J.H., 2016. Ensemble of deep convolutional neu- scenarios. So as future work, we aim to study how such ral networks for prognosis of Ischemic Stroke, in: In- meta-data (i.e. TICI score) could be incorporated in our ternational Workshop on Brainlesion: Glioma, Mul- architecture, to consolidate the impact of the clinical in- tiple Sclerosis, Stroke and Traumatic Brain Injuries, tervention and to further improve the 90-day lesion pre- Springer. pp. 231–243. diction. Coutts, S.B., Simon, J.E., Tomanek, A.I., Barber, P.A., Chan, J., Hudon, M.E., Mitchell, J.R., Frayne, R., Acknowledgement Eliasziw, M., Buchan, A.M., et al., 2003. Reliability of assessing percentage of di usion-perfusion mismatch. Adriano Pinto was supported by a scholarship from the Stroke 34, 1681–1683. Fundac ¸ao ˜ para a Ciencia ˆ e Tecnologia (FCT), Portugal (scholarship number PD/BD/113968/2015). This work El Tawil, S., Muir, K.W., 2017. Thrombolysis and was supported by FCT national funds, under the national thrombectomy for acute ischaemic stroke. Clinical support to R&D units grant, through the reference project Medicine 17, 161–165. UIDB/04436/2020 and UIDP/04436/2020. Gonzalez, R., Hirsch, J., Koroshetz, W., Lev, M., Schae- fer, P., 2007. Acute ischemic stroke: imaging and in- tervention. American Journal of Neuroradiology 28, References Badrinarayanan, V., Kendall, A., Cipolla, R., 2015. Gonzalez, ´ R.G., Hirsch, J.A., Koroshetz, W., Lev, M.H., Segnet: A deep convolutional encoder-decoder ar- Schaefer, P.W., 2011. Acute ischemic stroke. Springer. chitecture for image segmentation. arXiv preprint arXiv:1511.00561 . Grysiewicz, R.A., Thomas, K., Pandey, D.K., 2008. Epi- demiology of ischemic and hemorrhagic stroke: inci- Bauer, S., Gratz, P.P., Gralla, J., Reyes, M., Wiest, R., dence, prevalence, mortality, and risk factors. Neuro- 2014. Towards automatic MRI volumetry for treat- logic clinics 26, 871–895. ment selection in acute ischemic stroke patients, in: En- Higashida, R.T., Furlan, A.J., Roberts, H., Tomsick, T., gineering in Medicine and Biology Society (EMBC), Connors, B., Barr, J., Dillon, W., Warach, S., Broder- 2014 36th Annual International Conference of the ick, J., Tilley, B., et al., 2003. Trial design and report- IEEE, IEEE. pp. 1521–1524. ing standards for intraarterial cerebral thrombolysis for Berkhemer, O.A., Jansen, I.G., Beumer, D., Fransen, acute ischemic stroke. Journal of Vascular and Inter- P.S., Van Den Berg, L.A., Yoo, A.J., Lingsma, H.F., ventional Radiology 14, E1–E31. 16 Hinton, G.E., 2012. A practical guide to training restricted Bentley, P., Chen, L., et al., 2017. ISLES 2015-a Boltzmann machines, in: Neural networks: Tricks of public evaluation benchmark for ischemic stroke lesion the trade. Springer, pp. 599–619. segmentation from multispectral MRI. Medical image analysis 35, 250–269. Hochreiter, S., Schmidhuber, J., 1997. Long short-term memory. Neural computation 9, 1735–1780. McKinley, R., Hani, ¨ L., Gralla, J., El-Koussy, M., Bauer, S., Arnold, M., Fischer, U., Jung, S., Mattmann, K., Jenkinson, M., Beckmann, C.F., Behrens, T.E., Woolrich, Reyes, M., et al., 2017. Fully automated stroke tis- M.W., Smith, S.M., 2012. Fsl. Neuroimage 62, 782– sue estimation using random forest classifiers (faster). Journal of Cerebral Blood Flow & Metabolism 37, 2728–2741. Kamnitsas, K., Ledig, C., Newcombe, V.F., Simpson, J.P., Kane, A.D., Menon, D.K., Rueckert, D., Glocker, B., Memezawa, H., Smith, M.L., Siesjo, ¨ B.K., 1992. Penum- 2017. Ecient multi-scale 3d CNN with fully con- bral tissues salvaged by reperfusion following middle nected CRF for accurate brain lesion segmentation. cerebral artery occlusion in rats. Stroke 23, 552–559. Medical image analysis 36, 61–78. Menon, B.K., Qazi, E., Nambiar, V., Foster, L.D., Yeatts, Kemmling, A., Flottmann, F., Forkert, N.D., Minnerup, S.D., Liebeskind, D., Jovin, T.G., Goyal, M., Hill, J., Heindel, W., Thomalla, G., Eckert, B., Knauth, M.D., Tomsick, T.A., et al., 2015. Di erential e ect M., Psychogios, M., Langner, S., Fiehler, J., 2015. of baseline computed tomographic angiography collat- Multivariate dynamic prediction of ischemic infarction erals on clinical outcome in patients enrolled in the in- and tissue salvage as a function of time and degree terventional management of stroke iii trial. Stroke 46, of recanalization. Journal of Cerebral Blood Flow & 1239–1244. Metabolism 35, 1397–1405. Milletari, F., Navab, N., Ahmadi, S.A., 2016. V-net: Fully Kistler, M., Bonaretti, S., Pfahrer, M., Niklaus, R., convolutional neural networks for volumetric medical Buchler ¨ , P., 2013. The virtual skeleton database: an image segmentation, in: 3D Vision (3DV), 2016 Fourth open access repository for biomedical research and col- International Conference on, IEEE. pp. 565–571. laboration. Journal of medical Internet research 15. Mok, T.C., Chung, A.C., 2017. Deep adversarial net- Krizhevsky, A., Sutskever, I., Hinton, G.E., 2012. Ima- works for stroke lesion segmentation . genet classification with deep convolutional neural net- Monteiro, M., Oliveira, A.L., 2017. Fully convolutional works, in: Advances in neural information processing neural network for 3d stroke lesion segmentation . systems, pp. 1097–1105. Nair, V., Hinton, G.E., 2010. Rectified linear units im- Labeyrie, M.A., Turc, G., Hess, A., Hervo, P., Mas, J.L., prove restricted Boltzmann machines, in: Proceedings Meder, J.F., Baron, J.C., Touze, E., Oppenheim, C., of the 27th international conference on machine learn- 2012. Di usion lesion reversal after thrombolysis: a ing (ICML-10), pp. 807–814. mr correlate of early neurological improvement. Stroke 43, 2986–2991. Nielsen, A., Hansen, M.B., Tietze, A., Mouridsen, K., 2018. Prediction of tissue outcome and assessment Liebeskind, D.S., 2003. Collateral circulation. Stroke 34, of treatment e ect in acute ischemic stroke using deep 2279–2284. learning. Stroke 49, 1394–1401. Lucas, C., Heinrich, M.P., 2017. 2d multi-scale res-net Niu, Y., Gong, E., Xu, J., Pauly, J., Zaharchuk, G., for stroke segmentation. 2018. Improved prediction of the final infarct from Maier, O., Menze, B.H., von der Gablentz, J., Hani, ¨ L., acute stroke neuroimaging using deep learning, in: Heinrich, M.P., Liebrand, M., Winzeck, S., Basit, A., STROKE. 17 Oliveira, A., Pereira, S., Silva, C.A., 2018. Retinal vessel Robben, D., Boers, A.M., Marquering, H.A., Langezaal, segmentation based on fully convolutional neural net- L.L., Roos, Y.B., van Oostenbrugge, R.J., van Zwam, works. Expert Systems with Applications 112, 229– W.H., Dippel, D.W., Majoie, C.B., van der Lugt, A., 242. et al., 2020. Prediction of final infarct volume from native ct perfusion and treatment parameters using deep Pereira, S., Meier, R., McKinley, R., Wiest, R., Alves, learning. Medical image analysis 59, 101589. V., Silva, C.A., Reyes, M., 2018. Enhancing inter- Robben, D., Suetens, P., 2017. Dual-scale fully convo- pretability of automatically extracted machine learning lutional neural network for final infarct prediction, in: features: application to a RBM-Random Forest system Ischemic stroke lesion segmentation-ISLES challenge on brain lesion segmentation. Medical image analysis 2017, held in conjunction with MICCAI 2017, Date: 44, 228–244. 2017/09/10-2017/09/10, Location: Quebec City, Que- Pereira, S., Pinto, A., Amorim, J., Ribeiro, A., Alves, bec, Canada. V., Silva, C.A., 2019. Adaptive feature recombination Ronneberger, O., Fischer, P., Brox, T., 2015. U-net: and recalibration for semantic segmentation with fully Convolutional networks for biomedical image seg- convolutional networks. IEEE Transactions on Medi- mentation, in: International Conference on Medical cal Imaging 38, 2914–2925. doi:10.1109/TMI.2019. image computing and computer-assisted intervention, Springer. pp. 234–241. Pinto, A., McKinley, R., Alves, V., Wiest, R., Silva, C.A., Rose, S.E., Chalk, J.B., Grin, M.P., Janke, A.L., Chen, Reyes, M., et al., 2018a. Stroke lesion outcome pre- F., McLachan, G.J., Peel, D., Zelaya, F.O., Markus, diction based on MRI imaging combined with clinical H.S., Jones, D.K., et al., 2001. MRI based di usion information. Frontiers in Neurology 9, 1060. and perfusion predictive model to estimate stroke evo- lution. Magnetic resonance imaging 19, 1043–1053. Pinto, A., Pereira, S., Meier, R., Alves, V., Wiest, R., Silva, C.A., Reyes, M., 2018b. Enhancing clinical MRI Rumelhart, D.E., McClelland, J.L., 1986. Parallel dis- perfusion maps with data-driven maps of complemen- tributed processing: explorations in the microstructure tary nature for lesion outcome prediction, in: Medical of cognition. volume 1. foundations . Image Computing and Computer Assisted Intervention Scalzo, F., Hao, Q., Alger, J.R., Hu, X., Liebeskind, D.S., – MICCAI 2018, pp. 107–115. 2012. Regional prediction of tissue fate in acute is- chemic stroke. Annals of Biomedical Engineering 40, Pisov, M., Belyaev, M., Krivov, E., 2017. Neural networks 2177–2187. ensembles for ischemic stroke lesion segmentation . SMIR Online Platform, 2017. Ischemic stroke lesion Powers, W.J., Rabinstein, A.A., Ackerson, T., Adeoye, segmentation 2017. https://www.smir.ch/ISLES/ O.M., Bambakidis, N.C., Becker, K., Biller, J., Brown, Start2017. [Acessed: 2019-11-18]. M., Demaerschalk, B.M., Hoh, B., et al., 2018. 2018 guidelines for the early management of patients Vinh, N.X., Epps, J., Bailey, J., 2010. Information the- with acute ischemic stroke: a guideline for health- oretic measures for clusterings comparison: Variants, care professionals from the American Heart Associa- properties, normalization and correction for chance. tion/American Stroke Association. Stroke 49, e46–e99. Journal of Machine Learning Research 11, 2837–2854. Rekik, I., Allassonniere, ` S., Carpenter, T.K., Ward- Visin, F., Ciccone, M., Romero, A., Kastner, K., Cho, K., law, J.M., 2012. Medical image analysis methods Bengio, Y., Matteucci, M., Courville, A., 2016. Re- in MR/CT-imaged acute-subacute ischemic stroke le- seg: A recurrent neural network-based model for se- sion: Segmentation, prediction and insights into dy- mantic segmentation, in: Proceedings of the IEEE Con- namic evolution simulation models. a critical appraisal. ference on Computer Vision and Pattern Recognition NeuroImage: Clinical 1, 164–178. Workshops, pp. 41–48. 18 Winzeck, S., Hakim, A., McKinley, R., Pinto, J.A.A.D.S., Alves, V., Silva, C., Pisov, M., Krivov, E., Belyaev, M., Monteiro, M., et al., 2018. ISLES 2016 & 2017- benchmarking ischemic stroke lesion outcome predic- tion based on multispectral MRI. Frontiers in Neurol- ogy 9, 679. World Health Organization, et al., 2014. Global status re- port on noncommunicable diseases 2014. World Health Organization. Zeiler, M.D., Fergus, R., 2014. Visualizing and under- standing convolutional networks, in: European confer- ence on computer vision, Springer. pp. 818–833. Zivelonghi, C., Tamburin, S., 2018. Mechanical thrombectomy for acute ischemic stroke: the therapeu- tic window is larger but still “time is brain”. Functional neurology 33, 5.

Journal

Computing Research Repository – arXiv (Cornell University)

Published: Jan 2, 2021

There are no references for this article.