---
title: 'SimFA-python: Finite Element MRI Framework'
url: https://www.emergentmind.com/topics/simfa-python
type: topic
---

# SimFA-python: Finite Element MRI Framework

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 [1908.01719]. 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 [1908.01719].

## 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 [1908.01719].

For permeable membranes, the interface conditions on $\Gamma$ are

$$
\llbracket D \nabla U\cdot n^0\rrbracket = 0,\qquad \{D \nabla U\cdot n^0\} = -\kappa \llbracket U\rrbracket,
$$

with $\kappa$ denoting membrane permeability [1908.01719]. 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 $\mathbb{R}^3$:

$$
U_m = U_s e^{i\theta(t)},\qquad
D_m \nabla U_m\cdot n = (D_s \nabla U_s\cdot n)e^{i\theta(t)},
$$

with

$$
\theta(t)=\gamma (b_k-a_k) g_k \int_0^t f(s)\,ds
$$

for each periodic direction $k$ [1908.01719]. 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 $\Phi_h\in\{0,1\}$ is introduced, and the weak problem is written in the product space $V_h = V_{0h}\times V_{1h}$. The formulation is

$$
\langle \partial_t U_h, v\rangle_{\Omega_0\cup\Omega_1}
= -(i\gamma f(t) g\cdot x\, U_h + U_h/T_2, v)_\Omega
- (D\nabla U_h,\nabla v)_\Omega
+ \langle -\kappa \llbracket U_h\rrbracket,\llbracket v\rrbracket\rangle_\Gamma,
$$

for all test functions $v\in V_h$ [1908.01719]. 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 $\theta$-method on a partition $0=t^0<t^1<\dots<t^N=T$ with time steps $k^n=t^n-t^{n-1}$:

$$
\left(\frac{U_h^n-U_h^{n-1}}{k^n}, v\right)
=
\theta F(U_h^n,v,t^n)
+
(1-\theta)F(U_h^{n-1},v,t^{n-1}).
$$

The framework identifies $\theta=\tfrac{1}{2}$ with the second-order Crank–Nicolson scheme and $\theta=1$ with the backward Euler scheme [1908.01719]. 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` [1908.01719]. The library is complemented by auxiliary routines, batch scripts, and notebooks.

| Component | Function |
|---|---|
| `MRI_parameters` | Manages timing parameters $(\Delta,\delta,T)$, conversion between gradient strength $g$ and b-value, and symbolic definition of $f(t)$ via SymPy |
| `MRI_domain` | Wraps a DOLFIN mesh, function spaces, sub-domain and interface markers, diffusion tensor, $T_2$ field, and pseudo-periodic BCs |
| `MRI_simulation` | Sets initial conditions, solver parameters, time-stepping, executes `solve()`, and computes normalized signal $S(b)$ in `PostProcessing()` |
| `CreatePhaseFunc` | Generates the PUFEM phase function $\Phi_h$ |
| `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 $S(b)$ |

`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 $f(t)$, $\int f$, and the b–g conversion [1908.01719]. This makes arbitrary waveform specification part of the core parameter layer rather than an external precomputation.

`MRI_domain` wraps a DOLFIN mesh, continuous $P_1$ and discontinuous Galerkin function spaces, sub-domain and interface markers, and methods for imposing diffusion tensors, $T_2$ relaxation fields, and strong or weak pseudo-periodic boundary conditions [1908.01719]. 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(T)$ over $\Omega$ to compute the normalized signal $S(b)$ and exports `vtk` files for Paraview [1908.01719]. 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 [1908.01719]. 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 [1908.01719]. For local notebook-based work, the recommended mode is a Docker container based on `quay.io/fenicsproject/stable:${fenics_tag}` with a Jupyter server configured for browser access and Google Colab interoperability.

For HPC deployment, the documented workflow is Singularity plus MPI. A Singularity image is built from the Docker recipe, the preprocessing stage generates an `h5` input file, and `GCloudDmriSolver.py` is then launched with `mpirun` using specified parameters such as b-value, pulse timings, permeability, and gradient direction [1908.01719]. This deployment model is consistent with the framework’s description as portable and cloud-ready.

The intended usage pattern is correspondingly layered. Geometry and partition markers are prepared first; then `MRI_parameters` defines the sequence, `MRI_domain` imposes compartment-wise diffusion and relaxation properties, `MRI_simulation` advances the solution in time, and `PostProcessing()` returns $S(b)$ while exporting field data [1908.01719]. 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 [1908.01719]. In the validation examples, three-layered disk, sphere, and torus configurations showed agreement within 1%, and a cylinder example with $T_2=[\infty,40\ \mathrm{ms},40\ \mathrm{ms}]$ reproduced known signal attenuation [1908.01719].

A major capability is support for arbitrary waveforms. PGSE, OGSE, and double or trapezoidal pulses are all listed as matching reference signals [1908.01719]. This is a direct consequence of the symbolic `Piecewise` handling of $f(t)$ in `MRI_parameters`, which allows the waveform to be encoded at the parameter-definition stage rather than hard-coded into the PDE solver.

Performance results are given for a single neuron mesh of approximately $0.6$ million vertices and $2.5$ million cells with $\Delta t = 200\ \mathrm{ms}$. 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 [1908.01719]. 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 [1908.01719].

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 [1908.01719]. 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 $\Delta t=100$ instead of $\Delta t=10$ for artificial permeability, and that Crank–Nicolson should be used for second-order accuracy [1908.01719].

## 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 [1908.01719]. 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 [1908.01719].

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` [1810.00674]. That description emphasizes declarative and imperative APIs, approximately 120 predefined weak-form terms, and two-scale homogenization workflows rather than diffusion MRI [1810.00674].

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 [1908.01719]. In contrast, the SfePy description concerns a general finite-element package with a broader PDE and homogenization scope [1810.00674]. 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.

Source: https://www.emergentmind.com/topics/simfa-python