Stroke onset time determination using MRI relaxation times without non-ischaemic reference in a rat stroke model

Abstract

Background:

Objective timing of stroke in emergency departments is expected to improve patient stratification. Magnetic resonance imaging (MRI) relaxations times, ${\mathrm{T}}_2$ and ${\mathrm{T}}_{1\mathit{\rho }}$ , in abnormal diffusion delineated ischaemic tissue were used as proxies of stroke time in a rat model.

Methods:

Both ‘non-ischaemic reference’-dependent and -independent estimators were generated. Apparent diffusion coefficient (ADC), ${\mathrm{T}}_2$ and ${\mathrm{T}}_{1\mathit{\rho }}$ , were sequentially quantified for up to 6 hours of stroke in rats ( $\mathit{n}=8$ ) at 4.7T. The ischaemic lesion was identified as a contiguous collection of voxels with low ADC. ${\mathrm{T}}_2$ and ${\mathrm{T}}_{1\mathit{\rho }}$ in the ischaemic lesion and in the contralateral non-ischaemic brain tissue were determined. Differences in mean MRI relaxation times between ischaemic and non-ischaemic volumes were used to create reference-dependent estimator. For the reference-independent procedure, only the parameters associated with log-logistic fits to the ${\mathrm{T}}_2$ and ${\mathrm{T}}_{1\mathit{\rho }}$ distributions within the ADC-delineated lesions were used for the onset time estimation.

Result:

The reference-independent estimators from ${\mathrm{T}}_2$ and ${\mathrm{T}}_{1\mathit{\rho }}$ data provided stroke onset time with precisions of $\pm 32$ and $\pm 27$ minutes, respectively. The reference-dependent estimators yielded respective precisions of $\pm 47$ and $\pm 54$ minutes.

Conclusions:

A ‘non-ischaemic anatomical reference’-independent estimator for stroke onset time from relaxometric MRI data is shown to yield greater timing precision than previously obtained through reference-dependent procedures.

1. Introduction

The current pharmacological treatments of an ischaemic stroke, including recombinant plasminogen activator, can only be administered within strict time limits [6]. This is due to two key reasons, as (i) irreversible damage to brain parenchyma advances with time [26] and (ii) the benefits of pharmacotherapy are outweighed by adverse effects [4]. Estimates of the number of neural cells and functionality lost per minute of ischaemia [21] stress the earliest possible administration of thrombolytic and/or anti-ischaemic therapy. Thus, stroke duration is seen as a crucial proxy of both brain state and a predictor of functional recovery, and consequently will play a pivotal role in the objective determination of onset time in the clinical settings.

A growing body of evidence indicates magnetic resonance imaging (MRI) protocols provide clinically acceptable timing of stroke onset [9,20,27,29]. The difference between MRI relaxation times ${\mathrm{T}}_2$ and ${\mathrm{T}}_{1\mathit{\rho }}$ in ischaemic brain tissue, as defined by its low apparent diffusion coefficient (ADC), and contra-lateral non-ischaemic tissue has been observed to increase linearly with time, thus, providing the means for estimating the onset of ischaemic stroke [9,25]. For example, an application of the ${\mathrm{T}}_2$ difference approach for stroke timing has yielded promising results in patients [27], and a similar approach using ${\mathrm{T}}_2$ -FLAIR as a proxy of stroke duration [29] has been trialled in patients with variable success [3]. Of these, the latter apparently has an acceptable specificity for clinical exploitation, however, the sensitivity remains low [29] when compared with the use of quantitative ${\mathrm{T}}_2$ relaxation times [27].

The key difficulty in the application of MRI relaxation times in stroke timing clinically is that it is difficult to know the baseline value at each voxel of the ischaemic lesion before the onset of ischaemia, such that the deviation from baseline due to ischaemia cannot be calculated with certainty. For this reason, some authors have sought to select a set of non-ischaemic voxels as a ‘reference’ region, representative of the lesion prior to the onset of ischaemia [27]. However, this approach is confounded by the inherently varying MRI relaxation times across the brain [30], such that the ‘reference’ voxels may not faithfully represent those of the lesion prior to the onset of ischaemia. Given that the distribution of relaxation times in the human brain is broader than the change caused by ischaemia, the challenge of selecting an appropriate set of reference voxels, and the consequences of failing to do so, have presented a major barrier.

Motivated by these factors, we present here a mathematical procedure for ‘non-ischaemic reference’-independent stroke time estimation. Relaxometric MRI data from permanent focal ischaemia in rats was used to assess the capability of the mathematical model for timing of stroke onset.

2. Methods

2.1. Animal model

Male Wistar rats ( $\mathit{N}=8$ , approximately 300 g) underwent permanent middle cerebral artery occlusion (MCAO) to induce focal stroke [17]. Animals were euthanized by intravenous injection of saturated KCI at the end of the experiment. Animal procedures were conducted according to European Community Council Directives 86/609/EEC guidelines and approved by the Animal Care and Use Committee of the University of Eastern Finland.

2.2. Image acquisition

MRI data were acquired according to the methods described by Jokivarsi et al. [9]. The trace of the diffusion tensor ( $\mathrm{ADC}=1\mathrm{/}3·Trace[\mathrm{D}]$ ), ${\mathrm{T}}_2$ and ${\mathrm{T}}_{1\mathit{\rho }}$ MR images from a single 2 mm slice were obtained at 30–60 min intervals for 6 h after induction of focal brain ischaemia using a horizontal Varian Unity^Inova MRI system operating at 4.7T with a volume coil transmit/quadrature half-volume receive set-up (Rapid Biomedical GmbH, Rimpar, Germany). The ADC was determined from spin echo images with 4 bipolar gradients along each axis with 4 b-values ranging from 0 to 1,370 s mm² (field of view $2.56\times 2.56{\text{cm}}^2$ , $64\times 128$ data points, echo time 55 ms). Fast spin-echo (FSE) readout data ( $64\times 128$ pixels, echo spacing 10 ms, FOV $2.56\times 2.56{\text{cm}}^2$ ) was used for ${\mathrm{T}}_{1\mathit{\rho }}$ and ${\mathrm{T}}_2$ MRI. The on-resonance spin-lock ${\mathrm{T}}_{1\mathit{\rho }}$ MRI was acquired with a continuous wave ${\mathrm{T}}_{1\mathit{\rho }}$ approach. A contrast formation block of AHP-SL-AHP (AHP = adiabatic half passage, SL = spin-lock) segment was added before readout. Adiabatic spin-lock pulses ranging from 8 to 64 ms were used [SL amplitude (B_1,SL) 0.4 G, time to repetition (TR) 2.5 s and time to echo (TE) 6 ms]. ${\mathrm{T}}_2$ images were collected using the above FSE readout in front of a preparation block consisting of adiabatic pulses (AHP-TE/4-AFP-TE/2-AFP-TE/4-reverse; five TEs of 10, 20, 40, 60 and 80 ms). The adiabatic refocussing pulses in a single slice mode was used to minimise contributions of non-spin echo coherences [23]. Total acquisition time for $\mathrm{ADC}$ , ${\mathrm{T}}_2$ and ${\mathrm{T}}_{1\mathit{\rho }}$ was 24 min.

2.3. Image processing

All image post-processing was performed using in-house designed software written in the MATLAB programming environment (MathWorks, Natick, Massachusetts, USA).

${\mathrm{T}}_2$ and ${\mathrm{T}}_{1\mathit{\rho }}$ weighted MR images were filtered with a squared Hamming window function in the phase and frequency encoded dimensions of k-space prior to computing parametric maps. ${\mathrm{T}}_2$ and ${\mathrm{T}}_{1\mathit{\rho }}$ relaxation maps were computed from a voxel-wise monoexponential fit, where the data were first transformed to a logarithmic space then fitted via linear least squares. ${\mathrm{T}}_2$ and ${\mathrm{T}}_{1\mathit{\rho }}$ weighted images were generated by summing all echo and spin-lock time images, respectively.

Lesion detection was performed on reciprocal ADC images where regions of hyperintensity are suggestive of an ischaemic assault. More formally, we have defined a lesion as the largest connected cluster of voxels whose values all exceed one median absolute deviation above the whole brain median value [13]. The output from the detection procedure was a binary mask, where ischaemic voxels were flagged with the value of one and all others set to zero.

2.4. Statistical models for estimating onset time

At present, the state-of-the-art estimators for time from stroke onset are based upon measuring $\mathrm{\Delta }{\mathrm{T}}_2$ or $\mathrm{\Delta }{\mathrm{T}}_{1\mathit{\rho }}$ , i.e. the difference between the mean ${\mathrm{T}}_2$ or ${\mathrm{T}}_{1\mathit{\rho }}$ within the lesion and some user-defined reference region, for example, the contra-laterally mirrored lesion [20,25,27]. The mean value is typically estimated via the average of the voxel values within a given region. These averages are then used to train a linear regression model that is the main tool through which estimation is performed.

One of the primary disadvantages associated with these estimators is the arbitrary choice of the reference region, which if not chosen appropriately, in terms of the anatomical and/or tissue type match with respect to the ischaemic lesion, can dramatically affect their reliability. Furthermore, their robustness is governed by the ability to accurately estimate the mean ${\mathrm{T}}_2$ or ${\mathrm{T}}_{1\mathit{\rho }}$ values, which is challenging given that the number of ischaemic voxels can vary considerably between strokes and with time in stroked tissue. Also, voxel values can exhibit heavy-tailed empirical distributions for which the mean value is not well defined, making the average estimator potentially unreliable.

Motivated by the disadvantages mentioned above, we have explored the idea of designing robust reference-independent estimators that exploit information from the full empirical distribution of the ${\mathrm{T}}_2$ or ${\mathrm{T}}_{1\mathit{\rho }}$ data in ischaemic voxels only. For improved robustness with large as well as small regions, we underpin the linear regression estimators with a parametric statistical model of the empirical distribution that provides a low-dimensional structure with favourable regularity properties.

2.5. Distributional models for MR parameters

Appropriate parametric statistical models for ${\mathrm{T}}_2$ and ${\mathrm{T}}_{1\mathit{\rho }}$ voxel values within ischaemic lesions were identified by using the MATLAB statistical distribution fitting toolbox. In particular, histograms of ${\mathrm{T}}_2$ and ${\mathrm{T}}_{1\mathit{\rho }}$ within the lesion were computed, then area-normalised to represent probability densities. Individual (specified below) distributions were fitted to the data by maximum likelihood estimation and compared with their empirical distributions. The model that on average scored the best Akaike information criteria (AIC) value was then selected to design the regression method.

2.6. Training and testing of estimators

As with the state-of-the-art, our estimators were chosen to take the form of a multiple linear regression. Here, it is important to note that the assumption of linearity applies only to the coefficients of the estimators and does not necessarily imply that MR parameters vary linearly in time.

To mitigate the possibility of over-fitting we used a leave-one-out cross-validation approach to split the data used to train and to test the estimators. Precisely, our dataset of 8 rats (total of 45 sets of relaxometric data) was iteratively separated into a training set consisting of the scans of 7 rats, and a testing set consisting of the scans of the remaining rat. For each ${\mathrm{T}}_2$ scan in the training set, the parameters associated with the statistical model of ischaemic ${\mathrm{T}}_2$ were estimated by maximum likelihood and then used collectively to train a ${\mathrm{T}}_2$ polynomial regression estimator (this procedure was iterated over all rats, yielding a total of 8 ${\mathrm{T}}_2$ and 8 ${\mathrm{T}}_{1\mathit{\rho }}$ reference-independent estimators). A similar approach was used for signal intensities in ${\mathrm{T}}_2$ and ${\mathrm{T}}_{1\mathit{\rho }}$ weighted images. Finally, the performance of the regression estimates was assessed by measuring the root mean-squared-error (RMSE) with respect to the true onset times.

To serve as a comparison, we also computed linear regressions based on $\mathrm{\Delta }{\mathrm{T}}_2$ and $\mathrm{\Delta }{\mathrm{T}}_{1\mathit{\rho }}$ predictors, where the reference region was defined as the contra-lateral mirrored lesion as described previously [25].

3. Results

3.1. Lesion detection

ADC images, including lesion mask overlays, from a typical rat brain acquired over multiple time-points are given in Figure 1(A). As indicated in the images, the lesion masks form a single contiguous region and their volumes generally increase from time of stroke onset reaching $38\pm 8$ % of the MRI slice area by 6 hours. ADC values of voxels in the lesion were $49.7\pm 13.4$ % from those determined in the contralateral non-ischaemic brain parenchyma confirming presence of severe ischaemia [7].

Fig. 1.

(A) ADC images (green) with lesion masks (pink), (B) ${\mathrm{T}}_2$ and (C) ${\mathrm{T}}_{1\mathit{\rho }}$ maps from a typical rat across multiple time points (separated by approximately 1-hour intervals). (D) Distribution of ${\mathrm{T}}_2$ within ADC lesion and (E) distribution of ${\mathrm{T}}_{1\mathit{\rho }}$ within ADC lesion across multiple time-points.

3.2. Statistical models for MR parameters

Histograms are given for ${\mathrm{T}}_2$ (Figure 1(D)) and ${\mathrm{T}}_{1\mathit{\rho }}$ (Figure 1(E)) within the identified lesions for a typical rat (see Figure 1(A) for ADC, 1(B) for ${\mathrm{T}}_2$ and 1(C) for ${\mathrm{T}}_{1\mathit{\rho }}$ images). The histograms show that there was a right-shift in both ${\mathrm{T}}_2$ and ${\mathrm{T}}_{1\mathit{\rho }}$ distributions within the lesion as time from stroke onset increased. It is precisely this information that is exploited in the current $\mathrm{\Delta }{\mathrm{T}}_2$ and $\mathrm{\Delta }{\mathrm{T}}_{1\mathit{\rho }}$ estimators [9]. In addition to the value shift, the histograms also reveal that the general spread of the ${\mathrm{T}}_2$ and ${\mathrm{T}}_{1\mathit{\rho }}$ distributions increases with time from stroke onset. A similar investigation focusing now on the contra-lateral mirror reference was also performed and found that both ${\mathrm{T}}_2$ and ${\mathrm{T}}_{1\mathit{\rho }}$ distributions remain approximately stationary across time-points (not shown).

As mentioned earlier, in order to identify a low-dimensional representation of the empirical distribution, several statistical distributions were fitted to the data. Namely, the normal, lognormal, log-logistic [5] and scaled Burr distributions [1]. Of these models the 2-parameter log-logistic distribution, i.e. $$\begin{array}{c}\mathit{f}\mathrm{(}\mathit{x}\mid \mathit{\mu }\mathrm{,}\mathit{\sigma }\mathrm{)}={\displaystyle \frac{1}{\mathit{\sigma }}}{\displaystyle \frac{1}{\mathit{x}}}{\displaystyle \frac{exp\mathrm{(}\mathit{z}\mathrm{)}}{{[1+exp\mathrm{(}\mathit{z}\mathrm{)}]}^2}}\mathrm{,}\text{where}\mathit{z}={\displaystyle \frac{log\mathrm{(}\mathit{x}\mathrm{)}-\mathit{\mu }}{\mathit{\sigma }}}\mathrm{,}\end{array}$$ was found to be the most suitable for both ${\mathrm{T}}_2$ and ${\mathrm{T}}_{1\mathit{\rho }}$ values within the lesion, as determined by obtaining, on average, the best AIC value (Figure 1(D) and (E) display the distribution fits in red). Here, the log-logistic parameters $\mathrm{(}\mathit{\mu }\mathrm{,}\mathit{\sigma }\mathrm{)}$ respectively denote the log-median value and distributional shape for a given MR parameter within the lesion.

The log-logistic model was also found to be suitable for describing the distributions of ${\mathrm{T}}_2$ and ${\mathrm{T}}_{1\mathit{\rho }}$ signal intensities within the lesion. However, very little difference was observed between time-points, i.e. little to no drift or spread (Supplementary Information, Figure S1).

3.3. Regression estimators

The anatomical reference-independent regression estimators were designed using (functions of) the parameters associated with the log-logistic fits, i.e. the μ and σ parameters. For the ${\mathrm{T}}_{1\mathit{\rho }}$ -based estimator, a 3-parameter regression (including one constant term) was found to be most accurate: $$\begin{array}{c}{\mathit{t}}_{\mathrm{estimated}}={10}^3·[\mathrm{(}1.3249\pm 0.1668\mathrm{)}{\mathit{\mu }}_{{\mathrm{T}}_{1\mathit{\rho }}}+\mathrm{(}2.6680\pm 0.9232\mathrm{)}{\mathit{\sigma }}_{{\mathrm{T}}_{1\mathit{\rho }}}-\mathrm{(}5.8129\pm 0.7211\mathrm{)}]\end{array}$$ ( ${\mathit{R}}^2=0.9229\mathrm{,}\mathit{p}\ll 0.001$ ), with an associated RMSE of $\pm 27\text{mins}$ .

Fig. 2.

Estimated time vs. actual time of onset for several regression estimators. (A) Reference-dependent $\mathrm{\Delta }{\mathit{T}}_{1\mathit{\rho }}$ (B) Reference-dependent $\mathrm{\Delta }{\mathit{T}}_2$ (C) Reference-independent ${\mathit{T}}_{1\mathit{\rho }}$ multiple linear regression. (D) Reference-independent ${\mathit{T}}_2$ multiple polynomial regression. 95% prediction intervals are given by green curves. Zero residual error is indicated by red curves. Coloured dots correspond to estimates from individual rats. Both axes are in minutes.

For the ${\mathrm{T}}_2$ -based estimator, a 7-parameter estimator based on a third-order polynomial was chosen: $$\begin{array}{ll}\hfill {\mathit{t}}_{\mathrm{estimated}}=& {10}^6·[-\mathrm{(}1.4260\pm 1.6028\mathrm{)}{\mathit{\mu }}_{{\mathrm{T}}_2}-\mathrm{(}0.0426\pm 0.0455\mathrm{)}{\mathit{\sigma }}_{{\mathrm{T}}_2}\\ \hfill & +\mathrm{(}0.3475\pm 0.3929\mathrm{)}{\mathit{\mu }}_{{\mathrm{T}}_2}^2+\mathrm{(}1.0624\pm 1.0877\mathrm{)}{\mathit{\sigma }}_{{\mathrm{T}}_2}^2\\ \hfill & -\mathrm{(}0.0282\pm 0.0321\mathrm{)}{\mathit{\mu }}_{{\mathrm{T}}_2}^3-\mathrm{(}7.9923\pm 8.3537\mathrm{)}{\mathit{\sigma }}_{{\mathrm{T}}_2}^3+\mathrm{(}1.9490\pm 2.1789\mathrm{)}]\end{array}$$ ( ${\mathit{R}}^2=0.8918\mathrm{,}\mathit{p}\ll 0.001$ ), with an associated RMSE of $\pm 32\text{mins}$ .

Figures 2(C) and (D) show the fits for the ${\mathrm{T}}_2$ and ${\mathrm{T}}_{1\mathit{\rho }}$ estimators, respectively, where 95% prediction intervals are given by green curves. These results were notably more accurate than those associated with the 2-parameter regression based on $\mathrm{\Delta }{\mathrm{T}}_{1\mathit{\rho }}$ (Figure 2(A), $\mathit{p}\mathrm{<}0.001$ , pairwise Bartlett test) and $\mathrm{\Delta }{\mathrm{T}}_2$ (Figure 2(B), $\mathit{p}=0.01$ , pairwise Bartlett test), which only achieved RMSEs of $\pm 54\text{mins}$ and $\pm 47\text{mins}$ respectively. It is worth noting that uncertainties here are larger than previously reported [25]. This is primarily due to different uncertainty quantification approaches. Rogers et al. [25] used an analytical uncertainty propagation method that is liberal to approximate precision at a given time point, whereas here we assess precision empirically by cross-validation and report RMSE across all time points. An interval-wise breakdown of the results for all estimators is given in Table 1.

In addition to measuring the precision of reference-dependent estimators, we investigated their robustness through consideration of randomly-shifted reference regions (Supplementary Information, Figure S2 and S3). The results showed that overall precision is sensitive to reference placement, with RMSEs often larger than those reported for the mirrored contra-lateral reference.

Table 1

Interval-wise breakdown of estimated vs. actual onset times. Values are displayed as mean ± SD. 95% confidence intervals for estimated means were computed from $\text{mean}\pm {\mathit{t}}_{0.95\mathrm{,}\mathit{N}}·\mathrm{SEM}$ , where SEM is standard error of mean, and ${\mathit{t}}_{0.95\mathrm{,}\mathit{N}}$ is the Student t-value for the 95% interval and N degrees of freedom

${\mathrm{T}}_2$ reference-independent ( $\mathrm{RMSE}=\pm 32\text{mins}$ )
Timing (mins)		Parameters		Estimation (mins)
Interval	Mean time	${\mathit{\mu }}_{{\mathit{T}}_2}$ (log ms)	${\mathit{\sigma }}_{{\mathit{T}}_2}$ ( $\times {10}^2$ )	${\mathit{t}}_{\mathrm{estimated}}$	95% CI
000–060 ( $\mathit{N}=7$ )	41.3 ± 8.4	$3.99\pm 0.02$	$3.37\pm 0.54$	50.3 ± 14.5	37–63
060–120 ( $\mathit{N}=8$ )	97.0 ± 16.4	$4.05\pm 0.02$	$3.54\pm 0.70$	99.7 ± 22.4	81–118
120–180 ( $\mathit{N}=9$ )	157.8 ± 16.8	$4.09\pm 0.03$	$4.02\pm 0.67$	166.9 ± 31.6	143–191
180–240 ( $\mathit{N}=7$ )	218.4 ± 11.0	$4.12\pm 0.02$	$4.56\pm 0.86$	225.5 ± 30.0	199–252
240–300 ( $\mathit{N}=8$ )	271.9 ± 15.6	$4.14\pm 0.03$	$4.83\pm 0.94$	264.7 ± 39.4	233–297
300–360 ( $\mathit{N}=6$ )	333.8 ± 12.4	$4.17\pm 0.03$	$5.00\pm 0.78$	307.4 ± 25.6	282–333

${\mathrm{T}}_{1\mathit{\rho }}$ reference-independent ( $\mathrm{RMSE}=\pm 27\text{mins}$ )
Timing (mins)		Parameters		Estimation (mins)
Interval	Mean time	${\mathit{\mu }}_{{\mathit{T}}_{1\mathit{\rho }}}$ (log ms)	${\mathit{\sigma }}_{{\mathit{T}}_{1\mathit{\rho }}}$ ( $\times {10}^2$ )	${\mathit{t}}_{\mathrm{estimated}}$	95% CI
000–060 ( $\mathit{N}=8$ )	41.9 ± 10.4	$4.35\pm 0.02$	$3.31\pm 0.49$	37.5 ± 22.6	23–60
060–120 ( $\mathit{N}=8$ )	102.4 ± 11.6	$4.40\pm 0.01$	$3.68\pm 0.77$	115.3 ± 22.1	84–120
120–180 ( $\mathit{N}=8$ )	160.4 ± 11.2	$4.43\pm 0.02$	$4.11\pm 0.77$	172.2 ± 24.2	141–180
180–240 ( $\mathit{N}=8$ )	219.4 ± 12.9	$4.45\pm 0.02$	$4.93\pm 0.89$	218.0 ± 23.5	200–239
240–300 ( $\mathit{N}=7$ )	274.1 ± 11.5	$4.50\pm 0.03$	$4.93\pm 0.87$	275.2 ± 28.6	249–300
300–360 ( $\mathit{N}=6$ )	331.8 ± 12.4	$4.51\pm 0.02$	$5.35\pm 0.56$	305.3 ± 16.6	315–348

$\mathrm{\Delta }{\mathrm{T}}_2$ reference-dependent ( $\mathrm{RMSE}=\pm 47\text{mins}$ )
Timing (mins)		Parameters	Estimation (mins)
Interval	Mean time	$\mathrm{\Delta }{\mathrm{T}}_2$ (ms)	${\mathit{t}}_{\mathrm{estimated}}$	95% CI
000–060 ( $\mathit{N}=7$ )	41.3 ± 8.4	−2.92 ± 1.71	59.0 ± 41.6	22–96
060–120 ( $\mathit{N}=8$ )	97.0 ± 16.4	−0.02 ± 1.49	129.6 ± 36.3	100–159
120–180 ( $\mathit{N}=9$ )	157.8 ± 16.8	1.77 ± 1.97	173.1 ± 47.9	137–209
180–240 ( $\mathit{N}=7$ )	218.4 ± 11.0	3.16 ± 2.01	206.8 ± 48.9	163–251
240–300 ( $\mathit{N}=8$ )	271.9 ± 15.6	5.14 ± 2.07	255.0 ± 50.3	214–296
300–360 ( $\mathit{N}=6$ )	333.8 ± 12.4	6.29 ± 1.36	282.8 ± 33.0	250–316

$\mathrm{\Delta }{\mathrm{T}}_{1\mathit{\rho }}$ reference-dependent ( $\mathrm{RMSE}=\pm 54\text{mins}$ )
Timing (mins)		Parameters	Estimation (mins)
Interval	Mean time	$\mathrm{\Delta }{\mathrm{T}}_{1\mathit{\rho }}$ (ms)	${\mathit{t}}_{\mathrm{estimated}}$	95% CI
000–060 ( $\mathit{N}=8$ )	41.9 ± 10.4	2.17 ± 1.28	74.3 ± 23.9	55–94
060–120 ( $\mathit{N}=8$ )	102.4 ± 11.6	5.38 ± 1.79	134.1 ± 33.3	107–161
120–180 ( $\mathit{N}=8$ )	160.4 ± 11.2	7.42 ± 2.96	172.1 ± 55.1	127–217
180–240 ( $\mathit{N}=8$ )	219.4 ± 12.9	9.21 ± 2.76	205.4 ± 51.3	164–247
240–300 ( $\mathit{N}=7$ )	274.1 ± 11.5	11.66 ± 3.25	251.0 ± 60.4	197–305
300–360 ( $\mathit{N}=6$ )	331.8 ± 12.4	13.03 ± 1.63	276.4 ± 30.3	246–307

Estimators based on ${\mathrm{T}}_2$ and ${\mathrm{T}}_{1\mathit{\rho }}$ signal intensities, averaged over their respective echo and spin-lock times, were also considered. However, they performed poorly because the statistical distributions lacked any significant time-dependence (Supplementary Information, Figure S1). A reference-independent estimator that combined both ${\mathrm{T}}_2$ and ${\mathrm{T}}_{1\mathit{\rho }}$ log-logistic parameters was also constructed. However, our results indicated that these were highly correlated, and therefore no significant improvement in precision was observed.

4. Discussion

The current study has established a stroke onset timing procedure using quantitative relaxometric MRI data only from the ischaemic lesion delineated by low ADC. The use of low ADC for selection of tissue for relaxometric data assures that only ischaemic volumes [7] are used for timing. Low ADC eliminates non-ischaemic regions from timing analyses, because reperfusion in the early hours of ischaemia normalises ADC [11,16]. It should be noted that quantitative ADC alone provides no time information [15], neither does the rate at which ADC lesion expands [13] in the early hours of ischaemia in permanent MCAO. This ‘anatomical reference-independent’ procedure performs better in timing than the one using MRI data from a non-ischaemic reference region. As the quantitative ${\mathrm{T}}_2$ (which is commonly available in standard clinical MRI systems) provides excellent precision for stroke time estimates in the preclinical model, the presented procedure is thought to bear significant potential for stroke patient stratification.

The current ${}^1\mathrm{H}$ MRI-based timing estimation procedures rely on the differences in average values of distributions of relaxation times between ischaemic tissue and non-ischaemic reference, i.e. $\mathrm{\Delta }{\mathit{T}}_2$ and $\mathrm{\Delta }{\mathit{T}}_{1\mathit{\rho }}$ . A substantial body of evidence shows that the time dependence of both $\mathrm{\Delta }{\mathit{T}}_2$ and $\mathrm{\Delta }{\mathit{T}}_{1\mathit{\rho }}$ in the rat brain can be adequately described by a linear function [9,20,25] and $\mathrm{\Delta }{\mathit{T}}_2$ similarly in human brain [27]. The $\mathrm{\Delta }{\mathit{T}}_2$ -based estimators have been reported to classify stroke onset within 3 hours with a precision of 0.850 in rats [19] and 0.757 in humans [27] indicative of clinical utility. However, the precision of timing by these reference-dependent estimators strictly depends on the selection of an anatomically identical non-ischaemic reference region for the ischaemic lesion. This is a non-trivial task even for the simple rat brain (see Supplementary Information Figure S3). In humans, reference selection will become even more problematic due to inherent relaxation time differences in GM and WM [30], variation of ${\mathrm{T}}_2$ in WM due to microstructure [14], as well as the lack of ‘perfect anatomical symmetry’ of hemispheres to choose a suitable ‘mirror reference’ over the midline. Furthermore, stroke lesions in humans vary in location and size.

The stroke time estimators introduced in this paper use only the ${\mathrm{T}}_2$ or ${\mathrm{T}}_{1\mathit{\rho }}$ data from low ADC-defined ischaemic brain tissue. Our data indicate that information embedded within the distributions of ${\mathrm{T}}_2$ and ${\mathrm{T}}_{1\mathit{\rho }}$ can be exploited to build multiple linear regression estimators independent of MRI data from a non-ischaemic reference, provided relaxation data from multiple time points are available. In other words, the time dependency of the parameters in regression estimators needs to be determined similar to the $\mathrm{\Delta }{\mathit{T}}_2$ and $\mathrm{\Delta }{\mathit{T}}_{1\mathit{\rho }}$ time dependencies used in the recent estimators [9,27]. The pathophysiological basis underpinning the estimators is likely to be that the volume of brain tissue with abnormally elevated MRI relaxation times increases linearly with time [20]. This is due to the effects of cytotoxic oedema on both ${\mathrm{T}}_{1\mathit{\rho }}$ [10] and ${\mathrm{T}}_2$ relaxation times with macromolecular breakdown [13] followed by vasogenic oedema on hydrodynamics within ischaemic tissue [28]. In addition to cytotoxic and vasogenic oedema, altered chemical exchange processes may influence ${\mathrm{T}}_{1\mathit{\rho }}$ relaxation in ischaemia [8]. The spread of effects of ischaemia in the brain tissue in the early hours of the insult has been proposed to be the critical factor underpinning the time-dependent increase in ${\mathrm{T}}_2$ relaxation time due to cytotoxic oedema [13]. We speculate that this spread is captured by the log-logistic fits of the relaxation time distributions to provide accurate stroke time estimation.

Regarding the utility of the current time estimator in stroke imaging, several points should be considered. Firstly, while both ${\mathrm{T}}_{1\mathit{\rho }}$ and ${\mathrm{T}}_2$ MRI are in use in general radiology, ${\mathrm{T}}_{1\mathit{\rho }}$ is not part of brain MRI protocol in clinics due to specific absorption rate of energy concerns. Therefore, the estimator based on ${\mathrm{T}}_2$ data may bear clinical potentials. Secondly, it should be stressed that the time estimators are generated from the quantitative MRI data obtained during early hours of permanent ischaemia. Thus, the current estimators provide timing only from the tissue where residual perfusion is below the ischaemia threshold [2]. A previous study showed that reperfusion of ischaemic tissue, but not of that with compromised blood flow and normal ADC, leads to prolongation of tissue ${\mathrm{T}}_2$ [12]. Against this, reperfusion of previously ischaemic tissue during MRI data acquisition would invalidate the established $\mathrm{\Delta }{\mathrm{T}}_2$ vs stroke onset time relationships [9,20] and potentially also the ${\mathrm{T}}_2$ relaxation time distribution characteristics used by the current log-logistic-based estimator. In clinical settings reperfusion, either spontaneous or pharmacologically induced, is most likely to take place when the clot is small and/or the occlusions involves small vessels; instead, reperfusion of proximal branches of cerebral arteries is uncommon [22]. Thirdly, in the current data set the number of voxels in ADC lesions varied from 180 to 1,330. A human stroke lesion of 5 cm³ in volume would have $\mathrm{>}1\mathrm{,}600$ voxels in a ${\mathrm{T}}_2$ maps acquired at routine spatial resolution of 3 mm³, i.e. even small strokes would have enough ${\mathrm{T}}_2$ data for generation of reference-independent estimators. Finally, a matter to be considered is that the ischaemic rat tissue is almost entirely gray matter and our time estimator apparently holds only for this tissue type.

5. Conclusions

The current study indicates that relaxometric MRI data may provide stroke time estimators without the need for non-ischaemic reference data. Critical issues in the clinical implementation of such protocols include both scan time and data post-processing. A recent study demonstrated a comprehensive stroke MRI protocol can be acquired in six minutes [24]. Investigational MRI scan procedures, such as magnetic resonance finger printing [18], have been shown to further shorten scan times, thus paving the way for expedited clinical brain examination. As far as the post-processing is concerned, the current standard computers are capable of analysing MRI data sets in seconds. Overall, quantitative MRI-based stroke timing taking minutes is likely to become reality for clinical exploitation in the near future.