WO2014013476A1 - A control point interpolation method for the quantification of cerebral haemodynamics - Google Patents

Abstract

The present invention provides a method which involves the deconvolution of the concentration time curve in tissue with an arterial input function (AIF), of the following equation: Formula (I) where C(t) is the measured concentration time curve (CTC), (⊗ represents convolution, the proportionality constant alpha is a measure of brain tissue density and difference in hematocrit between capillaries and large vessels, CBF is the cerebral blood flow and R(t) is the residue function, the residue function describing the fraction of a tracer remaining in the tissue vasculature at a time t after its arrival, wherein Control Points (CPs) are used to parameterise the residue function, the Control Points being interpolated to produce a smooth continuous residue function curve.

Description

A CONTROL POINT INTERPOLATION METHOD FOR THE QUANTIFICATION OF

CEREBRAL HAEMODYNAMICS

TECHNICAL FIELD

The present invention generally relates to medical imaging and more particularly, relates to a method for providing clinicians with quantification of cerebral perfusion and residue function shape information to facilitate better diagnostic, prognostic or therapeutic decisions. The present invention is particularly suitable, but not limited, for use where a patient has suffered a stroke or other neuro-pathological condition which affects the cerebral perfusion.

BACKGROUND

Perfusion magnetic resonance imaging (MRI) is widely used in neurological diagnostic imaging for ischemia, stroke and brain tumour assessment. Dynamic susceptibility contrast magnetic resonance imaging (DSC-MRI) is frequently used in the measurement of cerebral perfusion in stroke and other pathological conditions. It has been shown that perfusion parameters like cerebral blood flow (CBF), cerebral blood volume (CBV) and mean transit time (MTT) can be used in acute stroke patients for quantification of cerebral ischemia. DSC-MRI requires the injection of a bolus of paramagnetic contrast agent followed by measurement of the MRI signal loss during its passage through the tissue. The measurement of these perfusion parameters is based on tracer kinetic theory or indicator dilution theory, that considers the tissue concentration time curve as the convolution of the arterial input function (AIF) with a CBF- scaled tissue residue function. The residue function describes the fraction of tracer remaining in the tissue vasculature at a time after its arrival. One major concern in inferring the crucial perfusion parameters from DSC-MRI is reliable and accurate voxel-wise deconvolution of the observed concentration time curve (CTC, C(t)) with respect to the measured AIF.

There are multiple methods to perform deconvolution. Initially an analytical model- dependent deconvolution technique was proposed but it has been superseded by more flexible nonparametric approaches like frequency domain deconvolution (Fourier transform) and later by singular value decomposition (SVD). Variants of the SVD approach such as tSVD (truncated SVD) and the time insensitive oSVD (circular-SVD) are widely used to estimate the tissue response function (TRF; defined as residue function multiplied with CBF), its maximum value being used to estimate CBF. One of the typical assumptions in the DSC-MRI quantification process is the absence of bolus delay and dispersion between the site of AIF measurement and tissue voxel; i.e. the measured AIF is assumed to reflect the exact input function to the tissue voxel being analysed. Partial volume effects in the AIF measurement arise due to the spatial resolution of a typical DSC-MRI acquisition, whereby both arterial and non-arterial components occupying AIF chosen voxels. To minimise partial volume effects, the AIF is typically measured at a major cerebral artery, such as the middle cerebral artery (MCA). This single AIF is then used as the input function for the whole brain and this strategy can contribute to bolus delay and dispersion. Such errors might cause inaccuracies in perfusion quantification, particularly in patients with cerebral ischemia, but they may even be present to a lesser extent in healthy volunteers.

Bolus delay has long been recognised as a source of systematic error in CBF quantification and deconvolution strategies have been designed to be delay insensitive.

The key issues with SVD-based methods are the underestimation of CBF and also the introduction of non-physiological oscillations in the resulting residue function, making the resulting residue function difficult to interpret. There have been attempts with various regularization approaches to reduce the oscillation in the residue function. Some of the widely used regularization methods involve: truncated SVD (threshold using Psvd), oSVD using oscillation index (OI) or adding of a penalty term in the optimization, for example, the Tikhonov regularization. Despite all these efforts the estimation of realistically smooth residue functions has not yet been achieved with SVD methods.

More recently, a vascular model (VM) has been proposed, this is implemented within a Bayesian framework as an alternative to SVD methods. The model is based on a vascular architecture where a gamma probability distribution function is assumed to model the underlying tracer transit times. The VM perfusion values are believed to produce more accurate estimates compared to the SVD method, but, being a model-based solution, it lacks the flexibility of SVD methods. In particular, the underlying gamma distribution in the VM is controlled by two coupled parameters that permit only a restricted family of shapes for the tissue response function which might not be appropriate in modelling altered vasculature in pathologies like stroke. The shape of residue function is also of interest because it contains information about tissue vascular integrity. For in vivo analysis the actual residue function shape is not known a priori and, particularly in pathology, might not be drawn from the set of functions currently assumed for typical residue functions.

Dispersion causes a smearing of the AIF function in time, altering the bolus profile and affecting the features of AIF prior to its arrival at the tissue voxel. If the temporal spread is not considered in the perfusion analysis modelling, these distortions can lead to considerable errors in the quantification of the perfusion parameters, which can have serious implications for diagnosis and management of patients with cerebral ischemia.

Accordingly there is a need to address the aforementioned deficiencies. The aim of the present invention is therefore to provide an algorithm that overcomes the deficiencies named above. The present invention is a deconvolution algorithm that produces accurate perfusion values over a wide range of physiological conditions.

A deconvolution method is described herein. This method can estimate cerebral perfusion along with a physiologically plausible residue function without requiring it to belong to a specific class of functional shapes. The resulting residue function shapes can be used to assess the residue function variability among different brain regions and its changes in pathology.

SUMMARY

The present invention provides a control point interpolation method for quantification of cerebral haemodynamics. In a preferred embodiment the method involves the deconvolution of the concentration time curve in tissue with an arterial input function (AIF), of the following equation:

C (t) = a^■ CBF · (Ca (t) <g) R (t))

where C(t) is the measured concentration time curve (CTC), <S> represents convolution, the proportionality constant a is a measure of brain tissue density and difference in hematocrit between capillaries and large vessels, CBF is the cerebral blood flow and R(t) is the residue function, the residue function describing the fraction of a tracer remaining in the tissue vasculature at a time t after its arrival, wherein Control Points (CPs) are used to parameterise the residue function, Control Points being interpolated to produce a smooth continuous residue function curve. In another embodiment, the quantification of cerebral haemodynamics is from dynamic susceptibility contrast magnetic resonance imaging (MRI).

In a further embodiment the optimisation method comprises a Bayesian Inference Scheme or a non-linear least-square scheme.

In yet a further embodiment the CPs has degrees of freedom in both amplitude and time.

In another embodiment each consecutive control point is related to its precursor by a ratio factor (Ω) and a time spacing (ζ), according to the following equation = __! * n_t where is the amplitude of the i"¹ CP and : = t-i + ζί In yet another embodiment, the first CP is fixed such that ti=0 with ri=l.

In another embodiment, the first CP is fixed at zero amplitude or allowed to vary in amplitude.

In yet a further embodiment, t is fixed at the last time point of the measured CTC.

In another embodiment the CTC is converted to a Signal Estimate S(t) using the following equation:

TE

S (t) = S (t₀) - exp— ^c «

where S(t) is the, S(t₀) is the baseline signal, TE is the echo time and κ the proportionality constant.

In a yet further embodiment, prior knowledge is incorporated for each parameter. In a further embodiment, the prior knowledge is in the form of a Gaussian distribution for each parameter with mean and standard deviation values.

BRIEF DESCRIPTION OF THE DRAWINGS

Many aspects of the disclosure can be better understood with reference to the following drawings. The components in the drawings are not necessarily to scale, emphasis instead being placed upon clearly illustrating the principles of the present disclosure. Moreover, in the drawings, like reference numerals designate corresponding parts throughout the several views. Figure 1 depicts the non-linear single input single-output system for the Control Point

Interpolation (CPI) method of the present invention.

Figure 2 depicts the CPI framework and its optimisation parameters like compensating for delay (δ) in AIF and estimating the shape of residue function with η control points (CP), such that CP i (t=0, r=l) and CP,_! (t=last time point of CTC, r=variable). Figure 3 depicts simulation results showing the estimated CBF plotted against the true

CBF values (SNR=20 and CBV 4%; the different scale of estimated CBF is shown in a, b and c).

Figure 4 depicts the mean percentage error in estimation of CBF (a) and mean absolute error in bolus arrival delay (b) for the present method with respect to simulated delay values of 0 sec and +5 sec with CBF=10-70ml/100gm/min, CBV 4% and SNR=20. Figures 5 (a, b, c) illustrate the results from one of the simulated datasets comparing the

TRF shape estimation for a CBF of lOml/lOOg/min (cerebral ischemia).

Figures 5(d, e, f) illustrate the mean shape estimated for the corresponding simulated residue functions with + σ bands for CP-interpolation.

Figure 6(a) illustrates the mean estimates of CBF calculated by the present method for a simulated residue function with a large MTT that could represent T2* leakage effect.

Figure 6(b) illustrates the simulated and estimated TRF for CBF=10ml/100g/min by the present method for simulated data with large MTT that may be consistent with T2* leakage effect. Figures 7(a, b, c) depict axial brain slices showing perfusion maps estimated by CP- interpolation method for a healthy participant.

Figure 7(d) shows a T2* gray scale image with two selected ROIs for residue function shape analysis.

Figure 7(e) shows the estimated residue function for two representative ROIs.

Figures 8(a, b, c) depict an axial slice showing perfusion estimates from CP-interpolation method for a patient with atherosclerotic diseases.

Figure 8(d) depicts the diffusion weighted image of the patient, showing an area of infarcted tissue in the left parietal lobe.

Figure 8(e) illustrates the obtained residue function for two ROIs.

Figure 9 depicts absolute simulated CBF against estimated CBF (SNR=100 and CBV 4%) for residue function observed using the CPI method in the analysis of a healthy participant.

Figure 10 depicts simulated TRF, dispersed TRF, mean fit achieved with oSVD, VM, CPI (a,c); VM+VTF, CPI (CPi=0), CPI (CPi=var.) (b,d). Gamma Dispersion Kernel (a,b) and Exponential Dispersion Kernel (c,d).

Figure 11 shows CBF, MTT and Delay estimates for a representative patient with CPI method. The DWI, R50 and RIO map of the same patient (bottom row), where R50 and RIO map clearly showed region of high R50 and RIO values in region of infarct surrounded by an area under ischemic stress in the left hemisphere compared to normal right hemisphere. Figure 12 (left) illustrates the Control Point Interpolation methods as elaborated in Figure

2; (right) top: dispersion variants of CPI method, CPi fixed to zero (CPi=0), bottom: CPI+VTF variant of CPI accounting for dispersion with Gamma Dispersion Kernel (GDK) as a function time-to-peak (p), sharpness (s) and time (t) as vascular transport function (VTF).

Figure 13 shows MTT maps (a) CPI; (b) CPI+VTF, showing ROI selected in healthy and in DWI positive ischemic tissue, c) CPI and CPI+VTF shows similar residue function in contralateral ROI with no dispersion, d) Infarcted ROI showed broader and slower decay in residue function with CPI compared to CPI+VTF. CPI+VTF also showed slower decay than ROI in healthy ROI as in c. Larger amount of dispersion was observed in ischemic ROI than normal ROI as seen with much broader residue function with CPIQ. Figure 14 shows DWI (a,e), RIO (b,c,f,g), S (dispersion) (d,h) maps for CPI (b,f) and CPI+VTF (c,d,g,h) showing area of infarction (from DWI) with higher values of RIO with both CPI and CPI+VTF. b) RIO with CPI also showed higher values in left MCA territory which improved post-CEA (f). c) RIO with CPI+VTF showed much reduced values in left MCA territory where high amount of dispersion was found (marked with arrow in d.). Dispersion was found to completely disappear post-CEA (h).

Figure 15 depicts an axial brain slice from a patient with atherosclerotic diseases showing MTT, CTTH and Skw both pre- and post-CEA with CPI and CPI+VTF. MTT and CTTH tend to reduce to normal values (arrow) and Skw tend to increase (arrow) to normal vales after CEA. Figure 16 depicts the mean residue function and distribution of transit times observed with CPI for normal, DWI-ring and DWI+ tissue pre and post-surgery.

DETAILED DESCRIPTION Having summarized various aspects of the present disclosure, reference will now be made in detail to the description of the disclosure as illustrated in the drawings. While the disclosure will be described in connection with these drawings, there is no intent to limit it to the embodiment or embodiments disclosed herein. On the contrary, the intent is to cover all alternatives, modifications and equivalents included within the spirit and scope of the disclosure as defined by the appended claims.

The theoretical concepts of tracer kinetic theory are described firstly; followed by the framework for the CP-interpolation method and its procedure. Subsequently, the simulation protocol that was used to evaluate the method will be elaborated. Finally, examples are given where the CP-interpolation method was applied to analyze data from a healthy participant and a cerebrovascular disease patient.

Theory

For a given tissue voxel it is assumed that the intravascular tracer delivery to the capillaries can be represented by an Arterial Input Function (AIF), denoted by Ca(t). According to indicator-dilution theory the concentration of the tracer in the capillaries at time t can be expressed as the convolution of Ca(t) with the tissue residue function R(t), scaled by the cerebral blood flow (CBF):

C (t) = a^■ CBF · (Ca (t) <g) R (t)) Equation 1

Where C(t) is the measured concentration time curve, (¾ represents convolution as defined on the right-hand side of the equation and the proportionality constant a is a measure of brain tissue density and difference in hematocrit between capillaries and large vessels (compensating for the fact that only plasma volume is accessible to contrast agent). Since the parameter a is generally indeterminable, it is usually replaced by a fixed value. The residue function is a monotonically decaying function with initial value of one, R (0) = 1. The analysis of Equation 1 typically involves deconvolution of the chosen AIF from the measured CTC. The resulting tissue response function contains the information about both CBF and the residue function shape. The estimated residue function shape itself contains information on microvascular flow heterogeneity. Another important perfusion parameter calculated from the CTC is the mean transit time (MTT). This signifies the average time for a contrast molecule to pass through the tissue vasculature following an ideal bolus injection and can either be calculated using the central volume theorem

CBV

MTT = Equation 2

CBF

where,

/ C(t) dt

CBV = Equation 3

/ C (t) dt or, estimated as the area under the deconvolved residue function :

Equation 4

The CTC to signal conversion relationship is included so that it is the measured DSC that is being modelled: S (t) = S (t₀) · exp- k Equation 5

where S (t) is the signal, S (to) the baseline signal, TE the echo time and k the proportionality constant.

CP-Interpolation Method

The present invention is a model-based approach to deconvolution, which can be used, for example, within an algorithm for the estimation of CBF from DSC-MRI perfusion data using the theory above. In the present invention, the model is structured so that it can provide the flexibility afforded by non-parametric approaches subject to suitable constraints on the smoothness of the residue function. The present CP-interpolation method models the residue function using a restricted number (η) of control points that have freedom in both amplitude and time. These control points are further parameters of the model to be determined from the data (Fig. 2). The control points act as virtual nodes, the full shape of the residue function being generated using an interpolation algorithm. In the present example, piecewise cubic spline interpolation is used, but alternatives such as truncated sine, windowed sine, nearest neighbour, natural neighbour, linear, quadratic, cubic, akima polynomial, B-spline, Lagrange and Gaussian interpolation algorithms could be used. Constraints on the control points may be implemented in the model such that the resulting residue function is plausible. Each consecutive control point may be related to the precursor by the ratio factor (Ω) and a time spacing (ζ):

n = - i * Equation 6 where is the amplitude of the i CP and :

U = ti-l + ζί Equation 7 where tj is the time point of the i CP The following is an example of the method and is not intended to suggest that these are the only parameters or values of the parameters that may be used. Other suitable parameters and their values may also be used.

The first control point may be fixed so that ti = 0 with ri = 1 to meet the fundamental definition of the residue function that R (0) = 1. The first control point may also be either fixed at zero amplitude or allowed to vary in amplitude. Additionally, ί_η may be fixed at the last time point of the measured CTC, but with free amplitude calculated from Equation 6, to avoid extrapolation errors. To encourage plausible residue function shapes and numerical convergence of the solution other constraints may be incorporated in the CPI method such as: CBF cannot be negative (CBF > 0); the maximum value permitted for Ω is 1.5 (Ω < 1.5) and a maximum of three control points can be put within a 77? (ζ > TR/3). The maximum permitted value of Ω may be selected to be 1.5 to accommodate the associated uncertainty in the real clinical data like bolus dispersion. Thus the total number of estimating parameters in the CP-interpolation method is 2η-1 (from parameters CBF (1), δ (1), Ω (η-l) and ζ (η-2)). The CTC estimate may be found using Equation 1 and converted to the signal estimate using Equation 5. In this specific DSC-MRI application CTC to signal conversion was required, however, generally the results of convolution are not required to be converted to other forms (signal), for example in arterial spin labelling perfusion MRI, the CTC could be used directly.

Figure 1 shows the non-linear single input-output system for the CP-interpolation method of the present invention. The tissue residue function may be estimated by setting control points (η) that form the basis of a smooth interpolation along with parameters CBF and δ.

Figure 2 shows CPI framework, compensating for delay (δ) in AIF and estimating the shape of residue function with η control points (CP), such that CPi (t = 0, r = 1) and CP_η (t = last time point of CTC, r = variable).

The convolution may be performed between the AIF and the parameterized residue function scaled with estimated CBF to generate a CTC. The parameters may be estimated by optimizing the resulting signal (from estimated CTC) with respect to the measured MR signal. The estimating parameters are Ω, ζ, δ and CBF. Methods

Implementation

For parameter estimation a non-linear optimization algorithm may be used such as nonlinear least squares. A Bayesian non-linear model fitting algorithm may be employed that permits prior information to be supplied in the form of a probability distribution.

In this example a total of 12 (η=12) control points are used, however, it is understood that more or less control points may be used depending on the SNR, prior knowledge of residue function shape or for optimising the algorithm processing time. The control points may be initialized according to an exponential residue function shape with a logarithmic temporal spacing to encourage more control points at earlier time points where the residue function may have larger amplitudes. Of course, if computational time is not important then the use of more control points is advantageous. The complete estimated residue function curve may be generated from the control points by interpolation using cubic Hermite piecewise splines. The final residue function curve generated with splines may be used in convolution with the measured AIF to form a predicted CTC after scaling with the CBF estimate. The AIF may be shifted according to the bolus delay estimate (δ), such that the AIF shape remains the same but shifted temporally by δ units. The CTC may be converted into a predicted signal using Equation 5. Ω, ζ may be estimated on a logarithmic scale, as this makes the model more approximately linear with respect to these parameters, aiding the convergence of the fitting algorithm, which depends upon a local linear approximation to the function using a Taylor expansion. Hence the estimated parameters (Fig. 2) in the CP-interpolation method are Θ (log (Ω), log (ζ), δ, CBF).

Prior Distribution

When using a Bayesian optimization method, prior knowledge may be incorporated in the form of a Gaussian distribution for each parameter with mean and standard deviation values given in Table 1. A non-informative prior may be implemented for CBF using a large standard deviation and zero mean, whereas δ may have a reasonably informed prior. For Ω and ζ the priors may be chosen so as to promote a smooth residue function whilst retaining the flexibility of the residue function shape to adapt to the data. Table 1 presents the priors that have been used here, priors may be varied according to requirement in different applications. Table 1: Priors for the CP-interpolation method.

7

Parameter Mean (μ) Variance (σ )

CBF (a.u) 0 10^b

Delay (S) (sec) 0 10

Log (Ω) Log (0.5)= -0.69 0.67

Log (9 Log (8) = 2.08 0.40

The method may be initialized with CBF and δ values from another method, such as the widely used singular value decomposition deconvolution, as these may provide a reasonable first estimate in the vicinity of the true solution.

Simulations

An AIF was simulated using a gamma-variate function which has been described previously:

Equation 8

where, to is the bolus arrival time, a, b and c are constants for which typical values were used: a=l, b=3, c=1.5. Simulations were be performed with CBV= 4% (4 ml/lOOg) and 2% (2 ml/lOOg) which are representative of normal gray and white matter. CBF values were varied from 10 to 70 ml/lOOg/min in increments of 10 ml/lOOg/min for CBV=4 and from 5 to 35 ml/lOOg/min in increments of 5 ml/lOOg/min for CBV=2 . Thus MTT varied in the range 3.43s to 24s. The actual tissue response function is usually not known a priori and in the prior art fairly simple residue functional shapes have been assumed. In the present invention three different types of residue functions were simulated:

1 ) Exponential Equation 9

2) Linear R (t) = Equation 10

3) Box-car R

⁽'Ho¹ Equation 11

The exponential residue function describes a well-mixed compartment model whereas the box-car function describes a vascular bed with "plug" flow where the capillaries are in parallel with equal length and mean transit times. In addition, a further simulation was carried out with a residue function that sought to model a CTC with very large MTT that could represent a T2* tracer leakage effect, where the CTC did not return to baseline and remained at higher values throughout the measurement period. It was generated as a special case of the exponential residue function (Equation 9) but with a considerably larger value of MTT of 2 minutes. Since, for this case, an MTT has been chosen and not calculated from CBF and CBV the central volume principle (Equation 2) will no longer apply.

A linear relationship may be assumed between tissue relaxivity (R₂*) and C(t). Concentration time curves were generated as described in Equation 1 and then converted to signal curves S(t) as in Equation 5 with S(t0)=l00 and 7E=65ms. The constant k was chosen such that 40% peak signal drop was achieved at a flow of 60 ml/lOOg/min with CBV=4%. The AIF signal curve was generated similarly but a value of proportionality constant k was chosen such that 60% peak signal drop was achieved. Zero mean Gaussian noise was added to the signal curves to produce a baseline SNR of 20 and 100 simulating typical DSC- MRI and the best possible image data acquired respectively. Values of +5sec, Osec and -5sec were used to evaluate the model for sensitivity with respect to bolus arrival delay (δ). For each combination of CBV, CBF, residue function and delay a total of 100 CTC were generated. Thus a total of 16,800 (2 x 7 x 4 x 3 x 100) concentration curves were analyzed during each simulation process with the method of the present invention.

In vivo Clinical Data

For clinical evaluation of the method, in vivo data was acquired from one healthy participant and one patient with a history of atherosclerotic disease. Images were acquired under Institutional Review Board (IRB) approved protocol for DSC-MRI study. The MRI acquisition was performed on a Siemens 3T Trio, GRE-DSC:77?/ E=1.5sec/30msec, 78 volumes, 128x128x22 matrix, 1.7x1.7x5mm voxels. Of course, any other suitable instrumentation may be used. An intra-venous (IV) bolus injection of 0.1 mmol/kg Magnevist® was performed followed by a 20 ml saline flush. DSC images were then analyzed with the method of the present invention.

The AIF was selected in the in vivo clinical data using an automated algorithm written in- house on a 30mm slab centred at the Circle of Willis. CTC features including the first moment (FM), peak height (Pmax), time to peak (TTP) and bolus arrival time (Td) were extracted and using a K-means clustering algorithm voxels were grouped into five classes (grey matter, white matter, vessel, CSF and background). The AIF was defined in the vessel subclass with highest Pmax and lowest FM and a venous output function (VOF) in the same class but with highest Pmax and greatest FM. The AIF was then scaled such that the mean value after the first passage of bolus (the tail of AIF) was equal to the tail of the VOF using the tail scaling technique.

Statistical Analysis

A regression coefficient (Y) was calculated for estimated vs. absolute simulated flow values and used to analyse the method of the present invention. Coefficients of determination (R ) were calculated using sum of squared error method and employed to evaluate the signal fit achieved from the method as well the fit of the estimated residue function against simulated. R may vary in range from 0 to 1, where zero represents a poor estimated values and values close to one represents a good agreement between simulated and estimated shape.

Results Simulations

Figure 3 shows the mean estimated CBF (10-70ml/100g/min) for the present method

(SNR=20 and CBV=4 ), regression coefficients being given in Table 2. Lower CBF values were estimated accurately by the method of the present invention.

The mean percentage error in estimation of CBF and absolute error in bolus arrival delay is shown in Figure 4 for simulated δ= Osec and +5sec. The method of the present invention is insensitive to the bolus arrival delay for estimation of CBF (fig. 4a), however, the present method of CPI shows a very low error in absolute delay. Thus the CP-interpolation method is able to estimate the actual delay values with great accuracy (fig. 4b).

Table 2: Regression coefficient (Y) for estimated vs. absolute simulated flow values. SNR=

20

CBV = 4 % CBV = 2 %

Linear Box Linear

Method Exp. Func. Exp. Func. Box Func.

Func. Func. Func.

CPI 0.96 1.17 1.55 0.94 1.15 1.43 Figures 5 (a, b, c) illustrate the results from one of the simulated datasets comparing the TRF shape estimation for a CBF of lOml/lOOg/min. The CP-interpolation method is generally found to be able to estimate a residue function shape in keeping with the original without any oscillations. Figures 5(d, e, f) show the mean and standard deviation (σ) over the 100 datasets of the estimated residue function from the CP-interpolation method of the present invention for data generated with exponential, triangle and box residue functions with CBF=10ml/100g/min and CBV= 4%. The mean estimate of the shape of residue function was found to closely resemble that simulated with a narrow variability across the 100 separate instances.

The present method exhibited a good correlation as measured by the coefficient of determination (R ) for the fitted signal against the true simulated signal, the CP-interpolation method having a correlation of (0.98+0.02).

Table 3 illustrates the corresponding coefficients of determination between the estimated residue function shapes and the simulated residue functions.

2

Table 3: Coefficient of determination (R ) for estimated residue function against the simulated residue function

Method R² (Exp. R² (Linear R² (Box R² (Bolus Leak Func.)

Func.) Func.) Func.)

CPI 0.96 + 0.03 0.95 + 0.04 0.82 + 0.10 0.75 + 0.11

Figure 6(a) illustrates the mean estimates of CBF calculated by the method of the present invention for a simulated residue function with large MTT or a T2* leakage effect. The corresponding regression coefficients (Y) are: Y (CPI) = 1.10.

Figure 6(b) shows the representative mean residue function solution obtained for this special case of residue function simulated with CBF=10 ml/lOOg/min (typical for cerebral ischemia). The CP-interpolation solution of the present invention matched the simulated residue function closely. Similar shapes of residue function were obtained for other CBFs simulated with this residue function. In vivo Data

Figures 7(a, b, c) show representative perfusion maps for a healthy participant estimated by the method of the present invention, along with residue function from two ROIs (region of interest of 3x3 voxels) (figure 7e). The CBF map, (figure 7a) shows a good gray/white matter contrast and a larger dynamic range for CBF estimates with the present CPI method. A large variation in MTT was observed with the CPI method. The method of the present invention showed a higher MTT in the ventricles as expected in the case of normal intact blood brain barrier. The residue function shapes estimated by CPI method of the present invention were found to be smooth as anticipated.

The data from a patient with atherosclerotic disease are shown in Figure 8, where the CPI methods shows high values of MTT in a small region of the left parietal lobe; the diffusion weighted image confirmed the presence of infarcted tissue within this lobe. Two ROIs were selected (3x3 voxels) (figure 8e), one in the region of infarcted tissue and another in a nearby normal appearing brain parenchyma. The CPI method demonstrated substantial change in the residue function shape between the two ROIs. With the CPI method the residue function from infarcted tissue ROI (Figure 8e) appeared to have had a slower decay consistent with high MTT observed in the region.

A CP-interpolation technique has been provided for the estimation of perfusion-related metrics in DSC MRI. The CP-interpolation method may use Bayesian optimisation. Singular value decomposition and its variants have been used extensively in DSC MRI analysis despite estimating physiologically unrealistic oscillatory residue functions.

The aim of the present CPI methodology is to improve the estimation of CBF by accurate estimation of the tissue residue function. In the present approach, the tissue residue function was estimated at a subset of points and then cubic spline interpolation was used to generate the complete smooth function. As mentioned above, other interpolation methods may also be used. The simulation study showed that the CP-interpolation method was accurate as shown by the larger values of the regression coefficients for simulated exponential residue functions.

The actual tissue response function is not known a priori and in the prior art fairly simple residue functional shapes (exponential or box) have normally been assumed. In an example of the present invention a more complex residue function was simulated with a large MTT that is consistent with T2* tracer leakage effect where CTC did not return to baseline but remained at a higher value throughout the measured time frame. A region like this should typically have higher values of transit time for a given CBF value, since the tracer will remain within the voxel for an extended time. As illustrated, the CP-interpolation method could estimate both the correct CBF and residue function shape. The CP-interpolation method of the present invention is able to accommodate a wider range of residue functions as demanded by the data.

The residue function shapes as estimated from CP-interpolation method in the clinical perfusion data (figures 7 and 8) were smooth and revealed that the residue function had both fast and slow components. The residue function as obtained from CPI method in normal healthy participant was averaged over a region of interest (4x4 voxels) and further simulated data were generated using this function and the similar methodology used above. The results showed that even with this complex shape of residue function the CPI method of the present invention estimated the CBF values correctly (Y (CPI) = 0.96), the absolute simulated CBF values (figure 9).

The results from the CPI method of the present invention suggest that in the in vivo data, the residue functions deviate from those assumed by the gamma distribution of transit times that form the basis for prior art methods. This may partially reflect the presence of contamination from contrast agent in larger arteries. Dispersion of the contrast bolus between the site of measurement of the AIF and the voxel is also a widely documented effect that will be represented by a change in the residue function shape if not accounted for elsewhere in the analysis. One of the strengths of a non-parametric deconvolution method of the present invention is that it is adaptable to such sources of variation that otherwise would violate the assumptions of a model-based approach such as the prior art methods. The CPI method can be viewed as effectively non-parametric, in that it retains the flexibility in the residue function shape. However, the estimation of the residue function is subject to constraints that regularize the shape avoiding the oscillations of prior art solutions. The present invention method provides a means by which well-behaved residue functions can be estimated.

The effectiveness of the CPI solution of the present invention was illustrated in a patient, where changes in the residue function were observable between healthy and infarcted tissue. The results show that the CPI method of the present invention may be able to offer residue function estimates that could be interpreted in terms of the underlying physiology. This would include the use of the residue function to probe the heterogeneity of the flow within the tissue. Furthermore, other comparisons and results are given below. BOLUS DISPERSION Bolus dispersion is a recognised effect in perfusion techniques that involve a global arterial input function (AIF) from a site distant that is used to characterize tissue signals, typically via a deconvolution operation; the problem is likely to be amplified in cases of cerebral ischemia. Dispersion introduces errors in the Cerebral Blood Flow (CBF) estimates and it is generally not possible to separate the effects of dispersion from the true residue function when using the common singular value decomposition methods.

Dispersion on two levels was investigated: 1) in terms of CBF accuracy, and 2) the ability to isolate dispersion effects in the residue function. Herein, a range of deconvolution approaches is investigated, including those specifically designed to handle bolus dispersion. Specifically, the control point interpolation (CPI) method is used to address the issues of dispersion.

Theory of Bolus Dispersion in Perfusion Imaging

The residue function scaled with CBF is referred to as the tissue response function (TRF (t) = CBF R(t)). Bolus dispersion affects the AIF, by altering its shape and amplitude while in transit between the measured AIF location and the tissue. This can be reconciled by means of dispersion parameters in the DSC model: for example, through the inclusion of a Vascular Transport Function (VTF) in convolution with the AIF. The VTF thus encapsulates the effects of dispersion in a simple mathematical function, giving a modified C(t):

C (t) = a · CBF · (Ca (t) <g) VTF(t) <g) R (t)) Equation 12 where Ca(t) is now the measured AIF.

The residue function is defined as a monotonically decaying function with an initial value of one, R(0)=1. Once dispersion is included in the perfusion model, the sharp AIF profile is now 'smeared out' over time by the VTF. Since convolution is associative, the VTF can instead be convolved with the residue function (keeping AIF unaltered) to examine the effects of dispersion on the residue function as would be estimated by perfect deconvolution. In such a scenario, the measured residue function may no longer start from one (i.e. R(0)≠l).

Vascular Transport Function Model

The VTF depends on several factors, such as vascular topology, tissue type, site of AIF measurement, CBV and the pathological condition of the vessels (vessel lumen, elasticity, stenosis). Depending on these pathophysiological conditions, an assumed VTF kernel functional form could be employed; for example, a gamma dispersion kernel has been considered as a reasonable approximation to model the effects of dispersion:

sl+sp

VTF t) = ^■ t^{sp ■} e^~st Equation 13

Γ (1 + sp)

where, s characterizes the "sharpness" of the kernel and p is the time-to-peak {ttp). When p=0 and as s -∞ this kernel approximates a delta function resulting in zero dispersion; whereas with larger value of p and smaller value of s the kernel will be associated with a greater amount of dispersion.

Methods

Simulations were performed with Cerebral Blood Volume=4ml/l OOg, CBF in the range 10-60 ml/lOOg/min, exponential tissue residue function (TRF), and two vascular transport functions (VTF) to model dispersion: 1) Gamma and 2) Exponential kernels. The Gamma kernel was parameterised by the following: time-to-peak (ttp) and sharpness (s); a range of low to high dispersion were evaluated using ttp = l-5sec; s = 2-0.5sec^_1. The Exponential kernel was parameterised by one parameter time-constant (tconst) that was varied from l-5sec. Concentration time curves were generated, converted to signal time course, to which Gaussian noise was added to achieve an SNR = 50 (100 simulations).

Quantification of residue function fitting achieved by the methods was performed by calculating the coefficient-of-determination, R for estimated residue function against simulated, and accuracy of CBF estimation with the regression coefficient {Y) for estimated CBF against simulated CBFs. In CPI, the residue function is estimated from a set of control points (CP) that form the basis of a smooth piecewise cubic spline interpolation. Each CP is allowed to vary in both amplitude and time, except for the first CP that is set to 1 (CPi=l). This CPi restriction could lead to errors when applied to dispersed data, where the initial value of the residue function subject to dispersion may no longer be unity.

Two possible solutions were considered: i) CPI (CPi=0): CPi fixed at 0, ii) CPI (CPi=var.): CPi has a variable amplitude (range 0-1) that is an additional parameter of the method.

Results

Figure 10 shows the simulated and dispersed tissue response function and mean fitting achieved with oSVD, VM, CPI (a,c) and with dispersion corrected methods (b,d), for gamma and exponential dispersion kernels used.

Table 4 show the results with a moderate degree of dispersion (ttp=1.5, s=0.7, tconst=2). oSVD showed poorest performance at CBF estimation; VM+VTF with dispersion correction was found to be superior to VM alone in residue function shape estimation but CBF estimates were still not accurate. Both dispersion corrected variants of CPI method performed better than other methods used. CPI (CPi=0) was found to be more precise at residue function shape estimation (R²=0.93-0.97) (fig. 13 b,d) whereas CPI (CPi =var.) was found to be more accurate in CBF estimates (F=0.93-0.99) in comparison to all other deconvolution methods used. Similar findings were obtained for the other simulated dispersion conditions.

Exponential Y 0.65 0.67 1.11 0.70 1.23 0.99 Dispersion R^z 0.83 0.75 0.90 0.73 0.93 0.77

Bolus dispersion is a practical problem in perfusion analysis, which is often ignored because of its unknown characteristics. It was hypothesised in the past that the true residue function and dispersion could be separated if a model-based approach is used (parameterising both the true residue function and dispersion). However, these simulations suggest that even a complete model-based deconvolution (VM+VTF) cannot clearly distinguish the two, leading to underestimation of CBF. Models with a more constrained relation between bolus arrival delay and dispersion have been proposed, but the existence of such a relationship still needs to be validated. Additionally, in vivo the actual residue function shape and VTF are both not known a priori, particularly in pathology. The CPI approach provides a method to avoid the strict model- based assumptions for the residue function and VTF, whilst still offering a smooth interpretable residue function estimate. Here a variant on the original CPI method was found to offer reasonable CBF estimation in the face of dispersion. However, CPI does not directly permit the separation of true residue function from the effects of bolus dispersion. Ultimately a further deconvolution of estimated dispersed residue function with an appropriate model of VTF might be required to separate these effects, where again some assumption has to be made for VTF.

DSC-MRI DSC-MRI analysis is an ill-posed inverse problem that involves the deconvolution of the

MR signal with an arterial input function. One aspect of quantification is the use of deconvolution techniques to estimate the residue function. The residue function is an important summary of capillary haemodynamics that might provide useful clinical information regarding evolution of stroke, for example, but has been hampered by non-physiological oscillatory signals when using Singular Value Decomposition (SVD) deconvolution or its variant methods.

Control Point Interpolation (CPI) deconvolution has recently been proposed with demonstrated ability to characterise the residue function from clinical data, as detailed above. The changes in residue function behaviour was investigated, along with flow heterogeneity in 17 patients with atherosclerotic disease by comparing results from SVD and CPI deconvolution methods.

A further aspect is to quantify the non-parametric residue function into physiologically interpretable and clinically useful parameters.

Material and Methods:

DSC data were acquired from five patients (median age: 65yrs [(47-85 yrs], M:F=3:2) with atherosclerotic diseases under an Institutional Review Board approved protocol. MRI data were acquired on a Siemens 3T Trio scanner with Diffusion Weighted Imaging (DWI) and GRE- DSC: TR/TE=1.5s/30ms. 78 volumes, 128x128x22 matrix, 1.7xl.7x5mm³ voxels. An intra- venous bolus injection of 0.1 mmol/kg Magnevist was performed followed by a 20 ml saline flush.

DSC images were analysed using the oSVD and with CPI deconvolution. In the CPI method, the tissue response function was estimated at a set of control points and then cubic spline interpolation was used to generate the complete smooth residue function.

A Region of Interest (ROI) analysis was performed both pre and post carotid endarterectomy on the DSC data. Residue function characteristics were evaluated from ROIs based on DWI and the perfusion weighted images under three criteria:

I) normal perfused regions 2x1 voxels (normal) in the contralateral hemisphere that was outside the hypoperfused and DWI infarct region;

II) region around DWI lesion (DWI-ring) but normal DWI; and

III) showing positive infarct tissue within a DWI lesion (DWI+).

A Maximum of four ROIs were selected from each patient in respective groups of normal, DWI-ring and DWI+ tissue. In total 20 ROIs were selected from patients under normal group, 20 in DWI-ring and 20 ROIs under the DWI+ group. Rigid body image registration was performed between pre- and post-surgery perfusion and DWI images using the FMRIB Linear Image Registration Tool (FLIRT) from the FSL toolbox before ROI selection. ROIs were selected in the post-surgical imaging data by creating a ROI mask in pre surgery images and transforming it to later using FLIRT algorithm. Residue functions obtained from CPI deconvolution among three ROI groups were analysed for variation in shape and heterogeneity. The residue function is a monotonic decreasing function; in order to characterise it, we measure the time taken by the residue to drop to 50% and 10% of its maximum value for three ROI groups, referred to as R50 and RIO respectively. Transit time distribution (h(t)) in the ROI was calculated from the residue functions (R(t)) using, h(t) = -dR(t)/dt. Full-Width at Half-Maximum (FWHM) of the distribution of transit times was used as a measure of heterogeneity. Non- parametric Kolmogorov-Smirnov test was used to measure the statistical significance (p<0.05).

Here, only 50% and 10% of the maximum value has been used, however, it is to be understood that the time taken by the residue function to decrease by any other suitable percentage, may be used to parameterise the residue function curve. Results:

The residue function shapes estimated by CPI method were found to be smooth, whereas using the SVD approach, the residue function shapes were not smooth. Using CPI, differences between normal and ischemic tissues were observed: the CPI residue function in DWI+ and DWI-ring ROIs showed a slower decay compared to normal ROIs. The residue function shapes were consistent, with high MTT in corresponding tissues.

Figure 11 shows the R50 and R10 map for a representative patient along with CBF, MTT, delay and DWI. R50 and R10 maps showed appreciable variation between normal, DWI- ring and DWI+ regions.

Estimation of residue function shape may be informative in clinical situations allowing flow heterogeneity and tissue oxygen delivery to be assessed. Using the CPI approach, smooth residue function shapes were obtained that might facilitate transit time heterogeneity analysis which in turn could offer more insight into cerebral haemodynamics than CBF or MTT. In vivo ROI analysis with CPI method for five patients with atherosclerotic disease showed appreciable difference in shape of residue function and transit time distribution among normal, DWI-ring and DWI+ tissues. The transit time distribution appeared to have more contribution from high transit times in ischemic tissue compared to normal. R50 and R10 values were considerably higher for ischemic tissue during pre-surgery analysis which dropped significantly post-intervention. However a wide variation in values was observed in DWI+ tissue group which could be a reflection of variation within infarct tissue type. The R50 and R10 could serve as clinically useful quantitative measures for non-parametrically estimated residue function using CPI. R50 and RIO parametric maps might provide information in addition to standard perfusion and diffusion imaging.

The R50 and RIO values obtained using the CPI technique could provide potential useful information on capillary heterogeneity and possibly an access to study pathophysiology of stroke.

Further Bolus Dispersion Methods

The blocked circulant SVD, (oSVD), vascular model (VM), and control point interpolation (CPI) methods were used in the present investigation. Both the VM and CPI in their original implementation, were designed to model the residue function as per its fundamental definition with (R(0)=1). To explicitly accommodate dispersion, both methods were thus modified, as described below.

Dispersion Modification for VM

A gamma dispersion kernel was implemented in conjunction with VM, referred to here as

VM+VTF. The VM has three parameters, CBF, δ and λ, where δ is the bolus arrival delay and λ is the shape parameter of the gamma probability density function to describes the transit times through the voxel. VM+VTF has five free parameters with the addition of s and p to describe the gamma dispersion kernel. Dispersion Modification for CPI

In the original CPI, the residue function is estimated from a set of control points (CP) that form the basis of a smooth piecewise cubic spline interpolation. Each CP is allowed to vary in both amplitude and time, except for the first CP, which is set to 1 (CPi=l). This CPi restriction therefore does not permit assessment of dispersion. Two modifications were considered to accommodate dispersion in the CPI method. Figure 12 shows two implementations of CPI with dispersion: i) CPI₀: CPi fixed at 0, ii) CPI+VTF: CPI method implemented with the VTF. CPI+VTF attempts to characterize both the true underlying residue function (via the CPI residue function) and the dispersion effects (via the VTF), whereas, CPIo variant retains the 'model-free' nature of the original CPI since no assumptions are made for residue function shape or for dispersion kernel. The trade-off is that it does not allow to separate capillary haemodynamics from dispersion; thus neither the true residue function nor the dispersion kernel can be estimated separately.

Implementation

The delay insensitive variant of SVD, oSVD (blocked circulant SVD), was implemented. The VM was implemented with priors supplied from the oSVD solution, the model-fitting being performed using a Bayesian non-linear model fitting algorithm. The CPI was implemented using the same Bayesian non-linear model fitting algorithm using priors. The six deconvolution methods (oSVD, VM, CPI, VM+VTF, CPI₀ and CPI+VTF) were implemented in MATLAB using software written in-house. Of course other software may be used.

The Bayesian non-linear model fitting algorithm used for optimisation also required prior information for the VTF. Priors for the VTF were provided in the form of a Gaussian distribution with high probability for low-to-medium level of dispersion (based on empirical preliminary analysis, data not shown); priors used on p and s were log(2)+2 (mean+SD) and the same values were used for initialisation of the optimisation routine. Parameters p and s were estimated on a logarithmic scale, since this makes the model more linear with respect to these parameters, aiding the convergence of the fitting algorithm, which depends upon a local approximation to the function using a Taylor expansion.

Simulations

An initial set of simulations was performed with no dispersion. Next, a range of values of bolus dispersion and delay were introduced. In all cases, accuracy was calculated as the ability to estimate CBF and to separate the effects of dispersion from the simulated tissue haemodynamics (residue function).

The AIF was simulated using a gamma-variate function. Simulations were performed with Cerebral Blood Volume (CBV) = 4ml/100g (4%), representative of normal gray matter, and CBF in the range 10-60 ml/lOOg/min, in increments of 10 ml/lOOg/min. Thus MTT varied from 4s to 24s. A well-known model of the vasculature bed as a single well-mixed compartment was assumed to model tissue haemodynamics: the residue function was hence modelled as an exponential decay function. The Gamma VTF (Gamma Dispersion Kernel, GDK) was used to model dispersion; a range of low to high dispersion values were evaluated using the parameters shown in Table 5. Table 5: Gamma Dispersion Kernel (GDK) coefficients used in simulation for low-high dispersion in the MR signal.

Level of Dispersion p s

Low (L) Ϊ 2

Medium (M) 3 1

High (H) 5 0.5

The increase in the level of dispersion is associated with an increase in the time-to-peak, decrease in the signal drop, and broadening of the tissue signal; the area under the corresponding CTC remains the same.

Concentration time curves were then generated. A linear relationship was assumed between tissue relaxivity (R₂*) and C(t). The CTC was converted to signal curves S(t) with S(t0) at arbitrary set tolOO and TR/TE=\sl 65ms. The linear proportionality constant was chosen such that 40% peak signal drop was achieved at a flow of 60 ml/lOOg/min with CBV=4%. The AIF signal curve was generated similarly, but with a constant that produced a 60% peak signal drop. Zero mean Gaussian noise was added to the signal curves to produce a baseline SNR of 50. Values of 0-5 s were used to evaluate the models for sensitivity with respect to bolus arrival delay (δ) and its effect when combined with dispersion. The simulations were thus performed under three conditions:

I) no dispersion and no delay

II) dispersion with no delay: low, medium and high dispersion

III) dispersion with delay: low, medium and high dispersion and delay (5=0-5s)

For each combination of CBF, dispersion coefficients and delay a total of 100 CTC were generated. Thus a total of 14,400 (6 x 4 x 6 x 100) concentration curves were analyzed during each simulation process with six deconvolution methods: a) delay insensitive SVD (oSVD); b) vascular model (VM); c) VM+VTF; d) CPI (21); e) CPI₀; and f) CPI+VTF. In vivo Clinical Data

For empirical evaluation, the various methods were retrospectively applied to DSC data acquired on a Siemens 3T Trio using gradient-echo EPI: ?/ E=1.5sec/30msec, 78 volumes, 128x128x22 matrix, 1.7x1.7x5 mm voxels. An intra-venous (IV) bolus injection of 0.1 mmol/kg Magnevist was performed, followed by a 20 ml saline flush. DSC perfusion data were acquired from one healthy participant (age= 40yrs, male) and two patients (76yrs male and 65yrs female) with a history of atherosclerotic disease. Both of the patients had left internal carotid artery (ICA) stenosis, one with 75% and the other with 90% occlusion. Patients underwent left carotid endarterectomy (CEA) and DSC-MRI data were acquired both pre- and post-CEA. Images were acquired under an Institutional Review Board (IRB) approved protocol for DSC-MRI study as part of a larger study (31,32). Diffusion weighted imaging (DWI) (TR/TE = 4.4s/93ms, b-values

= 0, 1000 s/mm 2 , 27slices, 1.6xl.6x3mm 3 voxels) was also acquired along with the DSC data.

The AIF for the in vivo data was selected using an automated algorithm written in-house, from a 30mm slab centered at the Circle of Willis. CTC features including the first moment (FM), peak height (Pmax), time-to-peak (ttp) and bolus arrival time (Td) were extracted and, using a K-means clustering algorithm, voxels were grouped into five classes (grey matter, white matter, vessel, CSF and background). The AIF was defined in the vessel subclass with highest Pmax and lowest FM, and a venous output function (VOF) in the same class but with highest Pmax and greatest FM. The AIF was then scaled such that the mean value after the first passage of bolus (the tail of AIF) was equal to the tail of the VOF using the tail scaling technique.

Analysis of Simulated data

The ratio of CBF (CBFratio) as estimated CBF / true CBF was calculated to quantify the error in CBF estimation. The accuracy of residue function estimation for each method was evaluated by calculating the sum of squared errors (SSE) between the estimated and true residue function shapes. Statistical significance at p<0.05 was calculated using the Student paired t-test for comparison among methods.

Analysis of in-vivo data

In addition to CBF, MTT and CBV parametric maps, the residue functions obtained from deconvolution methods were also analysed for spatial variations in shape for the patients both pre- and post-CEA. For visualisation purposes, residue functions from two ROIs (2x2x1 voxels), one in infarcted and another in normal tissue, will be shown. Residue function is a monotonic decreasing function; thus, the time taken for the residue function to drop to 10% of its maximum value was calculated and referred to as RIO maps. Larger values of RIO would represent more time taken by the residue function to decrease, which might be an indication of tissue under haemodynamic stress. To see the spatial variations in dispersion, p and S (S=l/s) parametric maps were also analysed; larger values of p and S are likely to be associated with a higher amount of dispersion. Time to peak from the residue function from CPIo was also recorded and treated in the same way as the dispersion coefficient (p) in CPI+VTF as a measure of dispersion. Residue functions obtained from CPI and CPI+VTF method were analyzed for variation in shape and capillary heterogeneity. Residue function is a monotonic decreasing function; transit time distribution could be calculated from the observed residue function. Transit time distribution (h(t)) in the tissue voxel was calculated by the numerical differentiation of the residue function as in:

d R(t)

h (t) = - Equation 14

d t

where, h(t) is the tissue capillary transit time distribution with residue function R(t). Standard Deviation of the transit time distribution was used as a measure of capillary transit time heterogeneity (CTTH). Smaller values of CTTH represent a homogeneous distribution whereas a larger value might represent a broader distribution with more heterogeneity. Skewness (Skw) was also measured to quantify the asymmetry in the transit time distribution. Higher positive value indicates more Skewness of the curve towards right and lower positive value suggests less asymmetry. Negative values of Skewness would be suggestive of a left Skewed curve, which is unlikely in this case.

Hence, the following parameters are evaluated in this clinical data for quantification of cerebral hemodynamics using CPI and its dispersion corrected CPI+VTF method:

1) MTT

2) CTTH

3) Skw Results

Simulations

No dispersion and no delay

CPI (CBFratio = 1.00+0.12) was the most accurate in case of no dispersion, followed by VM (CBFratio = 0.95+0.11) which significantly underestimated CBF (p<0.05). CPI+VTF (CBFratio = 1.1+0.14) significantly overestimated CBF in case of no dispersion; whereas CPIo and VM+VTF overestimated CBF by approximately 1.6-2 times. Table 6 shows the SSE in estimation of residue function; the methods having the lowest SSE are marked in bold where (*) indicates statistical significance (p<0.05) and § indicates methods with statistically significant lower SSE within the groups but which could not be differentiated among themselves. In the case of no dispersion and no delay CPI had the lowest error, but the error with the CPI, VM and CPI+VTF methods showed no significant difference.

Table 6: Sum of squared error (SSE) for simulated and estimated residue function with oSVD, VM, VM+VTF, CPI and CPI+VTF at no, low, medium and high dispersion with and without delay. Lowest SSE are marked in bold.

Dispersion oSVD VM VM + VTF CPI CPI + VTF

No Delay

No 0.91+0.28 0.17+0.16^§ 1.13+0.90 0.14+0.18^§ 0.16+0.16^§

Low 1.91+0.35 0.24+0.20 0.83+0.89 0.22+0.34 0.16+0.17*

Medium 4.28+0.51 0.83+0.42 0.63+0.47 1.19+1.26 0.47+0.62* High 7.25+0.16 2.69+0.72 2.33+0.74* 6.19+3.83 3.18+2.39

With Delay

Low 4.77+0.49 0.26+0.24 0.32+0.46 0.56+0.66 0.16+0.21*

Medium 6.81+0.82 0.85+0.42 0.58+0.42^§ 1.95+1.38 0.58+0.79^§

High 9.18+1.18 2.71+0.07* 2.36+0.75 3.32+2.38 3.06+2.36 * indicate the method with statistically significant lowest SSE (p<0.05) § indicate methods with statistically significant lowest SSE (p<0.05) but are no better than each other

Dispersion with no Delay

As expected, the value of CBFratio decreased when dispersion levels went from low to high. oSVD, VM and CPI underestimated CBF in the presence of dispersion with poorer performance at the highest level of dispersion. CPI+VTF was most accurate for CBF estimation with a small standard deviation at low (0.99+0.13) and medium (0.81+0.13) dispersion values. CPIo tended to overestimate CBF at low and medium levels of dispersion but was most accurate compared to other methods at the highest level (1.04+0.06) of dispersion. VM+VTF was found to overestimate CBF at a low dispersion level by 1.6 times and to underestimate at medium and high dispersion levels. VM+VTF also had the highest level of uncertainty in CBF estimation among all other methods (VM+VTF: Low=1.57+0.87, Med=0.94+0.50, High=0.60+0.27). Table 6 shows the sum of squared errors in the estimation of residue function shape with a low, medium and high degree of simulated dispersion. At low and medium dispersion levels CPI+VTF had the lowest error in residue function shape estimation, significantly lower than other methods. The accuracy in estimation of residue function decreases with an increase in the level of dispersion for all methods. At the highest level of dispersion all the methods were found to perform poorly, among which VM+VTF showed the lowest error.

Dispersion with Delay

The methods that attempt to account for dispersion were found to perform better in the case of dispersion combined with delay. Introduction of delay together with dispersion led to an impaired performance by both VM and CPI, which was found to exacerbate as dispersion and/or delay was increased. CPI+VTF was found to be the most accurate in CBF estimation at low (CBFratio = 0.99+0.13) and medium (CBFratio = 0.80+0.14) levels of dispersion; whereas CPI₀ was more accurate at the highest (CBFratio = 1.03+0.06) level of dispersion simulated. oSVD always underestimated the CBF (CBFratio = 0.77+0.07 to 0.57+0.5) and VM+VTF overestimated CBF at low dispersion (CBFratio = 1.30+0.69) and tended to underestimate at medium and high (CBFratio = 0.85+0.35 and 0.56+0.05 respectively) dispersion levels. At the medium dispersion levels, the mean CBFratio for VM+VTF was found to be close to CPI+VTF, but VM+VTF had three times higher variability than CPI+VTF. Table 6 shows the sum of squared errors from the estimation of residue function at low, medium and high degrees of dispersion along with delay. At low dispersion, CPI+VTF had the significantly lowest SSE (p<0.05) but at the medium level of dispersion CPI+VTF and VM+VTF both performed similar in residue function shape estimation.

In vivo data

The oSVD estimates of CBF were the lowest among all methods; in agreement with the simulation study, VM+VTF consistently estimated the highest CBF within any given voxel. MTT values estimated by oSVD were observed to be higher than other methods, however a larger variation in MTT was observed with other methods.

Figure 13 shows the MTT estimates from a representative patient calculated with CPI and CPI+VTF; higher values of MTT being observed in the left parietal white matter region. The DWI confirmed a region of ischemic tissue in parietal white matter. CPI and CPI+VTF also showed higher values of MTT in some parts of the left MCA territory surrounding DWI positive ischemic tissue. Figure 13 also shows the residue functions from within normal and infarcted tissue from one patient. The residue function with CPI (Figure 13d) was broader and decayed more slowly in infarcted ROI compared to healthy ROI (Figure 13c). The residue function from CPI+VTF showed a faster decay compared to CPI in the infarcted ROI but was still broader in infarcted ROI (Figure 13d) compared to the normal ROI (Figure 13c). Residue function with CPIo, dispersed residue function, also showed broadening (higher dispersion) in infarcted ROI (Figure 13c) compared to the normal ROI (Figure 13d).

RIO maps for CPI and CPI+VTF, and a map of the sharpness parameter (S) from the gamma vascular transport function for CPI+VTF are shown in Figure 14 for the patient before and after CEA. The CPI method showed the area of infarction with very high RIO values (Figure 14b), along with increased values elsewhere in the left hemisphere mainly in the MCA territory. CPI+VTF showed similar areas of higher values in RIO (Figure 14c) but values in left MCA territory were much reduced compared to the higher values seen in CPI results. The dispersion coefficient (5) indicated a higher amount of dispersion in left MCA territory (arrow in Figure 14d). Post-CEA perfusion analysis showed decrease in RIO across the brain for both CPI (Figure 14f) and CPI+VTF (Figure 14g) and S (Figure 14h) with CPI+VTF methods.

Figure 15, from a representative patient, shows improvement in hemodynamics after CEA. Pre-surgery images showed high values of MTT in the right MCA territory which appears to decrease to normal (same as contralateral hemisphere) post-surgery. C7TH was found to be higher in right MCA region before surgery compared to contralateral left which also reduced to normal after surgery. Skw was lower in the ischemic region which improved to higher normal values same as left contralateral hemisphere after surgery. Similar changes in hemodynamics were observed in analysis with CPI and CPI+VTF. Right-Left hemispherical differences were qualitatively more clearly visible in pre-surgery images with CPI compared to CPI+VTF. Figure 16a illustrates the mean residue function observed with CPI method showing variation between normal tissue and tissue with hemodynamic disturbance (pathology). The residue function observed in pre-surgery evaluation from ischemic tissue consistently showed a slower decay compared to ones observed in normal tissue; residue function from DWI+ tissue decay even slower than the DWI-ring tissue. Figure 16b shows the corresponding transit time distribution for normal, DWI-ring and DWI+ tissue type during pre-surgery perfusion analysis; larger contribution from the higher transit times is seen in DWI-ring and DWI+ tissue. Figures 16c, d also show that the variation in residue function shape and transit time distribution among the three tissue types decreases post-surgery. The hemodynamics in DWI-ring and DWI+ tissue appears to change towards normal tissue type post-surgery.

Bolus dispersion is challenging to quantify in DSC-MRI, and more problematic in clinical population such as cerebrovascular patients. In the current work we developed a means to assess dispersion and how it contributes to error in: 1) the CBF estimation and 2) the tissue residue function. It was found that in presence of dispersion the accuracy in CBF estimation is significantly improved with dispersion corrected CPI methods and also allowed to separate tissue residue function from effects of dispersion with greater accuracy. CPI+VTF method corrected the residue function for dispersion and showed variation in tissue residue function in in vivo data from patients with atherosclerotic diseases.

Some investigations have shown that dispersion reduces the accuracy in CBF estimation when using the standard SVD methods by as much as 70%. Similar findings were observed in this study, where oSVD was found to underestimate CBF by a large amount (CBFratio=0.57- 0.77). However, modified methods designed to specifically handle dispersion, (VM+VTF, CPI+VTF and CPI₀), were more accurate in CBF estimation (CBFratio=0.94-1.04). CPI+VTF was most accurate in CBF estimation at low to medium levels of dispersion whereas CPIo was better at high level of dispersion.

It is already known that the broadening in the DSC time course caused by dispersion might be interpreted as an increase in the MTT if the analysis does not account for dispersion. This source of error propagates to the CBF estimation, leading to underestimation. In this case, the tissue response function would appear to be broader and to decay more slowly (an effect of apparent increase in MTT) compared to the true response function. Such effects were observed in our study when the analysis was performed with the original implementations of VM and CPI that cannot accommodate dispersion, with the estimated residue functions having slower decay when compared to those simulated. VM+VTF and CPI+VTF were found to perform better in such cases, as the dispersion can now be accommodated in the estimated vascular transport function.

For the in vivo data, the residue function estimated by CPI in the infarcted ROI was found to be broader and with a slower decay for the patient with atherosclerotic disease compared to a healthy scenario. This was also seen in the RIO map (Figure 14). The residue function when assessed with CPI+VTF in the same infarcted ROI decayed faster than the CPI estimate; this would imply the presence of dispersion, as suggested by the simulation study. Higher values of S (dispersion coefficient) in the left MCA territory (Figure 14) with CPI+VTF were further evidence of dispersion in this region. The CPI+VTF results suggest that, while there was a contribution from dispersion, there was still some variation in residue function associated with the infarcted tissue, as could be clearly seen by a broader residue function with CPI+VTF in an infarcted ROI compared to a more rapid decay in the normal ROI (Figure 13c,d). This is reasonable since changes in the microvasculature and thus the residue function would be expected in infarcted tissue. The wider region of a slower residue function decay observed in the left MCA territory with CPI was absent from the results obtained from CPI+VTF, implying that there were no substantial observable changes in the microvasculature in this region, but that the effect was primarily due to dispersive effects (Figure 14). It was also later observed in the post carotid endarterectomy perfusion analysis (Figure 14f-h), where these dispersive effects were not observed with either of the methods (CPI and CPI+VTF). In vivo neither the actual residue function shape nor effects of dispersion (VTF) are known a priori, particularly in pathological conditions. The CPI+VTF approach provides a method to avoid the strict model-based assumptions for the residue function with a plausible approximation of the dispersion kernel whilst still offering a smooth interpretable residue function estimate. Using CPI+VTF a better estimation of CBF was achieved even in the presence of dispersion in simulations. CPI+VTF was also found to isolate the effects of dispersion and to estimate cerebral haemodynamics at low and medium levels of dispersion, which represents the most likely clinical scenario. All methods performed poorly at extreme levels of high dispersion, suggesting that this may well present a limit to what can be accommodated in perfusion analysis. In order to account for dispersion effects which cannot be parameterised, a non-parametric deconvolution method accommodating dispersion has also been proposed in this study in the form of CPIo. CPIo performed within acceptable limits in CBF estimation but the residue function and dispersion remained inseparable. The estimated residue function has the combined information of the residue function and dispersion similar to the SVD solution, although with the benefit of no oscillations. It might be expected that the CBF estimation could be improved, and the residue function and effects of dispersion separated, if a model-based approach for tissue vasculature is used, parameterising both the residue function and effects of dispersion. However, the results of the present simulations suggest that a model-based deconvolution (VM+VTF) cannot clearly distinguish the two, leading to errors in CBF estimation and a significant amount of uncertainty in both CBF and the residue function shape. It was observed that variations in the signal time course from the VM+VTF model were similar for changes either in the dispersion parameters or in the parameters of the VM transit time distribution, most probably because both are defined by gamma distributions. Thus, there would be substantial ambiguity when estimating either set of parameters.

A metric based on the time for the residue function to decay to 10%, RIO, has been used here to characterise the residue function in vivo. This is in recognition of the fact that whilst it might be able to estimate the residue function its interpretation still needs to be fully explored. RIO in particular measures changes in the tail of the residue function and is thus associated with long transit times through the voxel. Thus it may capture information beyond that already found in the MTT parameter. The clinical significance of this parameter is beyond the scope of this work.

A gamma kernel was used to characterise dispersion effects. This is mere example of possible implementation. A number of other functions have been used in the prior art methods to model dispersion, for example, an exponential dispersion kernel has previously been used for DSC. In practice, the exponential kernel represents a special case of the family of gamma dispersion kernels, with p=0. Therefore, the gamma kernel provides a fairly flexible model for the VTF. This is mere example of possible implementation. Many variations and modifications may be made to the above-described embodiments without departing from the principles of the present disclosure. An alternative to dispersion correction is the use of methods that derive local AIFs, these would bypass the errors introduced using analytical dispersion kernels. However, they suffer from issues due to partial volume effects and there is currently no solution to this. Additionally, it still remains to be validated whether these methods can reliably measure the local AIR In summary, described here is a method to improve the accuracy of CBF estimates when dispersion is present. The method developed did not fully resolve the ambiguity between the true residue function and dispersion; however, incorporating an additional vascular transport function model did show benefits. A reasonable separation of tissue haemodynamics and bolus dispersion can now be achieved with this approach for low to medium levels of dispersion.

Therefore the method of the present invention offers a viable means to measure residue function shape and transit time distribution that has not been available with the prior art methods.

It should be emphasized that the above-described embodiments are merely examples of possible implementations. Many variations and modifications may be made to the above- described embodiments without departing from the principles of the present disclosure. All such modifications and variations are intended to be included herein within the scope of this disclosure and protected by the following claims.

Claims

1. A method for the quantification of cerebral haemodynamics, the method involving the deconvolution of the concentration time curve in tissue with an arterial input function (AIF), of the following equation:

C (t) = a^■ CBF · (Ca (t) <g) R (t))

where C(t) is the measured concentration time curve (CTC), <S> represents convolution, the proportionality constant a is a measure of brain tissue density and difference in hematocrit between capillaries and large vessels, CBF is the cerebral blood flow and R(t) is the residue function, the residue function describing the fraction of a tracer remaining in the tissue vasculature at a time t after its arrival, wherein Control Points (CPs) are used to parameterise said residue function, said Control Points being interpolated to produce a smooth continuous residue function curve.

2. The method of claim 1 , wherein the quantification of cerebral haemodynamics is from dynamic susceptibility contrast magnetic resonance imaging (MRI).

3. The method of claim 1 , wherein the optimisation method comprises a Bayesian Inference Scheme

4. The method of claim 1 , wherein the optimisation method comprises a non-linear least- square scheme.

5. The method of claim 1 wherein the CPs have degrees of freedom in both amplitude and time.

6. The method of claim 5, wherein each consecutive control point is related to its precursor by a ratio factor (Ω) and a time spacing (ζ), according to the following equation

where is the amplitude of the i"¹ CP and :

5

7. The method of claim 6, wherein the first CP is fixed such that ti=0 with r₁=l.

8. The method of claim 6 wherein the first CP is fixed at zero amplitude or allowed to vary in amplitude.

9. The method of claim 7, wherein t is fixed at the last time point of the measured CTC.

10. The method of any preceding claim, wherein said interpolation is piece wise cubic spine interpolation, for generating the complete shape of the residue function from said control points.

11. The method of claim 1 , wherein the CTC is converted to a Signal Estimate S(t) using the following equation:

TE

S (t) = S (t₀) - exp— ^c «

where S(t) is the, S(t₀) is the baseline signal, TE is the echo time and κ the proportionality constant.

12. The method of claim 1 wherein prior knowledge is incorporated for each parameter.

13. The method of claim 12, wherein the prior knowledge is in the form of a Gaussian distribution for each parameter with mean and standard deviation values.

14. The method of claim 1 , wherein effects of a bolus dispersion are included in a Vascular Transport Function (VTF) in convolution with the AIF, in the form of the following equation: C (t) = a · CBF · (Ca (t) <g) VTF(t) <g) R (t))

15. The method of claim 14, wherein the VTF is defined by the following equation:

sl+sp

Γ (1 + sp)

where, s characterizes the "sharpness" of the kernel and p is the time-to-peak (ttp).

16. The method of any preceding claim, wherein the time taken by the residue function to decrease to X% of its maximum value is used to further parameterise said residue function curve, wherein X is any percentage between 0 to 100.

17. The method of any of claims 1 to 15, wherein the time taken by the residue function to decrease to a proportion of its maximum value is used to further parameterize said residue function curve.