WO2007051224A1

WO2007051224A1 - A method of characterising an area of interest within a body

Info

Publication number: WO2007051224A1
Application number: PCT/AU2006/001533
Authority: WO
Inventors: Paul R. Clark; Wanida Chua-Anusorn; Wenjie Pang; Sander J. Bangma
Original assignee: Resonance Health Analysis Services Pty Ltd
Priority date: 2005-11-01
Filing date: 2006-10-17
Publication date: 2007-05-10
Also published as: AU2006308776A1

Abstract

The present invention provides a method of characterising an area of interest within a body. The method comprises obtaining a frequency distribution of intensities of a plurality of magnetic resonance responses received from respective regions in the area of interest and characterising the obtained frequency distribution to obtain numerical parameters. The method also comprises incorporating at least one of the numerical parameters in a marker for a disease to determine the presence or degree of the disease.

Description

A METHOD OF CHARACTERISING AN AREA OF INTEREST WITHIN A

BODY

Field of the Invention

The present invention broadly relates to a method of characterising an area of interest within a body.

Background of the Invention Chronic liver diseases are amongst the most common diseases throughout the world, and encompass a diverse range of clinical manifestations including liver fibrosis, cirrhosis and cancer. Fibrosis can be classified as a wound-healing response to a variety of chronic stimuli. Progressive liver fibrosis or excessive deposition of extracellular matrix is the main cause of organ failure of any aetiology, from a variety of causes including viral, autoimmune, drug induced, cholestatic and metabolic disease and iron overload. The degree of scarring (fibrosis) of liver tissue is the single most important predictor of complications in patients with chronic liver diseases. Currently the gold standard for assessing fibrosis is examination of a specimen of liver obtained via an invasive biopsy procedure. The results of the assessment assist in a clinician's determination of patient treatment. The biopsy carries risks of complications ranging from pain through to internal bleeding or even death. In addition, poor patient acceptance makes use of the procedure problematic. Hence, there is a need to develop accurate and reliable non-invasive means to assess the severity of liver scarring. In assessment of fibrosis, liver biopsies are clearly- inadequate and have various disadvantages. The accuracy of the histological assessment of fibrosis in liver biopsy sample depends on the size of the specimens. Furthermore, interpretation of the biopsies carries intra-observer and inter-observer variation of 10-20% even among experienced pathologists.

Blood tests are one alternative where a combination of bio-chemical markers may be measured to provide an indirect assessment of fibrogenesis in the liver. None of the available biomarkers are ideal as they are not liver specific. It is also not possible to determine the degree of fibrosis.

Radiological techniques for the non-invasive measurement of liver fibrosis are presently in the early stages of clinical development. Ultrasound based transient elastography may be used to evaluate fibrosis by noninvasive measurement of tissue elasticity. However, ultrasound elastography is limited by sampling size, as it is measures only 0.5 to 1% of the total liver volume. Further, elasticity measurement may be perturbed by a number of factors such as fatty liver, ascites, blood flow, iron deposition, and narrow rib cages.

There is a clear need for an alternative non-invasive method of characterising an area of interest within a body, such as the human liver.

Summary of the Invention

The present invention provides in a first aspect a method of characterising an area of interest within a body, the method comprising: obtaining a frequency distribution of intensities of a plurality of magnetic resonance responses received from respective regions in the area of interest, characterising the obtained frequency distribution to obtain numerical parameters and incorporating at least one of the numerical parameters in a marker for a disease to determine the presence or degree of the disease.

The area of interest typically is an area of a liver and the disease typically is a liver disease, such as fibrosis, cirrhosis or cancer. In this case the method has the particular advantage that a non-invasive characterisation of the liver may be performed which typically has a higher sensitivity to fibrosis, cirrhosis or cancer than conventional characterisation methods such as those involving blood tests.

It is to be appreciated, however, that in variations of this embodiment the area of interest may not be associated with the liver but may be associated with any other body part.

The method may comprise the additional step of obtaining the marker for the disease, such as a marker for liver fibrosis, cirrhosis or cancer. The marker may be obtained by comparing numerical parameters associated with magnetic resonance responses for a group of individuals having the disea'se with those associated with individuals being disease free. For example, the marker may be obtained from a combination of numerical parameters determined by statistically analysing magnetic resonance characterisation results of the liver together with at least one other liver characterisation method. The statistical analysis may be a combination of multivariate logistic regression and receiver operating characteristic (ROC) curve analysis. The magnetic resonance responses may be responses detected by a magnetic resonance measuring instrument or may be obtained by processing the responses detected by the magnetic resonance measuring instrument. For example, the magnetic resonance response may be responses from iron or from iron related compounds in the liver. The magnetic resonance responses typically are associated with a relaxation property of respective regions, such as a transverse relaxation property (R₂) . For example, the transverse relaxation property may be determined by any one of: a mono-exponential analysis of the decay in transverse magnetisation; a bi-exponential analysis of the decay in transverse magnetisation; a multi- exponential analysis of the decay in transverse magnetisation; or a stretched-exponential analysis of the decay in transverse magnetisation.

The step of obtaining the frequency distribution typically comprises generating a magnetic resonance image of the area of interest . In this case the magnetic resonance image typically is generated for R₂ and each region may be associated with a voxel in the magnetic resonance image. For example, the intensity of each voxel may be a measure for the R₂ value recorded ^"for that voxel .

The step of characterising the obtained frequency distribution typically is computer software implemented and typically comprises displaying the obtained frequency distribution of the intensities. For example, the obtained frequency distribution of the intensities may be displayed in form of a graph or plot . The step of characterising the obtained frequency distribution to obtain the numerical parameters may comprise fitting a model to the distribution. For example, the model may comprise several component distributions, such as Gaussian distributions.

If the transverse relaxation response is determined via .. bi-exponential analysis of the decay in transverse magnetisation, the transverse relaxation response typically is any one of: a fast rate of decay in transverse magnetisation (R_2f) ; the population density of the fast rate of decay in transverse magnetisation (Pf) ; a slow rate of decay in transverse magnetisation (R_2s) ; the population density of the slow rate of decay in transverse magnetisation (P_s) ; or any combination of R₂f, Pf, R_2S, and/or P₃. To obtain the numerical parameters, the frequency distribution typically is characterised by statistical measures such as those associated with mean, standard deviation, coefficient of variation, skewness, kurtosis, Kolmogorov-Smirnov statistic, or Kuiper's statistic. Alternatively, the or each fitted component frequency distribution may have parameters derived from a scaled probability density function (PDF) model fitted to the obtained frequency distribution. For such a PDF the distribution may be a Gaussian or normal distribution, a Poisson distribution, a Weibull distribution, or any other appropriate statistical distribution. ' In one specific embodiment of the present invention the or each fitted component frequency distribution is associated with a Gaussian distribution. For example, a number of Gaussian frequency distribution components and their initial numerical parameter values may be determined by: generating a smoothed approximation to the obtained frequency distribution; subtracting the smoothed approximation of the obtained frequency distribution from the frequency distribution to obtain a difference distribution; for positive regions of the difference distribution, where the area of a given positive region is greater than a specified threshold value, designating a Gaussian component to be located in the vicinity of the bounds of that region on the frequency distribution to contribute to the multiple Gaussian model; and. determining initial estimates of the numerical parameter values for the Gaussian components.

For example, the initial estimates of the numerical parameters may include the mean, standard deviation or amplitude of the Gaussian components from measurements based on either the obtained frequency distribution, the smoothed approximation to the obtained frequency distribution, and/or the difference distribution. The method typically also comprises thresholding the obtained frequency distribution to retain only those distribution values that are above a predetermined percentage of the maximum distribution value. The step of characterising the obtained frequency distribution may comprise fitting multiple Gaussian components to the obtained frequency distribution. For example, this may include minimising a difference in a height of the multiple Gaussian distribution relative to a height of the obtained frequency distribution. The fitting of the multiple Gaussian components to the obtained frequency distribution may be performed with parameter constraints to determine optimum values for the mean, standard deviation, and amplitude of each Gaussian component .

The smoothed approximation to the obtained frequency distribution typically is determined as either an approximation to the distribution based on a single Gaussian component, an approximation based on neighbourhood averaged data of the measurements on the area of interest, an approximation based on neighbourhood averaged distribution values, or an approximation based on broader area-equivalent binning elements for the distribution.

The threshold value to determine which positive difference region will foster an associated Gaussian component typically is a percentage cut-off value of the area of the positive difference region relative to either the maximum positive region area, the total area of all positive regions, the total area of the difference distribution, or the area of the actual frequency distribution.

The initial estimate of the mean for the Gaussian component fostered from a respective positive difference region may be calculated as either the position of a maximum of the difference distribution between bounds of that region, the position of the maximum on the frequency distribution between the bounds of that region, or a weighted combination thereof.

The initial estimate of the amplitude of each Gaussian component may be set to a height of the obtained frequency distribution at the estimated mean position for the Gaussian component. The initial estimate of the standard deviation for the Gaussian component that is either a leftmost or rightmost Gaussian component may be calculated according to the expression

where σ is the standard deviation, a is the amplitude of the Gaussian component, μ is the mean of the Gaussian component, p is a percentage between 0% and 100%, and x_pa is the position of the frequency distribution where the height of the distribution is p percent of the amplitude a.

Further, an additional Gaussian component may be fitted that is neither the leftmost nor the rightmost Gaussian component. The initial estimate of the standard deviation for that Gaussian component may be calculated according to the expression

where σ is the standard deviation, μ is the mean of the Gaussian component, a is the amplitude of the Gaussian component, is the minimum frequency distribution value

between the mean positions of the two neighbouring Gaussian components, and is the position of the minimum frequency

distribution value.

In one specific embodiment of the present invention a multiple Gaussian model of more than two Gaussian components is reduced to an initial estimate of a two Gaussian model by progressively merging Gaussian components between the leftmost or rightmost Gaussian components (hereinafter referred to as inner Gaussian components) with either the leftmost or rightmost Gaussian components in a left to right or right to left fashion. In this case the method typically comprises : calculating a crossing point of an inner Gaussian component with the leftmost Gaussian component which occurs between the means of the two components (hereinafter referred to as the left-hand cross-point) , and determining the height of either Gaussian at the left-hand cross-point; calculating a proximity measure of either the position of the left-hand cross-point relative to the mean of the leftmost Gaussian,_, the height of the left-hand cross-point relative to the amplitude of the leftmost Gaussian_/ or any combination thereof, for which the measure increases in value as proximity of the left-hand cross-point to the peak of the leftmost Gaussian increases; calculating the crossing point of the inner Gaussian component with the rightmost Gaussian component which occurs between the means of the two components (hereinafter referred to as the right-hand cross-point) , and determining the height of either one of the Gaussians at the right-hand cross-point; calculating a proximity measure of either the position of the right-hand cross-point relative to the mean of the rightmost Gaussian, the height of the right-hand cross-point relative to the amplitude of the rightmost Gaussian, or any combination thereof, for which the measure increases in value as proximity of the right-hand cross-point to the peak of the rightmost Gaussian increases; and comparing the proximity measures for the left-hand and right-hand cross-points, and if the left-hand cross-point proximity measure is greater merging the inner Gaussian component and the leftmost Gaussian component into a new single Gaussian component that is reclassified as the leftmost Gaussian component, or if the right-hand cross- point proximity measure is greater merging the inner Gaussian component and the rightmost Gaussian component into a new single Gaussian component that is reclassified as the rightmost Gaussian component.

Merging of a leftmost or rightmost Gaussian component with an inner Gaussian component typically comprises fitting a single Gaussian component to the summation of the leftmost or rightmost Gaussian component with the inner Gaussian component .

The proximity measure of the left-hand or right-hand cross-point position relative to the mean of the leftmost or rightmost Gaussian component respectively may be equal to

respectively, where μ is the mean of the inner Gaussian component, μ_ι is the mean of the leftmost Gaussian component, and μ_r is the mean of the rightmost Gaussian component .

Alternatively, the proximity measure of the left-hand and/or right-hand cross-point height relative to the leftmost or rightmost Gaussian amplitude may be equal to

where x_c is the position of the left-hand or right-hand cross-point, μ is the mean of the leftmost or rightmost Gaussian component, σ is the standard deviation of the leftmost or rightmost Gaussian component, and a is the amplitude of the leftmost or rightmost Gaussian component. The initial estimate of the numerical parameters associated with the two Gaussian model may be refined by fitting the model to the frequency distribution with associated parameter constraints to determine optimum values for the mean, standard deviation, and amplitude of the two Gaussian components.

The step of fitting may also comprise fitting a first and a second Gaussian component and one further component associated with a mixture basis function, for example to account for mixing of diseased and healthy tissue below the length scale of a region within the area of interest. In this case the mixture basis function may be specified by the expression

where y_m(x) is the mixture basis function,

is the amplitude of the mixture basis function, are the

mean and standard deviation respectively of the first Gaussian component, and

and

are the mean and standard deviation respectively of the second Gaussian component.

The component associated with the mixture basis function typically is fitted to the obtained frequency distribution with associated parameter constraints to determine optimum values for the mean, standard deviation, and amplitude of the two Gaussian components, as well as the amplitude of the mixture basis function. For example, the amplitude of the component associated with the mixture basis function may be fixed to the height of a position at which the first and second Gaussian components cross each other between their means.

The present invention provides in a second aspect a computer system in which the above-described method is implemented.

The present invention provides in a third aspect a method of obtaining a marker for liver fibrosis and/or cirrhosis in which a numerical parameter can be incorporated to determine the presence of the liver fibrosis and/or cirrhosis, the numerical parameter being associated with a frequency distribution of intensities of a plurality of magnetic resonance responses received from respective regions in the area of interest, the method comprising: providing numerical characterisation parameters associated with liver characterisation methods, statistically analysing the numerical characterisation parameters to select those numerical characterisation parameters which have a dependency on a degree of the disease and which can be used in a formula for the marker and combining the selected numerical characterisation parameters in the formula for the marker.

The step of statistically analysing the numerical characterisation parameters typically comprises selecting those numerical characterisation parameters which have a correlation with a degree of the disease. The step of statistically analysing typically also comprises multivariate logistic regression and receiver operating characteristic (ROC) curve analysis.

The present invention provides in a fourth aspect a method of characterising a liver within a body, the method comprising: obtaining magnetic resonance responses from the liver, characterising the magnetic resonance responses to obtain selected numerical parameters and incorporating at least one of the selected numerical parameters in a marker for a liver disease to determine the presence or degree of the liver disease.

The invention will be more fully understood from the following description of specific embodiments of the invention. The description is provided with reference to the accompanying drawings .

Brief Description of the Drawings

Figure 1 shows a flow-chart illustrating a method of characterising an area of interest within a body according to an embodiment of the present invention; Figure 2 shows a flow-chart illustrating a method of characterising an area of interest within a body according to another embodiment of the present invention;

Figure 3 shows an example of R₂ frequency distribution characterisation for a patient confirmed to have liver fibrosis by needle biopsy;

Figure 4 shows a diagrammatic representation of the steps in determining a multiple Gaussian model to be fitted to the R₂ frequency distribution, where the number of Gaussian components to be fitted and their initial parameter values are determined by a distribution differencing technique;

Figure 5A shows an example R₂ frequency distribution of a patient with fibrosis for which the number and location of multiple Gaussian components has been determined from the distribution differencing technique and the initial parameter estimate methods;

Figure 5B shows the fit of the multiple Gaussian components of Figure 5A to the R₂ frequency distribution, to enable the determination of optimum values for the mean, standard deviation, and amplitude of each Gaussian component ;

Figure 6 illustrates the generation of a two Gaussxan model from a multiple Gaussian model via a process of merging an inner Gaussian component with either a leftmost or rightmost Gaussian component;

Figure 7A shows an example R₂ frequency distribution of a patient with fibrosis with the initial parameter estimates of a two Gaussian model following the merging of multiple Gaussian components;

Figure 7B shows the fit of the two Gaussian model of Figure 7A to the R₂ frequency distribution, to enable the determination of optimum values for the mean, standard deviation, and amplitude of each Gaussian component;

Figure 8 shows a two Gaussian fit of the R₂ frequency distributions for 3 patients with varying degrees of liver fibrosis; Figure 9 shows one possible MR marker, which is a combination of two numerical parameters characterising the frequency distribution, that distinguishes healthy individuals from cirrhotic individuals;

Figure 10 shows receiver operating characteristic (ROC) plots associated with three selected markers;

Figure 11 shows box plots associated with a first selected marker;

Figure 12 shows box plots associated with a second selected marker; and Figure 13 shows box plots associated with a third selected marker.

Detailed Description of Specific Embodiments

Referring initially to Figure 1 a method characterising an area of interest within a body according to a first aspect of the present invention is now described.

In broad terms embodiments of the present invention provide a method of characterising an area of interest in a body, such as a liver, to identify and/or quantify a disease such as liver fibrosis and/or cirrhosis.

In this embodiment this is achieved by characterising frequency distributions of processed magnetic resonance image intensities to obtain numerical parameters on the frequency distributions that are incorporated in MRI markers of fibrosis and/or cirrhosis.

Figure 1 shows the general flow of the method 100.

Initially one or more magnetic resonance images are acquired of an organ or tissue, which are then processed to give a quantitative image of magnetic resonance relaxation properties for a region of that organ or tissue (step 102) . A frequency distribution of magnetic resonance relaxation properties of that region is then generated (step 104) . The nature of the frequency distribution for that region is then characterised by a variety of numerical parameters (step 106) . The numerical parameters are correlated with measures of fibrosis and/or cirrhosis measured by any other method. In this embodiment, the correlated numerical parameters are formulated into a predictive MRI marker of liver fibrosis and/or cirrhosis by statistical analysis techniques to determine the presence of the disease (step 108) .

The flow chart shown in Figure 2 shows an illustration of a more specific example. Figure 2 illustrates method steps of method 200 for identifying and/or quantifying liver fibrosis and/or cirrhosis by characterising frequency distributions of transverse relaxation properties throughout a region of interest of the liver. Figure 2 outlines the process where bi-exponential transverse relaxation analysis has been considered performed. Bi-exponential transverse relaxation analysis is conducted as transverse relaxation behaviour in liver tissue predominantly arises from two distinct populations of hydrogen protons with two different transverse relaxation rates, one fast (R2f) and one slow (R2s) . The two different population components are postulated to arise from different hydrogen proton pools within the tissue. Liver fibrosis impacts on the hydrogen proton pools and is thus expected to leave a footprint on the frequency distributions that is characteristic of fibrosis and/or cirrhosis.

In the process outlined in Figure 2, a series of spin- echo images of the liver are acquired (step 202) , for which the decay in transverse magnetisation is modelled by a bi- exponential transverse relaxation process. Further, only the overall contribution from the sum of the fast and slow transverse relaxation rates is considered in this example, to give a single R₂ image (step 204) . A frequency distribution for the R₂ image is then generated (step 206) . The R₂ frequency distribution is then characterised by a variety of numerical parameters, which may include (step 208) : ^■ statistical parameters such as the mean, standard deviation, coefficient of variation, skewness, kurtosis, the Kolmogorov-Smirnov statistic for normalcy, and Kuiper' s statistic for normalcy (step 208a);

^■ parameters of scaled probability density functions (PDFs) fitted to the frequency distribution, for such PDFs as the Gaussian distribution, the generalised Poisson distribution, and the Weibull distribution (step 208b) ;

■ parameters of a multiple Gaussian model fitted to the frequency distribution, for which the parameters are the mean, standard deviation, and amplitude of each Gaussian component (step 208c) ;

■ parameters of a two Gaussian model fitted to the frequency distribution, where the initial parameter values for the two Gaussians are determined by a process of merging Gaussians from a multiple Gaussian model (step 208d) ;

■ parameters of a two Gaussian model with an accompanying mixture basis function of variable amplitude that is fitted to the frequency distribution (step 208e) ; or

■ parameters of a two Gaussian model with an accompanying mixture basis function with an amplitude fixed to the crossing point of the two Gaussians that is also fitted to the frequency distribution (step 208e) .

The numerical parameters are then statistically screened to determine those that best correlate with a known measure of fibrosis and/or cirrhosis (step 210) , and a statistical analysis is conducted to formulate an optimally predictive MR marker of liver fibrosis and/or cirrhosis (step 212) .

Examples of R₂ frequency distribution characterisation for the numerical parameters described above (steps 208 and 208a to 208e) are shown in Figure 3 (302, 302a to 302f respectively) for a patient confirmed to have liver fibrosis by needle biopsy.

For the multiple Gaussian model fitted to the R₂ frequency distribution, the number of Gaussian components to be fitted and their initial parameter values - the mean, standard deviation, and amplitude of each component - are determined by a distribution differencing technique. Method steps associated with this technique are illustrated in Figure 4 (inset A) . For the sake of description, the R₂ frequency distribution may be thought of as having been generated from a region of interest of a liver R₂ image. The R₂ frequency distribution 400 is shown with voxel count displayed as a function of R₂. The frequency distribution is then thresholded to retain only those distribution values that are above a given percentage level 402 of the maximum distribution value to remove long tails from the distribution. A "smoothed" frequency distribution 404 is then generated by either smoothing the R₂ image and generating a subsequently smoothed distribution, smoothing the frequency distribution directly by neighbourhood averaging, generating a smoothed distribution via broader area-equivalent binning elements to those used in construction of the frequency distribution, or by approximating the frequency distribution by a single Gaussian distribution. A difference distribution 406 is then generated by subtracting the smoothed frequency distribution from the original frequency distribution. Positive regions 408 of the difference distribution are then identified, and the number of Gaussian components 410 which will contribute to the multiple Gaussian model are determined based on the number of positive regions that have areas greater than a specified threshold value. The threshold value is specified as a percentage cut-off value for the area of the positive difference region relative to the maximum positive region area, the total area of all positive regions, the total area of the difference distribution, or the area of the actual frequency distribution. The Gaussian component fostered by a positive difference region is located in the vicinity of the bounds of that region on the frequency distribution. The initial estimates as to the mean, standard deviation, and amplitude of the component are determined from measurements based on either the frequency distribution, the smoothed approximation to the frequency distribution, the difference distribution, or any combination of all three distributions. In one specific embodiment, an initial estimate of the mean for a Gaussian component fostered from a positive difference region (Figure 4, inset B) may be calculated as either the position of the maximum on the difference distribution 412 between the bounds of the region, the position of the maximum on the frequency distribution 414 between the bounds of the region, or a weighted combination of the two different maximum position values. An initial estimate of the amplitude for a Gaussian component is then set to the height of the frequency distribution 416 at the estimated mean position for the Gaussian component.

To calculate an initial estimate of the standard deviation for a Gaussian component, two approaches are employed depending on whether the Gaussian component is a leftmost or rightmost Gaussian component (Figure 4, inset C) , or one in between (Figure 4, inset D) . An initial estimate of the standard deviation for a Gaussian component that is either the leftmost or rightmost Gaussian component is calculated according to the expression

where σ is the standard deviation, a is the amplitude 418 of the Gaussian component, μ is the mean 420 of the Gaussian component, p is a percentage between 0% and 100%, and x_pa 422 is the position of the frequency distribution where the height of the distribution is p percent of the amplitude a.

An initial estimate of the standard deviation for a Gaussian component that is not the leftmost or rightmost

Gaussian component is calculated according to the expression

where σ is the standard deviation, μ is the mean 424 of the Gaussian component, a is the amplitude 426 of the Gaussian component, is the minimum frequency distribution value

428 between the mean positions of the two neighbouring Gaussian components, and x_v is the position 430 of the minimum frequency distribution value .

An example of the number of Gaussian components obtained from the distribution differencing technique and the initial appearance of the components determined by initial estimates calculated as described above is shown in Figure 5A for the R₂ frequency distribution of a patient with fibrosis. The final fit of the multiple Gaussian model to the R₂ frequency distribution is shown in Figure 5B, for which the fitting was accomplished by minimising the difference in height of the multiple Gaussian model relative to the height of the frequency distribution with appropriate constraints on the Gaussian parameters, to enable the determination of optimum values for the mean, standard deviation, and amplitude of each Gaussian component.

For the two Gaussian model to be fitted to the R₂ frequency distribution (with or without the inclusion of a mixture basis function) , the initial estimates for the parameters of the two Gaussian components are determined from a multiple Gaussian model by progressively merging Gaussian components between the leftmost or rightmost Gaussian components (hereinafter referred to as inner Gaussian components) with the leftmost or rightmost Gaussian components in a left to right or right to left fashion.

The process of merging an inner Gaussian component with either a leftmost or rightmost Gaussian component is illustrated in Figure 6. The crossing point of the inner Gaussian component with the leftmost Gaussian component (hereinafter referred to as the left-hand cross-point 602) , which occurs between the means of the two components is calculated, and the height 604 of either Gaussian at the left-hand cross-point is determined. Similarly, the crossing point of the inner Gaussian component with the rightmost Gaussian component (hereinafter referred to as the right- hand cross-point 606) , which occurs between the means of the two components is calculated, and the height 608 of either Gaussian at the right-hand cross-point is determined. A proximity measure of the left-hand cross-point 602 with the peak of the leftmost Gaussian 610 and a proximity measure of the right-hand cross-point 606 with the peak of the rightmost Gaussian 612 is then determined. The left-hand proximity measure may be calculated as either the position of the left-hand cross-point 614 relative to the mean of the leftmost Gaussian 616, the height of the left-hand cross- point 604 relative to the amplitude of the leftmost Gaussian 618, or any combination thereof, provided that the measure increases in value as proximity of the left-hand cross-point to the peak of the leftmost Gaussian increases. Similarly, the right-hand proximity measure may be calculated as either the position of the right-hand cross-point 620 relative to the mean of the rightmost Gaussian 622, the height of the right-hand cross-point 608 relative to the amplitude of the rightmost Gaussian 624, or any combination thereof, provided that the measure increases in value as proximity of the right-hand cross-point to the peak of the rightmost Gaussian increases. The proximity measures are then compared, and if the left-hand cross-point proximity measure is greater the inner Gaussian component is merged with the leftmost Gaussian component into a new single Gaussian component that is reclassified as the leftmost Gaussian component, or if the right-hand cross-point proximity measure is greater the inner Gaussian component is merged with the rightmost

Gaussian component into a new single Gaussian component that is reclassified as the rightmost Gaussian component. The preferred method of merging involves fitting a single Gaussian component to the summation of the leftmost or rightmost Gaussian component with the inner Gaussian component. Following merging of the first inner Gaussian component, the remaining inner Gaussian components are considered for merging with either the leftmost or rightmost Gaussian components (in a left to right or right to left fashion) , until only the leftmost and rightmost Gaussian components remain.

In this specific embodiment the proximity measure of the left-hand cross-point position relative to the mean of the leftmost Gaussian is given by

where μ is the mean of the inner Gaussian component 626, μ_ι is the mean of the leftmost Gaussian component 616, and μ_r is the mean of the rightmost Gaussian component 622.

Similarly, the proximity measure of the right-hand cross- point position relative to the mean of the rightmost Gaussian is given by

where μ is the mean of the Gaussian component 626, μ_r is the mean of the rightmost Gaussian component 622, and μ_ι is the mean of the leftmost Gaussian component 616. In addition, the proximity measure of the left-hand or right-hand cross- point height relative to the amplitude of the leftmost or rightmost Gaussian is given by

where x_c is the position of the left-hand or right-hand cross-point,

is the mean of the leftmost or rightmost Gaussian component, σ is the standard deviation of the leftmost or rightmost Gaussian component, and a is the amplitude of the leftmost or rightmost Gaussian component. An example of the results of the merging technique to obtain the initial parameter estimates of the two Gaussian model is illustrated in Figure 7A for the R₂ frequency distribution of a patient with fibrosis. The final fit of the two Gaussian model to the R₂ frequency distribution is shown in Figure 7B, for which the fitting was accomplished with appropriate parameter constraints to enable the determination of optimum values for the mean, standard deviation, and amplitude of each Gaussian component. For the two Gaussian model which includes a mixture basis function to account for two material mixing below the length scale of resolution of the R₂ image, the mixture basis function is specified by the expression

where y_m(x) is the mixture basis function, a_m is the amplitude of the mixture basis function,

and

re the mean and standard deviation respectively of the first Gaussian component, and μ₂ and σ₂ are the mean and standard deviation respectively of the second Gaussian component.

The two Gaussian model incorporating the mixture basis function is fitted to the frequency distribution with associated parameter constraints to determine optimum values for the mean, standard deviation, and amplitude of the two Gaussian components, as well as the amplitude of the mixture basis function. The two Gaussian model with the mixture basis function may also be fitted where the amplitude of the mixture basis function is fixed to the height of the two Gaussian components where they cross between their means . Following calculation of the various numerical parameters that characterise the R₂ frequency distribution, the parameters are statistically screened to determine those that best correlate with a known measure of fibrosis and/or cirrhosis. (In the case of a bi-exponential transverse relaxation analysis of the magnetic resonance image data with unknown fast and slow transverse relaxation components, numerical parameters may have also been generated for the R_2ff P_f/ R_2s and P₃ frequency distributions.) Multivariate logistic regression and receiver operating characteristic (ROC) curve analysis is then performed to determine one or more predictive MRI markers of liver fibrosis and/or cirrhosis from a combination of the screened numerical parameters .

An illustration of the potential magnetic resonance frequency distribution characterisation of liver fibrosis is shown in Figure 8, which shows R₂ images and corresponding two Gaussian fits of the R₂ frequency distributions for 3 patients with varying degrees of liver fibrosis as determined by the Metavir scoring system. The Metavir scoring system is a well-established system for grading liver biopsy specimens in a semi-quantitative manner by a pathologist . The Metavir system is based on a 5 point scale: FO, no fibrosis; Fl, portal fibrosis alone; F2 , portal fibrosis with rare septae; F3, portal fibrosis with many septae; F4, cirrhosis. As shown in Figure 8 a-c, for increasing Metavir score of fibrosis, the frequency distribution changes shape from being characterised by predominantly one Gaussian component to being clearly comprised of increasing amounts of a second Gaussian component. Growth of a second Gaussian component on the R₂ frequency distribution may thus be associated with progression of liver fibrosis.

An example of one possible combination of numerical parameters from the characterisation of R₂ frequency distributions to give a potential MR marker to distinguish between healthy individuals and those with cirrhosis is shown in Figure 9. The particular MR marker shown is the ratio of Kuiper's statistic calculated on the R_2s frequency distribution (R2s KP-V) to the kurtosis calculated on the R_2f frequency distribution (R2f kurtosis) . This MR marker distinguishes a healthy group of individuals from two groups of cirrhotic individuals, one with normal hepatic iron concentration (HIC) , and the other with elevated HIC. The following presents the mathematical formulae for characterising the transverse relaxation distributions, which generate the sum total of MR numerical parameters . As described above, a frequency distribution of intensities of magnetic resonance responses is obtained for the liver and the obtained frequency distribution is characterised to obtain the numerical parameters which are then incorporated in the marker for a disease to determine the presence or degree of the disease.

Magnetic Resonance (MR) statistical parameters calculated on transverse relaxation distributions

The basis statistical measures calculated on the di-stributions of transverse relaxation properties as a function of each property x (where x is either R₂, R2f, R2s, P_f7 or P_s) , and for which there are N voxel measurements x_n of a given property for the region of interest over which the distribution is generated may include:

■ The mean μ(x) :

■ The standard deviation σ(x)

■ The coefficient of variation CV(x) :

^■ The skewness Skew(x) :

^■ The kurtosis Kurt(x) :

^■ The Kolmogorov-Smirnov statistic KS_D(x) :

where P(x) is the cumulative probability distribution for x, and P_norm(x) is the normal cumulative probability distribution for x based on and

■ Kuiper's statistic KP_v(x) :

The probability density functions fitted to the distributions of transverse relaxation properties as a function of each property x (where x is either R₂, R_2f, R_2s, P_f, or p_s) may include:

^■ The Gaussian distribution gss(x) :

where

is the amplitude of the Gaussian distribution, is the mean, and

is the standard deviation.

^■ The generalised Poisson distribution gpd(x) :

where is an amplitude scale factor, defines

the location, and

defines the shape of the generalised Poisson distribution (GPD) . The mean

of the GPD is:

and the standard deviation of the GPD is:

■ The standard Weibull distribution swb(x)

where

is an amplitude scale factor,

defines the location,

is a scale factor, and is a shape factor. The mean is ^•■

and the standard deviation is:

where is the gamma function.

^■ An effective Weibull distribution ewb(x) :

where

is an amplitude scale factor,

is the effective mean is the effective

standard deviation, and

is a shape factor,

Multiple Gaussian models fitted to the distributions of transverse relaxation properties as a function of each property x (where x is either R₂, R_2f, R_2s, P_{f >} or p_s) may- include :

" A multi-Gaussian model gsn(x) of N_gsn(x) Gaussians matching the number of significant peaks on the distribution:

where

is the amplitude of the «"' Gaussian distribution, is the mean, and is

the standard deviation.

" A 2 -Gaussian model:

where

is the amplitude of the 1" Gaussian distribution, is the mean, and is the

standard deviation, and is the amplitude

of the 2^nd Gaussian distribution,

is the mean, and

is the standard deviation, and where

" A 2-Gaussian model incorporating a mixture basis function to account for partial volume effect:

where is the amplitude of the 1^st Gaussian

distribution, is the mean, and is

the standard deviation, is the amplitude

of the 2^nd Gaussian distribution, is the

mean, and

is the standard deviation, and ^{is the} amplitude of the mixture basis

function, and where

^■ A 2-Gaussian model incorporating a mixture basis function of amplitude fixed to the height of the cross- point of the 2 Gaussians :

where

is the amplitude of the V Gaussian distribution,

is the mean, and is the

standard deviation,

is the amplitude of the 2^nd Gaussian distribution,

is the mean, and

is the standard deviation, and

is the amplitude of the mixture basis function, which may be given by:

where x_cis the cross-point of the 2 Gaussians, and where

The following MR parameter combinations were derived from the parameters above for incorporation into the logistic regression and ROC curve analysis to generate the predictive MR markers of fibrosis:

■ The product of the amplitude of the 1^st Gaussian component from a 2-Gaussian fit to the R_2s distribution with the amplitude of the 2^nd Gaussian component aussian fit to the R_2f distribution :

A logistic weighting of the peak position on the p_s histogram

about the mean peak position for

the entire volunteer group:

^■ The percentage area of the Gaussian component of minimum area from a 2-Gaussian fit to the p_s distribution, weighted by the difference of the means for the two components:

Table la-c lists the optimal MR parameters, from the above defined mathematical formulae, which are of greatest diagnostic value for potentially distinguishing liver fibrosis or cirrhosis. The parameters are ranked according to probability (Wilcoxon test) of distinguishing healthy volunteers (n = 12) from patients with liver disease as diagnosed by liver biopsy (Table Ia, n = 96) , and from subgroups of patients with normal (Table Ib, n = 36) and elevated (Table Ic, n = 60) liver iron concentration (p < 0.05) .

Table Ia

Table Ib

Table Ic

MR markers for discrimination of clinically significant fibrosis The following defines three examples of predictive MR markers for liver fibrosis. Logistic regression and ROC curve analysis of the optimal MR parameters resulted in several models from which three MR markers were selected. A first marker was determined from ROC curve analysis of the entire volunteer group, and included 8 of the optimal MR parameters listed in Table Ia plus age. Second and third predictive markers were determined by ROC curve analysis of the volunteer sub-groupings of individuals with normal LIC and elevated LIC, with the additional inclusion of the parameters listed in Table Ib and Ic, respectively. The predictive formulae for the two subgroups were then combined into an overall predictive model of fibrosis by logistic weighting according to the liver iron concentration determined from the mean R₂ value in the liver. The second predictive marker derived by this method included 10 MR parameters, whilst the third predictive marker included 7 MR parameters plus age . The three predictive markers are listed below, for which the relative weight of each parameter and the intercept are quoted with the standard error.

Marker 1:

where

and

Marker 2 :

where

Marker 3 :

Figures 10 shows ROC curves for the markers 1 to 3 that demonstrate their ability to distinguish clinically significant fibrosis (F2-F4) according to both METAVIR grading and percentage fibrosis which was obtained by computerised morphometric analysis of liver biopsy samples. The markers were modelled based on the ability to distinguish clinically significant fibrosis (F2-F4) according to the METAVIR system and thresholding of the percentage fibrosis measurement. For each marker, ROC curves were generated for the entire volunteer group, as well as the volunteer sub-groups with normal liver iron concentrations and elevated liver iron concentrations. Significantly, less variation in AUC was obtained for the fibrosis models between groups on the percentage fibrosis measurements, which was a phenomenon observed during fibrosis model development. In addition, the differences between all AUCs for the models with respect to percentage fibrosis were within the standard errors. For the ROC curves determined against METAVIR grading, Marker 1 and Marker 3, which included the age of the volunteers, had significantly higher AUCs for the entire volunteer group than did Marker 2, which did not include age. Marker 1 and Marker 3 also had significantly higher AUCs than Marker 2 against METAVIR grading for the sub-group of volunteers with normal liver iron concentrations. Furthermore, Marker 3 had a significantly higher AUC (outside of the standard error) than Marker 1 for the normal liver iron concentration group. For the sub-group of volunteers with elevated liver iron concentrations, there was no significant difference between markers as determined by the AUC against METAVIR grading. In addition, Marker 3 was also the only maker to show a significantly higher AUC for the normal liver iron concentration group when compared with the elevated liver iron concentration group.

Figures 11, 12 and 13 show box plots for the first, second and third marker respectively that distinguish healthy volunteers from volunteers with disease (Fig. 11a, 12a and 13a)_., volunteers with disease grouped according to clinically significant fibrosis (FO-Fl and F2-F4, Fig. lib, 12b and 13b) , and volunteers with disease according to distinct METAVIR grading (Fig. lie, 12c and 13c) . All markers showed a significant difference between healthy volunteers and the volunteers with disease according to the Wilcoxon test (p < 0.0001), with Marker 3 showing the greatest separation between the two groups. The markers also showed a significant difference between healthy volunteers and volunteers with disease of METAVIR grading FO-Fl (p = 0.0287 for Marker 1, p = 0.0003 for Marker 2, and p = 0.0019 for Marker 3), for which Marker 2 without age was the most statistically significant. A significant difference between groups comprising healthy volunteers and distinct METAVIR gradings was also observed (p < 0.0001) .

Although the invention has been described with reference to particular examples, it will be appreciated by those skilled in the art that the invention may be embodied in many other forms. For example, any other type of suitable marker for the disease may be used. Further, the area of interest may not be associated with the liver and the disease may not be associated with a liver disease and the marker may be a marker for a disease other than a liver disease.

Claims

What is claimed is :

1. A method of characterising an area of interest within a body, the method comprising: obtaining a frequency distribution of intensities of a plurality of magnetic resonance responses received from respective regions in the area of interest, characterising the obtained frequency distribution to obtain numerical parameters and incorporating at least one of the numerical parameters in a marker for a disease to determine the presence or degree of the disease.

2. The method of claim 1 wherein the area of interest is an area of a liver and the disease is a liver disease.

3. The method of claim 2 wherein the liver disease is at least one of liver fibrosis, liver cirrhosis and liver cancer.

4. The method of claim 2 or 3 wherein the marker is a marker for the liver disease.

5. The method of any one of the preceding claims comprising the additional step of obtaining the marker for the disease.

6. The method of claim 5 wherein the marker is obtained by comparing numerical parameters associated with magnetic resonance responses for a group of individuals having the disease with those associated with individuals being disease free.

7. The method of claim 6 wherein the marker also allows quantifying severity of the disease.

8. The method of any one of the preceding claims wherein the magnetic resonance responses are detected by a magnetic resonance measuring instrument .

9. The method of any one of claims 1 to 7 wherein the magnetic resonance responses are obtained by processing the responses detected by a magnetic resonance measuring instrument .

10. The method of any one of the preceding claims wherein each magnetic resonance response is associated with a relaxation property of the respective region.

11. The method of any one of the preceding claims wherein the magnetic resonance responses are associated with a transverse relaxation property (R₂) .

12. The method of claim 11 wherein the transverse relaxation property is determined by any one of : a mono- exponential analysis of the decay in transverse magnetisation; a bi-exponential analysis of the decay in transverse magnetisation; a multi-exponential analysis of the decay in transverse magnetisation; or a stretched- exponential analysis of the decay in transverse magnetisation.

13. The method of claim 11 or 12 wherein the transverse relaxation response is determined via bi-exponential analysis of the decay in transverse magnetisation and is one Of : a fast rate of decay in transverse magnetisation (R_2f) ; the population density of the fast rate of decay in transverse magnetisation (P_f) ; a slow rate of decay in transverse magnetisation (R_2s) ; the population density of the slow rate of decay in transverse magnetisation (P₃) ; or any combination of

14. The method of any one of the preceding claims wherein the step of analysing a plurality of magnetic resonance responses comprises generating a magnetic resonance image of the area of interest .

15. The method of claim 14 wherein the magnetic resonance image is generated for R₂ and each region is associated with a voxel in the magnetic resonance image.

16. The method of claim 15 wherein the intensity of each voxel is a measure for the R₂ value recorded for that voxel .

17. The method of any one of the preceding claims wherein the step of characterising the obtained frequency distribution is computer software implemented.

18. The method of any one of the preceding claims wherein the step of characterising the obtained frequency distribution comprises calculating statistical measures associated with at least one of a statistical mean, a standard deviation, a coefficient of variation, skewness, kurtosis, Kolmogorov-Smirnov statistic, or Kuiper's statistic.

19. The method of any one of claims 1 to 17 wherein the step of characterising the obtained frequency distribution comprises fitting a model having one or more component distributions .

20. The method of claim 19 wherein each fitted component frequency distribution has parameters derived from a scaled probability density function (PDF) model fitted to the obtained frequency distribution.

21. The method of claim 20 wherein the PDF is one of: a Gaussian distribution, a Poisson distribution, or a Weibull distribution.

22. The method of any one of claims 1 to 20 wherein the step of characterising the obtained frequency distribution comprises fitting at least one Gaussian distribution component and the initial numerical parameter values are determined by: generating a smoothed approximation to the obtained frequency distribution; subtracting the smoothed approximation of the obtained frequency distribution from the frequency distribution to obtain a difference distribution; for positive regions of the difference distribution, where the area of a given positive region is greater than a specified threshold value, designating a Gaussian component to be located in the vicinity of the bounds of that region on the frequency distribution to contribute to the multiple Gaussian model; and determining initial estimates of the numerical parameter values for the Gaussian components.

23. The method of claim 22 comprising thresholding the obtained frequency distribution to retain only those distribution values that are above a predetermined percentage of the maximum distribution value.

24. The method of claim 22 wherein fitting multiple Gaussian components comprises minimising a difference in a height of the multiple Gaussian component peaks relative to a height of the obtained frequency distribution.

25. The method of claim 22 wherein the smoothed approximation to the obtained frequency distribution is determined as one of: an approximation to the distribution based on a single Gaussian component, an approximation based on neighbourhood averaged data of the measurements on the area of interest, an approximation based on neighbourhood averaged distribution values, or an approximation based on broader area-equivalent binning elements for the distribution.

26. The method of claim 22 wherein the threshold value to determine which positive difference region will foster an associated Gaussian component is a percentage cut-off value of the area of the positive difference region relative to either the maximum positive region area, the total area of all positive regions, the total area of the difference distribution, or the area of the actual frequency distribution.

27. The method of claim 22 wherein the initial estimate of the mean for the Gaussian component fostered from a respective positive difference region is calculated as one of: a position of a maximum on the difference distribution between bounds of that region, the position of the maximum on the frequency distribution between the bounds of that region, or a weighted combination thereof.

28. The method of claim 22 wherein the initial estimate of the amplitude of each Gaussian component is set to a height of the obtained frequency distribution at the estimated mean position for the Gaussian component.

29. The method of claim 22 wherein the initial estimate of the standard deviation for the leftmost and rightmost Gaussian component is calculated according to the expression

where σ is the standard deviation, is the amplitude of

the Gaussian component, μ is the mean of the Gaussian component, p is a percentage between 0% and 100%, and x_pa is the position of the frequency distribution where the height of the distribution is p percent of the amplitude a.

30. The method of claim 22 wherein an additional Gaussian component is fitted that is neither the leftmost nor the rightmost Gaussian component and the initial estimate of the standard deviation for that Gaussian component is calculated according to the expression

where σ is the standard deviation, μ is the mean of the Gaussian component, a is the amplitude of the Gaussian component, y_nύn is the minimum frequency distribution value between the mean positions of the two neighbouring Gaussian components, and is the position of the minimum frequency

distribution value.

31. The method of claim 22 wherein a multiple Gaussian model of more than two Gaussian components is reduced to an initial estimate of a two Gaussian model by progressively merging an inner Gaussian component with either the leftmost or rightmost Gaussian components in a left to right or right to left fashion.

32. The method of claim 31 comprising: calculating a crossing point of an inner Gaussian component with the leftmost Gaussian component which occurs between the means of the two components (hereinafter referred to as the left-hand cross-point) , and determining the height of either one of the Gaussians at the left-hand cross-point; calculating a proximity measure of either the position of the left-hand cross-point relative to the mean of the leftmost Gaussian, the height of the left-hand cross-point relative to the amplitude of the leftmost Gaussian, or any combination thereof, for which the measure increases in value as proximity of the left-hand cross-point to the peak of the leftmost Gaussian increases; calculating the crossing point of the inner Gaussian. component with the rightmost Gaussian component which occurs between the means of the two components (hereinafter referred to as the right-hand cross-point) , and determining the height of either Gaussian at the right-hand cross-point; calculating a proximity measure of either the position of the right-hand cross-point relative to the mean of ths rightmost Gaussian, the height of the right-hand cross-point relative to the amplitude of the rightmost Gaussian, or any combination thereof, for which the measure increases in value as proximity of the right-hand cross-point to the peak of the rightmost Gaussian increases; and comparing the proximity measures for the left-hand and right-hand cross-points, and if the left-hand cross-point proximity measure is greater merging the inner Gaussian component and the leftmost Gaussian component into a new single Gaussian component that is reclassified as the leftmost Gaussian component, or if the right-hand cross- point proximity measure is greater merging the inner

Gaussian component and the rightmost Gaussian component into a new single Gaussian component that is reclassified as the rightmost Gaussian component .

33. The method of claim 32 wherein merging of the leftmost or rightmost Gaussian component with an inner Gaussian component comprises fitting a single Gaussian component to the summation of the leftmost or rightmost Gaussian component with the inner Gaussian component.

34. The method of claim 32 wherein the proximity measure of the left-hand cross-point position relative to the mean of the leftmost Gaussian component is equal to

where μ is the mean of the inner Gaussian component, μ_ι is the mean of the leftmost Gaussian component and μ_r is the mean of the rightmost Gaussian component .

35. The method of claim 32 wherein the proximity measure of the right-hand cross-point position relative to the mean of rightmost Gaussian component is equal to

where μ is the mean of the inner Gaussian component, μ_r is the mean of the rightmost Gaussian component and μ, is the mean of the leftmost Gaussian component.

36. The method of claim 32 wherein the proximity measure of the left-hand or right-hand cross-point height relative to the amplitude of either the leftmost or rightmost Gaussian component is equal to

where x_c is the position of the left-hand or right-hand

cross-point, μ is the mean of the leftmost or rightmost Gaussian component, σ is the standard deviation of the leftmost or rightmost Gaussian component, and a is the amplitude of the leftmost or rightmost Gaussian component.

37. The method of any one of claims 1 to 31 wherein the step of characterising the obtained frequency distribution comprises fitting at least two Gaussian component distributions and a component associated with a mixture basis function.

38. The method of claim 37 wherein the mixture basis function may be specified by the expression

where is the mixture basis function, is the

amplitude of the mixture basis function, and are the

mean and standard deviation respectively of the first Gaussian component, and

39. The method of claim 38 wherein the amplitude of the component associated with the mixture basis function is fixed to the height of a position at which the first and second Gaussian components cross each other between their means .

40. The method of claim 5 or 6 wherein the marker is a marker for liver fibrosis and/or cirrhosis and is obtained from a combination of numerical parameters determined by statistically analysing magnetic resonance characterisation results of a liver together with at least one other method of liver characterisation.

41. The method of claim 40 wherein the statistical analysis comprises multivariate logistic regression and receiver operating characteristic (ROC) curve analysis.

42. A computer system in which the method according to any one of the preceding claims is implemented.

43. A method of obtaining a marker for liver fibrosis and/or cirrhosis in which a numerical parameter can be incorporated to determine the presence or degree of the liver fibrosis and/or cirrhosis, the numerical parameter being associated with a frequency distribution of intensities of a plurality of magnetic resonance responses received from respective regions in the area of interest, the method comprising: providing numerical characterisation parameters associated with liver characterisation methods, statistically analysing numerical characterisation parameters to select those numerical characterisation parameters which have a dependency on a degree of the disease and which can be used in a formula for the marker and combining the selected numerical characterisation parameters in the formula for the marker.

44. The method of claim 43 wherein the step of statistically analysing numerical characterisation parameters comprises selecting those numerical characterisation parameters which have a correlation with a degree of the disease.

45. The method as claimed in claim 43 or 44 wherein statistically analysing comprises multivariate logistic regression and receiver operating characteristic (ROC) curve analysis.

46. A method of characterising a liver within a body, the method comprising: obtaining magnetic resonance responses from the liver, characterising the magnetic resonance responses to obtain selected numerical parameters and incorporating at least one of the selected numerical parameters in a marker for a liver disease to determine the presence or degree of the liver disease.