Validation and Robustness Analysis of Dynamic Contrast Enhanced MRI

(1)

UmeåUniversity

Radiation Sciences

Master’s Thesis in Engineering Physics

Validation and Robustness Analysis of Dynamic Contrast Enhanced MRI

Author:

Samuel Fransson

safr0031@student.umu.se

Supervisors:

Patrik Brynolfsson Tufve Nyholm

Examiner:

Jonna Wilén

August 17, 2015

(2)

CONTENTS CONTENTS

In Dynamic Contrast Enhanced MRI there are several steps from the initial signal to obtaining the pharmacokinetic parameters for tumor characterization. The aim of this work was to validate the steps in the flow of data focusing on T1-mapping, Contrast Agent (CA)-quantification and the pharmacokinetical (PK) model, using a digital phantom of a head. In the Digital Phantom tissues are assigned necessary values to obtain both a regular and contrast enhanced (using Parker AIF) representation and simulating an SPGR signal.

The data analysis was performed in a software called MICE, as well as the Digital Phantom developed at the department of Radiation Sciences at UmeåUniversity. The method of variable flip angles for the T1-mapping was analyzed with respect to SNR and number of flip angles, finding that the median value in each tissue is correct and stable. A "two point"

inversion recovery sequence was tested with optimal combination of inversion times for white matter and CSF and obtaining correct T1-values when the inversion times were close to the tissue T1, otherwise with large deviations seen. Three different methods for CA-quantification were analyzed and a large underestimation was found assuming a linearity between signal and CA-concentration mainly for vessels at about 60%, but also for other tissue such as white matter at about 15%, improving when the assumption was removed. Still there was a noticeable underestimation of 30% and 10% and the quantification was improved further, achieving a near perfect agreement with the reference concentration, taking the T₂^∗-effect into account. Applying Kety-model, discarding the vp-term, K^transwas found to be stable with respect to noise in the tumor rim but ve noticeably underestimated with about 50%.

The effect of different bolus arrival time, shifting the AIF required in the PK-model with respect to the CA-concentration, was tested with values up to 5 s, obtaining up to about 5%

difference in K^transas well as the effect of a vascular transport function obtained by the means of an effective mean transit time up to 5 s and up to about 5% difference in K^trans.

(5)

1 INTRODUCTION

1 Introduction

1.1 Dynamic Contrast Enhanced Magnetic Resonance Imaging

Magnetic Resonance Imaging (MRI) is an imaging modality using magnetic fields and gradients to reveal the presence of water in the body. When comparing to modalities such as CT it is superior in soft tissue contrast and since using no ionising radiation also associated with low risk.

Most MR-images are of a qualitative character for diagnostic purposes, imaging the volume of interest in the best possible way with respect to resolution, noise, etc. and viewed by a radiologist.

More recently quantitative images have come to be of interest producing so called parametric maps which are applicable when e.g. imaging volume with tumors to give additional insight.

Separating the tissues and determining the severity of the tumor may prove to be a hard task only using the diagnostic images, but since tissue may have very different response and properties one or several quantitative maps may be of help. There are different approaches to obtain the quantitative maps, revealing different properties, but here we are concerned with the Dynamic Contrast Enhanced (DCE) Magnetic Resonance Imaging (MRI), revealing information on uptake and washout in the imaged tissues of a substance called contrast agent (CA) injected into the patient. Where present the relaxation times which affects the signal are shortened. The Blood Brain Barrier, preventing molecules diffusing from blood to the brain, is broken in presence of tumors, hence the uptake of the CA is larger in tumor tissue and visible in the image. Information regarding vascular integrity and plasma fraction of e.g. tumors is obtained by looking at certain parameters such as the transfer constant between blood plasma and intercellular space. These biological properties determines the uptake and washout, favorable when classifying the severity of the tumor [1] Since the size of the tumor is not necessarily affected by treatment but rather the properties treatment response may be determined as well. It has also been shown to be a potential tool for individualizing radiotherapy [2]. Today the commonly used contrast agent is gadolinium based, e.g. Gd-DTPA, a strongly paramagnetic and stable complex. Even low concentrations will affect the signal by decreasing the relaxation time and since well tolerated suited for use in patients. The contrast agent mainly consist of the paramagnetic heavy metal and a chelate, where the latter one has as main purpose to protect the patient against the toxic heavy metal, but also has an effect on the relaxation [2].

When acquiring parametric maps the baseline image without CA is first obtained. The CA is injected and the same volume is scanned for a period of time obtaining a dynamic series from where one can see how the signal changes over time. A T₁-map is obtained by either the variable flip angle method or inversion recovery method and along with the dynamic series the CA-concentration is acquired. The flow of CA in the tissue is modelled with a pharmacokinetical model, applied to the CA-concentration. The most widely applied model is the one called Kety, describing flow of CA with two or three parameters of which the most commonly reported is K^trans, the transfer constant between the blood plasma and the extravascular extracellular space Fig. 1 represents the flow of data. The data analysis from image to the parameter maps are performed by using software with implemented equations and algorithms on pixel by pixel basis

But how accurate are all of the steps in the flow? How accurate is the T1-mapping, and does it depend on the SNR and number of flip angle or inversion recovery series? How accurate is the CA-quantification and what algorithm should one use? Do we obtain correct parameters in the pharmacokinetical (PK) model applied to the CA-quantification? Is it dependent on SNR or temporal resolution? All of these question are hard if not impossible to answer using real data,

(6)

1.2 Purpose 1 INTRODUCTION

Figure 1: The image flow from unprocessed data to parameter maps

since there is no way of knowing the actual T1-map, CA-quantity or parameters.

1.2 Purpose

By using simulated images obtained from already defined parameters all of the input to the signal is known and possible to compare with the calculated. A digital phantom was developed prior to this work as a digital representation of a head based on real MR-images from Brainweb [3, 4].

The possibility to add CA-concentration was implemented and simulated DCE-MRI images can hence be obtained . The purpose of this work is to validate the flow in Fig. 1 by using simulated images, focusing on the T1-mapping, CA-quantification and the parameter mapping, implemented in a software called MICE developed at the department of Radiation Sciences at UmeåUniversity.

If time is available development of the digital phantom implementing image artifacts is to be focused on.

1.3 Limitations

The analysis of the images is carried out in MICE and the validation hence limited to that software. Any clinical data is not taken into the result and the validation is hence limited to simulated data.

(7)

2 THEORY

2 Theory

2.1 MRI Physics

The MRI signal has its basis in the nuclear magnetic moment, arising from proton spin and the intrinsic angular momentum. Generally the spins are randomly distributed but when applying a static magnetic field the spins align with the field according to the Bloch equation, introducing a magnetization M(t)

dM(t)

dt = γM(t) × B(t) (1)

where γ called the gyromagnetic constant is equal to 2.7 · 10⁸ rads⁻¹T⁻¹. The spins aligned will also rotate with an angular frequency around the z-axis according to

ω₀= γB₀ (2)

where B(t) = B0ˆz and ω0 is called the Larmor frequency. If additionally an RF-field with the Larmor frequency is applied the spins are aligned with this field, making it possible to flip the spins into the xy’-plane, see Fig. 2 . When in the xy’-plane the spins will decay back to the equilibrium state, characterized by two time constants T1 and T2

The physical reason for the decay is different for T1 and T2. When the spins are in the xy’-plane they experience slightly different fields due to neighboring spins and molecular motion rotating slightly differently, see Eq. (2), going out of phase. This phenomenon is taken into account by T₂. Translation, vibration and rotation close to the Larmor frequency will cause loss of energy to the surroundings. This is characterized by T₁ and affected by the static field strength due to the frequency dependence.

The sequences to obtain signal from the magnetization are divided into two separate fami- lies called spin-echo and gradient-echo. In spin-echo sequences a 90^◦ -pulse is applied using an RF-field taking the magnetization to the xy’-plane. The spins are dephased due to applied gradients and T2-relaxation and after a time T E/2 a 180^◦ -pulse around the x’-axis shifts the spins location. If assuming that the surroundings does not change considerably the spins still experience the same magnetic field and continue the phase shift in the same direction. Applying a gradient rephasing signal produces an echo after the time TE which is measured. The time between consecutive 90^◦-pulses is denoted TR.

Gradient echo sequences also begin with a pulse flipping the magnetization to the xy’-plane, not necessarily 90^◦but rather the angle θ followed by a gradient field rapidly dephasing the signal.

After the time T E/2 another gradient field in the opposite direction is applied, rephasing the signal, creating an echo. However in this case relaxation due to T₂^∗is seen, shorter than T₂ since this gradient only rephases signal being dephased with the preceding gradient and not dephasing due to molecular motion, field inhomogeneitites etc. When applying short TR magnetization from the previous pulse remains and needs to be removed by applying a spoiler gradient dephasing the remaining magnetization.

For a spoiled gradient echo sequence, SPGR, the signal equation is S = M₀sin (θ) (1 − e^−T^R^/T¹)

1 − cos(θ)e^−T^R^/T¹ e

−TE/T₂^∗ (3)

(8)

2.2 The Contrast Agent 2 THEORY

Figure 2: The rotating frame of reference seen at frequency Ω with static field of magnitude B0 applied along the z’-axis. Applying an RF-field at the same frequency the spins will precess around the effective magnetic field caused by both the static field and the RF-field and if Ω equals the Larmor frequency the effective magnetic field is the RF-field, making it possible to flip the spins to the xy’-plane at an angle θ.

M₀ is the initial magnetization, θ the flip angle, TR repetition time and TE echo time. T₁ and T₂^∗ are tissue characteristic relaxation times.

2.2 The Contrast Agent

An exogenous paramagnetic contrast agent is sometimes used to change the contrast in different tissues by changing the relaxation times according to:

1 Ti

= 1 Ti,0

+ r_iCA_conc (4)

for i = 1,2 where Ti,0 is the initial relaxation time, ri relaxivity and CAconc concentration of the contrast agent. The same equation is valid for T₂^∗, however for blood an additional term is required, q₂^∗, such that

1

T_2,blood^∗ = 1

T_2,0,blood^∗ + r^∗₂CAconc+ q^∗₂CA²_conc (5) The contrast agent is normally injected intravenously, and is distributed from the blood stream to different organs and tissues in the body. The uptake is different between different tissues and sometimes between healthy and diseased tissues, which is useful when e.g. identifying malignancies.

(9)

2.3 SNR 2 THEORY

2.3 SNR

SNR, Signal to Noise Ratio, is defined as the mean pixel intensity in a region divided by the noise standard deviation

SN R = I

σ (6)

Although the mean of the noise is 0 in a normal distribution taking the magnitude of a complex signal generates a rician distribution, requiring a correction factor to Eq. (6) of 0.66, the Rayleigh factor. The practical SNR is dependent on the voxel size and acquisition time in the form

SN R ∝ ∆V

√N_tot

√

BW (7)

where ∆V is the voxel size, Ntotthe total number of acquisitions and BW the bandwidth related to the acquisition time through T_acq= N_tot/BW . Hence a clear correlation with the resolution and SNR where higher SNR requires lower resolution and/or longer imaging time.

2.4 Pharmacokinetical Modelling

2.4.1 Basic physiology

The blood brain barrier, BBB, is a diffusion barrier, preventing molecules from diffusing from the blood into the brain. Contrast agent molecules are too large to penetrate the BBB, so there is generally very little change in signal in healthy tissue in the brain. However, the vasculature around tumors are often coarse and broken, leading to a local disruptions in the BBB. which leads to local signal enhancements around the tumor. There are two important parameters in the study of perfusion. The blood flow carrying the CA-molecules and the permeability, the parameter related to the diffusion of molecules through the blood vessel. The permeability is generally assumed to be constant and the larger area the more diffusion can take place, hence the permeability surface area product, PS, is important.

2.4.2 Haematocrit

The haematocrit, HCT , in blood is defined as the volume ratio of red blood cells to the total volume, generally around 45%. This is an important parameter in the DCE-MRI study because of the CA distribution in the plasma, i.e. the fraction 1-HCT of the blood volume. The CA plasma concentration is however not directly measurable but rather the concentration in the whole arterial blood, ca,b and hence the HCT is required to convert it into the plasma concentration ca(t)

c_a(t) = ca,b(t)

(1 − HCTa) (8)

where HCTa is the haematocrit in arterial blood, and similarly for cp, tissue plasma concentration

cp= ca,b(t)

(1 − HCTt) (9)

(10)

2.4 Pharmacokinetical Modelling 2 THEORY

where HCTt is the haematocrit in tissue blood. [5] The haematocrit explains the necessity of vp

e.g. for the blood

vp= 1 − HCT (10)

2.4.3 Arterial Input Function

Arterial Input Function, AIF , describes the CA input to the tissue of interest [6] and is the concentration of CA in the blood pool [7] required for the pharmacokinetical model. By finding a large artery the AIF can be individually determined using an extraction protocol from e.g Parker et.al [7]. However there are simpler ways, assuming e.g. a functional form of the AIF useful when there is no artery large enough to obtain a patient specific AIF . This model however assumes no change between individuals, although it is probably dependent on e.g. heart rate and kidney function. There is a third way of obtaining the AIF proposed by Parker et.al [8]

in a functional form but with sufficient information to allow for accurate determination of pk- parameters. Performing AIF -measurements in a number of patient and fitting to the data an average AIF is obtained, applicable whenever the patient specific cannot be obtained.

Even if the AIF is obtained correctly at one location, it is not certain that this AIF is applicable at all locations. The dispersion of the contrast agent can be described by

AIF = AIF_d⊗ V T F (t) (11)

where ⊗ denotes convolution, AIF_d(t) is the distal AIF and VTF(t) the vascular transport function, taking into account the transport of the CA. Since this transport function is unknown it is preferable to measure the AIF close to the point for its application [2]. However if one assumes the volume of interest a well mixed compartment one can show the relation [9]

V T F (t) = βe^−βt (12)

where 1/β = M T Tef f, the effective mean transit time from AIFdto AIF . This dispersion model is tested in a DSC-MRI context [9] with values of M T Tef f up to 3 s

2.4.4 Kety model

Pharmacokinetics is the study of how a drug is distributed in the body , and in the MRI context the study of the distribution of the contrast medium. To model this distribution there are several different pharmacokinetical models and in this section the most common compartment (well-mixed space with uniform distribution of CA) models are mentioned called Kety and Brix. Each voxel in the image is considered to consist of a number of compartments depending on the model and a number of parameters are assigned to the model explaining the CA distribution.

The Kety model, a two compartment model visualized in Fig. 3(b) ,consist of a plasma compartment and an EES (Extravascular Extracellular Space) compartment with CA-concentration Cp(t) and Ce(t) and their respective fractional volume vp and ve, always in the range [0, 1]. The plasma compartment represent blood-plasma since the CA is distributed in the body through the blood, and the EES is the space surrounding the cells. K^trans and k21are the parameters describing the flow of CA between compartments.

(11)

Figure 3: A representation of the transport of CA-molecules from the blood vessels to the EES, charac- terized by the pharmacokinetic parameters.

In the Kety model the C_p(t)is assumed equal to the AIF , c_a(t), and hence known. The equation for the system is [10]

ve

dC_e(t)

dt =K^trans(Cp(t) − Ce(t)) and the total voxel concentration

C(t) = vpCp(t) + veCe(t) (13) Solving this yields the total concentration [10]

C(t) = v_pC_p(t) + K^transC_p(t) ⊗ e^−tk^ep (14) where kep= K^trans/ve. This is referred to as the extended Kety model. Since low vp is usually the case this is ignored

C(t) = K^transCp(t) ⊗ e^−tk^ep (15) containing only K^transand veas the unknown parameters and referred to as the Kety model.

2.4.5 Brix model

The Brix model is another two compartment model visualized in Fig. 3(a), differing since it does not require any prior knowledge of the AIF since no assumption of its equality to the plasma concentration is made. The CA-flow is characterized by k12and k21, flow rate between compartments, and Kinand kel, flow in and out from the plasma compartment and the equations become

dM₁(t)

dt =K_in− M1(t) (k_el+ k₁₂) + k₂₁M₂(t) (16) dM2(t)

dt =k12M1(t) − k21M2(t)

(12)

where M₁(t) and M₂(t) are the concentrations in the respective compartment. If relating k₂₁and k12through [11]

vpk12= vek21 (17)

and assuming that veis small, i.e that the peripheral EES compartment has only negligible effect on the central compartment, then

dC_p(t) dt = v_p

Kin

− kelC_p(t) (18)

dCe(t)

dt =k₁₂vp

v_eC_p(t) − k₂₁C_e(t) where v_pC_p(t) = M₁(t) and v_eC_e(t) = M₂(t) Solving yields [12]

C(t) = −AH k21− kel

(e^−k²¹^{(t−T A)}− e^−k^el^{(t−T A)}) (19)

where AH is the amplitude scaling constant, affected by the ratio of the EES, properties of tissue, scan parameters and infusion rate and TA the time of arrival of the CA. This is referred from now as the early Brix model. The Brix model assumes a linear relationship between CA-concentration and relative signal

Sr=S − S0

S₀ (20)

in the form

C(t) = δS_r (21)

where S is the contrast enhanced signal, i.e. from the time where the signal is affected by the CA, S0is the baseline i.e. the signal before contrast injection and δ the proportionality constant.

Hence the possibility to obtain a relation directly between CA and the signal in the form

Sr= −A k21− kel

e^−k²¹^{(t−T A)}− e^−k^el^{(t−T A)}

(22) where the new constant A contains both AH and δ.

If assuming that the flow between the compartments equals permeability surface area product PS and the in and outflow in the vascular compartment equals blood flow Fp, giving [13]

vp

dCp(t)

dt =Fp(ca(t) − Cp(t)) − P S(Cp(t) − Ce(t)) (23) v_edC_e(t)

dt =P S(C_p(t) − C_e(t))

These equations are nonetheless trickier to solve but the impulse response function h(t) is

(13)

derived

h(t) = Fp

e^tK⁺+K+− T_E⁻¹

K₊− K₋(e^−tK⁻− e^−tK⁺)

(24) K_± =1

2(T_p⁻¹+ T_E⁻¹± q

(Tp⁻¹+ T_E⁻¹)²− 4T_E⁻¹T_B⁻¹ Tp= vp

P S + Fp

TE = ve

P S T_B = v_p

Fp

The impulse response functions is interpreted as the concentration at an impulse input to the voxel, i.e. a deltafunction, so the total concentration is

C(t) = h(t) ⊗ ca(t) (25)

where c_a(t) is the AIF . This is referred from now as the open Brix model.

2.4.6 Properties

Having all the solutions the question regarding choice of model is evident. All models have pros and cons so let us list some properties of the models in Table 1

Table 1: Properties of models

Name Need AIF Need T1 #Parameters Plasma Concentration

Kety Yes Yes 2-3 Normal, biexponential

Early Brix No No 3 Fitted with single rate constant

Open Brix Yes No 4 Fitted with single rate constant

Neither of the Brix model require T1-mapping and Early Brix no AIF either, making it easy to implement but K^trans, containing information on permeability surface area product and flow cannot be obtained. However the similar k21 is obtained reflecting the diffusion between plasma and EES. [14]. The simplicity is useful when there is a need for a quick estimation but assuming a linear relationship between CA-concentration and relative signal is only approximately correct at low CA-concentration. The constant A contains information both from the sequence and the tissue, however the product Ak₂₁is related to the ratio permeability surface area product and the unit volume of vasculature [14]. The approximate value of A is [15]

A ≈ r₁TkCp(0)k^{P Sρ}/τ (26)

where r1 is the relaxivity, Tk is defined by Tk = 1

S

∂S

∂(1/T1)= ∂E

∂(1/T1) (27)

(14)

2.5 Contrast Agent quantification 2 THEORY

where

E = S(t)

S(0)− 1 (28)

and τ the duration of the CA-injection.

The more advanced open Brix model requires longer calculation time and possibly a higher noise sensitivity due to more parameters, although the information obtained is more extensive. The separation of PS and flow is possibly most important. However the flow obtained is not the actual flow in the venous but rather an apparent flow due to the necessity of an approximation since the venous output function is unknown. [2]. Looking at Eq. (23) the first part of the right hand side should be replaced by

F (ca(t) − Cp(t)) → F (ca(t) − Cv(t)) (29) where Cv is the venous contrast agent concentration. This parameter however cannot easily be obtained by measurements and therefore one assumes a relation between ca, Cv and Cp

through [13]

ca(t) − Cv(t) = 1

r(ca(t) − Cp(t)) (30)

with constant fraction 0 ≤ r ≤ 1. Hence the relation F =˜ F

r (31)

where the apparent flow obtained is larger than the actual flow. [13]

The Kety model is stable in the form of few parameters explicitly yielding information about the fractional volumes. The signal is not assumed to be proportional to the CA-concentration but both AIF , T1-mapping and r1 is required. As mentioned previously PS is an important parameter reflecting solidity of the BBB, and represented by k21 and K^trans implicitly in the early Brix and Toft model and PS in the open Brix model. The relation between PS, k21and K^transis [15]

k21=k^{P S}^v ve

(32)

where k^{P S}^v is the volume transfer constant related to PS through k^{P S}^v = P Sρ where ρ is the tissue density and [2]

K^trans= F_P

1 − e^{−P S/F}^P

(33) where FP is the flow. It has been shown that PS and ve are important parameters for the tumor and changes in these reflect the extent the chosen treatment has had on the tumor

2.5 Contrast Agent quantification

2.5.1 Linear "low CA"

The contrast agent affects both T1 and T2 and the signal through these parameters. One can

(15)

2.6 T₁-mapping 2 THEORY

and no T₂^∗-effect from Eq. (3) through [2]

C(t) = (r1T10)⁻¹ S(t) S(0)− 1

(34)

where S(t) is the signal at time t, S(0) the baseline signal, r1 relaxivity and T10 the initial T1

before CA-injection.

2.5.2 Linear

A more complete solution, not assuming a low concentration, is to solve directly from Eq. (3) again excluding the T₂^∗-term, [16]

C(t) = 1 T1

− 1 T10

/r1 (35)

E10=e^{−T R/T}¹⁰ (36)

B = 1 − E10

(1 − cos α ∗ E10) A =B ∗ S(t)

S(0) 1

T1(t)= − 1 T R

ln(1 − A) (1 − cos α ∗ A)

2.5.3 Non-linear

Another approach is by creating a least squares optimization procedure from Eq. (3) to the measured signal obtaining the CA-concentration, preferably using the linear "low CA" relation as an initial value, considering the T₂^∗-effect as well.

2.6 T

₁

-mapping

2.6.1 Variable Flip Angle

Through a set of images obtained at different flip angles the T₁-map is acquired through the algorithm [16]

x =S(θ)/ tan θ (37)

y =S(θ)/sinθ T10= − T R/ ln s

S(θ) is the vector containing the signal at all flip angles, the value s is the slope from a linear fit x to y. The algorithm is obtained from the signal Eq. (3) discarding the T₂^∗-term.

(16)

2.7 Temporal resolution 3 METHOD

2.6.2 Inversion Recovery

Another way is to use the method of inversion times and the relation S(T I) = S_e

1 − (1 − k) · e^{−T I/T}¹+ e^{−T R/T}¹

(38)

where TI is inversion time, Sethe signal intensity without the inversion sequence and k = cos(αef f) where αef f is the effective flip angle of the inversion pulse. Applying these sequences usually TR»T1 and the T1-map is acquired using a set of series at varying inversion times.

These sequences are sensitive to the effective flip angle αef f but instead measuring the baseline signal, i.e without any inversion applied, and two signals with different inversion times [17] an estimated T_1mis acquired in this "two point" inversion recovery

T_1m= T I2− T I1

ln_s

e−(−sIR1) s_e−s_IR2

(39)

where T I₁ and T I₂ are the inversion times, s_IR1 and s_IR2 their corresponding signals and s_e the baseline signal, insensitive to α_{ef f}. When calculating the T₁-map using this method the standard deviation of the error is [17]

SET 1=σT 1m

T1

1 2

 σ_s se

T1

T I2− T I1

pexp(−2 · T I₂/T₁) + exp(−2 · T I₁/T₁)

exp(−(T I1+ T I2/T1)) (40) where σT 1mis the standard deviation of the T1-map and σsthe standard deviation of the baseline signal.

2.7 Temporal resolution

The temporal resolution of the acquired dynamic series has to be sufficient to avoid any discrep- ancies between the calculated and actual parameters. Transforming Eq. (14), the Kety-model, to the frequency domain [18]

H(ω) = Ct(t)

C_p(t) = vp+ K^trans

iω + K^trans/v_e (41)

in which ω is the angular frequency one sees the singularity at iω = −K^trans/ve. The frequency range when acquiring a dynamic series is [2π/Tend, 2π/∆T ]where ∆T is the temporal resolution and Tend the acquisition time, where the highest frequency has to be at least two times the singularity frequency, called the cut-off frequency, according to the Nyquist criterion.

3 Method

3.1 MICE-studio

(17)

3.2 The Digital Phantom 3 METHOD

Figure 4: The data flow in MICE. Image data are processed and the result is displayed.

The images are obtain either from the MIQA (Medical Information Quality Archive, a database for medical images developed at the department of Radiation Sciences at UmeåUniversity) database or the local computer. In Fig. 4 an example of a flow is illustrated for DCE-MRI. There are many properties included in the software but here we are only concerned with those related to DCE- MRI, focusing on the T1-quantifier, CA-quantifier and Kety-estimator. The T1-quantifier use the method of Variable Flip Angle or Inversion Times, calculating a T1-map. The CA-quantifier takes as input the T1-map along with a dynamic series and a baseline, calculating the CA-concentration.

Both baseline and dynamic series are obtained from the initial dynamic series data, where the user separates the data into the two parts at the point where the signal enhancement due to CA is visible. The images before this time is taken as the baseline and the images after as the contrast enhanced signal. The Kety-estimator calculates K^trans, veand vp based on Parker AIF and the calculated CA-concentration.

3.2 The Digital Phantom

The digital phantom is essentially a MATLAB-based representation of a head with possibility to add tumors and a contrast agent. The phantom is based on images obtained from BrainWeb [3, 4]

with a 3D-model of a head with a label map describing the different tissues. Each tissue is assigned values for pd (proton density), T1, T2. T₂^∗, X(susceptibility), flow, K^trans, ve, vp, r1, r2, r₂^∗, q^∗₂ and CA (initial CA-concentration). See Appendix for the values applied to each tissue. There are two types of tumors available, called solid and granulated, where the solid is contained in one coherent piece whereas the granulated is more spread out. The size as well as the location must be specified and added to the phantom containing also the tissue values. The CA-concentration in the phantom is obtained with Parker AIF, Extended Kety model and the initial concentration with the tissue specific parameters defined to generate the image.

(18)

3.3 Workflow 3 METHOD

3.3 Workflow

3.3.1 T1-mapping

The validation of the T₁-mapping was performed using the variable flip angle method and the two-point inversion recovery method. For the variable flip angle method a number of images at different flip angles, without any CA-concentration, were obtained from the predefined values in Table A.1 for 1.5 T at different noise levels in the Digital Phantom. The images were then exported in the DICOM format and a T1-map obtained in MICE and the map exported back to MATLAB where it was analyzed. A reference T1-map was obtained from the predefined values generating the signal and the two maps were compared to see sensitivity with respect to noise and number of flip angles. Similarly for the two-point inversion recovery method were instead the inversion times were varied and the T1-map obtained in MICE.

3.3.2 CA-quantification

A dynamic series was simulated applying the Parker AIF and the predefined values for the Extended Kety-model in Table A.1 and exported to MICE. A reference T1-map was exported along with the dynamic series and the CA-concentration was obtained both in MICE and MATLAB using the different methods in section 2.5. The analysis was carried out in MATLAB where the timecurves between calculated and reference CA-concentration were compared. Noise was not applied in the context of testing the CA-quantification only and a reference T1-map applied to isolate the effect of the CA-quantification method applied.

3.3.3 Parameter mapping

The CA-concentration was obtained similarly to above but both with and without noise applied to the dynamic series. The T1-map was obtained using the variable flip angle method both with and without noise applied to the signal and different number of flip angles. In MICE the Kety-model and the Extended Kety-model was applied to the CA-concentration obtaining parameter maps of K^trans, veand vp which are compared to reference maps from the Digital Phantom used to generate the dynamic series.

The validation workflow was setup using the digital phantom and MICE in Fig. 5. Input parameter such as K^trans, veetc. are chosen and image series generated from the digital phantom applying the signal equation, exporting to MICE-studio performing calculations. The input parameters and those obtained from MICE are then compared. In some occasions the images obtained in MICE were exported to MATLAB and analyzed.

(19)

3.3 Workflow 3 METHOD

Figure 5: A validation workflow where images generated with known input are analyzed and the output from the analysis compared with the input.

(20)

4 RESULT

4 Result

4.1 T

₁

-mapping

The procedure described in the theory for T₁-mapping was tested using a range of SNR values and different flip angles with TR = 5.6 ms, TE = 2.63 ms and a resolution of 128x128x60. The map was separated into each tissue and analyzed, yielding Fig. 6. The number of voxels in each tissue at given resolution is seen in Table B.1 in appendix

When applying several flip angles in Fig. 6(a) the calculated T1 seems quite stable except with only minor deviations with higher noise present at 2-3%. Applying only 2 flip angles seen in Fig. 6(b) one sees clearly that the stability is dependent on the T1-value and at tissue with high T1, e.g. vessels, dura matter and CSF, deviations up to 15% is seen at the highest applied noise level. Since the added noise is equal all over the image, the SNR in tissue at higher T1 having a lower signal will also have a lower SNR and the standard deviation of e.g. vessels, dura matter and CSF having the largest T1will also be largest, as seen in Fig. 6, and the SNR given in fat will have the largest SNR of all tissue.

When looking at the procedure of the "two point" inversion recovery any pair of inversion times can be used but the optimal choice is when SN R_T₁= T_1m/SE_{T 1}is at maximum. According to Eq. (40) T I₁should be as short as possible while T I₂is to be optimized. Two series, one baseline series without the inversion sequence and one with T I₁ 50 ms, applying Eq. (38), both with TR 7000 ms and TE 45 ms, acquired from the Digital Phantom are the fixed points in optimization and by using fminsearch in MATLAB the optimal T I₂ is found for each tissue using 2000 ms as initial point. The result from the optimization is shown in Table 2. represents the optimal T I2 for a range of SNR values in the baseline series. Applying the optimal T I2 760 ms for white matter and the optimal 3700 ms for CSF both with a T I1of 50 ms and a baseline series yields Fig. 7 representing relative error in T1 at different SNR in the baseline series for each tissue.

The same level of the noise was applied to the inversion recovery series as the baseline series. In Fig. 7(a) the optimal T I2 for white matter was applied and the error for the tissue with T1 close to that for white matter is relatively small (max 20% in muscle with the highest noise) compared to that for tissue with T1 significantly higher. In dura matter and CSF there are large deviations up to about 60% independent of noise in the image. Fig. 7(b) shows the result when an optimal inversion time for CSF was applied and the result for the tissue with T1in this region is accurate, e.g. for dura matter and CSF, although overestimation up to 15% is seen in the CSF. For tissue with T₁ significantly lower the result is accurate only when no noise is applied, otherwise errors up to a factor of 2-6 is seen

4.2 Contrast Agent quantification

Testing both the linear "low CA" and the linear model in MICE and MATLAB as well as the non-linear in MATLAB, see Fig. 8, there is a clear discrepancy for the linear "low CA" relation mainly for vessels, underestimating with 60%, but also for other tissue as white matter with 15%. Without the assumption of linearity between the CA-concentration and signal the CA- concentration better coincides with the reference, but still an underestimation of 30% for vessels and 10% for white matter. Taking the T₂^∗-effect into account by applying the non-linear procedure (without the q₂^∗-term in Eq. (5)) further improves the quantification. It should however be noted

(21)

4.2 Contrast Agent quantification 4 RESULT

−0.5

−0.25 0 0.25

0.5 a)

Variable Flip Angle Variable flip angle

−0.5

−0.25 0 0.25 0.5

fat(200

ms) around

fat(500

ms) bone

marro w(500

ms) muscle/skin(569

ms) white

matter(656

ms) muscle(1078

ms) gra

ymatter(11

88 ms) vessels(1540

ms) dura

matter(2569

ms) csf(4070

ms)

b) RelativeerrorinT1

FA 2,5,10,15,20,25,30

SN R_{f at}= Inf SN R_{f at}= 61 SN R_{f at}= 31 SN R_{f at}= 20 SN R_{f at}= 15

FA 5,15

Figure 6: Relative error in T1 for each tissue at different SNR for (a) 2,5,10,15,20,25 and 30 degrees and (b) 5,15 degrees

(22)

4.2 Contrast Agent quantification 4 RESULT

−2 0 2 4 6 8 a)

Two point inversion recovery

−2 0 2 4 6 8

fat(200

ms) around

fat(500

ms) bone

marro w(500

ms) muscle/skin(569

ms) white

matter(656

ms) muscle(1078

ms) gra

ymatter(1

188 ms) vessels(1540

ms) dura

matter(2569

ms) csf(4070

ms)

b) RelativeerrorinT1

T I₂ 760 ms ^{SN R}_{SN R}^{f at}^{= Inf}

f at= 61 SN R_{f at}= 31 SN R_{f at}= 20 SN R_{f at}= 15 SN R_{f at}= 12

T I2 3700 ms

Figure 7: Relative error in T1 in each tissue at different SNR for the baseline series at (a) T I2 760 ms optimized for white matter and (b) T I2 3700 ms optimized for CSF

(23)

4.3 Parameter determination 4 RESULT

Table 2: Optimal T I2(ms) for maximizing SNRT₁ in each tissue at different SNR

56 28 19 14 11

Air 1000 1000 1000 1000 1000

CSF 3700 3700 3700 2100 2200

Gray Matter 1200 1200 1200 1200 1200

White Matter 760 760 760 770 770

Fat 250 250 250 250 250

Muscle 1300 1300 1300 1300 1300

Muscle/Skin 670 680 670 680 680

Skull 180 180 190 190 190

Vessels 1800 1800 1800 1800 1900

Around fat 600 600 590 600 600

Dura Matter 3000 3000 3000 3100 3100

Bone Marrow 590 590 600 600 600

time T E = 2.63 ms, flip angle 30^◦ , temporal resolution 5 s and resolution 128x128x60 and if not stated otherwise this setting is applied to all dynamic series

The explanation for the 60% underestimation of the CA-concentration for vessels compared to other tissue for the linear "low CA"-relation is found by looking at the signal relation to the CA-concentration using Eq. (3) in Fig. 9. For white and gray matter the normal CA- concentration peaks around 0.2 mM where a linear relation might be considered. However in vessels the concentration is often peaking at more than 5 mM and no linear relation in this range is seen.

4.3 Parameter determination

The pharmacokinetical model is last in the flow chain and errors introduced earlier in the process hence affects these. Generating images with a range of SNR-values and varying number of flip angles the stability and effects of noise were tested by obtaining a K^trans-map for each SNR.

Analyzed in MATLAB, masking the image giving only values in the tumor rim for about 5000 voxels the median and standard deviation for the volume could be obtained as seen in Fig. 10 representing the relative error in K^transusing the Kety model. The SNR is given as the maximum SNR in the dynamic series

4.3.1 Parameter Phantom

To thoroughly test what ranges the software and models handle a parameter phantom was developed, assigning each voxel a unique combination of the parameters K^trans, ve and vp. obtaining a CA-concentration from each of these combinations using Parker AIF . The Kety- model and Extended Kety-model was applied to these CA-concentrations, calculating K^transand taking the ratio calculated K^transand phantom K^transyielding Fig. 11 and Fig. 12 where the accuracy of the Kety-model alone was tested since applying the reference CA. A signal was also simulated from the CA-concentration using Eq. (3) from where CA-concentration was obtained in MICE with the linear relation, the Kety-model applied acquiring equal pattern as seen in Fig. 11 and Fig. 12.

(24)

0 0.1 0.2 0.3

a)

0 1 2 3 4 5 6 7

0 20 40 60 80 100 120 140 160

a)

b)

CA White Matter

Linear "low CA", Matlab Linear "low CA", MICE Non-linear, Matlab reference Linear, MICE Linear, Matlab

CA-concentration[mM]

Time [s]

CA vessels

Linear "low CA", Matlab Linear "low CA", MICE Non-linear, Matlab reference Linear, MICE Linear, Matlab

Figure 8: The CA-quantification in (a) white matter and (b) vessels using the linear "‘low CA" relation, the linear relation and the non-linear relation between signal and CA-concentration with the reference used to generate the image.

0 0.1 0.2

0 2 4 6 8 10

Magnitudesignal

CA-concentration [mM]

White Matter Gray Matter Vessels

Figure 9: A signal dependence on CA-concentration for white matter, gray matter and vessels using

(25)

−1

−0.5 0 0.5

1 a)

−1

−0.5 0 0.5 1

Inf

0 20 30 40 50 60 70 80 90

a)

b)

Relativeerror

K^trans FA 2,5,10,15,20,25,30 FA 2,30 FA 2,5

SNR

v_e FA 2,5,10,15,20,25,30 FA 2,30 FA 2,5

Figure 10: Relative error in (a) K^trans and (b) veat different SNR and T1-map generated from different set of flip angle applying the Kety model for the tumor rim

(26)

4.4 Partial Volume Effects 4 RESULT

a)

a) b)

c)

a) b)

c) d)

ve,phantom

vp,phantom = 0.01

0.2 0.4 0.6 0.8 1

vp,phantom= 0.087

ve,phantom

K_phantom^trans [min⁻¹] v_p,phantom = 0.174

0.2 0.4 0.6 0.8 1

K_phantom^trans [min⁻¹] v_p,phantom = 0.25

0.2 0.4 0.6 0.8 1

0.1 1 10

Figure 11: The ratio K_calc^trans/K_phantom^trans using the Kety model for different combinations of K_phantom^trans , ve,phantom and vp,phantom generating the CA-concentrations used for the Kety-estimation

The noticeably larger deviations following the increasing v_p up to a factor of more than 10 at some combinations at v_p0.25 in Fig. 11 is expected since the Kety-model Eq. (15) doesn’t take v_p into account, hence this pattern is not seen when applying the Extended Kety Model in Fig. 12, but even at low values there is a clear underestimation of K^trans and a factor of either 0 or close to 0 when at low ve. The temporal resolution 1 s in the dynamic series corresponds to 1 Hz as highest frequency and at the lowest veof 0.01 the Nyquist frequency spans 0.02-0.53 Hz for all K^trans according to Eq. (41), hence lower than 1 Hz. At temporal resolution 0.1 s the same pattern is seen, hence the temporal resolution should be sufficient and is no explanation to the result at low vein Fig. 11 and Fig. 12.

4.4 Partial Volume Effects

To this point the voxels only contain one tissue each. In reality voxels may obtain several tissues depending on the volume it conceals, introducing partial volume effects into the image.

Brainweb [3, 4] supplies weighted images for each tissue where the voxel values ranges between 0 and 1, corresponding to the volume fraction of the specific tissue in the voxel. The phantom was extended to handle these images and introduce partial volume effects by generating signal for each tissue, adding these together, and similarly for tumor rim with the possibility to weigh its occurrence in the voxel. Different images with a weighted rim were created and their corresponding

(27)

4.4 Partial Volume Effects 4 RESULT

a)

a) b)

c)

a) b)

c) d)

ve,phantom

vp,phantom = 0.01

0.2 0.4 0.6 0.8 1

vp,phantom= 0.087

ve,phantom

K_phantom^trans [min⁻¹] vp,phantom = 0.174

0.2 0.4 0.6 0.8 1

K_phantom^trans [min⁻¹] vp,phantom = 0.25

0.2 0.4 0.6 0.8 1

0.1 1 10

Figure 12: The ratio K_calc^trans/Kphantom^trans using the extended Kety model for different combinations of Kphantom^trans , ve,phantom and vp,phantomgenerating the CA-concentration used for the Kety-estimation

(28)

4.5 Bolus Arrival Time 4 RESULT

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

0 20 40 60 80 100 120 140 160

a)

0 0.5 1 1.5 2

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 a)

b)

CA-concentration[mM]

Time [s]

w = 1 w = 0.7 w = 0.5 w = 0.3 w = 0.2 w = 0.1 w = 0

weight

K_rel^trans ve,rel

Figure 13: (a)The CA-concentrations obtained for different weights of the tumor rim and (b), the relative K^transand ve with respect to the reference

since the reference CA-concentration is applied. Fig. 13 representing the CA-concentration for different rim weights and relative K^trans and ve. The overestimation of calculated K^trans at weight 1 of about 25% is expected since applying Kety-model discarding the vp term in the calculation, although it is present generating the image. Fig. 13 (a) clearly shows decreasing CA-concentration at lower weight and Fig. 13 (b) an approximate linearly decreasing trend for K^transbetween 0.2 and 1 in weight. A clear linearly decreasing trend for veis seen as well as an underestimation with 50% at weight 1.

4.5 Bolus Arrival Time

Since the calculation of parameters requires an input AIF a shift due to slightly different Bolus Arrival Time, BAT, might induce an error in the estimation of parameters. Even if the user supplies the correct division of timeframe of the initial data it is possible that the actual timepoint for this is somewhere in between the chosen and the neighboring timeframe, hence a shift of

±timeresolution is possible. Testing this shift with a time resolution of 5 s , Fig. 14 (a), shows stability for ve with respect to BAT but a noticeable shift up to 8% is seen for K^trans in the negative direction.

(29)

4.5 Bolus Arrival Time 4 RESULT

0.9 1 1.1

−5 −4 −3 −2 −1 0 1 2 3 4 5

a)

0.9 1 1.1

0 1 2 3 4 5

a)

b)

BAT [s]

K^trans_rel ve,rel

M T Tef f[s]

K^trans_rel v_e,rel

Figure 14: Relative K^trans and ve at different BAT (a) and MTTef f (b) to the reference corresponding to perfect alignment/no dispersion between the CA-concentration and the AIF applied to the Kety-model.

(30)

4.6 AIF -dispersion 4 RESULT

4.6 AIF -dispersion

As mentioned the dispersion of the AIF caused by the vascular transport function is difficult to model, but applying Eq. (12) it is possible to simulate dispersion by means of the M T Tef f. Simulations have been done with several values for M T Tef f creating a range of different CA- concentration curves using Parker AIF . The Kety model was applied with a nondispersed AIF trying to fit to the different CA-curves. Applying for the tumor rim yields Fig. 14 (b) where v_e only deviates at most about 3% at an M T T_{ef f} of 5 s but K^transup to about 8%.

Validation and Robustness Analysis of Dynamic Contrast Enhanced MRI

UmeåUniversity

Radiation Sciences

Master’s Thesis in Engineering Physics

Validation and Robustness Analysis of Dynamic Contrast Enhanced MRI

Author:

Samuel Fransson

safr0031@student.umu.se

Supervisors:

Patrik Brynolfsson Tufve Nyholm

Examiner:

Jonna Wilén

August 17, 2015

Contents

1 Introduction

1.1 Dynamic Contrast Enhanced Magnetic Resonance Imaging

1.2 Purpose

1.3 Limitations

2 Theory

2.1 MRI Physics

2.2 The Contrast Agent

2.3 SNR

2.4 Pharmacokinetical Modelling

2.5 Contrast Agent quantification

2.6 T

-mapping

2.7 Temporal resolution

3 Method

3.1 MICE-studio

3.2 The Digital Phantom

3.3 Workflow

4 Result

4.1 T

-mapping

4.2 Contrast Agent quantification

4.3 Parameter determination

4.4 Partial Volume Effects

4.5 Bolus Arrival Time

4.6 AIF -dispersion