Modelling MRS Spectra

A small language for writing down a spectrum, and an interpreter that draws what you wrote — on ¹H, ²H, ¹³C, ¹⁹F, ²³Na, ³¹P and ¹²⁹Xe

Almost everything an in vivo spectrum does can be written as one line of arithmetic. The free induction decay is a sum of complex exponentials, each one a resonance: an amplitude, a frequency, a starting phase, and a decay. The spectrum is its Fourier transform. Line broadening, zero filling, phase correction, echo time, decoupling, the difference between a doublet and a triplet — all of them are operations on that sum, applied before the transform or after it.

The honest version of that sentence is that the sum is what survives the physics. Getting from a spin Hamiltonian to a list of frequencies and amplitudes is the hard part, and for strongly coupled spin systems it cannot be done on the back of an envelope. What this page models is the easy part, which happens to be the part that describes most of what an in vivo spectrum looks like.Writing the sum out by hand is tedious in a way that hides the structure. So this page contains a small domain specific language for declaring it — nucleus, field, a list of resonances with their shifts and widths and couplings, an acquisition, a processing chain — together with a complete interpreter for that language: lexer, parser, compiler, and a simulator that synthesises the time-domain signal and transforms it. Nothing is fetched, nothing is uploaded, and the program you write is the whole of the model.

The reason to do this across nuclei rather than for protons alone is that the interesting differences between ³¹P, ¹³C and ¹²⁹Xe spectroscopy are almost entirely differences in three numbers, and a language makes that visible in a way a paragraph does not. The gyromagnetic ratio sets how many hertz a part per million is worth. The chemical shift range sets how many parts per million you have to spend. The coupling constants set whether a resonance arrives as one line or seven. Change those three things in the program below and a phosphorus spectrum becomes a carbon spectrum.

The workbench

The examples in the menu

The language in one page

A program is a sequence of statements in any order. Every statement has a fixed number of arguments, so several may share a line — which is what makes a one-line pool declaration legible — and an expression can never run past the end of a line. Comments run from # to the end of the line, except where a # is followed immediately by three or six hexadecimal digits, which is a colour.

The one place the grammar is stricter than it looks: a statement taking two numbers, such as couple, line and axis, reads each of them as a single term rather than as arithmetic, so that line -0.03 0.5 is two arguments rather than one subtraction. Put brackets around anything more complicated than a name or a number.Numbers may carry a unit — ppm, Hz, kHz, ms, us, s, T, deg — and the unit decides how the number is read, so width 8 Hz and width 0.15 ppm are both linewidths and mean different things at different fields. Times default to seconds, shifts to ppm, and linewidths and couplings to hertz, but it costs nothing to write the unit and the examples always do.

nucleus 31P
field 3 T
sweep 40 ppm
center -5 ppm
points 2048
noise 0.08

param pH = 7.05 [6.2, 7.6]
param Pi_ppm = 3.29 + 2.39 / (1 + pow(10, 6.75 - pH))

pool PCr  { shift 0        amp 1.00  width 10 Hz }
pool Pi   { shift Pi_ppm   amp 0.33  width 10 Hz }
pool bATP { shift -16.26   amp 0.85  width 10 Hz  couple 2 16.5 Hz }

process { lb 2 Hz  zerofill 4 }

report {
  "PCr/ATP" = area(PCr) / area(bATP)
  SNR       = snr(PCr)
}

The β-ATP triplet is the cleanest demonstration in the language that a “peak” and a “resonance” are different objects. One pool, one amplitude, one chemical shift — and three lines, in a 1:2:1 ratio, because the β phosphate is coupled to two others. Its integral is unchanged by the splitting; only its height falls.Every pool is a chemical species or, more precisely, a set of chemically equivalent nuclei. It has a shift, an amplitude, a width, and optionally a list of couplings, each of which splits every line the pool has so far into a first-order multiplet with binomial intensities. The amplitude is the integral of the whole multiplet and is conserved by splitting, which is the point.

Acquisition

These statements stand on their own at the top level, or may be gathered inside an acquire { … } block, which changes nothing but the tidiness.

Pools

A pool is a set of chemically equivalent nuclei — a metabolite, or one carbon of one. Everything inside is optional; an empty pool is a unit singlet at 0 ppm.

Processing

Everything in the process block happens after the signal exists and before it is drawn, which is the whole of the distinction between acquisition and post-processing.

Report

The distinction between area() and integral() is the point of the whole block. The first asks the model what it put there; the second measures the drawn spectrum between two ppm values, the way a person with a cursor does. Where peaks overlap the two disagree, and the size of the disagreement is the quantification bias you would have had.A report block is a list of named expressions, evaluated against the finished spectrum and printed under the FID. An optional string after the expression is a unit label.

report {
  pH        = 6.75 + log10((shiftmax(3.5, 6.5) - 3.29) / (5.68 - shiftmax(3.5, 6.5)))
  "PCr/ATP" = area(PCr) / area(bATP)
  linewidth = fwhm(PCr) "Hz"
  SNR       = height(PCr) / rms(-25, -20)
}

Expressions

Guards, shifts, amplitudes and report lines are all ordinary arithmetic over parameters and a few read-only quantities. Comparisons yield 1 and 0; and, or and not work on those.

What one part per million is worth

Gyromagnetic ratios in MHz/T, with the sign dropped; ¹²⁹Xe is negative, and the gas-phase value 11.777 is used rather than the bare-nucleus 11.860, because in vivo xenon work references itself to the gas peak in the same breath. Natural abundance is irrelevant for the hyperpolarised nuclei, where the polarisation is five orders of magnitude from thermal equilibrium and the isotope is enriched anyway.The chemical shift scale is defined in ppm precisely so that it does not depend on the field. The linewidth is not, and neither is the shim, and neither is anything else that happens in hertz. That mismatch is the whole of multinuclear spectroscopy's difficulty, and it is arithmetic:

Read the last two columns together. A 10 Hz line is 0.08 ppm on protons and 0.5 ppm on deuterium: the same shim, the same magnet, the same voxel, and a sixfold difference in what counts as resolved. ¹²⁹Xe has the opposite luck. Its gyromagnetic ratio is poor, so a dissolved-phase line 150 Hz wide is a perfectly ordinary 4 ppm — and the red cell and membrane resonances are 20.6 ppm apart, which is 730 Hz at 3 T and separable by a one-point Dixon with no fitting at all. The nucleus with almost the worst sensitivity has almost the best-resolved spectrum, because its electron cloud is enormous and polarisable and its chemical shift range runs to thousands of ppm.

Which is also why decouple 1H is not a convenience in ¹³C work but a precondition. Without it, the 30.2 ppm lipid methylene arrives as a triplet 125 Hz wide and every carbon in the spectrum overlaps its neighbour's multiplet.Carbon is the awkward case in the other direction. Its 200 ppm range looks generous until you notice that the one-bond ¹³C–¹H coupling is 125 to 170 Hz and does not shrink with field: at 3 T that splitting is 4 to 5 ppm wide, comparable to the entire separation between glutamate C4 and glutamine C4. Couplings are in hertz, shifts are in ppm, and going to a higher field is the only operation that improves the ratio between them.

Things worth trying

The examples are starting points rather than results. A few ways of pulling on them:

• Take the field away. Load the phosphorus example and drag B0 from 7 T down to 1.5 T. Nothing in the program changes except one number, and the phosphomonoester and phosphodiester pairs merge into the two anonymous humps that clinical 31P at 1.5 T actually delivers. Resolution in ppm is linear in field; nothing else in the model is.
• Turn off the decoupler. In the natural-abundance carbon example, drag dec to zero. Every carbon with attached protons explodes into a 125 Hz multiplet, the NOE evaporates, and a spectrum that was legible becomes a comb. Then put it back and watch the integrals change by more than the peaks do.
• Sweep the echo time. In the proton example, take TE from 30 ms to 288 ms. The lactate doublet flattens near 72 ms, inverts at 144, and comes back upright at 288, while everything else simply fades. The inversion is not relaxation; it is the coupling evolving through half a cycle, and the model shows the two effects separately because they are separate.
• Undersample on purpose. In the xenon lung example, reduce sweep until the enormous gas peak folds. Because the FID here is sampled rather than assembled in the frequency domain, aliasing happens for the right reason and lands in the right place, which is the cheapest way to convince yourself where a fold comes from.
• Look at what a huge peak does far away from itself. In the xenon lung example the gas resonance is a hundred times the two peaks that matter, and its line is given a Gaussian width because that is what shim imperfection produces. Change it to a pure Lorentzian of the same width and its wings alone — 200 ppm from the line that made them — lift the baseline under the membrane and red cell peaks. This is the same arithmetic that makes residual water a problem at 2 ppm in a proton spectrum.
• Truncate and apodise. In the lineshape example, set lb to zero and truncate to a few hundred points. The sinc wiggles at the foot of every peak are the Fourier transform of the rectangular window you just applied. Then raise lb until they vanish, and note what it cost in linewidth.
• Add a dead time. Set delay to 1.5 ms anywhere. The baseline rolls, the peaks go dispersive at the edges of the sweep, and phase1 set to 360 × delay × sweep degrees puts it all back. First-order phase is not a mystery; it is a time shift seen from the other side of the transform.
• Set the noise to zero and zero-fill. The peaks get smoother and the SNR in the report does not change. Zero filling interpolates; it is the one free operation in the chain, and it is free because it adds nothing.

What this does not model

If you need a basis set to fit real data with, this is not the tool: use a density-matrix simulation of the actual pulse sequence — GAMMA, FID-A, Vespa, or the vendor basis sets built with them. This page is for understanding and for teaching, and for sketching what an experiment might look like before you book the magnet.The model is first-order throughout, and it is worth being explicit about what falls outside it.

• Strong coupling. Multiplets here are binomial and symmetric. Real spin systems whose chemical shift difference is comparable to their coupling — citrate, glutamate, myo-inositol, most of the interesting proton metabolites — show roofing, extra lines, and intensities that no amount of arithmetic on couple will reproduce. Where the examples show those metabolites, the multiplets are drawn as explicit line lists chosen to look right, not derived.
• Localisation. There is no PRESS, no STEAM, no sLASER, no slice profile, and therefore no chemical shift displacement error, no outer volume contamination, and no J evolution during the pulses. TE enters as a single ideal weak-coupling factor.
• B₀ inhomogeneity. A real voxel contains a distribution of field offsets. Here that is represented by a Gaussian component in the lineshape, which is a decent approximation and not a model.
• Water suppression, eddy currents, sidebands, motion. None of them. The residual water in the proton example is a pool you declare, not something the simulation left behind.
• Concentration. Amplitudes are arbitrary units. Turning a peak area into millimolar requires relaxation correction, a reference signal, and a set of assumptions this page does not make for you.
• Hyperpolarised dynamics. The carbon and xenon examples show one spectrum from a series. Polarisation decay, flip-angle consumption and the kinetics that make those experiments worth doing are not simulated.

The chemical shifts and couplings in the examples are literature central values, gathered for the periodic chart for MR on this site and reused here. They are the middle of a reported range, not a measurement, and several of them move with pH, temperature, oxygenation and tissue.

How the interpreter works

Roughly twelve hundred lines of ordinary JavaScript in the source of this page, with no dependencies and no build step. View source and read from function lex.The pipeline is the textbook one. A lexer turns the source into tokens, keeping line and column so that errors can point at something, and emitting newlines as tokens, because one statement per line is what lets couple 2 16 Hz be a statement rather than a four-argument ambiguity. It also recognises nucleus names, which is the one place the grammar has to cope with an identifier that starts with a digit. A recursive-descent parser reads the statements and hands arithmetic to a precedence climber, producing a syntax tree.

The tree is then compiled to closures: every expression becomes a JavaScript function of one argument, the evaluation context. Identifier resolution happens once, at compile time, which is where an undefined parameter or a misspelled pool name is caught. Pool names inside area() and its relatives are resolved to indices at the same moment.

The simulator then does the only physics in the page. It evaluates the acquisition in dependency order, expands every pool into a list of lines by applying each couple in turn, and synthesises the complex FID one pool at a time by direct summation of decaying exponentials. Noise is added to the sum alone, so each component trace stays clean and the total does not. Apodisation, truncation, the first-point weight and zero filling are applied to all of them identically, then each is transformed by an iterative radix-2 FFT and phased. Keeping a separate spectrum per pool costs a little memory and buys the component traces, the measured linewidths, and an honest answer to the question of how much of that peak was actually the metabolite you named.