Semi-Parametric Arterial Input Functions for Quantitative Dynamic Contrast Enhanced Magnetic Resonance Imaging in Mice Torfinn Taxt, Rolf K. Reed, Tina Pavlin, Cecilie Brekke Rygh, Erling Andersen, Radovan Jiřı́k PII: DOI: Reference: S0730-725X(17)30227-8 doi:10.1016/j.mri.2017.10.004 MRI 8848 To appear in: Magnetic Resonance Imaging Received date: Revised date: Accepted date: 30 January 2017 15 September 2017 17 October 2017 Please cite this article as: Taxt Torﬁnn, Reed Rolf K., Pavlin Tina, Rygh Cecilie Brekke, Andersen Erling, Jiřı́k Radovan, Semi-Parametric Arterial Input Functions for Quantitative Dynamic Contrast Enhanced Magnetic Resonance Imaging in Mice, Magnetic Resonance Imaging (2017), doi:10.1016/j.mri.2017.10.004 This is a PDF ﬁle of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its ﬁnal form. Please note that during the production process errors may be discovered which could aﬀect the content, and all legal disclaimers that apply to the journal pertain. T ACCEPTED MANUSCRIPT RI P Semi-Parametric Arterial Input Functions for Quantitative Dynamic Contrast Enhanced Magnetic Resonance Imaging in Mice SC Torfinn Taxta,b , Rolf K. Reeda,c , Tina Pavlina,b , Cecilie Brekke Rygha , Erling Andersend , Radovan Jiřı́ke a Dept. of Biomedicine, University of Bergen, Jonas Lies vei 91, N-5020 Bergen, Norway of Radiology, Haukeland University Hospital, Jonas Lies vei 83, N-5020 Bergen, Norway c Centre for Cancer Biomarkers (CCBIO), University of Bergen, Jonas Lies vei 87, N-5021 Bergen, Norway d Dept. of Clinical Engineering, Haukeland University Hospital, Jonas Lies vei 83, N-5020 Bergen, Norway e Czech Academy of Sciences, Inst. of Scientific Instruments, Královopolská 147, 61264 Brno, Czech Rep. ED MA NU b Dept. Abstract PT Objective: An extension of single- and multi-channel blind deconvolution is presented to improve the estimation of the arterial input function (AIF) in quantitative dynamic contrast enhanced magnetic resonance imaging (DCE- AC CE MRI). Methods: The Lucy-Richardson expectation-maximization algorithm is used to obtain estimates of the AIF and the tissue residue function (TRF). In the first part of the algorithm, nonparametric estimates of the AIF and TRF are obtained. In the second part, the decaying part of the AIF is approximated by three decaying exponential functions with the same delay, giving an almost noise free semi-parametric AIF. Simultaneously, the TRF is approximated using the adiabatic approximation of the Johnson-Wilson (aaJW) pharmacokinetic model. Results: In simulations and tests on real data, use of this AIF gave perfusion values close to those obtained with the corresponding previously published nonparametric AIF, and are more noise robust. Conclusion: When used subsequently in voxelwise perfusion analysis, these semi-parametric AIFs should give more correct perfusion analysis maps less Preprint submitted to Magnetic Resonance Imaging October 20, 2017 T ACCEPTED MANUSCRIPT RI P affected by recording noise than the corresponding nonparametric AIFs, and AIFs obtained from arteries. Significance: This paper presents a method to increase the noise robustness SC in the estimation of the perfusion parameter values in DCE-MRI. Keywords: DCE-MRI, blind deconvolution, arterial input function NU 1. Introduction Quantitative DCE-MRI is an important tool in clinical disciplines such as MA oncology [1, 2, 3, 4], cardiology [5, 6, 7], and in biological research [8, 9]. DCEMRI is used in rodents to study perfusion changes in malignant tumors following treatment with new vasoactive drugs [10] and changes in perfusion of animal ED strains with knock out deletions of genes. To do quantitative DCE-MRI, the observed contrast tissue signals have to be converted to contrast concentration time sequences (perfusion sequences) [11]. PT Each perfusion sequence is a convolution of an AIF and a TRF multiplied by plasma flow [12]. Thus, the perfusion sequence has to be deconvolved with the AC CE AIF to find the TRF with its associated perfusion parameters. There are several approaches to estimate the AIF (a comprehensive review in [13]). The AIF can be derived from one or more artery signals close to the re- gion of interest (ROI). However, AIFs from arteries are degraded by saturation effects. For high peak concentration of the contrast agent in the arteries, the R1 relaxivity rate is no longer proportional to the contrast agent concentration [14]. High contrast agent concentration in arteries leads also to non-negligible T2* effects decreasing further the level of the measured arterial signal [15]. A measured AIF can also be degraded by flow artifacts and partial volume effects [16]. Furthermore, representative arteries are often difficult to obtain in the acquired image sequence, especially in preclinical applications. Finally, dispersion (change of the AIF shape due to the passage of the contrast-agent bolus from the site of the measured arterial signal to the tissue of interest) is ignored. A population based AIF can also be used [17]. This ignores the differences in 2 T ACCEPTED MANUSCRIPT RI P the vascular tree between different subjects and depends on the AIF acquisition methods used for creation of these population-based ”standards”. Another approach of AIF estimation is based on a reference tissue (e.g. mus- SC cle) [18]. The AIF is estimated from the tissue curve in the reference tissue and the presumably known perfusion parameters. This approach has been shown for the Tofts model. However, for more complex pharmacokinetic models described NU by more perfusion parameters, the complete set of perfusion parameters would have to be known, which is not very realistic. The AIF estimation approach used here is based on single- and multi-channel MA blind deconvolution [19, 20, 21]. The AIF is estimated simultaneously with the TRF from the measured tissue contrast agent concentration curves. This provides examination specific AIF estimates. These techniques avoid many of ED the problems described above. Because the contrast agent concentration in the tissue is much lower than in arteries, the saturation effects mentioned above for measured AIFs are avoided. PT The algorithm can be implemented as single- [22, 23] or multi-channel [24, 25] deconvolution, where the number of channels is the number of tissue region sig- AC CE nals processed simultaneously. Multi-channel algorithms rely on the assumption that the AIF is the same for all processed ROIs. This might be regarded as an additional prior information if the assumption is true, or as a source of error if the tissue signals differ in the dispersion term of the local tissue region specific AIFs. This also means that single-channel blind deconvolution provides the local AIF (accounting for dispersion) while the multi-channel technique leads to an AIF estimate specific for the whole region containing all channels. However, for channels selected from a small region, multi-channel blind deconvolution can also provide an estimate of a local AIF [26]. Blind deconvolution AIF estimation has been introduced in [24, 21, 25] for clinical DCE-MRI and extended later to preclinical DCE-MRI [10, 22, 27, 28, 23]. To obtain an AIF using blind deconvolution, a realistic pharmacokinetic model of the TRF for the tissue in the ROI has to be selected. Different pharmacokinetic models have a clear impact on the time course of the estimated 3 T ACCEPTED MANUSCRIPT RI P AIF [23]. The pharmacokinetic models of the TRF used in blind deconvolution are largely based on Tofts [24] and extended Tofts [21, 25] models. More advanced pharmacokinetic models used for blind deconvolution include the two SC compartment exchange model (2CXM [1]) [23], the adiabatic approximation to the tissue homogeneity model (aaJW [29]) [22, 23], the distributed capillary adiabatic tissue homogeneity model (DCATH [30]) [31] and the Gamma Capillary NU Transit Time (GCTT [32]) [28] models. Blind deconvolution has been applied with a nonparametric AIF [24, 10, 22, 23] or with a parametric AIF [21, 25, 28, 31]. Clinical applications of blind MA deconvolution have been based on Parker’s AIF model (sum of two Gaussian functions and a sigmoid multiplied by an exponential function) [31], Fluckiger’s AIF model (Parker’s model with the Gaussian functions replaced by gamma ED variate functions) [21], and its modifications where three gamma variate functions were used instead of two, while forcing some parameters to the same values [25]. The gamma variate AIF model of [25] has been used also for preclinical PT DCE-MRI [28]. A realistic AIF model is an important prior information that can improve AC CE the reliability of blind deconvolution, compared to use of a nonparametric AIF. On the other hand, the use of a less realistic AIF model might degrade the accuracy of blind deconvolution. This paper is focused on small animal DCEMRI using an advanced TRF model (namely the aaJW model). For these applications, estimation of intravascular perfusion parameters, such as plasma flow and volume is very sensitive to the initial part of the AIF and of the measured tissue concentration sequence. Hence, it is important to estimate the initial part of the AIF accurately. In this paper, a semi-parametric AIF formulation is proposed to keep a suffi- cient degree of freedom for the initial AIF part by its nonparametric formulation and to regularize estimation of the slowly decaying AIF part by its parametrization. This approach is proposed as an alternative to the existing blind deconvolution AIF estimation methods which are based either on a parametric or a nonparametric AIF. 4 T ACCEPTED MANUSCRIPT RI P Estimation of the semi-parametric AIF is a straightforward extension of the previously published method of blind deconvolution for estimation of a nonparametric AIF [22]. As in [22], the algorithm is based on the Lucy-Richardson SC expectation-maximization algorithm and provides estimates of the AIF and the TRF. Simulated perfusion sequences were studied first. The AIFs and TRFs were NU derived from the inner part of the masseter muscle and the temporalis muscles of normal mice and mice treated with the anaphylactic agent C48/80. This agent gives a huge increase in blood flow and the permeability surface area product MA of muscle, thus leading to a large increase in the signal to noise ratio (SNR) of the muscle DCE-MRI signals [33, 34, 35, 36, 23]. In the second part, real perfusion sequences from the masseter and tem- ED poralis muscle of normal mice and mice treated with the anaphylactic agent C48/80 were studied. Perfusion parameter estimates derived using the aaJW compared. PT model for the TRF and either the nonparametric or semi-parametric AIF, were In the simulations and using real data, the semi-parametric AIF formulation AC CE was more noise robust than the nonparametric AIF. On the other hand, the nonparametric AIF formulation was more accurate in high SNR cases. 2. Theory and Methods 2.1. Perfusion Analysis in DCE-MRI The observed perfusion sequence, Ct [n] (n is the time index), is a convolution of the AIF and TRF, Ct [n] = Cp [n] ∗ H[n]. (1) Here, Cp [n] is the AIF and H[n] is the TRF multiplied by plasma flow, Fp . To find the relevant perfusion parameters in DCE-MRI, H[n] has to be modeled and isolated from Cp [n] through deconvolution and parameter estimation. To estimate Cp [n] and H[n] simultaneously, blind deconvolution is used. 5 T ACCEPTED MANUSCRIPT RI P Here, an extension of a previously published single- and multi-channel blind deconvolution method [22] is used to solve this problem. The single-channel algorithm is an iterative algorithm (Fig. 1) consisting of two parts. First, the SC TRF is updated based on the current AIF and TRF estimates, and then the AIF is updated based on the current AIF and TRF estimates. In the first 200 iterations, only nonparametric AIF and TRF are used. In the following 90 NU iterations, updating of the nonparametric TRF is followed by its approximation with a parametric aaJW model and update of the nonparametric AIF is followed by its approximation with the semi-parametric model. 90 iterations were found MA to be sufficient (200 iterations were used in [22]). More iterations did not lead to any visible change in the AIF and TRF estimates, indicating that convergence has been reached. ED The standard Lucy-Richardson deconvolution algorithm is formulated for non-negative signals that start with zero values and end with zero values. This leads to the need for extrapolation of the observed trunctated perfusion se- PT quences [22] to the length of approximately 3 hours (wash out time expected for a small molecular weight gadolinium contrast agent). This is done as in [22] AC CE by approximation of the upper time tail of the tissue function with a gamma function. A least-mean-square error procedure is used. The multi-channel version of the blind deconvolution algorithm is a direct extension of the single-channel approach [22]. The first 200 iterations (nonparametric deconvolution part) are performed as single-channel, separately for each channel. The channel-specific AIF estimates are merged in the step ”SemiParametric approximation of nonparametric AIF” (Fig. 1) as follows. All AIFs are first delay corrected (see [22] for more details) and averaged prior to the semi-parametric approximation. 2.2. Nonparametric Updates of the TRF and AIF A modification of the basic algorithm [22] is introduced to give improved estimates of the nonparametric AIFs and TRFs. Extensive preliminary experiments showed less distortion of the decaying part of the AIF due to noise when 6 T ACCEPTED MANUSCRIPT RI P using the regularization term for the Lucy-Richardson algorithm of [37], page 1062, instead of the total variation regularization term used in the basic algorithm. For simplicity, a low-pass filtered version of the residual defined in [37] to the following update equations. Cps+1 [n] = { Ct [n]+R̄1s [n] Cps [n] ∗ Cps [−n]}H s [n], (2) ∗ H s+1 [−n]}Cps [n]. (3) NU H s+1 [n] = { SC is used as a special case of the wavelet transform approach in [37]. This leads Ct [n]+R̄2s [n] H s+1 [n] MA Here, s is the iteration index, H s [n] is the estimation of H[n] in the s-th iteration, Cps [n] is the estimation of the AIF in the s-th iteration and R̄is [n] is the low-pass filtered residual (filter implemented as a shifted sigmoid in the π 40 rad/sample). The residual is R1s [n] = ED frequency domain, cutoff frequency Ct [n] − Cps [n] ∗ H s [n] for Eq. (2) and R2s [n] = Ct [n] − Cps [n] ∗ H s+1 [n] for Eq. (3). The AIF is initialized to the tissue perfusion sequence, Ct [n], after addi- PT tional median filtering (length 50 samples here, found experimentally). The TRF is initialized as a constant for the whole time interval of the extrapolated AC CE tissue function (including the precontrast part). More details on the choice of initialization can be found in [22]. For some of the noisy perfusion sequences of the controls, the parametric TRF did not overlap well with the corresponding nonparametric TRF. Also, the fit of the estimated perfusion sequence to the observed perfusion sequence was suboptimal. Obviously, the convergence point of the algorithm was local and wrong. Additional systematic variation of the optimization parameters always resulted in a global optimum point with the best fit of the observed and estimated time sequences and a good overlap of the final nonparametric and parametric TRFs. The systematically changed parameters were the intravascular transit time, Tc (given in samples: 6 to 17) and the number of samples, N o, (100,200,300 and 400), of the tail of the observed signal sequence used to estimate the extrapolation of this sequence in the preprocessing step. For all pairs of Tc and N o values from the above given ranges, the whole perfusion 7 T ACCEPTED MANUSCRIPT RI P estimation algorithm was run. The pair giving the best fit of the observed and estimated time sequences was then selected. 2.3. Parametric Approximation of the TRF SC The nonparametric TRF is approximated with the aaJW model (described below). In the first iteration with approximation (s=201, Fig. 1), the delay of NU the parametric TRF is estimated. The vertical rising edge of the parametric TRF is aligned with the point of the steepest rise of the nonparametric TRF, estimated as the maximum of the difference-operator processed sequence (see MA [22] for more details). The delay corresponding to this alignment is then applied in all following iterations. The aaJW model of the TRF is selected here as one of the simplest models ED providing a complete description of the intravascular and extravascular distribution of the contrast agent as a function of time. Other alternatives include the 2CXM [1], the tissue homogeneity model (TH) [38], the distributed parameter PT model (DP) [39], the distributed capillary adiabatic tissue homogeneity model (DCATH) [30], and the Gamma Capillary Transit Time (GCTT) model [32]. AC CE The TH and DP models are computationally demanding options because of no time domain analytic solution of the TH model and the need for numeric integration for the DP model. The DCATH and GCTT models include an additional perfusion parameter describing the distribution of the mean transit times. This may lead to worse conditioning of the deconvolution problem. A recent comparison of the aaJW and the 2CXM models indicated that the aaJW model is more realistic for cross-striated muscle (used here for evaluation) [23]. For cross-striated muscles, the plug flow assumption induced to the intravascular space by the aaJW model is more realistic than the compartment assumption of the 2CXM. For the aaJW model, H[n] is modeled as [29]: H[n] = Fp , 0 ≤ nTs < Tc , H[n] = EFp · e−kep ·nT s , nT s ≥ Tc , 8 (4) EFp ve is the rate of RI P where Fp is plasma flow, Ts is the sampling interval, kep = T ACCEPTED MANUSCRIPT contrast agent transport from the extravascular extracellular space to plasma, Tc is the mean intravascular transit time, ve is the extravascular extracellular SC volume per unit volume of tissue (leakage space) and E ≤ 1 is the extraction fraction. The first row of (4) is the vascular phase and the second row is the parenchymal phase of H[k]. NU The extraction fraction, E, is the fraction of the contrast agent that diffuses unidirectionally from plasma into the surrounding tissue during a single contrast-agent transit through the capillary bed. It is related to Fp and the per- MA meability surface area product, P S, through [40] E = 1 − e−P S/Fp . The aaJW model is completely described by the four independent physiological parameters ED Fp , Tc , E and kep . 2.4. Semi-Parametric Approximation of the AIF The semi-parametric AIF formulation consists of a nonparametric initial PT part, as returned by the nonparametric update, followed by a decaying parametric part. This part is a sum of three decaying exponential functions with AC CE the same delay, i.e. the time from which the approximation by the model starts. Hence, the parametric part of the AIF is described by six parameters (scaling factor and time constant for each exponential) and the common delay. In the first iteration where the semi-parametric approximation of the AIF is done (s=201, Fig. 1), the common delay of the three exponential decaying functions is determined as the point of the steepest rise of the nonparametric AIF estimated as the maximum value of the difference operator processed nonparametric AIF. In all following iterations, this delay is calculated as the location of the maximum point of the nonparametric AIF. This approach gave the best results, based on preliminary experiments. In extensive preliminary experiments, only two decaying exponential func- tions with the same delay were used to approximate the decaying phase of the AIF [41, 42]. This approach always led to lower values of the extraction fraction than the corresponding nonparametric AIF, while the other perfusion parame9 T ACCEPTED MANUSCRIPT RI P ters remained fairly unchanged. Thus, two decaying exponential functions were considered inadequate for a valid approximation of the decaying part of the AIF. The six parameters of the semi-parametric AIF are initialized in the first SC semi-parametric AIF approximation step (s=201, Fig. 1). It is done using the successive subtraction method [43]. Then, these initial parameter estimates are used as starting values for the nonlinear least mean square error method of NU Levenberg-Marquardt (the algorithm mrqmin in Numerical Recipes in C [44]). In the following iterations, the parameter values of the previous iteration are MA used as starting values for the Levenberg-Marquardt method. 2.5. Scaling Blind deconvolution estimates the AIF and TRF except for a scaling factor. ED The scaling factor is estimated at the end of the blind deconvolution algorithm (Fig. 1) from a reference tissue, the masseter muscle (see Discussion for other approaches to AIF scaling). Based on literature values of ve [45, 46, 47] and PT the assumption that blood plasma volume, vp , is approximately 5% of ve , the reference value ve + vp = 0.13 ml/ml is used. According to the formulation of AC CE the TRF model (Eq. (4)), ve + vp equals the area under the TRF curve. From the definition of the convolution of two non-negative signals, it follows PN that the area under the curve of a result of a convolution (here n=1 Ct [n], where N is the number of samples of the extrapolated perfusion sequence) is the product of the areas under the curves of the convolution factors (here PN PN n=1 Cp [n] and n=1 H[n]). Hence, the sought scaling factor of the blind de- convolution result, i.e. the area under the curve of the AIF, is the area under the curve of the contrast agent concentration sequence measured in a masseter P muscle region, Ct ref [n], divided by 0.13 ml/ml. P To account for the variability in the real recordings, the median of the values Ct ref [n] obtained from all real recordings was used. This relies on the fact that the same contrast-agent dose and concentration were used for all mice and that approximately the same kidney function can be assumed for both normal and treated mice. 10 T ACCEPTED MANUSCRIPT RI P 2.6. Simulations Time sequences were made to simulate the perfusion pattern of the masseter and temporalis muscles of control mice and mice treated with C48/80 (see SC below), with and without added noise. Thus, ignoring noise, four different perfusion sequences were generated, mimicking the perfusion sequences obtained NU from the real data recorded on one control and one treated mouse (Section 2.7). The Four Basic Perfusion Sequences. Typical aaJW TRFs of the masseter and temporalis muscles of the control and treated mice were generated using the MA reference perfusion parameters (Table 1). These perfusion parameters were obtained from in vivo mice perfusion sequences (see below) processed with the single-channel blind deconvolution algorithm using the nonparametric AIF for- ED mulation. The length of the simulated TRFs was 16384 samples, i.e. the length of real perfusion sequences after extrapolation and additional zeropadding for use of the fast Fourier transform, see Section 2.7. PT The AIFs used for simulation were obtained from in vivo data after preprocessing (conversion to R1, median filtering and extrapolation, see below). AC CE Single arterial voxels selected manually in the lower part of both hemispheres were used. The AIFs of the control mice were all much degraded by noise. To obtain a usable AIF of these mice, AIFs from four different control mice were selected and averaged. The AIFs of the treated mice were clearly less degraded by noise and artifacts. Hence, only a single, typical AIF was selected for the treated mice. The two TRFs of the control mice were convolved with the AIF of the control mice. In the same way, the two TRFs of the treated mice were convolved with the AIF of the treated mice. No dispersion was modeled. The resulting perfusion sequences were almost noise free (the remaining noise was due to the noise in the AIFs) and consisted of 16384 samples (the length of real perfusion sequences after extrapolation and additional zeropadding for use of the fast Fourier transform) with sampling interval 1.079 s. These perfusion sequences were denoted ”inf-full” on the SNR axis (Figs. 3, 5). No extrapolation 11 T ACCEPTED MANUSCRIPT RI P was applied in processing of these sequences. To evaluate also the effect of extrapolation on the blind deconvolution performance, the simulated sequences were truncated to 896 samples (the length of SC the real recordings, see below) and processed as the real data. These perfusion sequences were denoted ”inf” on the SNR axis (Figs. 3, 5). NU Perfusion Sequences with Noise. The in vivo perfusion sequences had about 130 samples before the contrast arrived. By merging several of these precontrast signal parts obtained from several recordings (converted to Ct [n]) and truncating MA the result, a signal of 896 samples with no contrast agent contribution and the noise level of real perfusion sequences was generated. Two such noise signals for the simulated perfusion sequences were gener- ED ated, one for the control mice and one for the treated mice. The noise signals were added to the corresponding almost noise free signals, to simulate the real perfusion sequences as closely as possible. PT To generate perfusion sequences with lower SNRs, the noise signals were multiplied with a constant factor higher than 1.0 before they were added to the AC CE almost noise free signals. The noise corrupted simulated signals were median filtered as the real signals (see below) before the blind deconvolution. The SNR for a single perfusion sequence was defined as the mean value of the perfusion sequence divided by the standard deviation of the noise [22]. Only the signal corresponding to the recorded part of the signal was used to estimate the SNR. The SNR for the multi-channel algorithms was defined as the mean of the SNR over all channels. 2.7. Animals, Imaging and Relaxation Rates All animal procedures were approved by the National Animal Research Authority (Oslo, Norway). All experiments were done in accordance with these protocols. Female C57Bl6 mice were anesthetized with Isoflurane in oxygen, and the body temperature was kept at 37o C throughout the experiment. One group of mice (n=8) received an injection of C48/80 (50 µg in 0.1 ml saline) 12 T ACCEPTED MANUSCRIPT RI P and one group of mice (n=8) received saline only. The data from these 16 mice have been reported and used for analysis also in [22, 23]. C48/80 acts as a mast cell degranulating agent and induces an inflammation due to the release of his- SC tamine [36]. Imaging was done using the FLASH sequence of a 7T Pharmascan (Bruker Biospin, Germany) and a mouse head volume coil. Precontrast T1(0) maps were computed based on precontrast T1-weighted NU recordings with fixed TR/TE = 15.0/2.5 ms and flip angles α = 5◦ , 10◦ , 15◦ and 20◦ , and 10 repetitions per flip angle. The linear form of the FLASHequation [48] was used. This approach allows computation of T1(0) values MA without neglecting the T2∗ effect. In the following steps, the T2∗ effect was neglected, which is the standard approach, valid for short TE, as used here. The k · ρ map (k is a spatially ED invariant constant dependent on the scanner and its gain setting, ρ is spin density) was computed using the FLASH-equation, S = k·ρ·sin(α)(1−E1 )/(1− cos(α)E1 ). Here, S is the voxel intensity in the precontrast recordings with PT α = 5◦ , E1 = e−T R/T 1(0) and T 1(0) is the map calculated in the previous step. DCE-MRI (dynamic T1-weighted 2D FLASH sequence, TR/TE = 15.0/2.5 AC CE ms, flip angle 25◦ , number of frames K = 1000, temporal resolution 1.079 s, slice thickness 1 mm, one axial slice through the left masseter muscle, one oblique sagittal slice along the thoracic aorta and at least one arteria carotis communis) was performed immediately after administration of the drug/saline. The conversion from signal intensity to R1 time sequences was done using the FLASH-equation, where the term k · ρ was known from the analysis of the multiple flip angle precontrast images. The contrast agent (gadodiamide, 0.1 mmol (kg BW)−1 in saline) was injected manually through the tail vein after acquisition of 130 frames. Care was taken to use as similar injection speed as possible for all mice. To increase the SNR to acceptable levels, average signals in manually drawn regions of interest were used. The average masseter muscle signal was derived from 40-80 voxels in the left and right masseter muscles. The average temporalis muscle signal was derived from 13-26 voxels in the left and right temporalis 13 T ACCEPTED MANUSCRIPT RI P muscles. The resulting two time sequences were median filtered (window length 11 samples). This filtering was selected based on previous experiments [22] as SC a compromise allowing sufficient noise reduction while not degrading the highfrequency content of the underlying signal. The fairly long window of the median filter is acceptable for tissue perfusion sequences of muscles due to their slow NU changes in shape. For arterial signals, the smoothing effect of the median filter affected mainly the shape of the peak (Fig. 4). However, this slight distortion is acceptable since the arterial signals are used only for generation of synthetic MA perfusion sequences. The real perfusion sequences were extrapolated to 13800 samples (3 hours) and zeropadded to 16384 samples (4 hours 33 min). This size was chosen as ED the lowest power of two (to allow use of fast Fourier transform) giving a signal which is longer than the assumed minimum washout of 3 hours (see section 2.1). PT 2.8. Evaluation Methods The simulated and real sequences were processed by single- and multi- AC CE channel blind deconvolution using the nonparametric and semi-parametric AIF. Multi-channel blind deconvolution was applied to the masseter and temporalis muscle perfusion sequences, i.e. two channels were used. For simulated data, AIF scaling was done with respect to the true reference AIF, i.e. an ideal case of accurate AIF scaling was simulated. For real data, AIF scaling was done according to Section 2.5. Simulations. To quantify the fit between a simulated perfusion sequence, Ct [n], and the corresponding estimated perfusion sequence, Ĉt [n], the approximation difference between these two sequences was computed. First, both Ct [n] and Ĉt [n] were normalized to have area under the curve 1.0. Then, the sum of absolute values of all sample based differences between the two sequences was computed, and taken as the approximation difference. To quantify the fit between a simulated true AIF and an estimated AIF, the overlap of the two sequences was computed. First, both the simulated and 14 T ACCEPTED MANUSCRIPT RI P estimated AIFs were normalized to have area under the curve 1.0. Then the absolute value difference, absdiff, between the two AIFs was computed as the mean of absolute values of all sample based differences. The percent overlap SC was defined as (1−absdiff/2) · 100%. To evaluate the accuracy of the perfusion-parameter estimation, the measure e was introduced as a mean relative error. It is defined as [22] 1 |Fb − Fˆb | |Tc − Tˆc | |E − Ê| |kep − k̂ep | ( + + + ) · 100%, 4 Fb Tc E kep NU e= (5) MA where the symbols with accent are the estimated values and the values without accent are the known reference values. The parameters Tc , E and kep are independent of AIF scaling. The parameter Fp depends on both the shape and the scaling of the estimated AIF. As an ideal (accurate) scaling is assumed in the ED processing of the simulated data, only the effect of the estimated AIF shape on Tc , E, kep and Fp is evaluated by Eq. (5). As a simplification, all four perfusion PT parameters are weighted equally in Eq. (5), supposing equal diagnostic value and impact on the precision and accuracy of the derived perfusion parameters (e.g. P S, ve ). AC CE The paired replicate data were compared using Wilcoxon’s two-sided distri- bution free signed rank test [49]. The significance level was 0.05. In-Vivo Data. To quantify the fit between an observed perfusion sequence, Ct [n], and the corresponding estimated perfusion sequence, Ĉt [n], the approximation difference between the two sequences was computed in the same way as described for the simulated perfusion sequences. To quantify the accuracy of the estimated AIFs, an indirect method was used because no ground truth AIF was available. The measured AIFs were not considered as a ground-truth because of acquisition artifacts and because reasonable AIF measurements were possible only for a limited number of animals (due to acquisition artifacts). Hence, the percent overlap, as defined for the simulated data, was computed between the AIF estimated from the masseter and temporalis muscles for each animal. This provided a comparison of 15 T ACCEPTED MANUSCRIPT RI P the nonparametric and semi-parametric AIF formulations for the single-channel algorithm. For multi-channel deconvolution this comparison was not available because it estimated one AIF common for both muscles/channels. To compare SC the single- and multi-channel algorithms, the overlaps between the nonparametric and semi-parametric AIF estimates was compared. NU 3. Results MA 3.1. Simulations Perfusion Sequences. The approximation difference (see methods for definition) was in general clearly lower for the nonparametric AIF than for the semiparametric AIF. This was observed for both the single- and multi-channel al- ED gorithms (Fig. 3). This shows that the nonparametric AIF has clearly more degrees of freedom to adapt to the signal noise. It also reflects the noise in PT the true AIFs (Fig. 4), and that the decaying phase of even noiseless true AIFs is not exactly modeled by three exponential functions. Examples of simulated perfusion sequences and their approximations are shown in Fig. 2 (the ripples AC CE in the simulated curves are due to the added noise and the subsequent median filtering of the simulated signals). The approximation difference was higher for the multi-channel algorithms than the single-channel algorithms (Fig. 3, note the different scale on the vertical axes). The reason for this difference was that compared to single-channel deconvolution, multi-channel deconvolution had less degrees of freedom due to the assumption of the same AIF for all channels. Hence, the single-channel estimates of the perfusion sequences could adapt to the given noise realization more closely. AIFs. For both, single- and multi-channel deconvolution, with no noise, the nonparametric AIF gave higher overlap with the true AIF than the semi-parametric AIF (Figs. 4, 5). However, when SNR decreased, the semi-parametric AIF had 16 T ACCEPTED MANUSCRIPT RI P consistently higher overlap. This clearly shows that semi-parametric AIF estimation is more robust with respect to noise (Figs. 4, 5). Comparing the multi- and single-channel algorithms, with no noise, the over- SC lap of the estimated AIFs with the true AIF was about the same. When SNR was about 10 or lower, the multi-channel algorithm gave better overlap of the estimated and the true AIF than the single-channel algorithm (Fig. 5). This was NU the case for both the nonparametric and semi-parametric AIFs in controls and treated mice. This could be explained by the stabilizing function of the prior assumption in multi-channel blind deconvolution that the AIFs of all channels MA were the same. This assumption took effect mostly in low SNR conditions. TRFs. The percent error of the perfusion parameters (Eq. (5)) corresponded ED to the observations obtained from overlaps of the estimated and reference AIFs described above. However, for the percent error of the perfusion parameters, the PT effect of randomness was higher than for the AIF-overlap measure. Hence, the percent error of the perfusion parameters was treated as a supportive measure and the corresponding tables are not shown. AC CE Comparing the nonparametric and semi-parametric AIF formulations, the percent errors of the perfusion parameters were comparable (mostly under 10% for SNR>15 and in the range between 4% and 57% for SNR<15). For high SNR (SN R >15), the nonparametric AIF led to clearly better perfusion-parameter estimates. For low SNR (SN R <15), the semi-parametric AIF approach was slightly (but non-significantly) better. Comparing single- and multi-channel deconvolution, the obtained percent errors of the perfusion parameters were comparable. The multi-channel algorithm gave clearly less accurate perfusion-parameter estimates for SNR>5 and slightly better (but not significantly) for low SNR (SN R < 5). 3.2. Real Perfusion Time Sequences SNR and Perfusion Sequences. As expected from the known effects of C48/80 on muscle perfusion, the SNR of the perfusion sequences of the treated muscles 17 T ACCEPTED MANUSCRIPT RI P was clearly higher than that of the corresponding control muscles (Table 2). Example perfusion sequences are shown if Fig. 6 (noise is manifested as ripples remaining after median filtering). The SNR of the perfusion sequences of the SC control masseter and temporalis muscles did not differ significantly, while the SNR in the treated mice was significantly higher for the masseter muscle than for the temporalis muscle (Table 2). NU The approximation differences of the observed and estimated perfusion sequences for the single-channel algorithm with a nonparametric and semi-parametric AIF (Table 2) were comparable to the corresponding approximation differences MA of the simulations with a similar SNR (Fig. 3). These approximation differences were smaller for the treated mice than for the control mice (Table 2). For the multi-channel algorithm, the approximation differences were slightly higher for ED the real data (Table 2) than for the corresponding simulated data (Fig. 3). The approximation differences of the observed and estimated perfusion sequences for the single-channel algorithm with a semi-parametric AIF were larger PT than the corresponding differences of the single-channel algorithm with a nonparametric AIF (Table 2, p = 0.008 for both muscles). This was the same AC CE observation as for the simulated data. It showed that the use of the nonparametric AIF allowed more degrees of freedom of the model for the perfusion sequences than the use of the semi-parametric AIF. Compared to single-channel deconvolution, the multi-channel algorithm gave much larger approximation differences of the observed and estimated perfusion sequences (Table 2). This indicated that the AIF of the two muscles differed due to dispersion. Possible different delays were corrected for automatically in the algorithm. The multi-channel algorithm with a semi-parametric AIF gave approximation differences of the observed and estimated perfusion sequences that were similar to those of the multi-channel algorithm with a nonparametric AIF (Table 2). Again, the most likely explanation for this observation was that the AIFs of the two muscles differed due to dispersion. Noise was only a minor contributing factor. 18 T ACCEPTED MANUSCRIPT RI P AIFs. The overlap of the AIFs estimated from the masseter and temporalis muscles of the same mouse (Table 3, columns 2, 3) was clearly higher for the semi-parametric AIFs than for the nonparametric AIFs. For in vivo control mice SC recordings, this corresponded well to the results obtained on simulated data for low SNRs (including the SNRs of in vivo control mice, Fig. 5, left). However, the in vivo recordings of the treated animals had SNRs in the high SNR range NU of the simulated data (Fig. 5, left) where the semi-parametric AIF gave worse results than the nonparametric AIF (see the Discussion section below). The overlap of the nonparametric and semi-parametric AIF estimates for MA the same animal (Table 3, columns 4–6) showed a slightly better consistency of the multi-channel algorithm than the single-channel algorithm. Compared to the simulation results (Fig. 5), this corresponds to the low SNR range where ED the multi-channel algorithm is more reliable. This was true even for the treated group where the SNR was in the range where the accuracy of the single- and multi-channel algorithms was comparable in the simulations (see the Discussion PT section below). For the semi-parametric algorithms, the parameters of the three exponential AC CE functions used in the AIF approximation varied to a large extent. This was the case for each group of mice for both semi-parametric algorithms and between the masseter and temporalis muscles of the same mouse for the single-channel semi-parametric algorithm. TRFs. Perfusion parameters were compared for all real perfusion sequences for the single- and multi-channel algorithms with the nonparametric and semiparametric AIF formulations. Except for two cases (out of 32 cases: 4 perfusion parameters x 2 muscle tissues x 2 AIF formulations x 2 animal groups), the perfusion parameters (Fb , Tc , E and kep ) did not differ significantly. 19 T ACCEPTED MANUSCRIPT RI P 4. Discussion For both the nonparametric and semi-parametric AIFs, the estimated AIFs of the single- and multi-channel algorithms had a good overlap with the true SC AIFs when the SNR of the perfusion sequences was reasonable. The nonparametric and semi-parametric AIFs led to comparable perfusion parameter esti- NU mates. Thus, modeling the decay phase of real AIFs of mice with three exponential decaying functions with the same delay was adequate for finding an almost noise free semi-parametric AIF specific for each mouse. MA In contrast, as done previously [41, 42], using only two exponential decaying functions with the same delay was not sufficient in the present setup. The AIF model with two exponentials might be sufficient when a simpler TRF model is used, as in [41, 42] (extended Tofts model). However, using a more complete ED TRF model (such as the aaJW model used here) imposes higher demands on the accuracy of the AIF model. PT As expected, a semi-parametric AIF led to a better noise robustness than the nonparametric AIF. This conclusion follows from the overlaps of the estimated and true AIFs for low SNR in the simulation experiments (Fig. 5) and from the AC CE consistency of the AIFs estimated from real recordings (Table 3, rows 2–3). In low SNR conditions, the semi-parametric AIF formulation provides additional prior information which stabilizes the estimation process. On the other hand, the nonparametric AIF estimates adjusted to the noise and diverged from the true AIF for the low SNR case. This is illustrated by the fact that the approximation differences of the estimated and simulated/measured perfusion sequences showed a better fit for the nonparametric AIF also for low SNR (Fig. 3, Table 2). For real data, the superiority of the semi-parametric AIF formulation was shown by the consistency of the AIFs estimated independently from two muscles. The semi-parametric AIF formulation gave better results than the nonparametric AIF (Table 3) even for high SNR conditions (treated mice, see Table 2) which was not in line with the simulations (Fig. 5). This observation suggests 20 T ACCEPTED MANUSCRIPT RI P the effect of an additional source of ”’noise” present in the real data and not accounted for in the simulations. This could be the deviation of the theoretical pharmacokinetic model and the real pharmacokinetics of the muscles. SC Another source of inaccuracy could be conversion from signal intensity to contrast agent concentration, as described in [50]. A solution would be to adjust the acquisition parameters according to the recommendations in [50]. Also, the NU nonlinear relation of the contrast agent concentration to the R1 relaxation rate due to the finite water exchange rate between blood plasma and tissue interstitial space was not considered here. Taking this effect into account would require MA extension of the pharmacokinetic model at the expense of two additional free parameters to be estimated [51]. It has not been shown, so far, whether such a model extension can lead to more accurate and precise blind deconvolution. ED The approximation of a function with three decaying exponential functions is an ill-conditioned problem, known for a long time [43]. The ill-conditioning is thought to be related to the fact that decaying exponential functions do not PT form an orthogonal basis. The observed large variation of the estimated AIF parameters stated above is likely to have the same origin. AC CE Comparing the single- and multi-channel algorithms for low SNR, multi- channel deconvolution lead to better accuracy of the AIF estimation for simulated data (Fig. 5). The reason was most probably that the assumption of the same AIF for all channels provided additional prior information. The comparison in terms of approximation differences of the estimated and simulated / measured perfusion sequences (Fig. 3, Table 2) showed the best fit for single-channel deconvolution for all SNR cases. This is reasonable since the single-channel deconvolution algorithm has more degrees of freedom to adjust to the actual noise realization than the multi-channel approach. For real data, an additional drawback of the multi-channel approach has to be taken into account, namely the dispersion differences of the AIFs in the different channels (not modeled in the simulations). The choice of the aaJW pharmacokinetic model for the TRF was based on previous work on perfusion in mouse masseter muscle [22, 23]. The choice of 21 T ACCEPTED MANUSCRIPT RI P a realistic TRF model for a given tissue is critical in obtaining a realistic AIF estimate, since the TRF model will have a strong influence on the time course of the estimated AIF [23]. SC Muscle tissue is present in most perfusion images of the body. Here, the aaJW model seems adequate for the TRF. Hence, a realistic blind estimate for the AIF in the region of interest can be obtained, and subsequently used in NU non-blind voxel-by-voxel perfusion analysis of the whole image. However, the challenge of choosing a realistic TRF model for tissues other than muscle in the image remains. MA The presented AIF scaling method was based on a known value of ve + vp in a reference tissue and on the assumption that the area under the AIF curve is the same for all animals in the study, due to the same contrast agent dose and ED the same kidney function. These are certainly simplifications. In [31], the present AIF scaling was compared to other two approaches. The first method was scaling to the area under the curve of a measured AIF (an PT arterial perfusion sequence derived from the DCE-MRI recording). The second method was scaling to the area under the curve of the tail of a measured AIF AC CE (the part of the curve after the peak of signal intensity where the saturation artifacts are reduced, used previously in [21, 52]). The AIF scaling method presented here gave clearly the most consistent results [31]. In [53, 52], the use of one venous sample at a late time (when arterial and venous concentrations are the same) was mentioned. Another option, mentioned in [21, 52], is to apply an additional scan after the DCE-MRI acquisition which would provide an accurate measurement of contrast agent concentration in a nearby artery or vein. These techniques need to be implemented and compared to decide on the best AIF scaling method for a given application. This work was limited to DCE-MRI in small animals, where the decay phase is monotonously decaying [54, 15, 42]. In humans, the decay phase of the AIF usually has a bimodal or three-modal time course. Several fully parametric models for these AIFs exist [21, 17, 55]. A future challenge is to incorporate such an AIF model in the presented framework to allow estimation of almost 22 RI P noise free AIFs for single voxel perfusion analysis in humans. T ACCEPTED MANUSCRIPT The present blind deconvolution method was computationally demanding. For one perfusion sequence, the current implementation (Linux C-shell scripts SC and C files, not optimized for speed) gave computation time of approx. 10 min for single-channel and 12 min for a two-channel deconvolution, on a PC with Intel Core 2 Duo P8400 processor (2.53 GHz, 725 3 MB L2 cache), 3.5 GB NU RAM, Xubuntu Linux Xfce 4.4. By a speed optimized reimplementation of the algorithm, a very large reduction in compuation time is expected. MA 5. Conclusion For quantitative preclinical DCE-MRI perfusion analysis on simulated and real data, a semi-parametric model of the AIF was introduced in the estimation ED of the AIF. The performance of this model combined with the aaJW pharmacokinetic model for the TRF was studied using single- and multi-channel PT Lucy-Richardson blind deconvolution expectation-maximization algorithms. The algorithm consists of two steps. In the first part, nonparametric estimates of the AIF and TRF are obtained. In the second part, the TRF is AC CE approximated using the aaJW pharmacokinetic model and the decaying part of the AIF is approximated by three decaying exponential functions. This leads to almost noise free semi-parametric AIFs. The tests on simulated and real data showed that use of semi-parametric AIFs led to better noise robustness of the AIF estimation for low SNR, when compared to the corresponding nonparametric AIF approach. The use of such blind deconvolution semi-parametric AIFs in the subsequent nonblind voxelwise deconvolution should lead to more accurate perfusion parameter maps than use of the corresponding nonparametric AIFs. Acknowledgment This study was supported by grants from the the Czech Science Foundation (grant no. 16-13830S) and by the Ministry of Education, Youth, and Sports of 23 T ACCEPTED MANUSCRIPT RI P the Czech Republic (project No. LO1212). The MR scanning was performed at the Molecular Imaging Center (MIC) at the Department of Biomedicine, SC University of Bergen. References NU [1] G. Brix, F. Kiessling, R. Lucht, S. Darai, K. Wasser, S. Delorme, J. Griebel, Microcirculation and microvasculature in breast tumors: pharmacokinetic analysis of dynamic MR image series., Magn Reson Med 52 (2) (2004) MA 420–9. [2] P. L. Choyke, A. J. Dwyer, M. V. Knopp, Functional tumor imaging with dynamic contrast-enhanced magnetic resonance imaging., JMRI 17 (5) ED (2003) 509–20. [3] I. S. Gribbestad, K. I. Gjesdal, G. Nilsen, S. Lundgren, M. H. B. Hjel- PT stuen, A. Jackson, An introduction to dynamic contrast-enhanced MRI in oncology, in: A. Jackson, D. L. Buckley, G. J. M. Parker (Eds.), Dynamic Contrast-Enhanced Magnetic Resonance Imaging in Oncology, Springer, AC CE Berlin, 2005. [4] C. K. Kuhl, H. H. Schild, Dynamic image interpretation of MRI of the breast., JMRI 12 (6) (2000) 965–74. [5] J. Barkhausen, P. Hunold, M. Jochims, J. F. Debatin, Imaging of myocardial perfusion with magnetic resonance, JMRI 19 (6) (2004) 750–757. [6] F. H. Epstein, J. F. London, D. C. Peters, L. M. Goncalves, K. Agyeman, J. Taylor, R. S. Balaban, A. E. Arai, Multislice first-pass cardiac perfusion MRI: validation in a model of myocardial infarction., Magn Reson Med 47 (3) (2002) 482–91. [7] P. Germain, G. Roul, J. Baruthio, C. Jahn, P. M. Coulbois, B. Dumitresco, J. L. Dietemann, P. Bareiss, A. Constantinesco, Myocardial flow reserve 24 T ACCEPTED MANUSCRIPT RI P parametric map, assessed by first-pass MRI compartmental analysis at the chronic stage of infarction., JMRI 13 (3) (2001) 352–60. [8] L. Hermoye, L. Annet, P. Lemmerling, F. Peeters, F. Jamar, P. Gianello, SC S. Van Huffel, B. E. Van Beers, Calculation of the renal perfusion and glomerular filtration rate from the renal impulse response obtained with NU MRI., Magn Reson Med 51 (5) (2004) 1017–25. [9] R. Materne, A. M. Smith, F. Peeters, J. P. Dehoux, A. Keyeux, Y. Horsmans, B. E. Van Beers, Assessment of hepatic perfusion parameters with MA dynamic MRI., Magn Reson Med 47 (1) (2002) 135–42. [10] O. Keunen, M. Johansson, A. Oudin, M. Sanzey, S. A. A. Rahim, F. Fack, F. Thorsen, T. Taxt, M. Bartos, R. Jiřı́k, H. Miletic, J. Wang, D. Stieber, ED L. Stuhr, I. Moen, C. B. Rygh, R. Bjerkvig, S. P. Niclou, Anti-VEGF treatment reduces blood supply and increases tumor cell invasion in glioblas- PT toma., P Natl Acad Sci USA 108 (9) (2011) 3749–54. [11] D. L. Buckley, G. J. M. Parker, Measuring contrast agent concentration AC CE in T1-weighted dynamic contrast-enhanced MRI, in: A. Jackson, D. L. Buckley, G. J. M. Parker (Eds.), Dynamic Contrast-Enhanced Magnetic Resonance Imaging in Oncology, Springer, Berlin, 2005. [12] G. J. M. Parker, D. L. Buckley, Tracer kinetic modeling for T1-weighted DCE-MRI, in: A. Jackson, D. L. Buckley, G. J. M. Parker (Eds.), Dynamic Contrast-Enhanced Magnetic Resonance Imaging in Oncology, Springer, Berlin, 2005. [13] F. Calamante, Arterial input function in perfusion MRI: A comprehensive review, Prog Nucl Magn Reson Spectrosc 74 (2013) 1–32. [14] C. S. Landis, X. Li, F. W. Telang, P. E. Molina, I. Palyka, G. Vetek, C. S. Springer, Equilibrium transcytolemmal water-exchange kinetics in skeletal muscle in vivo, Magn Reson Med 42 (3) (1999) 467–478. 25 T ACCEPTED MANUSCRIPT RI P [15] M. Heilmann, C. Walczak, J. Vautier, J.-L. Dimicoli, C. D. Thomas, M. Lupu, J. Mispelter, A. Volk, Simultaneous dynamic T1 and T2* measurement for AIF assessment combined with DCE MRI in a mouse tumor SC model., Magn Reson Mater Phy 20 (4) (2007) 193–203. [16] A. Garpebring, R. Wirestam, N. Ostlund, M. Karlsson, Effects of inflow and radiofrequency spoiling on the arterial input function in dynamic contrast- Med 65 (6) (2011) 1670–9. NU enhanced MRI: a combined phantom and simulation study., Magn Reson MA [17] G. J. M. Parker, C. Roberts, A. Macdonald, G. A. Buonaccorsi, S. Cheung, D. L. Buckley, A. Jackson, Y. Watson, K. Davies, G. C. Jayson, Experimentally-derived functional form for a population-averaged high- ED temporal-resolution arterial input function for dynamic contrast-enhanced MRI., Magn Reson Med 56 (5) (2006) 993–1000. PT [18] M. Heisen, X. Fan, J. Buurman, B. M. t. H. Romeny, Effects of reference tissue AIF derived from low temporal resolution DCE-MRI data on pharmacokinetic parameter estimation, in: ISMRM-ESMRMB Joint Annual AC CE Meeting 2010. - Sweden, Stockholm, Stockholm, 2010, p. 4802. [19] P. Campisi, K. Egiazarian, Blind Image Deconvolution, CRC Press, New York, NY, 2007. [20] C. Yang, G. S. Karczmar, M. Medved, W. M. Stadler, Multiple reference tissue method for contrast agent arterial input function estimation., Magn Reson Med 58 (6) (2007) 1266–75. [21] J. U. Fluckiger, M. C. Schabel, E. V. R. Dibella, Model-based blind estimation of kinetic parameters in dynamic contrast enhanced (DCE)-MRI., Magn Reson Med 62 (6) (2009) 1477–86. [22] T. Taxt, R. Jiřı́k, C. B. Rygh, R. Grüner, M. Bartos, E. Andersen, F.R. Curry, R. K. Reed, Single-channel blind estimation of arterial input 26 T ACCEPTED MANUSCRIPT RI P function and tissue impulse response in DCE-MRI., IEEE T Bio-Med Eng 59 (4) (2012) 1012–21. [23] T. Taxt, T. Pavlin, R. K. Reed, F.-R. Curry, E. Andersen, R. Jiřı́k, Using SC Single-Channel Blind Deconvolution to Choose the Most Realistic Pharmacokinetic Model in Dynamic Contrast-Enhanced MR Imaging, Appl Magn NU Reson 46 (6) (2015) 643–659. [24] D. Y. Riabkov, E. V. R. D. Bella, Estimation of Kinetic Parameters Without Input Functions : Analysis of Three Methods for Multichannel Blind MA Identification, IEEE T Bio-Med Eng 49 (11) (2002) 1318–1327. [25] M. C. Schabel, J. U. Fluckiger, E. V. R. DiBella, A model-constrained Monte Carlo method for blind arterial input function estimation in dynamic ED contrast-enhanced MRI: I. Simulations, Phys Med Biol 55 (16) (2010) 4783– 4806. PT [26] J. U. Fluckiger, M. C. Schabel, E. V. DiBella, Toward local arterial input functions in dynamic contrast-enhanced MRI, JMRI 32 (4) (2010) 924–934. AC CE [27] R. Jiřı́k, K. Souček, M. Mézl, M. Bartoš, E. Dražanová, F. Dráfi, L. Grossová, J. Kratochvı́la, O. Macı́ček, K. Nylund, A. Hampl, O. H. Gilja, T. Taxt, Z. Starčuk, Blind deconvolution in dynamic contrast-enhanced MRI and ultrasound., in: Conf Proc IEEE Eng Med Biol Soc., Vol. 2014, 2014, pp. 4276–9. [28] I. Jacobs, G. J. Strijkers, H. M. Keizer, H. M. Janssen, K. Nicolay, M. C. Schabel, A novel approach to tracer-kinetic modeling for (macromolecular) dynamic contrast-enhanced MRI., Magn Reson Med 75 (3) (2016) 1142– 1153. [29] K. S. St. Lawrence, T.-Y. Lee, An Adiabatic Approximation to the Tissue Homogeneity Model for Water Exchange in the Brain: I. Theoretical Derivation, J Cerebr Blood F Met 18 (12) (1998) 1365–1377. 27 T ACCEPTED MANUSCRIPT RI P [30] T. S. Koh, V. Zeman, J. Darko, T. Y. Lee, M. F. Milosevic, M. Haider, P. Warde, I. W. Yeung, The inclusion of capillary distribution in the adiabatic tissue homogeneity model of blood flow., Phys Med Biol 46 (5) (2001) SC 1519–1538. [31] J. Kratochvı́la, R. Jiřı́k, M. Bartoš, M. Standara, Z. Starčuk, T. Taxt, Distributed capillary adiabatic tissue homogeneity model in parametric multi- NU channel blind AIF estimation using DCE-MRI, Magn Reson Med 75 (3) (2016) 1355–1365. MA [32] M. C. Schabel, A unified impulse response model for DCE-MRI, Magn Reson Med 68 (5) (2012) 1632–1646. [33] J. Boesiger, M. Tsai, M. Maurer, M. Yamaguchi, L. F. Brown, K. P. Claffey, ED H. F. Dvorak, S. J. Galli, Mast cells can secrete vascular permeability factor/ vascular endothelial cell growth factor and exhibit enhanced release PT after immunoglobulin E-dependent upregulation of fc epsilon receptor I expression., J Exp Med 188 (6) (1998) 1135–45. AC CE [34] M. E. Koller, K. Woie, R. K. Reed, Increased negativity of interstitial fluid pressure in rat trachea after mast cell degranulation., J Appl Physiol 74 (5) (1993) 2135–9. [35] A. Rosén, P. Strandberg, B. Uvnäs, Vascular Effects of Moderate Doses of Histamine and Histamine Liberator (Compound 48/80)., Acta Pharmacol Tox 13 (3) (2009) 267–276. [36] A. M. Rothschild, Mechanisms of histamine release by compound 48-80., Brit J Pharmacol 38 (1) (1970) 253–62. [37] J. L. Starck, E. Pantin, F. Murtagh, Deconvolution in Astronomy: A Review, Publ Astron Soc Pac 114 (800) (2002) 1051–1069. [38] J. A. Johnson, T. A. Wilson, A model for capillary exchange., Am J Phys 210 (6) (1966) 1299–303. 28 T ACCEPTED MANUSCRIPT RI P [39] W. C. Sangren, C. W. Sheppard, A mathematical derivation of the exchange of a labeled substance between a liquid flowing in a vessel and an external compartment, Bull Math Biophys 15 (4) (1953) 387–394. SC [40] P. S. Tofts, G. Brix, D. L. Buckley, J. L. Evelhoch, E. Henderson, M. V. Knopp, H. B. Larsson, T. Y. Lee, N. A. Mayr, G. J. Parker, R. E. Port, J. Taylor, R. M. Weisskoff, Estimating kinetic parameters from dynamic NU contrast-enhanced T1-weighted MRI of a diffusable tracer: Standardized quantities and symbols., JMRI 10 (3) (1999) 223–32. MA [41] D. M. McGrath, D. P. Bradley, J. L. Tessier, T. Lacey, C. J. Taylor, G. J. M. Parker, Comparison of model-based arterial input functions for dynamic contrast-enhanced MRI in tumor bearing rats., Magn Reson Med 61 (5) ED (2009) 1173–84. [42] D. K. Ragan, S. Y. Lai, J. A. Bankson, Fast, reproducible measurement of PT the vascular input function in mice using constrained reconstruction and cardiac sampling., NMR Biomed 24 (4) (2011) 373–84. AC CE [43] W. Wiscombe, J. Evans, Exponential-sum fitting of radiative transmission functions, J Comput Phys 24 (4) (1977) 416–444. [44] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical Recipes in C, Cambridge University Press, New York, NY, 1992. [45] R. K. Reed, H. Wiig, Compliance of the interstitial space in rats. I. Studies on hindlimb skeletal muscle., Acta Physiol Scand 113 (3) (1981) 297–305. [46] V. Goh, S. Halligan, J.-A. Hugill, C. I. Bartram, Quantitative assessment of tissue perfusion using MDCT: comparison of colorectal cancer and skeletal muscle measurement reproducibility., Am J Roentgenol 187 (1) (2006) 164– 9. [47] S. M. Galbraith, M. A. Lodge, N. J. Taylor, G. J. S. Rustin, S. Bentzen, J. J. Stirling, A. R. Padhani, Reproducibility of dynamic contrast-enhanced 29 T ACCEPTED MANUSCRIPT RI P MRI in human muscle and tumours: comparison of quantitative and semiquantitative analysis., NMR Biomed 15 (2) (2002) 132–42. [48] H.-L. M. Cheng, G. A. Wright, Rapid high-resolution T1 mapping by vari- SC able flip angles: accurate and precise measurements in the presence of radiofrequency field inhomogeneity., Magn Reson Med 55 (3) (2006) 566–74. NU [49] M. Hollander, D. A. Wolfe, Nonparametric Statistical Methods, John Wiley and Sons Inc., New York, NY, 1973. MA [50] M. C. Schabel, D. L. Parker, Uncertainty and bias in contrast concentration measurements using spoiled gradient echo pulse sequences, Phys Med Biol 53 (9) (2008) 2345–2373. ED [51] J. Zhang, S. Kim, Uncertainty in MR tracer kinetic parameters and water exchange rates estimated from T1-weighted dynamic contrast enhanced PT MRI, Magn Reson Med 72 (2) (2014) 534–545. [52] M. C. Schabel, E. V. R. Dibella, R. L. Jensen, K. L. Salzman, A ModelConstrained Monte Carlo Method for Blind Arterial Input Function Es- AC CE timation in Dynamic Contrast-Enhanced MRI: II) In Vivo Results, Phys Med Biol 55 (16) (2010) 4807–4823. [53] D. Y. Riabkov, E. V. R. Di Bella, Blind identification of the kinetic parameters in three-compartment models, Phys Med Biol 49 (5) (2004) 639–664. [54] D. Checkley, J. J. L. Tessier, S. R. Wedge, M. Dukes, J. Kendrew, B. Curry, B. Middleton, J. C. Waterton, Dynamic contrast-enhanced MRI of vascular changes induced by the VEGF-signalling inhibitor ZD4190 in human tumour xenografts., Magn Reson Med 21 (5) (2003) 475–82. [55] C. Yang, G. S. Karczmar, M. Medved, A. Oto, M. Zamora, W. M. Stadler, Reproducibility assessment of a multiple reference tissue method for quantitative dynamic contrast enhanced-MRI analysis., Magn Reson Med 61 (4) (2009) 851–9. 30 T ACCEPTED MANUSCRIPT RI P Tables Table 1: Parameters of the simulated TRFs.. tissue Fb Tc E [ml/min/100ml] [s] [-] controls, masseter m. 6.9 11.9 66.1 0.567 treated, masseter m. 29.4 6.5 51.3 0.661 controls, temporalis m. 4.8 7.6 59.1 0.472 treated, temporalis m. 17.4 [1/min] SC NU MA ED PT AC CE 31 kep 7.6 56.5 0.577 Approximation differences of the perfusion sequences (see Sec. 2.8 for definition). RI P Table 2: T ACCEPTED MANUSCRIPT Real data. Each difference multiplied by 100. Single-channel algorithm: column 2-7. Multichannel algorithm: column 8-10. Nonp: Nonparametric. Semi: Semi-Parametric. mouse SNR nonp temporalis m. semi SNR nonp controls 14.3 0.22 0.60 16.3 22 14.6 0.28 0.72 14.6 23 20.9 0.21 0.59 15.2 24 13.6 0.28 0.82 20.4 25 21.3 0.32 3.79 21.5 26 30.0 0.33 1.58 28 13.6 0.35 mean 18.3 0.28 SD 6.1 0.05 0.26 0.66 NU 19 semi mean of both muscles SNR nonp semi 15.3 2.67 2.54 SC masseter m. 0.71 14.6 6.78 6.38 0.19 0.38 18.1 6.50 6.23 0.31 1.33 17.0 5.92 5.99 0.18 0.72 21.4 14.80 14.65 24.2 0.36 1.33 27.1 2.13 2.86 0.73 13.5 0.41 1.12 13.6 6.24 6.17 1.26 18.0 0.25 0.89 18.2 6.43 6.40 4.1 0.12 0.37 4.7 4.15 3.99 ED MA 0.26 1.17 treated 35.1 0.18 1.14 30.6 0.21 0.75 32.9 5.75 5.77 6 86.2 0.07 0.48 34.1 0.11 0.58 60.2 1.92 1.99 7 83.8 0.07 0.27 42.7 0.08 0.47 63.3 8.04 8.18 AC CE PT 5 8 49.7 0.17 1.12 53.5 0.09 0.54 51.6 2.02 1.99 10 62.8 0.10 0.55 54.2 0.12 0.58 58.5 0.43 0.67 12 49.3 0.10 0.66 44.4 0.12 0.65 46.9 2.42 2.46 13 53.0 0.09 0.56 34.1 0.15 0.85 43.6 6.52 3.83 mean 60.0 0.11 0.68 41.9 0.13 0.63 51.0 3.87 3.55 SD 18.3 0.05 0.33 9.5 0.04 0.13 10.7 2.86 2.61 32 Overlap of AIFs in percent. Real data. Single-channel algorithm (single), overlap RI P Table 3: T ACCEPTED MANUSCRIPT of AIFs estimated from the masseter and temporalis muscles using nonparametric (nonp, column 2) and semi-parametric (semi, column 3) AIFs. Overlap of AIFs estimated using the nonparametric and semi-parametric AIFs, single-channel algorithm applied to the masseter SC (column 4) and temporalis (column 5) muscles and multi-channel algorithm (column 6). mass versus nonp versus semi multi 19 95.1 NU temp mouse single single single single nonp semi mass temp 94.8 92.0 94.7 94.4 22 94.1 95.6 96.4 94.4 96.5 23 93.1 96.1 96.6 94.8 96.0 24 91.9 96.8 95.7 95.2 96.2 ED MA controls 90.8 95.9 93.4 94.4 96.0 26 92.3 95.9 94.3 93.6 95.5 28 91.4 96.9 92.9 92.9 95.0 mean 92.7 96.0 94.5 94.3 95.7 SD 1.5 0.7 1.8 0.8 0.7 AC CE PT 25 treated 5 94.6 95.7 94.5 96.3 96.3 6 96.0 98.6 97.4 96.5 97.4 7 90.4 96.2 95.7 93.4 97.5 8 95.6 96.1 94.9 96.1 96.8 10 97.5 97.3 97.5 97.3 98.0 12 95.7 97.4 96.6 96.1 96.9 13 96.9 98.8 97.5 97.3 98.0 mean 95.3 97.2 96.3 96.1 97.3 SD 2.3 1.2 1.3 1.3 0.6 33 T ACCEPTED MANUSCRIPT RI P Figure Captions SC Figure 1. Flowchart of the single-channel blind deconvolution algorithm. NU Figure 2. Perfusion sequences. Two simulated sequences of the temporalis muscle and the corresponding estimated sequences with a semi-parametric AIF. Control-mouse conditions (SNR = 8.0) (a) and treated-mouse conditions (SNR MA = 23.3) (b). Figure 3. Approximation differences of the simulated perfusion sequences (see ED Sec. 2.8 for definition). Each difference multiplied by 100. Nonp: nonparametric estimate of the AIF. Semi: semi-parametric estimate of the AIF. Horizontal axis: PT ”inf” – noiseless sequences, ”inf-full” – noiseless sequences with no truncation Figure 4. Reference and estimated AIFs. Single-channel algorithm. The two AC CE reference AIFs used in the simulations (solid lines) and the corresponding semiparametric AIFs estimated by blind deconvolution of the simulated time sequences (dotted and dashed lines). (a) – The reference AIF for control mice and the AIFs estimated from masseter (Mass. AIF, SNR=22.7) and temporalis (Temp. AIF, SNR=16.0) muscles. (b) – The reference AIF for treated mice and the AIFs estimated from masseter (Mass. AIF, SNR=63.5) and temporalis (Temp. AIF, SNR=46.6) muscles. SNR values correspond to those of the real recordings. Insets show magnified AIF segments. Figure 5. Percent overlap between true and estimated AIFs (see Sec. 2.8 for definition). Nonp: nonparametric estimate of AIF. Semi: semi-parametric estimate of AIF. Horizontal axis: ’inf’ – noiseless sequences, ’inf-full’ – noiseless sequences with no truncation. 34 T ACCEPTED MANUSCRIPT RI P Figure 6. Observed and approximated perfusion sequences, real data, singlechannel algorithm, semi-parametric AIF. Temporalis muscle of control mouse SC 24 (a) (SNR = 20.4) and treated mouse 6 (b) (SNR = 34.1). Figure 7. Estimated AIFs of control mouse 24. Single-channel algorithm. Nonparametric and semi-parametric AIFs. Masseter muscle with SNR = 13.6 (a). NU Temporalis muscle with SNR = 20.4 (b). MA Figure 8. Estimated AIFs of treated mouse 6. Single-channel algorithm. Nonparametric and semi-parametric AIFs. Masseter muscle with SNR = 86.2 (a). AC CE PT ED Temporalis muscle with SNR = 34.1 (b). 35 T ACCEPTED MANUSCRIPT Preprocessing Initialization of AIF, TRF, iteration index s=1 S>200 - Parametric approximation of nonparametric TRF + Nonparametric updating of AIF Semi-Parametric approximation of nonparametric AIF MA s=s+1 + NU S>200 - SC Nonparametric updating of TRF RI P Figures S>290 Scaling of AIF and TRF AC CE PT ED Figure 1: Flowchart of the single-channel blind deconvolution algorithm. 36 T ACCEPTED MANUSCRIPT (a) RI P 1 Simulated Estimated 0.8 SC 0.6 0.4 0 0 5 NU 0.2 10 time [min] 15 (b) MA 1 0.8 0.6 Simulated Estimated ED 0.4 0.2 0 PT 0 AC CE Figure 2: Perfusion sequences. 5 10 time [min] 15 Two simulated sequences of the temporalis muscle and the corresponding estimated sequences with a semi-parametric AIF. Control-mouse conditions (SNR = 8.0) (a) and treated-mouse conditions (SNR = 23.3) (b). 37 T ACCEPTED MANUSCRIPT RI P 12 10 8 SC Approximation difference [a.u.] 3 1 20 40 Inf Inf−full ED Masseter, nonp Masseter, semi Temporalis, nonp Temporalis, semi PT 1 0.5 0 Figure 3: 20 40 60 80 SNR [−] Inf Inf−full Approximation difference [a.u.] Treated, single−channel dec. 1.5 0 60 80 SNR [−] MA 0 NU 2 0 Controls, multi−channel dec. Masseter, nonp Masseter, semi Temporalis, nonp Temporalis, semi 4 AC CE Approximation difference [a.u.] Approximation difference [a.u.] Controls, single−channel dec. Masseter, nonp Masseter, semi Temporalis, nonp Temporalis, semi 6 4 2 0 0 3.5 3 2.5 2 1.5 1 0.5 0 20 40 60 80 SNR [−] Treated, multi−channel dec. Masseter, nonp Masseter, semi Temporalis, nonp Temporalis, semi 0 20 40 60 80 SNR [−] Approximation differences of the simulated perfusion sequences (see Sec. 2.8 for definition). Each difference multiplied by 100. Nonp: nonparametric estimate of the AIF. Semi: semi-parametric estimate of the AIF. Horizontal axis: ”inf” – noiseless sequences, ”inf-full” – noiseless sequences with no truncation 38 Inf Inf−full Inf Inf−full T ACCEPTED MANUSCRIPT (a) 0.02 High SNR 0.01 0 5 0.5 1 1.5 10 time [min] 2 NU 0 SC 0.02 0.01 0 Low SNR RI P Reference 15 Low SNR High SNR MA (b) Reference 0.04 0.04 0.03 0.01 0 5 0 0.5 10 time [min] AC CE PT 0 0.02 ED 0.02 1 1.5 15 Figure 4: Reference and estimated AIFs. Single-channel algorithm. The two reference AIFs used in the simulations (solid lines) and the corresponding semi-parametric AIFs estimated by blind deconvolution of the simulated time sequences (dotted and dashed lines). (a) – The reference AIF for control mice and the AIFs estimated from masseter (Mass. AIF, SNR=22.7) and temporalis (Temp. AIF, SNR=16.0) muscles. (b) – The reference AIF for treated mice and the AIFs estimated from masseter (Mass. AIF, SNR=63.5) and temporalis (Temp. AIF, SNR=46.6) muscles. SNR values correspond to those of the real recordings. Insets show magnified AIF segments. 39 Controls, single−channel dec. SC 95 90 NU Masseter, nonp Masseter, semi Temporalis, nonp Temporalis, semi 85 0 20 AC CE 90 0 20 40 Masseter, nonp Masseter, semi Temporalis, nonp Temporalis, semi 60 80 SNR [−] Inf Inf−full Controls, multi−channel dec. 95 90 85 80 Inf Inf−full ED PT 95 Figure 5: 60 80 SNR [−] Treated, single−channel dec. 100 85 40 Nonp Semi 0 20 40 Inf Inf−full 95 90 85 Nonp Semi 0 20 40 Percent overlap between true and estimated AIFs (see Sec. 2.8 for definition). Nonp: nonparametric estimate of AIF. Semi: semi-parametric estimate of AIF. Horizontal axis: ’inf’ – noiseless sequences, ’inf-full’ – noiseless sequences with no truncation. 40 60 80 SNR [−] Treated, multi−channel dec. 100 Aif overlap [%] 80 Aif overlap [%] RI P 100 MA Aif overlap [%] 100 Aif overlap [%] T ACCEPTED MANUSCRIPT 60 80 SNR [−] Inf Inf−full T ACCEPTED MANUSCRIPT (a) −3 1.5 RI P x 10 SC 1 0 5 10 time [min] Measured Estimated 15 MA 0 NU 0.5 (b) ED x 10−3 1.5 PT 1 0.5 AC CE 0 0 Measured Estimated 5 10 time [min] 15 Figure 6: Observed and approximated perfusion sequences, real data, single-channel algorithm, semi-parametric AIF. Temporalis muscle of control mouse 24 (a) (SNR = 20.4) and treated mouse 6 (b) (SNR = 34.1). 41 T ACCEPTED MANUSCRIPT (a) −3 RI P 6 x 10 Nonparametric Semi−parametric 2 0 5 10 time [min] 15 MA 0 NU SC 4 (b) Nonparametric Semi−parametric ED x 10−3 6 PT 4 AC CE 2 0 0 5 10 time [min] 15 Figure 7: Estimated AIFs of control mouse 24. Single-channel algorithm. Nonparametric and semi-parametric AIFs. Masseter muscle with SNR = 13.6 (a). Temporalis muscle with SNR = 20.4 (b). 42 T ACCEPTED MANUSCRIPT (a) −3 RI P x 10 Nonparametric Semi−parametric SC 6 4 5 0 10 time [min] MA 0 NU 2 15 0.01 AC CE PT 0.005 ED (b) 0 0 5 Nonparametric Semi−parametric 10 time [min] 15 Figure 8: Estimated AIFs of treated mouse 6. Single-channel algorithm. Nonparametric and semi-parametric AIFs. Masseter muscle with SNR = 86.2 (a). Temporalis muscle with SNR = 34.1 (b). 43

1/--страниц