Research ArticleCHEMICAL PHYSICS

Quantum mechanical MRI simulations: Solving the matrix dimension problem

See allHide authors and affiliations

Science Advances  19 Jul 2019:
Vol. 5, no. 7, eaaw8962
DOI: 10.1126/sciadv.aaw8962

Figures

  • Fig. 1 Accuracy and run time statistics for a Stejskal-Tanner pulse sequence simulation followed by extraction of the known diffusion coefficient (18 × 10−10m2/s, dashed line) from the simulated data.

    The 1.5-cm-long sample contains a single type of protons with a chemical shift of 4.6 ppm at 11.74 T. The Stejskal-Tanner pulse sequence simulation uses ideal RF pulses; perfectly rectangular gradients are assumed with no stabilization delay. The duration of gradient pulses is δ = 2 ms, and the duration of the diffusion delay is Δ = 50 ms. Gradient amplitudes varied from 0 to 0.5 T/m. (Top) The diffusion coefficient extracted by fitting Eq. 34 to the simulated data for three finite difference stencil sizes as a function of the spatial grid size. (Bottom) Wall clock time (2 × Intel Xeon E5-2698) for the pulse sequence simulation for three finite difference stencil sizes as a function of the spatial grid size.

  • Fig. 2 Simulated Stejskal-Tanner diffusion attenuation profiles as a function of spatial grid size.

    The 1.5-cm-long sample contains a single type of protons with a chemical shift of 4.6 ppm at 11.74 T. The Stejskal-Tanner pulse sequence simulation uses ideal RF pulses; perfectly rectangular gradients are assumed with no stabilization delay. The duration of gradient pulses is δ = 2 ms, and the duration of the diffusion delay is Δ = 50 ms. The diffusion coefficient is 18 × 10−10 m2/s. (Top) Diffusion attenuation profiles for spatial grids and finite difference stencils of different sizes (the minimal grid that satisfies the spatial Nyquist condition in this system has 1280 points). (Bottom) Difference between the simulated diffusion attenuation profiles and the exact analytical answer for grids and finite difference stencils of difference sizes.

  • Fig. 3 DPFGSE pulse sequence element (67) setup to defocus the signals influenced by the pair of soft pulses.

    In our experience, this is the best solvent suppression method in existence. The sequence may be switched from selective suppression to selective excitation by turning off the hard 180° pulses.

  • Fig. 4 Examples of DPFGSE simulations.

    (A) Simulated 1H NMR spectrum of 20 mM GABA (a common metabolite) dissolved in water; water signal dominates the spectrum. (B) The same spectrum simulated in the presence of double pulsed field gradient spin echo (DPFGSE) water suppression. Amplitudes of first and second DPFGSE gradients are 0.1 and 0.15 T/m, respectively; gradient durations are all 1.0 ms; sample size is 15 mm; there are 500 points in the spatial grid; off-resonance selective inversion pulses at the frequency of the water signal are simulated with the Fokker-Planck method (51) and use 10-point Gaussian envelopes. (C) Same as (B), but including spatial diffusion with a diffusion coefficient of 2.6 × 10−9 m2/s. (D) Same as (C), but using DPFGSE signal selection rather than suppression mode. One of the GABA multiplets is selectively excited. (E) Same as (C), but also including spatial flow at 10 mm/s. All simulations run in minutes on a Tesla K40 graphics processing unit (GPU) using Spinach 2.4 and later, and are included in the example set.

  • Fig. 5 PRESS simulations for a 1D sample containing varying concentrations of three different two-spin systems.

    (Top) Concentration profiles of the three spin systems—system A, δ = ±1 ppm, J = 30 Hz; system B, δ = ±3 ppm, J = 10 Hz; system C, δ = ±2 ppm, J = 20 Hz. (Middle) Volume selection profiles excited by an off-resonance square pulse at three different frequencies. (Bottom) Magnitude mode PRESS NMR spectra of the three selected volumes. All simulations run in seconds on a Tesla K40 GPU using Spinach 2.4 and later, and are included in the example set.

  • Fig. 6 Simulation of PRESS excitation hotspots.

    The sample on the left has two dimensions (108 × 90 pixels), and the sample on the right has three (108 × 90 × 111 pixels). PRESS hotspots in two (108 × 90 pixels, left) and three (108 × 90 × 111 pixels, right) dimensions. All parameters and J-coupled spin systems are the same as in Fig. 5. Although its effect is not immediately visible, the 3D simulation includes isotropic diffusion with D = 2.6 × 10−9 m2/s to emphasize the ability of the formalism presented here to handle spatial dynamics in 3D. Gradients are tilted relative to the spatial grid by arbitrary angles. Both simulations (including volumetric visualization in the right panel) run in minutes to hours on a Tesla K40 GPU using Spinach 2.4 and later, and are included in the example set.

  • Fig. 7 PSYCHE (pure shift yielded by chirp excitation) strong coupling artefact simulation.

    The figure is zoomed into the central parts of the simulations of 500 MHz NMR (bottom) and PSYCHE (69) (top) spectra of rotenone (22 proton spins). The simulation runs in minutes on a Tesla K40 GPU using Spinach 2.4 and later. The standard NMR simulation is nearly instantaneous.

Tables

  • Table 1 Wall clock time benchmarks for matrix-vector multiplication using a single Xeon E5-2698 CPU core in Matlab on a computer equipped with 256 GB of RAM.

    Matrix-vector multiplication taskWall clock time, polyadic repWall clock time, explicit rep
    [AB]v with dim(A,B) ≤ 64, full0.37 ± 0.01 ms0.88 ± 0.12 ms
    [ABC]v with dim(A-C) ≤ 64, full1.8 ± 0.3 msOut of RAM
    [ABCD]v with dim(A-D) ≤ 64, full97 ± 14 msOut of RAM
    [AB]v with dim(A,B) ≤ 64, sparse0.21 ± 0.01 ms0.05 ± 0.01 ms
    [ABC]v with dim(A-C) ≤ 64, sparse2.1 ± 0.3 ms11.4 ± 1.6 ms
    [ABCD]v with dim(A-D) ≤ 64, sparse105 ± 16 msOut of RAM
  • Table 2 Matrix dimension and memory footprint statistics for the simulation of DPFGSE water suppression and selective excitation NMR spectrum shown in the bottom row of Fig. 4.

    Memory utilization is quoted as reported by Matlab using compressed column sparse format (74).

    Problem parameterValueNotes
    Liouville space dimension, full4,11246 (GABA) +
    16 (water)
    Liouville space dimension, reduced1,912IK-2 basis (45)
    Min. points in the spatial grid500Spatial Nyquist condition
    Space(x)Spin dimension, full2,056,000500 × (46 + 16)
    Space(x)Spin dimension, reduced195,500IK-2 basis (45), ZTE (42)
    Nonzeroes in evolution
    generator, matrix
    18,252,000446 MB
    Nonzeroes in evolution
    propagator, matrix
    >109>16 GB
    Nonzeroes in evolution
    generator, polyadic
    28,418640 kB
  • Table 3 Matrix dimension and memory footprint statistics for the simulation of 3D PRESS excitation profile in the right panel of Fig. 6.

    Memory utilization is quoted as reported by Matlab using compressed column sparse format (74).

    Problem parameterValueNotes
    Liouville space dimension, full409646
    Liouville space dimension, reduced483 × 42
    Spatial grid, points108 × 90 × 111
    Space(x)spin dimension, full4.4 × 10946 × 108 × 90 × 111
    Space(x)spin dimension, reduced5.2 × 106(3 × 42) × 108 × 90 × 111
    Nonzeroes in evolution generator, matrix425,094,48011 GB
    Nonzeroes in evolution generator, polyadic5,403,402233 MB
  • Table 4 Evolution generator dimension (spin × spatial), number of nonzeroes (nnz), and matrix-vector multiplication wall clock time (wct) statistics for a series of simulations of ultrafast “maximum-quantum” experiments (70) that include the effect of isotropic spatial diffusion using seven-point central finite difference operators.

    SystemParameterValueNotes
    2-SpinState space dimension16 × 500Sparse matrix representation fits into the L2 cache—polyadic processing is slower,
    memory footprint improvement is not significant
    nnz(L) and wct(Lv)a, sparse116,000, 0.5 ms
    nnz(L) and wct(Lv)a, polyadic12,000, 5 ms
    4-SpinState space dimension256 × 500Sparse matrix representation exceeds L2 cache, parity on the wall clock time, significant
    memory advantage for the polyadic format
    nnz(L) and wct(Lv)a, sparse2,600,000, 13 ms
    nnz(L) and wct(Lv)a, polyadic14,000, 18 ms
    6-SpinState space dimension4096 × 500Sparse matrix representation is in the gigabytes, polyadic processing is faster and
    offers a vast memory footprint advantage
    nnz(L) and wct(Lv)a, sparse37,000,000, 227 ms
    nnz(L) and wct(Lv)a, polyadic87,000, 185 ms

    aIntel Xeon E5-2698, single core.

    Navigate This Article