Comparison of two aspects of a PDE model for biological network formation

In space we consider the square domain $\Omega=[0,1]\times[0,1]$, that we discretize by a uniform Cartesian mesh with spatial step $h:=\Delta x=\Delta y$. We call $\Omega_{h}$ the discrete computational domain. The two variables conductivity and the pressure are ${\mathbb{C}_{ij}\approx\mathbb{C}(x_{i},y_{j})},m_{ij}\approx m(x_{i},y_{j})$ and ${p_{ij}\approx p(x_{i},y_{j})}$, defined at the center of the cell $(i,j)$, therefore the set of grid points is $x_{i}=(i-1/2)h,\,y_{j}=(j-1/2)h,\,(i,j)\in\{1,\dots,N\}^{2}$, $hN=1$.

In order to obtain second order accuracy in space, we use central differences for the computation of the space derivatives [12]. Discretizing in space Eq. (9), we have:

	$\displaystyle\displaystyle\frac{\partial m^{(1)}}{\partial t}$	$\displaystyle=$	$\displaystyle D^{2}\mathcal{L}\,m^{(1)}+c^{2}\left(\mathcal{D}_{x}p\right)^{2}\,m^{(1)}+c^{2}\mathcal{D}_{x}p\mathcal{D}_{y}p\,m^{(2)}-\alpha|m|^{2(\gamma-1)}m^{(1)}$		(12)
	$\displaystyle\displaystyle\frac{\partial m^{(2)}}{\partial t}$	$\displaystyle=$	$\displaystyle D^{2}\mathcal{L}\,m^{(2)}+c^{2}\left(\mathcal{D}_{y}p\right)^{2}\,m^{(2)}+c^{2}\mathcal{D}_{x}p\mathcal{D}_{y}p\,m^{(1)}-\alpha|m|^{2(\gamma-1)}m^{(2)}$		(13)

while the semi-discrete version of Eq. (4) is

	$\displaystyle\displaystyle\frac{\partial{C}^{(1,1)}}{\partial t}$	$\displaystyle=$	$\displaystyle D^{2}\mathcal{L}\,{C}^{(1,1)}+c^{2}\left(\mathcal{D}_{x}p\right)^{2}-\gamma|\mathbb{C}|^{\gamma-2}{C}^{(1,1)}$		(14)
	$\displaystyle\displaystyle\frac{\partial{C}^{(1,2)}}{\partial t}$	$\displaystyle=$	$\displaystyle D^{2}\mathcal{L}\,{C}^{(1,2)}+c^{2}\mathcal{D}_{x}p\mathcal{D}_{y}p-\gamma|\mathbb{C}|^{\gamma-2}{C}^{(1,2)}$		(15)
	$\displaystyle\displaystyle\frac{\partial{C}^{(2,2)}}{\partial t}$	$\displaystyle=$	$\displaystyle D^{2}\mathcal{L}\,{C}^{(2,2)}+c^{2}\left(\mathcal{D}_{y}p\right)^{2}-\gamma|\mathbb{C}|^{\gamma-2}{C}^{(2,2)}$		(16)

where $\mathcal{L}$ is the discrete Laplacian operator and $\mathcal{D}_{x}$ and $\mathcal{D}_{y}$ are, respectively, the discrete $x$ and $y$ first derivative operators, both using central difference approximation. The norm $|\cdot|$ we use is the Frobenius one. In this case we have three equations instead of four, because $\mathbb{C}$ is symmetric.

In order to have a fully implicit scheme, we define a compact form of the semi-discrete Eqs. (12-13) and Eqs. (14-16), as follows:

	$\displaystyle\displaystyle\frac{\partial{m}_{\rm comp}}{\partial t}$	$\displaystyle=$	$\displaystyle D^{2}\mathcal{L}\,{m}_{\rm comp}+c^{2}\mathcal{P}\,{m}_{\rm comp}-\alpha\mathcal{Q}^{m}(m){m}_{\rm comp}$		(17)
	$\displaystyle\displaystyle\frac{\partial{\mathbb{C}}_{\rm comp}}{\partial t}$	$\displaystyle=$	$\displaystyle D^{2}\mathcal{L}\,{\mathbb{C}}_{\rm comp}+c^{2}\mathcal{P}-\alpha\mathcal{Q}^{c}(\mathbb{C}){\mathbb{C}}_{\rm comp}$		(18)

where ${m}_{\rm comp}=[m^{(1)},m^{(2)}]^{T}$, ${\mathbb{C}}_{\rm comp}=[C^{(1,1)},C^{(1,2)},C^{(2,2)}]^{T}$ and $\mathcal{P}=\mathcal{D}_{\beta}p\,\mathcal{D}_{\eta}p$, with the suitable choice of $\beta,\eta\in\{x,y\}$. For the metabolic terms, we have

	$\displaystyle\mathcal{Q}^{m}(m)$	$\displaystyle=$	$\displaystyle|m|^{2(\gamma-1)}$		(19)
	$\displaystyle\mathcal{Q}^{c}(\mathbb{C})$	$\displaystyle=$	$\displaystyle|\mathbb{C}|^{\gamma-2}.$		(20)

Here, we discretize the Poisson equation for the pressure. To obtain a conservative scheme, we consider the following discretization: first, we extend the formula in Eq. (3), that becomes

$\displaystyle\partial_{x}\left(\left(r+C^{(1,1)}\right)\partial_{x}p\right)+\partial_{x}\left(C^{(1,2)}\partial_{y}p\right)+\partial_{y}\left(C^{(1,2)}\partial_{x}p\right)+\partial_{y}\left(\left(r+C^{(2,2)}\right)\partial_{y}p\right)=-S$

(21)

where we use the symmetry $C^{(1,2)}=C^{(2,1)}$.

Now, we discretize the components of the formula, one by one, since we use different discretizations. For simplicity we pose $\mathcal{C}^{1,1}=r+C^{(1,1)}$:

	$\displaystyle\partial_{x}\left(\mathcal{C}^{1,1}\partial_{x}\,p\right)_{i,j}\approx$	$\displaystyle\frac{1}{\Delta x}\left(\mathcal{C}^{1,1}_{i+1/2,j}\partial_{x}\,p_{i+1/2,j}-\mathcal{C}^{1,1}_{i-1/2,j}\partial_{x}\,p_{i-1/2,j}\right)$
	$\displaystyle=$	$\displaystyle\frac{1}{2\Delta x^{2}}\left(\left(\mathcal{C}^{1,1}_{i+1,j}+\mathcal{C}^{1,1}_{i,j}\right)p_{i+1,j}+\left(\mathcal{C}^{1,1}_{i-1,j}+\mathcal{C}^{1,1}_{i,j}\right)p_{i-1,j}\right)$
		$\displaystyle-\frac{1}{2\Delta x^{2}}\left(\mathcal{C}^{1,1}_{i+1,j}+\mathcal{C}^{1,1}_{i-1,j}+2\mathcal{C}^{1,1}_{i,j}\right)p_{i,j}$		(22)

where in the last line we consider the following approximations:

$$\partial_{x}\,p_{i+1/2,j}\approx\frac{p_{i+1,j}-p_{i,j}}{\Delta x},\quad\mathcal{C}^{1,1}_{i+1/2,j}\approx\frac{\mathcal{C}^{1,1}_{i+1,j}+\mathcal{C}^{1,1}_{i,j}}{2}.$$

We omit the term with both $y$-derivatives because it is analogue to the one with $x$-derivatives. Now we discretize the term with mix derivatives.

	$\displaystyle\partial_{x}\left(C^{(1,2)}\partial_{y}\,p\right)_{i,j}\approx$	$\displaystyle\frac{1}{\Delta x}\left(C^{(1,2)}_{i+1/2,j}\partial_{y}p_{i+1/2,j}-C^{(1,2)}_{i-1/2,j}\partial_{y}p_{i-1/2,j}\right)$		(23)
	$\displaystyle=$	$\displaystyle\frac{1}{8\Delta x^{2}}\left(C^{(1,2)}_{i+1,j}+C^{(1,2)}_{i,j}\right)(p_{i+1,j+1}-p_{i+1,j-1})$
		$\displaystyle-\frac{1}{8\Delta x^{2}}\left(C^{(1,2)}_{i-1,j}+C^{(1,2)}_{i,j}\right)(p_{i-1,j+1}-p_{i-1,j-1})$		(24)
		$\displaystyle+\frac{1}{8\Delta x^{2}}\left(C^{(1,2)}_{i+1,j}-C^{(1,2)}_{i-1,j}\right)(p_{i,j+1}-p_{i,j-1})$		(25)

again, we omit the term with $y,x$-derivatives because it is analogue to the one with $x,y$-derivatives.

At this stage, we describe the time discretization that we apply to the model. Since systems (8-6) and (3-4) are stiff in all their components, the choice of the time discretization is crucial for the efficiency.

We also need a high performing scheme in time, since at each time step we compute the solution of seven linear systems (two for the conductivity vector $m$, three for the conductivity tensor $\mathbb{C}$ and two Poisson equations for the pressure $p$).

As we shall see, for some values of the parameters, the well-posedness of the problem becomes weaker, which reflects the bad conditioning of the numerical problem. For such a reason, we adopt a symmetric scheme, which better preserves possible symmetries of the solution. In particular, we adopt a symmetric-ADI scheme for both the conductivity variables, which guarantees efficiency, second order accuracy and spatial symmetry.

Anyway, the scheme is not strictly second order accurate in time for two reasons. First, the pressure is computed at time $n$ rather than at an intermediate time $n+1/2$. Second, the metabolic term is treated partially explicitly and partially implicitly, thus destroying second order accuracy. Improvements of the order of accuracy in time are currently under investigation.

In Table 1 we define the tests we want to show in this paper, varying the parameters of the systems. This choice of parameters, a typical time scale is of the order of unit, while after time 15, the solution reached the steady state.

First, we check the accuracy of the schemes adopted. In Table 2 and Table 3 we see the error for the conductivity variables, calculated with Richardson extrapolation (see, e.g., [13]). We show the error for the module of the vector and of the tensor, and the parameters chosen are defined in TestA, TestB and TestC.

Here we define the initial conditions $m^{0}_{\rm comp}(\vec{x})$ and $\mathbb{C}^{0}_{\rm comp}(\vec{x})$, and the source function $S(\vec{x})$

		$\displaystyle m^{0}_{\rm comp}(\vec{x})=[1,1]^{T},\qquad\mathbb{C}^{0}_{\rm comp}(\vec{x})=[1,0,1]^{T},\qquad S(\vec{x})=E-\bar{E}$			(29)
		$\displaystyle E=\exp(-\sigma(\vec{x}-\vec{x}_{0})^{2}),\,\sigma=1000,\,\vec{x}_{0}=(0.1,0.1)$			(30)

where $\mathbb{I}$ is the identity matrix and $\bar{E}={\rm mean}(E)$.

In Fig. 1 we show three different quantities of the TestD: the module of the variables at final time (first column), the two components of the flux $|\mathbb{C}\nabla p|$ at final time (second column) and the energy as a function of time (third column). In the first row we have the results for the variable $m$ and in the second row the results for the variable $\mathbb{C}$. As expected, the energy decays in time for both variables, and it is very small at final time, which indicates that we are close to the steady state of the systems. The main difference between the two variables are the shape of the network, with a $Y-$shape for the conductivity vector and a $V-$shape for the tensor.

In Fig. 2 we show the results obtained when varying the parameter $D$ in Eqs. (8-9) and in Eqs. (3-6). The tests we consider are: TestG (first column), TestD (second column) and TestE (third column), with $D\in\{0.05,0.01,0.001\}$, for the variable $m$ in the first row and for $\mathbb{C}$ in the second row. The ramifications become more evident when decreasing the diffusivity, and they get thinner and thinner. For the first parameter chosen, $D=0.05$, we are not able to see those ramifications for the vector $m$ because the time-scale associated with the diffusion is too fast to capture the details.

In Fig. 3 we observe the dependence on the relaxation exponent $\gamma$. In the first column we report the results of TestH, in the second, those corresponding to TestD and in the third one, those corresponding to TestF. Again, in the first row we show the results for the variable $m$ and in the second row those for the variable $\mathbb{C}$. If $\gamma=1$ the results do not show the details of the network, and it seems that $\gamma=0.75$ is the parameter that better represents the leaf network.

In Fig. 4 we show the behaviour of the solution $\mathbb{C}$, when $\varepsilon\to 0$. We see the results for $\varepsilon\in\{10^{-2}(\rm left),10^{-3}(\rm center),10^{-4}(\rm right)\}$, and again we notice that for the largest value of $\varepsilon$ we are not able to see any ramification. While we see that for $\varepsilon$ smaller than $10^{-3}$ we are close to the asymptotic behaviour.

In this section we show some quantitative comparison between the two models. For this reason we consider well prepared initial data and we look for compatible parameters.

The goal of this part is to choose a convenient set of parameters in order to compare the two systems, trying to make them as close as possible. Now we distinguish the parameters $(D_{l},c_{l},\alpha_{l},\gamma_{l})$ with $l=1$, for the $\mathbb{C}-$model, and $l=2$, for the $m-$model.

For simplicity, the choice of parameter is performed by comparing the two models in one space dimension. In 1D the systems (3,6) and (8-9) read

	$\displaystyle C_{t}-D_{1}^{2}\,C_{xx}-c_{1}^{2}\,p_{x}^{2}+\alpha_{1}|C|^{\gamma_{1}-2}C$	$\displaystyle=0$		(31)
	$\displaystyle m_{t}-D_{2}^{2}\,m_{xx}-c_{2}^{2}\,p_{x}^{2}\,m+\alpha_{2}|m|^{2(\gamma_{2}-1)}m$	$\displaystyle=0.$		(32)

Now we suppose that $C$ has the following form

$$C=m^{2}+B,$$

(33)

where $B$ is a measure of the discrepancy between the two models, and we set the initial conditions so that $B(t=0)=0$. If we substitute Eq. (33) in Eq. (31), we have

$\displaystyle(m^{2})_{t}-D_{1}^{2}(m^{2})_{xx}-c_{1}^{2}p_{x}^{2}+\alpha_{1}|m^{2}+B|^{\gamma-2}m^{2}=-B_{t}+D_{1}^{2}B_{xx}-\alpha_{1}|m^{2}+B|^{\gamma-2}B.$

(34)

At this point we multiply Eq. (32) by a factor $(2m)$, and we obtain

$$2m\,m_{t}-2D_{2}^{2}m\,m_{xx}-2\,c_{2}^{2}\,p_{x}^{2}\,m^{2}+2\alpha_{2}|m|^{2(\gamma_{2}-1)}m^{2}=0.$$

(35)

After some manipulation, Eqs. (34,35) become:

	$\displaystyle 2m\,m_{t}-D_{1}^{2}(m^{2})_{xx}-c_{1}^{2}p_{x}^{2}+\alpha_{1}\,(m^{2})^{\gamma-1}=\underbrace{-B_{t}+D_{1}^{2}B_{xx}-\alpha_{1}|m^{2}|^{\gamma-2}B}_{:=\mathcal{R}}$		(36)
	$\displaystyle 2m\,m_{t}-D_{2}^{2}(m^{2})_{xx}+2D_{2}^{2}(m_{x})^{2}-2\,c_{2}^{2}\,p_{x}^{2}\,m^{2}+2\alpha_{2}(m^{2})^{\gamma_{2}}=0$		(37)

where $\mathcal{R}$ is the residual. We made use of the following identity in the equation for $m$

$$(m)^{2}_{xx}=2m\,m_{xx}+2\left((m_{x})^{2}\right).$$

(38)

Now we consider the difference between Eq. (36) and Eq. (37), and we obtain

$$(D_{2}^{2}-D_{1}^{2})(m^{2})_{xx}-2D_{2}^{2}(m_{x})^{2}+(2\,c_{2}^{2}m^{2}-c_{1}^{2})p_{x}^{2}+\alpha_{1}\,(m^{2})^{\gamma-1}-2\alpha_{2}(m^{2})^{\gamma_{2}}=\mathcal{R}$$

(39)

Since we want the residual to be small, in absolute value, a convenient choice for the sets of variables is the following

$$D_{1}=D_{2},\quad\alpha_{1}=2\alpha_{2},\quad\gamma_{1}=\gamma_{2}+1,\quad c_{1}=\sqrt{2}\,c_{2}|m|$$

(40)

and for the initial conditions we choose $m^{0}$ and $C^{0}$ such that, initially, we have

$$C^{0}=(m^{0})^{2}={\rm const}.$$

In this way the first derivative in space $m_{x}$ is also equal to 0 after one time step.

In order to show some results in 2D, we need to define the initial conditions for $m^{0}_{\rm comp}$ and $\mathbb{C}^{0}_{\rm comp}$ such that, at initial time, we have again

$$\mathbb{C}^{0}=m^{0}\otimes m^{0},$$

with $m^{0}_{\rm comp}=[\sqrt{2}/2,\sqrt{2}/2]^{T}$ and $\mathbb{C}^{0}_{\rm comp}=[0.5,0.5,0.5]^{T}$, while the values of the parameters are defined in TestM. In 2D, the equivalent expression of (33) is

$$\mathbb{C}=m\otimes m+\mathbb{B}.$$

In this subsection we comment on the solutions of the following tests

In Table 6 we show the time evolution of the norm of the difference between $\mathbb{C}$ and $m\otimes m$, to see how they move away from each other when we increase the time. In this table we see that the two solutions are very different, even after few time steps, and for this reason, they are difficult to be compared. The definition of $||\mathbb{B}||$ is the following

$$||\mathbb{B}||:=\frac{\Big{|}\Big{|}|\mathbb{C}|-|m\otimes m|\Big{|}\Big{|}}{\Big{|}\Big{|}|m\otimes m|\Big{|}\Big{|}}$$

after $n_{t}$ time steps, with ${\rm time}=n_{t}\Delta t$, such that $n_{t}=6(2^{k}),\,k=0,1,2,3,4,5,6$, where $|\mathbb{C}|:=\sqrt{\mathbb{C}_{11}^{2}+2\,\mathbb{C}_{12}^{2}+\mathbb{C}_{22}^{2}}$ and
$|m\otimes m|:=\sqrt{(m_{1}^{2})^{2}+2\,(m_{1}\,m_{2})^{2}+(m_{2}^{2})^{2}}$.

As it appears from the table, the two models move quickly far apart from each other, suggesting an intrinsically different behaviour.

Now we are interested in showing the solution of the $\mathbb{C}-$model, in the case of zero-diffusivity. Since the randomness of the network is common in nature but is also very effective in stabilizing the equations, we want to see if there is some analogy in considering the cases $D=0$ and $D\ll 1$. In Fig. 5 we show the agreement of the two solutions, with $D=0$ (left panel) and $D=10^{-5}$ (right panel). The other parameters are defined in TestN and TestO, while the initial condition is defined in Eq. (29) for Fig. 5.