Embedded network 1D-3D model (tissue perfusion)

We solve a tracer transport problem in a domain that consist of a porous medium block that has an embedded transport network. In this example the porous medium is brain tissue and the embedded network is the vasculature (blood vessels).

Table of contents. This description is structured as follows:

• Embedded network 1D-3D model (tissue perfusion)
  • Problem set-up
  • Network data and blood flow
  • Mathematical model
• Implementation
• Output
  • Folder layout and files
  • Part 1: Problem definition
  • Part 2: Spatial parameters
  • Part 3: Properties
  • Part 4: Main program

Problem set-up

In this example we simulate clearance of a substance present in the tissue through the blood. The tissue cube is assigned no-flux/symmetry boundary conditions assuming that identical cubes are mirrored on all sides. Therefore the tracer has to cross the vessel wall into the network (vessel lumen). In then get transported in the blood stream by advection and diffusion. In the network, inlet tracer mole fraction is zero and at the outlet the mole fraction gradient is zero, making the tracer being transported out by advection only. We write out VTK output and the total tracer concentration in the tissue as text file, in every time step.

Network data and blood flow

The domain consists of a small blood vessel network embedded in a porous tissue (pore space is the pore space of the extra-cellular matrix). The network is extracted based on the mouse cortical brain data from Blinder (2013). With the boundary condition estimated by Schmid (2017a) and the pressure and network data available from Schmid (2017b), blood flow simulations in DuMu^x give a pressure field in the entire network. A small part (200µm)³ has been extracted for this example. The example also contains a blood flow solver which uses the pressure boundary data to compute blood volume flux for every facet (vertex) in the network. Note: The blood flow solver is currently not documented in detail. This might be added in the future.

Mathematical model

We solve the following coupled, mixed-dimensional PDE system:

where the subscript T and B denote the tissue and the network (blood flow) compartment,

For the coupling, we use a formulation with line source and a perimeter integral operator, that is

All equations are discretized with a cell-centered finite-volume scheme (with two-point flux approximation) as spatial discretization method with tracer mole fraction as primary variables. The equations are based on DuMu^x's Tracer model, an advection-diffusion-reaction model for porous media. In time, we use an implicit Euler scheme. The arising linear system is solved with a stabilized bi-conjugate gradient solver (BiCGSTAB) with a block-diagonal zero-fill incomplete LU factorization (ILU0) preconditioner. For details on the spatial discretization scheme, we recommend the DuMu^x handbook or the DuMu^x paper.

Implementation

Below is an overview over the files in the example folder and links to the more detailed descriptions of some individual files. The main program flow can be found in the main function implemented in main.cc. This is the only source file that includes all the supplementary header files which contain important parts of the implementation. The CMakeLists.txt file instructs CMake to create an executable from main.cc and CMake will configure and figure out the necessary compiler command for compilation. (The executable is built in the build folder by typing make example_embedded_network_1d3d. See dumux-course basic exercise for more detailed instructions.)

Key classes that are implemented by the user are the problems (problem.hh), one for the tissue and one for the network. The problems define initial and boundary conditions. They also implement source terms. As the coupling between tissue and network is realized through point sources at integration points, the coupling condition is implemented in the problem classes. The problem classes are used by DuMu^x during assembly of the system matrix and the residual (discrete equations). For this to work, we have to tell the model which problem class to use. This is achieved through the property system in properties.hh where the Problem property is specialized for the main models created in this example (NetworkTransportModel and TissueTransportModel). These models are passed in the main function to the assembler.

Secondly, the spatial parameters (spatialparams.hh) are classes that specify (possibly) spatially varying parameter. On such parameter is the radius field for the network. In the class NetworkSpatialParams, the radius field is read from the grid file network.dgf (which is in the very simple, human-readable Dune Grid Format). As for the problem, spatial parameters have to be added to the model by specializing the SpatialParams property for the model in properties.hh.

Apart from problem and spatial params, the model (properties.hh) also has other configurable parameters. (In fact most of the inner workings of the assembler can be configured like this.) One example is the grid manager used for each model. Dune provides specialized implementations for certain grid types behind a common interface. In this exercise, we use a structured Cartesian grid (YaspGrid) for the tissue domain and an embedded network grid manager (FoamGrid) for the network. With the model configuration through the property system in mind, we can better understand the main program (main.cc) and how the boundary conditions and parameter setting make their way into the assembler.

Output

The simulation output VTK files at every time step that can be inspected in ParaView. It also writes a simple text file clearance_tracer_amounts.dat containing the total tracer amount at each time step. There is a small Python script plot.py, that visualizes this data. (Execute ./plot.py or python3 plot.py.)

Folder layout and files

└── embedded_network_1d3d/
    ├── CMakeLists.txt          -> build system file
    ├── main.cc                 -> main program flow
    ├── params.input            -> runtime parameter configuration file
    ├── plot.py                 -> small Python script to plot the tracer amount in the tissue
    ├── problem.hh              -> boundary & initial conditions and source terms
    ├── properties.hh           -> compile time model configuration
    ├── solver.hh               -> specialized solver for linear multi-domain problems
    ├── spatialparams.hh        -> spatially varying model parameters
    └── tracerfluidsystem.hh    -> fluid system with a single tracer

Part 1: Problem definition

Part 2: Spatial parameters

Part 3: Properties

Part 4: Main program