A sibling to the EEG simulator and the PET scanner simulator. Describe a recording in a few lines, record the magnetic field of current dipoles in the head with SQUIDs in a helmet or with OPMs on the scalp, let the elevator, the traffic, the heart and the eyes in, and then do what an MEG lab does: shield, compensate, project, separate, filter, average and localize. The simulator always knows the truth, so every answer can be checked.
When tens of thousands of pyramidal cells in a patch of cortex receive synaptic input at the same moment, the currents flowing inside their apical dendrites add up to a current dipole of a few tens of nanoampere-metres. Its magnetic field leaves the head undistorted by the skull and is about 100 femtotesla a few centimetres away: a billionth of the Earth's field, and a millionth of the field of a car passing in the street. Measuring it takes the most sensitive detectors there are, SQUIDs, coupled to pickup coils in a helmet; liquid helium in a dewar to keep them superconducting; a magnetically shielded room to keep the world out; and head position coils and a digitizer to know where the head is. Optically pumped magnetometers (OPMs) are a newer way to do the same thing at room temperature, on the scalp. Click a part of the system, or a stage below, to see how it works.
A fold of cortex. Pyramidal cells stand perpendicular to the cortical surface, so on the crown of a gyrus they point at the scalp and on the wall of a sulcus they lie along it. An excitatory synapse on the apical dendrite lets positive charge in; it flows down inside the cell (the primary current, blue arrows) and returns through the surrounding tissue (the volume current, grey). Seen from a few centimetres away the loop is a current dipole, and its magnetic field circles the primary current (dots: field out of the page, crosses: into it).
Magnetometer field map and EEG potential map of the active patch, from the same head model as the rest of the simulator. Synchronous cells add linearly; asynchronous ones add like a random walk, as the square root of their number. A radial patch on a gyral crown is nearly invisible to MEG; two opposed sulcal walls cancel for both methods.
A cross-section through a tangential dipole in the four-shell head: brain, cerebrospinal fluid, skull and scalp. Grey lines are the volume currents, computed for these conductivities; colour is the potential they set up. The skull conducts tens of times worse than the brain, so the current spreads sideways before crossing it and the scalp potential is smeared. In a spherical head the magnetic field outside does not depend on the conductivities at all: the volume currents' contribution is fixed by the primary current alone (Sarvas, 1987).
How each coil's response to a small magnetic source falls with distance, relative to a source 4 cm away. A magnetometer's falls as the cube of distance; a gradiometer subtracts two nearly equal readings of anything far away and its response falls one power faster. That is the gradiometer's distance filter: the brain is near, the world is far.
A DC SQUID is a superconducting ring broken by two Josephson junctions. Biased with a small current, the voltage across it oscillates with the magnetic flux through the ring, once per flux quantum, Φ0 = h/2e = 2.07 × 10-15 Wb. The electronics feed a current back into a coil beside the ring to hold the flux at the steepest point of the curve (the flux-locked loop); the feedback current is the output, linear over many flux quanta.
When the flux changes faster than the loop can follow, the SQUID passes through a whole period of its curve before the feedback catches up, and the loop locks again one flux quantum away. The output jumps by a fixed step and stays there: a flux jump or SQUID jump. Strong, fast fields (no shielded room, a nearby magnet) make them common.
SQUIDs and pickup coils are niobium, superconducting below 9 K, so they sit in liquid helium at 4.2 K inside a vacuum-insulated vessel shaped into a helmet at the bottom. A whole-head dewar holds about 70 litres and boils off several litres a day; many labs now recycle the gas. The price of the cold is distance: vacuum and insulation keep the coils about 2 cm from the scalp on top of the head and more where the head does not fill the helmet.
Layers of mu-metal, a nickel-iron alloy of very high permeability, give slowly varying fields an easier path around the room than through it. A layer of aluminium adds eddy-current shielding, which grows with frequency. Everything that must cross the walls does so through filtered connectors and waveguides; stimuli arrive by light, air tubes and plastic, and the participant is watched by camera. Opening the door, even an elevator next door, shows in the data.
Before the recording a 3D digitizer pen records three anatomical landmarks (the nasion and the two preauricular points), which define the head coordinate frame, the positions of four or five small head position indicator (HPI) coils taped to the scalp, and a hundred or more points on the scalp and the nose. In the helmet, each coil is driven with its own frequency, above the brain's band, and its field is fitted as a magnetic dipole: the coils' positions in the helmet give the head's position, continuously if the coils stay on.
Head position over this recording, as the HPI coils would track it (switch on Head movement on the Noise tab or in the program). With SQUIDs the helmet stays put and the head moves inside it, so every source's field map slides across the sensors.
Circularly polarized laser light pumps the spins of rubidium atoms in a small heated glass cell into one direction, and a pumped vapour stops absorbing the light. Any magnetic field makes the spins precess away and the cell absorbs more, so the light reaching the photodiode measures the field. The effect is sharpest near zero field: within about a nanotesla. That is why OPMs need a well-shielded room and coils that null the remaining field, and why moving them through a gradient of a few nanotesla per metre makes artifacts.
| SQUID magnetometer | OPM | |
|---|---|---|
| Operating temperature | 4 K (liquid helium) | about 150 °C inside the cell, skin-safe outside |
| Distance to the scalp | 2 to 4 cm, fixed helmet | about 6 mm, on the head |
| Noise | 2 to 5 fT/√Hz | 7 to 20 fT/√Hz |
| Dynamic range | very large, with flux-locked loop | a few nT; needs nulling coils |
| Bandwidth | DC to kHz | DC to about 150 Hz |
| Head movement | the head moves under fixed sensors | sensors move with the head, through any residual field |
| Fits | one adult head size | any head, children, movement |
Every MEG study can be written down as a handful of decisions: which sensors, how fast to sample, what the online filters do, which shielding and noise rejection to use, what the participant does and when, and what is in and around the head. The program on the left is exactly that, in a small language. Change a line and press Build recording: the simulator records a new session and every other tab updates. The examples are a good place to start; the language is described under the editor.
A program is a few blocks in braces. Inside a block each statement is a keyword and its values. A number takes a unit
(Hz, ms, s, fT, pT, nT, fT/cm, nAm,
cm, mm, deg, bpm) or a %, and the unit says what the number is, so
blinks 0.2 Hz 2 pT and blinks 2 pT 0.2 Hz mean the same. In a run of numbers the last unit applies to all:
epoch -100 400 ms. Line breaks are only spacing, and # starts a comment.
| Statement | Meaning |
|---|---|
| recording { ... } | |
| sensors magnetometers | axial | planar | opm | The array: 105 SQUID magnetometers, 105 axial gradiometers (5 cm baseline), 105 triplets of one magnetometer and two planar gradiometers (315 channels, as in the 306-channel systems), or 105 OPMs on the scalp. |
| rate 500 Hz | Sampling rate, 200 to 2000 Hz. |
| highpass 0.1 Hz | off lowpass 150 Hz | Online filters (first-order high-pass, fourth-order Butterworth low-pass). Default low-pass 0.3 × rate. |
| mains 60 Hz | Line frequency, 50 or 60 Hz. |
| noise 3 fT | Sensor noise density per √Hz (fT/cm/√Hz for planar gradiometers). Defaults: SQUIDs 3, axial 4, OPMs 15. |
| background 150 fT | off | RMS of spontaneous brain activity on the magnetometers (random cortical dipoles with 1/f time courses). |
| hpi on | Head position coils, each driven at its own frequency above the band of interest. |
| seed 1 | Random seed. The same program and seed always produce the same recording. |
| shielding { ... } | |
| room on | off | The magnetically shielded room. Off means a recording in an ordinary lab. |
| references on | Gradient compensation from the reference magnetometers high in the dewar (or on the wall, for OPMs). |
| ssp 3 | Signal space projection with this many empty-room vectors. |
| sss on | tsss | off order 8 3 | Signal space separation, or its temporal extension, with the inside and outside expansion orders. Needs sensors planar. |
| nulling on | For OPMs: coils around the head that cancel the room's residual field. |
| head { ... } | |
| position 0 0 -1 cm | Where the head sits in the helmet, relative to its centre (x right, y forward, z up). |
| source dipole | pair | auditory | somatosensory | alpha | spike name { ... } | |
| at C3 depth 2.5 cm | Place a dipole 2.5 cm under the scalp beneath a 10-10 position or a sensor (MLT31). Or pos x y z cm. |
| orient tangential | radial | inward | 40 deg dir 90 deg | Tilt from pointing straight out (0 radial, 90 tangential), and the direction of the tangential part: 0° toward the nose, 90° toward the right ear. toward C4 also works. |
| amp 20 nAm signal sine | burst | pink 10 Hz | dipole: moment and time course (a steady sinusoid, a waxing and waning rhythm, or 1/f noise). |
| dipole A { at C3 ... } rhythm 20 Hz coupling 0.6 lag 10 ms | pair: two dipoles sharing a fraction of a common rhythm, the second lagging the first. |
| side both | left | right amp 30 nAm latency 100 ms jitter 8 ms | auditory: the M50, M100 and M200 from the supratemporal plane; somatosensory: the N20m and later peaks from area 3b of the hand area. |
| peak 100 ms 35 ms 1 condition deviant 1.4 | Replace the default peaks (latency, full width at half maximum, gain), and set the gain for each condition. |
| beta 20 Hz 18 nAm erd 60% rebound 120% | somatosensory: a motor-cortex beta rhythm that drops after each stimulus and rebounds. |
| freq 10 Hz closed 100% open 20% erd 40% | alpha: two parieto-occipital generators; amplitude with the eyes closed and open, and the drop after each stimulus. |
| rate 0.2 Hz width 70 ms wave 50% | spike: interictal epileptiform spikes at random times, with an aftergoing slow wave. |
| timeline { ... } | |
| block closed | open 20 s | A rest block with the eyes closed or open. |
| trials 120 isi 1 s jitter 0.2 s { condition a 80% ... } | A run of stimuli, one every isi plus a random 0 to jitter, in random, alternate or sequence order. |
| artifacts { ... } | |
| environment 50 nT line 100 nT | Field changes from an elevator, traffic and the city, and the mains field, as they would be without the shielded room. |
| cardiac 70 bpm 3 pT blinks 0.2 Hz 2 pT saccades 0.15 Hz 1 pT muscle 0.05 Hz 2 pT | The heart, the eyes and the muscles, with their largest field on the magnetometers. |
| dental 400 pT | Magnetized dental work, moved by the jaw. |
| jumps 0.03 Hz 6 pT | SQUID flux jumps (ten times as frequent without the room). |
| bad MRT31 noisy | flat | jumpy | One misbehaving sensor. |
| movement 5 mm | Head movement in the helmet: drift, sinking, and sudden shifts. |
| … off | Declare an artifact but start with it switched off. Every artifact can be toggled on the Noise and Recording tabs. |
| analysis { ... } | |
| highpass 1 Hz lowpass 40 Hz notch 60 Hz order 4 phase zero | causal | Starting offline filters. |
| epoch -100 400 ms baseline -100 0 ms reject 4 pT 1500 fT/cm | Epoch window, baseline and peak-to-peak rejection thresholds (field for magnetometers, gradient for planar gradiometers). |
| ica 25 compensate on | Number of ICA components, and head movement compensation through SSS. |
The current inside the cells (the primary current) drives volume currents through the whole head. EEG measures the potential those volume currents set up on the scalp; MEG measures the magnetic field of both kinds of current. In a spherical head the field outside depends only on the primary current's position and its tangential part: a radial dipole makes no field outside at all, and a deep one makes very little, because the field falls with distance and the radial share grows toward the centre. The EEG potential map of a tangential dipole has its two poles on either side of the dipole, along it; the magnetic field circles the current, so the MEG map's two poles lie across it, rotated by 90 degrees. The skull smears the potential but not the field. Move the dipole below and watch both maps.
Direction is measured in the plane of the scalp above the dipole: 0° toward the nose, 90° toward the right ear.
Colour is the strength of the field in this plane on a logarithmic scale; lines follow its direction. The current flows out of the page at the dipole (dot). Outside the head the field is that of the whole head, primary and volume currents together.
Colour is the potential (red positive, blue negative); grey lines are the volume currents, leaving the head of the dipole and returning to its tail through the brain, the CSF, the skull and the scalp. EEG electrodes sit on the outer circle.
Largest signal across the array as the dipole turns from pointing out (0°) through tangential (90°) to pointing in (180°), each method relative to its own maximum. MEG falls to zero at radial; EEG barely changes.
Largest signal for the current orientation as the dipole moves along its radius toward the centre of the head, each relative to its value at the shallowest depth. Gradiometers and OPMs favour shallow sources most; EEG reaches deepest. At the centre of the sphere every orientation is radial and the MEG signal vanishes.
Sensitivity of one channel to a dipole anywhere on a sphere of cortex, in its best orientation at each point, on the same map projection as the field maps (the cross marks the sensor). A magnetometer sees a broad ring and is blind right under itself, where every current is either radial or makes a field along the surface; a planar gradiometer sees best what is right under it; an OPM on the scalp sees a tighter patch than a SQUID magnetometer at the same site, and more strongly. The EEG electrode, against the average reference, is most sensitive to radial currents right under it. The number is the area of cortex above half the channel's largest sensitivity.
The field of every source in the head reaches every sensor, so each channel records a weighted sum of all of them. The weights, one column per source and one row per sensor, form the mixing matrix A in x = A s. They depend on where each source is, its depth and its orientation, and on where the sensors are; for MEG, unlike EEG, not on the conductivities. Two sources whose field maps overlap are hard to tell apart, from the sensors and from any inverse method. Select a source to move it and watch its map, the mixing matrix and the recording change together.
The map is the source's column of the mixing matrix for the array in use, per 10 nAm of moment. A tangential source makes an outgoing and an ingoing pole on either side of itself, with the zero line right above it; the planar gradiometers instead peak right above it.
The top trace is what the sensor records, before noise rejection and without offline filters. Below it, each signal's share, all on the same scale; they add up to the top trace. Choose another sensor and the same sources appear with different weights.
Each column is one generator's field pattern, scaled to its largest value. Hover for values. Below, how alike the patterns are: the cosine of the angle between two columns. Two sources with similar patterns mix almost indistinguishably.
An evoked field of a few hundred femtotesla has to be measured next to the Earth's 50 microtesla, field changes of tens of nanotesla from moving steel and currents in the building, the mains, the heart's picotesla and the eyes and muscles. No single trick removes all of it. The shielded room takes out most of the world; gradiometers ignore what is far; reference sensors measure the remaining field and subtract it; signal space projection removes the directions in sensor space that noise occupies; signal space separation splits the field into a part from inside the helmet and a part from outside and keeps only the inside. Each step costs something. Switch them on one at a time and watch the budget below.
Fields at the head. On the left, typical values; on the right, the largest value of each signal in a 10 s stretch of this recording on the magnetometers (or the array's own sensors): what reaches the sensors (open circle) and what is left after the noise rejection chosen on the left (filled).
Each signal's own contribution, RMS over the main sensors in a 10 s stretch of the recording, before and after the noise rejection chosen above; all steps are linear, so the parts add up to the whole. The brain's signals are at the top: whatever removes the noise should leave them alone.
Traces are grouped by region from front to back, left before right. The same data can be shown as they reach the sensors, after the noise rejection chosen on the Noise tab, or after the offline filters and ICA as well. Switch signals on and off on the right to learn each one's signature: where on the helmet it appears, how big, how fast and how often. Click in the traces to put a cursor there and see the field map at that instant.
MEG spectra are usually drawn as an amplitude spectral density, in femtotesla per root hertz, because that is how sensor noise is specified: a good SQUID has a flat floor of a few fT/√Hz with a rise at the lowest frequencies. Brain activity falls roughly as 1/f above that floor and rhythms stand out as peaks: alpha near 10 Hz over the back of the head with the eyes closed, beta near 20 Hz over the motor strip. Line noise is a needle at 50 or 60 Hz and its harmonics; the environment and the heart crowd the lowest frequencies. A time-frequency map shows when power changes: alpha appearing as the eyes close, or beta dropping after a stimulus and rebounding.
Noise-rejected data. Power is integrated over the band from each channel's spectrum and interpolated between sensor sites.
After the noise rejection come the steps an analyst chooses. Filters remove what you do not want and distort what you do, more so the steeper and closer to the signal they are, and a filter that does not shift latencies must look into the future. ICA can separate the heart, the eyes and a bad sensor from the brain only to the extent that they have fixed field patterns and independent time courses. Head movement blurs every average; with an array that supports SSS the field can be rebuilt as if the head had stayed still. The Recording, Spectra and Averaging tabs all show the output of this pipeline.
ICA assumes the recording is x = A s: a fixed mixture of signals that are statistically independent. It finds an unmixing matrix that makes the recovered signals as non-Gaussian as possible, then shows each component's field pattern (a column of A) and its time course. Channels are divided by their noise first so magnetometers and gradiometers can be mixed. Mark components to remove and their back-projection is subtracted from the data. Here you can check each decision against the truth.
ICA has not been run on this recording.
Select a component to see its time course and spectrum.
The head's position over the recording, from the HPI coils: displacement of the centre of the head from where it started, in millimetres, along each axis of the helmet.
The evoked field's map, averaged over the trials at the head positions where they happened: as if the head had stayed still, as recorded, and after SSS rebuilt each trial's inside field for the starting position. Compensation needs an array that supports SSS (planar gradiometers with magnetometers). OPMs move with the head, so for them the brain does not move at all.
An evoked field of a hundred femtotesla is buried in hundreds of femtotesla of everything else. If the response is the same on every trial and the rest is not time-locked to the stimulus, averaging N trials leaves the response untouched and shrinks the rest by a factor of √N. Signal-to-noise ratio grows as √N, so doubling it costs four times the trials. The catch is in the assumptions: if the latency wanders from trial to trial the average is smaller and wider than any single response, and artifacts that are time-locked or huge survive. Epochs are cut around each event, baseline corrected, screened for artifacts, and then averaged.
Noise is the RMS of the average in the baseline (what you can measure) and the RMS of the average minus the true evoked field over the whole epoch (what only a simulation can measure). Signal is the true field's largest absolute value after the stimulus. On a log-log plot √N is a straight line of slope one half.
One row per epoch of the selected sensor, colour for field. A single trial is mostly noise, so each row is usually averaged with its neighbours, as in the ERP images of EEG and MEG papers; rejected epochs are left raw and marked on the left. Sorting by the true latency jitter, which the simulator knows and an experimenter does not, shows how the average blurs a response that wanders in time.
The forward problem, from a known source to the field at the sensors, has exactly one answer, and physics supplies it. The inverse problem, from the field back to the sources, has infinitely many: Helmholtz showed in 1853 that infinitely many current distributions inside a conductor produce the same field outside. Every inverse method adds an assumption to pick one. Dipole fitting assumes one or a few point sources. Minimum norm spreads the current over the whole brain and picks the distribution with the least total power. A beamformer assumes the sources are uncorrelated in time. Before any of that the head must be put in the same frame as an MRI, because a source is only as well placed on the anatomy as the coregistration that put it there.
Loading the template MRI; a synthetic head stands in until it arrives. The three landmarks were digitized on the participant before the recording; marking the same three points on the MRI gives the transform between the two frames. The dots are the digitized head shape placed with that transform: where they float off the scalp or sink into it, the transform is wrong. The preauricular points lie just in front of each ear and are the hardest to mark the same way twice. Every source estimate below is drawn on the anatomy through this transform, so its error moves them all.
Press Estimate to fit, minimum-norm and beamform the data above. The first time takes a few seconds, to compute the field of every point of the source grid; after that the estimates follow the controls.
The tabs follow the life of an MEG study. The Instrument tab is the hardware, from a patch of cortex to a SQUID in liquid helium. The Experiment tab is the study, written as a short program; everything else is computed from it. MEG & EEG fields is the physics of one dipole, and Sources the sources of this recording. Noise & shielding is what MEG labs spend most of their effort on. Recording, Spectra, Processing and Averaging are the analysis; Source estimation runs the problem backwards. The concept boxes at the top of each tab open and close, so a projected screen can stay uncluttered.
tools/build_template.py; the source grid is its grey and white matter, the
digitizer and ICP use its scalp, and the drawings of the head follow its outline. Sources and licences are in
data/README.md.Built for NBL 425/625, Methods in Human Neuroimaging. See also the EEG simulator and the course's EEG versus MEG and wearable OPM-MEG demos.