{'fraction_non_alphanumeric': 0.033791005508509414, 'fraction_numerical': 0.012853188632190524, 'mean_word_length': 4.344171060339777, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bolus arrival time estimation in dynamic contrast-enhanced magnetic resonance imaging of small animals based on spline models

Dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) is used to quantify perfusion and vascular permeability. In most cases a bolus arrival time (BAT) delay exists between the arterial input function (AIF) and the contrast agent arrival in the tissue of interest which needs to be estimated. Existing methods for BAT estimation are tailored to tissue concentration curves which have a fast upslope to the peak as frequently observed in patient data. However, they may give poor results for curves that do not have this characteristic shape such as tissue concentration curves of small animals. In this paper, we propose a novel method for BAT estimation of signals that do not have a fast upslope to their peak. The model is based on splines which are able to adapt to a large variety of concentration curves. Furthermore, the method estimates BATs on a continuous time scale. All relevant model parameters are automatically determined by generalized cross validation. We use simulated concentration curves of small animal and patient settings to assess the accuracy and robustness of our approach. The proposed method outperforms a state-of-the-art method for small animal data and it gives competitive results for patient data. Finally, it is tested on in vivo acquired rat data where accuracy of BAT estimation was also improved upon the state-of-the-art method. The results indicate that the proposed method is suitable for accurate BAT estimation of DCE-MRI data, especially for small animals.


Introduction
Dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) is used to study perfusion and vascular permeability in tissues such as tumor vasculature (Sourbron and Buckley, 2012). This makes DCE-MRI especially interesting for the characterization of the vascular status in clinical studies (Padhani and Miles, 2010) and in preclinical studies employing small animals (Glowa et al., 2017;Hectors et al., 2018;Keunen et al., 2011;Schreurs et al., 2017).
DCE-MRI is based on acquiring MR-images at high temporal resolution after intravenously administering a contrast agent (CA) bolus. Concentration time curves describing the CA accumulation are extracted from the tissue of interest. Quantitative analysis of these tissue concentration curves (TCCs) is often performed by fitting pharmacokinetic models to them (Sourbron and Buckley, 2012). Many of these models require an arterial input function (AIF) describing the arrival of the CA in an artery close to the tissue of interest (Sourbron, 2010;Tofts et al., 1999). A common non-invasive technique for patient individual AIF extraction is to derive it directly from the acquired DCE-MRI images from a vessel close to the tissue of interest. Due to limited image resolution and partial volume effects it is not always possible to choose an artery in the direct vicinity of the tissue. Hence, the AIF is often extracted from a larger artery further upstream. However, due to the spatial distance, the CA bolus arrives later at the tissue of interest than suggested by the upstream measured AIF. Various studies have shown that neglecting this temporal delay decreases accuracy of the pharmacokinetic parameter estimates (Koh et al., 2011;Mehrtash et al., 2016;Nadav et al., 2017).
There are two standard approaches to determine the delay time: The first approach includes the time delay as an additional free parameter in the fitting process of the pharmacokinetic model (Kershaw and Buckley, 2006;Meyer, 1989). The second type estimates the bolus arrival times (BATs) separately for the TCC and for the AIF; that is, the time points where the concentration curves start to rise. The time delay is the difference between the BAT of the TCC and that of the AIF. The advantage of this second type is that it does not require to include the time delay as an additional fit parameter which is likely to decrease fit stability. A comparative study by Kershaw and Buckley (2006) showed that the second type gives more accurate results than the first one.
The gold standard for BAT estimation is the method of Cheong et al. (2003) which is based on fitting a linear or convex quadratic function to the sample points in the neighborhood of the BAT. An extension was proposed by Singh et al. (2009) which gives comparable results. These existing methods were mainly developed for patient data which often exhibit fast upslopes (i.e. the peak of the signal is close to the BAT). However, applying the method of Cheong et al. (2003) to TCCs which do not exhibit this characteristic slope, e.g. TCCs of small animals or weakly perfused tissues, may lead to poor results. Another limitation is that estimated BATs have to coincide with one of the acquisition time points. This can lead to inaccurate estimates, in particular when working with low temporal resolution data.
Even though many DCE-MRI studies are performed in a preclinical environment employing small animals (Glowa et al., 2017;Hectors et al., 2018;Keunen et al., 2011;Schreurs et al., 2017), to our knowledge, there is no suitable approach for BAT estimation of DCE-MRI data of rats and mice. TCCs tend to vary dramatically depending on the underlying physiology and vascular architecture of the respective tissue ranging from steeply increasing and subsequently rapidly decreasing TCCs in e.g. intact vessels, to monotonously and shallowly increasing TCCs in necrotic areas. Hence, a more flexible model is required, able to adapt to various kinds of CA uptake curves.
In this work we propose a novel method for BAT estimation of DCE-MRI TCCs that do not have a fast upslope -as often observed in small animal data and patient data of tissues without a significant intravascular fraction. The method is based on approximation of concentration curves by splines which allow to adapt to a large variety of curve shapes after the BAT. The BAT itself and the involved secondary model parameters are estimated automatically by generalized cross validation. Furthermore, the method is able to estimate BATs that do not coincide with one of the sampling time points. An extensive study using simulated data and in vivo rat data shows that the proposed method gives more accurate BAT estimations than the state-of-the art method.

The proposed method
2.1.1 TCCs and BATs DCE-MRI measurements are performed by acquiring a series of images at time points t n (n = 1, 2, . . . , N ). After a number of images, the CA is given as a bolus. In this paper, it is assumed that the MR-signal intensities are already converted into concentrations of the CA. Thus, the given data c n are noisy samples of a true TCC C(t): for n = 1, . . . , N.
Here, the n represent error terms which are assumed to be independent and identically distributed Gaussian distributed random variables with zero mean and unknown variance. For simplicity, we focus on the often used case of uniformly sampled time points, i.e., t n = (n−1)∆t, noting that the proposed method could also be extended to the more general case employing unevenly sampled time points. Our goal is to estimate the BAT t BAT from the noisy data c n which is the time point where the true TCC, C(t), starts to rise. Note that t BAT does not necessarily coincide with one of the sample time points t n .
The BAT naturally divides a TCC into two parts. Prior to bolus arrival, the TCC is constantvalued at zero (unless previous boluses were injected), whereas after the BAT, it rises and adopts different shapes depending on the underlying physiology. A common type of TCCs reaches their maximum almost immediately after the bolus arrival which is exploited in the method of Cheong et al. (2003). But often the time between the BAT and the curve maximum is much longer, e.g. in small animal data or patient data of tissues without a significant intravascular fraction. Thus, building on the above assumption of a steep initial signal slope can lead to poor estimates of the BAT. In the following, we present a non-parametric approximation model based on smoothing splines that allows for adaption to various curve shapes,in particular for TCCs without fast upslopes. We note that the proposed method is also applicable to the raw MR data before conversion to concentration.

Smoothing splines
Smoothing splines are a classical signal estimation method which have been applied to various biomedical problems (Craven and Wahba, 1978;Unser, 2002). We give a brief description here (for further details, the reader is referred to the references De Boor (1978); Reinsch (1967);Unser (1999); Wahba (1990)). First, recall that an interpolating spline of order k is a function u that passes through each data point while minimizing its smoothness costs, as represented by the square integral of its k-th derivative  (2). The given noisy TCC data (blue crosses) are approximated by a continuous piecewise defined function u * . It is a constant function before the BAT t BAT and afterwards a smoothing spline which is able to adapt to the shape of the curve by adjusting its stiffness parameters α and k. An estimate for the BAT, t BAT , is determined automatically along with suitable spline stiffness parameters α and k by generalized cross validation. Note that t BAT does not need to coincide with one of the sample times t 1 , . . . , t N .
Instead of interpolating, a smoothing spline approximates the data by optimizing a tradeoff between smoothness and fidelity to the data, as represented by the residual sum of squares. More precisely, a smoothing spline of order k ∈ N is the solution of the following optimization problem: Find a sufficiently smooth function u on [t 1 , t N ] that minimizes the cost function where u (k) denotes the k-th derivative of the function u. The parameter α > 0 adjusts the relative weight of data fidelity and smoothness: the larger the value of α, the smoother the resulting spline. To understand the role of the order k, it is instructive to look at the solutions for very large α: For α → ∞, the solution of (1) yields an ordinary polynomial of degree k − 1 on [t 1 , t N ]. This is because the smoothness cost (the second term in (1)) is equal to zero exactly for polynomials of maximum degree k − 1. Thus, for α → ∞, the minimizer of (1) coincides with the least squares approximation by a polynomial of degree k − 1.

Proposed approximation model
The above discussion on the characteristics of a concentration curve motivates to approximate a TCC by a continuous function that is constant on the interval [t 1 , t BAT ] and that is a smoothing spline on [t BAT , t N ]. A graphical illustration is given in Figure 1. Thus, the approximating function u * is defined by the minimizing condition subject to u is constant on [t 1 , t BAT ]. (2) As the approximation u * is not required to be smooth at the joint it typically has a kink at t BAT .
Since u is constant on [t 1 , t BAT ], problem (2) can be reduced tõ where N is the index of the last time frame before t BAT ; that is, N = max{n : t n ≤ t BAT }.
There is a simple relation between the solutions of the constrained variant (2) and that of the unconstrained variant (3): Our principal interest is to estimate the parameter t BAT from the concentration curve c. Alongside the estimation of t BAT , the spline parameters α and k have to be determined as well because these are needed for adapting the model to the curve shape after the BAT. All three parameters are estimated automatically as explained in detail in Section 2.1.5.

Discretization and numerical computation of the approximating function
To numerically compute the approximating functions, we discretize the optimization problem (3). First, a finite difference approximation of the k-th derivative u (k) is employed. A natural choice for the set of discretization nodes q are the BAT and the subsequently sampled time points, i.e., q = (t BAT , t N +1 , t N +2 , . . . , t N ). Letṽ be the vector that contains the ordinates of u at these discretization nodes; that is,ṽ i =ũ(q i ) for i = 1, . . . , I with I = N − N + 1. A first order finite difference approximation of the k-th derivative is given bỹ where the ω (k,t BAT ) ij ∈ R are weights which we specify in the following. First note that the weights depend on t BAT because the first discretization node coincides with t BAT . Determining the weights ω (k,t BAT ) ij is a standard procedure; see e.g. Fornberg (1996). Further, we approximate the integral in (3) by a Riemannian sum. In summary, we obtain the discretized problem Just asṽ * is a discrete approximation to the continuous functionũ * , the expanded vector v * = (ṽ * 1 , . . . ,ṽ * 1 ,ṽ * 2 , . . . ,ṽ * I ) (of length N ) is a discrete approximation to u * .
Equation (5) is a least squares problem inṽ. It is convenient to reformulate it in matrix notation,ṽ where v 2 denotes the squared Euclidean norm, v 2 = n v 2 n , and the involved matrices are as follows. W is a diagonal matrix with the entries where 1 N −1 = (1, . . . , 1) T denotes a vector of (N − 1) ones, e T 1 = (1, 0, . . . , 0) is the first unit vector in R N −1 and E I is the identity matrix of size I. Further, Ω implements the discrete approximation to the k-th derivative given in (4). It is explicitly given by Ω i,i+j = ω (k,t BAT ) ij for i = 1, . . . , I, j = 0, . . . , k − 1, and zero otherwise. For readability, we omit the dependence of the involved matrices on t BAT and k in our notation.
A standard computation (differentiation of the functional in (6) with respect toṽ and setting it equal to zero) reveals thatṽ * is the solution of the linear system Solving this equation forṽ * and using the relation v * = Bṽ * gives the explicit solution where H (t BAT ,α,k) ∈ R N ×N is the influence matrix given by Note that H depends on t BAT , on α, and on k because the right hand side in (8) does so. Thus, for fixed values of t BAT , α, and k, the discrete approximation v * of u * is given by multiplication of the data vector c with the matrix H (t BAT ,α,k) .

Parameter estimation by generalized cross validation
Our principal interest is to determine the BAT by estimating the parameter t BAT . For this, the spline parameters α and k need to be estimated as well because these affect the adaption of the model to the curve shape after the BAT. A widely used approach for determining such parameters is generalized cross validation (GCV) (Craven and Wahba, 1978;Golub et al., 1979), which has shown to provide good performance for spline-based estimators (Aydin and Omay, 2006;Wahba, 1985). The GCV score is a ratio of the mean squared error (MSE), given by v * − c 2 2 /N, and a quantity based on the degrees of freedom (DOF) of the estimator. As the proposed model leads to a linear estimator, the DOF coincide with the trace of the influence matrix H (t BAT ,α,k) . Thus the generalized cross validation score for the proposed model has the explicit form The optimal set of parameters (t * BAT , α * , k * ) in the sense of GCV is the one that provides the minimum GCV score: and k in a set of positive integers {k min , . . . , k max }. A minimizer t * BAT is the proposed estimate of the BAT.
It is possible that (9) has multiple minima meaning that there are multiple equally good candidates for t * BAT in the sense of the proposed model. In such a scenario -which is the exception rather than the rule -the utilized minimization algorithm implicitly takes the decision for one of the minima.
2.1.6 Computational minimization of the GCV score and implementation All computational work was performed in MATLAB R2017b (MathWorks, Inc. Natick, Massachusetts, USA). The finite difference weights in (4) are computed using the implementation "fdweights" of T. Driscoll available at the MATLAB File Exchange. 1 Regarding the minimization of (10) with respect to the spline order k, we first discuss the role of that parameter: By construction, a spline of order k receives no smoothness penalty when approximating data by a polynomial of degree k − 1. This means that long-term polynomial trends of order k − 1 are particularly well captured. (However, note that splines are more flexible than polynomials.) The long-term trends of the considered types of TCCs are beyond linear which suggests to choose k min > 2. On the other hand, approximation with high-order polynomials is prone to overfitting, especially in the case of high noise levels, and may get illconditioned, in particular for longer signals. For the considered types and lengths of TCCs, we found that k max ≤ 6 leads to stable solutions. To improve the stability in case of even higher noise or even longer TCCs, one can consider reducing k max to 5 or 4. Summing this up, we found k min = 3 to k max = 6 to be a reasonable range for the signals considered in this work.
To minimize the GCV score with respect to the two continuously defined variables t BAT and α, MATLAB's solver "fminunc", an iterative solver based on Newton's method, is used with the logarithm of the GCV score as target function. (This is equivalent to minimizing the original GCV score as the logarithm is strictly monotone.) We use the variant based on numeric differentiation since deriving analytic expressions for the Jacobian and the Hessian of the GCV score is involved. To avoid getting trapped in a local minimum, we first perform a coarse search over a grid of α and t BAT and use the optimum values as starting point for the iterative solver. To save computing time, we limit the coarse search interval for t BAT to the interval between the time points t min and t max which correspond to the minimum and the maximum of the signal. Although this restriction did not negatively affect the results for the TCCs considered in this work, one could consider dropping this restriction when processing data with even stronger noise so that the coarse search space does not get restricted in a disadvantageous way. Note that the fine search always uses the full search space [t 1 , t N −k ]. If the injection time of the CA is known, it can be used to obtain an improved lower boundary for t BAT . For the evaluation in this work, we do not assume that the injection time is known, so that the method is applicable even if the injection time has not been recorded. To prevent the iterations to run out of the search domain for (t BAT , α) barrier functions are added. These are functions that are zero in the interior of the search domain and grow fast and smoothly to infinity when approaching the boundaries. Further, computing the influence matrix H based on the normal equation in (8) can be numerically instable when working with higher spline orders and with longer signals. To avoid this we perform the computation based on the least squares formulation in (6) using the QR solver implemented in MATLAB. We also observed that very small values for α often led to overfitting. To prevent this, we limit the search range for α by a minimum value α min which was set to α min = 1 throughout this work. A pseudocode of the proposed method is given in Appendix A. Our reference implementation is provided online. 2

Validation using simulated data
For validation, the proposed method was applied to simulated TCCs and the respective AIFs with known true BAT.

Data simulation
Data were simulated using an in-house developed tool within the "Medical Imaging Interaction Toolkit (MITK)" Nolden et al., 2013). TCCs were generated by forward convolution of an AIF with the respective tissue response function of a pharmacokinetic compartment model. We used the tissue response functions of two established pharmacokinetic models. The two compartment exchange model (2CXM) (Sourbron and Buckley, 2012) is characterized by four parameters, the plasma flow F p , the permeability P S, and the fractional volumes of the extravascular extracellular compartment v e and the plasma compartment v p . In the Extended Tofts Model (ETM) (Sourbron and Buckley, 2011;Tofts et al., 1999) the parameters F p and P S are combined to the transfer constant K trans .
High resolution TCCs were simulated for an acquisition time of 360 s with a temporal resolution of 0.25 s and BAT at 34.75 s (Figure 2). BAT was defined as the last point on the concentration curve's baseline.

Estimation of BAT on simulated data for different temporal resolutions and SNRs
To evaluate the robustness of the proposed method, simulated TCCs were altered for four temporal resolutions (∆t = 1 s, 2 s, 5 s, 7 s) and four noise levels defined by their signalto-noise ratio (SNR) (SNR = 100, 50, 25, 10). SNR was defined as the ratio between the maximum of the concentration curve to the standard deviation of the Gaussian noise added to the concentration curve. Different temporal resolutions were realized by respectively down-sampling the original data. BATs were estimated by the proposed method for 1000 noise realizations for each configuration (curve type, SNR, temporal resolution) to obtain statistical information.
Results were compared to those obtained by the method of Cheong et al. (2003) who are using a linear-quadratic piecewise continuous function for approximating the concentration-time curve  and estimating the BAT. (A more detailed description of Cheong's method and its mathematical relation to the proposed method is given in Appendix B and C). To our knowledge this method is still the current gold standard for BAT estimation.

Estimation of BAT on in vivo rat data
To validate our approach for in vivo data, DCE-MRI data of five tumor-bearing Copenhagen rats (Charles River Laboratories Inc., MA, USA) were investigated. All experiments were approved by the governmental review committee on animal care , and animals were kept under standard laboratory conditions. Briefly, animals were transplanted with fresh fragments of the syngeneic Dunning prostate adenocarcinoma subline R3327-HI (Isaacs et al., 1978) in both distal thighs. When tumors reached a size of 10 × 10 mm 2 they were imaged at a 1.5 T MRI system (Symphony, Siemens, Erlangen, Germany) using an in-house built animal coil. Animals were anaesthetized by inhaling a mixture of 2 % isoflurane (Baxter, Unterschließheim, Germany) and 2 l/min oxygen. T2-weighted images were acquired with a turbo spin echo sequence (TR 3240 ms, TE 82 ms, slice thickness: 1.5 mm, pixel size: 0.35 mm × 0.35 mm) for anatomy references. Two DCE slices (slice thickness: 4.5 mm, pixel size: 0.99 mm × 0.99 mm) were positioned transversely in the rat: one slice in the heart for acquisition of the AIF, the other one through the largest diameter of the tumors. DCE-images were acquired with a T1-weighted TURBO-FLASH sequence (TR 272 ms, TE 1.67 ms, flip angle 20 • ) for a total acquisition time of 380 s with a temporal resolution of 0.75 s. 0.1 mmol/kg Gd-DTPA (Magnevist, Bayer Healthcare Pharmaceuticals Berlin, Germany) was manually injected i.v. approximately 30 s after starting the DCE sequence. Even though the experimenters followed a strict routine for CA injection, some temporal variation for the time point of injection might have occurred between animals, as our experimental set-up only allowed to record the time point of entering the MR room and not the actual time point of injection.
DCE-data was converted to concentration curves in MITK by means of absolute signal enhancement (C(t) = κ · (S(t) − S(0)), where S is the MR signal at time t and κ delineates a proportionality which was set to 1). The AIF was extracted from one voxel in the left ventricle of the rats' hearts with the steepest rise, highest signal amplitude, and least noise. The tumors were segmented based on the acquired T2-weighted images. TCCs were extracted voxel-wise for each tumor. BATs were estimated using the proposed method as well as the method developed by Cheong et al. (2003). Additionally, we compared our results to an adapted version of Cheong's method, which only considers sample points after the CA injection time for BAT estimation (Cheong et al. adapted). Because our experimental set-up did not allow for exact recording of the injection time, we set a plausible common lower boundary. As we can ensure that no CA was injected before 25 s, earlier time points were not considered for BAT estimation. Thus the adapted version of Cheong's method considered the sample point at 25.37 as the first one possible.

BAT estimation on simulated rat data
Exemplary fits by the proposed method to one TCC type for all temporal resolutions and noise levels are displayed in Figure 3. Figure 4 displays the results of BAT estimation for various TCC types and AIFs for various temporal resolutions and SNRs representing different data quality.
Cheong's method provides accurately estimated BATs for AIFs with high temporal resolutions. The accuracy is diminished for AIFs with low temporal resolutions. For curves with a less significant peak in the vicinity of the BAT, such as TCC type 3 for the ETM and TCC type 1 for the 2CXM, the estimated BATs are relatively far off the ground truth and are scattered over a wide range. Cheong's method gives poor results for estimated BATs of the remaining TCCs.
The proposed method gives comparable results for the AIFs at high temporal resolutions and performs better for the AIFs at the lowest temporal resolution. For the TCCs, the proposed method could precisely estimate BATs of all shapes and in most cases did not exhibit much variability as indicated by the small range between the 5 % and 95 % percentiles. A slight overestimation of BATs was observed for TCCs simulated with the 2CXM. This could be explained by the shallow initial rise in concentration observed for TCCs simulated with this model. Furthermore, the results are robust towards increased noise levels and low temporal resolutions of the input data. In almost all configurations the proposed method significantly improves the accuracy of BAT estimation upon the state-of-the-art method for the tested TCCs.

BAT estimation on simulated patient data
Human AIFs are characterized by a higher and steeper initial peak which might lead to different TCC shapes as compared to those of small animals (Figure 2). Results of BAT estimation by both methods are displayed in Figure 5.
Cheong's method gives accurately estimated BATs with only small variation for low temporal resolution TCCs. Some major variation in the obtained results can only be observed for data simulated with the 2CXM at the lowest SNR. Cheong's method performed equally well for all simulated patient TCCs. The proposed method gives results close to the ground truth with a slight overestimation of BATs for low SNR TCCs and data simulated with the 2CXM.
Thus, for patient data, the proposed method gives competitive results, yet Cheong's method performs slightly better.

BAT estimation on in vivo rat DCE-MRI data
The proposed method was also tested on in vivo data that was acquired from 10 tumors in 5 rats. BATs were estimated for each tumor voxel and the results are displayed in Table 1.
The compared methods estimated identical BATs for the rats' AIFs (Table 1, Figure 6 E-F). However, for the TCCs the BATs estimated by Cheong's method are implausible as almost all of them were estimated to occur even before the lower bound of injection times. The adapted version of Cheong's method which excludes implausible sample points below the lower bound of 25.37 s for BAT estimation (see Section 2.3) resulted in this lower bound in most cases. However, visual inspection reveals that the TCCs do not rise significantly before approximately  35 s. This means that the results of the adapted version remain uninformative in most cases. The BATs estimated by the proposed method are very close to each other within one animal and appear to be at reasonable time points by visual inspection (Table 1, Figure 6 A-D). DCEcurves of animal 4 with estimated BATs for both tumors and AIF are displayed as exemplary results in Figure 6. Despite the noticeable variation in TCC shapes occuring even within one tumor, the proposed method fits all curves equally well and estimates reasonable BATs. We emphasize that the proposed method does not require a lower bound of injection times as a priori information.

Computation time
On a laptop computer with 2.2 GHz Intel Core i7, 16 GB RAM, MATLAB R2017, the computation time for the simulated TCCs with N = 361 time samples was as follows: The computation time of the baseline approach (Cheong et al., 2003) is around 0.15 seconds per signal curve. The proposed method took around 108 seconds for processing a single TCC. Processing 100 TCCs simultaneously took around 347 seconds, i.e. around 3.5 seconds per TCC in average. The reason for the more favourable runtimes in simultaneous processing is that a series of expensive matrix operations do not have to be recomputed for every signal. Finally, we found that the optimal order k * of the smoothing splines resulted to be 5 or 6 in most cases.

Discussion
Estimating the time delay between the arrival of the CA in the upstream artery and the tissue of interest accurately is an important step to obtain correct parameters from pharmacokinetic modelling of DCE-MRI data (Mehrtash et al., 2016). One way to reliably obtain the delay is to take the time difference between the estimates of the BATs of the TCC and of the AIF, respectively (Kershaw and Buckley, 2006). The gold standard method for estimating BATs has been developed for TCCs with a fast upslope to the peak near the BAT as often observed in patient data. However, we found that this method does not give accurate BAT estimations in small animal data and shallowly rising concentration curves and to our knowledge, there is currently no commonly accepted method addressing this problem. To close this gap, we have developed a novel method for BAT estimation in such data. The proposed method approximates the TCC by a constant function before the BAT and by a smoothing spline after the BAT. Determination of the BAT and adaption of the spline to the curve shape is done automatically by generalized cross validation. We used simulated data with known ground truth to evaluate the accuracy of our approach and compared our results to the method by Cheong et al. (2003). The results showed that Cheong's method works very well on patient data with fast upslopes but gives unsatisfactory results for rat TCCs. In contrast, the proposed method estimates BATs that are close to the ground truth. The main reason is that the method of Cheong et al. (2003) heavily relies on a fast upslope from the BAT to the peak of the concentration curve which is frequently not the case in small animal data. A central advantage of the proposed method is that it does not need such a special curve characteristic. Instead, the smoothing spline is able to adapt to a variety of curve shapes after the BAT.
The arrival of the CA does typically not coincide with one of the time points of the image acquisition. The state-of-the-art method, however, limits estimation of BATs to those time points which imposes a systematic upper limit to the estimation accuracy. The proposed method estimates BATs on a continuous scale; i.e. the estimated BAT can lie in between the time points of the image acquisition. This improves the accuracy of the results especially for data with low temporal resolution.
Even though the proposed method was developed for small animal data, we found that it is flexible enough to be used for BAT estimation in patient data as well: For the higher temporal resolutions and higher SNRs, the experimental results are competitive with those obtained by the method of Cheong et al. (2003). Yet, for lower temporal resolutions and lower SNRs Cheong's method gave more accurate estimates of the BATs in patient data.
The proposed method was also applied to real in vivo acquired rat data. As the ground truth is not known for this data we can only make qualitative statements on the results. Cheong's approach gave reasonable results for the AIFs which fulfill the criterion of a fast upslope from the BAT to the concentration maximum. For the TCCs, however, the estimated BATs can be considered as implausible because they laid before the actual injection time. Adapting Cheong's method to consider only BATs after a plausible lower time point mostly resulted in that lower bound. In either case, the estimates for the BATs are mostly not informative. In contrast, the BATs estimated for the TCCs by the proposed method were found to be very close to each other, for each tumor and each animal, respectively. When evaluating the location of the BATs within the TCCs, they appear to be set at reasonable positions.
Due to the complexity of the model, the computational effort of the proposed method is higher than that of the method of Cheong, especially for long signals. Generally, the average computation time per signal can be dramatically reduced when processing multiple signals simultaneously (for example, all voxels of all the imaged objects). Computational time can be reduced further at a potential trade-off in accuracy: First the search space for the order k can be restricted to the values 5 and 6 as we observed them to be the optimal order k * of the smoothing splines in most cases. Second, the signal tails could be cropped to reduce the length of the signal.

Conclusion
We proposed a novel method for BAT estimation in DCE-MRI data which is able to adapt to various types of TCCs; in particular it does not require the TCCs to have a fast upslope from the BAT to the peak. The method is based on a flexible spline-based approximation model and automatic parameter estimation by generalized cross validation. We evaluated the method on simulated data (rat and patient) with known ground truth as well as on in vivo rat data. The proposed method was found to be more accurate in rats than the method of Cheong et al. (2003), which has been developed for the analysis of patient data with fast upslopes. Furthermore, the proposed method even gave competitive results in the considered patient data. We propose using this method for reliable BAT estimation of DCE-MRI data that does not have a fast upslope -as often observed in small animal data and patient data of tissues without a significant intravascular fraction -to correct the delay between CA arrival in the AIF and the tissue of interest. the BAT and the peak time: u(t n ) = β 0 , if t 1 ≤ t n ≤ t BAT , β 0 + β 1 · (t n − t BAT ) + β 2 · (t n − t BAT ) 2 , if t BAT < t n ≤ t p . (11) Here, n ∈ {1, . . . , p − 2} and β 0 , β 1 , β 2 are the model parameters. We note that also a linear variant was proposed which is obtained by setting β 2 equal to zero. We found that in average imposing this restriction did not lead to improved estimates, thus we focus on the quadratic method. The method proceeds as follows. For each BAT candidate t BAT ∈ {t 1 , . . . , t p }, the optimal model parameters are computed by a least squares fit using non-negativity constraints.
The candidate that provides the least residual sum of squares is considered as the best estimate for the BAT. Note that with this method the estimated BAT is necessarily one of the sampled time points.
C Mathematical relation of the proposed method to the method of Cheong et al. (2003) As mentioned earlier, in the limit α → ∞, the proposed spline model (2) reduces to least squares approximation by a polynomial of degree k − 1 after the BAT. Hence, the unconstrained parametric model of Cheong et al.,i.e. model (11) without non-negativity constraints on β, appears as the limit case α → ∞ of the proposed non-parametric model of order k = 3. Furthermore, minimizing the GCV for this case reduces to minimizing the residual sum of squares. This is because the denominator in (9) is constant for polynomial approximation of fixed degree. Also recall that Cheong's method limits the search space for the BAT to the sample time points t 1 , . . . , t p−2 . To summarize, the unconstrained variant of Cheong's method can be formulated in terms of the proposed method (9) as the limit α → ∞ of the optimization problem argmin t∈{t 1 ,...,t p−2 } GCV(c 1:p ; t, α, 3) where c 1:p = (c 1 , . . . , c p ) denotes the signal until the time of maximum. Therefore the proposed approach can be seen as generalization of the unconstrained variant of the method of Cheong et al. (2003).