M³ – Monologues, MRI, Miscellanea

View-angle tilting MRI (part 3): Fewer approximations

This is a third (and presumably final) text about view-angle tilting (or VAT) MRI concluding the explanations posted on 2020-09-10 and on 2020-09-15. I don’t really like the mathematical approximations that I used in my first text and, therefore, I would like to try again; the simple graphical visualization should correspond to an equally simple mathematical derivation …

Let’s start as before with the MRI signal equation for an excited slice without any distortions

with frequency-encoding direction x and frequency-encoding gradient G_x (phase encoding is generously omitted). The x integral is over the complete imaged object (but for a compact object that fits into our field of view, we can integrate from  – ∞ to ∞). The slice-selection direction is denoted by ẑ, where I use the hatted (“^”) variable to differentiate the undistorted (original) slice position ẑ from the distorted position z = z(x) analyzed below.

I will ignore the slice thickness, since the essence of the VAT technique can pretty well be explained without it. And everything gets really complicated if we consider varying slice thicknesses – or even slices split in several parts as illustrated in the following figure.

It’s useful to keep in mind the resonance relation between an excitation pulse with frequency ω_rf applied after switching on the slice-selection gradient G_z,slc and the slice position ẑ (without field inhomogeneities):

If we now consider inhomogeneity-induced distortions of a slice (due to a field inhomogeneity, ΔB(x,z)), then each original slice position ẑ (corresponding to the excitation frequency ω_rf) is “moved” to a new position, z(x), which is, in general, different for each x. The new position, z(x), is the solution of the equation

for each (fixed) value of x. The relation to the original slice position, ẑ, can be expressed (after division by γG_z,slc) as

which can be re-arranged to

(Of course, the previous equation is not an explicit solution describing z(x), since z(x) still appears also on the right-hand side of the equation.)

If we go back to the MRI signal equation and include the distorted slice geometry, z(x), the integral changes to

where, ω(x,z) is used to describe the Larmor frequencies during readout. Without field inhomogeneities, this is simply (and independent of z) ω(x,z) = γG_x x. The inhomogeneity, ΔB(x,z) changes the Larmor frequencies to ω(x,z(x)) = γ (G_x x + ΔB(x, z(x))). If we also include the additional VAT readout gradient, G_z,VAT, we get

We can now insert the expression for z(x) derived above into the third addend of this last expression, which gives

So, by setting G_z,VAT = G_z,slc (as proposed for VAT MRI), the inner parenthesis vanishes and neither ΔB nor z(x) (which implicitly contains ΔB as well) appears in the expression for the Larmor frequency. This means that the artifacts in readout direction (due to changed readout frequencies) are removed:

Other artifacts caused by the changed slice geometry, z(x), obviously remain (as indicated by the term ρ(x, z(x)) in the integral) and are not corrected by the VAT approach.

Since I’ve omitted the integration (∫dz …) over the finite slice thickness, the meaning of the slice coordinate ẑ in the final results is not completely obvious. Actually, the appearance of ẑ during readout corresponds to a tilted readout (that’s why it’s called view-angle tilting), which becomes clearer if one includes the omitted integration over the slice direction in the resulting MRI signal equation.

View-angle tilting MRI (part 2): Explanations and visualizations

Here I would like to extend the introduction to view-angle tilting (or VAT) MRI posted on 2020-09-10. Using the same field inhomogeneity and slice geometry as in the previous text, I’m trying to provide some more (visual) explanations about what’s going on.

Let’s consider two sample slices as shown in the following figure. Both slices have the same thickness Δz and are really meant to be at the same z position (with respect to the field inhomogeneity that I’m going to add later on); they are only plotted at different z positions to combine them into one diagram.

The top slice is a homogeneous slice with the same signal (magnetization density, color-coded in gray) for all x positions (except for the periphery without signal). The bottom slice has alternating signals with two low-signal areas; there is a change from dark to bright signal exactly at the center of the slice. Below each slice, the signal intensity, I(x) – obtained by a readout in x direction – is plotted.

Using the field inhomogeneity defined in the previous text, both slices can be plotted with inhomogeneity-induced distortions in z direction, which is illustrated in the next figure. If we assume (which is generally not the case in real life) that the signal of the underlying object doesn’t vary in slice (z) direction, then the signal intensity I(x) remains exactly the same as before (i. e., the projection of the signals onto the x axis is unchanged). (This is true at least as long as we ignore that the thickness of the excited slice might also change as function of x.)

However, as previously explained the distorted slice geometry is not the only effect of the field inhomogeneity. Field inhomogeneities during readout influence the spatial encoding in readout (x) direction. This is illustrated in the next figure, in which the slice distortion effect in z direction is omitted for better visualization and only in-slice effects are visualized.

The voxels are scaled in x direction (some are wider, others narrower than in the reference plot above). Since the number of protons in each voxel remains the same, the signal, which is proportional to the proton density, is scaled inversely proportionally to the voxel width, so narrow voxels appear brighter. The strongly increased signal intensity, I(x), in the right half of the plot is called “signal pile-up”.

In addition, the in-slice geometry is distorted which is nicely illustrated in the blue intensity curve: the change from low to high signal is clearly shifted from the exact center to the right-hand side of the plot; this is, of course, a direct consequence of the wider voxels in the left half of the slice.

Both effects, the distorted slice geometry as well as the in-slice signal artifacts are shown together in the following figure. The spatial shift in x direction becomes very obvious as a shift of the minimum of the slice shape to the right. The signal intensity after readout, I(x) remains the same as in the previous figure.

This situation now is the starting point for artifact correction with the view-angle tilting approach. The artifacts shown in the signal intensity curves, I(x), can – perhaps surprisingly – be corrected by a simple change of the view angle as illustrated below.

Just by projecting the distorted spin density onto a tilted axis (shown at the bottom of the following plot), the signal artifacts (areas with low signal or signal pile-up) are corrected. In this example, the projection is numerically evaluated and yields the almost perfectly straight intensity curve in the bottom plot.

A second visual “proof” of the VAT approach are the gray projection lines (one for each voxel) that appear equidistant in the tilted projection, but very non-equidistant in the conventional projection shown at the top of the figure.

In addition, the sinc-smoothing effect of the VAT pulse sequence is visualized at the periphery of the projected slice, where the signal increases (and decreases) in a relatively long ramp instead of a sharp step.

Finally, it’s interesting to see how the strength, G_z, of the slice-selection gradient influences the VAT approach (assuming that the readout gradient, G_x, is not changed). A lower slice-selection gradient results in more severe through-slice distortions (since the field inhomogeneities are now stronger relative to the gradient). But the increased slice distortion has a surprising advantage in terms of in-slice artifacts as illustrated in the following figure: increased through-slice (z) distortions shift the slice away from the maximum of field inhomogeneities and, thus, reduce the remaining field inhomogeneity in the slice (along the x direction) after excitation.

The reduced in-slice artifacts are shown in the next figure together with the uncorrected signal intensity (top) and VAT-corrected intensity (bottom).

Obviously, the voxels are less distorted along the readout (x) direction, and the uncorrected signal demonstrates reduced signal pile-up. The angle of the tilted axis at the bottom is smaller than above; it’s always arctan(G_z / G_x) which gives 45° in the first case and about 26.565° now. The numerical projection on this axis shows strongly reduced artifacts (compared to the uncorrected readout) as well as reduced sinc-smoothing compared to the VAT readout above.

View-angle tilting MRI: Theory and approximations

This text is first of all an explanation and note to myself – but some parts of it might be of more general interest (if this is not the case for you, please just stop reading …).

Metallic implants can cause severe image artifacts in MRI. The problem of these implants in MRI is their magnetic susceptibility, χ, that influences the static magnetic field B₀ in their neighborhood. Consequently, the (ideally) perfectly homogeneous static field B₀(x,y,z) = B₀ e_z becomes slightly (or not so slightly) inhomogeneous: B₀(x,y,z) = B₀ + ΔB(x,y,z) (considering only the z component from here on). A simplified two-dimensional example of such an inhomogeneity ΔB(x,z) is shown in the following figure (together with another kind of desired inhomogeneity, namely a linear gradient field).

MR imaging usually relies on the assumption that B₀ is homogeneous; the presence of ΔB(x,y,z) results in image artifacts such as voxels moved to wrong positions, signal pile-up, and signal voids.

To simplify things, I ignore the phase-encoding direction (y) and restrict this discussion to the frequency-encoding direction (x, with frequency-encoding gradient G_x; displayed horizontally) and the slice direction (z, with slice-selection gradient G_z; displayed vertically). Slice selection in MRI works by applying a slice-selection gradient G_z during excitation. This means that each position z corresponds to a unique Larmor frequency ω(z) = γz G_z, which is used to excite the spins at the desired z positions. A slice is selected by applying a radio-frequency (rf) pulse with a certain frequency (ω_rf) and bandwidth (Δω_rf). Thus, all spins with Larmor frequencies ω_rf ±Δω_rf / 2 are excited. For a perfectly linear gradient G_z (i. e., a gradient over a homogeneous background field), these spins form perfect planes (illustrated at the left hand side of the following figure for three different excitation frequencies). In the presence of inhomogeneities, spins are no longer excited in linear planes, but in often very nonlinear structures (illustrated by the central slices on the right hand side).

Now look at the central, most severely deformed slice: During slice selection, all Larmor frequencies in this slice (independent of x) were exactly centered around ω_rf (only these Larmor frequencies were in resonance with the excitation pulse). However, after switching off the slice-selection gradient (as shown in the next figure), the Larmor frequencies in this slice are no longer the same at all positions: Obviously, the frequencies close to the center are higher due to the field inhomogeneity.

So, a first effect of field inhomogeneities are deformed slice geometries. A second effect results from the remaining field inhomogeneity in the excited slice shown above. This second effect becomes relevant during frequency encoding, when the spin localizations, x, are encoded by Larmor frequencies that are assumed to depend linearly on x: ω(x) = γx G_x. However, in the selected slice, the position x is no longer proportional to the Larmor frequency, ω(x), as shown in the following figure.

The deviation from the linear dependency shifts voxel locations in readout direction resulting in distortion artifacts. This effect (but not the deformed slice geometry) is (at least approximately) corrected by the view-angle tilting (VAT) MRI pulse sequence suggested by Z. H. Cho et al. (1988).

Mathematically, these effects can be analyzed using the basic signal equation of MRI. Without field inhomogeneities, the signal S(t) of the magnetization density ρ(x,z) during the readout is:

Frequently ignored, but here made explicit, is the integral over the slice thickness Δz from z₁ to z₂ = z₁ + Δz. As written above, the position and thickness of the slice are defined (in the case without field inhomogeneities) by the slice-selection gradient G_z and by the frequency (ω_rf) as well as bandwidth (Δω_rf) of the excitation rf pulse:

The signal equation also shows that each position x is uniquely encoded by a corresponding Larmor frequency ω(x) = γx G_x.

With additional field inhomogeneities, the Larmor frequencies change proportional to ΔB(x,z): during excitation, the frequencies are ω(x,z) = γ(z G_z + ΔB(x,z)) and, during readout, they are ω(x,z) = γ(x G_x + ΔB(x,z)). In both cases, the field inhomogeneity can be directly transformed into spatial shifts:

and

In both cases, the distortions (Δz(x,z) and Δx(x,z)) decrease when stronger gradients (G_z, G_x) are used.

So, in the presence of field inhomogeneities, the signal equation must be changed accordingly by adding an x-dependency to the slice positions z₁ and z₂ (describing, e. g., the curvy slice shape shown above) as well as using the modified readout frequency:

The slice positions z_1,2(x) – and, more generally, the position(s) z(x) corresponding to an excitation frequency ω_rf – are given by the following implicit equation (for each value of x):

which can give very nonlinear solutions for z(x), since ΔB(x,z) can be a complicated function of x and z. The slice can even split into several components centered at different z positions, which means that the z integral is not over a single interval, but over several intervals.

The view-angle tilting sequence (shown in the following figure) adds an additional readout gradient in slice-selection direction with the same magnitude G_z as the original slice-selection gradient.

The additional gradient during readout must be included into the signal equation as additional component γG_z z of the Larmor frequency:

For the following, a significant approximation is made: The complicated implicit equation for the slice position, z(x), given above is substantially simplified by replacing the x-dependent slice position z(x) as argument inside ΔB(x, z(x)) with the original slice position, ẑ:

The last part of this equation can be easily solved, giving:

where the relation between excitation frequency, ω_rf and original slice position, ẑ, was used in the final step: γẑG_z = ω_rf, i. e., ẑ = ω_rf / (γG_z).

So, in this approximation, the slice is simply shifted by  – ΔB(x, ẑ) / G_z and the slice thickness remains unchanged, i. e.

(A slightly less approximated calculation is presented by K. M. Koch et al. (2010), who introduce a first-order Taylor series:

where D(x, z) = ∂ΔB(x,z) / ∂z is the first derivative of the inhomogeneity function. Then

yielding

This describes shifting and scaling (by 1 + D(x,ẑ) / G_z) of the slice variable, z, and, thus, includes now scaling of the slice thickness.)

Using the simpler (zeroth-order) version z_1,2(x) = ẑ_1,2 – ΔB(x, ẑ) / G_z, the signal equation reads

and this can be transformed (by substituting z′ = z – ΔB / G_z) to

The innermost parentheses now contain the desired correction term that compensates the shifts in readout (x) direction:

This last equation describes a readout in a tilted readout direction (not longer along G_x, but along the vector (G_x, G_z)), without any distortions (which were related to the canceled term ΔB / G_x) along the readout direction.

Assuming that ρ(x, z′ – ΔB / G_z) does not vary significantly over the the z′ integration interval (say, a thin slice), it can be replaced by ρ(x, ẑ – ΔB / G_z). The integral over z′ can now be performed independently of ρ, resulting in a sinc-shaped modulation of the signal over the time t or, after reconstruction, a certain amount of blurring.

If additional corrections of the slice distortions along the z axis are required, techniques such as “Slice Encoding for Metal Artifact Correction” (SEMAC) as proposed by Lu et al. (2009) can be used in combination with the VAT approach (however, at the cost of substantially prolonged scan times).

Why are T₁ and T₂ of water what they are?

Relaxation is one of the most basic and obvious effects in nuclear magnetic resonance (NMR) experiments. Relaxation describes how the nuclear magnetization returns to its equilibrium state after an excitation by a pulsed radio-frequency field. The transverse magnetization (precessing in the plane perpendicular to the external magnetic field B₀ and yielding the measurable NMR signal) decays exponentially with a time constant called T₂; the longitudinal magnetization (parallel to B₀) recovers exponentially with a time constant called T₁.

For pure (distilled) water, these relaxation times are in the order of 1 to 5 seconds as established by numerous experiments. However, the calculation of these relaxation time constants is known to be notoriously complicated and difficult (requiring detailed quantum mechanical analysis of the spin ensembles and the spin interactions); a typical result of such an analysis is shown in the figure below.

So I wondered if the NMR relaxation time constants T₁ or T₂ of pure water can be roughly estimated using only some simple assumptions and considerations such that the resulting quantities are at least in the correct order of magnitude.

The main reason for transverse and longitudinal relaxation is the magnetic dipole-dipole interaction between spins (other effects such as interaction with the thermal radiation field or electric interactions with the electrons can be neglected for the water protons). The dipole field of a magnetic moment m is

or, which is enough for our order-or-magnitude estimations,

where μ₀ = 4π × 10^–7 T² m³ / J is the vacuum permeability.

The magnetic moment m of the proton can be easily expressed by the gyromagnetic ratio γ ≈ 267.5 × 10⁶ rad / (sT), which is defined as the ratio between magnetic moment and angular momentum I, m = γI. Approximating the angular momentum I of the proton by I ~ ℏ ≈ 1.055 × 10^–34 Js leads to |m| ~ γℏ ≈ 3 × 10^–26 J / T and

A first idea about T₂ or transverse relaxation (introduced already 1946 in the seminal NMR paper by F. Bloch) is to consider spins (protons) in a rigid lattice, where the magnetic field of the nearest neighboring proton induces a phase shift at the spin of interest leading to signal dephasing. Using the shortest distance that is relevant for the nuclei in water molecules, namely the intramolecular proton-proton distance d_pp ≈ 0.15 × 10^–9 m, the additional, dephasing dipole field of a single neighboring proton is δB = B_dipole(d_pp). This yields a dephasing angular frequency δω = γδB = γB_dipole(d_pp) (in addition to the Larmor frequency due to the external B₀ field). Assuming that transverse relaxation requires an additional phase angle of the order of 1, the time constant T₂ can be estimated to be of the order

Well, this value for T₂ ~ 5 μs is obviously way too small for water – by about 6 orders of magnitude! So we are not even close … but there is still something to learn from this: The very short T₂ is in fact quite exactly what we find in solids (as e.g. in water ice with T₂ ≈ 4 μs as measured by T. G. Nunes et al.), which shouldn’t be too surprising since we have assumed a rigid spin configuration above. Consequently, we’ve learned now that the constant random motion of water molecules in liquid water must play a very important role for the quantitative understanding of relaxation.

To obtain at least an order-of-magnitude estimation of the actual T₂ relaxation time (and in fact also of T₁) of liquid water, we have to include some additional information about the random motion of water molecules. A handy physical parameter to describe the effects of this random motion is the correlation time τ_c, which can be used to describe, e. g., a fluctuating magnetic field component ⟨B(t)B(t + τ)⟩ = ⟨B(t)²⟩ exp(–|τ| / τ_c), where ⟨ · ⟩ denotes the ensemble average over all spins and the result is assumed to be independent of the time t. So basically, after the correlation time τ_c the fluctuating field strength becomes statistically independent of its former values.

For the protons of water, this correlation time is extremely short, in the order of picoseconds: τ_c ≈ 5 × 10^–12 s. This value is in fact about the diffusion time of a water molecule over a region with diameter d_pp, which is d_pp² / (6D) ≈ 2 × 10^–12 s for D ≈ 2 × 10^–9m² / s. It is important to note that there is almost no dephasing over the correlation time, i. e., τ_c δω ≈ 10^–6 ≪ 1. This means that the static dephasing as assumed above cannot take place since the protons move too quickly and the static effects are averaged out. After about τ_c, the motion has completely changed the orientation of the vector distance r between the protons and, hence, also the magnetic dipole field direction and the magnitude of each field component.

To obtain an estimation of T₂ using this correlation time, we now consider the phase angles ϕ(t) within our ensemble of spins after a time t – more exactly, the additional phase angles (as seen in a rotating frame of reference) due to the fluctuating random magnetic fields δB(t). These phase angles exhibit a random distribution with mean value ⟨ϕ(t)⟩ = 0, but the width of this distribution ⟨ϕ(t)²⟩ increases proportional to the time ⟨ϕ(t)²⟩ = αt. Similarly as above, we estimate our relaxation time T₂ to be the time required for this standard deviation to become of order 1; i. e., T₂ ~ (1)/(α) = (t)/(⟨ϕ(t)²⟩).

What can we say about the value of α? The factor α must have the dimension of 1/time; it should increase with the local magnetic field strength δB = δω / γ and also with the correlation time τ_c (the longer the correlation time, the less effective becomes the averaging). So, the simplest approach to express α by the given quantities with the correct physical dimension is to set α = δω²τ_c. With this guess, we find

in surprisingly good agreement with the experimental data.

To obtain a very similar result by an actual calculation, we can quantitatively determine the width of the distribution ⟨ϕ(t)²⟩. The relation between the fluctuating magnetic field component δB(t) seen by an individual spin and the phase evolution of this spin is

yielding

or, by substituting τ = t″ – t′,

We can now use the correlation time relation from above for the fluctuating random magnetic field component δB(t), namely ⟨δB(t) δB(t + τ)⟩ = ⟨δB(t)²⟩ exp(–|τ| / τ_c) with the average squared field strength ⟨δB(t)²⟩ ~ B_dipole(d_pp)² ~ δω² / γ² to obtain

To solve this double integral without long calculations (which are nevertheless straightforward), we need the fact that τ_c ≪ t and therefore almost always τ_c ≪ τ (the second integration variable), so the exponential function in the integral is almost always 0.

Hence, the second (inner) integral can be approximated for almost all values of t′ as ∞∫–∞exp(–|τ| / τ_c) dτ = 2∞∫0exp(–τ / τ_c) dτ = 2τ_c. With this approximation, we find

and, consequently,

which is again an estimation that agrees well with the observed values.

A final short comment on T₁: In non-viscous liquids such as water, T₁ and T₂ are approximately the same. While loss of phase coherence is required for T₂ relaxation, energy transfer from individual spins to the liquid is causing T₁ relaxation. This energy transfer occurs due to transverse fluctuating magnetic fields δB_x,y(t) with spectral components agreeing with the Larmor frequency ω = γB₀. Since these frequencies are much lower than the inverse correlation time, i. e., γB₀ ≪ 1 / τ_c, they basically occur with the same probability as the dephasing T₂ effects considered above. A more detailed analysis of T₁ and T₂ requires quantitative estimations of the spectral density functions J(ω) of the dynamic processes involved.

Update (2019-02-29)

A concise introductory summary about some approaches and results of quantitative relaxation theory can be found in the text book “Principles of Nuclear Magnetic Resonance Microscopy” by Paul T. Callaghan, section 2.5.

What’s special about 1.256 431 208 626… in MRI?

Consider the following simple question in magnetic resonance imaging (MRI) or spectroscopy (MRS): Given a fixed total measurement time, T_total, (e. g., a typical breath-hold duration of 16 s or a maximum accepted sequence duration of 10 min) and the possibility to fit into this duration an acquisition with several repeated read-outs (that are to be averaged to increase the data quality), what is the optimum balance between the repetition time (TR or T_R) and the number of averaged signals (also called simply the number of averages, N, or the number of excitations, NEX or N_ex)? (For simplicity, let’s consider only pulse sequences with 90 ° excitations (e. g., spin-echo sequences), in which all available magnetization is flipped to the transverse plain in each repetition.)

The solution (for 90° excitations)

In general, a single MRI acquisition (without averaging) requires M repeated excitations separated by T_R (e. g., M could be the number of phase-encoding steps for conventional spin-echo acquisitions or M = 1 for single-shot EPI or one-dimensional MR spectroscopy). Hence, the acquisition time without data averaging is T_acq = MT_R, and, neglecting for now that the number of averages must be integer, this number of averages is N = T_total / T_acq = T_total / (MT_R). In the following, the actual relevant time parameter is, thus, T_avail = T_total / M, the time available for each required excitation, and the number of averages is N = T_avail / T_R. We can also express the repetition time by these parameters as:

On the one hand, averaging data increases the obtainable signal-to-noise ratio (SNR or Ψ) proportional to the square root of the number of averages:

On the other hand, increasing the number of averages (at fixed total measurement duration!) results in a shortened time for longitudinal relaxation: T_R = T_avail / N and, thus, in less SNR:

The resulting total SNR is proportional to the product of both factors, i. e.,

The SNR (as a function of T_R at constant total measurement time) has, thus, a maximum in an intermediate range; for very short T_R (and, consequently, many averages), the longitudinal magnetization cannot relax sufficiently, which reduces the available signal; for very long T_R, the relaxation is approximately complete anyway, but the number of averages decreases.

To calculate the optimum T_R (or N), T_R and T_avail are best expressed in terms of T₁ using τ = T_R / T₁ and T = T_avail / T₁, which gives

The maximum of this expression with respect to τ is obtained by setting its derivative to zero, i. e.

which gives

The solution of this equation is

where W_–1 is the lower branch of the Lambert W-function (or product-log or Omega function, which gives the solution for W in z = Wexp(W), i. e., it is the inverse function of f(W) = Wexp(W)). Hence, the optimal choice for T_R is (at least theoretically):

Fortunately, the function to be maximized has a rather broad maximum, and one obtains values between 99 % and 100 % of the maximum SNR for T_R between about 0.9937 T₁ and 1.5773 T₁ and values between 95 % and 100 % of the maximum SNR for T_R between about 0.7293 T₁ and 2.0882 T₁ (values determined by numerical evaluation and illustrated in the first figure at the top).

So, practically, choosing T_R to be T₁ gives still a nearly optimal SNR; if a larger number of signal acquisitions is preferred, e. g., for statistical evaluation, then T_R may even be shortened to 0.75 T₁, and if – on the other hand – the base SNR of every single acquisition becomes very low (but is to be reconstructed before signal averaging) then T_R may be chosen to be 2 T₁ to increase the quality of each single (not averaged) data set.

The presented result is not new; in fact, it has been derived before (at least numerically) in several publications. An early example is, e. g., a publication by R. R. Ernst and R. E. Morgan (1973).

The presented analysis is valid only for simple 90° excitations; FLASH sequences or steady-state pulse sequences which refocus the transverse magnetization as well will exhibit a different behavior.

Of course, optimizing T_R as described above is not an option if T₁-weighted or T₂-weighted images are to be acquired. In these cases, the optimal value of T_R depends on the desired contrast and not on the total SNR.

(This is a shortened version of a slightly longer pdf document that contains some more details and discusses also a few special cases.)

Looking for reference values of water diffusion?

This is a short advertisement for my web document “Diffusion Coefficients of Water”. If you have ever looked for reference values of the self-diffusion coefficient of water at different temperatures, then chances are that you might like to bookmark the link above.

(If you have never heard of molecular (self-)diffusion, then it’s probably not so relevant for you. As an ultra-short explanation: Self-diffusion is a term for the thermal motion of molecules. If we could look at the individual molecules of a liquid such as water (say, with a really big microscope), then we would see all molecules in constant random motion. This molecular motion is known as the cause of the Brownian motion, which has been observed already in 1827 by the botanist Robert Brown. The extent of this motion is directly related to the temperature of the liquid or – which is essentially the same – to the thermal kinetic energy of the molecules. This motion can be described by a physical quantity called the (self-)diffusion coefficient D, which has units of m² / s. The meaning of a diffusion coefficient of, e. g., D = 2 × 10^–3 mm² / s is that after a time t = 1 s the space “visited” by an average diffusing molecule (in a certain statistical sense) has a size s (i. e., a radius) of s = √( 6Dt ) = √( 12 × 10^–3 mm² ) ≈ 0.1 mm.)

Why is it good to know the exact relationship between the temperature and the diffusion coefficient of water? Typical applications are:

• Using these data, we can determine temperatures by measuring diffusion coefficients of water, e. g., using diffusion-weighted MRI.
• We can check how well our MRI system is calibrated by measuring the temperature conventionally (i. e., with a thermometer) and compare the MRI-determined diffusion coefficient with a reference value.
• An elegant variant of the previous approach is to use an ice-water phantom (i. e., a sample of liquid water floating in a mixture of ice and water such that the temperature is known to be 0 °C) as suggested in a publication by T. L. Chenevert et al. (2011).

The web document introduced above provides an interactive interface to calculate self-diffusion coefficients of water at different temperatures (or, alternatively, to calculate the temperature corresponding to a given diffusion coefficient). This calculation is based on the results taken from several published articles with exactly this kind of data (namely, measured diffusion coefficients of water at different temperatures).

The web page has actually existed for quite some time (since about 2002), but it was never made public. No (hyper)links to it existed and it was used only by me (and by some colleagues and our Master or PhD students). Thus, I was somewhat surprised when I recently found out that the page is indexed by google and has already been cited in several reports and theses. Consequently, I decided to slightly update its content (and to pretty up its appearance). The most relevant change is that I added a non-linear (quadratic) fit that describes the data in the Arrhenius plot much more accurately (and can be used for data over a greater temperature range) than the previously used linear fit. Note that using the (recommended) quadratic fit, calculated diffusion coefficients for intermediate temperatures around 15 to 30 °C will differ from those of the earlier (linear) version by a few tenths of a percent.

Double field strength gives double SNR – or not?

Some time ago at lunch, we had a discussion about the advantages of high magnetic field strengths B₀ in MRI. We happily agreed that higher field strengths result in higher signal-to-noise ratios (SNR). But then several opinions surfaced about the exact quantitative relation between the SNR and B₀ – ranging from linear to quadratic and including some very specific exponents in between such as 7/4. It turns out that more than one correct answer exists … and there are some surprising technical subtleties.

As starting point of a quantitative discussion, we have to define what we consider as SNR in MRI: The SNR (sometimes given the symbol Ψ) is defined as the amplitude of the image signal (at some point or region of interest) divided by the standard deviation of the noise signal. (There are other SNR definitions used in other disciplines that involve squared amplitudes (signal power) and logarithms (resulting in decibels), but in MRI we generally use the very simple ratio of the signal amplitude and the noise standard deviation.)

To simplify things, we can divide the analysis of field strength and SNR into two parts: the discussion of the signal and of the noise. The signal part is rather straight forward: If the field strength B₀ increases, then

• the Larmor (precession) frequency ω = γB₀ of the nuclei increases proportionally with B₀ (the proportionality constant γ is the gyromagnetic ratio)
• and the nuclear magnetization M_N increases; for realistic field strengths and “normal” temperatures (around 300 K), the nuclear magnetization is approximately proportional to the field strength M_N ≈ χ_v B₀, where χ_v is the nuclear magnetic (volume) susceptibility given by χ_v = ((N / V) γ² ℏ² I (I + 1))/(3kT) (with N / V: spin density, I: nuclear spin quantum number, k: Boltzmann constant, T: sample temperature).

The measured signal S is the voltage induced in the radio-frequency (rf) receiver coil by the precessing nuclear magnetization. As known from the theory of electromagnetic induction, this induced voltage is proportional to both the (angular) frequency ω = γB₀ and the magnetization M_N = χB₀, and taking both factors together, we find S ∝ B₀².

So, if we stop at this point, we may hope for four times the SNR at double field strength. But we haven’t considered the noise yet. And there the trouble begins …

There are several sources of noise in MRI and depending on the setup, different sources can dominate the noise generation. Generally, the thermal noise voltage U_noise can be expressed as U_noise = √( 4kT R Δf ), where Δf is the signal bandwidth and R is a resistance associated to the rf receiver coil. For very small samples (milliliters), this resistance R is dominated by the actual coil resistance R_coil. R_coil is proportional to the square root of the Larmor frequency R_coil ∝ √( ω ) ∝ √( B₀ ) because of the rf skin effect, which reduces the penetration depth and, thus, the conducting area of the coil wires proportional to 1 / √( ω ). Therefore, U_noise ∝ R_coil^1 / 2 ∝ ω^1 / 4 ∝ B₀^1 / 4 and the resulting SNR is proportional to S / U_noise ∝ B₀^7 / 4 in this case.

However, a different looking result is presented by A. Abragam in his classic monograph “The principles of nuclear magnetism” (1961). There (on p. 83) we find for the SNR Ψ ∝ √( Q ) B₀^3 / 2 with the coil quality factor Q. This is derived from the shunt resistance R = QLω of a circuit (here the receive coil) with inductance L. This apparently contradicting result can be explained if the frequency dependence of the quality factor Q ∝ √( ω ) (for a Q-optimized solenoid coil) is considered, which then gives the same exponent of 7 / 4 as before.

Finally, for larger samples – such as human subjects in clinical MRI – additional noise is generated because of the inductive (or magnetic) losses associated with the electrical conductivity of the sample or tissue. Electrical power is dissipated in the sample proportional to the squared induced voltage in the sample, i. e., proportional to ω² and proportional to the squared current I² in the coil. This dissipation of power can therefore be expressed as an additional apparent coil resistance R_sample that is also proportional to the squared Larmor frequency R_sample ∝ ω². Consequently, considering only this apparent coil resistance associated to the sample, we find U_noise ∝ R_sample^1 / 2 ∝ ω ∝ B₀, and now the resulting SNR is proportional to S / U_noise ∝ B₀.

These results (and detailed derivations) can be found in two classic MRI papers by D. I. Hoult and R. E. Richards (1976) and by D. I. Hoult and P. C. Lauterbur (1979).

An obvious question now is: Can’t we simply measure the SNR dependence on the field strength to find out what’s actually going on in MRI? Well, we can try, but there are several complications. One major problem is that we cannot use the same rf coils at different field strengths. There will generally be differences of the coil design for different field strengths and these must be expected to influence the measured SNRs. Another complication is the influence of the electromagnetic wavelengths in tissue that approach the dimensions of the samples at about 3 to 7 T. Furthermore, the relaxation times (T₁, T₂, T₂^*) change with B₀ and may influence SNR measurements.

Nevertheless, there are some publications reporting SNRs at different field strengths, e. g. by D. I. Hoult et al. (1986), J. T. Vaughan et al. (2001), C. Triantafyllou et al. (2005), and recently by R. Pohmann et al. (2016). They all show a clear increase of SNR with B₀, but there are still some discrepancies with respect to the exact dependency. Particularly the latest paper by R. Pohmann et al. demonstrates and discusses a slightly better than linear increase of SNR at high fields. But in conclusion and as rule of thumb, assuming an approximately linear relationship of SNR and B₀ appears still justified for large (clinical) MRI systems.

Is T₂ necessarily shorter than T₁?

If you are familiar with MRI (or NMR in general), then probably also with the relaxation time constants T₁ and T₂. These tissue-specific (or substance-specific) constants describe how fast the nuclear magnetization returns to its equilibrium value, M₀, after excitation by a pulsed radio-frequency (rf) field. Shortly summarized, T₁ describes the exponential recovery of the longitudinal magnetization M_L (i. e., of the magnetization parallel to the external static magnetic field B₀), and T₂ describes the exponential decay of the transverse magnetization M_T (that is precessing in the plane orthogonal to B₀).

Typically (as shown in the first figure), T₂ values of tissue are considerably lower than T₁ values, i. e., the transverse magnetization decays quicker than the longitudinal relaxation needs for recovery. For most tissues in vivo, T₁ varies between about 300 ms and 3 s, while T₂ varies between about 10 ms and 200 ms. Longer T₂ relaxation times (up to about 3 s as well) are found for liquids.

So one may ask if there are good physical reasons for T₂ values being shorter than or – at most – equal to T₁. As a physicist, I’d start with checking some extreme cases, e. g., assuming that T₂ is much longer than T₁, i. e. T₂ ≫ T₁. Then the longitudinal magnetization can fully recover while at the same time some transverse magnetization would be preserved. As a result, the magnitude of the total magnetization √( M_L² + M_T² ) would become greater than the equilibrium value M₀ – which is physically impossible.

To analyze these properties of T₁ and T₂ in more detail, longitudinal and transverse relaxation can also be plotted together in a diagram showing the transverse magnetization on the horizontal axis and the longitudinal magnetization on the vertical axis. These diagrams show the evolution of the magnetization (for experts: in a rotating frame of reference) after a 90° rf pulse (all trajectories start at the lower right corner of the diagram). The left-hand side of the following figure shows three cases T₂ = T₁ / 2 (blue), T₂ = T₁ (green), and T₂ = 2 T₁ [sic!] (cyan); and all three curves show a “benign” behavior in that they lie in the shaded area below the black circle segment M_L² + M_T² = M₀². This means that the magnitude of the total magnetization vector is always smaller than M₀.

However, the right-hand side of this figure shows what’s happening if T₂ becomes greater than 2 T₁ – in this case, T₂ = 3 T₁ (red curve): Now the magnitude of the total magnetization vector increases above the physical limit of M₀, i. e., the red line crosses the black border of physically benign behavior!

In fact, it can be shown (see appendix if interested) that the maximum T₂ value, for which the red curve stays always below the black line, is exactly T₂ = 2 T₁. And as so often, almost everything that is physically possible is also realized in nature (although the case T₁ < T₂ < 2 T₁ is really extremely rare), as described by Malcolm H. Levitt in his highly recommendable NMR text book “Spin dynamics” (2nd ed., section 11.9.2, note 13):

The case where T₂ > T₁ is encountered when the spin relaxation is caused by fluctuating microscopic fields that are predominantly transverse rather than longitudinal. One mechanism which gives rise to fields of this form involves the antisymmetric component of the chemical shift tensor (not to be confused with the CSA). […] Molecular systems in which this mechanims is dominant are exceedingly rare (see F. A. L. Anet, D. J. O’Leary, C. G. Wade and R. D. Johnson, Chem. Phys. Lett., 171, 401 (1990)).

So, the answer to the title question is: No, T₂ can in fact be greater than T₁ in very special circumstances, but it can never be greater than 2 T₁.

Appendix

The maximum T₂ value, for which the red curve stays always below the black line, is exactly T₂ = 2 T₁. This can be seen by analyzing the inequality

First, we divide by M₀² and set T₂ = αT₁ as well as β = exp(–t / T₁), yielding

which is (after subtraction of 1 and division by β)

β is by definition (for positive t) between 0 and 1, so the right-hand side of the last inequality is a linear function descending from 2 to 1 (i. e. always  ≤ 2). Its left-hand side has very different shapes depending on α: it is increasing from 0 to 1 for 0 < α ≤ 2 (since then the exponent 2 / α – 1 ≥ 0); but it is going to infinity for β → 0 if α > 2 (since then the exponent 2 / α – 1 < 0). So, the last inequality will not hold in the latter case for sufficiently small values of β, which means that non-physical behavior occurs if α > 2 or, using the definition from above, if T₂ > 2 T₁.

Gadolinium does not increase the signal, but the speed of MRI

You may have heard that gadolinium-based MRI contrast agents can enhance or increase the signal of tissue. This is generally a good description of what’s going on. Here, however, I would like to argue why this is – strictly speaking – not true: contrast agents cannot really increase the signal available for MRI.

Some basic facts first: gadolinium-based contrast agents are very frequently used in clinical MRI to improve the image contrast as illustrated in the following example.

An important fact about (conventional) MRI contrast agents is that it’s never the contrast agent itself that is visible in MR images. Instead, the contrast agent changes the behavior of the atomic nuclei in its neighborhood – in clinical MRI, these are the nuclei of hydrogen, i. e. the protons. As a consequence, these protons now appear brighter in T₁-weighted MRI than protons which are not influenced by the contrast agent.

So, apparently the proton signal is increased by gadolinium? Yes, apparently … Actually, there is always a maximum signal that is available for MRI and that depends on three major factors:

• the number of available protons, which is related to the proton density ρ of the tissue: the more protons (per voxel), the higher the signal;
• the magnetic field strength B₀: the higher the field strength, the higher is also the (thermal) nuclear magnetization and, hence, the measured signal;
• the receiver coil: the more efficient the receive system (the radio-frequency coil), the higher the signal.

But the presence of a contrast agent does not increase this maximum signal.

Instead, we are cheating: First, we artificially decrease the MRI signal by choosing short repetition times (TR). And only afterwards, a certain part of this suppressed signal is recovered due to the influence of the contrast agent! This is illustrated in the following diagram:

Obviously, the maximum signal, S_max, with contrast agent is exactly the same as the maximum signal without contrast agent – but we can obtain this maximum signal considerably faster (i. e, at shorter TRs). That’s why gadolinium can be described to increase the speed, but not the MRI signal. In agreement with this observation, no additional gadolinium-induced signal enhancement can be found in proton-density-weighted MR images (with very long TRs). But, hypothetically, if contrast agent could be distributed homogeneously in the tissue (which in reality is not possible), then PD-weighted MRI could be accelerated by using shorter TRs without changing the contrast.

Simultaneous multi-slice excitation done differently

Last week, Koray Ertan and colleagues published an MRM paper (doi: 10.1002/mrm.27031) (expanding on an ISMRM abstract presented earlier this year) with an interesting approach for simultaneous multi-slice excitation. So, what’s it about?

In MRI, we are frequently interested in data from a single two-dimensional slice of the imaged subject or object – or actually in data from several such slices. These slices are then displayed as conventional 2D MR images. A procedure called slice selection is used to restrict our data acquisition to each single slice. To understand slice selection, one has to know that MRI is based on the resonant excitation of spins in a static magnetic field B₀. Spins have a characteristic, so-called Larmor frequency ω = γB depending on the magnetic field B (the constant of proportionality γ is the gyromagnetic ratio). By applying a radio-frequency (rf) field with exactly this Larmor frequency, spins can be excited (i. e., they can be made to generate a measurable signal).

The basic idea of slice selection is to excite only spins in a single slice by (first) superposing a linear magnetic field gradient g_z e. g. in z direction, resulting in the spatial field distribution (shown in the first figure):

(These linear magnetic fields or gradient fields are one of the basic ingredients of MR imaging. Each MR imager comes with built-in coils to apply gradient fields in all possible spatial orientations.)

Then, by choosing an rf frequency ω_slc = 2πf_slc, we can select a plane in space where

and only spins in this plane centered around z_slc = (ω_slc / γ – B₀) / g_z are excited (see first figure).

More advanced MRI techniques can excite several spatially separated slices at once by applying rf fields with more than one frequency – the difficult part is then to separate data acquired from these slices for reconstruction.

The nice idea realized in the paper by Koray Ertan et al. is to excite multiple slices not by applying rf fields with a mixture of several frequencies, but instead by modifying the gradient field to a spatially non-linear magnetic field. Consequently the mapping between frequencies ω and spatial positions (z) is no longer one-to-one, but several positions can correspond to a single frequency:

So, depending on the shape of the magnetic field variation two or more slices can be excited using a single excitation frequency. This has the advantage that we can use short and simple standard rf pulses for slice excitation. The obvious disadvantage is that it requires additional gradient hardware providing the non-linear magnetic fields, which is currently not available at existing MRI systems.

However, this may change in future, since there are currently some promising approaches for non-linear encoding fields. In particular, I’m thinking about an impressive MRM paper (doi: 10.1002/mrm.26700) by Sebastian Littin and colleagues published also this year, in which an 84-channel matrix gradient coil is presented, which is capable of providing very flexible linear or non-linear field configurations.