WO2008132469A1 - Estimation of a measure of contrast agent concentration using an analytical mathematical model - Google Patents

Abstract

A vascular input function, to estimate a measure of contrast agent concentration, which is decomposed into a function component representing the bolus of contrast-enhancing agent and a function component representing equilibration of the bolus by leakage from the circulatory system of the subject in to other parts of the body of the subject (referred to as the 'body transfer function'). Restricting the components to suitable analytical (e.g. algebraic) mathematical functions ensures that the calculations can be carried out analytically.

Description

ESTIMATION OF A MEASURE OF CONTRASTAGENT CONCENTRATION USING AN ANALYTICAL MATHEMATICAL MODEL

The present invention relates to methods and apparatus for estimating a measure relating to the amount or concentration of image contrast-enhancement agent within an imaged subject.

Radiological imaging of a subject may enable a determination of valuable physiological parameters when an image contract-enhancing agent is employed. Such an agent is a substance to be administered to a subject which, when the subject is subsequently imaged, appears as a change of the contrast (e.g. brightness) of the parts of the imaged subject at which the substance is present. By monitoring the temporal and/or spatial development of the concentration of the substance in a subject, via a sequence of radiological images over a time period, one may enable an estimate to be made of certain physiological parameters (often referred to as "kinetic parameters") of the subject.

Image contrast-enhancing agent is typically administered into the blood stream of a subject to be imaged. The quantity of agent immediately after administration takes the form of a relatively highly concentrated bolus of substance which disperses, or attenuates in concentration, over time as it is distributed through the circulatory system of the subject. The manner and mechanisms by which this attenuation of concentration occurs may be approximately represented, or modelled, mathematically using a "Vascular Input Function (VIF)". By using an appropriate VIF one may describe the temporal development of a measure of the concentration of the agent within a specified component of the subject's circulatory system. For example, applying the vascular input function to further mathematical processing, according to separate models of the manner in which the agent propagates from the circulatory system to other parts (e.g. tissues) of the subject's body, may provide kinetic parameters.

Obtaining accurate quantitative estimates of tissue kinetic parameters using Dynamic Contrast Enhanced radiological imaging (e.g. Magnetic Resonance Imaging (DCE-MRI) ) requires an accurate characterisation of the concentration of contrast agent in the blood plasma.

Many studies have demonstrated that inaccuracies in the vascular input function may propagate through to the estimated tissue parameters, and so it is important to accurately characterise the vascular input function. Using mathematically-defined models for the vascular input function has advantages over directly measuring agent concentration levels, as the estimation procedure is more stable, and parameter estimates tend to be more accurate if the vascular input function is appropriate.

However, models and parameter values proposed in current methods are often only applicable when there is a negligible proportion of plasma present in the imaged region of the subject under study (often referred to as the "plasma fraction") . When the plasma fraction is non- negligible such approaches tend to produce over-estimates of kinetic parameters (e.g. the volume-transfer constant trans

K ) . Furthermore, many existing methods are numerically intensive (e.g. requiring complex numerical processing of mathematical operations) and therefore slow and complex to implement.

The invention aims to provide an efficient and robust methodology that can be used in the case of fast temporal sampling of image data, and/or where a non-negligible plasma fraction is present.

At its most general, the invention proposed is to employ a vascular input function to estimate a measure of contrast agent concentration, which is decomposed into a function component representing the bolus of contrast- enhancing agent and a function component representing equilibration of the bolus by leakage from the circulatory system of the subject in to other parts of the body of the subject (referred to hereafter as the "body transfer function") . Restricting the components to suitable analytical (e.g. algebraic) mathematical functions ensures that the calculations can be carried out analytically. Even within these constraints, realistic representations can be produced. These may describe a peak in contrast-enhancing agent concentration at a region where the bolus first passes, observed after injection of a bolus of contrast agent. A subsequent recirculation peak may be described to represent concentration rises as the bolus passes a subsequent time, and the later phases.

The invention may implement a general framework for generating functional forms to define vascular input functions which are efficient to implement. By appropriately specifying the components of a vascular input function it is possible to generate functions that are realistic, and that ensure the analytical curves representing a concentration of contrast agent can be analytically calculated. This means that the computations necessary to estimate kinetic parameters from measured data are efficient (and, therefore, quick) , which is important when used in clinical practice.

The methodology may give analytic solutions for both the vascular input function and the representation of the concentration of contrast agent in tissue. The latter may be derived using the vascular input function. This avoids the need for costly numerical convolutions when fitting the functions or representations to measured data, and thus increases the potential applicability of quantitative radiology (e.g. DCE-MRI) in a clinical setting.

Accordingly, in a first of its aspects, the present invention may provide a method for estimating a measure of the concentration of an image contrast-enhancing agent within (e.g. the body or circulatory system) of an imaged subject including; providing an analytical mathematical function (e.g. a VIF) comprising function components and including a sum of a first function component (e.g. representing the shape of an injected bolus of a contrast agent) and a second function component in which the second function component includes a convolution of the first function component with a third function component (e.g. a body- transfer function) , wherein each function component includes one or more adjustable parameters; adjusting the numerical value of one or more adjustable parameters of the function to fit the function to a sequence of values of an image pixel (or pixels) representing a location in the imaged subject at each of a succession of times thereby to determine values for said parameters representative of the sequence; using the representative parameters to estimate a measure of a concentration of the image contrast- enhancing agent within the imaged subject (e.g. within the body, or circulatory system such as within blood plasma, and/or within tissue) .

Thus, the observed (imaged) temporal changes in the contrast (e.g. brightness) of a pixel (or group of pixels) representing a point or region of interest in a subject's body, due to the presence of an amount of contrast-enhancing agent, may used as the basis for estimating the optimal representative values for adjustable parameters of an analytical expression for the VIF. These values, and the analytical expression for the

VIF, may be used to estimate physiological (e.g.

"kinetic") parameters of the subject or any other measure of, or related to, a concentration of contrast agent

(e.g. the concentration itself) within the subject. The function may provide an analytical (e.g. algebraic) representation for a vascular input function which does not require lengthy computations to implement in practice. This is to be contrasted with existing purely numerical (or non-analytical) expressions for vascular input functions which require extensive (e.g. time- consuming) computation to implement (such as lengthy numerical integrations) in practice.

Furthermore, by requiring that the second component of the VIF includes a convolution of the first component, one may provide a methodology by which the first function component may be representative of the bolus of a contrast-enhancing agent ("contrast agent") and the third function component may be representative of the equilibration of the bolus by leakage of the agent from the circulatory system to other parts of the subject's body.

This enables a simple application and interpretation of the adjustable parameters of the function in use. The representative adjustable parameters may provide a measure of the amount of contrast agent in the bolus, the rate of attenuation/leakage of the agent in the bolus, and other measures representative of (or related to) concentration of contrast agent in an imaged region (e.g. blood plasma) . Furthermore, by intimately relating the first and second function components in this way one may help ensure that some adjustable parameters are common to both the first and second function components. This constraint enables the procedure of fitting the functions to image data to be more stable, robust and accurate. It greatly reduces the likelihood of unrealistic or unphysical parameter values arising from the fitting procedure, which may often occur in existing methodologies in which adjustable parameters are not sufficiently constrained.

The representative parameters, once determined, may themselves be used directly to provide a measure of (or relating to) a concentration of contrast agent within e.g. the blood plasma of the imaged subject at the imaged region. Additionally, the representative parameters may be used to determine an estimate of a concentration of the contrast agent in another part of the body of the imaged subject (e.g. within imaged tissue). This may be done by using the representative parameters in a model, expression or procedure for determining levels of concentration of contrast agent in tissue from an estimate of contrast agent concentration levels in plasma. The fitting procedure may include minimising the value of a measure of the difference between the values of the image pixel (s) in the sequence and the corresponding values of the function. The fitting procedure preferably aims to find values for the adjustable parameters which result in a function having a shape most closely matching the overall shape of the temporal sequence of the pixel values. This aims to enable the function to be used as a continuous analytical representation of the temporal development of the value of the image pixel (s) representing the location in the imaged subject in question.

It will be understood that the independent variable of the function is the time variable (t) such that the function may represent a time variation of contrast agent concentration, with a set of adjustable parameters.

The method may include providing the integral of the function, determining values for a cumulative sum of the sequence of pixel values, adjusting the numerical value of adjustable parameters of the function to fit the integral of the function to the cumulative sum. The method may include determining the representative parameters to be those which optimise the fit of the integral of the function to the cumulative sum. The fitting procedure may include a least-squares fitting of the integral of the function to a cumulative sum of the pixel values in the sequence.

It has been found that fitting the integral of the function to the cumulative sum of the sequence of values of a pixel tends to lead to much more stable and realistic representative parameter values. Noise in the image data is inherently smoothed by the process. Also, those parameters of the function which are representative of an injected bolus of contrast agent are forced to best match the cumulative sum of pixel value data rather than the individual values of the pixel themselves. This can be particularly beneficial when the sequence of values of the pixel contains (as is often the case) a few values (e.g. 2 or 3) much larger than others representing a brief surge in contrast agent concentration levels due to a passage of the bolus directly through, or close to, the imaged position. With such few data points being available to represent such a relatively significant data feature, standard curve-fitting methods often simply adjust the values of the adjustable parameters of the function to values to produce a good fit of the function to the few large pixel values, but in doing so often provide implausibly large values for times interpolated between the large pixel values (e.g. the function may "overshoot" between the large pixel values) . This typically results in representative values for the parameters which are unrealistic or even physiologically implausible .

The proposed method forces the fitting procedure to take account of all data and so apply an important constraint which is much more likely to result in accurate, meaningful and practically useful parameter values for use in analysing image data.

The integral of the function may most preferably be an analytical mathematical function (e.g. an algebraic function) since the function itself is analytical. This enables the fitting procedure to be performed quickly and accurately in avoiding the need to employ time-consuming or error-prone numerical integrations, which would otherwise be the case. Rather, merely the values of a relatively simple analytical expression need be evaluated.

The method may ^' include generating a cost function and adjusting the numerical value of adjustable parameters φ of the function to minimise the numerical value of the cost function, where the cost function includes the term χ{φ) given by

where

and C₁Q₁,φ) is the numerical value of the integral of the function for the time /, within the sequence, and y is the j^th pixel value in the sequence. The integral of the function is most preferably an analytical mathematical expression. The integral of the function (denoted cAt)) may be given by

and is preferably used/provided in analytical form.

The cost function may employ not only a cumulative sum of measured data and the integrated function, but may also employ the values of the differences between the pixel values and the corresponding function values, such as a sum (e.g. cumulative) of the sequence of such differences.

The cost function may be a weighted sum of a measure of difference between the cumulative sum of pixel values and the function integral, and a measure of difference between the pixel values and the function. For example, the cost function may be given by

The weights W₁ and w_p may be adjusted to adjust the relative importance of the two parts of the cost function.

The image pixel values may represent a location containing blood plasma. The method may include using the representative parameters to estimate a measure of the concentration of contrast-enhancing agent within blood plasma of the imaged subject at said location.

At least some of the image pixel values may represent a location within a bolus of contrast-enhancing agent. The first function component containing said representative parameters may represent a concentration of contrast- enhancing agent within the bolus .

The first function component may be arranged to so represent a smooth and continuous peak in a changing concentration of contrast agent at an imaged location. In this way, a rise, and subsequent fall, of a concentration of contrast agent may be represented by the first function component so as to represent an effect of a passage of the bolus through (or near to) the imaged region represented by the pixel (s) in respect of which the first function component is to be representative. In particular, the first function component may be arranged to represent the smooth, but often rapid, rise and fall in concentration due to the first-pass of the bolus through (or near to) an imaged location.

The first function component may include the term c_B(t) given by c_B(t) = a_Bt^a exp(—μ_Bt) where t represents time, a_B and μ_B are adjustable parameters, and a may have integer value e.g. a value of 1 (one), or any other integer (e.g. 2, 3, 4 or 5 etc), or may be zero.

The first function component may include the term c_B(t) given by

where t represents time, a_B and μ_B are adjustable parameters. The time variable is preferably limited to a range of values within one period of the term cos(μ_Bt) , e.g. mt_B ≤ t ≤ (m + \)t_B where t_B=2πμ_B ^~ , and m = 0, 1, 2, 3 ... etc. The third function component may include the term c_G(t) given by c_c(t) -a_G exp(-μ_Gt) where t represents time, a_G and μ_G are adjustable parameters. This function component may be representative of the body transfer function discussed above. The function may be given by c_p(t)= c_B(t)+ c_B(t)®G(t) where Θ represents a convolution operation.

The function c_p{t) and/or the representative parameters may be used to estimate a measure of the concentration c_t(f) of the image-contrast enhancing agent within tissue in the body of the imaged subject. For example, function and/or parameters may be used to describe the leakage of contrast agent into the extracellular-extravascular space (EES) of the imaged tissues. The function may be used as a component of a further function c_t(t) representative of the concentration of image-contrast enhancing agent within tissue in the body of the imaged subject. The further function may be of the form:

where h (t) is preferably an analytical mathematical function representative of the mechanism(s) for propagation of contrast agent within the imaged tissue (e.g. a tissue residue function). The parameter v_p may be an adjustable parameter. The function h (t) may contain adjustable parameters, and may take the form:

h(t) = K'^mseκp(rk_tpt) trans -1 where K (min ) may be representative of the volume transfer constant between the blood plasma and the EES,

-1 and k_ep (min ) may be representative of the rate constant between the EES and the blood plasma. The adjustable parameter v_P may be representative of the proportion of plasma present (the "plasma fraction") . The fraction of trans contrast agent in the EES may be given by v_e=K /k_ep. The method may include estimating any one or more of the above kinetic or physiological measures using the value of the associated representative parameter. The method

may include fitting the further function c_t(t — τ_o) to a

temporal sequence of values of an image pixel of imaged tissue of the subject (e.g. consecutive pixel brightness values associated with a tissue region undergoing changing levels of contrast agent concentration) in order to determine the representative parameter values (e.g. being those associated with the optimal ,fit of the

further function to observed data) . The parameter τ_o may have a value to represent the time of arrival of contrast in the subject at the location represented in the image

data by the image pixel in question. The value of τ_o may be determined by the fitting procedure with τ_o an adjustable parameter.

It will be appreciated that the method described above may be implemented on suitably arranged apparatus, such as apparatus including suitably programmed computer means, and the invention encompasses such apparatus.

In a another of its aspects, the invention may provide apparatus for estimating a measure of the concentration of an image contrast-enhancing agent within an imaged subject (e.g. within the body) including computer means for providing a representation of an analytical mathematical function comprising function components and including a sum of a first function component and a second function component in which the second function component includes a convolution of the first function component with a third function component, wherein each function component includes one or more adjustable parameters, and for receiving a sequence of values of an image pixel representing a location in the imaged subject at each of a succession of times; the computer means being arranged to adjust the numerical value of adjustable parameters of the function to fit the function to the sequence thereby to determine values for said parameters representative of the sequence; the computer means being further arranged to use the representative parameters to estimate a measure of the concentration of the image-contrast enhancing agent within (e.g. the body of) the imaged subject.

The computer ^'means may be arranged to provide a representation of the integral of the function, to determine values for a cumulative sum of the sequence of pixel values, to adjust the numerical value of adjustable parameters of the function to fit the integral of the function to the cumulative sum. The computer means may be arranged to determine the representative parameters to be those which optimise the fit of the integral of the function to the cumulative sum.

The computer means may be arranged to generate a cost function and to adjust the numerical value of adjustable parameters φ of the function to minimise the numerical value of the cost function, where the cost function includes the term χ(φ) given by

where

and C,(t_t,φ) is the numerical value of the integral of the function for the time t, within the sequence, and y_} is the j^th pixel value in the sequence. The integral of the function c_p(t) may be given by

and the computer means may provide a representation of the integral in analytical form. The computer means may be arranged to determine the cost function as a weighted sum of a measure of difference between the cumulative sum of pixel values and the function integral, and a measure of difference between the pixel values and the function. For example, the cost function may be determined as

The weights W₁ and w_p may be adjusted by the computer (or user) to adjust the relative importance of the two parts of the cost function.

The integral of the function is preferably an analytical mathematical expression.

The image pixel values may represent a location containing blood plasma. The computer means may be arranged to use the representative parameters to estimate a measure of the concentration of contrast-enhancing agent within blood plasma of the imaged subject at said location.

At least some of the image pixel values may represent a location within a bolus of contrast-enhancing agent. The first function component containing said representative parameters may represent a concentration of contrast enhancement agent within the bolus.

The first function component may include the term c_B(t) given by c_B(t) = a_Bt^a exp(—μ_Bt) where t represents time, a_B and μ_B are adjustable parameters, and a may have an integer value e.g. of 1 (one), or any other integer (e.g. 2, 3, 4 or 5 etc) , or may be zero.

The first function component may include the term c_B(t) given by c_B(t) = a_B(l-cos(μ_Bt)) where t represents time, a_B and μ_B are adjustable parameters. The time variable is preferably limited to a range of values within one period of the term cos(//_βt) , e.g. mt_B ≤ t ≤ {m + l)t_B where t_B=2πμ_B ^Λ, and m = 0, 1, 2, 3 ... etc. The third function component may include the term c_G(t) given by c_G(t) = a_G exp(-μ_Gt) where t represents time, a_G and μ_G are adjustable parameters.

The function may be given by c_p(t)= c_B(t)+ c_B(t)®G(t) where <8> represents a convolution operation. The computer means may be arranged to use the function c_p(t) and/or the representative parameters to estimate a measure of the concentration c_t{t) of the image-contrast enhancing agent within tissue in the body of the imaged subject. For example, the computer means may be arranged to enable the function and/or parameters may be used to describe the leakage of contrast agent into the extracellular- extravascular space (EES) of the imaged tissues. The computer means may be arranged to provide a representation of a further function c_t(t) , and may be arranged to use the function as a component of the further function. The computer means may be arranged to use the further function c,(t) in a representation of the concentration of image-contrast enhancing agent within tissue in the body of the imaged subject. The further function may be of the form: c,(0 ="^,(0+^(0®A(O where h (t) is preferably an analytical mathematical function representative of the mechanism (s) for propagation of contrast agent within the imaged tissue (e.g. a tissue residue function) . The computer means may

be arranged to provide the parameter v_p may be an

adjustable parameter. The function h (t) may contain adjustable parameters, and may take the form:

trans 1 where K (min^~ ) may be provided as representative of the volume transfer constant between the blood plasma and

-1 the EES, and k_ep (min ) may be provided as representative of the rate constant between the EES and the blood plasma. The computer means may be arranged to provide a value of the adjustable parameter v_P which may be representative of the proportion of plasma present (the "plasma fraction") . The computer means may be arranged to determine an estimate of the fraction of contrast agent trans in the EES by determining v_e=K /k_ep. The computer means may be arranged to estimate any one or more of the above kinetic or physiological measures using the values of the associated representative parameters. The computer means

may be arranged to fit the further function c_t(t) to a

temporal sequence of image pixel values of imaged tissue of the subject (e.g. consecutive pixel brightness values associated with a tissue region undergoing changing levels of contrast agent concentration) in order to determine the representative parameter values (e.g. being those associated with the optimal fit of the further function to observed data) .

In any of the above aspects, the pixel values to which functional forms may be fitted may be determined according to the acquired radiological (e.g. NMR) signal intensity converted to contrast agent concentration using the method of Wang et al . Magn. Reson. Med. 5(5), (1987) pp399-416. Methods other than that of Wang et al may be used for this purpose. The aspects described herein are suitable for use in analysis of images acquired by nuclear magnetic resonance imaging, CT imaging, PET imaging or any other image acquisition method in which image contrast-enhancing agents may be' used.

In any of its aspects, the third function component (e.g. BTF) may be given by:

where i,j=l, 2, 3... denotes the i^th or j^th term in the sum of terms, and a^ , τ£⁽ and /4° are parameters adjustable to be representative of recirculation of the bolus, and

OQ and μ^ are parameters adjustable to be representative of equilibration, renal excretion or any other processes. Any of the amplitude terms a^ or a^ may be set to zero as desired.

There now follow examples of the invention described with reference to the accompanying drawings of which:

Figure 1 illustrates graphically the form of VIFs (left panel) , and both the integral of the VIFs over time (right panel) ; Figure 2 illustrates graphically the form of VIFs (left panel) , and both the integral of an analytical VIF together with a cumulative sum of pixel values (right panel) ;

Figure 3 illustrates four tissue kinetic parameter maps for image data of a subject containing a primary bladder tumour;

Figure 4 illustrates apparatus for determining kinetic parameters from image pixel data.

Vascular Input Function:

After a bolus of image contrast-enhancing agent

("contrast agent") is injected in to a peripheral vein of a subject to be imaged, the bolus passes through the cardio-pulmonary system. The bolus is then transported to the arterial side of all organs, including the imaged region, so that an element of the shape of the required input function is the bolus shape as it enters the imaged region, denoted c_B(t) . As the bolus travels through the tissues, some contrast agent generally leaks from the bolus, and the network of capillaries tends to delay and disperse the bolus.

This process mixes the bolus of contrast agent ^• with the entire blood pool due to continual recirculation. Transfer of contrast agent from the plasma into the whole body leakage space eventually leads to an equilibrium concentration in the blood plasma. These processes of recirculation, mixing and leakage predominantly occur over timescales below 10-20 minutes. At the same time, renal excretion permanently removes contrast agent from the plasma pool. However, this process generally happens more slowly, e.g. at a rate equivalent to a concentration attenuation with a half life of around an hour. All these effects are represented by a body transfer function (BTF) denoted G(t). The function driving this process is the shape of the bolus of contrast agent for the imaged region, c_B(t). The concentration of contrast agent in the blood plasma c_P{t) may then be represented by a superposition of the bolus shape and its shape after some modification by the body transfer function, that is:

c_p(t)= c_B(t)+ c_B(t)<8>G(t) (1)

where <8> represents a convolution operation.

Three different examples for the vascular input function

(Cp(t)) are described in detail below, and of great importance is that all three give relatively simple analytic forms for the plasma concentration, as are the resulting functions representing tissue concentration derived. from the VIFs. A brief discussion of more complex models that have greater realism, while maintaining analytic tractability is also given.

VIF #1:

Consider using c_B{t)= a_Bexp(-μ_Bt) and G(t)= a_Gexp(-μ_Gt) . This describes the bolus with an exponential, and the BTF can be interpreted as a two-compartment model (the blood plasma and extracellular-extravascular space (EES) of the whole body) with transfer rate μ_G. The resulting input function is:

Cp(0 = A_B exp(-/V) + A_G exp(-/V) ( 2 ) with A_B = a_B-a_Ba_G/ ( μ_B-μ_G) and AQ = a_Ba_G/ ( μ_B-μ_G) •

VIF #2 :

Another choice is to use c_B(t)= a_Btexp(-μ_Bt) for the bolus model, with the property that c_B(0) = 0, so the input function is continuous at t = 0. With G(t)= a_Gexp(-//_Gt) for the BTF this leads to:

with A_B = a_B-a_Ba_G/ { μ_B-μ_c) and AQ = a_Ba_G/ { μ_B-μ_G) .

VIF #3 :

A limitation of the previous model is that although the curve is continuous, the gradient is not continuous at t = 0. A solution is to use a raised cosine to represent the shape of the bolus, of the form:

With G(t) = a_Gexp(-μ_Gt) for the BTF this leads to:

where f { t, μ_G) is defined by :

The components of the function are relatively simple, and the resulting equations only involve simple mathematical functions that are easy to compute.

Other VIFs: The three VIFs (c_P(t)) described above have the advantage that they have relatively simple algebraic forms. A generalisation of the first function component c_B(t) describing the shape of the bolus is to use a gamma- a variate function of the form c_B(t)= a_Bt exp(-μ_Bt) . The additional α term gives increased flexibility in choosing the representation of the bolus shape, particularly the initial rise of the contrast agent concentration through the bolus. Algebraic solutions exist for c_P(t) when the term a is an integer, and these involve sums containing at least (α+1) terms to be summed, while for general values (including non-integer) of the term α, the solution for c_P{t) involves using special mathematical functions, which are available as (computationally expensive) library functions. The advantage of the raised cosine in VIF #3 is to achieve a similar bolus shape to this model with α>l, but to avoid the need for large summations or special functions.

Generalisations of the BTF can be used to model the other processes mentioned above.

Examples of analytical mathematical functions are described below in terns of cosine bolus models and gamma- function bolus models, such models may use the following step function explicitly wherever it is suitable to do so:

According to a raised cosine bolus model, the bolus is modelled with one complete cycle of a raised cosine function, as follows: B(t;μ) = (1-cosOO)^•(s(t)-s(t-t_B)), where t_B = 2πμ^~λ .

The convolution of a bolus model and an exponential function may serve as the main template function used to construct the input function and tissue function. It can be evaluated for small t and with κ = 0 , which is ensured by the form used for the helper function /|() given below:

The function /JQ is given by:

This form is preferable to avoid problems with numerical underflow. The first case for /J(-) listed above is the basic definition of the function, and the other three cases listed below the first case are formed from combinations of the following taylor-series

approximations , The

threshold ε, is preferably set according to the numerical precision being used for the computations, and for floating point calculations (e.g. using IDL), e,= 2xlO^~2 is appropriate. When κ = 0 and |μ/|>£_j, the second listed of the above-listed cases of the expression for _/j(-) may be used. When this is evaluated it may be reduced to t—μ^~ι sin(μt) , which will give the correct numerical result.

The convolution of a raised cosine model with gamma form may serve as an alternative template function used to construct the tissue function when one of the input function parameters or is (or is substantially) equal to one of the leakage parameters. The convulusion is:

10 where

The function /₂(-) is given by :

25 This form is preferable to avoid problems with numerical underflow, and mirrors the definition of _/J(-) . For floating point calculations (e.g. using IDL), ε₂=8xl(T² is appropriate. This function also reduces to the correct form when K = 0. A St. Lawrence & Lee model may be employed whereby the input function is comprised as follows.

The bolus model may be c_B(t) = a_BB(t;μ_B) and the body transfer

function may be G(t) = a_Ge ^μ° -s(t) , so the input function is given by

c_p(t) = c_B{t) +c_B{t)®G(O = a_BB(t;μ_B)+a_Ba_cE(t;μ_B,μ_c).

For the St Lawrence & Lee model the tissue residue function may be given by

R(t) = F-(s(t)-s(t-T_c))+E_fF-c^'k'^Λ'-^τ^ ._S(t-T_c),

Where F is an adjustable parameter quantifying blood flow rate (ml/lOOg/min) , E_f is a parameter quantifying the extraction fraction (no units), k_ep is the return rate constant (min^"1) and T_c is the capilliary transit time (min) . The tissue residue function can be considered as a sum of three exponentials with rates 0 , 0 and k respectively and delays 0 , T_c and T_c . In general, if the tissue residue function is R(t) , then the equation for the tissue curve may be given by

c_t(t) = c_p{t)®R{t) =c_B(t)®R{t)+c_B(t)®G{t)®R(t).

The first term in the equation for the tissue curve can be written:

c_B(t)®R(t) = a_BF-E(t;μ_B,O)-a_BF-E(t-T_c;μ_B,O)+a_BE_fF-E(t-T_c;μ_B,k).

The second term of the equation for the tissue curve can be considered in two stages.

The second part of the second term in the equation for the tissue curve can be written:

The whole of the second term in the equation for the tissue curve can then be written:

These terms may be combined, and similar template functions collected together to give:

By writing the function in this way it is possible to determine when numerical underflow may occur. The first two lines compute the difference between two template functions whose arguments differ only by a delay of T_c . If T_c is very small then this will cause underflow problems, which can be avoided by placing a lower-bound on T_c . For floating point calculations (e.g. IDL) a lower limit of lCr³ is appropriate. The third line will experience underflow problems when k_ep∞μ_G. These terms are generated

by the convolution of e ""' and e ^μ°' which come from the tissue residue function and the body transfer function respectively. If μ_G ⁼k_ep then the convolution of these

terms gives te ^v , which can be used to provide an alternative form that avoids numerical underflow. By transforming the parameters to the dimensionless quantities μ_ct and k_ept , it can be shown that the appropriate threshold is

, when the third line should be replaced with

-a_Ba_GE_fF■D(t -T_c\μ_B,k_ep) .

The following describes calculations involving a gamma function bolus representation according to an example of the invention.

The bolus model is a gamma variate function with integer exponent (the integer exponent m is analogous to the integer exponent α referred to elsewhere herein) ,

Convolution of this with an exponential has an analytic form,

Computation of E{t;μ,κ) may be performed as follows.

To accurately compute this function for a range of input parameters it is preferable to construct various alternative functions that will compute correctly for large or small values. To reduce the effect of large values one may initially compute the logarithm of the function, where for t>0 we obtain

logE(t;μ,κ) = log(m!)+(m+1)log(t)-μt+f_m(t(μ-K)).

The helper function f_m(-) is defined by

and the ancillary parameters ε_m, L_1n and Z_1n are defined in the following table. The values in this table are designed to give a result that has a fractional accuracy of 10^~6 , and the result may be accurate over a large dynamic range of t, μ and K. Importantly, the special case μ=κ may be accurately computed.

computation of E(t;μ,κ)

The form of the helper function is preferred to avoid problems with numerical underflow, overflow and division by zero. A particularly important case when this could occur is with μ = κ, so z = 0. In this case /_m(-) may be evaluated via

, so the final computation may be correctly computed without division by zero. Another case when z = 0 is when t = 0, but this is easily handled as the final result can be directly set to zero.

The thresholds ε_m are preferably chosen to ensure that the direct form of the equation, α_m(z) , is only used when the result will not suffer from underflow. This is caused by the difference between e^∑ and the summation in the expression for ot_m(z) underflowing the machine precision, which happens when IzI is small. The summation limits L are preferably chosen to ensure the alternative form β_m(z) uses sufficient terms to be accurate -- the exact, but impractical choice is L_1n = oo . The upper limits Z_m are preferable so that the e^z term in oc_m(z) does not overflow for large positive z , when the alternative 7_m(z) expression is used. For large negative z , the expression for ^a _m(^z) ^may preferaby be used, and although the computation of e^z may suffer from underflow, the result may still be accurate.

Firstly one may preferably choose the thresholds to ensure ce_m(z) is accurately computed. Underflow may affect the computation of the second logarithm term in the definition of ot_m(z) , which for large positive z tends towards e^z , and for large negative z tends towards z^mlm\ . These limiting functions may be accurately computed for suitably large \z\, so thresholds ε_m preferably exist such that oc_m{z) is accurately computed for \z\>ε_m. To choose the thresholds note that for |z|<l,

The approximation may be justified because when |z|<l, |z"/n!| is a strictly decreasing series for all n>0, and the higher order terms in the infinite summation may be neglected. Numerical underflow can occur when the result of the calculation is small with respect to the two subtracted terms, that is when | z\^m+x /(m+l)!< e^z (or the summation). This may occur when |z|<l, which implies e^'ssl, so numerical underflow may be identified simply by |z|^m+I /(m+1)!. The thresholds can be defined by

, where ε_f is preferably chosen to be larger than the floating-point precision of the computations. The following table uses ε_/=l(r⁵ , and gives thresholds and values of a_m(z) for z = ±ε_m.

Since all ε_m are less than 1 the decreasing property of the series given earlier holds. The computed values of α.(±εj are all close to ±ε_f as intended, so these may all be accurately computed if the machine precision is better than ε_f . The modulus operations in the definition of cc_m(-) are preferable because μ<κ gives z<0, so the logarithms may otherwise not be correctly computed. However, note that

since the argument of the final logarithm is always positive for all real z .

With these thresholds one may define upper limits for the summation in β_m(z) that prefereably ensure the truncation error is negligible. Since all the thresholds are less than unity, the summation series is decreasing, so L_1n may be defined to ensure that the largest truncated term of the series is smaller than the largest term in the series by an appropriate fraction. Thus;

for \z\<ε_m, where ε_t is the chosen truncation error. The left-hand side (LHS) of this inequality is largest when Z | == ε_m , so using this value of z in the expression, and substituting ε_f , this equation becomes

One may tabulate the LHS for different L_1n and select the smallest L_1n that satisfies the inequality for a given ε, .

The following table gives values of

as before. For ε, =10^~6, tabulated values below 10^~" are therefore preferably deemed as having acceptable L_1n .

The upper thresholds Z_m may be be derived by finding solutions to the equation

as this is the fractional truncation error committed by neglecting the summation in the expression for &_m{z) . As before one may employ ε_r=10^~6, and the solutions are given in table 1. A formula for the computation of the convolution of the bolus with te " is .

To accurately compute D{t;μ,κ) for a range of input parameters it is preferable to construct various alternative functions that will compute correctly for large or small values. To reduce the effect of large values one may compute the logarithm of the function, where for t>0 we obtain

logD(t;μ,κ) = log(m!)+(m+2)log(t) -μt+g_m(t(μ-K)).

The helper function g¹ (•) is defined by

and the ancillary parameters are defined in the following table. The values in this table are designed to give a result that has a fractional accuracy of 10^~6 , and the result is accurate over a large dynamic range of t, μ and K. Importantly, the special case μ = κ will be accurately computed.

Table 2: Values of ancillary parameters for the computation of G(t;μ,κ)

The thresholds and summation limits are deduced in the same way as described above. The bracketed term in the expression for <5_m(-) can be manipulated to give

The terms in the right-hand side (RHS) have decreasing magnitude when |z|<l since /n>0 and n>m+2. So for |z|<l this expression can be approximated by the largest term in the series, that is z^m+2/(w+2)! . Underflow problems may occur when the result is much smaller than the subtracted terms, i.e. when Since this may

occur when the thresholds can be defined via

and this results in the following thresholds.

These thresholds are all less than one, so the decreasing property of the series holds meaning that the various approximations are appropriate. From the design of these thresholds 17_m(±ε_m) \∞(m+\)ε_f , and the above table demonstrates that this approximation is accurate. With these thresholds it is possible to specify the summation limits in ξ_m(z) . The ratio of the largest term to the largest truncated term in cω IS

and if this ratio is less

than ε, at z = ε_m, then this is equivalent to

{L_m-m)ε_m £^m_+!l(L_m-\-\)\<ε_tε_f . The following table gives the LHS

of this inequality for

and so tabulated values below 10 give acceptable L_1n .

The upper thresholds Z_m can be derived by finding solutions to the equation

10 as this is the fractional truncation error committed by neglecting the summation in the expression for δ_m{z) . As before one may use ε, =10^~6, and the solutions are given in

table 1.

Contrast Agent Concentration in Tissue:

With the VIF provided in a representative analytical form, it may be used to generate a measure of a concentration of a contrast agent in a tissue of the imaged subject.

The extended-Kety model (Kety S. Pharmacol . Rev. 3,

(1951), ppl-41; Tofts P. J. Magn . Reson. Imag. 7, (1997), pp91-101) may be used to describe the leakage of contrast agent into the extracellular-extravascular space (EES) of the imaged tissues, and takes the form

c, (/) = v_pc_p (0 + c_p (0 <8> {K'^rαm exp(-k_ept)}

trans 1 where K (min^~ ) is the volume transfer constant between

the blood plasma and the EES, k_ep (min ) is the rate constant between the EES and the blood plasma and v_P is the proportion of plasma present (the "plasma fraction") . The fraction of contrast agent in the EES is given by trans v_e=K /k_ep < 1. For the input functions described above, the tissue concentration is given explicitly for VIF #1 by:

The equations for VIF #2 are given by: c Xt) = v_p[A_BtQxp{-μ_Bt) + 4_;(exp(-/i_G0 - exp(-μ_Bt))]

The equations for VIF #3 are given by (equation 8):

and, for t>t_E

where

and where c_p(t) is given by equation (4) . Special forms of these expressions which arise when k_ep μ_B or μ_G are given by the following.

When k_ep = μ_B then:

VIF #1:

VIF #2:

VIF #3: Unchanged.

When k_ep = μ_G then: VIF #1:

VIF #2:

VIF # 3 :

VIF Evaluation :

VIF #2 and #3 are more appealing than VI F # 1 because for both the initial phase of the bolus rises from zero , rather than j umping abruptly ( see figure 1 ) . However , it is instructive to quantitatively determine how this difference affects the tissue parameter estimates , and this is done by means of a simulation-based experiment .

Evaluation via Simulated Data:

Data were simulated using a population-based input function, obtained from a cohort of patients undergoing DCE-MRI examinations (Parker et al. Magn. Reson. Med. 45, (2006), pp993-1000) shown in figure 1. Test data were generated from these curves using: y. = c,(nT- t_o) + ε_H

for n =0, 1, ..., N, where c_t(...) is the simulated tissue curve, T is the sampling interval, to is the bolus arrival

time, and ε_n is a Gaussian-distributed random sample with

2 variance σ . Eight simulation scenarios were considered, trans , using combinations of K e{0.1, 0.4} min , v_ee{0.2, 0.6} and v_pe{0, 0.05}, and 2000 data sets were simulated for each case. The other parameters were σ =0.02 mM, T=3/60 min, N=60, and to=13.5/6O min, which were chosen to reflect those in a typical fast imaging sequence capable of observing the first-pass of an injected bolus.

For each data set, estimates of the tissue parameters were obtained for each of the proposed vascular input functions using equations (6), (7) and (8) in a least- squares fitting routine, which also included the onset time to as a fit parameter. The value of the parameter v_p was constrained within the fitting routine to be above 0 (zero) since this parameter cannot physically be

2 negative. In each case a p-value was computed using a χ - statistic of the form:

where y_n are the predicted measurements using the least-

squares estimate, and there are N-4 degrees of freedom.

2 If the residuals are genuinely Gaussian with variance σ , then the expected p-value is 0.5, and typically a threshold of p<0.001 is used to reject implausible VIFs. Parameters for the proposed input functions were obtained by least-squares fitting the time-integral of the population input function to the time-integral of each of the proposed input function models over the same range as the test data, that is te [0, 3] . The time-integral of the population input function was obtained by numerical integration, and the corresponding integrals for the proposed VIFs by analytic integration of equations (2), (3) and (4) . Included in the fit was an additional time- offset parameter to account for the fact that the peak of the bolus shape in the input function is not at t=0. These curves are shown in figure 1, along with the time- integrals, and the estimated parameter values detailed in table 1. The table also details the corresponding amplitude parameters (a_B and a_G) for the first function component (bolus shape) and the third function component (BTF) of the two models, a_B and a_G. The BTF parameters A₆ and μ_G are very similar for all three models, which is not surprising since the form of the BTF is the same in each case. The total area under the bolus is a_B/μB=0.863 for VIF #1, a_B/(μ_B) =0.848 for VIF #2 and a_B/μ_B=0.782 for VIF #3. If the first function component (c_B(t)) is interpreted as a distribution over arrival times, then

-1 the mean arrival time is μ + t₀ = 0.187 for VIF #1,

2μ + to = 0.179 for VIF #2, and 0.5μ + t0 = 0.141 for VIF #3, which again are similar. These similarities are reflected in the very close agreement between ^' the time- integrals as shown in the right-hand panel of figure 1.

In figure 1, in both panels the black dash-dot line is the population derived input function, the gray dashed line is VIF #1, the solid black line is VIF #2, and the solid gray line is VIF #3. The left panel shows the curves c_p(t), the right panel shows the time-integral of the curves, and in both panels all three are virtually identical after t=l sec.

Table 1. Estimated Input Function Parameters, including to, the time-offset for the fit. Note that the units on a_B and A_B are different because of the different model forms.

The overall objective of fitting the VIF to the image pixel sequence is to fit a given VIF (vascular input function) to noisy data. The fitted curves and pixel data each consist of an initial period with zero concentration followed by a short peak (due to the bolus) and a longer period of gradually decreasing concentration, as shown in figures 1 and 2.

In principle, standard nonlinear least-squares fitting methods can be used to give estimates of the model parameters by directly fitting the VIF curve to the data. However, with current time-resolution limits on the data acquisition hardware (particularly with MR imaging systems) , the bolus is often only observed at one or two time points. Since the VIF model typically has two (or more) parameters to define the bolus, if there are two data points for the bolus there will be a parameter set that matches the data. However, the resulting curve is often implausible - usually the peak of the curve comes between the two data points, and is much too high.

The solution employed here is to fit the integral of the VIF model curve to the cumulative sum of the pixel data using a least-squares method. This tends to give more plausible results because the fitting is driven by the area and duration of the part of the data dominated by the bolus, rather than the amplitude and duration of that part. The procedure is mathematically described below.

If c_p{t,φ) is the "true" VIF parameterised by φ, then the data y_n acquired at times t_n are modelled with

where ε_n are (unknown) error terms and n = 1, 2, ... N. The cumulative sums of the data are given by

for n = 1 , 2 , ... N, and the integrated VI F is given by

Depending on the choice of c_p{t,φ), C,(t,φ) may be available analytically. The fitting process uses a cost function of the form

Standard algorithms are used to find the φ that minimises χ(φ) . Another possibility is to use a cost function that combines both the original data and the cumulative data, that is

The weights W₁ and w_p are used to adjust the relative importance of the two parts of the cost function.

Simulation Results:

Table 2 shows mean values for the parameter estimates over the 2000 data sets for each scenario described in the previous section. Uncertainties are given as plus/minus two standard deviations, and the average p-

2 value for the χ -statistic is given for each case. For the cases with v_p=0 the parameter estimates have very little bias, while the average p-values indicate that the VIF fits are indistinguishable from what one would expect with the true VIF (i.e. using the population input function).

Table 2. Mean parameter estimates and p-values for eight test scenarios. Uncertainties are plus/minis two standard deviations .

When Vp=O.05, vif #1 gives biassed estimates, particularly for V_p, and the low p-values indicate that in this context a bi-exponential input function is distinguishable from the population input function. This is not surprising since the shape of the bi-exponential input function is very different from the population input function. However, for VIF #2 the bias is much reduced - the two standard deviation interval contains the true value in every case. The average p-value is distinctly less than 0.5, but is sufficiently large to suggest that statistically, VIF #2 is little different from the population input function in this context. For VIF #3 the bias is similar to VIF #2, though the variance is generally slightly smaller. However, the p-values suggest that VIF #3 is statistically indistinguishable from the population input function in this context.

The average time taken to compute the least-squares estimates of the tissue curves (c_t(t)) for the three VIFs was also measured - VIF #1 was the fastest, VIF #2 took 10% longer than VIF #1 and VIF # 3 took 50% longer than VIF #1. However, the average time taken to find the least- squared estimates with v_p =0.05 was fastest for VIF #2, VIF #3 took 3% longer than VIF #2, and VIF #1 took 30% longer than VIF #2. This seemingly paradoxical result may be because VIF #1 is not very similar to the true model, so the least-squares cost function is likely to be non- quadratic and therefore requires more iterations to find the minimum. For VIF #3, since it is very similar to the true VIF, the least-squares cost function is likely to be nearly quadratic, and so the minimum will be found in very few iterations, which offsets the additional time required to compute the tissue curves.

In-v±vo Example :

The simulations suggest that VIF #2 is a reasonable compromise between statistical accuracy and model complexity. Therefore, an example is given of the application of VIF #2 to some in-vivo data taken from a patient with a bladder carcinoma. The dynamic data for this study were obtained using a Siemens Avanto, and consisted of 70 dynamic measurements acquired every 5.6s using a 3D spoiled gradient-echo sequence during shallow breathing. Magnevist contrast agent (relaxivity of 4.26

-1 -1 s mM ) was administered using a power injector with a dose of O.lmmol/kg body weight, and the imaging parameters were TR/TE = 4.36/1.34ms, α = 24°, 256x256

interpolated matrix size, 12 slices, 5mm slice thickness, 1 acquisition. The acquired signal intensity was converted to contrast agent concentration using the method of Wang et al. Magn. Reson. Med. 5(5), (1987) pp399-416.

Parameters for the input function #2 were obtained using data taken from the femoral artery. A sequence of the values of the pixel with the highest peak concentration was used, and this sequence is shown with the fitted input function (c_p(t)) in the left panel of figure 2. As with the input functions for the simulation, the input function parameters (including an onset-time parameter) were derived by least-squares fitting the integral of the input function to the cumulative sum of the data.

In figure 2, the left panel shows the data (dots) taken from the femoral artery, and the fitted input function curve (solid line) . For reference the population input function is also shown (dashed line) . The right panel shows the cumulative sum of the data (dots) and the curve used to generate the fits, which is the integral of the input function #2.

The fitted curve and cumulative sum are shown in the right panel of figure 2, and the parameters given in table 3.

Table 3. Estimated Input Function Parameters

Fitting to the cumulative sum of the data tends to be more stable since the noise is smoothed by this process, and also because the parameters for the first function component c_B(t) (bolus model) are forced to match the corresponding area under the data, rather than the data points themselves. This is particularly beneficial in this example as there are two parameters in the model that define the bolus, and only two data points acquired during the bolus passage.

The parameters that best fit the original data will therefore be able to match the two first-pass data points, but the resulting input function is then unrealistic. In this example, fitting the input function directly to the data in the left panel of figure 2 gives an input function with a bolus area 90% larger than the fit shown in the figure, and a peak bolus amplitude of 52.9 mM, which is physiologically implausible. This fit is not shown on the figure as the rest of the curves would not be visible. For reference, the figure also shows the population input function, where it is clear that for this patient, the overall input function amplitude is about twice the population average.

Figure 3 shows estimated representative parameter maps for four regions of interest (ROI) each indicated by an arrow and superimposed on an anatomical Tl-weighted image. Together, the four ROIs contained 3456 pixels, and the combined execution time was 41 seconds using a Pentium 4, 3.4GHz processor with IGB RAM running WindowsXP Professional. The centrally placed ROI is a primary bladder carcinoma, the ROI in the top right of the image is a malignant or involved lymph-node, and the two gluteal muscles are included to give some estimates that can be compared to literature values. The proportion of pixels from which estimates were successfully obtained are 463/474 = 98% for the tumour ROI, 237/253 = 94% for the lymph-node ROI and 2637/2729 = 97% for the muscle ROIs. An unsuccessful fit was defined as one where v_e>l or where the least-squares minimisation routine (e.g. a Levenburg-Marquardt algorithm implemented in IDL) failed to converge.

In figure 3, tissue kinetic parameter maps are shown for a data set containing a primary bladder tumour and a malignant or involved lymph-node. The two ROIs near the bottom of the images are the gluteal muscles, and the tone-scale bars at the right-hand edge of each figure panel show the numerical tone scaling for each parameter. Table 4 gives summary of statistics for the representative parameter estimates for the successful fits in the four ROIs. Of particular note is the plasma fraction in the two tumour regions with median values of 4% for the bladder tumour and 2% for the lymph-node. The simulations indicate that (contingent on the input function being appropriate) these values are reliable with VIF #2 (and VIF #3), but that VIF #1 would give biassed estimates. The plasma fraction is very small in the muscle, and this is as expected for resting muscle.

Table 4. Statistics for in-vivo parameter estimates from three regions using VIF #2 and the parameters in table 3.

The median v_e for the two muscle ROIs is 0.084, which

agrees well with expected values and the median K is

-1

0.027 min which is within the expected range for muscle. The spread of v_e estimates in both tumours is quite large, which is to be expected due to the chaotic nature of tumour anatomy.

Figure 4 schematically illustrates apparatus for estimating a measure of concentration of contrast agent within the body of an imaged subject 11. An image acquisition device (such as an NMR scanner) 10 acquires image data of a region of the body of the subject 11 as a time sequence of successive images of the region within which concentration levels of contrast agent vary. The image data 12 is input to a computer 13 for analysis.

The computer 13 includes a storage device 14 arranged to store image data generated by the imaging apparatus 10. The image data received by the computer may be received directly from the image generator or may be received as pre-stored image data 16 input to the computer via an external data storage device 20 in which previously acquired image data 22 was stored for subsequent input to the computer.

The computer includes a processing unit 15 arranged or programmed to implement the method described above. The processing unit is operably connected to the data store 14 via a data transfer link 18 via which image data is passed to the processor unit. Parameter values determined by the processor unit, according to the methodology described above, are transmittable from the processor unit to the memory store via a data transfer link 17. Stored parameter values calculated by the processor unit 15 may then be passed to a user input/output device and/or graphical interface unit 21 in order to be communicated to the user (e.g. graphic biological data 19) . In addition to, or as an alternative to, the input of image data 16 via a remote storage device 20, representative parameter values associated with a vascular input function may be input to the computer as input data for use, according to a methodology described above, in generating estimates of kinetic parameters associated with contrast agent concentration levels in tissue of the imaged subject. Accordingly, the processor unit 15 may be arranged or programmed to implement such methods as described above, and as illustrated in Figures 1 to 3 herein.

The processor unit is programmed, or arranged, to provide a desired VIF and/or a function c_t(t) representative of a concentration of contract agent in imaged tissue, as desired, with adjustable parameters. The processor unit is arranged to fit the VIF or c_t(t) to input image data as described above to generate representative parameter values and/or estimates of kinetic parameters associated with the imaged subject.

The examples described above are intended to be non- limiting and modification or variants such as would be readily apparent to the skilled person are encompassed by the invention.

Claims

CLAIMS :

1. A method for estimating a measure of the concentration of an image contrast-enhancing agent within an imaged subject including; providing an analytical mathematical function comprising function components and including a sum of a first function component and a second function component in which the second function component includes a convolution of the first function component with a third function component, wherein each function component includes one or more adjustable parameters; adjusting the numerical value of adjustable ^' parameters of the function to fit the function to a sequence of values of an image pixel representing a location in the imaged subject at each of a succession of times thereby to determine values for said parameters representative of the sequence; using the representative parameters to estimate a measure of the concentration of the image contrast-enhancing agent within the imaged subject.

2. A method according to any preceding claim including providing the integral of the function, determining values for a cumulative sum of the sequence of pixel values, adjusting the numerical value of adjustable parameters of the function to fit the integral of the function to the cumulative sum, and determining the representative parameters to be those which optimise the fit of the integral of the function to the cumulative sum.

3. A method according to Claim 2 including generating a cost function and adjusting the numerical value of adjustable parameters φ of the function to minimise the numerical value of the cost function, where the cost function includes the term χ(φ) given by

where

and C₁(J₁,φ) is the numerical value of the integral of the function for the time t_t within the sequence, and y is the j^th pixel value in the sequence.

4. A method according to Claim 2 or 3 in which the integral of the function is an analytical mathematical expression.

5. A method according to any preceding claim in which the image pixel values represent a location containing blood plasma and the method includes using the representative parameters to estimate a measure of the concentration of contrast-enhancing agent within blood plasma of the imaged subject at said location.

6. A method according to any preceding claim in which at least some of the image pixel values represent a, location within a bolus of contrast-enhancing agent and the first function component containing said representative parameters represents a concentration of contrast enhancement agent within the bolus.

7. A method according to any preceding claim in which the first function component includes the term c_B(f) given by c_B{t) = a_Bt^a exp(-μ_Bt) where t represents time, a_B and μ_B are adjustable parameters, and a has integer value or is zero.

8. A method according to any of preceding claims 1 to 6 in which the first function component includes the term c_B(t) given by c_B(t) = a_B(l-cos(_ju_Bt)) where t represents time, a_B and μ_B are adjustable parameters.

9. A method according to any preceding claim in which the third function component includes the term c_G(t) given by c_G(t) = a_c exp(-μ_Gt) where t represents time, a_G and μ_G are adjustable parameters.

10. Apparatus for estimating a measure of the concentration of an image contrast-enhancing agent within an imaged subject including; computer means for providing a representation of an analytical mathematical function comprising function components and including a sum of a first function component and a second function component in which the second function component includes a convolution of the first function component with a third function component, wherein each function component includes one or more adjustable parameters; and for receiving a sequence of values of an image pixel representing a location in the imaged subject at each of a succession of times; the computer means being arranged to adjust the numerical value of adjustable parameters of the function to fit the function to the sequence thereby to determine values for said parameters representative of the sequence; the computer means being further arranged to use the representative parameters to estimate a measure of the concentration of the image-contrast enhancing agent within the imaged subject.

11. Apparatus according to Claim 10 in which the computer means is arranged to provide a representation of the integral of the function, to determine values for a cumulative sum of the sequence of pixel values, to adjust the numerical value of adjustable parameters of the function to fit the integral of the function to the cumulative sum, and to determine the representative parameters to be those which optimise the fit of the integral of the function to the cumulative sum.

12. Apparatus according to Claim 11 in which the computer means is arranged to generate a cost function and to adjust the numerical value of adjustable parameters φ of the function to minimise the numerical value of the cost function, where the cost function includes the term χ{φ) given by

where

and C_j(t_t,φ) is the numerical value of the integral of the function for the time I₁ within the sequence, and y_} is the j^th pixel value in the sequence.

13. Apparatus according to Claim 11 or 12 in which the integral of the function is an analytical mathematical expression.

14. Apparatus according to any of claims 10 to 13 in which the image pixel values represent a location containing blood plasma and the computer means is arranged to use the representative parameters to estimate a measure of the concentration of contrast- enhancing agent within blood plasma of the imaged subject at said location.

15. Apparatus according to any of claims 10 to 14 in which at least some of the image pixel values represent a location within a bolus of contrast- enhancing agent and the first function component containing said representative parameters represents a concentration of contrast enhancement agent within the bolus.

16. Apparatus according to any of claims 10 to 15 in which the first function component includes the term c_B(t) given by c_B(t) = a_Bt^a exp(-μ_Bt) where t represents time, a_B and μ_B are adjustable parameters, and a has integer value or is zero.

17. Apparatus according to any of claims 10 to 16 in which the first function component includes the term c_B(t) given by c_B{t) = a_B(l - cos(μ_Bt)) where t represents time, a_B and μ_B are adjustable parameters .

18. Apparatus according to any of claims 10 to 17 in which the third function component includes the term c_G(t) given by c_G(t) = a_G exp(-μ_Gt) where t represents time, a_G and μ_G are adjustable parameters .

19. Apparatus according to any one of preceding claims 10 to 20 including a computer means programmed to perform the method according to any one of claims 1 to 9.

20. A computer means programmed to perform the method of any one of claims 1 to 9.

21. A computer program product containing a computer program for performing the method of any one of claims 1 to 9.

22. A computer program for performing the method of any one of claims 1 to 9.

23. A method substantially as described in any one embodiment hereinbefore with reference to the accompanying drawings .

24. Apparatus substantially as described in any one embodiment hereinbefore with reference to the accompanying drawings .