Finite element discretization of a biological network formation system: a preliminary study

TL;DR

We address numerical discretization of the Cai-Hu biological network formation model, a coupled elliptic-parabolic system for pressure $p$ and conductivity tensor $\mathbb{C}$. The authors propose a low-order finite element discretization in space combined with a semi-implicit backward Euler time-stepping, implemented in FreeFEM, and validate the scheme via accuracy tests and energy decay observations. The work demonstrates that the method captures steady-state network patterns and reveals how parameters such as the background permeability $r$, diffusivity $D^2$, and metabolic exponents $\alpha,\gamma$ shape ramification and stability. The findings provide a practical numerical framework for exploring multi-scale, stiff network formation dynamics and guide future developments toward more precise and efficient solvers, including monolithic formulations, IMEX time integrators, adaptivity, and parallelization.

Abstract

A finite element discretization is developed for the Cai-Hu model, describing the formation of biological networks. The model consists of a non linear elliptic equation for the pressure $p$ and a non linear reaction-diffusion equation for the conductivity tensor $\mathbb{C}$. The problem requires high resolution due to the presence of multiple scales, the stiffness in all its components and the non linearities. We propose a low order finite element discretization in space coupled with a semi-implicit time advancing scheme. The code is {verified} with several numerical tests performed with various choices for the parameters involved in the system. In absence of the exact solution, we apply Richardson extrapolation technique to estimate the order of the method.

Figures (3)

• Figure 1: In this figure we show three results considering different values of the background permeability $r$. From left to right we have the results ofTest-r1 with $r = 10^{-1}$,Test-r2 with $r = 10^{-2}$ andTest-DD1with $r = 10^{-3}$. We plot two different quantities for each test: in the first row, the norm of the solution $\mathbb{C}$ at final time; in the second row, the time dependent plot of the energy $\mathcal{E}[\mathbb{C}]$.
• Figure 2: In this figure we show three results considering different values of the squared diffusivity $D^2$. From left to right we have the results of Test-DD1with $D^2 = 2\cdot 10^{-3}$,Test-DD2with $D^2 = 10^{-3}$ andTest-DD3with $D^2 = 5\cdot 10^{-4}$. Same format as Fig. \ref{['fig_testr']}.
• Figure 3: In this figure we show three results considering different values for $\alpha$ and $D^2$. From left to right we have the results ofTest-Al1with $\alpha = 0.75$ and $D^2 = 10^{-3}$,Test-Al2with $\alpha = 1.5$ and $D^2 = 5\cdot 10^{-4}$ andTest-DD3with $\alpha = 0.75$ and $D^2 = 5\cdot 10^{-4}$. Same format as Fig. \ref{['fig_testr']}.