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

Clarissa Astuto  King Abdullah University of Science and Technology (KAUST),
4700, Thuwal, Saudi Arabia  Daniele Boffi  King Abdullah University of Science
and Technology (KAUST), 4700, Thuwal, Saudi Arabia  Jan Haskovec  King
Abdullah University of Science and Technology (KAUST), 4700, Thuwal, Saudi
Arabia  Peter Markowich  King Abdullah University of Science and Technology
(KAUST), 4700, Thuwal, Saudi Arabia  Giovanni Russo  Department of Mathematics
and Computer Science, University of Catania, Viale Andrea Doria 6, 95125,
Catania, Italy

###### Abstract

We compare the solutions of two systems of partial differential equations
(PDE), seen as two different interpretations of the same model that describes
formation of complex biological networks. Both approaches take into account
the time evolution of the medium flowing through the network, and we compute
the solution of an elliptic-parabolic PDE system for the conductivity vector
$m$, the conductivity tensor $\mathbb{C}$ and the pressure $p$. We use finite
differences schemes in a uniform Cartesian grid in the spatially two-
dimensional setting to solve the two systems, where the parabolic equation is
solved by a semi-implicit scheme in time. Since the conductivity vector and
tensor appear also in the Poisson equation for the pressure $p$, the elliptic
equation depends implicitly on time. For this reason we compute the solution
of three linear systems in the case of the conductivity vector
$m\in\mathbb{R}^{2}$, and four linear systems in the case of the symmetric
conductivity tensor $\mathbb{C}\in\mathbb{R}^{2\times 2}$, at each time step.
To accelerate the simulations, we make use of the Alternating Direction
Implicit (ADI) method.

The role of the parameters is important for obtaining detailed solutions. We
provide numerous tests with various values of the parameters involved, to see
the differences in the solutions of the two systems.

##  1 Introduction

We study two elliptic-parabolic systems of partial differential equations
(PDE) describing the formation of biological network structures. Both systems
are derived as gradient flows of an energy functional, consisting of a
diffusive term, activation term and metabolic cost term. The energy functional
can be seen as a continuum version of its discrete counterpart introduced in
[1]. Here the authors consider a total energy consumption function for a
general class of biological transport networks (seen also in [2]), including a
material cost function of the network. In the papers [1, 3], the authors adapt
the dynamics of local information for biological transport network and their
relation with the optimization principle, which is known to be very common in
nature.

Assuming the validity of Darcy’s law for slow flow in porous media (see, for
instance, [4, 5]), the energy functional is constrained by a Poisson equation
for the fluid pressure. We use two modes of description of the network
conductivity: first, in terms of a conductance vector $m$, and, second, in
terms of a symmetric positive definite conductance tensor $\mathbb{C}$. Taking
the $L^{2}$-gradient flow with respect to $m$ (see, for instance, [6, 7, 8, 9,
10]) and, resp., with respect to $\mathbb{C}$ (see [11]), leads to two
structurally similar elliptic-parabolic PDE systems.

In the parabolic equation, a reaction term appears, representing the metabolic
cost of the network, while the diffusion term describes the randomness in the
material structure. The third term, called activation term, describes the
tendency of the network to align with the principal direction of the material
flow. The elliptic Poisson equation describes local mass conservation and is
equipped with a right-hand side describing the distribution of sources and
sinks of the material. This distribution is supplied as a datum and is
supposed to be time independent.

The aim of this paper is to present several numerical simulations of both the
vector-valued $m$-model and the tensor-valued $\mathbb{C}$-model, with the
goal of comparing their solutions in various parameter settings. The initial
datum is chosen such that $\mathbb{C}=m\otimes m$, i.e., initially, the
principal direction (eigenvector corresponding to the largest eigenvalue) of
$\mathbb{C}$ is aligned with $m$. Note that, in the spatially two-dimensional
setting, the second eigenvalue of $m\otimes m$ is zero, with eigenspace
orthogonal to $m$. We shall observe that the two PDE systems develop
structurally and qualitatively similar solutions, however, they differ in
quantitative details, for instance, the number and location of branches.

##  2 Mathematical model

We introduce the energy functional $\mathcal{E}_{\mathrm{tens}}$ for the
tensor-valued model,

$\mathcal{E}_{\mathrm{tens}}[\mathbb{C}]:=\int_{\Omega}\frac{D^{2}}{2}|\nabla\mathbb{C}|^{2}+c^{2}\nabla
p[\mathbb{C}]\cdot\mathbb{P}[\mathbb{C}]\nabla
p[\mathbb{C}]+M(|\mathbb{C}|)\mathrm{d}x,$ (1)

where $\mathbb{C}=\mathbb{C}(x)\in\mathbb{R}^{2\times 2}$ is the conductivity
tensor, the diffusivity parameter $D\in\mathbb{R}$ measures the effect of
random fluctuations in the network, and the activation parameter $c^{2}>0$
controls the strength of the network formation feedback loop. The total
permeability tensor is of the form
$\mathbb{P}[\mathbb{C}]:=r\mathbb{I}+\mathbb{C}$, where the scalar function
$r=r(x)\geq r_{0}>0$ describes the isotropic background permeability of the
medium. The scalar pressure $p=p[\mathbb{C}]$ of the fluid transported within
the network is the unique solution (up to an additive constant) of the Poisson
equation

$-\nabla\cdot\left(\mathbb{P}[\mathbb{C}]\nabla p\right)=S,$ (2)

subject to homogeneous Neumann boundary condition on $\partial\Omega$. The
source/sink distribution $S=S(x)$ in the mass conservation Eq. (2) is to be
supplemented as an input datum and is assumed to be independent of time. The
metabolic cost function $M:\mathbb{R}^{+}\to\mathbb{R}^{+}$ describes the
dependence of the metabolic expenditure of maintaining the network on its
transportation capacity, see [1]. The expression $|\mathbb{C}|$ denotes the
Frobenius norm $|\mathbb{C}|:=\sqrt{\sum_{i=1}^{2}\sum_{j=1}^{2}{C}_{ij}^{2}}$
and
$|\nabla\mathbb{C}|:=\sqrt{\sum_{i=1}^{2}\sum_{j=1}^{2}\sum_{k=1}^{2}(\partial{C}_{ij}/\partial
x_{k})^{2}}$.

Taking the $L^{2}$-gradient flow of the energy (1) constrained by (2), we
obtain the parabolic-elliptic system

$\displaystyle-\nabla\cdot\left(\left(r\mathbb{I}+\mathbb{C}\right)\nabla
p\right)$ $\displaystyle=S$ (3)
$\displaystyle\frac{\partial\mathbb{C}}{\partial
t}-D^{2}\Delta\mathbb{C}-c^{2}\nabla p\otimes\nabla
p+\frac{M^{\prime}\left(|\mathbb{C}|\right)}{\left|\mathbb{C}\right|}\mathbb{C}$
$\displaystyle=0,$ (4)

see [11] for details of the derivation. In our paper, we shall make the
generic choice for the metabolic cost function (see, e.g., [9, 6, 7, 8]),

$M(s):=\frac{\alpha}{\gamma}s^{\gamma}\qquad\mbox{for }s\geq 0,$ (5)

where $\alpha>0$ is the metabolic constant and $\gamma>0$ the metabolic
exponent. For instance, to model leaf venation in plants, one chooses
$1/2<\gamma<1$, see [1]. Then Eq. (4) becomes

$\frac{\partial\mathbb{C}}{\partial t}-D^{2}\Delta\mathbb{C}-c^{2}\nabla
p\otimes\nabla p+\alpha|\mathbb{C}|^{\gamma-2}\mathbb{C}=0$ (6)

To derive the vector-valued model, we make the ansatz $\mathbb{C}:=m\otimes
m$. Inserting this into (1) with $D=0$ and noting that $|\mathbb{C}|=|m|^{2}$,
we obtain

$\mathcal{E}=\int_{\Omega}c^{2}\nabla p[m]\cdot\mathbb{P}[m]\nabla
p[m]+M\left(|m|^{2}\right)\mathrm{d}x$

with $\mathbb{P}[m]=r\mathbb{I}+m\otimes m$ and, with a slight abuse of
notation, we now denote by $p=p[m]$ the solution of the Poisson equation (2)
with $\mathbb{P}[\mathbb{C}]$ replaced by $\mathbb{P}[m\otimes m]$. Re-
introducing the Dirichlet integral $D^{2}\int_{\Omega}|\nabla
m|^{2}\mathrm{d}x$, we arrive at the energy functional

$\mathcal{E}_{\mathrm{vect}}[m]:=\int_{\Omega}D^{2}|\nabla m|^{2}+c^{2}\nabla
p[m]\cdot\mathbb{P}[m\otimes m]\nabla p[m]+M(|m|^{2})\mathrm{d}x.$ (7)

Note that this functional is different from the one obtained by replacing
$C=m\otimes m$ in Eq. (1), the two functionals mainly differencing in the
contribution of the diffusion term. As we shall see, we are mainly interested
in studying the behavior of the system for very small diffusion coefficients,
so that we expect such difference not being so crucial.

Taking the $L^{2}$-gradient flow with respect to the vector variable $m$ and
using the explicit form (5) for the metabolic function, we obtain the system

$\displaystyle-\nabla\cdot\left(\left(r\mathbb{I}+m\otimes m\right)\nabla
p\right)$ $\displaystyle=S$ (8) $\displaystyle\frac{\partial{m}}{\partial
t}-D^{2}\Delta m-c^{2}(\nabla p\otimes\nabla p)m+{\alpha}|m|^{2(\gamma-1)}m$
$\displaystyle=0.$ (9)

From now on, we shall denote system (3)–(6) as the $\mathbb{C}-$system, while
we call (8)–(9) the $m-$system. Let us now shortly discuss the differences
between the two systems.First, the activation term $c^{2}\nabla p\otimes\nabla
p$ in (4) is quadratic and depends on $\mathbb{C}$ only through the solution
$p=p[\mathbb{C}]$ of the Poisson equation (3). In contrast, the activation
term $c^{2}(\nabla p\otimes\nabla p)m$ in (9) is cubic. Moreover, the
metabolic term $\alpha|\mathbb{C}|^{\gamma-2}\mathbb{C}$ in (4) becomes
singular at $\mathbb{C}=0$ if (and only if) $\gamma<1$. This obviously causes
difficulties for both the analytical treatment of the system (see [11] for
details) and its numerical resolution. We shall discuss in Section 3.2.1 how
the numerical difficulties can be overcome. In contrast, the metabolic term
${\alpha}|m|^{2(\gamma-1)}m$ in (9) only becomes singular at $m=0$ if
$\gamma<1/2$. Taking into account that for typical applications in biology
only the parameter range $\gamma\geq 1/2$ is relevant (see [1]), the metabolic
term in (9)) does not require any special treatment for the $m-$system.

We define the equations on a domain $\Omega\subset\mathbb{R}^{2}$, and we
define the boundary conditions on $\partial\Omega$ for the pressure and the
conductivity. We choose homogeneous Dirichlet boundary conditions for $m$ and
$\mathbb{C}$, and homogeneous Neumann conditions for $p$.

$m(t,\vec{x})=0,\quad\mathbb{C}(t,\vec{x})=0,\quad\mathbb{P}[\mathbb{C}]\nabla
p(t,\vec{x})\cdot\nu=0,\quad\vec{x}\in\partial\Omega,\,t\geq 0$ (10)

where $\nu$ is the outgoing normal vector to $\partial\Omega$. In all our
numerical test, the numerical support of $\mathbb{C}$ or $m$, after a long
time, is well within $\Omega$, therefore the solution should not be
particularly sensitive to the boundary conditions on the two variables. In any
case, different boundary conditions for $m$ and $\mathbb{C}$ are currently
under investigation.

To close the system, we prescribe an initial condition for the conductivity
vector and tensor

$m(t=0,\vec{x})=m^{0}(\vec{x}),\quad\mathbb{C}(t=0,\vec{x})=\mathbb{C}^{0}(\vec{x}),\quad\text{
in }\Omega.$ (11)

A direct consequence of the boundary condition defined for the pressure, that
is also a necessary condition for the solvability of the Poisson equation, is
that the source function has to be vanishing mean, i.e.,
$\displaystyle\int_{\Omega}S(\vec{x})d\Omega=0$.

##  3 Numerical schemes

In this section we define the numerical schemes used to discretize in space
and time the Eqs. (8 \- 9) and (3-4). We adopt semi-implicit second order
schemes.

###  3.1 Space discretization

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.

###  3.2 Time discretization: symmetric-ADI method

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.

####  3.2.1 Time discretization for the conductivity vector

Here we focus on the time discretization for the Eqs. (8-9).

Given $m^{n}\approx m(t^{n})$, we compute $p^{n}$ by solving the Poisson
equation

$-\mathcal{L}\left(m^{n}\otimes m^{n}\right)\,p^{n}=S,$ (26)

where $\mathcal{L}\left(m^{n}\otimes m^{n}\right)\in\mathbb{R}^{N^{2}\times
N^{2}}$ is the discrete elliptic operator, in both directions ($x$ and $y$),
with variable coefficients and corresponding to zero Neumann conditions.

The symmetric-ADI method to solve the Eq. (17) works as follows. We start with
the $y$-direction implicit and $x$-direction explicit, and then we consider
the opposite order in the second step of the ADI scheme

$\displaystyle(y-{\rm impl},x-{\rm expl})\qquad\tilde{m}$ $\displaystyle=$
$\displaystyle m^{n}+\frac{\Delta t}{2}\mathcal{L}_{y}\tilde{m}+\frac{\Delta
t}{2}\mathcal{L}_{x}m^{n}+\Delta t\,c^{2}\mathcal{P}_{xy}^{n}(m^{n})$

and we solve for $\tilde{m}$,

$\displaystyle(x-{\rm impl},y-{\rm expl})\qquad m^{n+1}_{y}$ $\displaystyle=$
$\displaystyle\tilde{m}+\frac{\Delta
t}{2}\mathcal{L}_{x}{m^{n+1}_{y}}+\frac{\Delta
t}{2}\mathcal{L}_{y}\tilde{m}-\Delta
t\,\alpha\mathcal{Q}^{m}\left({m^{n}}\right)m^{n+1}_{y}$ $\displaystyle+\Delta
t\,c^{2}\mathcal{P}_{y}^{n}\,m_{y}^{n+1}$

and we solve for $m^{n+1}_{y}$.

The second time we apply the ADI method, we first consider the $x$-direction
implicit and $y$ explicit, and then we exchange the order. Thus we have

$\displaystyle(x-{\rm impl},y-{\rm expl})\qquad\hat{m}$ $\displaystyle=$
$\displaystyle m^{n}+\frac{\Delta t}{2}\mathcal{L}_{x}\hat{m}+\frac{\Delta
t}{2}\mathcal{L}_{y}{m}^{n}+\Delta t\,c^{2}\mathcal{P}_{xy}^{n}(m^{n})$

here we solve for $\hat{m}$,

$\displaystyle(y-{\rm impl},x-{\rm expl})\qquad m^{n+1}_{x}$ $\displaystyle=$
$\displaystyle\hat{m}+\frac{\Delta
t}{2}\mathcal{L}_{y}{m^{n+1}_{x}}+\frac{\Delta
t}{2}\mathcal{L}_{x}\hat{m}-\Delta
t\,\alpha\mathcal{Q}^{m}\left({m^{n}}\right)m_{x}^{n+1}$ $\displaystyle+\Delta
t\,c^{2}\mathcal{P}_{x}^{n}\,m_{x}^{n+1}$

and now we solve for $m^{n+1}_{x}$, where
$\mathcal{Q}^{m}\left(m^{n}\right)=|m^{n}|^{2(\gamma-1)}$ and
$\displaystyle\mathcal{L}_{\beta}$, with $\beta=x,y$, are the discrete
operators for the second derivatives in $x$ and $y$ directions, respectively,
with $\displaystyle\mathcal{L}_{\beta}\in\mathbb{R}^{N\times N}$. For the
pressure term we have $\mathcal{P}_{x}^{n}=(\mathcal{D}_{x}p^{n})^{2}$ for the
first component of the vector $m^{(1)}$, and
$\mathcal{P}_{y}^{n}=(\mathcal{D}_{y}p^{n})^{2}$ for the second component
$m^{(2)}$. While the force term is
$\mathcal{P}^{n}_{xy}(m^{n})=[\mathcal{D}_{x}p^{n}\mathcal{D}_{y}p^{n}\,m^{(2)n},\mathcal{D}_{x}p^{n}\mathcal{D}_{y}p^{n}\,m^{(1)n}]^{T}$.

At the end, we calculate the average of the two solutions $m^{n+1}_{y}$ and
$m^{n+1}_{x}$ to obtain the conductivity vector at time $t^{n+1}$

$\displaystyle m^{n+1}$ $\displaystyle=$
$\displaystyle\frac{1}{2}{m^{n+1}_{x}}+\frac{1}{2}{m^{n+1}_{y}}.$ (27)

####  3.2.2 Time discretization for the conductivity tensor

As we said in the description of the two systems, for the conductivity tensor
$\mathbb{C}$, the reaction term is very stiff because of the exponent $\gamma$
that belongs to the interval $(0.5,1)$. After a finite time, we are basically
dividing by zero at each time step. For this reason we introduce a small
regularizing parameter $\varepsilon$ in the equation, as follows

$\mathcal{Q}^{c}(\mathbb{C})=|\mathbb{C}+\varepsilon|^{\gamma-2}\mathbb{C}$
(28)

and we study the behaviour of the system as $\varepsilon$ becomes smaller and
smaller.

Given the pressure at time $t^{n}$ from Eq. (26), we apply the symmetric-ADI
method to solve the Eq. (18). As explained before, we first choose the
$y$-direction implicit and then we exchange the two directions. The scheme
reads

$\displaystyle(y-{\rm impl},x-{\rm expl})\qquad\tilde{\mathbb{C}}_{1}$
$\displaystyle=$ $\displaystyle\mathbb{C}^{n}+\frac{\Delta
t}{2}\mathcal{L}_{y}\tilde{\mathbb{C}}_{1}+\frac{\Delta
t}{2}\mathcal{L}_{x}\mathbb{C}^{n}+\Delta t\,\mathcal{P}^{n}$

and we solve for $\tilde{\mathbb{C}}_{1}$,

$\displaystyle(x-{\rm impl},y-{\rm expl})\qquad\mathbb{C}^{n+1}_{y}$
$\displaystyle=$ $\displaystyle\tilde{\mathbb{C}}_{1}+\frac{\Delta
t}{2}\mathcal{L}_{x}{\mathbb{C}^{n+1}_{y}}+\frac{\Delta
t}{2}\mathcal{L}_{y}\tilde{\mathbb{C}}_{1}-\Delta
t\,\alpha\mathcal{Q}^{c}\left({\mathbb{C}^{n}}\right)\mathbb{C}^{n+1}_{y}$

and we solve for $\mathbb{C}^{n+1}_{y}$.

Now, as before, we apply the ADI method for the second time for $\mathbb{C}$
starting with the $x$-direction implicit, and we have

$\displaystyle(x-{\rm impl},y-{\rm expl})\qquad\tilde{\mathbb{C}}_{2}$
$\displaystyle=$ $\displaystyle\mathbb{C}^{n}+\frac{\Delta
t}{2}\mathcal{L}_{x}\tilde{\mathbb{C}}_{2}+\frac{\Delta
t}{2}\mathcal{L}_{y}{\mathbb{C}}^{n}+\Delta t\,\mathcal{P}^{n}$

here we solve for $\tilde{\mathbb{C}}_{2}$,

$\displaystyle(y-{\rm impl},x-{\rm expl})\qquad\mathbb{C}^{n+1}_{x}$
$\displaystyle=$ $\displaystyle\tilde{\mathbb{C}}_{2}+\frac{\Delta
t}{2}\mathcal{L}_{y}{\mathbb{C}^{n+1}_{x}}+\frac{\Delta
t}{2}\mathcal{L}_{x}\tilde{\mathbb{C}}_{2}-\Delta
t\,\alpha\mathcal{Q}^{c}\left({\mathbb{C}^{n}}\right)\mathbb{C}_{x}^{n+1}$

and here we solve for $\mathbb{C}^{n+1}_{x}$.

Finally we calculate $\mathbb{C}^{n+1}$ from the two solutions
$\mathbb{C}^{n+1}_{y}$ and $\mathbb{C}^{n+1}_{x}$

$\displaystyle\mathbb{C}^{n+1}$ $\displaystyle=$
$\displaystyle\frac{1}{2}{\mathbb{C}^{n+1}_{x}}+\frac{1}{2}{\mathbb{C}^{n+1}_{y}}.$

The quantities above have the following expressions:
$\mathcal{Q}^{c}\left(\mathbb{C}^{n}\right)=|\mathbb{C}^{n}+\varepsilon|^{\gamma-2}$,
and for the pressure,
$\mathcal{P}^{n}=\mathcal{D}_{\beta}p^{n}\,\mathcal{D}_{\eta}p^{n}$, with the
suitable choice of $\beta,\eta\in\\{x,y\\}$ for the four components of the
conductivity tensor.

##  4 Numerical results

In this section we perform several simulations with the aim of studying the
effect of the various parameters. In particular, we check the agreement of the
two models for the $m-$system in Eqs. (8-9) and for the $\mathbb{C}-$system in
Eqs. (3-6).

###  4.1 Accuracy tests and qualitative agreements

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.

|  | $\alpha$ | $c$ | $D$ | $\varepsilon$ | $\gamma$ | $r$ | $t_{\rm fin}$
---|---|---|---|---|---|---|---|---
Accuracy $m$ |  TestA: | 0.5 | 1 | 0.01 | - | 0.75 | 0.1 | 1
Accuracy $\mathbb{C}$ |  TestB: | 1 | 1 | 0.01 | 0.1 | 1.75 | 0.1 | 1
Accuracy $m$ |  TestC: | 0.5 | 5 | 0.01 | - | 0.75 | 0.01 | 1
$D=0.05$ |  TestG: | 0.75 | 5 | 0.05 | $10^{-3}$ | 0.75 | 0.005 | 15
$D=0.01$ |  TestD: | 0.75 | 5 | 0.01 | $10^{-3}$ | 0.75 | 0.005 | 15
$D=0.001$ |  TestE: | 0.75 | 5 | 0.001 | $10^{-3}$ | 0.75 | 0.005 | 15
$\gamma=1$ |  TestH: | 0.75 | 5 | 0.01 | $10^{-3}$ | 1 | 0.005 | 15
$\gamma=0.75$ |  TestD: | 0.75 | 5 | 0.01 | $10^{-3}$ | 0.75 | 0.005 | 15
$\gamma=0.5$ |  TestF: | 0.75 | 5 | 0.01 | $10^{-3}$ | 0.5 | 0.005 | 15
$\varepsilon=10^{-2}$ |  TestI: | 0.75 | 5 | 0.01 | $10^{-2}$ | 0.75 | 0.005 | 15
$\varepsilon=10^{-3}$ |  TestD: | 0.75 | 5 | 0.01 | $10^{-3}$ | 0.75 | 0.005 | 15.
$\varepsilon=10^{-4}$ |  TestL: | 0.75 | 5 | 0.01 | $10^{-4}$ | 0.75 | 0.005 | 15
Table 1: In this table we define all the tests that we show in Section4.1. The
first three rows show the parameters for the accuracy tests for $m$ and
$\mathbb{C}$, and the results are summarized in Table 2-3. The second three
rows define the parameters that we use in Fig. 2, where we compare the results
changing the diffusivity. The third three rows are the tests showed in Fig. 3,
varying $gamma$ and the results of the last three rows are in Fig. 4, where we
change the stabilization parameter $\varepsilon$. For the accuracy tests, the
number of points of the discretization is specified in Table 2-3, while, for
all the other tests, the number of points is fixed and it is $N=600$.

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)$.

N | error | order
---|---|---
20 | - | -
40 | 0.036030 | -
80 | 0.0492860 | -0.4520
160 | 0.01454106 | 1.7610
320 | 0.00690830 | 1.0737
640 | 0.001529779 | 2.1750

N | error | order
---|---|---
20 | - | -
40 | 0.036012 | -
80 | 0.0493010 | -0.4531
160 | 0.01456192 | 1.7594
320 | 0.00691103 | 1.0752
640 | 0.001528055 | 2.1772

Table 2: Accuracy test of the $m$-system (8-9): we show the $L^{2}-$norm of the relative error for $|m|$, with the parameters defined in TestA (left) and TestC (right). $N$ | error2 | order
---|---|---
25 | - | -
50 | $9.066\times 10^{-2}$ | -
100 | $4.625\times 10^{-2}$ | 0.97
200 | $1.571\times 10^{-2}$ | 1.56
400 | $4.149\times 10^{-3}$ | 1.92
800 | $7.347\times 10^{-4}$ | 2.50
Table 3: Accuracy test of the $\mathbb{C}$-system (3-4): we show the
$L^{2}-$norm of the relative error for $|\mathbb{C}|$, with the parameters
defined in TestB.

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.

Figure 1: In this figure we show three different quantities of the same
computations, with the parameters defined in TestD: the module of the
variables at final time (left panels), the flux also at final time (central
panels) and the energy as a function of time (right panels). The first row is
about the variable $m$ and the second one is for the variable $\mathbb{C}$.

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.

Figure 2: In this figure we show the difference of the results on varying the
diffusivity $D$. In the first column we have the results of the TestG (D =
0.05), in the second column we choose the parameters of the TestD (D = 0.01)
and in the third one the  TestE (D = 0.001). The first row shows the results
of the variable $m$ and the second row those of the variable $\mathbb{C}$.

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.

Figure 3: In this figure we show the difference of the results on varying the
relaxation exponent $\gamma$. In the first column we have the results of the
TestH ($\gamma$ = 1), in the second column we choose the parameters of the
TestD ($\gamma=0.75$) and in the third one the  TestF ($\gamma=0.5$). The
first row shows the results of the variable $m$ and the second row of the
variable $\mathbb{C}$.

Figure 4: In this figure we show the results for the variable $\mathbb{C}$
varying the parameter $\varepsilon$. On the left we have the results of the
TestI ($\varepsilon=10^{-2}$) in the central panel TestD
$(\varepsilon=10^{-3})$ and on the right  TestL $(\varepsilon=10^{-4})$.

###  4.2 Quantitative agreement

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

|  | $\alpha_{1},\alpha_{2}$ | $c_{1},c_{2}$ | $D_{1}=D_{2}$ | $\varepsilon$ | $\gamma_{1},\gamma_{2}$ | $r$ | $N$
---|---|---|---|---|---|---|---|---
set of parameters |  TestM: | 1, 0.5 |  $\sqrt{2}$, 1 | 0.1 | $10^{-1}$ | 1.75, 0.75 | 0.1 | 600
Table 4: In this table we define the two set of parameters in Eq. (40).  |  | $\alpha$ | $c$ | $D$ | $\varepsilon$ | $\gamma$ | $r$ | $t_{\rm fin}$ | $N$
---|---|---|---|---|---|---|---|---|---
$D=0$ |  TestN: | 0.75 | 5 | 0 | $10^{-3}$ | 0.75 | 0.005 | 15 | 600
$D=10^{-5}$ |  TestO: | 0.75 | 5 | $10^{-5}$ | $10^{-3}$ | 0.75 | 0.005 | 15 | 600
Table 5: In this table we define the two tests showed in Figs.5-6, with $D=0$
(TestN) and $D=10^{-5}$ (TestO).

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}}$.

time | 0.01 | 0.02 | 0.04 | 0.08 | 0.16 | 0.32 | 0.64
---|---|---|---|---|---|---|---
$||\mathbb{B}||$ | 0.0348 | 0.0538 | 0.0697 | 0.0982 | 0.1509 | 0.2611 | 0.5320
Table 6: Here we show the values of $\mathbb{B}$, after $n_{t}$ time steps,
where $n_{t}=2^{k},k=0,1,2,3,4,5,6$.

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.

In Fig. 6 we illustrate the long time solution for the $\mathbb{C}$-model
obtained with the following space dependent initial condition:

$\mathbb{C}^{0}=f(x,y)\mathbb{I},\quad f(x,y)=(2-|X+Y|)\exp(-10(|X-Y|))$ (41)

Figure 5: Comparison between TestN and TestO, with initial condition defined
in Eq. (29). The main difference is that, on the left we have $D=0$, while, on
the right, $D=10^{-5}$.

Figure 6: Comparison between TestN and TestO. Here the initial condition is
defined in Eq. (41), with $D=0$ on the left and $D=10^{-5}$ on the right.

We also note that, if we set the diffusion coefficient equal to zero, the
model for $m$ reduces to a reaction equation for the conductivity. This means
that if $\gamma>1/2$ the support of the unknown $m$ remains unchanged. In
particular, it cannot extend, while in some regions the numerical support
(i.e. the region in which $|m|$ becomes lower than a given small threshold)
may shrink.

##### Alternative boundary conditions

Furthermore, in the case of zero-diffusivity for the $m-$model, we observe an
anomalous behavior of the solution near the boundaries. In order to overcome
such a problem we propose an ad hoc boundary condition as illustrated below.

Let us consider the equations for $m$ with $D=0$ and in the limit of steady
state. We have

$\displaystyle c^{2}\nabla p\otimes\nabla p\,m-\alpha|m|^{2(\gamma-1)}m$
$\displaystyle=0\quad{\rm in}\,\partial\Omega$

that we can write as follows

$\displaystyle(m\cdot\nabla p)\nabla
p=\frac{\alpha}{c^{2}}|m|^{2(\gamma-1)}m\quad{\rm in}\,\partial\Omega$ (42)

Thus, we can deduce that $m\propto\nabla p$. This means that there exists a
constant $\beta$, such that,

$m=\beta\nabla p.$ (43)

If we substitute the Eq. (43) in Eq. (42), we obtain

$\beta|\nabla p|^{2}\nabla p=\frac{\alpha}{c^{2}}\beta^{2(\gamma-1)}|\nabla
p|^{2(\gamma-1)}\beta\nabla p$ (44)

and if we solve it for $\beta$, we have the boundary condition for $m$

$\displaystyle\left.m\right|_{\partial\Omega}=\underbrace{\left(\frac{c^{2}}{\alpha}|\nabla
p|^{4-2\gamma}\right)^{\frac{1}{2(\gamma-1)}}}_{\beta}\nabla p.$ (45)

Analogously, we can find also a boundary condition for the vector
$\mathbb{C}$. Again, starting with zero-diffusivity and at steady state, we
have

$\nabla p\otimes\nabla
p=\frac{\alpha}{c^{2}}|\mathbb{C}|^{\gamma-2}\mathbb{C}\quad{\rm
in}\,\partial\Omega$ (46)

that means that the conductivity is proportional to the tensor product of the
pressure gradient, i.e. $\mathbb{C}\propto\nabla p\otimes\nabla p$. Again, we
look for a constant $\tilde{\beta}$, such that,

$\mathbb{C}=\tilde{\beta}\nabla p\otimes\nabla p.$ (47)

Now we substitute the Eq. (47) in Eq. (46), and we solve for $\tilde{\beta}$.
In this way, as before, we find the expression for the conductivity tensor at
the boundary

$\left.\mathbb{C}\right|_{\partial\Omega}=\underbrace{\left(\frac{c^{2}}{\alpha}|\nabla
p|^{-2(\gamma-2)}\right)^{\frac{1}{\gamma-1}}}_{\tilde{\beta}}\nabla
p\otimes\nabla p.$ (48)

Conditions (45,48) might be a reasonable choice in the case of zero
diffusivity. This treatment has the drawback of introducing additional non-
linearity to the system. Alternative boundary conditions are currently under
investigation.

Another aspect we want to focus on, is the steady state for the $m-$model. In
two different cases: $\gamma<1$ and $\gamma>1$ (as the authors show in [7]).
For this reason we define different initial conditions for the vector $m$,
such that

$\displaystyle m^{0,1}_{1}=1,\quad m^{0,1}_{2}=\sqrt{2};\quad
m^{0,2}_{1}=5,\quad m^{0,2}_{2}=5;$ $\displaystyle
m^{0,3}_{1}=(2-|X+Y|)\exp{(-10|X-Y|)},\quad m^{0,3}_{2}=m^{0,3}_{1};$

and in Fig. 7 we plot the following quantity

${\rm diff}=\frac{||m^{\beta}-m^{\eta}||}{||m^{\eta}||}$ (49)

where $||\cdot||$ is the Frobenius norm, and $\beta,\eta\in\\{1,2,3\\}$. In
this way we see the difference between the solutions (with initial condition
$m^{0,1}$ and $m^{0,2}$ in the left panel and with $m^{0,1}$ and $m^{0,3}$ in
the right panel), as function of time. In this way, we support with numerical
evidences that, for the $m-$model and for $\gamma>1$, the steady state is
unique and it does not depend on the initial conditions, as expected (see e.g.
[7]).

In Fig. 8, that is the case for $\gamma<1$, we see that we reach two different
steady states, when choosing two different initial conditions, suggesting that
the steady state solution is not unique when $\gamma<1$.

Figure 7: In this figure we show the difference between two different
solutions choosing, as initial condition, $m^{0,1},m^{0,2}$ (on the left) and
$m^{0,1},m^{0,3}$ (on the right). We plot the expression defined in Eq. (49),
as a function of time, with $\gamma=1.75>1$, and the others parameters are
defined in TestD.

Figure 8: In this figure we show the steady states when $\gamma=0.75<1$. On
the left the initial condition is $m^{0}=m^{0,1}=1$ while, on the right, the
initial condition is a function of space, $m^{0}=m^{0,3}$, and the parameters
are defined in TestD.

##  5 Conclusions

In this paper we use an elliptic-parabolic model to study the formation of
biological network, and, in particular, of leaf venation networks. Throughout
the paper we compare the solutions of two different systems, that derive from
the Cai-Hu model, one for the conductivity vector and one for the conductivity
tensor, and we explore the dependence of the solution on the parameters of the
two models.

As we said before, all the components of the two systems are very stiff. In
particularly, the $\mathbb{C}$-system is more challenging because of the
negative exponent in the reaction term. For this reason, we add a
regularization parameter $\varepsilon$, and we compute numerical solutions for
smaller and smaller values of $\varepsilon$. This parameter prevents the
instability coming from the division by zero in the reaction term.

We make use of finite differences scheme to compute the two solutions, with
central differences for the space discretization and a symmetric-ADI method in
time. The convergence rate is calculated numerically and denotes the second
order accuracy.

At the end we added some quantitative comparison between the two systems,
choosing more suitable sets of parameters. We see that the two solutions
differ significantly, even after few time steps. This aspect makes it
problematic any kind of direct comparison between the two systems.

Then, we showed some results in the case of zero diffusivity for the
conductivity tensor. The last tests we consider are in agreement with the
results achieved in [7], where the authors prove that there is a unique steady
state for the $m-$system when the metabolic exponent $\gamma$ is greater than
1, and provide evidence that this is not true when $\gamma<1$.

## Acknowledgments

G.R. thanks ITN-ETN Horizon 2020 Project ModCompShock, Modeling and
Computation on Shocks and Interfaces, Project Reference 642768, and the
Italian Ministry of Instruction, University and Research (MIUR) to support
this research with funds coming from PRIN Project 2017 (No.2017KKJP4X entitled
”Innovative numerical methods for evolutionary partial differential equations
and applications”.

## References

  * [1] Dan Hu and David Cai.  Adaptation and optimization of biological transport networks.  Physical review letters, 111(13):138701, 2013.
  * [2] Eleni Katifori, Gergely J Szöllősi, and Marcelo O Magnasco.  Damage and fluctuations induce loops in optimal transport networks.  Physical review letters, 104(4):048704, 2010.
  * [3] Dan Hu.  Optimization, adaptation, and initialization of biological transport networks.  Notes from lecture, 1:3–1, 2013.
  * [4] Henry Darcy.  Les fontaines publiques de la ville de Dijon: exposition et application… Victor Dalmont, 1856.
  * [5] Shlomo P Neuman.  Theoretical derivation of darcy’s law.  Acta mechanica, 25(3):153–170, 1977.
  * [6] Di Fang, Shi Jin, Peter Markowich, and Benoit Perthame.  Implicit and Semi-implicit Numerical Schemes for the Gradient Flow of the Formation of Biological Transport Networks.  The SMAI journal of computational mathematics, 5:229–249, 2019\.
  * [7] Jan Haskovec, Peter Markowich, and Benoit Perthame.  Mathematical analysis of a pde system for biological network formation.  Communications in Partial Differential Equations, 40(5):918–956, 2015.
  * [8] Jan Haskovec, Peter Markowich, Benoît Perthame, and Matthias Schlottbom.  Notes on a pde system for biological network formation.  Nonlinear Analysis, 138:127–155, 2016.  Nonlinear Partial Differential Equations, in honor of Juan Luis Vázquez for his 70th birthday.
  * [9] Giacomo Albi, Marco Artina, Massimo Foransier, and Peter A. Markowich.  Biological transportation networks: Modeling and simulation.  Analysis and Applications, 14(01):185–206, 2016.
  * [10] Jan Haskovec, Peter Markowich, and Simone Portaro.  Emergence of biological transportation networks as a self-regulated process, 2022.
  * [11] Jan Haskovec, Peter Markowich, and Giulia Pilli.  Tensor pde model of biological network formation.  arXiv preprint arXiv:2111.03889, 2021.
  * [12] P. Wesseling.  Principles of computational fluid dynamics.  Springer, 2001.
  * [13] Lewis Fry Richardson.  Ix. the approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam.  Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character, 210(459-470):307–357, 1911.


# Self-sustained price bubbles driven by Bitcoin innovations and adaptive
behavior

Misha Perepelitsa and Ilya Timofeyev  Department of Mathematics
University of Houston
4800 Calhoun Rd.
Houston, TX. email<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

We show that infinite divisibility of a trading commodity leads to a self-
sustained price bubble when traders use adaptive investment strategies. The
adaptive strategy can be viewed as a psychological response of a trader to the
situation when the trader’s estimation of future prices does not match the
actual, realized price.

We use a multi-agent model to illustrate the price bubble formation and to
quantify its main statistical properties such as the return, the volatility,
and the systematic risk of the price bubble to crash. We discuss the
plausibility for bubbles to drive prices of digital currencies.

##  1\. Introduction

With an introduction of Bitcoin in 2010, digital currencies became a global
social phenomenon and they are gaining popularity with investors. One of the
distinctive features of digital currencies is their ultimate divisibility, as
virtually any fraction amount of it can be traded. Currently, the minimum
trading amount for Bitcoin is set at $10^{-6}$ units, but the design of the
currency allows it to be lowered by any order of magnitude. This, technically,
solves the problem of liquidity, as no matter how high the price is, one can
always buy some of it, say, with 1 dollar. And it seems that the divisibility
indeed comes handy, as the price for Bitcoin skyrocketed in the last decade.
It’s growth is attributed to many factors such as the novelty effect, the
demand from the black market, global access to the currency through the
internet, the usage of bitcoin as storage in markets with hyper-inflated
national currencies, and, not least, the fact that bitcoin became a sort of a
social cult of a decentralized currency, see Dodd (2017). Apart from the black
market demand, which alone can not contribute to this phenomenal growth, all
other factors point at a price bubble, which seems to be a consensus among the
economists, see an article by Wolff-Mann (2018).

In this paper we show that this unique property of a traded commodity to be
infinitely divisible makes possible the emergence of self-sustained price
bubbles. By this we mean that a group of investors can “pump” the price up at
a positive rate, without ever running out of cash. More than that, we show it
may happen without any collusion between the investors, if each of them
follows a psychologically appealing, adaptive investment strategy.

It is hard to estimate if a speculative bubble of this kind is behind the
growth of the Bitcoin. Several mechanisms a likely to be at work here. For
example, the bubble might have started as an adaptive speculative bubble that
later morphed into a spontaneously organized, global Ponzi scheme, or the
other way around. Other mechanisms of price bubbles in digital currencies are
possible, see for example Griffin and Shams (2020), or, for generic models of
bubbles, Sornette (2003). It is certainly not the goal of this work to fit the
parameters of the model to the time series of Bitcoin prices and to estimate
its statistical significance. Rather than that, our intention was to show the
potentiality of this sort of bubble to develop in the market of digital
currencies.

Let us mention that the adaptive behavior of investors is not limited to
trading of digital currencies but rather is an ubiquitous phenomenon. Thus the
model for the price bubble dynamics, considered in this paper, has a wider
applicability. In particular, it shows that for “classical commodities” there
must be a continuous in-flow of investments for a price bubble to persist, and
provides a analytical tool to obtain quantitative information on the operation
of a price bubble.

The paper organized as follows. In section 2 we describe the model and its
main properties and provide numerical estimates of the return and the
volatility of the price.

In section 3, we give numerical estimates of the systematic risk for the price
bubble to crash, measured by two different metrics, one based on the low cash-
to-stock portfolios of traders, and the other, on low current-cash-to-initial-
cash ratio. We obtain two different estimates of the expected stock price
adjusted by the risk of a crash.

We argue that the second metric is a more relevant measure of risk for the
adaptive speculative bubbles. This gives another evidence of a bubble being
self-sustained. On the other hand, we argue that overestimation of the risk by
the first metric and the high volatility of the returns may contribute to
longevity of speculative bubbles generated by this mechanism in efficient
markets.

Additionally, we describe a relation between the volatility of the price and
the systematic risk, see formula (7) on page 7. It allows for the direct
estimation of the systematic risk from an empirical values of the volatility
of the price time series. Several interesting statistical laws of a price
bubble in a fixed “psycho-state” are derived from this formula, such as the
distribution of the times of a crash, and the maximum expected price that can
be attained in a life-time of a bubble.

In the section 4 we give a psychological interpretation of two main parameters
of our model and consider a possibility that these parameters have cyclic
dynamics in the domain of the psychological variables, thus defining an
intrinsic life span of a speculative bubble. This view of a price dynamics as
an output of a cyclic model suggests that the price dynamics generated by this
model can be in agreement with the no-arbitrage principle.

Models of adaptive economic behavior were developed by many authors in many
different contexts. In strategic decision making they were advanced as models
of convergence by Cross (1983), models of learning in games by Fudenberg and
Levine (1998), and reinforcement learning models by Roth and Erev (1995), to
mention just a few contributions. Our work is close to the methodology of
multi-agent dynamic simulations of Stigler (1964), Kim and Markowitz (1989),
Arthur et al. (1997), Levy et al. (1994, 1995, 2000), Lux (1995, 1998), Lux
and Marchesi (1999, 2000). Price bubbles in efficient markets is a well-known
phenomenon. They were discussed, for example, by Mandelbrot (1966), Blanchard
and Watson (1982) and Sornette (2003).

##  2\. Asynchronous stochastic price pump

Adaptive behavior is an universal paradigm. Not surprisingly it is widely used
as an investment strategy. Consider for a example an investor who has
investment in a stock and has certain cash reserve and would like to keep the
stock-to-cash ratio at certain level, so that the investor is regularly buying
or selling stocks to be at that target level. We let our investor, perhaps a
bitcoin enthusiast, to be somewhat uncertain about the investment goals, and
in particular, to have infinite investment horizon. This might be rather
natural assumption in the case of the digital currencies, since they have been
introduced just recently. The investor might consider the market as a complex,
exogenous environment with constantly changing characteristics. As an adaptive
strategy the investor may increase his investment in stock if the market
outperformed the investor’s expectations as measured by the stock-to-cash
ratio and decrease the investment, otherwise. For example, a trader may
increase the target stock-to-cash ratio by a factor of $\alpha>1,$ in the
former case, and decrease it by the factor of $\beta<1,$ in the latter.
Parameters $\alpha$ and $\beta$ express a degree of the psychological response
of the investor to the market performance and can be thought of as measures of
investor’ optimism and pessimism, or, greed and fear of losing.

What will happen if all other investors in the market follow similar
strategies, trying to adapt to the changing market environment that they
themselves create? One might expect some sort of anomaly to arise. The
scenario was considered by Perepelitsa and Timofeyev (2019) based on the
analysis of a multi-agent model, called asynchronous stochastic price pump
(ASPP) that we describe below.

###  2.1. 2-agent model

To introduce the model we start with a market of only two agents trading
shares of a single asset. The agents, do not know the asset fundamental value
(and maybe there isn’t), but expect that the price can grow at moderate rates
for long periods of time, and they always prefer increasing values to
decreasing.

The agents meet at time intervals that we label with $n=0,1,2..$ At period
$n,$ the first agent’s portfolio has $(s^{n}_{1},b^{n}_{1}),$ dollar value of
stock and bond investments, and the second: $(s^{n}_{2},b^{n}_{2}).$ The stock
price per share at period $n$ is $P^{n}.$ The next period investment portfolio
$(s^{n+1}_{i},b^{n+1}_{i})$ and price $P^{n+1}$ determined trough the
following steps.

Evaluating the last change in their portfolios, agents set target stock-to-
bond ratios for the next period, that we denote by $k^{n+1}_{i}.$ The new
price $P^{n+1}$ is set in such a way that the dollar amount of funds that
agent 1 wishes to move from stocks to bonds, to be at the target ratio
$k^{n+1}_{1},$ equals the dollar amount agent 2 wants to move from bonds to
stock to be at ratio $k^{n+1}_{2}.$ This balance is expressed as

(1)
$\frac{P^{n+1}}{P^{n}}s^{n}_{1}{}-{}k^{n+1}_{1}b^{n}_{1}{}={}k^{n+1}_{2}b^{n}_{2}{}-{}\frac{P^{n+1}}{P^{n}}s^{n}_{2}.$

Here we assume that any fractional amount of a share can be exchanged, and for
simplicity of the exposition, consider the case with zero interest in the bond
account. With the newly set price the agents re-balance portfolios and
evaluate the market performance by comparing ratios
$P^{n+1}s^{n}_{i}/(P^{n}b^{n}_{i})$ to target values $k^{n+1}_{i}.$ If
$P^{n+1}s^{n}_{i}/(P^{n}b^{n}_{i})>k^{n+1}_{i},$ from the agent $i$ point of
view, the market performed better than expected, and if
$P^{n+1}s^{n}_{i}/(P^{n}b^{n}_{i})<k^{n+1}_{i}$ it performed worse. Since
while one agent is selling shares the other one is buying, they have different
opinions on the market performance. The agents update the target portfolio
ratios according to the rule

(2) $k^{n+2}_{i}{}={}\left\\{\begin{array}[]{ll}\alpha
k^{n+1}_{i}&\dfrac{P^{n+1}s^{n}_{i}}{P^{n}b^{n}_{i}}>k^{n+1}_{i},\\\ \\\
k^{n+1}_{i}&\dfrac{P^{n+1}s^{n}_{i}}{P^{n}b^{n}_{i}}=k^{n+1}_{i},\\\ \\\ \beta
k^{n+1}_{i}&\dfrac{P^{n+1}s^{n}_{i}}{P^{n}b^{n}_{i}}<k^{n+1}_{i}.\end{array}\right.$

If the agent identifies a growing market she increases her stock-to-bond ratio
by a fixed amount, while the other agent decreases her ratio. The feedback
reflects the fact that when faced with a series of bad investments the agent
will reduce the equity part of the portfolio, while with the investment
growing better than expected, the agent will take riskier position.

It easy to see that there is an equilibrium when no shares are traded and
price doesn’t change: $s^{n}_{i}=s^{0}_{i},$ $b^{n}_{i}=b^{0}_{i},$
$k^{n}_{i}=s^{0}_{i}/b^{0}_{i},$ $P^{n}=P^{0}.$ Small deviations will set off
a non-trivial dynamics. In this process, agents will alternatively change
their preferences by factors $\alpha$ and $\beta,$ so that, on average, the
stock-to-bond ratios change at the rate $(\alpha\beta)^{1/2}$ per period. The
gross return, $P^{n+1}/P^{n},$ after a short transient period will settle at
the average value of $(\alpha\beta)^{1/2}$ per period, but will keep
oscillating around that value, as in Figure 1. The stock price goes up when
agents have bias toward the increase of the risk preferences, expressing the
optimistic outlook, quantified by the condition $\alpha\beta>1.$

The story of this process is simple: agents are trading stock and and move
their stock-to-bond ratios up and down with an upward average trend. This
drives up the average demand and consequently the price.

Figure 1. Gross return for first 25 trading periods in a 2-agent ASPP model,
with $\alpha=3.01,$ and $\beta=0.36.$ Circles indicate the values of gross
returns for alternating ups and downs of stock price, with a net return per
period approaching value $\sqrt{\alpha\beta}=1.012,$ (solid line).

###  2.2. Properties of $N$-agent ASPP

For a mulit-agent model, we populate the market with many copies of agents
considered above and let them trade anonymously, by submitting their buy-sell
orders to a market maker. Thus, we consider a set of $N$ agents, described by
their portfolios $(s_{i},b_{i}),$ $i=1\ldots N,$ of dollar values of a stock
and cash accounts, and let $k_{i}$ stand for agent $i$ stock-to-bond ratio,
either actual, or a future, target value. $P_{0}$ will denote a current price
per share and $P$ the new price determined by agents’ demand. Let
$\\{i_{l}\,:\,l=1\ldots m\\}$ be the set of “active” agents, i.e. the ones
setting the new price. The set of active traders is determined each trading
period by a random draw from the population.

If $x_{i_{l}}$ is the dollar amount that agent $i_{l}$ wants to invest in
stock, then

$\frac{\frac{P}{P_{0}}s_{i_{l}}+x_{i_{l}}}{b_{i_{l}}-x_{i_{l}}}{}={}k_{i_{l}}.$

The demand-supply balance is

$\sum_{l=1}^{m}x_{i_{l}}{}={}0,$

which can be solved for $P:$

$\frac{P}{P_{0}}{}={}\left(\sum_{l=1}^{m}\frac{k_{i_{l}}b_{i_{l}}}{1+k_{i_{l}}}\right)\left(\sum_{l=1}^{m}\frac{s_{i_{l}}}{1+k_{i_{l}}}\right)^{-1}.$

Once the price is set, agents move corresponding amounts between cash and
stock accounts, re-balancing their portfolios. The update mechanism for new
stock-to-bond ratios is given by (2). The interaction is repeated the
following trading periods with new, randomly selected sets of active agents.
In this trading scheme the total amount of cash in the system and number of
shares of stock remain constant.

There is an equilibrium state, when all agents have balanced portfolios,
$s_{i}/b_{i}=k_{i},$ $i=1,\ldots,N.$ In this case stock price doesn’t grow,
$P_{n}=P_{n-1}.$ When the parameters are out of the equilibrium, even by a
small degree, the system exhibits non-trivial, divergent dynamics, such as
shown in Figure 3.

In this divergent regime, after a certain transient period, the logarithm of
the price start fluctuating as normally distributed white noise with linearly
growing in time mean, see Perepelitsa and Timofeyev (2019). The distribution
of returns depends on the mechanism for selecting a random group of active
agents. If the number of active agents, $m,$ is fixed the returns in the
stationary regime follow log-normal distribution. If number $m$ is itself
random, for example chosen from an uniform distribution, the returns have
significantly more mass at the tails and at 1, than a corresponding log-normal
distribution.

To estimate the values of $\sigma$ in (4) we performed numerical simulation of
ASPP model, with $N=500$ and $m=10,20,40,$ on 200,000 trajectories for each
pair of test values of $(\alpha,\beta).$ The values of $\sigma$ for the case
with $m=10$ active traders are reported in Table 1. Figure 2 shows the scatter
plot of those values against the values of $\alpha-\beta$ and the
corresponding linear regression line. Linear regression coefficients are
listed in Table 2. It can be seen from that table that formula
$c(m)(\alpha-\beta)$ is a good approximation for volatility $\sigma,$ as a
function of $\alpha$ and $\beta.$ In all three cases, there is a positive
correlation of about $0.7$ between the mean log-return and the volatility.

$\alpha\,\backslash\,\beta$ | 0.8 | 0.825 | 0.85 | 0.875 | 0.9 | 0.925 | 0.95 | 0.975
---|---|---|---|---|---|---|---|---
1.1 | 0.0775 | 0.0679 | 0.0590 | 0.0508 | 0.0438 | 0.0380 | 0.0340 | 0.0318
1.2 | 0.0889 | 0.0817 | 0.0753 | 0.0700 | 0.0659 | 0.0629 | 0.0609 | 0.0596
1.3 | 0.1062 | 0.1066 | 0.0962 | 0.0927 | 0.0899 | 0.0878 | 0.0864 | 0.0854
1.4 | 0.1253 | 0.1211 | 0.1117 | 0.1150 | 0.1129 | 0.1112 | 0.1101 | 0.1093
1.5 | 0.1446 | 0.1413 | 0.1385 | 0.1363 | 0.1346 | 0.1333 | 0.1324 | 0.1318
Table 1. $\sigma$ as a function of $\alpha$ and $\beta$ for ASPP model with $N=500$ and $m=10,$ in a stationary regime reached after 2 years of trading time.  | c | d
---|---|---
m=10 | 0.2066 | 0.0081
m=20 | 0.1457 | 0.0055
m=40 | 0.1084 | 0.0027
Table 2. Linear regression coefficients $(c,d)$ in $\sigma=c(\alpha-\beta)+d$
for several values of $m.$ Figure 2.  Values of sigma from Table 1 plotted
against $\alpha-\beta.$

Combining these findings and assuming for the simplicity of presentation that
there is one trade per day (total 360 per year), and using $t_{n}=n/360$ to
label $n^{th}$ trading period on the time axis measured in years, we obtain
that daily prices $P(t_{n})$ are independently distributed with

(3) $\ln(P(t_{n}))\in\mathcal{N}(n\ln r,\,\sigma/\sqrt{2}),\quad
r{}={}(\alpha\beta)^{\frac{m}{2N}},\quad\sigma=c(m)(\alpha-\beta).$

In this formulas, $N$ is the total number of investors and $m$ is the number
of investors participating in one trade per day, and $c(m)$ is a positive
number, depending on $m.$ Accordingly, for the daily returns,
$R(t_{n})=P(t_{n})/P(t_{n-1}),$ we have

(4) $\ln(R(t_{n}))\in\mathcal{N}(\ln r,\,\sigma)$

where $r$ and $\sigma$ as in (3).

Let’s illustrate this behavior on a particular model setup, in which there is
total of $N=500$ investors, with $m=10$ randomly chosen agents to participate
in a daily trade, and the value of optimism is set at $\alpha=1.2$ and
pessimism at $\beta=0.96$ Figure 3 illustrates typical characteristic features
of a stochastic price pump. The top panel shows time series of the logarithm
of the price of a share with a “healthy” growth rate $1.66$ per year (66%
increase). The panel below shows that the returns are not constant but
fluctuate as a white noise with volatility of about $0.07$ (7% daily price
changes). This price growth is achieved by a modest amount of trading, with
daily trading volume not exceeding $0.3\%$ of the total cash in the system,
see the third panel from the top in Figure 3. Combined with the exponential
growth of the price, this means that the daily trading volume measured in the
units of the shares progressively decreases, see the bottom panel in Figure 3.
If, initially, trades evolved exchanges of amount of shares of order 1, then
by year 10, the typical daily volume is of order $10^{-2}$ and by year 20, of
order $10^{-4}.$ This property is a consequence of the exponential growth of
the price combined with the fact that the trading conserves the amount of cash
in the system.

This example shows that investors acting adaptively can generate a long-
lasting, positive price dynamics, and a necessary condition for an ASPP to
function for prolonged time is an infinite divisibility of traded commodity,
which is one of the main characteristic properties of digital currencies.

Our main interest in this paper will be the systematic risk of ASPP and its
relation to the volatility, which we discuss in the next section.

Figure 3. Four typical time series describing one simulation of ASPP for a
period of 20 years. From top to bottom: logarithm of the stock price,
logarithm of the daily gross returns, daily trading volume as percentage of
the total cash in the system, daily trading volume in units of shares.

##  3\. Systematic risk

Multi-agent models offer a convenient tool for estimating systematic risk of
sharp price drops (crashes), since the information about investment portfolio
of each agent is a part of the description of the model. Here, by systematic
risk we understand the risk of a breakdown of the system that generates the
price dynamics, in this case, the breakdown of ASPP. We consider two metrics
for measuring this risk. The first is based on the proportion of agents with
low cash-to-stock ratios. Let us denote by $q(t)$ the proportion of agents
with cash-to-stock ratio less than, say 10%, at time $t.$ Higher values of
$q(t)$ mean higher likelihood of a crash, as a larger proportion of agents may
spontaneously initiate stock sellout. We will use a hazard rate function

(5) $h_{q}(t){}={}\frac{q(t)}{1-q(t)}$

to quantify the probability of a crash. According to this formula, the hazard
rate changes monotonically from $0,$ when $q(t)=0,$ to $+\infty,$ when
$q(t)=1.$ Recall that the hazard rate function is defined so that for an
infinitesimal time interval of length $dt,$ $h_{q}(t)dt$ is the probability of
a crash to occur during interval $[t,t+dt],$ given that it didn’t occur before
time $t.$ Then, the probability that no crash occurs up to time $t,$ equals

(6) $S_{q}(t)=e^{-\int_{0}^{t}h_{q}(s)\,ds}.$

Thus, the “discounted” stock price at time $t$ equals $P(t)S_{q}(t),$
assuming, for simplicity, that a crash brings price to zero.

Alternatively, the systematic risk can be computed based off the proportion of
agents with low ratio of cash reserves at time $t$ to the initial cash
reserves, say less than 60%. We denote this proportion by $r(t),$ and the
corresponding hazard rate and the survival functions by $h_{r}(t)$ and
$S_{r}(t),$ defined by (5) and (6) where $h_{q}$ is replaced by $h_{r}.$

Figure 5 shows plots of the expected price, and the prices discounted by two
hazard rates, $h_{q}$ and $h_{r},$ for the same model as used in the
simulation in Figure 3. Measured by the first metric, the price collapses to
zero, while in the second metric the price continues to grow (eventually it
will go down as well).

It might seem that the first metric seems to be more relevant measure of a
sudden sellout. Note however, that a low cash-to-stock ratio is not caused by
traders spending all cash buying stock, but rather on the anomalous price
growth. For this reason the second metric is more likely to be used by
adaptive speculators with long investment horizon. As explained in the
previous section the latter is model of main interest in this paper.

Another difference between two metrics, that we would like to emphasis here,
is that the first depends on the initial investment portfolios, and can be
arbitrary high at time $t=0$ because of that. Whereas the second is always
zero at $t=0$ and increases gradually in the process of trading, thus being a
characteristic of a fully developed stationary regime of the price pump.

We can conclude that the persistence of the price bubbles in the digital
currencies may result from overestimation of the risk by either the hazard
rate $h_{q}(t),$ whereas the actual risk in an adaptive speculation model
given by the hazard rate $h_{r}(t)$ with significantly smaller values, or the
high volatility which is positively correlated with the return, see formulas
for $r$ and $\sigma$ in (3). The latter might be the reason why short term
investors stay away from them.

As can be seen from the previous examples, for fixed values of $\alpha$ and
$\beta,$ the volatility of the price is constant. At the same time, the hazard
function increases, implying that the volatility can not be a correct measure
of the systematic risk. However, they are not completely unrelated. To show
this, let us define a cumulative volatility at $t$ years as

$V(t)=\int_{0}^{t}\sigma^{2}(s)\,ds.$

For a time series of daily stock price $P(t_{n}),$ the cumulative volatility
$V(t_{n})=\sum_{t_{m}\leq t_{n}}[\ln R(t_{m})-\ln R(t_{m-1})]^{2}/720,$ where
$R(t_{n})=P(t_{n})/P(t_{n-1})$ is the gross return.

Figure 6 shows a scatter plot of daily values of $V(t_{n})$ and the hazard
rate $h_{r}(t_{n}),$ and the best fit linear approximation (linear
regression). This numerical evidence is strongly suggestive of a linear
relation between $V(t)$ and $h_{r}(t),$ or, equivalently, a linear relation

(7) $\sigma^{2}=\gamma\frac{dh_{r}}{dt}.$

Table 3 lists the statistical values of parameter $\gamma$ for several values
of $(\alpha,\beta).$ This formula provides an important practical tool for
estimation of the systematic risk from the volatility, since the latter can be
computed from time series of the price.

From this formula we can find the hazard rate $h_{t}(t)=\sigma^{2}t\gamma,$
and the PDF for the time of a crash. The latter has a Weibull distribution

$\frac{\sigma^{2}t}{\gamma}{\rm
exp}\left({-\frac{\sigma^{2}}{2\gamma}t^{2}}\right).$

Whith that, the expected discounted stock price becomes equal

$P_{d}(t)=P(0){\rm
exp}\left(r_{0}t-\frac{\sigma^{2}}{2\gamma}t^{2}\right),\quad
r_{0}=\frac{360m}{2N}\ln(\alpha\beta).$

It can be seen that there is a maximum value $\max P_{d}(t)$ of the price,
attained when $t=2r_{0}\gamma/\sigma^{2}.$ It places a limit on the price that
a bubble can reach in a fixed “psycho-state” $(\alpha,\beta).$

Figure 4. Proportions of traders with low cash-to-stock portfolios (left) and low cash-to-initial-cash ratio (right) as functions of time in one simulation of the model for the period of 20 years. Figure 5. Expected stock price discounted by the risk of a crash. Top graph is the expected price without discounting. The graph in the middle is the expected price discounted by hazard rate $h_{r}(t).$ The graph at the bottom is the expected price discounted by hazard rate $h_{q}(t).$ Figure 6. Scatter plot for cumulative volatility $V(t_{n})$ and hazard rate $h_{r}(t_{n}),$ in one simulation of the model for the first 20 years. $\alpha\,\backslash\,\beta$ | 1 | 0.99 | 0.96
---|---|---|---
1.1 | 0.2909 | 0.3932 | 0.8343
1.2 | 0.2591 | 0.2763 | 0.3577
1.3 | 0.3644 | 0.3732 | 0.4084
Table 3. Values of $\gamma$ in formula (7). Each value is an average of the
linear regression coefficients of $V(t_{n})$ on $h_{r}(t_{n}),$ for 1000
20-year path simulations of ASPP model with $N=500$ agents in the system,
$m=10$ active agents.

##  4\. Conclusions

In the last section we listed several factors that can contribute to a
prolonged life span of the price bubble. But will the bubble perpetuate its
growth forever? What will happen in an ideal case when no crashes due to
systematic risk occur, when there are no arbitrageurs out there, and investors
have very long, or even infinite, investment horizon? The answer to these
questions lies in the interpretation of parameters $(\alpha,\beta)$ of the
model, that we identify as two psychological motivators, greed and fear,
receptively. Here, we are assuming that a coherent psychological state in the
population is achieved through the copying behavior and inter-connectedness of
the investors by information networks, so that all investors share the same or
almost the same values of $\alpha$ and $\beta.$ These parameters are combined
together in the formula for the return (3), which shows that the price pump
can work both ways, as a price inflator or a deflator.

Being of psychological nature $(\alpha,\beta)$ will change according to some
intrinsic mechanisms, that combine personal psychological response of each
investor to the changes in the composition of his/her portfolio and the degree
of influence of this response on the behavior of other investors. As an
approximation, one can assume that $(\alpha,\beta)$ co-evolve in such a way
that $\alpha$ (greed) is self-reinforcing and generates $\beta$ (fear), and
fear, in its turn, inhibits greed, by increasing investor’s sensitivity to
idiosyncratic noise.

According to this model, the evolution of a price bubble will proceed in a
“predator-prey” type cycle in the domain of psychological variables, triggered
by an exogenous shock. Along the cycle, the return and volatility of the price
will changes according to formulas from (3), producing a pattern of increasing
returns with low volatility, followed by a period of growing volatility,
during which the return reaches its maximum and starts to decrease, which
eventually ends with sharp decrease on return and low volatility. At the end
of the cycle the system will reach some “background level” of
$(\alpha,\beta),$ at which they balance each other. The system will remain in
this “no bubble” phase until another external factor pushes it on a road to a
new bubble.

This type of dynamics can produce complex patterns of price variations. Since
the period of a bubble cycle is a non-constant (random) variable, there can be
little predictability of the future price, even more so with crashes due to
systematic risk. This suggests an intriguing idea that the price fluctuations
due to adaptive speculative market behavior might be consistent with the no-
arbitrage principle.

## References

  * [1] Arthur, W.B., Holland, J.H., LeBaron, B, Palmer, R.G., Tayler, P. (1997). Asset Pricing under Endogenous Expectations in an Artificial Stock Market. in W.B. Arthur, S. Durlauf and D. Lane. eds. The Economy as an Evolving Complex System II. Addison-Wesley.
  * [2] Blanchard. O.J., Watson, M.W. (1982). Bubbles, rational expectations and speculative markets. in Crisis in Economic and Financial Structure: Bubbles, Bursts and Shocks, P. Wachtel, editor, Lexington Books, Lexington, MA.
  * [3] Cross, J.G. (1983). A theory of adaptive economic behavior. Cambridge Univ. Press.
  * [4] Dodd, N. (2017). The social life of bitcoin. Theory, Culture & Society, 35 (3). 35–56.
  * [5] Fudenberg, D., Levine, D.K. (1998). The theory of Learning in Games. The MIT Press.
  * [6] Griffin, J.M., Shams, A. (2020). Is Bitcoin Really Untethered? J. of Finance, vol. LXXV, no. 4.
  * [7] Levy, M., Levy, H., Solomon, S. (1994). A microscopic model of the stock market: Cycles, booms, and crashes. Economics Letters, 45:103–111.
  * [8] Levy, M., Levy, H., Solomon, S. (1995). Microscopic simulation of the stock market: The effect of microscopic diversity. Journal Physique I, 5:1087–1107.
  * [9] Levy, M., Levy, H., Solomon, S. (1997). New evidence for the power law distribution of wealth. Physica A, 242:90–94.
  * [10] Levy, M., Levy, H., Solomon, S. (2000) Microscopic Simulation of Financial Markets: From Investor Behavior to Market Phenomena. Academic Press.
  * [11] Lux, T. (1995). Herd behavior, bubbles and crashes. Economic Journal, 105:881–896.
  * [12] Lux, T. (1998). The socio-economic dynamics of speculative markets: Interacting agents, chaos, and the fat tails of return distributions. Journal of Economic Behavior and Organization, 33:143–165.
  * [13] Lux, T., M. Marchesi, M. (1999). Scaling and criticality in a stochastic multi-agent model of a financial market. Nature, 397:498–500, 1999.
  * [14] Lux, T., Marchesi, M. (2000). Volatility clustering in financial markets: A micro- simulation of interacting agents.International Journal of Theoretical and Applied Finance, 3:67–702.
  * [15] Kim., G., Markowitz, H.M. (1989). Investment rules, margin and market volatility. Journal of Portfolio Management. 16:45–52.
  * [16] Mandelbrot, B. (1966). Nonlinear forecasts, rational bubbles, and martingales. Journal of Business 39, 242–255.
  * [17] Perepelitsa, M., Timofeyev, I. (2019) Asynchronous stochastic price pump, Physica A: J. Stat. Mech. Appl. 516, (2019), 356–364.
  * [18] Roth, A.E., Erev, I. (1995). Learning in extensive-form games: experimental data and simple dynamics models in the intermediate term. Games and Economic Behavior, 8 164–212.
  * [19] Sornette, D. (2003). Why stock markets crash: critical events in complex financial systems. Princeton University Press.
  * [20] Stigler, G.J. (1964). Public regulation of the securities market. Journal of Business, 37.
  * [21] Wolff-Mann, E. (2018). Only good for drug dealers: More Nobel prize winners snub bitcoin. Yahoo Finance. https://web.archive.org/web/20180612141128if_/https://finance.yahoo.com/news/good-drug-dealers-nobel-prize-winners-snub-bitcoin-184903784.html