Skip to content

DCE-MRI Analysis Tutorial

DCE-MRI workflow: data loading, T1 mapping, signal-to-concentration conversion, model fitting, and BIDS export.

Prerequisites

  • Completed Getting Started tutorial
  • DCE-MRI data with known acquisition parameters
  • VFA data for T1 mapping (or pre-computed T1 map)

Using the CLI? Generate a config with osipy --dump-defaults dce > config.yaml, edit it, then run osipy config.yaml data.nii.gz. See How to Run Pipeline from YAML. The tutorial below covers the Python API for step-by-step control.

Background

DCE-MRI measures tissue perfusion by tracking the passage of a gadolinium-based contrast agent. The analysis pipeline involves:

  1. T1 Mapping: Estimate pre-contrast T1 relaxation times
  2. Signal Conversion: Convert MRI signal to contrast agent concentration
  3. Model Fitting: Fit pharmacokinetic models to extract parameters

Key parameters include:

Parameter Symbol Units Description
Transfer constant Ktrans min⁻¹ Rate of contrast transfer from plasma to EES
EES volume fraction ve fraction Extravascular extracellular space volume
Plasma volume fraction vp fraction Blood plasma volume in tissue

For more theory, see Understanding Pharmacokinetic Models.

Step 1: Load Your Data

import numpy as np
import osipy

# Load DCE-MRI data (returns PerfusionDataset)
dce_dataset = osipy.load_nifti("dce_4d.nii.gz")

# Load VFA images for T1 mapping
# These should be acquired at different flip angles
vfa_data = osipy.load_nifti("vfa_images.nii.gz")

# Define acquisition parameters
flip_angles = np.array([2, 5, 10, 15])  # degrees
tr = 5.0  # Repetition time in ms

# DCE timing
n_timepoints = dce_dataset.shape[-1]
temporal_resolution = 3.5  # seconds between frames
time = np.arange(n_timepoints) * temporal_resolution

print(f"DCE data shape: {dce_dataset.shape}")
print(f"VFA data shape: {vfa_data.shape}")
print(f"Time points: {n_timepoints}, spanning {time[-1]:.1f} seconds")

Data Organization

DCE-MRI data should be a 4D array with shape (x, y, z, time). VFA data should be a 4D array with shape (x, y, z, flip_angles).

Step 2: Compute T1 Map

Accurate T1 mapping is essential for converting signal to concentration:

# Set acquisition parameters on VFA dataset so compute_t1_map has what it needs
from osipy.dce import DCEAcquisitionParams
vfa_data.acquisition_params = DCEAcquisitionParams(
    tr=tr, flip_angles=flip_angles
)

# Compute T1 map using Variable Flip Angle method
t1_result = osipy.compute_t1_map(vfa_data, method="vfa")

# The result is a T1MappingResult with t1_map (ParameterMap) and quality_mask
t1_values = t1_result.t1_map.values
t1_quality = t1_result.quality_mask

print(f"T1 map shape: {t1_values.shape}")
print(f"T1 range: {t1_values[t1_quality > 0].min():.0f} - {t1_values[t1_quality > 0].max():.0f} ms")

T1 Quality Check

Typical tissue T1 values at 3T:

  • White matter: ~800-900 ms
  • Gray matter: ~1200-1400 ms
  • Blood: ~1600-1900 ms
  • CSF: ~4000+ ms

If your values are outside these ranges, check your flip angles and TR.

Step 3: Create a Tissue Mask

# Simple threshold-based mask
# Use the mean DCE signal to identify tissue (.data accesses the raw array)
mean_signal = dce_dataset.data.mean(axis=-1)
signal_threshold = mean_signal.max() * 0.1  # 10% of max

# Combine with T1 quality
mask = (mean_signal > signal_threshold) & (t1_quality > 0)

# Optionally: exclude CSF (very high T1)
mask = mask & (t1_values < 3000)

print(f"Tissue voxels: {mask.sum()} / {mask.size}")

Step 4: Convert Signal to Concentration

from osipy.dce import DCEAcquisitionParams

# Define acquisition parameters
acq_params = DCEAcquisitionParams(
    tr=5.0,            # ms
    te=2.0,            # ms
    flip_angles=[15],  # degrees (DCE flip angle)
    baseline_frames=5,
    relaxivity=4.5,    # mM^-1 s^-1 (typical for Gd-DTPA at 3T)
)

# Convert signal to concentration
# dce_signal_to_concentration wraps signal_to_concentration from osipy.dce
concentration = osipy.dce_signal_to_concentration(
    signal=dce_dataset.data,
    t1_map=t1_result.t1_map,
    acquisition_params=acq_params,
)

print(f"Concentration shape: {concentration.shape}")

# Check concentration range
valid_conc = concentration[mask]
print(f"Concentration range: {valid_conc.min():.3f} - {valid_conc.max():.3f} mM")

Negative Concentrations

If you see negative concentrations, check:

  1. Baseline selection (should be pre-contrast)
  2. T1 map accuracy
  3. Flip angle and TR values

Step 5: Select an Arterial Input Function

The AIF describes contrast agent concentration in blood plasma. You have several options:

# Parker AIF - widely used population average
aif = osipy.ParkerAIF()(time)

# Alternative: Georgiou AIF (different shape)
# aif = osipy.GeorgiouAIF()(time)

print(f"AIF peak: {aif.concentration.max():.2f} mM at {time[aif.concentration.argmax()]:.1f} s")

Option B: Measured AIF

# If you have manually selected an arterial ROI:
# aif_roi_indices = [...]  # Indices of arterial voxels
# aif_signal = dce_dataset.data[aif_roi_indices].mean(axis=0)

# Convert to concentration using the same acquisition parameters
# aif_conc = osipy.dce_signal_to_concentration(
#     signal=aif_signal,
#     t1_map=None,  # Pass None if no T1 map for AIF voxels
#     acquisition_params=acq_params,
#     t1_blood=1600,  # Blood T1 at 3T in ms
# )

# Create AIF object
# aif = osipy.ArterialInputFunction(time=time, concentration=aif_conc)

See How to Use a Custom AIF for detailed instructions.

Step 6: Fit the Pharmacokinetic Model

# Fit Extended Tofts model
result = osipy.fit_model(
    "extended_tofts",
    concentration=concentration,
    aif=aif,
    time=time,
    mask=mask
)

# Extract parameter maps (result is a DCEFitResult object)
ktrans = result.parameter_maps["Ktrans"].values
ve = result.parameter_maps["ve"].values
vp = result.parameter_maps["vp"].values
r_squared = result.r_squared_map
quality_mask = result.quality_mask

# Statistics on valid voxels
valid = quality_mask
print(f"\nFitting Results (valid voxels: {valid.sum()}):")
print(f"  Ktrans: {ktrans[valid].mean():.4f} ± {ktrans[valid].std():.4f} min⁻¹")
print(f"  ve:     {ve[valid].mean():.4f} ± {ve[valid].std():.4f}")
print(f"  vp:     {vp[valid].mean():.4f} ± {vp[valid].std():.4f}")
print(f"  R²:     {r_squared[valid].mean():.4f} ± {r_squared[valid].std():.4f}")

See Understanding Pharmacokinetic Models for model comparison and selection guidance.

Compare with Standard Tofts

# Compare with Standard Tofts
result_standard = osipy.fit_model(
    "tofts",
    concentration=concentration,
    aif=aif,
    time=time,
    mask=mask
)

print(f"\nStandard Tofts R²: {result_standard.r_squared_map[valid].mean():.4f}")
print(f"Extended Tofts R²: {r_squared[valid].mean():.4f}")

Step 7: Visualize Results

Parameter Maps

import matplotlib.pyplot as plt

# Select a representative slice
slice_idx = dce_dataset.shape[2] // 2

fig, axes = plt.subplots(2, 3, figsize=(12, 8))

# Row 1: Parameter maps
im0 = axes[0, 0].imshow(ktrans[:, :, slice_idx], cmap='hot', vmin=0, vmax=0.5)
axes[0, 0].set_title('Ktrans (min⁻¹)')
plt.colorbar(im0, ax=axes[0, 0])

im1 = axes[0, 1].imshow(ve[:, :, slice_idx], cmap='viridis', vmin=0, vmax=1)
axes[0, 1].set_title('ve')
plt.colorbar(im1, ax=axes[0, 1])

im2 = axes[0, 2].imshow(vp[:, :, slice_idx], cmap='plasma', vmin=0, vmax=0.2)
axes[0, 2].set_title('vp')
plt.colorbar(im2, ax=axes[0, 2])

# Row 2: Quality metrics
im3 = axes[1, 0].imshow(r_squared[:, :, slice_idx], cmap='gray', vmin=0, vmax=1)
axes[1, 0].set_title('R²')
plt.colorbar(im3, ax=axes[1, 0])

im4 = axes[1, 1].imshow(quality_mask[:, :, slice_idx], cmap='binary')
axes[1, 1].set_title('Quality Mask')
plt.colorbar(im4, ax=axes[1, 1])

# Anatomy reference
im5 = axes[1, 2].imshow(mean_signal[:, :, slice_idx], cmap='gray')
axes[1, 2].set_title('Mean Signal')
plt.colorbar(im5, ax=axes[1, 2])

for ax in axes.flat:
    ax.axis('off')

plt.tight_layout()
plt.savefig('dce_parameter_maps.png', dpi=300, bbox_inches='tight')
plt.show()

Visualize Time Curves

# Plot concentration time curves for sample voxels
fig, ax = plt.subplots(figsize=(10, 5))

# AIF
ax.plot(time, aif.concentration, 'r-', linewidth=2, label='AIF')

# Sample tissue voxels
sample_indices = np.where(valid)
n_samples = min(5, len(sample_indices[0]))
for i in range(n_samples):
    x, y, z = sample_indices[0][i*10], sample_indices[1][i*10], sample_indices[2][i*10]
    ax.plot(time, concentration[x, y, z, :], alpha=0.5, label=f'Voxel {i+1}')

ax.set_xlabel('Time (s)')
ax.set_ylabel('Concentration (mM)')
ax.set_title('DCE-MRI Concentration Time Curves')
ax.legend()
ax.grid(True, alpha=0.3)

plt.savefig('dce_time_curves.png', dpi=150)
plt.show()

Step 8: Export Results

Save your results in BIDS-compliant format:

Experimental Feature

BIDS derivative export is partially implemented and may not produce fully compliant output for all use cases.

Export to BIDS derivatives

# Export to BIDS derivatives
# export_bids signature: (parameter_maps, output_dir, subject_id, session_id, metadata)
# parameter_maps must be dict[str, ParameterMap]
osipy.export_bids(
    parameter_maps=result.parameter_maps,
    output_dir="derivatives/osipy",
    subject_id="01",
    session_id="01",
)

print("Results exported to derivatives/osipy/")

The output structure will be:

derivatives/osipy/
└── sub-01/
    └── ses-01/
        └── perf/
            ├── sub-01_ses-01_Ktrans.nii.gz
            ├── sub-01_ses-01_ve.nii.gz
            ├── sub-01_ses-01_vp.nii.gz
            ├── sub-01_ses-01_R2.nii.gz
            ├── sub-01_ses-01_quality-mask.nii.gz
            └── sub-01_ses-01_dce.json  # Provenance metadata

Complete Example

import numpy as np
import osipy
from osipy.dce import DCEAcquisitionParams

# 1. Load data (returns PerfusionDataset objects)
dce_dataset = osipy.load_nifti("dce_4d.nii.gz")
vfa_data = osipy.load_nifti("vfa_images.nii.gz")

# 2. Acquisition parameters
acq_params = DCEAcquisitionParams(
    tr=5.0, te=2.0, flip_angles=[2, 5, 10, 15],
    baseline_frames=5, relaxivity=4.5,
)
temporal_resolution = 3.5

# 3. Time array
n_timepoints = dce_dataset.shape[-1]
time = np.arange(n_timepoints) * temporal_resolution

# 4. T1 mapping (requires DCEAcquisitionParams on the dataset)
vfa_data.acquisition_params = DCEAcquisitionParams(
    tr=acq_params.tr, flip_angles=acq_params.flip_angles,
)
t1_result = osipy.compute_t1_map(vfa_data, method="vfa")

# 5. Create mask
mean_signal = dce_dataset.data.mean(axis=-1)
mask = (mean_signal > mean_signal.max() * 0.1) & (t1_result.quality_mask > 0)

# 6. Signal to concentration
concentration = osipy.dce_signal_to_concentration(
    dce_dataset.data, t1_result.t1_map, acq_params
)

# 7. AIF
aif = osipy.ParkerAIF()(time)

# 8. Fit model
result = osipy.fit_model("extended_tofts", concentration, aif, time, mask=mask)

# 9. Export (expects dict[str, ParameterMap])
osipy.export_bids(result.parameter_maps, "derivatives/osipy", "01", "01")

Next Steps

Troubleshooting

Low R² Values

  • Check T1 map quality (typical R² > 0.9 for VFA fitting)
  • Verify AIF timing matches bolus arrival — try fit_delay=True to estimate arterial delay automatically
  • Consider simpler model (Standard Tofts)

Unrealistic Ktrans Values

  • Ktrans > 1 min⁻¹ often indicates fitting failure
  • Check time units (osipy expects seconds)
  • Verify relaxivity value for your contrast agent

Memory Issues

  • Use masking to reduce voxel count
  • Process in chunks for very large datasets
  • Consider GPU acceleration for large-scale fitting