Computation of frequency- and time-domain Jacobians in optical tomography with Monte Carlo simulations

This paper presents a complete theoretical framework and open-source Monte Carlo implementation for computing frequency- and time-domain Jacobians in optical tomography, demonstrating their necessity for accurate modeling in low-scattering regimes and the benefits of realistic detector modeling for short source-detector separations.

The Gist

The Big Picture: Mapping the Brain with Light

Imagine trying to see inside a thick, foggy forest (your brain) using only a flashlight. You can't see the trees clearly because the fog scatters the light in every direction. This is what happens when scientists try to image the human brain using near-infrared light. The brain is full of "fog" (tissue) that bounces light around.

To figure out what's happening inside—like if a part of the brain is active or if there is a tumor—scientists use a technique called Optical Tomography. They shine light in at one spot and measure how much light comes out at another spot. By doing this many times, they try to build a 3D map of the brain's insides.

The Problem: The "Gold Standard" is Slow and Incomplete

To make this map accurate, scientists need a mathematical guide called a Jacobian. Think of a Jacobian as a "sensitivity map." It answers the question: "If I change the fog density in this tiny spot, how much will the light coming out at the detector change?"

For a long time, the most accurate way to calculate these maps was using Monte Carlo (MC) simulations. This is like running a massive video game where you simulate billions of individual photons (particles of light) bouncing around the brain to see where they end up. It's the "gold standard" because it's incredibly accurate.

However, there were two big gaps in this method:

1. Missing Tools: While scientists could simulate simple light measurements, they couldn't easily simulate more advanced measurements (like light that oscillates at a specific radio frequency or light that arrives at different times) using this gold-standard method.
2. The "Foggy" Shortcut: Because simulating billions of photons takes a supercomputer a long time, many scientists use a shortcut called the Diffusion Approximation (DA). This is like assuming the fog is perfectly uniform and smooth. It's fast, but it breaks down in "clear" spots in the brain (like the fluid-filled spaces around the brain) where the light doesn't behave like smooth fog.

What This Paper Did

The authors, working with a powerful software called MCX (Monte Carlo eXtreme), did three main things:

1. Built New Tools for the Simulation

They wrote new mathematical formulas to let the simulation calculate Jacobians for Frequency-Domain (light that wiggles like a radio wave) and Time-Domain (light that arrives in a specific time sequence) measurements.

• The Analogy: Imagine you were previously only able to count how many raindrops hit a bucket. Now, they gave you tools to also measure the speed of the raindrops and the pitch of the sound they make when they hit. This gives you much more information about the storm.

2. Created a "Realistic" Detector

In many simulations, the detector is treated like a magical black hole that catches any light hitting a specific circle on the skin. In reality, the detectors are fiber-optic cables with glass prisms that only catch light coming from specific angles.

• The Analogy: Imagine trying to catch rain with a bucket.
  • Old Model: The bucket is a giant, wide funnel that catches rain from any angle.
  • New Model: The bucket is a narrow straw. It only catches rain falling straight down.
  • The Result: The authors added a "post-processing" step to their simulation. After the light hits the skin, they check: "Did this photon hit the straw at the right angle?" If not, they discard it. They found this changes the sensitivity map, especially for short distances between the light source and the detector.

3. Proved the Shortcut is Flawed in "Clear" Areas

They compared their new, super-accurate Monte Carlo maps against the "shortcut" (Diffusion Approximation) maps using models of newborn babies' heads.

• The Finding: In areas where the brain is very "foggy" (high scattering), the shortcut works great. But in areas with Cerebrospinal Fluid (CSF)—which is like clear water compared to fog—the shortcut fails. It predicts that the light is much more sensitive to changes than it actually is.
• The Takeaway: If you are studying the brain, you cannot trust the shortcut near the fluid-filled spaces. You need the heavy-duty Monte Carlo simulation to get the right answer.

Why This Matters (According to the Paper)

• Better Maps: By using these new formulas, scientists can now build more accurate 3D maps of the brain, especially for newborns who have different brain structures than adults.
• Short Distances: For measurements taken very close together (short distances), the realistic detector model (the "straw" vs. the "funnel") matters. It reduces the sensitivity to the very surface of the skin and slightly increases the sensitivity to the deeper brain tissue.
• Validation: The paper proves that when you remove the "clear fluid" from the model, the fast shortcut matches the slow, accurate simulation. This confirms that the difference they saw earlier was indeed caused by the fluid, not a mistake in their math.

Summary

The authors upgraded the "gold standard" simulation software to handle more complex types of light measurements and added a realistic model for how the detector "sees" the light. They proved that while fast shortcuts work well in thick fog, they fail in clear fluid, and that realistic detector models are crucial for getting accurate readings, especially when the light source and detector are close together.

Technical Summary

1. Problem Statement

Optical Tomography (OT) is a non-invasive imaging technique used to reconstruct the spatial distribution of optical properties (absorption $\mu_a$ and scattering $\mu_s$) in biological tissues. Accurate image reconstruction relies on solving an ill-posed inverse problem, which requires a forward model and a Jacobian matrix (sensitivity profile) describing how measurements change with respect to tissue parameters.

While Monte Carlo (MC) simulations are the "gold standard" for modeling light transport (solving the Radiative Transfer Equation, RTE) due to their accuracy in complex geometries and heterogeneous media, the methodology for computing Jacobians in Frequency-Domain (FD) and Time-Domain (TD) modalities has been incomplete. Specifically:

• Existing MC methods often lack direct formulations for FD scattering Jacobians and TD Jacobians (intensity and mean time-of-flight).
• The Diffusion Approximation (DA), commonly used for speed, fails in low-scattering regions (e.g., cerebrospinal fluid in neonatal heads) and near sources/detectors.
• Standard MC detector models often assume isotropic reception, ignoring the physical constraints of real-world probes (e.g., prism-terminated optical fibers and Numerical Aperture limitations), leading to inaccuracies in short source-detector separations (SDS).

2. Methodology

The authors developed a comprehensive theoretical framework and implemented it in the high-performance, GPU-accelerated Monte Carlo eXtreme (MCX) simulator.

A. Theoretical Derivations (Perturbation MC)

Using the Microscopic Beer–Lambert Law (mBLL) and the Perturbation Monte Carlo (pMC) "replay" mode, the authors derived analytical formulas for Jacobians:

• Frequency-Domain (FD): Derived formulas for the complex intensity ($I_{FD}$), log-amplitude, and phase shift with respect to both $\mu_a$ and $\mu_s$. This involves weighting photon paths by $e^{i\omega t}$.
• Time-Domain (TD): Derived formulas for time-integrated intensity and Mean Time-of-Flight (TOF) with respect to $\mu_a$ and $\mu_s$. This utilizes the first temporal moment of the Time Point Spread Function (TPSF).
• Generalization: The formulas were generalized to handle split-voxel models (hybrid mesh-voxel) to accurately model curved boundaries and prism terminals.

B. Implementation in MCX

The new Jacobian types were implemented in the MCXLAB (v2025.10) environment. The computation follows a two-step "replay" approach:

1. Forward Simulation: A baseline simulation records photon trajectories that reach the detector.
2. Replay Simulation: The same trajectories are replayed to accumulate sensitivity metrics (path lengths, collision counts, TOF) without re-sampling random numbers.

• New Output Metrics: The authors introduced three new replay output types to support these calculations:
  • rfmus: Weighted scattering counts.
  • wltof: TOF-weighted path lengths.
  • wptof: TOF-weighted scattering counts.

C. Realistic Detector Modeling

A post-processing function was developed to model prism-terminated optical fiber detectors (common in Aalto University's FD instruments).

• Criteria: Photon packets are accepted only if they exit the tissue within a specific area, refract through a glass prism, and enter the fiber bundle within the Numerical Aperture (NA) limits.
• Split-Voxel Boundary (SVMC): To handle the curvature of the head surface under the prism, the authors utilized the SVMC approach, where boundary voxels are split into sub-volumes with tilted planes, allowing for more accurate modeling of photon exit angles.

D. Validation

• Internal Consistency: Validated against analytical linearization of FD and TD relationships (e.g., phase $\approx$ mean TOF $\times$ frequency).
• External Validation: Compared MC-derived Jacobians against a Finite Element Method (FEM) solver using the Diffusion Approximation (DA) in neonatal head models.
• Test Cases: Simulations were performed on full-term and preterm neonatal head models (segmented into scalp/skull, CSF, grey matter, white matter) at 798 nm.

3. Key Contributions

1. Complete Jacobian Framework: Provided the first complete theoretical derivation and MC implementation for FD scattering Jacobians and TD mean TOF Jacobians (both absorption and scattering).
2. Realistic Detector Modeling: Introduced a rigorous post-processing method to simulate prism-terminated fiber detectors with NA limitations, moving beyond simple isotropic reception models.
3. Split-Voxel Jacobians: Extended the pMC "replay" mode to work with Split-Voxel Monte Carlo (SVMC) for accurate sensitivity profiling on curved surfaces.
4. Open-Source Software: Integrated these capabilities into the open-source MCX simulator, making them available to the biophotonics community.

4. Results

• MC vs. DA Agreement:
  • In high-scattering regimes (e.g., grey/white matter), MC and DA Jacobians showed excellent agreement.
  • In low-scattering regimes (specifically the CSF layer), DA significantly overestimated sensitivity compared to MC. This confirms that DA is invalid for CSF-rich regions (common in neonates), necessitating MC for accurate reconstruction.
• FD vs. TD Consistency: The FD phase-shift Jacobians and TD mean TOF Jacobians showed high consistency, validating the theoretical link between the two modalities at typical modulation frequencies (100 MHz).
• Impact of Detector Modeling:
  • The realistic prism/NA model rejected ~98% of photon packets that would be accepted by a simple isotropic model.
  • Sensitivity Shift: The realistic model reduced sensitivity to the most superficial tissues (scalp) and slightly increased sensitivity to deeper brain tissue at short source-detector separations (< 2 cm). This is crucial for separating extracerebral noise from brain signals.
• Convergence: Scattering Jacobians required significantly more photon packets ($10^{12}$–$10^{13}$) than absorption Jacobians to achieve visual smoothness due to the subtraction of similar-magnitude terms in the derivative formula.

5. Significance

• Improved Reconstruction Accuracy: By providing accurate Jacobians for FD/TD data and realistic detector models, this work enables more precise image reconstruction, particularly for neonatal brain imaging where CSF layers are prominent and DA fails.
• Short-Channel Optimization: The study highlights that realistic detector modeling is essential for short-separation channels (< 2 cm), which are critical for removing superficial physiological noise in functional brain imaging.
• Versatility: The framework supports complex geometries and heterogeneous tissues, making it applicable to various biomedical contexts (e.g., breast imaging, muscle imaging) and non-biological turbid media.
• Future Directions: The authors suggest using DA for fast initial estimates and MC for refinement, particularly for short channels and non-diffusive regions. The availability of scattering Jacobians opens new avenues for reconstructing scattering changes (e.g., related to neuronal activation or cell swelling) alongside absorption changes.

In summary, this paper bridges a critical gap in optical tomography by providing a robust, high-performance computational tool for generating accurate sensitivity profiles across all major measurement modalities (CW, FD, TD) while accounting for the physical realities of modern imaging probes.