We included patients registered between March 2014 and June 2016 in the Multicenter Randomized Clinical Trial of Endovascular Treatment for Acute Ischemic Stroke in the Netherlands (MR CLEAN) Registry. The MR CLEAN Registry is an ongoing, national, prospective, open, multicenter, observational monitoring study covering all 18 stroke intervention centers that perform EVT in the Netherlands, of which 16 participated in the MR CLEAN trial (10). The registry is a continuation of the MR CLEAN trial collaboration and includes all patients undergoing EVT (defined as entry into the angiography suite and receiving arterial puncture) for acute ischemic stroke in the anterior and posterior circulation. In the current analysis we included those patients who adhered to the following criteria: age 18 years and older, treatment in a center that participated in the MR CLEAN trial, and LVO in the anterior circulation (internal carotid artery (ICA), internal carotid artery terminus (ICA-T), middle (M1/M2) cerebral artery, or anterior (A1/A2) cerebral artery), shown by CT angiography (CTA) or digital subtraction angiography (DSA) (11).

#### Clinical Baseline Characteristics

We assessed the following clinical characteristics at admission: National Institutes of Health Stroke Scale (NIHSS), Glasgow Coma Scale, medical history (TIA, ischemic stroke, intracranial hemorrhage, subarachnoid hemorrhage, myocardial infarction, peripheral artery disease, diabetes mellitus, hypertension, hypercholesterolemia), smoking, laboratory tests (blood glucose, INR, creatinine, thrombocyte count, CRP), blood pressure, medication (thrombocyte aggregation inhibitors, oral anticoagulant drugs, anti-hypertensive drugs, statins), modified Rankin Score (mRS) before stroke onset, administration of intravenous tPA (yes/no), stroke onset to groin time, transfer from another hospital, and whether the patient was admitted during weekend or off hours.

### Radiological Baseline Parameters

All imaging in the MR CLEAN Registry was assessed by an imaging core laboratory (11). On non-contrast CT, the size of initial lesion in the anterior circulation was assessed by the Alberta Stroke Program Early CT Score (ASPECTS). ASPECTS is a 10 point quantitative topographic score representing early ischemic change in the middle cerebral artery territory, with a scan without ischemic changes receiving an ASPECTS of 10 points (12). In addition, presence of leukoaraiosis and old infarctions, hyperdense vessel sign, and hemorrhagic transformation of the ischemic lesion were assessed on noncontrast CT.

On CTA, the core lab determined clot burden score, clot location, collaterals, and presence of intracranial atherosclerosis. The clot burden score evaluates the extent of thrombus in the anterior circulation by location scored on a 0–10 scale. A score of 10 is normal, implying clot absence; a score of 0 implies complete multi-segment vessel occlusion (12). Presence of intracranial carotid artery stenosis, atherosclerotic occlusion, floating thrombus, pseudo-occlusion, and carotid dissection were scored on CTA of the carotid arteries. Collaterals were assessed using a 4 point scale, with 0 for absent collaterals (0% filling of the vascular territory downstream of the occlusion), 1 for poor collaterals (>0% and ≤50% filling of the vascular territory downstream of the occlusion), 2 for moderate collaterals (>50% and <100% filling of the vascular territory downstream of the occlusion), and 3 for excellent collaterals (100% filling of the vascular territory downstream of the occlusion) (13).

#### Treatment Specific Variables

Variables collected during EVT were type of sedation during the procedure (general or conscious), use of a balloon guiding catheter, carotid stent placement, performed procedure (DSA only or thrombectomy), and type of EVT-device (stent retriever, aspiration device, or a combination of both). In addition, data were collected on adverse events during the procedure (perforation, dissection, distal thrombosis on DSA).

Interventional DSA parameters in our dataset were occluded vessel segment (ICA: origin, cervical, petrous, cavernous, supraclinoid, M1-M4, A1, A2), arterial occlusive lesion (AOL) recanalization score before and after EVT (14), evidence of vascular injury (i.e., perforation, or dissection, vasospasm, new clot in different vascular territory or distal thrombus confirmed with imaging), and modified Thrombolysis in Cerebral Infarction (mTICI)-score before and after EVT. The mTICI-score grades the following categories of cerebral reperfusion: no reperfusion of the distal vascular territory (0), minimal flow past the occlusion but no reperfusion (1), minor partial reperfusion (2a), major partial reperfusion (2b), and complete reperfusion (3) (15). Further variables analyzed were time from stroke onset to recanalization, post-EVT stay on intensive care, high care or stroke care, NIHSS after EVT (<48 h), delta NIHSS (pre-treatment NIHSS subtracted from NIHSS <48 h after EVT) and hemicraniectomy or symptomatic intracranial hemorrhage <48 h after EVT.

### Outcome

The primary radiological outcome was good reperfusion defined as modified TICI-score directly post-procedure (post-mTICI) ≥ 2b (15). The primary clinical outcome was functional independence at 3 months after stroke (mRS ≤ 2). We excluded patients in whom any of the main outcomes (3-months mRS and post-mTICI) were missing.

To investigate the full potential of Machine learning compared with conventional methods in different settings after stroke we defined two prediction settings:

First, we assessed the probability of good reperfusion and good 3-months functional independence in our cohort of patients that underwent EVT based only on variables that were available on admission before entry into the angiography suite. With this baseline prediction setting we are able to investigate the added value of machine learning for models that could potentially support future clinical decision making regarding the performance of EVT yes or no.

Second, we tested the models for predicting 3-months functional independence in patients after EVT was performed. For this analysis we used all variables collected up to 48 h after the end of the endovascular procedure (baseline and treatment variables).

## Machine Learning Algorithms

The machine learning algorithms used in our study were Random Forests, Artificial Neural Network and Support Vector Machine, because they are among the algorithms that are currently most widely and successfully used for clinical data (2–7). Each one of them represents a different algorithm "family," each with radically different internal algorithm structures (16). Since it was not known beforehand which kind of algorithm would perform best, we chose algorithms with different internal structures to increase the probability of good discriminative performance. We also used Super Learner, which is an ensemble method that finds the optimal weighted combination of predictions of the Random Forests, Artificial Neural Network and Support Vector Machine algorithms used in this study. Ensemble methods, such as Super Learner have been shown to increase predictive performance by increasing model flexibility (17). For the implementation of all machine learning algorithms we used off-the-shelf methods in the Python module Scikit-Learn (18).

#### Super Learner

Super Learner is a stacking algorithm using cross-validated predictions of other models (i.e., a machine learning algorithm and logistic regression) and assigning weights to these predictions to optimize the final prediction. Super Learner's predictive performance has been found to surpass individual machine learning models in various clinical studies (17, 19, 20).

#### Random Forests

Random Forests consists of a collection of decision tree classifiers that are fit on random subsamples of patients and variables in the dataset. The variation of the subsampled variables creates a robust classifier. In the decision trees, each node represents a variable and splits the input data into branches based on an objective function that determines the optimal threshold for separating the outcome classes. The predictions from each tree are used as "votes," and the outcome with the most votes is considered the predicted outcome for that specific patient (6, 21). From the Random Forests algorithm variable importances can be derived, which are the sum of weights of nodes of the trees containing a certain variable, averaged over all trees in the forest (22).

#### Support Vector Machine

Support Vector Machine (SVM) is a kernel-based supervised machine learning classifier which can also be used to output probabilities. The SVM works by first mapping the input data into a high dimensional variable space. For binary classification, a hyperplane is subsequently determined to separate two classes such that the distance between the hyperplane and the closest data points is maximized (23).

#### Artificial Neural Network

In this study we use the multilayer perceptron, a popular class of artificial neural network architecture composed of one or more interconnected layers of neurons that process data from the input layer into predictions for the output layer. The algorithm computes a weight for each neuron based on input activation. These weights are updated by backpropagation and stochastic gradient descent (24, 25).

### Logistic Regression

For logistic regression, generally a set of variables has to be selected to be included in the model. Since model performance can rely heavily on selecting the right variables, we tested five different variable selection methods prior to logistic regression. We first selected variables based on prior knowledge, a still widely used method that could be considered "classical" (26). We selected 13 variables available at baseline that were included in a previous study for a similar purpose (27) (**Supplementary Table Ia**). In addition, from baseline and treatment variables we selected 15 variables based on expert opinions of vascular neurologists and radiologists (**Supplementary Table Ib**).

We further considered four automated variable selection methods: (i) backward elimination, which is also considered to be a more classical approach (26), and three state-of-the-art variable selection methods: (ii) least absolute shrinkage and selection operator (LASSO) (28), (iii) Elastic Net, which is a modification of the LASSO found to outperform the former while still having the advantage of a similar sparsity of representation (29), and (iv) selection based on Random Forests variable importance.

#### Analysis Pipeline

We imputed missing values using multiple imputations by chained equations (MICE) (30). Variables with 25% missing values or more were discarded from further analysis. All remaining variables used in this study are listed in **Supplementary Tables II**, **III**. In total, 53 baseline variables and 30 treatment variables were used as input for machine learning algorithms and automated variable selection methods for logistic regression.

The ordinal clinical (NIHSS) and radiological (clot burden and ASPECTS) scores were presented as continuous scores in all models to increase model efficiency, and we assumed linear trends underlying the ordinal scores.

We used nested cross-validation (CV), consisting of an outer and an inner CV loop. In the outer CV loop we used stratified CV with 100 repeated random splits resulting in a training set including 80% and a test set including 20% of all patients. Each training set was used as input for the inner CV loop, consisting of 10-fold CV (31, 32). In the inner CV loop we selected variables for the logistic regression models using the different variable selection methods, and optimized hyperparameters of all machine learning models. Hyperparameters are tuning parameters specific to each machine learning algorithm whose values have to be preset and cannot be directly learned from the data. We optimized hyperparameters with the random grid search module from Scikit-Learn (18). We selected those with highest area under the receiver operating characteristic (AUC) across all internal CV folds to find the best set of selected variables and hyperparameters. **Figure 1** shows a schematic representation of our nested CV methodology.

For all Random Forests models of both prediction settings we ranked variables by decreasing variable importance. For each variable we assessed the frequency of being among the 15 most

\*National Institutes of Health Stroke Scale score.

\*\*Direct oral anticoagulant drugs.

important variables in a Random Forests model for each of the 100 external CV folds (**Table 3**).

#### Model Performance

We assessed model discrimination (the ability to differentiate between patients with good and poor outcome) with receiver operating characteristic (ROC) analyses. Because of our outer CV loop with 100 repeated random splits, we obtained 100 different AUCs from every model. We computed the average ROC-curve and mean AUC with 95% confidence intervals (CI) for all models. We evaluated differences between mean AUCs of the best performing machine learning model and best performing logistic regression model by computing the difference of means including the associated 95% CI.

#### RESULTS

Of the 1,627 patients registered between March 2014 and June 2016, we excluded 244 patients for this analysis because of age <18 (n = 2), posterior circulation stroke (n = 79), missing MR CLEAN trial center (n = 20), and missing mRS or postmTICI (n = 143). Mean age was 69.8 years (SD ± 14.4) and 738 (54%) of the 1,383 included patients were men. In total, 531 (38%) patients had good reperfusion after EVT and 525 (38%) were functionally independent (mRS ≤ 2) 3 months after stroke. Baseline characteristics are shown in **Table 1**.

### Prediction of Good Reperfusion After EVT in Patients at Time of Admission

Discrimination between good and poor reperfusion of the best machine learning algorithm (Support Vector Machine, mean AUC: 0.55) and the best logistic regression model (using backward elimination, mean AUC: 0.57) was similar (difference of mean AUCs: 0.02; 95% CI: 0.01–0.03). Discrimination was poor for all models, with a mean AUCs ranging from 0.53 to 0.57 (**Table 2**). Variable selection using LASSO or Elastic Net was not possible likely because the signal-to-noise ratio was insufficient (18).

### Prediction of 3-Months Functional Independence in Patients at Time of Admission

Discrimination of good functional outcome of the best machine learning algorithm (Super Learner, mean AUC: 0.79) and the best logistic regression model (using LASSO, mean AUC: 0.78) was similar (difference of mean AUCs: 0.01; 95% CI: 0.01–0.01).

Discrimination was moderate for all models, with a mean AUCs ranging from 0.77 to 0.79.

### Prediction of 3-Months Functional Independence in Patients After Performance of EVT

Discrimination of good functional outcome of the best machine learning algorithm (Random Forests, mean AUC: 0.91) and the best logistic regression model (using prior knowledge, mean AUC: 0.90) was similar (difference of mean AUCs: 0.01; 95% CI: 0.00–0.01).

Discrimination was good for all models, with mean AUCs ranging from 0.88 to 0.91.

We performed a post-hoc analysis in patients with good reperfusion as defined by post-mTICI ≥ 2b, predicting 3 months functional outcome both at time of admission and after performance of EVT. We did not find significant differences in performance between machine learning algorithms and logistic regression models in this patient subset (data not shown).

In **Table 3** we show the top 15 variables based on the frequency of being among the 15 most important variables in a Random Forests model for each of the 100 external CV folds.

### DISCUSSION

We found no difference in performance between best performing machine learning algorithms and best performing logistic regression models in predicting radiological or clinical outcome in stroke patients treated with EVT. For prediction of good reperfusion using variables available at baseline, all models showed a poor discriminative performance. This could indicate that reperfusion after EVT depends on characteristics not present in our variables available at time of admission, such as vascular anatomy or interventionalist related factors. Prediction TABLE 2 | Discrimination of machine learning algorithms and logistic regression models across the various prediction settings.

\*Model discrimination is assessed by calculating mean Area Under the Curve (AUC) of the receiver operating characteristic across all outer cross-validation folds.

\*\*Logistic regression using automated variable selection methods.

UVariable selection not possible, likely due to insufficient signal-to-noise ratio.

‡Logistic regression using variables based on prior knowledge.

TABLE 3 | Variable importance of Random Forests for various prediction settings (used variables: predicted outcome).

NCCT, non-contrast CT; CRP, C-Reactive Protein; RR, blood pressure; NIHSS, National Institutes of Health Stroke Scale score.

\*Frequency of being among the 15 most important variables in a Random Forests model for each of the 100 external CV folds.

\*\*Location of intracranial occlusion on CTA.

of 3-months functional independence using variables known at baseline was moderate, predicting 3-months functional independence using baseline and treatment variables resulted in good performance.

We hypothesized that machine learning would outperform logistic regression models due to simultaneous assessment of a large number of variables, and more efficient processing of nonlinear relations and interactions between them. Although a large number of variables (83 in total, see **Supplementary Tables II**, **III**) were available for analysis in the MR CLEAN Registry database, performance of best machine learning algorithms and best logistic regression models were similar. This could indicate that interactions and non-linear relationships in our dataset were of limited importance.

To interpret our results, several methodological limitations have to be considered. First, due to their great flexibility machine learning algorithms are prone to overfitting, which results in optimistic prediction performance. To account for overfitting we used nested CV, which is considered to be an effective method for this aim (33). Second, our outer CV loop resulted in 100 AUCs per model leading to relatively small confidence intervals of mean AUCs. Although this increases the probability of statistically significant differences between mean AUCs of various models, the clinical relevance of these mean AUC differences is difficult to interpret. Because in our study mean AUC differences between models are minimal, clinical relevance of these differences is also negligible. Third, we used data from a registry. Registries might be prone to selection bias. However, we expect that selection bias in our study was minimal because the MR CLEAN Registry in principle covers all patients treated with EVT in the Netherlands. In addition, in all centers patients were treated according to national guidelines, and registration of treatment was a prerequisite for reimbursement (11).

Strong points of this study include the large sample size and standardized collection of patient data. Moreover, because of extensive hyperparameter tuning and state-of-the art variable selection methods, machine learning and logistic regression models were compared at their best performance. In several other studies that compared machine learning algorithms with only logistic regression methods using variables based on prior knowledge, machine learning outperformed logistic regression (6, 7, 34). Variable selection based on prior knowledge has the major drawback that predictive patterns in the data may be missed, as variable selection is strictly based on the literature and expert opinion (26). In our study however, logistic regression using variables based on prior knowledge performed similarly to logistic regression using automated variable selection methods.

The distinction between machine learning and "classical" regression methods is largely artificial. However, a clear distinction between various machine learning algorithms and logistic regression exists in terms of model transparency, which could be seen as the understanding of the mechanism by which the model works (35). Logistic regression has the advantage of transparency at the level of individual variable coefficients, since from these coefficients odds ratios can be derived. However, variable importances derived from the Random Forests algorithm also offer insight in the importance of individual variables for prediction performance (22). These variable importances take interaction between variables into account and have a similar interpretation for continuous and discrete variables, unlike odds ratios which constitute an effect per unit change of a predictor. Hence, Random Forests could be used as an efficient screening tool to pick up predictive patterns in the data that could potentially lead to further hypothesis-driven research. In **Table 3** we show the top 15 variables from either the baseline or baseline and treatment variable set, based on Random Forests variable importance. The majority of variables in **Table 3** do not overlap with the selection of variables based on prior knowledge, potentially providing researcher with additional information.

In this dataset we found no clinically relevant differences in prediction of reperfusion and 3-months functional independence across all models. However, since it is generally not known on beforehand which type of model will result in the best predictive performance in a new dataset, our methodology could be of importance in future studies. We present an analysis pipeline with both machine learning algorithms and logistic regression models including state-of-the-art variable selection methods. Assessing predictive performance of all models simultaneously enables the researcher to make the proper trade-off between predictive performance and model transparency. As our analysis pipeline is fully automated and input variables and outcome label can be altered at will, it is relatively easy to reuse in future studies. The Python code of our pipeline has been made publicly available

#### REFERENCES

1. Goyal M, Menon BK, van Zwam WH, Dippel DW, Mitchell PJ, Demchuk AM, et al. Endovascular thrombectomy after large-vessel ischaemic stroke: a meta-analysis of individual patient data from five randomised trials. Lancet (2016) 387:1723–31. doi: 10.1016/S0140-6736(16) 00163-X

in an online repository (https://github.com/L-Ramos/MrClean\_ Machine\_Learning).

### ETHICS STATEMENT

The central medical ethics committee of the Erasmus Medical Centre Rotterdam, the Netherlands, evaluated the MR CLEAN Registry protocol and granted permission to carry out the study as a registry. All subjects gave written informed consent in accordance with the Declaration of Helsinki.

### AUTHOR CONTRIBUTIONS

HvO lead author, study design, analysis and interpretation, critical revision manuscript for important intellectual content. LR study design, analysis and interpretation, critical revision of manuscript for important intellectual content. AH, MvL, ES, HL, SO, KZ, EV, and HM study design, critical revision of manuscript for important intellectual content. NK data acquisition, critical revision of manuscript for important intellectual content. DD, MvW, IvdS, WS, and CM data acquisition, critical revision of manuscript for important intellectual content. MW supervisor of lead author, data acquisition, study design, critical revision of manuscript for important intellectual content.

### FUNDING

MW was supported by a personal ZonMw VIDI grant (91717337), a Dekker Clinical Established Investigator Grant from the Netherlands Heart Foundation (2016T068) and a Fellowship grant from the Netherlands Brain Foundation (F2014(1)-22). LR was supported by an ITEA 3 grant (14003 Medolution). The MR CLEAN Registry is partially funded by unrestricted grants from Toegepast Wetenschappelijk Instituut voor Neuromodulatie, Twente University (Twin), Erasmus MC, AMC, and MUMC.

#### ACKNOWLEDGMENTS

See supplementary file for full details of the acknowledgements (MR CLEAN Registry acknowledgements).

#### SUPPLEMENTARY MATERIAL

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

hepatitis c by incorporating longitudinal data. Hepatology (2015) 61:1832–41. doi: 10.1002/hep.27750

**Conflict of Interest Statement:** DD reports grants from the Dutch Heart Foundation, AngioCare, Medtronic/Covidien/EV3, MEDAC/LAMEPRO, Penumbra, Top Medical/Concentric, and Stryker during conduct of the study; grants from Stryker European Operations BV, Medtronic, Dutch Heart Foundation, Brain Foundation Netherlands, The Netherlands Organisation for Health Research and Development, Health Holland Top Sector Life Sciences and Health, and consultation fees from Stryker, Bracco Imaging, and Servier, received by the Erasmus University Medical Centre, outside the submitted work. CM reports grants from TWIN, during the conduct of the study and grants from CVON/Dutch Heart Foundation, European Commission and from Stryker outside the submitted work (paid to institution), and is shareholder of Nico.lab. HM is founder and shareholder of Nico-lab.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 van Os, Ramos, Hilbert, van Leeuwen, van Walderveen, Kruyt, Dippel, Steyerberg, van der Schaaf, Lingsma, Schonewille, Majoie, Olabarriaga, Zwinderman, Venema, Marquering, Wermer and the MR CLEAN Registry Investigators. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Decision Criteria for Large Vessel Occlusion Using Transcranial Doppler Waveform Morphology

Samuel G. Thorpe1 \*, Corey M. Thibeault 1 , Nicolas Canac1 , Seth J. Wilk 1 , Thomas Devlin2 and Robert B. Hamilton1

*1 Neural Analytics, Inc., Los Angeles, CA, United States, 2 Erlanger Medical Center, Chattanooga, TN, United States*

#### Edited by:

*David S. Liebeskind, University of California, Los Angeles, United States*

#### Reviewed by:

*Qing Hao, Icahn School of Medicine at Mount Sinai, United States Marek Czosnyka, University of Cambridge, United Kingdom*

\*Correspondence: *Samuel G. Thorpe<EMAIL_ADDRESS>
#### Specialty section:

*This article was submitted to Stroke, a section of the journal Frontiers in Neurology*

Received: *04 May 2018* Accepted: *21 September 2018* Published: *17 October 2018*

#### Citation:

*Thorpe SG, Thibeault CM, Canac N, Wilk SJ, Devlin T and Hamilton RB (2018) Decision Criteria for Large Vessel Occlusion Using Transcranial Doppler Waveform Morphology. Front. Neurol. 9:847. doi: 10.3389/fneur.2018.00847* Background: The current lack of effective tools for prehospital identification of Large Vessel Occlusion (LVO) represents a significant barrier to efficient triage of stroke patients and detriment to treatment efficacy. The validation of objective Transcranial Doppler (TCD) metrics for LVO detection could provide first responders with requisite tools for informing stroke transfer decisions, dramatically improving patient care.

Objective: To compare the diagnostic efficacy of two such candidate metrics: Velocity Asymmetry Index (VAI), which quantifies disparity of blood flow velocity across the cerebral hemispheres, and Velocity Curvature Index (VCI), a recently proposed TCD morphological biomarker. Additionally, we investigate a simple decision tree combining both metrics.

Methods: We retrospectively compare accuracy/sensitivity/specificity (ACC/SEN/SPE) of each method (relative to standard CT-Angiography) in detecting LVO in a population of 66 subjects presenting with stroke symptoms (33 with CTA-confirmed LVO), enrolled consecutively at Erlanger Southeast Regional Stroke Center in Chattanooga, TN.

Results: Individual VCI and VAI metrics demonstrated robust performance, with area under receiver operating characteristic curve (ROC-AUC) of 94% and 88%, respectively. Additionally, leave-one-out cross-validation at optimal identified thresholds resulted in 88% ACC (88% SEN) for VCI, vs. 79% ACC (76% SEN) for VAI. When combined, the resultant decision tree achieved 91% ACC (94% SEN).

Discussion: We conclude VCI to be superior to VAI for LVO detection, and provide evidence that simple decision criteria incorporating both metrics may further optimize.

Performance: Our results suggest that machine-learning approaches to TCD morphological analysis may soon enable robust prehospital LVO identification.

Registration: Was not required for this feasibility study.

Keywords: transcranial doppler, ultrasound, ischemic stroke, large vessel occlusion, decision tree, diagnostic biomarker

## INTRODUCTION

Acute Ischemic Stroke (AIS) is the leading cause of longterm disability in the United States, and fifth leading cause of death (1). Current treatment for AIS includes the use of intravenous tissue Plasminogen Activator, and endovascular mechanical thrombectomy with a clot extraction or aspiration device. Although these therapies provide effective treatment options for Large Vessel Occlusion (LVO), their use is still limited by short time windows from symptom onset during which they are optimally effective (2–4). Indeed, only a small fraction of candidate patients who could ultimately benefit from endovascular treatment currently receive it (5). Early LVO identification is key to enabling rapid triage and transfer to comprehensive stroke centers, thus facilitating access to appropriate care. Computed Tomography Angiography (CTA) is the current gold standard for stroke diagnosis, but is limited to inhospital use, or a low number of prohibitively expensive mobile stroke ambulances. Unfortunately, current prehospital stroke assessment scales lack reliability due to training requirements and low inherent accuracies (6, 7); causing delays in triage, transfer, and treatment.

Transcranial Doppler (TCD) ultrasound is a reliable diagnostic tool for assessing the presence and severity of LVO (8–11), which has the additional advantages of being noninvasive, inexpensive, and portable. Because it directly measures Cerebral Blood Flow Velocity (CBFV), TCD is a strong candidate technology for prehospital diagnosis and assessment of LVO. Indeed, bedside TCD examinations to detect stenosed and/or occluded intracranial vessels are routinely conducted as standard of care at many comprehensive stroke centers (12). Numerous studies have been published comparing TCD diagnosis of arterial LVO with CTA imaging; reporting sensitivity (SEN) and specificity (SPE) ranging between 79 and 98% depending on occlusion location (13–17). A limiting factor of these studies is the TCD operator's ability to locate and interpret the CBFV waveform. Such challenges have contributed to TCD being critically underutilized for stroke assessment.

For stroke diagnosis, specialized training is required to inspect flow velocity and morphology across multiple vessels and depths in both cerebral hemispheres. One of the most cited papers for stroke diagnosis using TCD was published by Demchuk et al. (10), which instructs the operator to categorize waveforms according to evidence of stroke-related pathology; namely dampened, blunted, minimal, or absent signal. A number of additional TCD exam methodologies with different criterion for LVO assessment have been published (15, 17). Typically, CBFV and power M-mode (PMD) waveforms are obtained for flow through the Middle, Anterior, and Posterior Cerebral Arteries (MCA, ACA, and PCA) in each cerebral hemisphere, as well as the Internal Carotid Arteries (ICA). Heuristic assessments are then made based on numerous features, including relative velocities, collateral flow, PMD resistance signatures, and the presence of pathological waveform morphologies.

Assessment of these categories relies heavily on qualitative interpretation by specialists which cannot be replicated by less formally trained personnel. The challenge of moving LVO detection to the prehospital setting thereby obviates the need for objective metrics by which first responders might reliably evaluate TCD signals. An intuitive first candidate for such a metric is CBFV asymmetry, as it is already well established that velocity disparity, both between homologous vessels in opposite hemispheres (10, 18) as well as adjacent vessels in an occluded hemisphere (10, 15), can be indicative of vascular occlusion. One promising metric for LVO detection based on velocity disparity was published by (15); showing area under the Receiver Operating Characteristic curve (ROC-AUC) of 92.6%. However, their metric also relied on PMD resistance signatures as a predictive feature, which were not objectively computed, and was limited in application to occlusions of the MCA.

However intuitive, assessment of velocity asymmetry also comes with the inherent concern that velocity estimates in adjacent vessels and opposite hemispheres can be greatly impacted by anatomical variability (incident angle of the vessel and ultrasound beam), as well as by intrarater measurement inconsistency (19). Moreover, reliance on mean velocity disparity inherently discards the morphological information currently utilized in routine stroke assessment protocols. Such assessments incorporate morphological information explicitly, but in a subjective manner which requires expert interpretation. However, a number of recent studies have observed morphological changes associated with various medical conditions which are both objectively quantifiable (20–22), and independent of significant changes in mean velocity (23). Pulsatility Index is an example of a well known and widely clinically utilized morphological TCD variable (24); one which evidence suggests is not useful for detecting LVO (15).

Toward the aim of quantifying TCD waveform morphology for the purpose of LVO identification, we have recently proposed a diagnostic biomarker termed Velocity Curvature Index (VCI) (25, 26). Mathematically, it is a straightforward extension of the concept of graph curvature; one which is sensitive to the morphological structure of the pathological waveforms first described by Demchuk et al. (10). In this work we retrospectively compare the diagnostic utility of VCI to that of a standard Velocity Asymmetry Index (VAI) for the detection of LVO in a clinical subject population collected in-hospital. Additionally, we evaluate a simple decision tree classifier designed to incorporate complimentary information from both metrics. Decision trees are a commonly used diagnostic methodology in several areas of medicine (27, 28), which have previously been used with TCD variables to screen for cervical vascular injury (29). To these ends, we employ leave-one-out cross validation and subsequent sensitivity analysis to assess performance as diagnostic thresholds are weighted toward detection of true positives. Our goal is the validation of TCD-based decision criteria which are objective, intuitive, and easily communicated; allowing physicians and first responders alike a common language for LVO assessment.

### MATERIALS AND METHODS

### Subject Examination and Imaging

We acquired TCD waveforms from two clinical populations enrolled consecutively at Erlanger Health System's Southeast Regional Stroke Center in Chattanooga, TN, from October 2016 through September 2017. The LVO group was comprised of patients with CTA-confirmed occlusion of the M1/M2 branches of the MCA and/or ICA (proximal extracranial or terminal intracranial segments); these occlusion locations being selected since they are the large cerebral arteries most amenable to neurovascular intervention. The In-Hospital Control group (IHC) consisted of patients who arrived at the hospital presenting with stroke symptoms, but were later confirmed negative for LVO by imaging. Patients in both groups received TCD examinations in addition to standard care (patient history, monitoring, pharmaceuticals, and CT/A/perfusion imaging). Patients were eligible for enrollment if a complete TCD exam was acquired within 4 h of CTA, and none of the following exclusions applied: (1) Head CT findings consistent with acute primary intracranial hemorrhage, (2) Hemodynamically unstable patients requiring pharmacological support for hypotension; (3) Subjects who underwent partial or full craniotomy; (4) Additional intracranial pathologies present (tumor, hydrocephalus, etc.); (5) Anticipated insufficient time to acquire a complete set of scan as described by the protocol; (6) Significant hemodynamic pharmacological agent (cocaine, amphetamine, etc.); (7) Subjects who are under arrest for a felony.

CTA examinations were performed using a GE Lightspeed VCT 64-section multidetector scanner (GE Healthcare, Milwaukee, WI) with a slice thickness of 0.625 mm, and bolus injection of 70–150 mL of Omnipaque 350 (GE Healthcare, Milwaukee, WI) contrast material (4.0 mL/s). CTA images were reformatted in the coronal and sagittal plane, and 10-mm maximum intensity projection reconstructions were rendered and sent to PACS for review. Occlusion location was determined by the radiologist on call, who was blinded to any results of the TCD examination.

Complete TCD examinations included (at minimum) one pair of left/right MCA scans at depths between 45 and 60 mm, each containing 15 or more individual beat waveforms (see **Figure 1**). Subjects for whom complete examinations were not obtained were counted as missing/indeterminate data and excluded from analysis. Data was acquired during available intervals between patient testing/treatment, and in no way impacted patient care. The TCD technician was often present during initial evaluation of the subject, and so was not entirely blinded to all clinical information or imaging results. Sample size was not predetermined for this feasibility study, being established pragmatically as the maximum number of subjects attainable in the allotted time frame. Experiment protocols were approved by University of Tennessee College of Medicine Institutional Review Board (ID: 16-097). Reporting in this manuscript is structured in accordance with the Standards for Reporting of Diagnostic Accuracy Studies [(30); see **Appendix** for detailed criteria checklist].

#### Waveform Processing and Feature Extraction Recording

TCD scans were acquired by a trained technician using 2 MHz hand-held ultrasound probes. CBFV signals associated with the left/right MCA were identified via insonation through the transtemporal window. CBFV envelopes were digitally sampled at 125 Hz and recorded throughout the entire exam. Once the CBFV signal was identified and optimized at a specific depth, waveform recordings were then made in 30-s intervals. The technician was instructed to obtain recordings for as many depths as possible between 45 and 60 mm in both the left/right cerebral hemispheres. Start times for each interval were marked by the technician using a custom event remote, which prompted a 30 s countdown to a corresponding stop event. TCD envelopes and event information were aligned using custom software (Python 2.7; Kivy 1.9) running on Windows 10.

#### Processing

Average beat waveforms from each recorded depth interval were extracted using a combination automated beat identification algorithm with manual checking/editing. In this procedure, individual beats within each interval were first identified automatically using an internally developed beat extraction tool, and displayed to the user for manual confirmation/editing. Detected beats which lacked clear pulsatile structure and/or deviated anomalously from the group average (usually due to probe displacement during recording), were excluded. The remaining beats were then aligned and averaged, resulting in a single representative beat waveform for each recorded depth interval (see examples in **Figure 1**).

Since Doppler velocities scale with the cosine of the incident angle between the ultrasound beam and underlying blood flow (31), TCD waveforms for a given vessel with the highest measured velocities are assumed to most accurately reflect reality. In line with this reasoning, for each subject we selected a single bilateral (left/right) pair of average beat waveforms for analysis consisting of those with maximal mean velocity across all recorded depths for each hemisphere.

#### VCI

Curvature is a well-defined mathematical property of space curves which quantifies the degree to which a curve deviates from being "straight" at a given point. VCI is an application of the curvature metric specific to TCD which quantifies the degree to which a beat is blunted and/or dampened. Since curvature is a nonlinear function sensitive to small inflections associated with high frequency noise, we first smooth the average beat waveform via convolution with the Hanning window (9 ms). Moreover, we elect to consider only curvature associated with the beat systolic complex, where the signal-to-noise ratio is presumably greatest. The systolic complex, or "beat canopy," comprises the proportion of the beat with the highest velocities and richest morphological structure.

To compute VCI for a given TCD beat waveform, curvature is first computed for each time point (ti) of the smoothed beat (denoted x(ti) below) via the discretized equation for graph curvature (equation 1) expressed in terms of finite differences. 1 and δ 2 in equation 1 denote the first order backward (equation 2) and second order central (equation 3) finite difference equations. VCI, defined by equation 5, is computed as the sum of curvature taken over all individual time points comprising the beat canopy (C). The beat canopy is defined in equation 4

as the set of time points wherein velocity exceeds one quarter of its total diastolic-systolic range (td, and ts denoting the time points corresponding to diastolic minimum and systolic max, respectively). Since the hypothesized effect of occlusion on the TCD waveform is to lower VCI in the occluded vessel, when assessing a bilateral pair of waveforms we take VCI as the minimum computed for each member of the pair. VCI is a positive metric which, in principle, has no upper bound, but is clearly bounded in practice (see **Figure 2A**).

$$k\left(t\_i\right) = \frac{\left|\delta^2\left[\mathbb{X}\right](t\_i)\right|}{\left(1 + \left(\Delta\left[\mathbb{X}\right](t\_i)\right)^2\right)^{\frac{1}{2}}}\tag{1}$$

$$
\Delta\left[\boldsymbol{x}\right](t\_i) = \boldsymbol{x}(t\_i) - \boldsymbol{x}(t\_{i-1})\tag{2}
$$

$$
\delta^2\left[\mathbf{x}\right](t\_i) = \mathbf{x}\left(t\_{i+1}\right) - 2\mathbf{x}\left(t\_i\right) + \mathbf{x}\left(t\_{i-1}\right) \tag{3}
$$

$$C = \left\{ i \colon \mathfrak{x}\left(t\_i\right) \geqslant \mathfrak{x}\left(t\_d\right) + \frac{\mathfrak{x}\left(t\_s\right) - \mathfrak{x}\left(t\_d\right)}{4} \right\} \tag{4}$$

$$VCI = \sum\_{i \in C} k\left(t\_i\right) \tag{5}$$

#### VAI

Velocity Asymmetry Index is a metric which quantifies the degree to which average CBFV observed for a vessel in a given cerebral hemisphere differ from that observed in the corresponding vessel in the opposite hemisphere (see LVO example in **Figure 1**). The hypothesis that CBFV in an occluded vessel may be lower than that of the corresponding unaffected hemisphere is intuitive, but also supported by previous work (18). For a bilateral pair of left/right average beat waveforms, denoted xL(t) (with NL total time points), and xR(t) (with NR time points) in equations 6 and 7, respectively, VAI (defined in equation 8) is computed as the minimum average velocity across hemispheres divided by the corresponding maximum. By definition, VAI is a positive definite metric bounded on the closed interval [0, 1].

$$
\mu\_L = \frac{1}{N\_L} \sum\_{i=1}^{N\_L} \varkappa\_L(t\_i) \tag{6}
$$

$$
\mu\_{\mathcal{R}} = \frac{1}{N\_{\mathcal{R}}} \sum\_{i=1}^{N\_{\mathcal{R}}} \varkappa\_{\mathcal{R}}(t\_i) \tag{7}
$$

$$VAI = \frac{\min\left(\{\mu\_L, \mu\_R\}\right)}{\max\left(\{\mu\_L, \mu\_R\}\right)}\tag{8}$$

#### Feature Statistical Analysis

For both VCI and VAI, resultant group samples were not normally distributed. Accordingly, we tested significance of observed differences in group distributions for each feature using the Mann-Whitney U test statistic. The U-statistic is directly proportional to the common language effect size (by a factor of the product of the group sample sizes under comparison), which is equivalent to the area under the Receiver Operating Characteristic curve (ROC-AUC). Additionally, we computed ROC curves detailing separability of subject group distributions (LVO vs. IHC) for each feature. Specifically, the ROC curves give

True Positive Rate (TPR) as a function of False Positive Rate (FPR) for each possible feature threshold. To quantify binary classification performance, we computed SEN, SPE, and accuracy (ACC) at the thresholds which maximized Youden's J-statistic (32):

$$J = TPR - FPR\tag{9}$$

Additionally, we bootstrapped 95% confidence intervals (CI) on group means for each feature; accomplished by iteratively resampling each group distribution with replacement (10,000 iterations), and each time taking the mean (CI given by the 2.5th and 97.5th percentile of the resultant empirical distribution) (33). Statistical tests and ROC curves were computed using standard python libraries; SciPy version 1.0, and Scikit-learn version 0.19.1 (34), respectively.

### LVO Classification

#### Decision Tree Classifier

We sought to combine VCI and VAI into a single binary classifier with simple and intuitive decision criteria. The approach adopted here is to augment the bilateral VCI assessment such that subject pairs with VCI less than some low critical threshold are classified as LVO, whereas pairs exceeding some high critical threshold are classified as IHC. Pairs observed to fall between these thresholds are deemed uncertain and decided then based on VAI. This procedure effectively partitions the 2D decision space into two subspaces with piece-wise linear boundaries (see **Figure 2**). The procedure for fitting the thresholds based on a given set of training data were as follows. First, the low VCI threshold (VCIMIN) was fit using all the training data. Subjects with paired VCI below the threshold were predicted as LVO, and set aside. Next, the high VCI threshold (VCIMAX) was fit from the remaining data. Subjects with supra-threshold VCI were predicted as IHC, and set aside. Finally, the remaining data was used to fit the VAI threshold (VAICRIT), with sub/supra-threshold subjects predicted as LVO/IHC, respectively. Specifically, each of the three thresholds (VCIMIN, VCIMAX, and VAICRIT) were fit as the threshold which maximized Youden's J statistic for the data applicable to each decision (see also the sensitivity weighted Jstatistic used to determine thresholds for sensitivity analysis in section Sensitivity Analysis).

#### Model Cross-Validation

For both individual diagnostic metrics (VAI and VCI), as well as the decision tree model, leave-one-out cross-validation was performed to assess generalization of classification performance near decision boundaries. In this iterative procedure, a single subject is removed from the pooled data, and the predictive model is derived via training (i.e., threshold optimization) with the remaining subjects. The excluded subject's data is then predicted using the trained model, and this procedure is repeated for each subject to obtain a complete set of crossvalidated predictions from which to assess binary classification performance metrics (SEN/SPE/ACC).

#### Sensitivity Analysis

For many clinical problems, including LVO detection, diagnostic net benefit is optimized by increased detection of true positives at the cost of missing true negatives (i.e., SEN is prioritized over SPE). However, poor diagnostics for which SEN is maximized often have no clinical value, as SPE may plummet and overall ACC approaches chance. In order to assess how performance characteristics of our classifiers changed when priority is weighted toward increased sensitivity, we performed a sensitivity analysis wherein we iterated cross-validation of each model, each time incrementing classification thresholds away from the starting point of Youden's maximal J, toward increasing sensitivity. The procedure can be conceptualized as simply adjusting the thresholds up along the associated ROC curves toward increased true positive rate. Mathematically, this was accomplished by introducing a parameter (α) to modify the formula for Youden's J statistic as given in formula 9, and choosing thresholds to maximize the resultant index (Jα). Classifier performance was assessed by cross-validating each model with α ranging from 0.5 (threshold equivalent to Youden's maximal J) to 1 (maximal sensitivity) in steps of 0.01.

$$\mathbf{J}\_{\alpha} = \alpha \mathbf{TPR} - (1 - \alpha) \mathbf{FPR} \tag{10}$$

#### RESULTS

#### Subject Demographics

The current analyses included 33 LVO subjects (16 female), and 33 IHC subjects (13 female), with average ages of 66.9 years (SD = 15.7), and 56.4 years (SD = 16.3), respectively. A total of 156 subject screenings were attempted at Erlanger Medical Center, of which 68 were excluded due to screening failures (time required to complete exam, subject compliance, etc.). Of subjects with sufficient initial screenings, 50 and 38 were initially enrolled in the LVO and IHC groups, respectively. Of the LVO subjects, 3 were discontinued (subject either expressed desire to discontinue, or was transferred or died before enrollment could be completed). An additional 14 LVO and 5 IHC subjects were subsequently excluded due to the presence of disqualifying criteria unknown at the time of enrollment. In the LVO group, there were 20 subjects with M1 occlusions, 3 with M2 occlusions, and 8 with ICA occlusions. An additional subject had dual occlusions of both the M1 and ICA (same hemisphere), and another additional subject had bilateral occlusions of both ICA, in addition to an M2 occlusion. TCD exams were performed an average of 33 min (SD = 20) post-CTA for IHC subjects compared to 43 min (SD = 44) post-CTA for LVO subjects (difference not significant between groups; t = −1.15, p = 0.26). At time of admittance, LVO subjects were more physiologically/cognitively impaired as assessed by National Institute of Health Stroke Scale (NIHSS), with average scores of 16.8 (SD = 6.6), compared to 1.9 (SD = 2) for IHC (differences strongly significant between groups; t = −12.2; p << 0.001). No adverse events were reported for any subjects as a result of TCD examination.

#### Individual Metric Statistical Comparisons

**Figure 3** shows VAI and VCI metric distributions for LVO and IHC groups (A, C), and associated ROC curves depicting separability of group metrics (B, D). VAI means were greater for IHC subjects (0.89, CI = 0.86–0.92) relative to LVO (0.65, CI = 0.58–0.72). Associated ROC-AUC was observed at 88.4%, with significant group distribution differences confirmed by Mann-Whitney testing (p << 0.001). Similarly, VCI means were greater for IHC subjects (4.95, CI = 4.55–5.36) relative to LVO (2.66, CI = 2.38–2.97); with associated ROC-AUC observed at 94.2%. Significant group distribution differences were likewise confirmed by Mann-Whitney testing (p << 0.001). SEN/SPE/ACC at thresholds corresponding to Youden's maximal J are detailed in **Table 1** for both metrics.

#### Sensitivity Analysis

**Figure 4** shows SEN, SPE, and ACC dependence on the alpha weighting parameter using leave-one-out cross-validation for each classifier. By definition, the sensitivity of each classifier increases with increased alpha. Interestingly, performance

FIGURE 3 | Group feature distributions (A,C) were significantly different for both metrics (*p* << 0.001). Associated ROC curves (B,D) confirm both VAI and VCI provide diagnostically relevant information concerning the presence of LVO, with the greater separability observed for VCI suggesting it more information rich.

TABLE 1 | Descriptive information and performance indicators comparing LVO and IHC groups for VAI and VCI metrics at thresholds corresponding to Youden's maximal J statistic.

*In column 3, CI refers to 95 percent Confidence Intervals around the respective means given in column 2.*

indicator trajectories vary substantially between the VAI classifier and the other two (VCI and decision tree). For VAI, observed SPE and ACC are optimal near the maximal J (alpha = 0.5), and rapidly degrade with increased prioritization of SEN corresponding to alpha greater than 0.6. In contrast, for VCI and the decision tree, a stable range of alpha exists away from the maximal J for which SPE, and ACC are optimized and all performance indicators are stable. This range, roughly 0.6–0.8, is indicated in gray in **Figure 4**. Above this range, SEN increases for VCI and decision tree are accompanied by precipitous decreases in ACC and SPE. Together these results suggest a natural optimal alpha range for prioritization of sensitivity for the VCI and decision tree classifiers.

**Figure 5** shows cross-validated confusion matrices for each classifier with alpha specified at 0.6, which represents the start of the optimal range for the VCI and decision tree classifiers, and the tail end of the optimal range for the VAI classifier. For VAI, we observed an overall accuracy of 79%, with SEN/SPE of 76%, and 82%, respectively. For VCI, we observed an overall accuracy of 88%, with SEN/SPE of 88%, and 88%, respectively. Finally, for the decision tree we observed an overall accuracy of 91%, with SEN/SPE of 94%, and 88%, respectively. Together, these results demonstrate the superiority of the VCI classifier relative to VAI. However, within the framework of the decision tree, VAI helped to increase SEN of LVO identification relative to VCI alone. **Figure 5** results are summarized in **Table 2**.

### DISCUSSION

To our knowledge, this work represents the first published LVO decision criteria based on TCD variables which can be computed algorithmically and interpreted objectively. Results from all classifiers fall into the range observed in previous TCD studies using complex multi-vessel recording protocols (13, 17). Moreover, previous studies using predictive variables amenable to ROC analysis have not been subject to cross-validation in the manner we have performed here. Most importantly, these metrics substantially outperform stroke severity scales currently in prehospital use, which recent reviews suggest are unlikely to predict LVO with both high sensitivity and specificity (6, 7). Specifically, (6) published performance indicators for 5 such stroke assessment scales (RACE, 3ISS, LAMS, CPSSS, and PASS); reporting ACC and SEN capped at 74 and 64% across all scales. A sense of the potential for improvement upon these numbers can be gleaned from our simple decision tree which achieved cross-validated ACC and SEN exceeding 90%.

Previous work assessing TCD efficacy in detecting LVO align well with our current results. Tsivgoulis et al. (17) detected occlusions and stenoses based on the presence of the pathological waveforms described by Demchuk et al. (10) with SEN/SPE of 79 and 94%, respectively. A similar exam protocol was

FIGURE 4 | Cross-validated performance indicators for the VAI (A) and VCI (B) metrics as well as combined decision tree classifier (C). Sensitivity increases with the alpha weighting parameter as specificity decreases. VAI specificity decreases rapidly with increased sensitivity, whereas VCI and the decision tree display a stable range (indicated in gray) wherein specificity and accuracy are optimal.

TABLE 2 | Performance indicators for leave-one-out cross-validated classifiers comparing LVO and IHC groups, with classification thresholds specified at alpha equal to 0.6.

used by Brunser et al. (13), but with additional power Mmode criteria also considered, which achieved SEN/SPE of 90 and 94% detecting occlusion of any non-specific artery. The sensitivity and specificity observed for our metrics compare reasonably well with these previous results, which is especially encouraging considering our features were extracted from bilateral examination of a single vessel. Reliance on the MCA signals is pragmatic, but also represents a notable opportunity for improvement upon our current experimental paradigm. The MCA possesses the longest expected segment of probeable depths, and is thus most easily insonated and reliably located, but there is clearly more diagnostic information available in other vessels in the form of relative morphology and collateral flow which may improve performance in future experiments.

It is notable that the performance indicator curves we observed for VCI and the decision tree were extremely similar. This is a natural consequence of the decision criteria, which dictates that VAI is used only to decide "uncertain" subjects, thus serving mainly to improve upon VCI sensitivity to the degree allowable by the training data. Further work is needed to determine if there are specific occlusion types or patient demographics for which each metric works particularly better or worse, which should help to optimize decision criterion. In this configuration, VCI is doing the "heavy lifting" in our decision tree. It is an effective diagnostic metric because it is sensitive to the morphological structure of Demchuk's minimal, blunted, and dampened flows (10). The blunted waveform, for example, possesses an inherently smooth (i.e., low curvature) systolic complex, and is thus readily quantified by VCI. Forthcoming work will analyze in specific detail the manner in which VCI captures the subtle morphological variations associated with pathological LVO waveforms.

Some limitations of our study and directions for future work should be noted. The primary factor which could potentially inhibit generalisability of our results is the small sample size of our study. Much further data is required to refine estimates of morphological variability inherent in LVO and clinical control patient populations. Additionally, numerous important subgroup analyses are required to determine if/how TCD morphology depends on demographic and clinical factors (age, gender, occlusion type, etc). It should also be noted that the TCD technician's exposure to patient clinical information represents a potential source of bias which should be mitigated in future work by more thorough blinding. Finally, the relationships between curvature, heart rate, and stroke pathology require further investigation. In our sample, LVO subjects had an average heart rate of 87.9 (SD = 22.2) beats per minute (BPM), vs. 71.2 (SD = 11) BPM for IHC; which was significantly different between groups (t = 3.8, p < 0.001). However, there was no significant correlation between heart rate and curvature within either subject group (r = −0.006, p = 0.9 for IHC; r = −0.29, p = 0.1 for LVO). Moreover, when we use heart rate itself as a predictor to distinguish between groups, we observe an AUC of 73%, considerably underperforming both the VAI (88% AUC), and VCI (94% AUC) metrics. So, while it is possible that heart rate accounts for some degree of variance between groups, it remains unclear whether the effect is causal or correlative. It is certainly plausible that the lack of blood supply characteristic of LVO causes heart rate to increase; meaning heart rate is effectively part of the diagnostic signal. Further work is needed to establish how elevated heart rate might affect VCI when occlusion is not present.

One strength of the current approach is simplicity of data acquisition and communicability of decision criteria. However, as the amount of data we acquire increases, the subtle variations we will ultimately wish to detect will undoubtedly require more complicated and abstracted models. Given the efficacy of initial results, the road map to such models is encouraging. Incorporation of depth dependent and interhemispheric morphological dynamics across multiple vessels might ultimately allow precise prehospital localization of occlusion, and distinction between occlusive and hemorrhagic strokes. Moreover, digital rendering of individual subject vasculature, currently possible with emergent technologies such as 3D MRA time of flight imaging, could facilitate ultra rapid mapping and scanning across multiple vessels, as well as development of anatomically realistic mathematical models of cerebral hemodynamics. Such models could dramatically increase our understanding of the fluid mechanics involved in vascular occlusion, and the associated impact on morphological biomarkers like VCI.

### CONCLUSIONS

Our results suggest both VAI and VCI contain robust information concerning the presence of intracranial occlusion. Both are objective and real-time computable, and thus represent promising candidate metrics for the development of TCD-based prehospital LVO diagnostic systems. The feature distributions and classification performance indicators we observed suggest VCI may be superior to VAI for LVO detection, but even simple approaches to feature combination such as the decision tree analyzed herein may serve to further increase diagnostic accuracy. More data is needed to determine how well these decision criteria scale and generalize to a wider range of subject demographics and pathologies. Nonetheless, these results demonstrate the foundational potential for machine-learning approaches to TCD morphological analysis to enable faster and more widespread access to life saving medical intervention in the future.

#### ETHICS STATEMENT

This study was carried out in accordance with the recommendations of Declaration of Helsinki and The University of Tennessee College of Medicine Institutional Review Board with written informed consent from all subjects. All subjects gave written informed consent in accordance with the Declaration of Helsinki. The protocol was approved by the University of Tennessee College of Medicine Institutional Review Board (ID: 16-097).

### AUTHOR CONTRIBUTIONS

ST had access to all data used in the study and takes responsibility for data integrity and accuracy of data analysis. SW, RH, and TD study concept and design. ST and CT analysis. ST, CT, and NC interpretation of data. ST drafting of the manuscript. CT, NC, SW, RH, and TD critical revision of the manuscript for

### REFERENCES

important intellectual content. ST statistical analysis. CT, NC, RH, SW, and TD technical or material support. RH and TD study supervision.

### ACKNOWLEDGMENTS

We acknowledge and thank Leonardo Martinez-Torres, James LeVangie, and Professor Yaser Abu-Mostafa for their profound contributions to this work. In addition, the supporting efforts of Ben Delay, Michael O'Brien, Mina Ranjbaran, Shankar Radhakrishnan, Leo Petrossian, and Kian Jalaleddini were greatly appreciated.

#### SUPPLEMENTARY MATERIAL

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

response to intravenous thrombolysis for stroke. Stroke (2007) 38:948–54. doi: 10.1161/01.STR.0000257304.21967.ba

waveform morphology in healthy volunteers. Clin Sci. (2006) 111:47–52. doi: 10.1042/CS20050365

**Conflict of Interest Statement:** At the time that this research was conducted, ST, CT, NC, SW, and RH were employees of Neural Analytics, Inc., and TD was a paid consultant. All authors either hold stock or stock options in Neural Analytics, Inc.

Copyright © 2018 Thorpe, Thibeault, Canac, Wilk, Devlin and Hamilton. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Machine Learning in Acute Ischemic Stroke Neuroimaging

#### Haris Kamal\*, Victor Lopez and Sunil A. Sheth

Department of Neurology, University of Texas at Houston Health Science Center, Houston, TX, United States

Machine Learning (ML) through pattern recognition algorithms is currently becoming an essential aid for the diagnosis, treatment, and prediction of complications and patient outcomes in a number of neurological diseases. The evaluation and treatment of Acute Ischemic Stroke (AIS) have experienced a significant advancement over the past few years, increasingly requiring the use of neuroimaging for decision-making. In this review, we offer an insight into the recent developments and applications of ML in neuroimaging focusing on acute ischemic stroke.

Keywords: stroke, neuroimaging, machine learning (artificial intelligence), neurosciences, support vector machina (SVM), stroke management, stroke diagnosis

#### Edited by:

David S. Liebeskind, University of California, Los Angeles, United States

#### Reviewed by:

Jens Fiehler, Universitätsklinikum Hamburg-Eppendorf, Germany Maurizio Acampa, Azienda Ospedaliera Universitaria Senese, Italy

\*Correspondence:

Haris Kamal<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Stroke, a section of the journal Frontiers in Neurology

Received: 26 June 2018 Accepted: 22 October 2018 Published: 08 November 2018

#### Citation:

Kamal H, Lopez V and Sheth SA (2018) Machine Learning in Acute Ischemic Stroke Neuroimaging. Front. Neurol. 9:945. doi: 10.3389/fneur.2018.00945

Machine Learning (ML), considered a branch of artificial intelligence, is a field of computer science and engineering that facilitates extraction of data based on pattern recognition. A computer learns from previous mistakes after repeated analysis of data and masters tasks that were previously considered too complex for a machine to process (1). The development of these systems to interpret data in neuroimaging has provided valuable information for research in matters of the interaction, structure, and mechanisms of the brain and behavior in certain neurological disorders (2, 3).

Machine learning systems are now being implemented in the clinical neurosciences to devise imaging-based diagnostic and classification systems of neoplasms of the brain (4–6), certain psychiatric disorders (7–11), epilepsy (12, 13), neurodegenerative disorders (14–20), and demyelinating disorders (21–23). In this review, we discuss the present-day role of ML focusing on acute ischemic stroke (AIS), discussing its potential and limitations.

### MACHINE LEARNING IN THE CLINICAL NEUROSCIENCES

The use of neuroimaging in the evaluation of many neurological diseases such as dementia, epilepsy, demyelinating diseases, depression, and schizophrenia has grown tremendously. This burgeoning interest has been met with an expansion of ML algorithms in neurosciences (1, 24).

Oliveira et al. (14) evaluated an unsupervised ν-One-Class Support Vector Machine (ν-OC-SVM) trained with neuroimaging variables, such as cortical thickness and cerebral volume of the brain, from healthy subjects to calculate an abnormality index and compare it with patients diagnosed with mild cognitive impairment (MCI) and Alzheimer's disease (AD). The method correctly classified AD subjects as outliers with an accuracy of 84.3%, and the brain abnormality index was directly associated with the group diagnosis, clinical data, biomarkers, and risk of future conversion to AD.

In schizophrenia, Greenstein et al. (9) used Random Forest (RF), a machine learning algorithm, to discriminate between childhood-onset schizophrenia and healthy patients based on brain magnetic resonance imaging (MRI) measurements of regions of interest (ROI): left temporal lobes, bilateral dorsolateral prefrontal regions, and left medial parietal lobes. The algorithm correctly classified groups with 73.7% accuracy, and a greater brain-based probability of illness

**93**

was associated with a statistically significant worse functioning and fewer developmental delays. Machine learning can also help distinguish between subsets of a certain disease. Bleich-Cohen et al. (7) utilized Searchlight Based Feature Extraction (SBFE), a data-driven multi-voxel pattern analysis (MVPA) approach, to search for activation clusters of cognitive loads in brain functional Magnetic Resonance Imaging (fMRI). This ML method helped to identify the two subgroups of schizophrenic patients with and without Obsessive-Compulsive Disorder (OCD) with a 91% accuracy, successfully delineating between symptom severity and a psychiatric comorbidity.

An et al. (12) compared whole-brain white matter changes in patients with mesial temporal epilepsy and matched healthy controls, evaluating tract-based spatial statistics and fractional anisotropy with an ML approach. This ML-based approach discriminated each group accurately and demonstrated high sensitivity to changes in fractional anisotropy in mesial temporal epilepsy patients, which may be beneficial when no lesion can be identified on neuroimaging. Moghim et al. (13) introduced a predictive model for seizure occurrence in a single patient. This approach was based on a multi-class support vector machine (SVM) and 14 selected features of an electroencephalogram in patients with epilepsy. The predicted time of seizure with a window between 20 and 25 min was reported with an average sensitivity of 90.15, 99.44% specificity, and 97% accuracy.

Lesion burden estimation in traumatic brain injury (TBI), AIS, dementia, and multiple sclerosis serves to identify the affected regions, the extent of damage, and therefore, the functional outcome in such patients. Kaminatas et al. (25) proposed an approach for lesion segmentation using a multimodal brain MRI based on an 11-layers deep, multi-scale, 3D Convolutional Neural Networks (CNN) called Deep Medic. Their proposed novel training scheme is based on two main components, a 3D CNN that produces accurate soft segmentation maps and a connected Conditional Random Field that imposes regularization constraints on the CNN output and produces the final hard segmentation labels. This allows for a deeper and more discriminative delimitation of lesion burden, with the highest reported accuracy observed in a cohort of patients with severe TBI.

### CHALLENGES IN ACUTE ISCHEMIC STROKE

Stroke is the leading cause of serious long-term disability and the fifth leading cause of death in the United States, with its prevalence increasing with advancing age in both males and females, as each year ∼795,000 Americans experience a new or recurrent stroke (26). This burden is coupled with a direct medical expense of an estimated \$23.6 billion according to the last annual report of 2014 (26). With the increasing complexity of the acute ischemic stroke therapy and the rising of per-person costs, there is a real and urgent need for a technological solution to aid in the streamlined care of patients and selection of the appropriate therapeutic intervention.

Present treatments for AIS revolve around rapid reperfusion of ischemic tissue, using intravenous (IV) thrombolytic medications such as tissue plasminogen activator (tPA) and/or endovascular techniques to mechanically remove the obstruction to blood flow. Contemporary clinical trials are now implementing a higher complexity of neuroimaging modalities to define treatment standards, resulting in an increased economic as well as logistical burden on healthcare. The WAKE-UP multicenter clinical trial (27) used magnetic resonance imaging (MRI) in patients that presented with an unknown time of onset of symptoms to identify brain regions that exhibited a restricted diffusion on diffusion-weighted imaging (DWI) scan and no T2-signal hyperintensity on fluid-attenuated inversion recovery (FLAIR) sequence, estimating the onset of the infarct to be <4.5 h and thus guiding stroke therapy. Previous to this study, non-contrast head CT, an imaging modality that is widely and readily available, was the only imaging screen used to assess for tPA eligibility. The new, expanded tPA indication requiring MRI poses challenges for a majority of centers, which do not have ready access to this type of imaging emergently and 24/7.

The growing dependence on neuroimaging in determining treatment options for acute ischemic stroke is observed as well for endovascular stroke therapy (EST), which has shown to improve outcome when used in combination with standard medical care (28). In 2015, numerous clinical trials demonstrated a clear benefit of endovascular treatment over medical management alone for a select group of patients with acute ischemic stroke seen within 6 h of the onset of stroke (29–33), and relied on imaging modalities including NCHCT, CT/MR angiography (CTA/MRA) and CT/MRI perfusion (CTP/MRP) scans. Results derived from these trials showed an advantage in using advanced imaging modalities in identifying patients with a higher likelihood of better outcomes from EST. Two additional clinical trials (34), DAWN and DEFUSE3, published in 2018 evaluated a much larger population of stroke patients, those presenting up to 24 h after their symptoms, and required the use of perfusion imaging with CT or MRI.

This increased reliance on neuroimaging has led to a tremendous improvement in our ability to care for patients with AIS but has been coupled with a number of challenges. Specifically, limited availability of these imaging modalities, a shortage of specialists to promptly interpret these studies, as well as inter-observer variability have limited the implementation of the above findings. Indeed, studies evaluating inter-observer performance on Alberta Stroke Program Early CT Score (ASPECTS), a 10-region imaging grading system in stroke, showed significant variability (35–39). Further adding to the complexity of acute stroke treatment is that while the need to perform and interpret advanced neuroimaging has recently increased, the urgency with which such evaluation is being performed has remained the same. For every minute that a patient with a large vessel occlusion fails to be treated, an estimated 1.9 million neurons and 14 billion synapses are lost in the brain (40). Trials evaluating efforts to promptly assess and treat patients with AIS have demonstrated superior outcomes and decreased morbidity. In patients treated with

intravenous thrombolysis, reducing treatment times by 15 min was associated with reduced in-hospital mortality, reduced incidence of symptomatic intracranial hemorrhage, and a greater likelihood of independent ambulation at discharge (41, 42). In patients treated with endovascular therapy, for every 15-min reduction of onset to recanalization of the occluded artery, 34 per 1,000 treated patients had improved disability outcome (43). As such, there is an urgent need for systems to rapidly and precisely interpret neuroimaging data in acute ischemic stroke.

### IMPLEMENTATIONS OF MACHINE LEARNING IN ACUTE ISCHEMIC STROKE

Machine learning algorithms have been used to assist in the diagnosis and individualized treatment decisions in acute ischemic stroke. The implementations of machine learning are numerous, from early identification of imaging diagnostic findings (44), estimating time of onset (27, 45), lesion segmentation (46), and fate of salvageable tissue (47, 48), to the analysis of cerebral edema (49, 50), and predicting complications (51–53) and patient outcomes (54–57) after treatment. A summary of the most recent articles investigating the applications of machine learning for automated diagnosis and outcome prediction in acute ischemic stroke is given in **Table 1**.

One of the most relevant clinical criteria to decide if a patient with an acute ischemic stroke is eligible for IV thrombolysis with tPA is a time from symptom onset of <4.5 h, but in medical practice, stroke symptom onset is usually unknown. Ho et al. (45) developed a deep learning algorithm based on an autoencoder architecture to extract imaging features in perfusion-weighted images (PWI) in MRI to determine the time elapsed since stroke onset.

Lesion estimation and identification of salvageable tissue are essential in the acute decision making in stroke, but the expense and resources involved present a challenge for physicians. Chen et al. (46) used a framework with two CNNs to segment stroke lesions using DWI in MRI. One CNN was a combination of two DeconvNets (EDD Net), and the second CNN was a multiscale convolutional label evaluation net (MUSCLE Net) to help reduce the potential false positives detected by the EDD Net. The dataset was built with clinical acquired DWI from 741 subjects, exhibiting a high lesion detection rate, and accuracy.

Measurement of the perfusion-diffusion mismatch and calculation of infarction probability using MRI-based approaches for tissue-at-risk evaluation can be applied in stroke treatment decisions. Bouts et al. (47) analyzed the ability of five algorithms to depict potentially salvageable tissue using MRI imaging from rats subjected to a right-sided MCA occlusion without subsequent reperfusion, and with spontaneous or thrombolysis-induced reperfusion. The highest accuracy of riskbased identification of acutely salvageable ischemic tissue that could recover on subsequent reperfusion was observed using a generalized linear model (Dice's similarity index = 0.79 ± 0.14). Similarly, Huang et al. (48) used an SVM to predict infarct on a pixel-by-pixel basis using acute cerebral blood flow (CBF) and apparent diffusion coefficient (ADC) on MRI data. Serial images were collected during the acute phase up to 3 h and again at 24 h from 12 rats in each of the stroke groups exposed to a 30-min, 60-min, or permanent middle cerebral artery (MCA) occlusion. The accuracy observed for this approach was high in all groups and was enhanced by adding neighboring pixel information and spatial infarction incidence.

Takahashi et al. (44) designed a method to identify a hyperdense MCA, also known as the MCA dot sign, an important evaluation in an NCHCT as it represents a thrombus in a vessel. The authors created ROIs around the Sylvian fissure region and identified MCA dots based on the morphologic top-hat transformation, and classified images using an SVM with four features. Two hundred and ninety-seven CT images from seven patients with an MCA dot sign were classified by an SVM system, which exhibited a maximum sensitivity of 97.5% at a false positive rate of 1.28 per image and 0.5 per hemisphere while assessing the MCA dot sign.

Another application of ML in AIS is predicting factors that will contribute to neurological deterioration and increased morbidity, such as cerebral edema. Chen et al. (49) proposed a machine learning algorithm using serial CT scans of stroke patients to delineate and measure cerebrospinal fluid (CSF) volume over time, as it may represent a sensitive biomarker of cerebral edema progression. The initial cohort consisted of 155 subjects and preliminary processing using a generalized estimating equations (GEE) model top to calculate CSF volumes over time, adjusting for age, demonstrated that a reduction in CSF volume from baseline to final CT was correlated with infarct volume, the presence of cerebral edema, and the degree of midline shift. Comparatively, Dhar et al. (50) validated an automated technique for intracranial CSF segmentation by an ensemble of RFbased machine learning with a geodesic active contour (GAC) segmentation. CSF spaces were outlined on scans performed within 6 h of stroke onset and then closest to 24 h later in 38 patients. This method accurately tracked changes in CSF volume with an average DSC > 0.7. Pearson correlation coefficients between the changes in CSF and the ground truth were found to be statistically significant. These algorithms represent a potential for future research and may serve as a biomarker of cerebral edema severity.

The outcome of acute ischemic stroke patients is dependent on therapy, and risks for complications should be considered when deciding for stroke therapy. Yu et al. (53) established a method to predict the location and extent of hemorrhagic transformation (HT) in stroke, the most severe complication following reperfusion therapy. PWI and DWI of 165 patients treated with reperfusion therapy in a stroke center were collected and analyzed using five machine learning approaches, with Kernel spectral regression exhibiting an accuracy of 83.7 ± 2.6%. A multi-center retrospective study (52) assessed the predictive power for hemorrhagic transformation of PWI in MRI. Dynamic T2- weighted perfusion MR images from 263 patients from four medical centers were collected and served as input for linear and nonlinear predictive models, the latter having an average accuracy >85% in predicting HT. In one study, Nielsen et al. (54) ran a deep learning convoluted neural network (CNNdeep) with 9 biomarkers as input to calculate lesion volume in patients

\*The results displayed for each article are the most accurate or relevant in matter of the machine learning approach utilized according to the author.

CNN, Convoluted Neural Network; GAC, Geodesic Active Contour; HT, Hemorrhagic transformation; MCA, Middle Cerebral Artery; mRS, modified Rankin Scale; RF, Rain Forest; sICH, symptomatic Intracranial Hemorrhage; SR-KDA, Spectral Regression Kernel Discriminant Analysis.

treated with IV tPA. Input data from 29 untreated patients and 35 patients that received IV tPA were compared. This model predicted final infarct volume with 88% accuracy, being superior to other models in this study. Bentley et al. (51) predicted the risk of symptomatic intracerebral hemorrhage (sICH) after IV thrombolysis therapy. CT images of 116 patients who were treated with IV tPA, 16 of which had sICH, were entered as inputs into an SVM along with clinical severity. They found a better prognostication of the SVM when compared to the traditional clinician-based prognostication tools such as Hemorrhage after thrombolysis (HAT), and Sugar, Early Infarct signs, Dense cerebral artery sign, Age, and NIHSS scores (SEDAN).

Machine learning algorithms based on structural and functional MR images as input may assist in predicting motor deficits in stroke patients. Forkert et al. (56) applied 12 SVM classification models in calculating the corresponding 30-day mRS score of ischemic stroke patients through parameters including lesion overlap from different brain regions, stroke laterality, and other optional features such as infarct volume, NIHSS at admission, and patient age. Superior mRS prediction was observed by integrating the optional features and providing stroke location information, with a multi-value mRS prediction accuracy of 56%, and a dichotomized mRS (0–2 vs. 3–5) prediction accuracy of 85%. In a study by Rondina et al. (57), a proposed model to predict upper extremity motor deficit in 50 stroke patients was developed from data on structural MRI instead of functional MRI. Lesion probability images were derived using patterns of voxels and was then compared to lesion load per ROI in predicting outcomes, with the former providing better results when multiple regions of interest such as a range of cortical and subcortical motor areas and corticospinal tract were analyzed.

### CURRENT CHALLENGES AND FUTURE DIRECTIONS IN MACHINE LEARNING FOR ACUTE ISCHEMIC STROKE

Early promising results have demonstrated that ML techniques may be useful as decision support tools in treatment choices for AIS. To improve the generalizability of the findings discussed above, however, there are a number of limitations in currently existing architectures that need to be addressed. The first limitation is that of sample size. Deep learning algorithms using medical imaging often require datasets of tremendous magnitude, the types of which may not be readily available. For example, an ML algorithm demonstrated superior performance at differentiating skin cancer lesions from their benign corresponding equivalent when compared against 21 board-certified dermatologists, using a dataset of nearly 130,000 images (58). A dataset of this size in AIS for public use does not currently exist. This shortcoming, however, has been recognized as a problem that can and ought to be solved, and multiple calls for the creation of such a repository have been made (59). The obstacles in inter-institutional data sharing, as well as a lack of funding to correctly pre-process and curate these images, along limitations to host such a dataset account for some of the delays in the creation of this repository.

Another limitation encountered in neuroimaging-based ML techniques is the need for labeling regions of interest or "gold standard" findings on the images. That is to say, beyond collecting the images, the images and the findings on the images would need to be identified for the question being evaluated. For example, a study evaluating the presence or absence of a hyperdense MCA would need each image to be tagged with the true result, to train

the algorithm. Without foresight, this degree of manual curating could be required for each individual project.

#### CONCLUSION

Machine learning applications are expanding in the medical field for diagnostic and therapeutic purposes, and the rapidly expanding and increasingly neuro-imaging reliant field of AIS is proving to be fertile ground. There is a particular need for ML solutions in this field, which is faced with the challenge of increasingly complex data, with limited human

### REFERENCES

expert resources. Future directions in ML for AIS may require collaborative approaches across multiple institutions to build a robust dataset for efficient training of ML networks.

## AUTHOR CONTRIBUTIONS

HK and SS contributed equally to manuscript conception, design, revisions, and approved the submitted version. VL contributed to the manuscript design, revisions and edits and approved the submitted version.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Kamal, Lopez and Sheth. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Learning to Predict Ischemic Stroke Growth on Acute CT Perfusion Data by Interpolating Low-Dimensional Shape Representations

Christian Lucas 1,2 \*, André Kemmling3 , Nassim Bouteldja1 , Linda F. Aulmann4 , Amir Madany Mamlouk 5 and Mattias P. Heinrich1

1 Institute of Medical Informatics, University of Lübeck, Lübeck, Germany, 2 Graduate School for Computing in Medicine and Life Sciences, University of Lübeck, Lübeck, Germany, 3 Department of Clinical Radiology, University Hospital Münster, Münster, Germany, 4 Institute of Neuroradiology, University Medical Center Schleswig-Holstein, Lübeck, Germany, 5 Institute for Neuro- and Bioinformatics, University of Lübeck, Lübeck, Germany

#### Edited by:

Fabien Scalzo, University of California, Los Angeles, United States

#### Reviewed by:

David Robben, KU Leuven, Belgium Ulas Bagci, University of Central Florida, United States Aalpen A. Patel, Geisinger Health System, United States

> \*Correspondence: Christian Lucas

<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Stroke, a section of the journal Frontiers in Neurology

Received: 03 May 2018 Accepted: 02 November 2018 Published: 26 November 2018

#### Citation:

Lucas C, Kemmling A, Bouteldja N, Aulmann LF, Madany Mamlouk A and Heinrich MP (2018) Learning to Predict Ischemic Stroke Growth on Acute CT Perfusion Data by Interpolating Low-Dimensional Shape Representations. Front. Neurol. 9:989. doi: 10.3389/fneur.2018.00989 Cerebrovascular diseases, in particular ischemic stroke, are one of the leading global causes of death in developed countries. Perfusion CT and/or MRI are ideal imaging modalities for characterizing affected ischemic tissue in the hyper-acute phase. If infarct growth over time could be predicted accurately from functional acute imaging protocols together with advanced machine-learning based image analysis, the expected benefits of treatment options could be better weighted against potential risks. The quality of the outcome prediction by convolutional neural networks (CNNs) is so far limited, which indicates that even highly complex deep learning algorithms are not fully capable of directly learning physiological principles of tissue salvation through weak supervision due to a lack of data (e.g., follow-up segmentation). In this work, we address these current shortcomings by explicitly taking into account clinical expert knowledge in the form of segmentations of the core and its surrounding penumbra in acute CT perfusion images (CTP), that are trained to be represented in a low-dimensional non-linear shape space. Employing a multi-scale CNN (U-Net) together with a convolutional auto-encoder, we predict lesion tissue probabilities for new patients. The predictions are physiologically constrained to a shape embedding that encodes a continuous progression between the core and penumbra extents. The comparison to a simple interpolation in the original voxel space and an unconstrained CNN shows that the use of such a shape space can be advantageous to predict time-dependent growth of stroke lesions on acute perfusion data, yielding a Dice score overlap of 0.46 for predictions from expert segmentations of core and penumbra. Our interpolation method models monotone infarct growth robustly on a linear time scale to automatically predict clinically plausible tissue outcomes that may serve as a basis for more clinical measures such as the expected lesion volume increase and can support the decision making on treatment options and triage.

Keywords: ischemic, stroke, prediction, growth, learning, shape, CT, perfusion

### 1. INTRODUCTION

Cerebrovascular diseases, in particular strokes, are one of the leading global causes of death in developed countries (1). Acute stroke, which is usually caused by the blockage of cerebral blood flow due to a blood clot, is often diagnosed through CT or MR perfusion imaging (beside others, such as CTA). The derived perfusion parameter maps, e.g., Cerebral Blood Volume (CBV) or Time To Drain (TTD), provide spatio-temporal distributions of a contrast medium bolus within brain tissue. In contrast to native CT or standard MR sequences, such as T2 or FLAIR, perfusion images with their apparent functional signals enable the delineation of the potential infarct area even in the early acute phase and allow to visually assess the expected stroke severity, which helps the radiologist to come to a final therapy decision as early as possible.

In order to decide for a treatment the doctor has to weigh the risk of a therapy such as thrombolysis or thrombectomy against the expected outcome. For instance, (2) describe hemorrhages, such as symptomatic intracerebral hemorrhage, as a typical risk of intravenous thrombolysis therapy. For large vessel occlusions, mechanical thrombectomy improves functional outcomes but is logistically challenging. It is of major importance to consider the immediate availability of a therapy option, since the expected outcome strongly relies on the onset-to-treatment time (3). Depending on the expected time until revascularization, the radiologist has to estimate if further progression of the stroke can be avoided so that substantial parts of the tissue-at-risk within the penumbra could be salvaged. In the infarct core, however, as evident by a decrease in CBV, severe tissue injury and permanent vascular collapse have occurred.

Since stroke lesions vary widely in shape or size, and also evolve spatially heterogeneously over time, it is challenging for the radiologist to estimate growth or the size of the potentially stroke-affected tissue. For this reason, it is difficult to derive a time window in which a specific therapy path may be beneficial over another. Deep learning with CNNs has become popular in medical image analysis over the past recent years by clearly exceeding the so far state-of-the-art results, potentially capable of modeling this complex relationship.

#### 1.1. Objective

We present a novel tool for automatic stroke tissue outcome estimation using a CNN with a convolutional auto-encoder (CAE) that incorporates learned stroke shapes of core and penumbra. In a proof-of-concept, the trained model is able to predict the stroke lesion growth for patients with successful recanalization based on a given time-to-treatment for the thrombectomy (**Table 1**) and the CTP imaging parameter maps CBV and TTD. An evaluation of the method shows the practicality in principle on a limited dataset and the discussion provides pros and cons that suggest to further investigate this approach for clinical use.

#### 1.2. Outline

In order to gain a fundamental understanding of the method and its design choices, we provide a methodological overview in TABLE 1 | Inclusion criteria of the dataset for evaluation.

the following section: First, a review on the established stroke image analysis methods in clinical research literature is given; Second, the foundations of different image representations that can help solving higher-level tasks for image analysis as well as machine learning methods that have been investigated for stroke imaging are described; Third, we explain the use of CAEs for regularization in image segmentation by learning shape representations and how our work is based upon it. Subsequent, a detailed description of the assumptions and components of our method is provided—this third section also explains how to reproduce the method, that is, how to train the shape space and predict follow-up lesions based on noisy shape estimates. The fourth section lists the materials for a comparative evaluation and discusses its results, before we provide a final conclusion in the last section.

#### 2. IMAGE ANALYSIS

Classic thresholding methods for stroke image analysis have the drawback of only modeling a single univariate hard decision border between affected and non-affected tissue. Even when splitting into subgroups of different admission times (4) or distinguishing core against penumbra, the result will be a binary map of the affected vs. unaffected tissue. Further, purely voxel-based methods can produce irregular and physiologically implausible shapes. Statistical models, e.g., a linear regression model as used by Kemmling et al. (5), can cope with the variances in the data according to the complexity of the model and proper parameterization, while simple models will usually show a strong bias when used for high-dimensional problems.

#### 2.1. Representation and Spaces

In general, noisy images make it difficult to operate in image space, for instance, to apply a threshold to images for extracting regions that define the outcome. The input representation is not always suitable to detect the complex input patterns that determine the output. As known from signal processing and analysis, transforming the input into another representation (extract features) can often make it easier to perform classification or regression tasks. There are many transformations, e.g., non-linear kernel methods can bring a representation into a higher dimensionality where it may become linear separable. As this is a vast field that shall not be described here in detail, we emphasize that the input image data often needs to be fit into a regularized model or transformed to another representation first, on which the high-level task becomes easier to solve.

Kolouri et al. (6) propose transport spaces based on optimal transport theory to model biomedical problems such as tumor growth. The idea of looking at images as mass particle distributions is related to tissue distributions in biology, for example, when learning the transformations of some sample images onto a mean non-pathological image to extract the main modes of variance ideally representing the change from benign to pathological tissue (7). The modes could be extracted by principle component analysis (PCA) or other machine-learning approaches, e.g., auto-encoders (8). The formulation of the transport space and described applications suggest to make use of it for modeling other biological growth processes. However, to our knowledge, this has not yet been investigated for the cerebrovascular domain and applicability remains unclear for stroke tissue prediction.

While for all the above-mentioned methods their transformations to acquire a new representation are predefined, there are other models whose parameters can be estimated from samples. For instance, this could be learning a statistical distribution, like the shape and appearance of point representations (9), where the parameterization of the probabilistic distribution is learned from a training set. With a suitable representation at hand, there are several ways to machine-learn rather than fitting the input representation with its outcome into a statistical model. This often leads to more accurate results: Before Deep Learning has become popular in the last years, the medical image community had investigated Decision and Regression Forest models extensively, and they have shown good performance over statistical linear models or boosting approaches (10). However, these methods rely on previously specified or separately learned feature representations that need to be extracted from the image data first.

Opposite to the prior definition of the representations used, one can also machine-learn the representation and the classification (or regression) both at once using nonlinear artificial neural networks that are capable of learning sufficiently complex models without the need of tuning the right parameterization by hand. The review paper of Lee et al. (11) shows that the methods used with deep learning are still new to the field of stroke imaging and analysis. Some attempts with other machine learning methods have been made for diagnosis and prognosis, however, those models usually predict disease scores or specific clinical outcomes but not tissue outcome.

#### 2.2. Deep Learning for Stroke Imaging

Deep Learning with artificial neural networks is based on the idea of perceptrons where the output of a perceptron is computed by the weighted sum of its inputs x followed by an activation function σ (e.g., rectifier as in Equation 1). The power of such networks has been proven in the early work of Hornik et al. (12) by the fact that even a single hidden layer perceptron network with a proper activation function is capable of approximating any mathematical function. However, estimating that relation between input and output requires a lot of data and proper regularization since we have an underdetermined system when learning the coefficients wi (neuron weights for incoming connections), otherwise.

$$z = \sigma(\sum\_{i} w\_i \chi\_i), \quad \sigma(\boldsymbol{\wp}) = \max(0, \boldsymbol{\wp}) \tag{1}$$

Although known for a long time, training and regularization of such networks is difficult. As a consequence of this, there have only been few attempts to utilize them for spatial data such as medical images, e.g., as proposed by Huang et al. (13) for predicting tissue fates of stroke on acute image data, but their performance could not be tweaked to exceed other former stateof-the-art approaches. For image data, the breakthrough came with deep (i.e., many layers) CNNs automatically learned through the back-propagation algorithm. Their layers form a feature hierarchy of increasingly complex features detected by the single layers through convolving the input of shared-weights kernel neurons, which themselves can be simulated by perceptrons. Interspersing pooling layers with spatial strides allows to learn texture or in general global features of the input images. See Schmidhuber (14) for an accurate explanation of these principles.

One of the first approaches modeling stroke tissue outcome with a deep learning CNN has been presented by Stier et al. (15). They trained a 2D-patch-based architecture with respect to the Tmax feature from MR perfusion observed for acute ischemic stroke patients and a follow-up segmentation on FLAIR about 4 days later. The patch-based method clearly outperformed voxel-based approaches. As with other typical black box-like deep learning models, there are no further hyperparameters or constraints that can be set to control the prediction, e.g., for estimating the effect of time.

There are two major challenges in deep learning for stroke analysis tasks with regard to the data: First, there exists a general lack of accessible medical (ground truth) data and, second, the data is of irregular temporal nature. That makes it difficult to apply regular sequence models, such as Markov chains, recurrent neural networks, or the Long-Short-Term-Memory of Hochreiter and Schmidhuber (16). The data is temporally scattered: The points in time tOnset, tImaging, and tTreatment are sampled as patients rush into the hospital's emergency room and cannot be collected in a regular manner.

#### 2.2.1. U-Net Architecture

The U-Net architecture of Ronneberger et al. (17) has been successfully and widely used for biomedical applications by producing semantic segmentations through a fully-convolutional CNN (18) that additionally incorporates skip connections between the context encoding and the refining decoding path for each scale level. The encoding-decoding pattern has established well especially for fully-convolutional networks and is also known from auto-encoders, as used in our proposed method of this paper. Considering different scales is usually a good approach to capture context and details, and this works already well with just two pathways as in the DeepMedic architecture of Kamnitsas et al. (19) who won the sub-acute ischemic stroke lesion segmentation task of the first ISLES challenge in 2015 (20).

At the 2017 edition of the ISLES challenge (21) we presented a robust network on perfusion image data to predict an average lesion outcome and ranked second overall for the binary segmentation output. Many of the top-ranked methods exploited a U-Net architecture, such as the challenge winner (22) who used a 3D U-Net within an ensemble along with other networks and focused on its hyperparameter optimization. In our 2D network instead, we added further skip connections within the encoding path to enhance sensitivity in particular for the difficult smaller lesions in comparison to a standard U-Net (23).

We did not observe advantages when providing clinical variables (e.g., disease scores, time points) as constant input features along with the perfusion images to predict the followup lesion, although they are known to be good predictors for the outcome. In fact, the 2D U-Net performance on the ISLES data could not benefit from the additional information and so we only used spatial perfusion maps to train on. The visual results suggested that rather the robust image features for detecting some highly probable necrotic stroke tissue were learned (cf. also our experiments later in this paper: **Table 4**, **Figure 9**). This makes such a network suitable for segmenting present perfusion lesions, but requires a new strategy to make use of clinical variables for predicting follow-up lesions. It remains an open question in literature, how to ideally incorporate clinical data in a U-Net-only architecture.

#### 2.2.2. Biomedical Shape Regularization

Auto-encoders (AE) are one of several unsupervised methods to learn meaningful features from a data representation by typically encoding the input data x ∈ into a lower-dimensional representation (bottleneck) and decoding this representation to get the reconstruction z ∈ of the input x (Equation 2). This can be achieved through classical fully-connected layers or also by shared-weights convolutional layers for image data. In a convolutional auto-encoder (CAE) as introduced by Masci et al. (8) encoder E(x) usually consists of a typical convolutional feature hierarchy (akin to CNNs) that results in a discriminative latent code y ∈ M, which could be a feature vector or map. Decoder D(y) computes a reconstruction z back in input space . During training, the weights of both are optimized such that a loss L(x, z), e.g., the mean squared error 1 n P n (x − z) 2 for n training samples, is minimized. If used with volumetric segmentations, one can learn shape embeddings on a low-dimensional manifold M (**Figure 1**) with its dimensions representing some main modes of the shapes by optimizing the CAE:

$$z = (D \circ E)(\mathfrak{x}).\tag{2}$$

The principle of shape-constrained segmentation learning was proposed by Ravishankar et al. (24), whose cascaded architecture includes a U-Net and a CAE for shape regularization. While the U-Net follows the same encoder-decoder principle like the convolutional auto-encoder, it does not learn local geometry and shape but produces rather noisy predictions through its skip connections that skip its inner bottleneck. The authors combine both sub-tasks of segmentation and reconstruction in an overall loss to utilize the anatomically regularizing bottleneck of the

auto-encoder for completing noisy kidney segmentations, which improves the segmentations by about 5% compared to U-Net only.

l ).

Oktay et al. (25) presented an anatomically constrained neural network (ACNN) approach to also incorporate shape constraints of anatomical labels as prior knowledge. Their generic training scheme can be applied to various image analysis tasks and was documented for image segmentation and super-resolution. By using a CAE that is trained on ground truth shapes, they constrain the predicted image segmentations to lie close to the learned latent representation of the ground truth. In the end the decoder produces an anatomically constrained reconstruction of the segmentation from the learned shape space, because the segmentation has indeed been forced during training to lie close to the anatomy shape ground truth.

#### 2.3. Our Contribution

In this paper – based on the robust results that U-Nets achieve on perfusion imaging data and the shape-constrained network idea of Ravishankar et al. (24)—we present a novel methodology that:

## 3. METHODS

Our idea is to estimate a time-to-treatment-dependent tissue outcome based on CBV and TTD perfusion images. We hypothesize that the minimum and maximum extents of the potential final stroke lesion can be approximated by delineating the core and penumbra area on the perfusion maps. This includes the following assumptions:

It should be noted that training a model normalized to the acute stroke phase time range of 24 h is possible (cf. results of our evaluation later in this paper: **Table 3**, **Figure 8**) and recommended, if enough follow-up lesion data is available to sample roughly the entire space between core (0 h) and core+penumbra (24 h) representations to avoid areas of uncertainty.

#### 3.1. Architecture

The method consists of a two-phase neural network that combines three main components for automatic shape-constrained follow-up lesion prediction (**Figure 2**):

First, the perfusion images ICBV and ITTD are processed by a U-Net U to compute the segmentation estimates Sˆ c and Sˆ cp. Second, the encoder E of the CAE transforms each segmentation into a low-dimensional shape embedding yˆc and yˆcp of a shape space that must be learned beforehand. Linear interpolation (yˆi , Equation 3) between the latent core and core+penumbra codes yˆc and yˆcp is conducted according to the expected tImaging→Treatment time, which must be normalized by the remaining time to reach 10 h after onset (corresponding to the total core+penumbra).

$$\hat{\jmath}\_i = \hat{\jmath}\_\varepsilon + \eta(\hat{\jmath}\_{cp} - \hat{\jmath}\_\varepsilon), \quad \eta = \frac{t\_{\text{Imaging} \to \text{Treatment}}}{10 - t\_{\text{Onset} \to \text{Imaging}}} \tag{3}$$

This linear interpolation in shape space is crucial, as it corresponds to a non-linear interpolation of the reconstructed shapes on the manifold (**Figure 1**). The decoder D of the CAE is required to compute that reconstruction of the interpolated code yˆi in the voxel space with a segmentation Sˆ l for the final lesion as result (**Figure 3** illustrates a binarized 3D segmentation).

#### 3.1.1. Cascading Networks for Prediction

The construction of a two-phase network targets different subtasks that constrain the learning of the high-dimensional and complex overall task of follow-up lesion prediction. When discriminating only the final lesion binary label from background through high-dimensional multivariate input, it makes it hard for a machine learning algorithm to properly generalize the relation between input and outcome. Instead of directly estimating a follow-up segmentation from the input data, we guide the first sub-network (U-Net) to segment core and penumbra correctly, which is—as explained before—of major importance for predicting final lesion tissue outcome. Once this data is provided, the second sub-network (CAE) can learn the most salient shape features on a rather simplified representation with respect to the task (shape probability maps vs. different physiological CTP parameters) along with clinical data to estimate the follow-up lesion.

#### 3.1.2. U-Net

Instead of taking the 2D U-Net that we have used before at the ISLES 2017 challenge, we employ a smaller standard 3D U-Net U to reduce computational and memory demand while it can cope well with the three-dimensional nature of stroke volume data. Furthermore, instead of forwarding the full 128 × 128 × 28 input CTP images the U-Net is trained on randomly positioned cubic patches of size 64 × 64 × 28 (with additional padding of 20 voxels in each direction) and thus needs to forward 4 patches for segmentation of one single case. It receives patches from ICBV and ITTD as input and estimates Sˆ c , Sˆ cp = U(ICBV, ITTD). The U-Net is build of doubleconvolutional blocks as known from Ronneberger et al. (17), while each 3 × 3 × 3 convolution is preceded by a batch normalization layer of Ioffe and Szegedy (26) for data whitening. The blocks are spread over three resolution levels (two Max-Poolings) with 16, 32, and 64 channels, respectively. This sums up to a total of about 355.000 network parameters.

#### 3.1.3. Convolutional Auto-Encoder

Focusing on the subsequent CAE, we had to ensure some minimum number of layers (**Figure 4**) to detect the most salient and descriptive abstract stroke shape features while being wellregularized in order to reconstruct a good general estimate of the shape. This requires a bottleneck layer between E and D to produce low-dimensional latent codes that must have a limited but sufficient dimensionality. Consequently, only the main modes should be represented in the code with their major variances describing different stroke shapes without noise or overfitting of training samples. Regarding overfitting, our linear approach for the interpolation on the low-dimensional codes is in principle also robust against noise in time.

The input of the CAE is forwarded akin to the U-Net through double-convolutional blocks including batch normalization

layers. Instead of two Max-Pooling operations, three additional 2-stride convolutional layers intersperse those blocks, while the final block is a single 3 × 3 × 3 convolution that convolves the feature map to the 10 × 10 × 1 bottleneck size. Like for the U-Net, 3 × 3 × 3 filter kernels throughout all convolutional layers of the CAE, as introduced by Simonyan and Zisserman (27), are exclusively used to decrease complexity compared to networks with bigger kernels while the same receptive field sizes can be reached by just stacking more layers of the smaller kernels. This results in more non-linearities and usually generalizes better.

The decoder is a mirrored encoder, with deconvolutional layers replacing the 2-stride convolutions. With this architecture, a single 128×128×28 shape segmentation image can be encoded into a low-dimensional code representation, and then decoded back to a segmentation image. Note that although the U-Net is fed with both CBV and TTD images as separate channels at once, the CAE's encoder has to forward their segmentations independently to get two separate latent codes, which will be then interpolated and forwarded in one step through the decoder to get the final lesion prediction.

Since there are several global scalar predictors (time, age, sex, clinical scores) that might be required to model the space properly, we also tried to map the shape segmentation input of the CAE to a vector representation (instead of spatial feature maps) in the bottleneck by using a fully-connected or convolutional layer with appropriate kernel size. This would have allowed us to add an arbitrary number of scalar values directly as additional dimensions to the latent code in the bottleneck and to easily quantify the dimensionality of such shape representations in the latent space. However, regularization is difficult – even when using dropout (28)—and the reconstructions are less accurate than with latent spatial feature map representations. Therefore, the CAE remains convolutional-only and no clinical variables other than the combination η of both time predictors for the interpolation are used as per definition in Equation (3).

#### 3.2. Training

The training is conducted in three consecutive steps (see **Figure 5** for illustration of steps 2 and 3) that are characterized by different objectives formulated in their corresponding losses:

The U-Net is initially trained beforehand using the stochastic gradient variant ADAM of Kingma and Ba (29) for optimization. The ground truth segmentations Sc and Scp of core and core+penumbra are used to penalize their predictions with a bigger loss for less overlap in the SoftDice measure, which is defined for all voxel positions i in a segmentation A with ground truth B and a small constant ǫ as follows:

$$\text{SoftDice}(A, B) = \frac{2 \cdot \sum\_{i} (A\_i B\_i) + \epsilon}{\sum\_{i} (A\_i A\_i + B\_i B\_i) + \epsilon} \tag{4}$$

#### 3.2.1. Shape Space Learning

First, a low-dimensional shape space is learned that embeds the ground truth shapes of core and core+penumbra segmentations, Sc and Scp. This is enforced by loss LShape in Equation (5), that consists of three parts: (1) Reconstruction Rc , Rcp, and Rl of core, core+penumbra and final lesion shape, (2) the property of the reconstructed core/lesion volume to still be a subset of the core+penumbra volume, and (3) a L1 loss for keeping the latent code yl of the lesion shape close to the linear interpolation yi of the core and core+penumbra codes:

$$L\_{\text{Shap\epsilon}} = \sum\_{s \in \{c, cp\}} L\_{\text{SoftDice}}(R\_s, S\_s) + \sum\_{s \in \{c, l\}} L\_{\text{mono}}(R\_s, R\_{cp}) + \alpha L1(y\_l, y\_i) \tag{5}$$

For the first 25 epochs α = 0 holds, otherwise α = 1. We observed that this helps the non-convex optimization function to first learn to reconstruct the shape into the correct brain hemisphere and was found to robustly prevent the network from getting trapped in implausible minima at the beginning of the training. Since SoftDice ∈ [0, 1] with 0 indicating non-overlap and 1 for full overlap, we need to define

FIGURE 4 | Representational complexity of the CAE and reduced dimensionality in the bottleneck. The numbers given top to bottom indicate: Quadratic size of the first two spatial dimensions, size in the third spatial dimension (axial), and number of feature channels. There is one input channel for the segmentation image and one output channel for its reconstruction.

LSoftDice(S, Sˆ) = 1 − SoftDice(S, Sˆ). To force the CAE to learn that the interpolation is only growing when interpolating along the (time) trajectory from core until reaching the core+penumbra segmentation at maximum, we define a constraint Lmono for all voxel positions i in two segmentation images A, B so that the reconstructions of the core segmentation and the intermediate lesion interpolation are monotone increasing to the total core+penumbra segmentation:

$$L\_{mnoo}(A,B) = \sum\_{i} \max(A\_i - B\_i, 0) \tag{6}$$

#### 3.2.2. Noisy Shape Interpolation

In the second training phase, the encoder and decoder weights from the shape space learning phase before will be fixed. A second encoder is then learned for the U-Net predictions Sˆ c and Sˆ cp of core and core+penumbra to compute latent representations yˆc and yˆcp that are located close to the shape embeddings yc and ycp of the corresponding ground truth segmentations in terms of L1 norm. LPrediction in Equation (7) further enforces the monotone properties for the reconstructed segmentations and the main goal of high overlap for the prediction Sˆ l decoded from the interpolated code yˆi to the actual follow-up ground truth Sl :

$$L\_{\text{Prediction}} = L\_{\text{SoftDice}}(\hat{\mathbf{S}}\_l, \mathbf{S}\_l) + \sum\_{s \in \{c, l\}} L\_{\text{mono}}(\hat{\mathbf{S}}\_s, \hat{\mathbf{S}}\_{cp}) + \sum\_{s \in \{c, l, i\}} L1(\hat{y}\_s, y\_s) \tag{7}$$

This way, the decoder D of the first phase ideally decodes an approximate representation from the shape space so that the reconstruction of the core and core+penumbra estimates should be close to the ground truth reconstruction of core and core+penumbra, (D ◦ E2)(Sˆ c) ≈ (D ◦ E1)(Sc) and (D ◦ E2)(Sˆ cp) ≈ (D ◦ E1)(Scp). Moreover, the main goal remains to achieve an interpolation as close as possible to the true lesion segmentation:

$$D\left(\eta \cdot E\_{\mathbf{2}}(\hat{\mathbf{S}}\_{\mathbf{C}^{\mathbf{p}}}) + (1 - \eta) \cdot E\_{\mathbf{2}}(\hat{\mathbf{S}}\_{\mathbf{c}})\right) \approx (D \odot E\_{\mathbf{1}})(\mathbf{S}\_{\mathbf{l}}) \tag{8}$$

The individual loss terms have not been weighted by further manual parameters as we found no benefit for other than the uniformly weighted individual loss parts for both LShape and LPrediction. Further, we tried to apply both losses beforehand in an alternating and joint manner but this did not let the optimizer find proper minima; in particular, learning the exclusive occurrence of a stroke on either hemisphere could not be learned, which is what the network basically learns during the first epochs of training.

### 4. EXPERIMENTS

We run a 5-fold test on a 29 subjects dataset, because the time demand for a full cross-validation was too high. In order to test each of the 5 folds, we thus had to train 5 models with the remaining 4 folds. Four of the folds consist of 6, and 1 fold consists of 5 cases. About one fourth of the training samples were used as validation set, so that a training set consists of 17 or 16 cases and is validated in each epoch on 6 other cases (disjoint with the test fold).

We prevent overfitting of the model by training until the validation loss converges and choose the model with the lowest validation loss. Since the model is eventually tested on a different fold of patients not used for the training and validation, we also avoid that the evaluation results could be tuned on the validation loss optimum. Due to the huge number of parameters in our 3D sub-networks, the memory demand for gradient computation increases rapidly, so a batch size of 4 had to be chosen to fit the training data into 11 GB of GPU memory.

TABLE 2 | Characteristics for subjects with manually segmented core, penumbra and follow-up lesion of the retrospectively collected data.

### 4.1. Data

We used a dataset of 29 subjects from the Neuroradiology department at the University Hospital Schleswig-Holstein formerly collected for the TRAVESTROKE project for which one rater had created manual segmentations on the CT perfusion (CTP) modalities CBV (Cerebral Blood Volume) and TTD (Time-To-Drain) at the time of admission for core and core+penumbra, as well as a lesion segmentation on follow-up CT after treatment. The data was acquired with a Siemens Somatom Definition AS 40 (Siemens Healthcare GmbH Forchheim, Germany) and the raw data was deconvolved using the vendor algorithm to get CT perfusion parameters such as CBV or TTD. All patients of the dataset had been treated successfully with thrombectomy (TICI score 2b or 3). See baseline characteristics of the subjects included in the evaluation in **Table 2**.

The dataset was pre-processed with FSL-FLIRT (30) for affine registration to correct tilted heads and transform them into common space. A downsampling was applied so that input size of the CBV and TTD maps was 128 × 128 × 28 voxels for reducing the computational demand. Additional clinical data for each subject was given. For the evaluation we only used the two durations tOnset→Imaging and tImaging→Treatment, which were normalized according to Equation (3). All shape segmentations have been elastically deformed in each epoch to augment the limited available training data for learning generic features of the CAE such that the representations in shape space can be robustly reconstructed.

### 4.2. Comparison

In section 2.2.1, we referred to the ISLES 2017 task of predicting a follow-up lesion based on MR perfusion data. We participated in this challenge with a single U-Net directly predicting the final follow-up lesion, as many of the teams were using unconstrained CNNs [as presented in (21)]. Since this does not lead to accurate predictions of a progressed stroke when facing acute image data, we compare our proposed method with this simple U-Net approach with and without clinical time

While the first row shows an oracle prediction overlap for the theoretically best-fit per case interpolation within our model (could be before or after the true time-to-treatment), the second row lists the result for our proposed approach using the actual time-to-treatment ground truth to predict the follow-up lesion from ground truth. The third row indicates that the effect of choosing a different normalization value from within the time range of the acute stroke phase can rather be neglected.

TABLE 4 | Experimental results from 5-fold test data based on the U-Net's Sˆ c, Sˆ cp, and lesion estimates: The average values for the Dice overlap with the ground truth segmentations Sc, Scp, and Sl are presented.

The baseline U-Net-only approach has the same architecture as used in combination with the CAE above, but instead of predicting core and penumbra (2 channel output), it directly predicts the follow-up lesion (1 channel output) as described in section 4.2. We input CBV and TTD maps only (2in), or additionally both tOnset→Imaging and tImaging→Treatment as constant image channels (4in), whereas the latter performed worse coincident with our previous experience (see section 2.2.1 and Figure 9). The highest overlaps for core, core+penumbra, and lesion are highlighted in bold.

points as input. Unfortunately, the ISLES dataset consists of MR perfusion data without appropriate core and penumbra segmentations, so we cannot directly compare on the same dataset.

Our sub-task of linearly interpolating along the trajectory between core and core+penumbra requires representations of such in a suitable non-linear shape space that has to be learned before. In order to show the advantage of conducting this in such a shape space, we compare with the naïve way of linearly interpolating the shape segmentations in . This can be simply computed with the same η as defined before in Equation (3):

$$
\hat{S}\_l = \hat{S}\_c + \eta (\hat{S}\_{cp} - \hat{S}\_c) \tag{9}
$$

#### 5. RESULTS

The reconstruction results (**Table 3**) demonstrate the capability of our learned model to make time-dependent predictions and present the advantages of our CAE approach using core

FIGURE 7 | Visual comparison of linear interpolation in image and shape space for an axial slice of a single case from the 5-fold test (Ground truth denoted as GT). Compare the predictions at 2h with the actual follow-up at 1.7h: With linear interpolation of the segmentations, the core area Sˆ c remains unchanged while the huge penumbra area is faded in. Note that for a normalization of 24h, the fading over time would progress even slower. With the CAE and its interpolation in shape space, the shape grows non-linearly, first locally, and then quickly into the surrounding tissue-at-risk areas segmented by the U-Net in Sˆ cp.

FIGURE 8 | Three example cases (as in Figure 9) from the 5-fold test data with their results for the reconstruction. The top row shows a case with fast admission and treatment times, where the lesion appears about the same size as the initial core. The middle and bottom row are interpolated at η located one third on the trajectory from core to core+penumbra, while one case is an early and the other is a late admission; However, the case with bigger core progresses slowly with respect to the non-linear reconstruction compared to the smaller core case at the bottom. Reference overlays for manual segmentations Sc (red), Scp (orange), and Sl (green) are shown as outlines. Note that models trained with a normalization of 10 or 24 h output similar predictions.