SimFA-python: Finite Element MRI Framework

• SimFA-python is a portable, open-source finite-element framework that simulates diffusion MRI by solving the Bloch–Torrey equation with support for multi-compartment media and permeable interfaces.
• It integrates advanced finite element techniques, including Partition-of-Unity methods and pseudo-periodic boundary conditions, to handle arbitrary diffusion-encoding waveforms.
• It offers container-ready workflows for local, cloud, and HPC deployments, enabling reproducible research and scalable parallel execution.

SimFA-python is a portable finite-element framework for diffusion MRI simulation that solves the Bloch–Torrey equation in Python, builds on FEniCS for mesh handling, variational forms, and solvers, and adds diffusion-MRI-specific support for multi-compartment media, permeable interfaces, arbitrary diffusion-encoding waveforms, and pseudo-periodic boundary conditions (Nguyen et al., 2019). It is described as fully open-source and container-ready, with a workflow spanning a lightweight Python library, command-line pre-/post-processing scripts, and high-level Jupyter notebooks for Google Colaboratory or Singularity/Docker environments. Within computational diffusion MRI, its stated purpose is to reduce the gap between finite-element software capabilities and the simulation requirements of the MRI community, while supporting reproducible science, cloud execution, and MPI parallelization (Nguyen et al., 2019).

1. Mathematical model for diffusion MRI

SimFA-python is organized around the Bloch–Torrey equation in a domain $\Omega$ divided into two sub-domains $\Omega_0,\Omega_1$ with interface $\Gamma$. In strong form, the transverse magnetization $U(x,t)$ satisfies

$$\partial_t U = -i\gamma f(t) g\cdot x \, U - U/T_2(x) + \nabla\cdot(D(x)\nabla U),$$

where $\gamma$ is the proton gyromagnetic ratio, $f(t)$ is the temporal waveform of the diffusion-encoding gradient, $g$ is the gradient vector, $D(x)\in\mathbb{R}^{3\times 3}$ is the piecewise-constant diffusion tensor, and $T_2(x)$ is the local transverse relaxation time (Nguyen et al., 2019).

For permeable membranes, the interface conditions on $\Omega_0,\Omega_1$0 are

$\Omega_0,\Omega_1$1

with $\Omega_0,\Omega_1$2 denoting membrane permeability (Nguyen et al., 2019). This interface treatment is central to the framework’s support for multi-compartment media.

The framework also implements pseudo-periodic boundary conditions for domains considered tiled in $\Omega_0,\Omega_1$3:

$\Omega_0,\Omega_1$4

with

$\Omega_0,\Omega_1$5

for each periodic direction $\Omega_0,\Omega_1$6 (Nguyen et al., 2019). In diffusion MRI terms, this allows periodic microstructural models to be coupled to general gradient waveforms without abandoning the finite-element formulation.

2. Variational formulation and time integration

To impose permeability, SimFA-python uses the Partition-of-Unity FEM (PUFEM). A two-phase indicator $\Omega_0,\Omega_1$7 is introduced, and the weak problem is written in the product space $\Omega_0,\Omega_1$8. The formulation is

$\Omega_0,\Omega_1$9

for all test functions $\Gamma$0 (Nguyen et al., 2019). The role of this form is to encode diffusion, relaxation, gradient-induced phase accumulation, and membrane exchange within a single variational statement.

Time discretization is performed with the implicit $\Gamma$1-method on a partition $\Gamma$2 with time steps $\Gamma$3:

$\Gamma$4

The framework identifies $\Gamma$5 with the second-order Crank–Nicolson scheme and $\Gamma$6 with the backward Euler scheme (Nguyen et al., 2019). In practice, this places SimFA-python within the standard finite-element time-stepping tradition while retaining direct control over diffusion-MRI-specific source terms and interface operators.

3. Software architecture and principal abstractions

SimFA-python is centered on the Python file DmriFemLib.py, with three main classes encapsulating the workflow: MRI_parameters, MRI_domain, and MRI_simulation (Nguyen et al., 2019). The library is complemented by auxiliary routines, batch scripts, and notebooks.

Component	Function
MRI_parameters	Manages timing parameters $\Gamma$7, conversion between gradient strength $\Gamma$8 and b-value, and symbolic definition of $\Gamma$9 via SymPy
MRI_domain	Wraps a DOLFIN mesh, function spaces, sub-domain and interface markers, diffusion tensor, $U(x,t)$0 field, and pseudo-periodic BCs
MRI_simulation	Sets initial conditions, solver parameters, time-stepping, executes solve(), and computes normalized signal $U(x,t)$1 in PostProcessing()
CreatePhaseFunc	Generates the PUFEM phase function $U(x,t)$2
GetPartitionMarkers	Reads Gmsh “physical groups” into a DOLFIN MeshFunction
PreprocessingOneCompt.py / PreprocessingMultiCompt.py	Build and save h5 input for one- and multi-compartment domains
GCloudDmriSolver.py	Reads h5, unpacks domain and parameters, runs the simulation, and writes out $U(x,t)$3

MRI_parameters manages timing parameters, the conversion between gradient strength and b-value, and the symbolic definition of the encoding waveform through SymPy. Its member fs_sym holds the Piecewise symbolic expression, and calling mp.Apply() builds callables for $U(x,t)$4, $U(x,t)$5, and the b–g conversion (Nguyen et al., 2019). This makes arbitrary waveform specification part of the core parameter layer rather than an external precomputation.

MRI_domain wraps a DOLFIN mesh, continuous $U(x,t)$6 and discontinuous Galerkin function spaces, sub-domain and interface markers, and methods for imposing diffusion tensors, $U(x,t)$7 relaxation fields, and strong or weak pseudo-periodic boundary conditions (Nguyen et al., 2019). It supports both single- and multi-compartment domains through IsDomainMultiple, and periodic versus non-periodic setups through IsDomainPeriodic.

MRI_simulation handles the initial condition, linear solver parameters, and time-stepping. Its method solve(domain, parameters, linsolver, initial_condition) executes the time loop and stores snapshots, while PostProcessing() integrates $U(x,t)$8 over $U(x,t)$9 to compute the normalized signal $$\partial_t U = -i\gamma f(t) g\cdot x \, U - U/T_2(x) + \nabla\cdot(D(x)\nabla U),$$0 and exports vtk files for Paraview (Nguyen et al., 2019). The architecture therefore separates sequence definition, geometric/physical domain definition, and numerical solution.

4. Execution workflow, containers, and computing environments

SimFA-python is explicitly designed for portability across local, cloud, and HPC environments (Nguyen et al., 2019). The prerequisites are Docker or Singularity, Python 3, and optionally conda.

For hosted execution, the framework supports Google Colaboratory, where FEniCS is installed through apt, after which DmriFemLib.py and the example notebooks can be loaded directly (Nguyen et al., 2019). For local notebook-based work, the recommended mode is a Docker container based on quay.io/fenicsproject/stable:$${fenics_tag}</code> with a Jupyter server configured for browser access and Google Colab interoperability.</p> <p>For HPC deployment, the documented workflow is Singularity plus MPI. A Singularity image is built from the Docker recipe, the preprocessing stage generates an <code>h5</code> input file, and <code>GCloudDmriSolver.py</code> is then launched with <code>mpirun</code> using specified parameters such as b-value, pulse timings, permeability, and gradient direction (<a href="/papers/1908.01719" title="" rel="nofollow" data-turbo="false" class="assistant-link" x-data x-tooltip.raw="">Nguyen et al., 2019</a>). This deployment model is consistent with the framework’s description as portable and cloud-ready.</p> <p>The intended usage pattern is correspondingly layered. Geometry and partition markers are prepared first; then <code>MRI_parameters</code> defines the sequence, <code>MRI_domain</code> imposes compartment-wise diffusion and relaxation properties, <code>MRI_simulation</code> advances the solution in time, and <code>PostProcessing()</code> returns$$\partial_t U = -i\gamma f(t) g\cdot x \, U - U/T_2(x) + \nabla\cdot(D(x)\nabla U),$1 while exporting field data (Nguyen et al., 2019). This suggests a workflow in which mesh generation, physics specification, and solver execution can be decoupled without leaving Python.

5. Validation, scaling, and representative problem classes

The framework is reported to have been verified against matrix-formalism references, compared to random-walk methods, and parallelized with MPI for HPC and cloud platforms (Nguyen et al., 2019). In the validation examples, three-layered disk, sphere, and torus configurations showed agreement within 1%, and a cylinder example with $\partial_t U = -i\gamma f(t) g\cdot x \, U - U/T_2(x) + \nabla\cdot(D(x)\nabla U),$2 reproduced known signal attenuation (Nguyen et al., 2019).

A major capability is support for arbitrary waveforms. PGSE, OGSE, and double or trapezoidal pulses are all listed as matching reference signals (Nguyen et al., 2019). This is a direct consequence of the symbolic Piecewise handling of $\partial_t U = -i\gamma f(t) g\cdot x \, U - U/T_2(x) + \nabla\cdot(D(x)\nabla U),$$3 in <code>MRI_parameters</code>, which allows the waveform to be encoded at the parameter-definition stage rather than hard-coded into the <a href="https://www.emergentmind.com/topics/prompt-driven-exploration-pde" title="" rel="nofollow" data-turbo="false" class="assistant-link" x-data x-tooltip.raw="">PDE</a> solver.</p> <p>Performance results are given for a single neuron mesh of approximately$$\partial_t U = -i\gamma f(t) g\cdot x \, U - U/T_2(x) + \nabla\cdot(D(x)\nabla U),$4 million vertices and $\partial_t U = -i\gamma f(t) g\cdot x \, U - U/T_2(x) + \nabla\cdot(D(x)\nabla U),$5 million cells with $\partial_t U = -i\gamma f(t) g\cdot x \, U - U/T_2(x) + \nabla\cdot(D(x)\nabla U),$6. On this problem, the reported scaling is near-ideal up to 32 cores on Tegner and good speedup up to 500 cores across 25 nodes (Nguyen et al., 2019). Reported wall-clock figures include roughly 30 minutes per b-value on Google Cloud using 8 cores and about 7 minutes on Tegner using 20 cores (Nguyen et al., 2019).

The framework also includes examples that extend beyond conventional full-3D compartment models. A 1D manifold discretization of thin dendrites is reported to accelerate computation by two orders of magnitude versus a full 3D mesh, and an extracellular-space mesh with 0.9 million tetrahedra computes three principal directions in roughly 30 minutes per b-value on Colab (Nguyen et al., 2019). In addition, the recommendations state that periodic meshes together with the transformed PDE in Equation (5) in the appendix should be used for pseudo-periodic boundary conditions to allow larger time steps, for example $\partial_t U = -i\gamma f(t) g\cdot x \, U - U/T_2(x) + \nabla\cdot(D(x)\nabla U),$7 instead of $\partial_t U = -i\gamma f(t) g\cdot x \, U - U/T_2(x) + \nabla\cdot(D(x)\nabla U),$8 for artificial permeability, and that Crank–Nicolson should be used for second-order accuracy (Nguyen et al., 2019).

6. Reproducibility, open-source positioning, and relation to adjacent Python FEM frameworks

SimFA-python is presented as an open-source framework intended to support reproducible science in computational diffusion MRI (Nguyen et al., 2019). The documented assets include mesh scripts, preprocessing and post-processing scripts, DmriFemLib.py, and example notebooks, all hosted on GitHub and runnable inside a single Docker or Singularity image. The stated consequence is end-to-end reproducibility: one can clone the repository, launch the container, and reproduce the figures, timings, and signals in the paper or extend the framework to new geometries and sequences (Nguyen et al., 2019).

Within the Python FEM landscape, the provided literature also contains a distinct description labeled “SimFA-python (SfePy),” referring to SfePy as an open-source, Python-based finite-element package with modules such as sfepy/mesh, sfepy/discrete, sfepy/terms, sfepy/solvers, sfepy/postproc, and sfepy/homogenization (Cimrman et al., 2018). That description emphasizes declarative and imperative APIs, approximately 120 predefined weak-form terms, and two-scale homogenization workflows rather than diffusion MRI (Cimrman et al., 2018).

This suggests a potential terminological ambiguity in secondary summaries. In the diffusion MRI literature, SimFA-python denotes the FEniCS-based framework organized around DmriFemLib.py, PUFEM interface handling, symbolic waveform specification, and diffusion-MRI post-processing (Nguyen et al., 2019). In contrast, the SfePy description concerns a general finite-element package with a broader PDE and homogenization scope (Cimrman et al., 2018). For diffusion MRI specifically, the defining characteristics of SimFA-python are therefore not merely its use of Python and finite elements, but its domain-specific treatment of Bloch–Torrey dynamics, permeable interfaces, pseudo-periodic boundary conditions, arbitrary encoding waveforms, and portable containerized execution.