# Requires Python 3.10 or higher
pip install MRzeroCore

MR-zero is a framework for easy MRI sequence optimization and development of self-learning sequence development strategies. The vision is documented in this paper. These goals are backed by a modern Bloch simulation (see this paper). More material can be found in the literature.

Quick start

For a quick introduction, look at the MR-zero Playground!

An ever growing collection of Jupyter Notebook demonstrate core capabilities of MR-zero. These notebooks run on Google Colab - directly in your browser, no installation needed!

Alternatively, install MR-zero Core locally:

# Requires Python 3.10 or higher
pip install MRzeroCore

Simulate any Pulseq sequence in one line

MR-zero makes MRI sequence simulation easy - just one line for simulating .seq files:

import MRzeroCore as mr0

# That's it - automatic phantom download and simulation!
signal, ktraj_adc = mr0.util.simulate('your_sequence.seq')

Even simpler with PyPulseq - no need to worry about writing .seq files:

import pypulseq as pp
# Create sequence with PyPulseq
seq = build_your_sequence()
# ... build sequence ...
signal, ktraj_adc = mr0.util.simulate(seq)
# ... reconstruct, e.g.: with a NUFFT of signal and ktraj_adc

Further documentation

MR-zero uses PyTorch for computation, enabling GPU compute automatic differentiation with backpropagation. This means you can easily develop your own loss functions, sequence building code and more - and then optimize any input parameters efficiently with gradient descent. In the example above, this could mean to simply extend build_your_sequence() by sequence parameters like flip angles and writing an image based loss function for the reconstruction. Then you can lie back and let the computer find the best pulse train.

How this can be done will be explained in the following pages, where this documentation lists the most important ideas and applications of MR-zero. If you think that something is missing or misleading, pleas open an issue on GitHub or directly via the button on the top right on every documentation page.

Playground

Welcome to Playground MR0, a playground to share, vary and simulate MR sequences. MR sequences are written in the Pulseq standard using the pypulseq library. Pulseq files are simulated with the efficient Phase Distribution Graph Bloch simulation. Here we share links to example colabs that contain various MR sequences or let you upload your own seq file for simulation.

Many of the examples are build using PyPulseq and simulate the resulting .seq files with MR0. These .seq files could also be measured on any MRI scanner using a Pulseq interpreter.

Code and simulate PyPulseq

Notebook
Free Induction Decay	mr0_FID_seq.ipynb
Spin Echo CPMG	mr0_SE_CPMG_seq.ipynb
Stimulated Echo 3 pulses - 5 echoes	mr0_STE_3pulses_5echoes_seq.ipynb
FLASH 2D sequence	mr0_FLASH_2D_seq.ipynb
GRE EPI 2D sequence	mr0_EPI_2D_seq.ipynb
DWI SE EPI 2D sequence	mr0_DWI_SE_EPI.ipynb
Diffusion prepared STEAM	mr0_diffusion_prep_STEAM_2D_seq.ipynb
RARE 2D sequence	mr0_RARE_2D_seq.ipynb
TSE 2D sequence	mr0_TSE_2D_multi_shot_seq.ipynb
Interactive GRE to FLASH	mr0_GRE_to_FLASH.ipynb
balanced SSFP sequence	mr0_bSSFP_2D_seq.ipynb
DREAM STE for B0, B1, TxRx mapping	mr0_DREAM_STE_seq.ipynb
DREAM STID for B0, B1, TxRx mapping	mr0_DREAM_STID_seq.ipynb
Pulseq with RF shimming	pulseq_rf_shim.ipynb
Compare runtimes and accuracy settings for MR0 simulation	mr0_time_vs_acc_bspBSSFP.ipynb

Plot and simulate predifined .seq files

Notebook
Simulate own uploaded seq files	mr0_upload_seq.ipynb

MR-zero optimization

Gradient descent optimizations using automatic differentiation by backpropagation. Some notebooks use pulseq-zero for optimizable sequence definitions with PyPulseq.

Notebook
IR FLASH 2D sequence for T1 mapping using a fit	mr0_opt_FLASH_2D_IR_Fit_T1.ipynb
IR FLASH 2D sequence for T1 mapping using a NN	mr0_opt_FLASH_2D_IR_voxelNN_T1.ipynb
FLASH flip angle opt. for PSF (with pulseq-zero)	Pulseq_zero_FLASH_FAopt_PSFtask.ipynb
TSE flip angle opt. for SAR (with pulseq-zero)	Pulseq_zero_TSE_FAopt_SARtask.ipynb
DESC with pulseq-zero	pulseq_zero_DESC_demo.ipynb

MR-double-zero optimization

Gradient-free optimization with nevergrad

Notebook
Ernst angle optimization	mr00_FLASH_2D_ernstAngle_opt.ipynb

MR plot wall of fame

famous historic plots recreated

MR0 example notebooks

The following sequences are examples of how to realize various tasks in MR-zero rather than demonstrations of specific MRI sequences.

Notebook
Pure MR0 FLASH	flash.ipynb
pulseq FLASH	pulseq_flash.ipynb
pulseq pTx FLASH	pulseq_sim_pTx.ipynb

Integration

MRI simulations are usually not used in isolation. Instead, simulated sequences are written with various tools, optimized flip angle trains are measured, quantified phantoms are viewed in different programs. The following chapters list currently existing integrations with MR-zero and show how to utilize them.

Pulseq Integration

MRzero Core provides seamless integration with the Pulseq standard, making it easy to simulate any Pulseq sequence file in just one line of code.

Get Started Immediately

Simplest possible - one line simulation:

import MRzeroCore as mr0

# That's it - simulate any .seq file with automatic phantom download!
signal, ktraj_adc = mr0.util.simulate(seq)  # Downloads standard phantom if not found

seq.plot(plot_now=False)
mr0.util.insert_signal_plot(seq=seq, signal =signal.numpy())  # show sequence plots with simulated signals inserted in ADC plot.
plt.show()

Full control with prepass-mainpass:

import MRzeroCore as mr0
obj_p = mr0.VoxelGridPhantom.brainweb('subject05.npz') #load a phantom object from file
obj_p = obj_p.interpolate(64, 64, 32).slices([15])
obj_p = obj_p.build()

seq0 = mr0.Sequence.import_file('your_sequence.se')  # Read in the Pulseq seq file

graph = mr0.compute_graph(seq0, obj_p, 2000, 1e-4)  # pre-pass
signal = mr0.execute_graph(graph, seq0, obj_p, print_progress=False) # main-pass

seq.plot(plot_now=False)
mr0.util.insert_signal_plot(seq=seq, signal =signal.numpy())
plt.show()

Try now: mr0_upload_seq.ipynb

Simplified Workflows

Utily functions

MRzero now includes streamlined utility functions for common simulation tasks:

1. One line - automatic everything

   # Quickest 2D brain phantom simulation 
   signal, ktraj_adc = mr0.util.simulate(seq)  # Downloads standard phantom automatically
   

2. Two lines - custom phantom parameters

   # Quick 2D brain phantom sim with modifications
   obj = mr0.util.load_phantom([128,128], B0_polynomial=(0,0,-150,1,1,0))  # Custom B0 map
   signal, ktraj_adc = mr0.util.simulate(seq, obj)
   

3. Three lines - full control

   # Custom phantom from URL/file with arbitrary modifications and high accuracy
   obj = mr0.util.load_phantom([128,128], B0_polynomial=(0,0,-150,1,1,0), url='numerical_brain_cropped.mat')  
   obj.B1 = obj.B1**3     # Very inhomogeneous B1+ field
   signal, ktraj_adc = mr0.util.simulate(seq, obj, accuracy=1e-5)
   

Direct PyPulseq Integration

import pypulseq as pp

# Create sequence with PyPulseq
seq = pp.Sequence()
# ... build sequence ...

# Simulate directly - no need to write/read .seq file!
signal, ktraj_adc = mr0.util.simulate(seq)  # seq can be PyPulseq object

For MATLAB Pulseq Users

See MATLAB integration

Examples - Ready-to-run in Google Colab

Notebook	Complexity
Simulate own uploaded seq files	Beginner	mr0_upload_seq.ipynb
FLASH 2D sequence	Beginner	mr0_FLASH_2D_seq.ipynb
GRE EPI 2D sequence	Intermediate	mr0_EPI_2D_seq.ipynb
DWI SE EPI 2D sequence	Advanced	mr0_DWI_SE_EPI.ipynb
TSE 2D sequence	Advanced	mr0_TSE_2D_multi_shot_seq.ipynb

Browse All 20+ Examples →

Further Features

Visualization

# Plot PyPulseq sequence plot with MR-zero signals
seq.plot(plot_now=False)
mr0.util.insert_signal_plot(seq=seq, signal =signal.numpy())
plt.show()

# Plot the detected k-space used by MR-zero
seq.plot_kspace_trajectory()

#Plot the phase distribution graph
graph = mr0.compute_graph(seq0, obj_p, 2000, 1e-4)  # pre-pass
graph.plot()

Simulation

# Simple simulation
signal = mr0.util.simulate(seq)

# With custom phantom
phantom = mr0.util.load_phantom(size=(64, 64))
signal, kspace = mr0.util.simulate(seq, phantom)

Getting Help

Quick help:

• Common issues: FAQ & Troubleshooting
• GitHub Issues: Report bugs or ask for features
• Github Discussions: Ask Questions

MATLAB Pulseq ↔ MR-zero Integration Guide

This guide shows how to seamlessly integrate MATLAB Pulseq with Python MR-zero for powerful sequence development and simulation workflows.

Quick Start for MATLAB Users

Already have MATLAB Pulseq sequences? Get started in 3 steps:

1. Export your sequence: seq.write('my_sequence.seq') in MATLAB
2. Upload to Colab: Use our (mr0_upload_seq.ipynb )
3. Simulate instantly: No Python installation required!

Integration Workflows

Workflow 1: MATLAB seq creation → Python Simulation

Use Case: Leverage MR-zero’s fast simulation while keeping your MATLAB development workflow

% In MATLAB: Create your sequence as usual
seq = mr.Sequence();
% ... build your sequence ...
seq.write('my_flash.seq');

Then switch to python and load the seq file.

# In Python/Colab: Simulate with one line
import MRzeroCore as mr0
seq = mr0.Sequence.import_file('my_flash.seq')
signal = mr0.util.simulate(seq)

You can just upload your seq file in mr0_upload_seq.ipynb

Workflow 2: MATLAB seq creation → Call Python from Matlab using sys command

see Pulseq demo simMR0

Need Help? Try our interactive examples first, or check the GitHub issues for community support.

Pulseq-zero

Frequently asked questions

Common questions and solutions for MRzero Core and Pulseq integration.

What’s the easiest way to simulate a sequence?

A: Code and export your sequence as Pulseq file. One line:

import MRzeroCore as mr0

# Automatic phantom download and simulation
signal, ktraj_adc = mr0.util.simulate('your_sequence.seq')

For complete Pulseq integration features, see: Pulseq Integration Guide →

Can I simulate PyPulseq sequences directly without writing .seq files?

A: Yes! This is new functionality:

import pypulseq as pp
import MRzeroCore as mr0

# Create sequence with PyPulseq
seq = pp.Sequence()
# ... build sequence ...

# Simulate directly - no .seq file needed!
signal, ktraj_adc = mr0.util.simulate(seq)  # seq is PyPulseq object

Note, this is will still generate a temporary seq file. For avoiding seq files completely, check out Pulseq-zero

How do I use custom phantoms with B0 field modifications?

A:

# Load phantom with polynomial B0 augmentation
obj = mr0.util.load_phantom(
    [128, 128], 
    B0_polynomial=(0, 0, -150, 1, 1, 0),  # (const, lin_x, lin_y, quad_x, quad_y, quad_xy)
    url='numerical_brain_cropped.mat'     # Custom phantom URL or local file
)

# Further customize
obj.B1 = obj.B1**3     # Inhomogeneous B1+ field
obj.T2 = 80e-3         # Custom relaxation times

# High-accuracy simulation
signal, ktraj_adc = mr0.util.simulate(seq, obj, accuracy=1e-5)

Does MRzero support ISMRMRD raw data format?

A: Current status: ISMRMRD raw data support is not yet available, but development is underway as a first pull request exists. So it can be expected in upcoming releases. Stay updated:

• Watch the GitHub repository for ISMRMRD support announcements
• Check release notes for new features

My .seq file won’t load / import fails

A: Pulseq 1.4.2 is used and tested by us the most. Pulseq 1.3 and 1.5 should work too, but might need some different load and plot functions.

See also in the comments in mr0_upload_seq.ipynb

Getting version compatibility errors

A: Most issues are resolved by:

• Using matching Pulseq versions between creation and simulation
• Updating to latest MRzeroCore: pip install --upgrade MRzeroCore
• Checking GitHub Issues for known compatibility problems

Simulation is too slow / runs out of memory

A: Reduce sequence resolution, go to single slice, reduce phantom resolution, reduce accuracy.

# Reduce phantom size for testing
phantom = mr0.util.load_phantom(size=(32, 32))  # Instead of (128, 128)

# Reduce simulation accuracy for speed
signal = mr0.util.simulate(seq, accuracy=1e-2)  # Default: 1e-3

# Use smaller sequence for testing
# Remove unnecessary repetitions during development

GPU not being used / CUDA errors

A:

import torch
print(f"CUDA available: {torch.cuda.is_available()}")
print(f"CUDA device: {torch.cuda.get_device_name(0) if torch.cuda.is_available() else 'None'}")

# MR-zero uses GPU when seq0 and obj is using cuda
```python
graph = mr0.compute_graph(seq0, obj, 200, 1e-4)
# Minimum emitted and latent signal of 0.01 (compared per repetition to strongest state)
signal = mr0.execute_graph(graph, seq0.cuda(), obj.cuda(), 0.01, 0.01)

How does MRzero handle slice-selective pulses?

A: Current limitation: MRzero uses instantaneous pulses and currently ignores slice selection gradients. This means:

# In PyPulseq - slice selective pulse
rf, gz, gzr = pp.make_sinc_pulse(
    flip_angle=90*np.pi/180, 
    slice_thickness=5e-3,  # This is ignored in MRzero
    duration=1e-3
)

To have correct gradient behavior (spoiling, diffusion), our interpreter puts the instantanoues pulse in the center of the pulse, and correctly plays out the gradients. Thus, if you have assymetric pulses, you have to be careful, and better replace them by block pulses.

Workarounds:

1. For 2D simulations: Use thin phantom slices to approximate slice selection
2. For comparison with Scans: Consider this limitation when comparing to real measurements, apparent flip angles might be smaller due to slice profile effects.

How do I choose the right accuracy parameter?

A: The accuracy parameter controls simulation precision vs speed. When in doubt rerun with higher accuracy.

# High accuracy (slow, precise)
signal = mr0.util.simulate(seq, accuracy=1e-8)  # Research-quality

# Fast accuracy (quick testing)
signal = mr0.util.simulate(seq, accuracy=1e-2)  # Development/testing

# Very fast (rough estimates)
signal = mr0.util.simulate(seq, accuracy=1e-1)  # Initial debugging only

What does the accuracy parameter actually control?

A: It sets the min_emitted_signal and min_latent_signal thresholds in the PDG simulation:

• Lower values = more computation but higher precision
• Higher values = faster simulation but potential artifacts
• For sequences with small signals (like multi-echo), use lower values

Does MRzero support magnetization transfer effects?

A: Current status: MRzero currently does not model MT effects. The simulation assumes:

• Single pool (free water) only
• No bound pool interactions
• No MT saturation effects

How accurate is diffusion simulation in MRzero?

A: MRzero handles currently only isotropic diffusion, using the EPG equations from Weigel et al.:

How are relaxation times handled?

A:

# Set realistic values in phantom
phantom.T1 = 1.0    # seconds
phantom.T2 = 80e-3  # seconds  
phantom.T2dash = 50e-3  # T2' dephasing effects with 1/T2* = 1/T2 + 1/T2'

# All relaxation effects are included in simulation
# including T2* decay during readout

Tips:

• Use literature values for tissue types
• T2’ does NOT include intravoxel B0 inhomogeneity effects, this can be achieved by high enough B0 map resolution.

Sequence created in MATLAB doesn’t work in Python

A: Check Pulseq version compatibility

# Try different import approaches
seq = mr0.Sequence.import_file('file.seq')  # Standard since Pulseq 1.4, shoudl work for 1.5
seq = mr0.Sequence.from_seq_file('file.seq')  # Alternative for Pulseq 1.3 

Resources:

• GitHub Issues: Report bugs
• Examples: Interactive examples
• API Docs: Complete reference
• Pulseq Integration: Integration guide

Still stuck? Try the interactive examples first - they solve most common issues!

NIfTI phantoms

In previous versions, MR-zero used ad-hoc formats for its simulation phantoms. It stored the right combination of maps in numpys compressed file format. Starting 2026, we fully switch to a new file format with a precise specification. It is a combination of the widely used NIfTI file format an a JSON configuration file. Exact definitions of the contents of this configuration file makes writing loaders for different applications / programming languages easy. At the same time, care was taken to allow great flexibility and human readability of the files.

You can see the specification here. A closer look at the coordinate system is shown here.

To see how to use phantoms, you can either look at examples or view the API documentation.

BrainWeb NIfTI phantoms

Phantoms are built from BrainWeb data and stored in the NIfTI format. This data is not included directly. Instead, a BrainWeb downloader is provided by mr0, which automatically downloads the necessary data and combines it with literature values for the physical properties at different field strenghts. It is provided by `mr0.generate_brainweb_phantoms:

import MRzeroCore as mr0

# Stores the downloaded data in a `cache` folder in the given directory.
# All generated phantoms are emitted as a series of folders.
mr0.generate_brainweb_phantoms("output/brainweb")

These phantoms can be loaded as tissue dictionaries, where individual tissues can be modified before converting it into simulation data. For more information look at the NIfTI API documentation.

tissues = mr0.TissueDict.load("output/brainweb/subj05/subj05-3T.json")
# you can modify individual tissues - its a dictionary of `VoxelGridPhantom`s:
tissues.pop("fat")
tissues["gm"].T1 *= 1.2

# either merge all tissues into one VoxelGridPhantom with `tissues.combine()`,
# which can be useful for plotting the maps or slightly faster simulation
tissues.combine().plot()
# ...or build simulation data directly to get partial volume effects:
sim_data = tissues.build()

Custom phantoms for experimentation

Sometimes it is useful to generate phantoms consisting of individual voxels. This can help with educational tasks, analysing effective resolution and PSFs and more.

For a single voxel:

obj_p = mr0.CustomVoxelPhantom(
    pos=[[-0.25, -0.25, 0]],
    PD=[1.0],
    T1=[3.0],
    T2=[0.5],
    T2dash=[30e-3],
    D=[0.0],
    B0=0,
    voxel_size=0.1,
    voxel_shape="box"  # Can also be "sinc" or "gauss"
)

For a more complex structure try an L-shape:

obj_p = mr0.CustomVoxelPhantom(
    pos=[[-0.4, -0.4, 0], [-0.4, -0.2, 0],
         [-0.3, -0.2, 0], [-0.2, -0.2, 0], [-0.1, -0.2, 0]],
    PD=[1.0, 1.0, 0.5, 0.5, 0.5],
    T1=[1.0, 0.5, 0.5, 0.5, 2],
    T2=0.1,
    T2dash=0.1,
    D=0.0,
    B0=0,
    voxel_size=0.1,
    voxel_shape="box"
)

Predefined phantoms

Older versions of MR-zero Core stored phantoms as MATLAB .mat files or as NumPy .npz. These are deprecated in favour of the new NIfTI phantom standard, but still available online.

A simple 2D brain phantom as MATLAB .mat file

This phantom comes from a quantification measurement. It is 2D, 1Tx (single channel B1+) and contains the following data:

• PD (proton density)
• T1 relaxation
• T2 relaxation
• T2’ dephasing
• ADC (apparent diffusion coefficient)
• dB0 inhomogeneities
• B1+ transmit field
• B1- recieve field

The phantom is stored on GitHub: numerical_brain_cropped.mat. Load it with MR-zero:

# Load file with a shell command - works in interactive notebooks (e.g.: on Google Colab):
!wget https://github.com/MRsources/MRzero-Core/raw/main/documentation/playground_mr0/numerical_brain_cropped.mat &> /dev/null

# After downloading, load with MR-zero
import MRzeroCore as mr0
obj_p = mr0.VoxelGridPhantom.load_mat('numerical_brain_cropped.mat')
# You can modify the phantom first:
# change resolution, remove dB0 inhomogeneities, increase diffusion...
brain_phantom_res = 64
obj_p = obj_p.interpolate(brain_phantom_res, brain_phantom_res, 1)
obj_p.B0[:] = 0
obj_p.D *= 1.5
# Convert the phantom into sparse SimData
obj_p = obj_p.build()

Predefined 3T BrainWeb phantom

The old, pre-generated phantoms from BrainWeb data are replaced by the new NIfTI phantoms. There is still one available for backwards compatibility, stored on GitHub: subject05.npz It contains all maps except for recieve fields:

• PD (proton density)
• T1 relaxation
• T2 relaxation
• T2’ dephasing
• ADC (apparent diffusion coefficient)
• dB0 inhomogeneities (generated empirically)
• B1+ transmit field (generated empirically)
• B1- recieve field

# Load file with a shell command - works in interactive notebooks (e.g.: on Google Colab):
!wget https://github.com/MRsources/MRzero-Core/raw/main/documentation/playground_mr0/subject05.npz &> /dev/null

sz = [64, 64]
# (i) load a phantom object from file
obj_p = mr0.VoxelGridPhantom.brainweb('subject05.npz')
obj_p = obj_p.interpolate(sz[0], sz[1], 32).slices([15])
obj_p.size[2]=0.08

# ========================================================
# EXAMPLE: insert rectangular "Tumor" in the diffusion map
# ========================================================

# typical brain tumor ADC values are around ~1.5 * 10^-3 mm^2/s,
# which lies between GM/WM and CSF (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3000221)
# mr0 uses D in units of 10^-3 * mm ^2/s  this is the same as µm^2/ms

# construct tumor border region
for ii in range(15, 25):
    for jj in range(15, 25):
        obj_p.D[ii, jj] = torch.tensor(0.75)

# construct tumor filling
for ii in range(16, 24):
    for jj in range(16, 24):
        obj_p.D[ii, jj] = torch.tensor(1.5)

# Convert the phantom into sparse SimData
obj_p = obj_p.build()

NIfTI Phantom Specification

The “NIfTI phantom specification” describes the storage of data suitable for MR imaging simulations. A phantom consists of one or more NIfTI files for per-voxel data and a single JSON file that defines the phantom and references the NIfTI files.

1. Storage

The naming of files is a convention: implementations can but are not required to reject non-conforming names. The supported convention is:

📂 subj42
├ 📄 subj42.nii.gz
├ 📄 subj42_dB0.nii.gz
├ 📄 subj42_B1+.nii.gz
├ ...
├ 📄 subj42-3T.json
└ 📄 subj42-7T.json

The phantom name is subj42, used as directory name and file prefix. It is available in two variants: subj42-3T.json and subj42-7T.json. Contents of the JSON file are described in Section 4.

Per-voxel data is stored in <name>_<property>.nii(.gz) files; the density map omits the property postfix. See Section 2 for properties and Section 3 for the NIfTI format requirements. The density property is required for every tissue loaded from NIfTI (as opposed to a constant value).

2. Phantom Definition

Units

The JSON file contains a "units" field that explicitly specifies the unit of each property. Currently, only the default units listed below are supported; automatic conversion may be added in the future.

System

The JSON file contains a "system" field that describes the physical system of the MRI experiment.

The gyro value implicitly defines the nucleus: 42.5764 MHz/T corresponds to ^1H. Tissue properties (particularly relaxation times and dB0) are specific to the system they were measured at. A phantom defined at 3T should not be used in a 7T simulation without adjusting tissue values accordingly.

Tissue Properties

The following properties are defined by the NIfTI phantom specification. If stored as per-voxel NIfTI data, the file postfix should match the key (e.g.: subj42_dB0.nii.gz).

3. NIfTI Data

Per-voxel tissue properties are stored in .nii files following the NIfTI v1.1 specification, optionally gzip-compressed (.nii.gz).

• Each file contains a single property for all tissues
• Data must be 4-dimensional (use singleton dimensions for non-3D data)
  • Dimensions 1-3: spatial (size 1 if unused)
  • Dimension 4: tissue index
• All NIfTI files must share the same resolution and orientation
• Spatial data should follow the RAS+ convention (index 0: R, 1: A, 2: S, growing towards positive) to ensure correct orientation for tools ignoring the affine matrix
• The affine matrix must transform data into RAS+ using mm as units (as per NIfTI spec)

4. JSON Phantom File

The phantom is defined in a .json file with the following structure:

{
  "$schema": "https://github.com/mrx-org/bifti-phantoms/bifti-phantom-v1",
  "units": {
    "gyro": "MHz/T",
    "B0": "T",
    "T1": "s",
    "T2": "s",
    "T2'": "s",
    "ADC": "10^-3 mm^2/s",
    "dB0": "Hz",
    "B1+": "rel",
    "B1-": "rel"
  },
  "system": {
    "gyro": 42.5764,
    "B0": 3.0
  },
  "tissues": {
    "gm": {
      "density": "subj42.nii.gz[0]",
      "T1": 1.56,
      "T2": 0.083,
      "T2'": 0.32,
      "ADC": 0.83,
      "dB0": "subj42_dB0.nii.gz[0]",
      "B1+": [
        "subj42_B1+.nii.gz[0]",
        "subj42_B1+.nii.gz[1]",
        "subj42_B1+.nii.gz[2]",
        "subj42_B1+.nii.gz[3]",
        "subj42_B1+.nii.gz[4]",
        "subj42_B1+.nii.gz[5]",
        "subj42_B1+.nii.gz[6]",
        "subj42_B1+.nii.gz[7]"
      ]
    },
    "wm": {
      "density": "subj42.nii.gz[1]",
      "T1": 0.83,
      "dB0": "subj42_dB0.nii.gz[0]"
    },
    "fat": {
      "density": "subj42.nii.gz[4]",
      "dB0": {
        "file": "subj42_dB0.nii.gz[0]",
        "func": "x - 420"
      }
    }
  }
}

• $schema: identifies the format, must end in nifti-phantom-v1 or bifti-phantom-v1 (optionally followed by a file extension).
• units: units for all properties (see Units).
• system: physical system parameters (see System).
• tissues: keys are arbitrary tissue names; values define tissue properties.

Each tissue must have a density of type file_ref (it defines the tissue’s shape). All other properties are optional and fall back to their default values. B1+ and B1- are lists of maps to support multi-channel data (a single-element list is valid).

Tissue Property Values

Each property in a tissue can be set via one of three methods:

file_ref

A file reference has the form <file_name>[<index>]. The file must be a .nii(.gz) in the same directory as the JSON file. The index selects a volume along the 4th dimension of the NIfTI file, allowing multiple tissues’ data to share a single file.

mapping_func

A mapping function is an equation evaluated per voxel. It supports basic arithmetic (+ - * / ( )) and the following variables:

Examples:

Coordinate System

• We always use RAS+ in a subject-aligned coordinate system
• NIfTI’s can store two orientations at once and do not specify which one to use
• MITK uses a LPS+ coordinate system and negates the xy affine entries on loading
• The scanner says data is in the SCANNER coordinate system, but this changes with sequence settings.
• Phantom z direction should always point in $B_0$ direction

Note

We assume that in measurement and FOV, MRI sequences are aligned to the subject.

When storing phantoms, we should always orient them to the subject-aligned RAS+ system (origin best at center of FOV but can be arbitrary). Correctly stored with sform_code == 2 and qform unused (qform_code == 0), which is the default for nibabel but check for simpleITK!

MR-zero

The MR-zero coordinate system is equivalent to NIfTI’s RAS+ with voxel coordinates at their centers. MR-zero (currently) ignores the affine matrix, except to extract the physical size of the NIfTI data. It generates the voxel coordinates with the following code:

pos_x, pos_y, pos_z = torch.meshgrid(
    size[0] * torch.fft.fftshift(torch.fft.fftfreq(shape[0])),
    size[1] * torch.fft.fftshift(torch.fft.fftfreq(shape[1])),
    size[2] * torch.fft.fftshift(torch.fft.fftfreq(shape[2])),
    indexing="ij"
)

There is no corner vs center here - voxels are points. The grid follows the typical FFT definition torch.fft.fftshift, which means they are placed as follows:

# Even number of samples, e.g.: N = 32
pos = [-16, -15, ..., -1, 0, 1, ..., 15] / 32 * size
# Odd number of samples, e.g.: N = 21
pos = [-10, ..., -1, 0, 1, ..., 10] / 21 * size

This is matching typical FFT impls and k-space encodings: a plain FFT reco will return a pixel-perfect image without half-voxel shift. (This also means that for even sample counts, the phantom extends are not exactly centered. On the other hand, there is always a voxel exactly at (0, 0, 0)).

NIfTI coordinate system

The “true” specification of NIfTI seems to be the C header: nifti1_h.pdf - search for section “3D IMAGE (VOLUME) ORIENTATION AND LOCATION IN SPACE”.

The specification says that NIfTI’s are in RAS+:

the continuous coordinates are referred to as (x,y,z). The voxel index coordinates […] are referred to as (i,j,k)
[…]
The (x,y,z) coordinates refer to the CENTER of a voxel. In methods 2 and 3, the (x,y,z) axes refer to a subject-based coordinate system, with +x = Right +y = Anterior +z = Superior. This is a right-handed coordinate system.
[…]
The i index varies most rapidly, j index next, k index slowest.

The 3 Methods this refers to are:

1. qform_code == 0 && sform_code == 0: compat for ANALYZE files, only supports scaling
2. qform_code > 0: use a quaternion for additional rotation
3. sform_code > 0: full affine matrix specifies rotation, scale, offset (and even skewing)

Problem
NIfTI’s allow qform and sform to co-exist and doesn’t specify which one to use. It states:

In this scheme, a dataset would originally be set up so that the Method 2 coordinates represent what the scanner reported. Later, a registration to some standard space can be computed and inserted in the header. Image display software can use either transform, depending on its purposes and needs.

nibabel defaults to use sform and in normal use never sets different transforms for sform and qform, expect when set manually.

Mapping

The affine matrices map $[i,j,k] \mapsto (x, y, z)$. Even though the standard says that $(x, y, z)$ always refer to the subject-based RAS+ system, the origin of the coordinate system is specified twice, by the qcode_form and the scode_form:

Conclusion

NIfTI’s written with nibabel always use (unless forced) sform_code == 2: affine matrix with arbitrary coordinates, since the scanner system is typically not known when creating phantoms from code. Similarly, when reading files, we usually have no way of transforming between scanner and subject coordinate systems, since this information is not stored in the NIfTI…

Warning

…except if both sform and qform are used and one of them selects scanner coordinates and the other one subject coordinates. We should not use this in our generated phantoms, but have to take care when loading NIfTI files that use this property:

Different viewers could decide differently which of the both systems to use - double check when debugging a wrong orientation!

This information is also found in the NIfTI FAQ - Q17 and Q19

State selection

Documentation on PDG State Parameters and Configuration

Overview

The Phase Distribution Graph (PDG) states are frequently initialized with default parameters. This document outlines the essential considerations and configurations required for optimal performance, especially when dealing with exotic sequences.

Standard PDG Configuration

Typically, PDG states are invoked using the following standard parameters:

graph = mr0.compute_graph(seq0, obj_ph)
signal = mr0.execute_graph(graph, seq0, obj_ph)

The standard values used then are here written explicitly:

graph = mr0.compute_graph(seq0, obj_ph, max_state_count=200, min_state_mag=0.0001)
signal = mr0.execute_graph(graph, seq0, obj_ph, min_emitted_signal=0.01, min_latent_signal=0.01)

Importance of State Selection Thresholds

For arbitrary sequences, careful examination of state selection is crucial. This process remains important for all users employing PDG for scientific purposes.

As outlined in the main paper, there are four main parameters to adjust:

1. Threshold Number of States (max_state_count)
2. Threshold Magnetization (min_state_mag)
3. Threshold Emitted Signal (min_emitted_signal)
4. Threshold Latent Signal (min_latent_signal)

To resolve issues of inaccurate simulations one can adjust the parameters as follows to get a more accurate simulation:

graph = mr0.compute_graph(seq0, obj_ph, max_state_count=5000, min_state_mag=1e-12)
signal = mr0.execute_graph(graph, seq0, obj_ph, min_emitted_signal=0.001, min_latent_signal=0.001)

This adjustment will result in a longer simulation duration, but more accurate results

We recommend periodically validating results with more precise simulations.

Each sequence may require unique adjustments to balance simulation speed and accuracy. It’s often beneficial to start with a higher state count for a more accurate baseline.

For sequences with poor timing or spoiling, state explosion can occur (Hennig 1991). The state selection thresholds will filter out these states, potentially leading to inaccurate simulation results.

API overview

All available functions (snake_case naming) / classes (CamelCase) are listed here. In the following pages, more information is listed. For detailed descriptions of each items, look at their Python docstrings (as written in the source code or shown by your Python IDE) as well.

import MRzeroCore as mr0

• Sequence building blocks
  • mr0.Sequence
  • mr0.chain
  • mr0.Repetition
  • mr0.Pulse
  • mr0.PulseUsage
• Simulation data
  • mr0.SimData
  • mr0.VoxelGridPhantom
  • mr0.TissueDict
  • mr0.CustomVoxelPhantom
• NIfTI phantoms
  • mr0.PhantomUnits
  • mr0.PhantomSystem
  • mr0.NiftiRef
  • mr0.NiftiMapping
  • mr0.NiftiTissue
  • mr0.NiftiPhantom
• Simulation
  • mr0.isochromat_sim
  • mr0.compute_graph
  • mr0.compute_graph_ext
  • mr0.Graph
  • mr0.execute_graph
• Reconstruction
  • mr0.reco_adjoint
• util module
  • mr0.util.get_signal_from_real_system
  • mr0.util.insert_signal_plot
  • mr0.util.pulseq_plot
  • mr0.util.pulseq_plot_142
  • mr0.util.pulseq_plot_pre14
  • mr0.util.imshow
  • mr0.util.load_phantom
  • mr0.util.simulate
  • mr0.util.simulate_2d
• mr0.sig_to_mrd
• mr0.pulseq_write_cartesian
• mr0.generate_brainweb_phantoms

MRI Sequence definitions

The MRzeroCore Sequence is a Python list of Repetitions, each of which is starting with an instantaneous Pulse. To construct sequences, you can

• load from .seq files: Pulseq Integration
• construct with pulseq-zero: Pulseq-zero
• defined in the MR-zero directly:

import MRzeroCore as mr0
import numpy as np

seq = mr0.Sequence()

# Iterate over the repetitions
for i in range(64):
    rep = seq.new_rep(2 + 64 + 2)
    
    # Set the pulse
    rep.pulse.usage = mr0.PulseUsage.EXCIT
    rep.pulse.angle = 5 * np.pi/180
    rep.pulse.phase = i**2 * 117*np.pi/180

    # Set encoding
    rep.gradm[:, :] = ...

    # Set timing
    rep.event_time[:] = ...

    # Set ADC
    rep.adc_usage[2:-2] = 1
    rep.adc_phase[2:-2] = np.pi - rep.pulse.phase

Sequence

Defines an MRI sequence. Derived from a Python list.

Sequence(repetitions, normlized_grads=True)

• Sequence.cpu(): move sequence data to CPU
• Sequence.cuda(device=None): move sequence data to selected or default CUDA device
• Sequence.device (property): the device of the repetition data
• Sequence.clone(): make a copy of the repetition with cloned
• Sequence.new_rep(event_count): appends a Repetition.zero(event_count) to the sequence
• Sequence.get_kspace(): Returns the kspace trajectory of the measured signal as a single tensor (ADC events only, concatenated repetitions)
• Sequence.get_full_kspace(): Returns the full kspace trajectory as a list of per-repetition tensors (containing all events)
• Sequence.plot_kspace_trajectory(): Plot the result of get_full_kspace()
• Sequence.get_contrast_mask(contrast): Returns a bool mask for the selected contrast which can be applied to the simulated signal or kspace
• Sequence.get_contrasts(): Returns a list of all used contrast IDs (apply for all repetitions)
• Sequence.shift_contrasts(offset): Shift all contrast IDs by offset (apply for all repetitions)
• get_duration(): Return the total duration of the sequence in seconds

Sequences can be imported from Pulseq .seq files or Siemens .dsv files:

Sequence.import_file(file_name, exact_trajectories, print_stats, default_shim, ref_voltage, resolution)

chain

Helper function to combine multiple sequences into a multi-contrast sequence.

chain(*sequences, oneshot=False)

Returns a new sequence object. If oneshot is set to:

• True, sequences are simply concatenated.
• False, the adc_usages are shifted so that every sub-sequence has their own set of IDs.

Repetition

Part of a Sequence, containing an RF pulse and a list of event_count events.

Repetition(pulse: Pulse, event_time, gradm, adc_phase, adc_usage)

• Repetition.event_count (property): number of events in this repetition
• Repetition.cpu(): move repetition data to CPU
• Repetition.cuda(device=None): move repetition data to selected or default CUDA device
• Repetition.device (property): the device of the repetition data
• Repetition.zero(event_count): create a zero’d repetition with the given number of events
• Repetition.clone(): make a copy of the repetition with cloned
• Repetition.get_contrasts(): return a list of all used adc_usage values (except zero)
• Repetition.shift_contrasts(offset): add offset to all adc_usage values (except zero)

Pulse

Contains the definition of an instantaneous RF Pulse.

Pulse(usage: PulseUsage, angle, phase, shim_array, selective)

• Pulse.cpu(): move pulse data to CPU
• Pulse.cuda(device=None): move pulse data to selected or default CUDA device
• Pulse.device (property): the device of the pulse data
• Pulse.zero(): alternative constructor which zero-initializes everything
• Pulse.clone(): make a copy of the pulse with cloned

PulseUsage

Used for computing the k-space trajectory or to select the pulse type in legacy exporters.

Simulation Phantoms

Phantoms based on BrainWeb data can be downloaded and generated easily, as explained in Generating Phantoms. The data required by the simulation defined by SimData. It should not be created directly, but rather by defining or loading a VoxelGridPhantom or CustomVoxelPhantom.

Example for loading a BrainWeb phantom, scaling it and using it as SimData on the GPU:

phantom = mr0.VoxelGridPhantom.load("subject05.npz")
phantom = phantom.interpolate(128, 128, 32).slices([16])
data = phantom.build().cuda()

Tip

A new type of phantoms is now officially supported! These are based on NIfTIs and allow partial volume effects, segmented phantoms and easy phantom reconfiguration without code changes. Take a look at their description and API.

Voxel Shape

Simulation data (and therefore all Phantoms) define a shape for their voxels. The following options are available:

Most cases should use "sinc":
Input maps, signals, and reconstructed images are all sampled as points. Extending this bandwidth-limited data to continuous functions means that sinc-shaped voxels are the correct form in basically all cases. Using box-shaped voxels means that a round trip of putting maps into the simulation and then quantifying them from the reconstructed images will inevitably blur them. Box shaped voxels do not have flat signal response like sinc shaped voxels have (up to the nyquist frequency). Point shaped voxels would be a good fit too, except that they would prohibit gradient spoiling (and are a bad description of actual tissue).

SimData

Contains the data used for simulation (together with the Sequence). The values are flattened tensors, which only contain voxels with a non-zero proton density.

# Usually not constructed directly but via the `.build()` function on a phantom
SimData(PD, T1, T2, T2dash, D, B0, B1, coil_sens, size, voxel_pos, nyquist, dephasing_func, recover_func=None, phantom_motion=None, voxel_motion=None, tissue_masks=None)

• SimData.cpu(): move data to the CPU
• SimData.cuda(): move data to the GPU
• SimData.device (property): the device of the repetition data
• SimData.recover(): reconstruct the VoxelGridPhantom or CustomVoxelPhantom that was used to build this SimData

VoxelGridPhantom

Phantom format where tissue parameters are cartesian grids of voxels. The data is stored in 3D or 4D tensors with indices [(channel), x, y, z].

VoxelGridPhantom(PD, T1, T2, T2dash, D, B0, B1, coil_sens, size, phantom_motion=None, voxel_motion=None, tissue_masks=None)

Parameters of this function are identical to those of SimData, except that they are stored in 3D or 4D tensors. They are indexed [(channel), x, y, z].

• VoxelGridPhantom.load(file_name): Loading phantom from a .npz file
• VoxelGridPhantom.load_mat(file_name, ...): Loading phantom from a .mat file
• VoxelGridPhantom.slices(slices): Removing all but the specified list of slices (in z direction)
• VoxelGridPhantom.scale_fft(): Resizing the phantom by doing FFT -> truncate -> FFt
• VoxelGridPhantom.interpolate(): Resizing the phantom by using torch.nn.functional.interpolate(..., mode='trilinear')
• VoxelGridPhantom.plot(): Plotting all maps of the phantom (center slice)
• VoxelGridPhantom.build(PD_threshold): Convert this phantom into SimData, sparsifying the arrays (flattening and masking where pd < PD_threshold)

TissueDict

A tissue dictionary can store multiple, overlapping tissues (which lead to partial volume effects). It is a dictionary of VoxelGridPhantoms and provides additional helper functions.

TissueDict is modelled to hold the data loaded from a NiftiPhantom. To load such a file, you can call TissueDict.load() with the config .json of the NIfTI phantom. Alternatively, you can create a NiftiPhantom and then load the phantom (and the volume data) with TissueDict.load() - this time passing a folder and the NiftiPhantom.

TissueDict(dict[str, VoxelGridPhantom])

• TissueDict.load(path, config=None): either load from .json or from a dir + NiftiPhantom
• TissueDict.save(path, gyro, B0): save as NIfTI phantom
• TissueDict.interpolate(): apply [`VoxelGridPhantom.interpolate()] to all tissues
• TissueDict.slices(): apply [`VoxelGridPhantom.slices()] to all tissues
• TissueDict.combine(): merge all tissues into a single VoxelGridPhantom
• TissueDict.build(PD_threshold, voxel_shape): convert into SimData

CustomVoxelPhantom

Phantom format where tissue parameters are 1D lists of voxels with arbitrary positions.

CustomVoxelPhantom(pos, PD=1.0, T1=1.5, T2=0.1, T2dash=0.05, D=1.0, B0=0.0, B1=1.0, voxel_size=0.1, voxel_shape="sinc")

The size of this phantom is computed from the extends of the voxel positions.
Useful for testing reconstruction, sub-voxel positioning, teaching and more.

Parameters of this function are identical to those of SimData, except that they can also be python lists (and are converted to tensors internally). In addition, it contains a convenience method for plotting the phantom by rendering the voxels (and their selected shape) at any resolution.

• CustomVoxelPhantom.generate_PD_map(): Voxels can be sinc- box- or gauss-shaped and at arbitrary positions. This function renders them into a fixed-resolution proton density map.
• CustomVoxelPhantom.plot(): Render and plot all maps of this phantom
• CustomVoxelPhantom.build(): Convert into SimData instance for simulation

NIfTI Phantoms

The classes shown here reflect the content of the .json config file structure. These files define how various NIfTIs are combined into a phantom. You can load a complete phantom directly with mr0.TissueDict.load. But you can also load only the configuration with mr0.NiftiPhantom.load. This allows you to modify it as needed before loading the phantom or saving the changes.

All classes are dataclasses, augmented with additional functions.

PhantomUnits

Stores the physical units used by the phantom. These are currently fixed to the values used by MR-zero (default values) and are stored for documentation (no conversion is implemented). However, different units and automatic conversion might be implemented in the future.

PhantomUnits(gyro, B0, T1, T2, T2dash, ADC, dB0, B1_tx, B1_rx)

• PhantomUnits.default(): return PhantomUnits with the default units given above
• PhantomUnits.from_dict(config): load from a dictionary and check if units are valid
• PhantomUnits.to_dict(): convert to a dictionary

PhantomSystem

Describes the physical properties of the MRI experiment.

PhantomSystem(gyro, B0)

• PhantomSystem.from_dict(config): load from a dictionary
• PhantomSystem.to_dict(): convert to a dictionary

NiftiRef

Reference to a NIfTI file, including the index to the tissue. NIfTIs can contain multiple 3D volumes, stacked in the 4th dimension. For a description look at the specification.

NiftiRef(file_name, tissue_index)

As given by the spec, NiftiRefs are stored as "path/to/nifti.nii.gz[<index>]" strings:

• NiftiRef.parse(config): parse from a string
• NiftiRef.to_str(): convert to a string

NiftiMapping

Combination of NiftiRef and a mapping function for simple modifications. For a description look at the specification.

NiftiMapping(file, func)

• NiftiMapping.parse(config): parse from a dictionary (calling NiftiRef.parse)
• NiftiMapping.to_dict(): convert to a dictionary

NiftiTissue

Definition of a single tissue. Each scalar property is given by a float, NiftiRef or NiftiMapping. B1_tx and B1_rx are lists to support multi-channel data. Units are given by PhantomUnits.

NiftiTissue(density, T1, T2, T2dash, ADC, dB0, B1_tx, B1_rx)

where Prop = float | NiftiRef | NiftiMapping

• NiftiTissue.default(): return a NiftiTissue with the default units given above - density must be passed
• NiftiTissue.from_dict(config): load from a dictionary
• NiftiTissue.to_dict(): convert to a dictionary

NiftiPhantom

Represents the configuration given by a phantom .json file.

NiftiPhantom(units, system, tissues)

The class attribute schema_uri is always "https://github.com/mrx-org/bifti-phantoms/bifti-phantom-v1" and is not a constructor parameter.

• NiftiPhantom.default(gyro=42.5764, B0=3): return a NiftiPhantom without tissues
• NiftiPhantom.load(path): load from .json file, given by the path
• NiftiPhantom.save(path): save to a .json file, given by the path
• NiftiPhantom.from_dict(config): load from a dictionary
• NiftiPhantom.to_dict(): convert to a dictionary

Simulation

MR-zero provides two simulations:

• a purposefully simple isochromat based simulation to check for correctness
• a state-of-the-art PDG simulation, built for speed and accuracy.

Both simulations are differentiable and capable to run on any device supported by PyTorch (all CPUs and NVIDIA GPUs). While the isochromat simulation is rather slow and only there for confirmation purposes, PDG is one of the fastest currently available simulations. You can find more information about it in the Paper.

It is strongly adviced to use the PDG simulation - isochromat_sim only exists for testing purposes.

Phase Distribution Graph (PDG) simulation

Phase Distribution Graphs split the magnetization into multiple states, forming a Graph over the duration of the sequence. More details can be found in the PDG Paper, once published.

The simulation runs in two passes:

• The pre-pass (compute_graph) simulates a large number of states, but only roughly estimates the signal. This information is used to assess which states are important for the signal, and which can be skipped because they don’t contribute anything.
• The main-pass (execute_graph) uses the exact signal equation and takes all voxels into account, but simulates as few states as possible. It does so by using the information generated by the pre-pass.

The precision / speed tradeoff for both passes can be chosen as follows:

# Default parameters: max 200 states per repetition
# with at leat 0.0001 of maximum possible magnetization
graph = mr0.compute_graph(seq, data, 200, 1e-4)
# Minimum emitted and latent signal of 0.01 (normalized)
signal = mr0.execute_graph(graph, seq, data, 0.01, 0.01)

The pre-pass limits the state count by both the maximum number (higher = more states / slower + more precise simulaiton) of states and the minimum magnetization (lower = more states). In practice, it turns out that the second parameter is less surprising in how it influences the quality of the simulation, which means that the first parameter should be set to a high number and the second parameter should be increased if the pre-pass is too slow (usually the main-pass is more likely to be the bottleneck).

The main-pass has two parameters which are used as threshold to determine which states are used in the signal equation and which states are simulated. Usually they are both set to the same number. Lower numbers mean that more states are above the thresold and will be simulated. It is not useful to set min_emitted_signal to a value lower than min_latent_signal, as this is trying to include to states into the signal equation that are not simulated.

Examples:

# Very precise simulation
graph = mr0.compute_graph(seq, data, 10000, 1e-7)
signal = mr0.execute_graph(graph, seq, data, 1e-7, 1e-7)

# Average precision
graph = mr0.compute_graph(seq, data, 10000, 1e-4)
signal = mr0.execute_graph(graph, seq, data, 1e-3, 1e-4)

# FID only
graph = mr0.compute_graph(seq, data, 10000, 1e-3)
signal = mr0.execute_graph(graph, seq, data, 1, 0.1)

Optimal settings can always depend on the sequence itself. Some sequences are more demanding and require a very precise simulation, some are already decently described by only a few states.

compute_graph

This utility function computes the parameters for compute_graph_ext from the passed data. Returns the computed Graph.

compute_graph_ext

Compute the Phase Distribution Graph for the given sequence and settings. Runs a simulation of seq for a single voxel with the given properties. Discards states with magnetization below min_state_mag and the states with the lowest magnetization if surpassing max_state_count. Returns the computed Graph.

execute_graph

Takes the Sequence, the SimData and the Graph and calculates the measured ADC signal. Because of the work done by the pre-pass, only the minimal work needed in order to achieve the desired precision is executed. This precision can be tuned by the min_emitted_signal and min_latent_signal thresholds. Higher values lead to less states being simulated, which improves speed and reduces accuracy. A value of 1 will mean that 2 states will be simulated (z0 and one +), resulting in the FID signal. A value of 0 means that everything will be simulated that somehow contributes to the signal, even if vanishingly small.

Graph

A wrapper for the retuned graph, which provides two additional methods:

Graph.plot

Creates a matplotlib scatter plot of all PDG states. Does not create or show the figure, which allows to modify it, for example:

graph = compute_graph(seq, data)

plt.figure()
graph.plot(dephasing="k_x")
plt.grid()
plt.show()

Graph.plot_data

Returns the data used in Graph.plot.

isochromat_sim

Isochromat based Bloch simulation. The code is implemented for readability, not speed. This allows to verify its correctness, bugs are not expected.

Warning

This simulation

• is slow
• does not support diffusion
• always uses "box" - voxels, which is usually a bad choice and will introduce slight blurring. See Voxel Shape for more information. This choice was made for simplicity.

Reconstruction methods

NUFFT

Note

Currently, only one single reconstruction method exists here - despite more (like NUFFT or GRAPPA) being used in many scripts since a long time. We could integrate them here with MR-zero.

Nufft reconstruction is not integrated into mr0 directly, but can be realized with torchkbnufft and code similar to the following:

import torchkbnufft as tkbn

# Construct a sequence "seq" and simulate it, resulting in "signal"

kdata = signal.view(1, 1, -1)
ktraj = seq.get_kspace()[:, :2].T / (2 * np.pi)
im_size = (64, 64)

adjnufft_ob = tkbn.KbNufftAdjoint(im_size)
dcomp = tkbn.calc_density_compensation_function(ktraj, im_size=im_size)
reco = adjnufft_ob(kdata * dcomp, ktraj).view(im_size)

reco_adjoint

Adjoint reconstruction builds a very simple backwards version of the encoding / measurement operation. It essentially resembles a DFT and is thus slower and consumes more memory than an FFT, but can handle any readout trajectory.

util module

Additional helper functions.

import MRzeroCore as mr0
mr0.util.<function>()

get_signal_from_real_system

Wait for scanner TWIX files on the given path, and return if ready. Useful if a shared network folder between scanner and research PC is available.

Note

This is developed for internal use and might not be useful for others.

insert_signal_plot

Insert a measured signal into a currently open pypulseq plot. This is only supported by newer pulseq versions, which do not show (and close) the figure immediately if plot_now=False.

seq.plot(plot_now=False)  # option doesn't exist for older pulseq versions
mr0.util.insert_signal_plot(seq, signal)
plt.show()

pulseq_plot

Wrapper around insert_signal_plot, pulseq_plot_142 and pulseq_plot_pre14 which selects the appropriate plotting function based on the installed pulseq version.

pulseq_plot_142

Use only with pypulseq versions 1.4.x / tested with 1.4.2

Identical to pulseq plot(), but returns the generated axes! This is needed to insert the signal into the ADC plot. Modified from pypulseq 1.4.2 sequence.py - plot() - see pulseq documentation for details.

pulseq_plot_pre14

Use only with pypulseq versions before 1.4

Pulseq sequence plot with some aesthetical changes. Most important, it can insert the simulated signal into the ADC plot. Modified from pypulseq 1.2.0post1 sequence.py - plot() - see pulseq documentation for details.

imshow

Modified version of matplotlib imshow(), tailored for showing MR reconstructions. reco is expected to be an array with dimensions [x, y], [x, y, z] or [coil, x, y, z]

• plots images in RAS+
  • x is left-to-right
  • y is bottom-to-top
• plots 3d data (z) as a side-by-side grid of its slices
• applies quadratic coil-combination for multi-coil data

# Replace
import matplotlib.pyplot as plt
plt.imshow(reco, *args, **kwargs)

# with
import MRzeroCore as mr0
mr0.util.imshow(reco, *args, **kwargs)

load_phantom

Download a phantom from the MR-zero repository (or a custom URL) for getting started quickly.

DEFAULT_PHANTOM_URL = "https://github.com/MRsources/MRzero-Core/raw/main/documentation/playground_mr0/numerical_brain_cropped.mat"

simulate

Helper function to simulate using PDG.

simulate_2d

Caution

Deprecated - use load_phantom + simulate instead. For parameter descriptions, look at these functions.

Other methods

Miscellaneous methods provided by MR-zero. These are currently exported directly, but might be sorted into sub-modules if the list grows.

generate_brainweb_phantoms

Builds MR-zero phantoms with the data provided by BrainWeb. These phantoms are not stored in the repository; call this function once to generate them locally.

The configuration files with the used tissue parameters can be found here: brainweb_data.json

sig_to_mrd

Write the simulation output to an ISMRMD file; supported by many reconstruction tools.

pulseq_write_cartesian

Warning

Outdated - we suggest writing sequences with Pulseq(-zero) instead. Exporters might be removed in the future.

Export an MR-zero sequence to Pulseq. Since MR-zero sequences have no limitations in timing, slew rates, amplitudes and more, this exporter has to try to convert this into a sequence that can actually be measured. The exporter only supports cartesian trajectories due to its use of trapezoidal gradients. There is no guarantee that it produces the expected results.

MR-zero sign conventions

(from version 1.0.0 on)

1. Scope

This note defines the sign convention used in MR-zero for simulated MR signals, k-space encoding, FID formation, spectral reconstruction, and phase interpretation.

The central convention is:

$$ M_{xy}(t) = M_x(t) + iM_y(t) $$

For a proton with positive gyromagnetic ratio in a positive \(B_0\) field, the transverse magnetization rotates clockwise when viewed along \(+B_0\). In this convention, free precession is written as:

$$ M_{xy}(t)=M_{xy}(0)e^{-i\omega_0 t}, \qquad \omega_0=\gamma B_0 . $$

Thus, MR-zero uses the physical Bloch-equation sign convention:

$$ \boxed{ \text{physical MR signal evolution uses a negative complex exponent} } $$

This does not mean that every reconstruction step should numerically call fft.

It means that the signal formation / encoding kernel is FFT-like, while the reconstruction from raw samples back to object space is iFFT-like.

2. FFT terminology

Using the usual NumPy/MATLAB convention, the forward discrete Fourier transform is:

$$ A_k=\sum_{m=0}^{N-1}a_m \exp\left(-2\pi i,\frac{mk}{N}\right). $$

The inverse discrete Fourier transform is:

$$ a_m=\frac{1}{N}\sum_{k=0}^{N-1}A_k \exp\left(+2\pi i,\frac{mk}{N}\right). $$

In this note:

• FFT-type kernel means an \(e^{-i(\cdot)}\) kernel.
• iFFT-type kernel means an \(e^{+i(\cdot)}\) kernel.
• Scaling factors such as \(1/N\), \(2\pi\), fftshift, and ifftshift are treated separately.

Therefore, with the MR-zero physical convention,

$$ s(t)\sim e^{-i\omega t}, $$

a spectrum with a peak at \(+\omega\) is obtained using:

$$ \rho(\omega)=\int s(t)e^{+i\omega t},dt . $$

In NumPy/MATLAB-style code, this corresponds to a scaled ifft, not fft.

3. Imaging and spectroscopy have the same sign structure

The same physical sign governs both spatial encoding and spectral/FID evolution.

3.1 Imaging

Gradient-induced dephasing gives:

$$ s(\mathbf k) = \int \rho(\mathbf r) \exp(-i,\mathbf k\cdot\mathbf r),d\mathbf r . $$

Equivalently, if using cycles rather than radians:

$$ s(\mathbf k) = \int \rho(\mathbf r) \exp(-i2\pi,\mathbf k\cdot\mathbf r) ,d\mathbf r . $$

This is a forward Fourier encoding of the image into k-space.

Therefore image reconstruction is the inverse transform:

$$ \rho(\mathbf r) = \int s(\mathbf k) \exp(+i,\mathbf k\cdot\mathbf r),d\mathbf k . $$

So, under this convention:

$$ \boxed{ \text{k-space}\rightarrow\text{image uses an iFFT-type reconstruction} } $$

3.2 Spectroscopy

For spectroscopy, a spin ensemble with frequency distribution \(\rho(\omega)\) gives an FID:

$$ s(t) = \int \rho(\omega) \exp(-i\omega t),d\omega . $$

This is the same structure as imaging:

$$ s(\mathbf k) = \int \rho(\mathbf r) \exp(-i,\mathbf k\cdot\mathbf r),d\mathbf r . $$

The correspondence is:

$$ (\mathbf k,\mathbf r)\leftrightarrow(t,\omega). $$

Thus the reconstruction is also the inverse-sign transform:

$$ \rho(\omega) = \int s(t)\exp(+i\omega t),dt . $$

So, under the MR-zero convention:

$$ \boxed{ \text{FID}\rightarrow\text{spectrum uses an iFFT-type reconstruction} } $$

In code:

spec = np.fft.fftshift(np.fft.ifft(fid, axis=time_axis), axes=time_axis)
spec *= fid.shape[time_axis]   # optional: remove NumPy's 1/N ifft scaling

4. Fat, water, and \(B_0\) phase signs

Assume water is the reference frequency.

Water rotates clockwise at the reference frequency. Fat is more shielded and therefore precesses more slowly than water. Relative to water, fat has a negative frequency offset:

$$ \Delta f_{\text{fat}} < 0 . $$

Under the MR-zero convention,

$$ s(t)=e^{-i2\pi\Delta f t}. $$

Therefore, for fat,

$$ \Delta f_{\text{fat}}<0 \quad\Rightarrow\quad s_{\text{fat}}(t) = e^{+i2\pi|\Delta f_{\text{fat}}|t}. $$

So fat accumulates positive phase relative to water.

For a general \(B_0\) perturbation,

$$ \Delta f(\mathbf r)=\frac{\gamma}{2\pi}\Delta B_0(\mathbf r), $$

and the phase evolves as:

$$ \phi(\mathbf r,t) = -2\pi\Delta f(\mathbf r)t = -\gamma \Delta B_0(\mathbf r)t . $$

Therefore, under the MR-zero physical convention:

$$ \boxed{ \phi \text{ and } \Delta B_0 \text{ are anticorrelated} } $$

Positive \(\Delta B_0\) gives negative phase accumulation.

Negative frequency offsets, such as fat relative to water, give positive phase accumulation.

5. Relation to textbook and MRS-processing conventions

5.1 Keeler: physical precession versus processing convention

There is an apparent sign change in common NMR teaching material.

In the vector-model description, the physical Larmor precession is commonly written as:

$$ \omega_0=-\gamma B_0 . $$

For nuclei with positive gyromagnetic ratio, such as protons, this means the Larmor frequency is negative under that coordinate convention. The transverse magnetization rotates clockwise when viewed along \(+B_0\).

With:

$$ M_{xy}=M_x+iM_y, $$

this corresponds to:

$$ M_{xy}(t)=M_{xy}(0)e^{-i\omega t} $$

for a positive physical precession magnitude \(\omega>0\).

However, in Fourier-processing discussions of NMR spectroscopy, the complex FID is often written differently. The two receiver channels may be defined as:

$$ S_x(t)=S_0\cos\Omega t, \qquad S_y(t)=S_0\sin\Omega t. $$

Then the complex signal is combined as:

$$ S(t)=S_x(t)+iS_y(t). $$

This gives:

$$ S(t) = S_0\cos\Omega t+iS_0\sin\Omega t = S_0e^{+i\Omega t}. $$

In MR-zero terminology, the distinction is:

$$ \boxed{ \text{physical magnetization convention: } e^{-i\omega t} } $$

but:

$$ \boxed{ \text{receiver/processing convention: } e^{+i\Omega t} } $$

MR-zero intentionally follows the first convention, not the second.

6. Why MR-zero uses the physical sign

MR-zero simulates the MR signal directly in the rotating frame.

It does not simulate the complete analog/digital receive chain, including:

• coil voltage induction,
• mixer phase,
• receiver reference phase,
• quadrature-channel ordering,
• vendor-specific channel combination,
• storage-format convention.

Therefore MR-zero should not silently assume a particular receiver convention such as:

$$ S(t)=I(t)+iQ(t) $$

versus:

$$ S(t)=I(t)-iQ(t). $$

Those two choices differ by complex conjugation:

$$ S^*(t) \quad\Longleftrightarrow\quad e^{-i\omega t}\leftrightarrow e^{+i\omega t}. $$

The MR-zero convention is instead:

$$ \boxed{ \text{simulate the physical complex transverse magnetization } M_x+iM_y } $$

Thus, for a positive frequency offset,

$$ s(t)=e^{-i\omega t}. $$

This gives one consistent sign rule for both imaging and spectroscopy.

7. Relation to MRS literature saying “FFT”

Some MRS literature states that conversion of the discrete time-domain FID signal into a spectrum is performed using a discrete Fourier transformation such as the FFT, and that conversion from the spectral domain back to the time domain uses the inverse Fourier transform.

MR-zero deliberately differs from that numerical statement.

In MR-zero:

$$ s(t)=e^{-i\omega t}. $$

Therefore, a peak at \(+\omega\) is obtained using:

$$ \rho(\omega) = \int s(t)e^{+i\omega t}dt, $$

which is an iFFT-type operation.

So the MR-zero documentation should state explicitly:

$$ \boxed{ \text{MR-zero differs from the common MRS “FID}\rightarrow\text{spectrum = FFT” wording} } $$

This is not a disagreement about the physics of MR.

It is a disagreement about which complex signal is called “the FID.”

The MRS-processing convention usually assumes that the stored complex FID already follows a receiver/storage convention like:

$$ S_{\text{stored}}(t)\sim e^{+i\Omega t}, $$

so that a standard forward FFT places positive offsets at positive spectral frequencies.

MR-zero instead simulates:

$$ S_{\text{MR-zero}}(t)\sim e^{-i\Omega t}. $$

Thus, one possible conversion between the two conventions is:

$$ S_{\text{stored}}(t) = S_{\text{MR-zero}}^*(t). $$

8. Receiver and file-format conventions can still change the sign

Scanner data can differ from the MR-zero physical convention.

A receiver chain, reconstruction pipeline, or file converter may store either:

$$ s(t) $$

or:

$$ s^*(t). $$

Conjugation flips the exponent sign:

$$ s(t)=e^{-i\omega t} \quad\Rightarrow\quad s^*(t)=e^{+i\omega t}. $$

This changes the apparent reconstruction rule.

A truly iFFT-type reconstruction can look like an FFT-type reconstruction if the stored data have already been conjugated.

Thus conjugation both mirrors the frequency/spatial axis and inverts phase.

9. Explicit conversion knob

Users who want the common MRS receiver/storage convention can apply:

sig_receiver = np.conj(sig_mrzero)

The important design decision is that MR-zero does not build the receiver convention into the physical simulation. Instead, the conversion is explicit:

signal = np.conj(signal)

10. Imaging consequences of conjugated k-space

Let the physically expected k-space be:

$$ K(\mathbf k) =\mathcal F{\rho(\mathbf r)}. $$

Then the physically consistent reconstruction is:

$$ \rho(\mathbf r) =\mathcal F^{-1}{K(\mathbf k)}. $$

If the stored scanner data are instead:

$$ K_{\text{stored}}(\mathbf k)=K^*(\mathbf k), $$

then:

$$ \mathcal F^{-1}{K_{\text{stored}}} $$

gives a flipped and conjugated image, while using the opposite transform can fix the apparent orientation but leaves the phase convention changed.

In practical terms:

ifft(kspace)       -> correct reconstruction if kspace follows the physical convention
fft(kspace)        -> mirrored version of ifft(kspace), same phase after mirroring
ifft(conj(kspace)) -> mirrored image with inverted phase
fft(conj(kspace))  -> orientation may look fixed, but phase is conjugated

So if a scanner stores conjugated k-space, using fft instead of ifft can make the image orientation look correct, while the phase remains opposite to the MR-zero physical phase convention.

This should be treated as a vendor-, sequence-, and file-format-dependent issue rather than a universal law.

11. Spectroscopy consequence

For an MR-zero FID simulation with water at \(0\) Hz and fat at, for example,

$$ \Delta f_{\text{fat}}=-420\ \mathrm{Hz}, $$

the simulated FID should be:

$$ s(t) = A_{\text{water}} + A_{\text{fat}}e^{-i2\pi(-420)t}. $$

The correct MR-zero spectral reconstruction is then:

spec = np.fft.fftshift(np.fft.ifft(fid))
spec *= len(fid)

This gives fat at \(-420\) Hz relative to water.

Using fft directly would mirror the spectral axis and place fat at \(+420\) Hz unless the frequency axis is also reversed or the FID is conjugated.

Equivalent options are:

# MR-zero physical convention
spec = np.fft.fftshift(np.fft.ifft(fid)) * len(fid)

# Equivalent, but using conjugated data
spec = np.fft.fftshift(np.fft.fft(np.conj(fid)))

# Equivalent if one also flips the frequency axis
spec = np.fft.fftshift(np.fft.fft(fid))
freq = -np.fft.fftshift(np.fft.fftfreq(len(fid), dt))

Thus:

$$ \boxed{ s(t)=e^{-i2\pi ft} \quad\Rightarrow\quad \text{use iFFT if } f \text{ should appear at } +f} $$

12. Practical MR-zero rule

For MR-zero simulations, use the physical convention:

$$ M_{xy}=M_x+iM_y, \qquad M_{xy}(t)\sim e^{-i\omega t}. $$

Then the consistent reconstruction rules are:

The core conclusion is:

$$ \boxed{ \text{In MR-zero, both k-space and FID formation use the same }e^{-i(\cdot)} \text{ physics} } $$

Therefore:

$$ \boxed{ \text{Both image reconstruction and FID-to-spectrum reconstruction are iFFT-type operations} } $$

If scanner data require fft to look correct, this is evidence of a stored-data convention change, most commonly conjugation and/or axis reversal, not a different physical sign law.

13. Suggested validation experiments

The following experiments can reveal which convention is being used by a simulator, scanner export, or reconstruction package.

13.1 Spectroscopy FID test

Simulate:

$$ s(t)=1+ae^{-i2\pi(-420)t}. $$

Expected MR-zero result:

• ifft(fid) gives fat at \(-420\) Hz.
• fft(fid) gives fat at \(+420\) Hz unless the frequency axis is reversed.
• fft(conj(fid)) gives the same frequency placement as ifft(fid) up to scaling and shift conventions.

13.2 \(B_0\) phase test

Simulate a positive \(B_0\) offset:

$$ \Delta B_0>0. $$

Expected MR-zero phase:

$$ \phi(t)=-\gamma\Delta B_0t. $$

So positive \(\Delta B_0\) should produce negative phase.

13.3 Fat/water GRE test

For a GRE acquisition with water and fat:

$$ \Delta f_{\text{fat}}<0. $$

Expected MR-zero phase:

$$ \phi_{\text{fat-water}}(TE) =-2\pi\Delta f_{\text{fat}}TE>0. $$

Thus fat should accumulate positive phase relative to water.

13.4 Imaging k-space conjugation test

Compare:

img_ifft      = ifftn(kspace)
img_fft       = fftn(kspace)
img_ifft_conj = ifftn(np.conj(kspace))
img_fft_conj  = fftn(np.conj(kspace))

Expected behavior:

• fft(kspace) mirrors the ifft(kspace) image.
• Conjugating k-space inverts phase.
• fft(conj(kspace)) can appear to fix orientation while preserving the conjugated phase convention.

These tests separate three effects that are often confused:

$$ \boxed{ \text{transform sign} \quad \text{vs.} \quad \text{axis reversal} \quad \text{vs.} \quad \text{complex conjugation} } $$

14. Final takeaway

The MR-zero sign change is consistent and should be kept.

$$ \boxed{ \text{MR-zero simulates the physical } M_x+iM_y \text{ signal directly} } $$

Therefore:

$$ \boxed{ s(t)=e^{-i\omega t} } $$

and:

$$ \boxed{ \text{FID}\rightarrow\text{spectrum uses an iFFT-type transform} } $$

This is the same sign logic as spatial MRI encoding:

$$ s(\mathbf k)=\int\rho(\mathbf r)e^{-i\mathbf k\cdot\mathbf r}d\mathbf r \quad\Rightarrow\quad \rho(\mathbf r)=\int s(\mathbf k)e^{+i\mathbf k\cdot\mathbf r}d\mathbf k . $$

In one sentence:

$$ \boxed{ \text{MR-zero keeps the physics sign; receiver, vendor, and MRS-software signs are explicit conversions} } $$

References

• Keeler, Understanding NMR Spectroscopy:
  https://faculty.washington.edu/seattle/physics541/%202010-reading/nmr-spectroscopy.pdf

• Near et al., MRS preprocessing recommendations, NMR in Biomedicine:
  https://analyticalsciencejournals.onlinelibrary.wiley.com/doi/10.1002/nbm.4257

• NumPy FFT conventions:
  https://numpy.org/doc/stable/reference/routines.fft.html

• MRI signal equation notes:
  https://larsonlab.github.io/MRI-education-resources/MRI%20Signal%20Equation.html

• NIfTI-MRS specification:
  https://github.com/wtclarke/mrs_nifti_standard/blob/master/specification.MD

Pulseq-zero

Pulseq-zero allows to define MRI sequences with the Python¹ port of Pulseq²: PyPulseq³, and use them within MR-zero⁴. This way they are deeply integrated in a differentiable digital twin, enabling not only the simulation of the defined sequence but also the efficient optimization of any sequence parameter and any loss function, using the power of PyTorch⁵ and gradient-descent with backpropagation.

Pulseq-zero uses PDG⁶, a fast, analytical and physically exact simulation model that calculates signals that are comparable to in-vivo measurements within seconds. At the same time, the required changes to the sequence code are minimal; Pulseq-zero exports the optimized sequence works by simply using the installed PyPulseq without any interference.

Table of contents

1. General Information
2. Usage
3. Development
4. API
5. References

1. General Information

Pulseq-zero can be cloned from this repository but is also hosted on PyPI, install it locally with:

pip install pulseqzero

Note

Pulseq-zero does not declare any runtime dependencies, but it expects pypulseq, torch, MRzeroCore, numpy, and matplotlib to already be in the environment. Pulseq-zero 1.0 is compatible with PyPulseq 1.4.x and 1.5.x. Version-specific kwargs (freq_ppm, no_signal_scaling, use_block_cache, etc.) are detected at import time via inspect.signature, so the same pulseq-zero installation adapts automatically to whichever PyPulseq version is in the environment.

Migration from 0.x: the mode-switching facade has been removed. Replace import pulseqzero; pp = pulseqzero.pp_impl with import pulseqzero as pp, and drop any with pulseqzero.mr0_mode(): wrappers — seq.to_mr0() and seq.write() both work unconditionally now.

Pulseq-zero was displayed at ESMRMB 2024! You can view the abstract here: abstract/abstract.md. This project is affiliated with MR-zero and PDG but none of the other technologies. It relies on the following amazing projects:

• Python is the programming language used for Pulseq-zero
• Pulseq is a vendor-agnostic library and file format for sequence definition and transfer to real systems
• PyPulseq is the port of Pulseq to Python
• MR-zero is a digital twin of the full measurement and reconstruction pipeline for sequence optimization and discovery
• PyTorch is an ecosystem of tools for efficient tensor math with GPU accelleration, autograd through backpropagation and a wide variety of optimizers
• PDG (short for Phase Distribution Graphs) is a state-of-the-art Bloch simulation that produces accurate MRI signals for any sequence, orders of magnitude faster than other approaches

2. Usage

Pulseq-zero is a drop-in replacement for PyPulseq: any existing PyPulseq script runs under pulseq-zero by swapping the import — or, if you don’t want to touch the script at all, by installing a two-line sys.modules hijack in the driver (see §2 · Reuse an unmodified PyPulseq script). The same script then writes .seq files and is consumed differentiably by MR-zero — no context managers, no mode flags.

import pulseqzero as pp

# Build the sequence exactly like a PyPulseq script.
seq = pp.Sequence()
seq.add_block(pp.make_delay(10e-3))

Define the sequence as a function

Wrap the sequence code in a function so the same definition can drive both .seq export and MR-zero simulation / optimization:

def my_gre_seq(TR, TE):
    seq = pp.Sequence()
    # ... create your sequence ...
    seq.add_block(pp.make_delay(TR - 3e-3))
    # ... more sequence creation ...
    return seq

Reuse an unmodified PyPulseq script

You don’t have to swap the import at all. If a sequence script still reads import pypulseq as pp, a two-line hijack in the driver redirects every pypulseq attribute lookup to pulseq-zero:

import sys
import pulseqzero
sys.modules["pypulseq"] = pulseqzero

from my_existing_pypulseq_script import main  # its `import pypulseq as pp` now resolves to pulseqzero
seq = main()
seq.to_mr0()   # differentiable — no edits to the sequence script

The demo/ workspace uses exactly this pattern: demo/write_tse.py is a near-verbatim copy of the PyPulseq 1.5 upstream TSE example and still imports pypulseq; demo/main.py installs the hijack before importing it and gets a differentiable sequence for free.

Caveats. The hijack must run before the downstream script is first imported. from pypulseq.submod import X (submodule access, e.g. pypulseq.convert) will fail unless pulseq-zero mirrors submod — top-level import pypulseq and from pypulseq import foo both work. Calls to entry points pulseq-zero deliberately doesn’t wrap (adiabatic pulses, sigpy, SLR, etc.) still raise NotImplementedError with a named workaround (see §4).

Application

• Export a .seq file and plot (goes through PyPulseq under the hood; a one-off translation warning is emitted so you notice if it fires inside a hot loop):

  seq = my_gre_seq(14e-3, 5e-3)
  seq.plot()
  seq.write("tse.seq")
  

• Simulate with MR-zero:

  import MRzeroCore as mr0
  
  seq = my_gre_seq(14e-3, 5e-3).to_mr0()
  graph = mr0.compute_graph(seq, sim_data)
  signal = mr0.execute_graph(graph, seq, sim_data)
  reco = mr0.reco_adjoint(signal, seq.get_kspace())
  

• Optimize sequence parameters with PyTorch:

  TR = torch.tensor(14e-3, requires_grad=True)
  TE = torch.tensor(5e-3, requires_grad=True)
  optimizer = torch.optim.Adam([TR, TE], lr=0.001)
  
  for _ in range(100):
      optimizer.zero_grad()
      seq = my_gre_seq(TR, TE).to_mr0()
      loss = my_loss(seq)
      loss.backward()
      optimizer.step()
  
  # After optimization: export using the same script.
  my_gre_seq(TR, TE).write("tse_optim.seq")
  

3. Development

The recommended dev toolchain is uv. Install it once (install instructions), then run the demos straight from the repo root:

uv run demo/main.py        # end-to-end optimization demo (needs demo/brain.npz)
uv run demo/write_tse.py   # build a TSE sequence, plot it, and emit tse_pypulseq.seq

uv run resolves the demo workspace member defined in pyproject.toml, installs the pinned demo deps (PyPulseq 1.5.0.post1, torch, MRzeroCore, matplotlib) into .venv/, and picks up the editable pulseqzero checkout — no manual pip install -e . step.

Good to know:

• No unit test suite, no linter. CI (.github/workflows/adapter-check.yml) runs demo/test_pypulseq_sequence_all_funcs.py on every push/PR, an API completeness/parity probe against the coverage table below — not a numerical-correctness test suite. The two demo scripts are the closest thing to an acceptance gate: demo/main.py must complete 30 Adam iterations with non-NaN data loss and monotonically-decreasing SAR, and demo/write_tse.py must produce a .seq file that round-trips byte-for-byte against the pypulseq reference (the unified adapter guarantees this; see to_pypulseq() in wrapper/sequence.py).
• PyTorch CUDA pin. demo/pyproject.toml references the CUDA 12.6 wheel index (download.pytorch.org/whl/cu126). Swap that index URL (or remove it) if you’re on CPU-only or a different CUDA version.
• Headless plotting. seq.plot() forwards to pypulseq’s plot and expects an interactive matplotlib backend. Export MPLBACKEND=Agg to run headless (Agg will render but not show — useful for CI-style runs).
• Falling back to plain pip. If you’d rather skip uv, pip install --editable . from the root still works — but the root declares no runtime deps, so you need pypulseq==1.5.0.post1, torch, MRzeroCore, numpy, and matplotlib already in the env (a venv created with --system-site-packages and an existing MR-zero install is the path of least resistance).
• Warnings are single-fire. seq.to_pypulseq() (and seq.write() / seq.plot() / other forwarders) emit a UserWarning the first time they run from a given call site, then stay quiet — Python’s default warning filter dedupes by (message, module, lineno). If the warning fires inside an optimization loop, move the call out of the loop.

4. API

Pulseq-zero provides the whole pypulseq 1.5 API, with some notable exceptions: Some functions are not differentiable (like plotting) and just re-exports. A few methods are not provided (like sigpy pulse optimization) as they are not compatible with the approach taken by pulseq- zero.

The full API can be found in TOC.md.

Differentiable rounding

PyPulseq aligns many events to the block / gradient / ADC raster, which requires rounding — and rounding kills gradients. Pulseq-zero ships pp.round / pp.ceil / pp.floor that match PyTorch semantics but act like the identity function on the backward pass:

my_param = torch.tensor(1.5, requires_grad=True)
some_calc = pp.round(torch.sin(my_param))
some_calc.backward()

assert some_calc == 1
assert my_param.grad == torch.cos(my_param)

Use these whenever you round a timing (or any sequence quantity) that flows from an optimization parameter. For plain numeric rounding outside optimization, np.round / torch.round are fine.

See MATH.md for the full set of differentiable math helpers, including round_half_up and the internal interp used for gradient waveform resampling.

seq.to_mr0() and seq.write()

Every pulseqzero.Sequence supports both paths unconditionally:

• mr0_seq = seq.to_mr0() — build an MRzeroCore.Sequence for PDG simulation / optimization.
• seq.write("out.seq") — translate the internal event graph through PyPulseq and emit a .seq file. A one-time warnings.warn is raised per call so you notice if it fires inside a hot loop (move it out of the optimizer).

If you need a native PyPulseq Sequence for a one-off exotic call, seq.to_pypulseq() is the explicit escape hatch.

System (Opts) is configured once, on the Sequence

Pulseq-zero applies one system to the whole sequence: the system= you pass to pp.Sequence(system=...) (or Opts.default when you pass none). That system is authoritative — it is the one used for every event at to_mr0(), write(), and to_pypulseq().

A system= handed to an individual make_* call is honored only for that call’s construction-time math (e.g. deriving a slew-limited rise_time); it is not stored on the event and not carried to conversion/export. Consequently, pulseq-zero does not support mutating the system between calls, nor using a different/derated Opts for individual events: at export the Sequence’s system wins, so an event built under different limits can silently disagree with it (or fail PyPulseq’s re-validation). If you genuinely need per-event limits, do that one-off through seq.to_pypulseq().

This matches the universal pulseq idiom — create one Opts at the top, pass it to the Sequence and (by convention) to each make_* call; since it is the same object, everything stays consistent. Every script in the PyPulseq example suite checked (write_gre, write_epi, write_tse, write_haste, write_radial_gre, write_ute, write_epi_se_rs) creates exactly one Opts and reuses it throughout.

Differentiability

Gradients flow through the following quantities end-to-end (set requires_grad=True and they thread through to seq.to_mr0()):

• RF flip_angle, phase_offset, freq_offset, delay
• ADC phase_offset, freq_offset, delay, dwell
• Gradient amplitude (trapezoidal and arbitrary), rise_time, flat_time, fall_time, delay
• Block / repetition / TR / TE durations

The following are not differentiable today (they affect pulse shape, which is materialized eagerly via PyPulseq):

• Pulse shape parameters: duration when used to shape the envelope, time_bw_product, apodization, center_pos, slice_thickness (as it feeds shape generation), dwell for pulses
• Gradient waveform samples for arbitrary gradients (the scale is differentiable, the samples aren’t)
• Opts fields (max_grad, rasters, dead times) — intentionally numeric

Pulse-shape autograd can be added back per-factory via an opt-in flag if it ever becomes load-bearing.

pTx shim_array (optional)

Every RF factory (make_block_pulse, make_gauss_pulse, make_sinc_pulse, make_arbitrary_rf) accepts an optional shim_array kwarg: a per-channel transmit shim, carried through unchanged to to_mr0() / to_pypulseq(). It is only meaningful together with Martin’s pTx fork of pypulseq (which adds set_tx_mode); pulseq-zero re-exports set_tx_mode when that fork is detected at import time, and is a no-op for everyone else. Regular single-channel sequences can ignore this parameter entirely.

5. References

1. python programming language: https://www.python.org/ ↩

2. Layton K et al: Pulseq: A rapid and hardware-independent pulse sequence prototyping framework. MRM 2017, doi: 10.1002/mrm.26235 ↩

3. Keerthi SR et al: PyPulseq: A Python Package for MRI Pulse Sequence Design. JOSS 2019, doi: 10.21105/joss.01725 ↩

4. Loktyushin A et al: MRzero - Automated discovery of MRI sequences using supervised learning. MRM 2021, doi: 10.1002/mrm.28727 ↩

5. Paszke A et al: PyTorch: An Imperative Style, High-Performance Deep Learning Library. arxiv 2019, doi: 10.48550/arXiv.1912.01703 ↩

6. Endres J et al: Phase distribution graphs for fast, differentiable, and spatially encoded Bloch simulations of arbitrary MRI sequences. MRM 2024, doi: 10.1002/mrm.30055 ↩

pulseqzero.math: differentiable math helpers

This documents src/pulseqzero/math.py, the small set of differentiable helpers pulseq-zero uses internally (and exposes to users) wherever plain torch/numpy would otherwise break the optimization pipeline: rounding operations that would zero out gradients, and gradient/waveform interpolation that pypulseq itself does with numpy.

ceil, floor, round, and round_half_up are re-exported at the top level (pulseqzero.ceil, etc. — see the README’s “Differentiable rounding” section). interp is an internal helper used by wrapper/grad_funcs.py and is not re-exported; import it explicitly with from pulseqzero.math import interp if you need it directly.

Straight-through rounding: ceil, floor, round

• Forward pass: behaves exactly like the corresponding torch/numpy rounding operation.
• Backward pass: returns the incoming gradient unchanged (a straight-through estimator) — these functions act as the identity for autograd, even though the forward value is discrete.
• Use case: rounding a timing (or any other sequence quantity) that is derived from an optimized requires_grad=True parameter, without killing its gradient.
• Plain Python numbers / numpy arrays are also accepted: pulseq-zero falls back to np.ceil / np.floor / np.round when the input isn’t a tensor autograd can attach to.

import pulseqzero as pp

y = pp.ceil(x)   # differentiable version of torch.ceil
y = pp.floor(x)  # differentiable version of torch.floor
y = pp.round(x)  # differentiable version of torch.round

For plain numeric rounding outside of an optimization (nothing needs a gradient through it), np.round / torch.round are simpler and fine to use instead.

round_half_up(n, decimals=0)

Differentiable rounding that rounds halves away from zero instead of using torch.round’s / np.round’s banker’s rounding (round-half-to-even). Implemented as floor(abs(n) * 10**decimals + 0.5) / 10**decimals, so it inherits floor’s straight-through backward pass.

Known bug (matches pypulseq): because the sign is stripped before rounding and never restored, round_half_up returns the wrong sign for negative inputs (e.g. round_half_up(-2.5) does not give -3). This mirrors a bug in pypulseq’s own round_half_up, kept here intentionally for parity — do not rely on this function for negative values.

interp(x, xp, fp, left=None, right=None, tol=None)

Autograd-compatible 1D linear interpolation mirroring numpy.interp — differentiable in x, xp, and fp (unlike numpy.interp/torch, which have no autograd support for this at all). Unlike ceil/floor/round above, this is not a straight-through estimator: the forward value and its gradient both come from the real linear-interpolation formula.

• Linear within [xp[0], xp[-1]]; outside that range it returns left (for x < xp[0]) and right (for x > xp[-1]), defaulting to the edge values fp[0] / fp[-1] like numpy. Pass left=right=0 to zero-fill outside the support.
• xp must be sorted ascending.
• tol is a slack on the boundary test, in addition to what numpy.interp does: the default (a few ULP relative to xp’s magnitude) lets a query landing exactly on a boundary survive float rounding, e.g. when xp and x are computed in different dtypes. Pass tol=0 for a strict boundary when the query grid is built from the exact xp.

This is what lets pulseq-zero superimpose and re-sample gradient waveforms (arbitrary gradients, add_gradients, split_gradient_at, …) while keeping gradients flowing through both the waveform amplitudes and the timing.

Literature

This is a list of all papers and conference abstracts (known to us) which reference MR-zero. If you think something is missing, please open an issue or pull request to extend the list.

Papers

T1-weighting in Steady-State FLASH MRI–Diffusion Is Not Only Supportive but Mandatory for the Contrast
Weinmüller S, Charles D, Endres J, Leupold J, Gritsch F, Wagner F, Schneider R, Zaiss M.
Magn Reson Med. 2026; https://doi.org/10.1002/mrm.70443

Agentic MR sequence development: leveraging LLMs with MR skills for automatic physics-informed sequence development
Zaiss M, Aly A, Endres J, Dornstetter T, Weinmüller S, Maier A.
arXiv 2026; https://arxiv.org/abs/2604.13282

MR Simulation with Phase Distribution Graphs: Off-Resonant Pulse Response and Slice-Selection
Dietz F, Endres J, Weinmueller S, Link J, Zaiss M.
arXiv 2026; https://arxiv.org/abs/2606.15164

MR sequence design to account for nonideal gradient performance.
West D, Glang F, Endres J, Leitão D, Zaiss M, Hajnal J, Malik S.
Magn Reson Med. 2026; https://doi.org/10.1002/mrm.70093

Controlling sharpness, SNR, and specific absorption rate for 3D fast-spin echo at 7T by end-to-end learning
Dawood P, Blaimer M, Herrler J, Liebig P, Weinmüller S, Malik S, Jakob P, Zaiss M.
Magn Reson Med. 2025; https://doi.org/10.1002/mrm.30533

Phase distribution graphs for fast, differentiable, and spatially encoded Bloch simulations of arbitrary MRI sequences.
Endres J, Weinmüller S, Dang HN, Zaiss M.
Magn Reson Med. 2024; https://doi.org/10.1002/mrm.30055

MR-zero meets FLASH – Controlling the transient signal decay in gradient- and rf-spoiled gradient echo sequences.
Weinmüller S, Endres J, Dang HN, Stollberger R, Zaiss M.
Magn Reson Med. 2024; https://doi.org/10.1002/mrm.30318

MRzero - Automated discovery of MRI sequences using supervised learning.
Loktyushin A, Herz K, Dang HN.
Magn Reson Med. 2021; https://doi.org/10.1002/mrm.28727

MR-zero meets RARE MRI: Joint optimization of refocusing flip angles and neural networks to minimize T2-induced blurring in spin echo sequences.
Dang, H, Endres, J, Weinmüller S et al.
Magn Reson Med. 2023; https://doi.org/10.1002/mrm.29710

Slice-Profile-Enabled Phase Distribution Graphs for MRI Simulation
Hussain S, Hoinkiss DC, Huber J, Kuhlen V, Günther M.
arXiv 2026; https://arxiv.org/abs/2606.09233

Conference abstracts

Dietz F, Weinmüller S, Endres J, Zaiss M.
Advancing Phase Distribution Graphs: Time-effective simulations including shaped pulses and slice-selective excitation
In: DACH-ISMRM Annual Meeting. Bamberg, Germany; 2026.

Weinmüller S, Dang HN, Endres J, Glang F, Loktyushin A, Zaiss M.
A blurring-free 3D snapshot readout for fast CEST- or relaxation-prepared MRI
In: 24. Jahrestagung Der Deutschen Sektion Der ISMRM. ; 2022:6-7.

Dang H, Endres J, Weinmüller S, Glang F, Loktyushin A, Dörfler A, Maier A, Zaiss M.
Joint optimization of flip angle design and DenseNet parameters for reduced T2 blurring in TSE sequences
In: 24. Jahrestagung Der Deutschen Sektion Der ISMRM. Aachen, Germany; 2022:14-15.

ISMRM

Direct Encoded Signal Control with Phase Distribution Graphs for readout-tailored multipulse pTx
Weinmüller S, Endres J, Dawood P, Freudensprung M, Glang F, Nagel A, Zaiss M. ISMRM 2025, abstract 2753; https://archive.ismrm.org/2025/2753.html

Analytical simulation of motion in MRI using phase distribution graphs
Endres J, Zaiss M.
ISMRM 2025, abstract 4633; https://archive.ismrm.org/2025/4633.html

Non-idealized system (NIS) optimization of EPI sequences at ultra-high field
West D, Glang F, Endres J, Leitão D, McElroy S, Zaiss M, Hajnal J, Malik S.
ISMRM 2024, abstract 0528; https://archive.ismrm.org/2024/0528.html

Overcoming System Imperfections Using End-to-End MR Sequence Design.
West D, Glang F, Endres J, Zaiss M, Hajnal J, Malik S.
ISMRM 2023, abstract 1242; https://archive.ismrm.org/2023/1242.html

MRzero sequence generation using analytic signal equations as forward model and neural network reconstruction for efficient auto-encoding.
Weinmüller S, Dang HN, Loktyushin A, Glang F, Doerfler A, Maier A, Schölkopf B, Scheffler K, Zaiss M.
ISMRM 2021, abstract 1761; https://archive.ismrm.org/2021/1761.html

Advances in MRzero: supervised learning of parallel imaging sequences including joint non-Cartesian trajectory and flip angle optimization.
Glang F, Loktyushin A, Herz K, Dang HN, Deshmane A, Weinmüller S, Zaiss M et al.
ISMRM 2021, abstract 4200; https://archive.ismrm.org/2021/4200.html

MRzero – Automated invention of MRI sequences using supervised learning
Loktyushin A, Herz K, Glang F, Schölkopf B, Scheffler K, Zaiss M.
ISMRM 2021, abstract 4199; https://archive.ismrm.org/2021/4199.html

MRzero with dAUTOMAP reconstruction: automated invention of MR acquisition and neural network reconstruction
Dang H, Weinmüller S, Loktyushin A, Glang F, Dörfler A, Maier A, Schölkopf B, Scheffler K, Zaiss M.
ISMRM 2021, abstract 0333; https://archive.ismrm.org/2021/0333.html

ESMRMB 2025

Book of Abstracts ESMRMB 2025 Online 41st Annual Scientific Meeting 8–11 October 2025. https://doi.org/10.1007/s10334-025-01278-8

Breast tissue segmentation at 7T using Chemical Exchange Saturation Transfer (CEST)
Magda Duarte, Katharina Breininger, Katharina Tkotz4, Sebastian Bickelhaupt, Moritz Zaiss
abstract 041

Agentic MR sequence development - Leveraging LLMs with MR tools and tests for physics-informed sequence development
Moritz Zaiss, Jonathan Endres, Simon Weinmüller
abstract 046

Breast Digital Twin: simulating with MR-zero
Magda Duarte, Felix Dietz, Tobias Dornstetter, Jonathan Endres, Simon Weinmüller, Sebastian Bickelhaupt, Moritz Zaiss
abstract 055

End-to-End Optimization of Variable Flip Angle Schemes for Enhanced Apparent Resolution in 7T 3D Fast Spin Echo MRI
Peter Dawood, Martin Blaimer, Angelika Mennecke, Fraticelli Laura, Simon Weinmüller, Jonathan Endres, Peter Jakob, Moritz Zaiss
abstract 066

Highly accelerated snapshot readout for high-resolution CEST imaging
Martin Freudensprung, Patrick Liebig, Simon Weinmüller1, Magda Durate, Moritz Zaiss
abstract 076

Simulation of dynamic B0 with Phase Distribution Graphs
Jonathan Endres, Simon Weinmüller, Moritz Zaiss
abstract 119

Phase graph-based MRI simulation including off-resonant pulse response
Felix Dietz, Simon Weinmüller, Jonathan Endres, Moritz Zaiss
abstract 130

Towards B0 and B1 Robust Spiral Turbo Spin Echo Sequence at 7T
Simon Weinmüller, Jürgen Hennig, Felix Dietz, Peter Dawood, Jonathan Endres, Moritz Zaiss
abstract 233

Patient specific quantitative MR twins for synthetic previews and protocol planning
Deepak Charles, Simon Weinmüller, Fabian Wagner, Rainer Schneider, Jonathan Endres, Moritz Zaiss
abstract 359

Exact motion simulation vs. Retrospective application on static data
Jonathan Endres, Simon Weinmüller, Moritz Zaiss
abstract 445

Optimal control pulses and MIMOSA for CEST preparation at 7 T
Martin Freudensprung, Clemens Stilianu, Simon Weinmüller, Rudolf Stollberger, Moritz Zaiss
abstract 504

ESMRMB 2024

Book of Abstracts ESMRMB 2024 Online 40th Annual Scientific Meeting 2–5 October 2024. https://doi.org/10.1007/s10334-024-01191-6

Quantalizer: a quantitative digital twin to create and run tailored sequences
Simon Weinmüller1, Jonathan Endres, Hoai Nam Dang, Junaid Rajput, Felix Glang, Andreas Maier, Moritz Zaiss
abstract 011

MR-zero-tailored dummy refocusing pulses in TSE sequences reduce artifacts
Jacob Malich, Simon Weinmüller, Jonathan Endres, Peter Dawood, Moritz Zaiss
abstract 112

kT-points optimization using Pulseq and MR-zero at 7T
Martin Freudensprung, Simon Weinmüller1, Jonathan Endres,,Armin Nagel, Moritz Zaiss
abstract 422

Joint reconstruction of EPI and FLASH using Phase Distribution Graphs as a general MRI reconstruction model
Magda Duarte, Jonathan Endres, Junaid Rajput, Simon Weinmüller1, Moritz Zaiss
abstract 456

ESMRMB 2023

Book of Abstracts ESMRMB 2023 Online 39th Annual Scientific Meeting 4–7 October 2023. https://doi.org/10.1007/s10334-023-01108-9

Dynamic DREAM MRI: B0, B1 and Tx/Rx-phase mapping for assisting motion tracking systems
Tim Baum, Simon Weinmüller, Armin Nagel, Martin Vossiek, Moritz Zaiss
LT52

A simple pTx Pulseq extension for pulse-specific B1 shimming
Martin Freudensprung, Simon Weinmüller1, Jonathan Endres, Armin Nagel, Moritz Zaiss
LT53

A magnetization-prepared DREAM sequence for CEST imaging with an intrinsic dynamic B0 and B1 reference (CEST-MP-DREAM)
Tim Baum, Simon Weinmüller, Jonathan Endres, Patrick Liebig, Moritz Zaiss
LT68

DREAM-zero – Optimized variable flip angles for decreased image blurring in magnetization-prepared DREAM sequences
Simon Weinmüller, Tim Baum, Hoai Nam Dang, Jonathan Endres, Moritz Zaiss
P198

ESMRMB 2019

MRI zero: Fully automated invention of MRI sequences using supervised learning
Loktyushin A, Herz K, Glang F, Schölkopf B, Scheffler K, Zaiss M.
Magn Reson Mater Phy. 2019;32(Suppl 1):S119-S120. https://doi.org/10.1007/s10334-019-00754-2

ESMRMB 2020

Autoencoding T1 using MRzero for simultaneous sequence optimization and neural network training
Dang H, Loktyushin A, Glang F, Herz K, Dörfler A, Schölkopf B, Scheffler K, Maier A, Zaiss M.
Magn Reson Mater Phy. 2020;33(Suppl 1):S27-S28. https://doi.org/10.1007/s10334-020-00874-0