Time to peak and full width at half maximum in MR perfusion: valuable indicators for monitoring moyamoya patients after revascularization

Moyamoya disease (MMD) is a chronic, steno-occlusive cerebrovascular disorder of unknown etiology. Surgical treatment is the only known effective method to restore blood flow to affected areas of the brain. However, there are lack of generally accepted noninvasive tools for therapeutic outcome monitoring. As dynamic susceptibility contrast (DSC) magnetic resonance imaging (MRI) is the standard MR perfusion imaging technique in the clinical setting, we investigated a dataset of nineteen pediatric MMD patients with one preoperational and multiple periodic DSC MRI examinations for four to thirty-eight months after indirect revascularization. A rigid gamma variate model was used to derive two nondeconvolution-based perfusion parameters: time to peak (TTP) and full width at half maximum (FWHM) for monitoring transitional bolus delay and dispersion changes respectively. TTP and FWHM values were normalized to the cerebellum. Here, we report that 74% (14/19) of patients improve in both TTP and FWHM measurements, and whereof 57% (8/14) improve more noticeably on FWHM. TTP is in good agreement with Tmax in estimating bolus delay. Our study data also suggest bolus dispersion estimated by FWHM is an additional, informative indicator in pediatric MMD monitoring.

Recently, DEFUSE-3 18 and DAWN 19 trials have suggested that deconvolution-based time to maximum (Tmax) parameter may be a more effective biomarker for identifying salvageable ischemic tissue from infarct core for selecting acute ischemic stroke patients. Using deconvolution with an arterial input function (AIF) selected from the contralateral middle cerebral artery, Tmax is regarded as an AIF-normalized bolus arrival delay time with an optimal threshold of 6 s for identifying critically hypoperfused tissue across studies and patients 20 . It is natural that Tmax has been examined in several MMD studies 6,21,22 . However, to the best of our knowledge, there has been no prior study of pediatric MMD longitudinal changes using Tmax. One possible reason is that the deconvolution process is sensitive to AIF selection 23 . It is highly challenging for clinicians to select the optimal AIF at any DSC MRI scans consistently given the hemodynamic complexities in MMD 24 .
In this study, we analyzed the hemodynamic transition in cerebral cortex tissue of 19 pediatric patients who were treated with indirect vascularization procedures and had one preoperational and multiple periodic DSC MRI examinations after operation over a time span of 4 to 38 months. Cortical areas of blood perfusion improvement by collateral development (arteriogenesis) 25 were evaluated by DSC MRI bolus delay and dispersion at a variety of time thresholds. Shortened delay indicates that new faster blood supplying routes are established while narrowed dispersion indicates that more direct routes, possibly via larger arteries, are supplying blood. We took a rigid gamma variate model (GVM) fitting approach 26 to derive two nondeconvolution-based perfusion parameters: time to peak (TTP) and full width at half maximum (FWHM) 27 for estimating bolus delay and dispersion respectively. In addition, GVM-derived bolus characteristics also allowed fully automated AIF selection for estimating Tmax 26 . Here we present the bolus delay and dispersion changes in pediatric MMD in terms of TTP, FWHM, and Tmax quantitatively and the correlation agreement analyses between TTP, Tmax, and FWHM by intraclass correlation coefficient (ICC) 28 .

Results
DSC-MRI data. The study data consisted of 19 patients (12 boys, 7 girls, no older than 19 years, mean age 10.7 ± 4.6), each with 3 to 8 scans (mean, 6.3 ± 1.4) for periods of time ranging from 4 to 38 months (mean, 16.3 ± 9.2). By excluding 11 low quality scans from statistical analysis, there were in total 108 scans (mean per person, 5.7 ± 1.5). Each patient had at least one of three indirect revascularization procedures: left encephaloduroarteriosynangiosis (EDAS), right EDAS, and multiple burr hole operations. All patients with Suzuki stages range from I to IV 29 . There were 14 ischemic, 5 no-stroke, and 0 hemorrhage cases in MRI imaging findings. Table 1 lists patient information with the chronological sequence of both the DSC-MRI scans and surgical operations. Note that most patients received 2 surgery treatments in different time order. The shortest latest effective scan was 3 months.
Nondeconvolutional perfusion parameter maps. Nondeconvolution-based MR perfusion parameters TTP and FWHM were estimated for all 119 DSC MRI scans (including low quality scans) listed in Table 1. For computational simplicity and consistency, the estimations were done only on an image slice near convexity as removal of cerebrospinal fluid was not required in this region. Figure 1 illustrates the image processing procedure and the resultant TTP and FWHM maps.
The computation in Fig. 1 is summarized as follows. Assuming a linear relationship between the tissue contrast agent concentration and the change in T2* relaxation rate, contrast concentration time curves C(t) were computed on a voxel-by-voxel basis without smoothing by: where S(t) denoted MR signal intensity time curves, S 0 the precontrast signal intensity, and TE the echo time. Let the maximal S(t) be S 0 and TE = 1 for computational simplicity. Level C 0 and bolus arrival time t 0 of C(t) were decided by a linear-linear model 30 . The level-adjusted concentration G(t) was modeled by: a GVM with four parameters A, t 0 , α, and β, where (TTP -t 0 ) = αβ by definition. The number of free parameters in G(t) were reduced to 2 by fixing the values of t 0 and TTP 31 . G(t) was solved efficiently by linear regression using data points between t 0 to t 2 (G(t 2 ) = half maximum) 26,31 .
TTP and FWHM were derived in high precision using the fitted model G(t) defined by Eq. (2) (see Supplementary Information for estimation error analysis). Note that in order to compare inter-scan results, all the TTP and FWHM were normalized by referring to the mean TTP c and FWHM c (subscription c for cerebellum) of the green boxes at cerebellum manually selected by operators. TTP and FWHM normalization formulas were empirically selected such that TTP, FWHM, and Tmax perfusion color maps were in good agreement. (Refer to Methods Section for further details). Longitudinal monitoring examples. Figure 3 shows an example of a 12-year-old boy (patient 17 in Table 1 and Fig. 2) who underwent right EDAS procedure first and left EDAS five months later. While the  www.nature.com/scientificreports/ FWHM/2 (rescaled, half-FWHM) by sequence appears consistent, the normalized FWHM n clearly indicates that tissue bolus time curves on the right hemisphere become narrower gradually but those on the left become wider at month + 3, and narrower thereafter. Tmax and TTP also show similar longitudinal patterns. Figure 4 is another example of an 11-year-old boy (patient 12 in Table 1 and Fig. 2) who underwent left EDAS procedure first and right EDAS five months later. TTP and Tmax maps at + 0 (a few days before the first operation) and + 17 months appear similar as if the delay time parameters only fluctuate without major changes. However, FWHM n maps illustrate that the left hemisphere progresses consistently throughout the observation time window with observable lowering FWHM n even 17 months after the initial left EDAS operation. Digital subtraction angiography (DSA) images of the same patient performed in months + 0 (before left EDAS procedure), + 5 (after right EDAS procedure), and + 17 (follow-up) are illustrated in Fig. 5. The DSA patterns suggest that the FWHM n parameter has an inverse relationship with arteriogenesis stimulated by the left EDAS procedure.
Correlation agreement between TTP and Tmax. We conducted a set of eight correlation agreement analyses between the proposed rigid GVM-derived TTP and eight different Tmax derivations using all 108 clean DSC MRI scans. The full ICC results of the TTP and Tmax areas at time thresholds 2, 3, 4, 5, and 6 s at convexity are listed in Table 2 respectively.
The deconvolution for estimating Tmax was implemented by using Singular Value Decomposition (SVD) in time domain 32 and Fast Fourier Transform (FFT) in frequency domain 33 and noted as FFT-, FFT + , SVD-, and SVD + , where ' + ' and '−' representing with ( +) and without (−) using the proposed GVM preprocessing. Three top-ranked AIFs were automatically selected by in-house software 26 to derive three Tmax estimates for each FFT and SVD derivation. These three estimates were then used to derive the average (A) and median (M) Tmax at each voxel. In total, Table 2 lists eight Tmax's which are noted as FFT -A , FFT -M , FFT +A , FFT +M , SVD -A , SVD -M , SVD +A , and SVD +M according to their derivation approaches.
The level of agreement was interpreted by the guidelines proposed by Koo et al. 28 : ICC < 0.50 (poor), 0.5-0.75 (moderate), 0.75-0.90 (good), and > 0.90 (excellent). Poor to moderate agreement was found between all eight Tmax's and TTP at thresholds 2 and 3 s for most computation settings in Table 2. Only one case had moderate to good agreement with TTP at threshold 3 s. It was SVD +M with ICC 0.67 (95% confidence intervals (CI): 0.53, 0.77). Overall, Tmax estimated by SVD +A and SVD +M were in best agreement with TTP. SVD +M was above the average at thresholds 4 (ICC 0.79; 95% CI: 0.70, 0.85), 5 (ICC 0.85; 95% CI: 0.78, 0.89), and 6 (ICC 0.90; 95% CI: 0.84, 0.93) seconds, which were moderate to good, good, and good to excellent respectively. Figure 6a shows the correlation scatterplots for TTP and Tmax (derived by SVD +M ) at thresholds 1, 2, 3, 4, 5, and 6 s. Figure 6b shows the correlation scatterplots for TTP and FWHM n at thresholds 1, 2, 3, 4, 5, and 6 s respectively. The data points are more scattered in Fig. 6b as compared with Fig. 6a. The ICC agreement evaluations between perfusion delay time areas estimated by FWHM/2, FWHM n , and TTP at thresholds 2, 3, 4, 5, and 6 s are listed in the end of Table 2, where FWHM is divided by 2 (noted FWHM/2) for scaling the time measurements down to a range roughly between 0 and 6 s. The highest ICC was 0.57 (95% CI: 0.43, 0.69) at time threshold 6 s for FWHM/2 and 0.69 (95% CI: 0.58, 0.78) at time threshold 3 s for FWHM n respectively. FWHM n is in poor to moderate agreement, but FWHM/2 is prominently in poor agreement with TTP. www.nature.com/scientificreports/

Discussion
TTP and Tmax are two popular perfusion parameters in monitoring cerebral perfusion changes after revascularization in patients with MMD using DSC MRI. Both parameters are shown to be sensitive enough to monitor transitional bolus delay time changes in several studies 6,14,15,21,22 . Recent large-scale multi-center trials have suggested that the deconvolution-based Tmax parameter may be a more effective biomarker with an optimal threshold of 6 s for identifying critically hypoperfused tissue across studies and patients [18][19][20] . However, deconvolution-based Tmax is known to be sensitive to AIF selection. It is particularly challenging for longitudinal MMD monitoring as the contralateral side is rarely clear given a MMD patient with multiple operations at different time.
To alleviate the uncertainty of AIF selection and to serve as an alternative viewpoint, we proposed using a rigid GVM approach to derived cerebellum-normalized bolus delay (TTP) and dispersion (FWHM n ) parameters for longitudinal MMD monitoring, assuming no vertebral arterial abnormality existed to affect cerebellum circulation. The resultant TTP was statistically found to be in good agreement with Tmax in delay time map area measurement (Figs. 3, 4, 6a). The resultant FWHM n was shown to be in less agreement with TTP (Figs. 3, 4, 6b) but equally sensitive to detect bolus contrast time curve improvement. In Fig. 2b Fig. 2. Therefore, our method is likely available for Suzuki grades between I and IV. Regarding patients 13 and 18, who had no improvement in either TTP or FWHM n , we notice that both patients had operation performed within 3 months right before the latest DSC-MRI scans. We need further scans to evaluate these two cases more accurately. (See Supplementary Information for the perfusion maps of patients 13 and 18).  Table 1) who underwent a right EDAS operation first and a left EDAS five months later. Normalized delay time areas at thresholds of 2, 4, and 6 s are given under each map. Tmax (first row, derived by SVD +M ) shows modest perfusion improvement while TTP (second row) indicates more improvement. FWHM/2 (third row) shows consistent pattern in green color. FWHM n (fourth row) becomes narrower gradually on the right hemisphere (arrow head) but wider at month + 3 on the left hemisphere (arrow), and narrower thereafter. www.nature.com/scientificreports/ Note that a variety of Tmax implementations exists and FWHM n does not always change in accordance with TTP (Table 2), perfusion delay time parameters for MMD monitoring should be implemented and interpreted with a few caveats discussed in the following.
First, using a rigid GVM preprocessing approach to reduce scan noise-caused error, the fitted model in Eq. (2) is more robust and precise for deriving both bolus delay (TTP) and dispersion (FWHM) estimations. In addition, GVM can remove the recirculation effect from DSC MRI data for deconvolution computation. As a  Table 1) who underwent a left EDAS operation first and a right EDAS five months later. Normalized delay time areas at thresholds of 2, 4, and 6 s are given under each map. Both Tmax (derived by SVD +M ) and TTP show that perfusion on the left hemisphere improved between months + 5* to + 8 but returned to the pre-surgery condition (+ 0*) after month + 11. However, FWHM n indicates that the contrast time curves on the left hemisphere have become narrower (arrow) compared to pre-surgery (arrow head).  Table 2 shows that Tmax estimates by FFT + and SVD + are in better agreement with TTP than FFTand SVDrespectively. Nonetheless, we recommend applying SVD + only for numerical stability reasons because the resultant G(t) defined by Eq. (2) could cause wiggles in some FFT + computations. Second, using an empirical formula: TTP = TTP 0 -TTP c + 1, we were able to normalize raw TTP 0 to the similar time range of Tmax with an optimal threshold of 6 s by referring to cerebellar TTP c estimated from manually selected cerebellar patches. The " + 1" adjustment was empirically decided to achieve best ICC scores for all tested Tmax computations. Tmax estimated by SVD +M was in moderate to excellent agreement with TTP at delay time thresholds between 4 and 6 s that was of clinical importance. The normalization formula of FWHM n was empirically decided in a similar fashion.
Third, the clinical indications and interpretations of Tmax or TTP should not be conducted alone without comparing with other perfusion parameters. As a reminder, the scans (Fig. 4) of patient 12 at scanning time + 0 and + 17 months may appear similar to human eyes. Without reviewing other perfusion parameters, the radiological diagnosis might state that the TTP and Tmax delay time maps were fluctuating but stable without major changes. However, FWHM n illustrated a complete different scenario that the left hemisphere progressed consistently throughout the observation time window with observable FWHM n narrowing. By reviewing DSA images (Fig. 5), we found that FWHM n narrowing agreed with local blood flow improvement in an apparently inverse relationship. Table 2. ICC agreement analysis between Tmax/FWHM and TTP using the normalized delay time areas at thresholds of 2, 3, 4, 5, and 6 s. * MRP stands for the estimated MR perfusion parameters Tmax and FWHM. Tmax estimated by different deconvolution methods are noted with FFT, SVD with ( +) and without (−) GVM preprocessing; the characters ' A' and 'M' stand for the average and median of Tmax estimates derived from the top three AIFs respectively. www.nature.com/scientificreports/ Our study had two main limitations. First, we did not validate the accuracy of the automatically detected AIFs because accurate determinations of true AIF were challenging for MMD patients who were in different disease stages and underwent a variety of revascularization procedures. Instead, we conducted a series of ICC analyses between TTP and 8 different Tmax derivations with 3 automatically top-ranked AIFs. The statistical ICC analyses indicated that the automatically selected AIFs were relatively consistent and accurate for Tmax estimations. Second, our data size was too small for statistical grouping analysis with age, sex, associated disease stages, operative methods, and clinical results. In spite of the aforementioned limitations, our results indicated that nondeconvolution-based time to peak and peak width changes were observable perfusion characteristics by using cerebellar normalization.
In conclusion, Good to excellent agreement between time-to-maximum and time-to-peak MR perfusion delay time evaluations of pediatric moyamoya disease can be achieved by using a rigid gamma variate model fitting approach. Our study data also suggest bolus dispersion estimated by properly normalized FWHM may be an additional, informative indicator in pediatric MMD monitoring. The innovation of this study can be applied to clinical practice directly without changing current DSC-MRI scanning protocols.

Methods
This study and all experimental protocols were approved by the institutional review board of National Taiwan University Hospital, Taipei, Taiwan. Written informed consent was obtained from all participants and the children's parents. All methods were carried out in accordance with relevant guidelines and regulations. Pharma, Berlin, Germany) was injected 5 s after the DSC MRI commenced at 3 cm 3 /s injection rate and followed by a 20-cm 3 saline flush. The raw DSC MRI data was converted to contrast concentration using Eq. (1) and preprocessed by conventional image segmentation procedures to remove non-brain tissues.

Materials.
Rigid GVM fitting. It is possible to do least-squares fitting of Eq. (2) directly, however, it requires a multiple linear regression strategy which is computationally more expensive and numerically less stable as exponential functions are involved. To simplify the computation and make the GVM fitting more robust, here we present a 2-stage rigid GVM computation strategy, which is illustrated as C(t) and G(t) in Fig. 1.
The precontrast level C 0 and the bolus arrival time t 0 in Eq. (2) were estimated by a linear-linear-model LLM(t). In this model, the portion of contrast concentration time curves C(t) for 0 ≤ t ≤ TTP ′ (red circled points in C(t) with the maximum at t = TTP', Fig. 1) were approximated by a linear-linear piecewise function 26,30 : such that the sum of squared error SSE(t 0 ) was minimized. Substituting t by x in Eq.
(3) can be simplified as : Parameters C0 and C1 in Eq. (5) can be solved very efficiently by a linear least-squares estimation. Therefore, Eq. (4) can be minimized by linearly searching along t 0 from 0 to TTP' with an appropriate small increment.
Stage 2. The level-adjusted concentration G(t) in Eq. (2) was commonly modeled by four parameters A, t 0 , α, and β: where (TTP − t 0 ) = αβ by definition. The number of free parameters in G(t) can be reduced to 2 by fixing the values of t 0 and TTP 26,31 . Rearrange G(t) with t' = (t-t 0 )/(TTP − t 0 ) and G max = A(TTP-t 0 ) α exp(-α) at t = TTP, www.nature.com/scientificreports/ Take the natural logarithm of both sides of G(t'): The new formula has the form of y = b + ax, where G max and α can be determined from linear regression of the natural logarithm of the observed values. In this model, the portion of contrast concentration time curves G(t) for t 0 < t ≤∼ t 2 (red circled points between t 0 and roughly estimated t 2 in G(t), Fig. 1) were approximated. We name this approach as rigid GVM because t 0 and TTP are fixed in each regression model. Perfusion parameter estimation. MR perfusion parameter estimation was done on an image slice near convexity as removal of cerebrospinal fluid was not required in this region. On this convexity image slice, absolute areas were assessed at delay time thresholds 1, 2, 3, 4, 5 and 6 s respectively and normalized by the total cortical area. Tissue bolus-derived delay time parameters Tmax, TTP, and dispersion parameter FWHM are delineated as follows.
Tmax was estimated by Fast Fourier Transform (FFT) in frequency domain 33 and Singular Value Decomposition (SVD) in time domain 32 , with and without GVM preprocessing. The resultant residue curves were up-sampled at 0.5 s time resolution to find Tmax 33 . The top three AIF candidates detected by a fully automated algorithm 26 were used to compute three individual Tmax estimates at each pixel. The mean (average) and median values of these three estimates were also computed and analyzed.
TTP was determined as the best fitting result to the rigid GVM, ln(G(t')) defined by Eq. (7). For computational speed, we kept the LLM-solved t 0 fixed and only conducted linear research for TTP in the neighborhood of the maximal C(t) with a time step of 0.1 s. The raw result, TTP 0 , was normalized by subtracting the average result, TTP c , estimated at manually selected cerebellum areas: The adjustment of 1 s was determined empirically to correlate well with Tmax.
FWHM, peak width measured at half maximum (FWHM = t 2 -t 1 as illustrated in Fig. 1), was derived from the resultant GVM. In addition, a normalized version FWHM n referring to FWHM estimated in the cerebellum was also defined empirically as All GVM preprocessing and delay time computation were performed by in-house software implemented on the Matlab platform (version 2019b; Mathworks, Natick, MA).

Data availability
The datasets for this study are protected patient information. Some data may be available for research purposes from the corresponding author upon reasonable request.