Murray's law for discrete and continuum models of biological networks

TL;DR

This paper validates Murray's law, a key scaling relation for biological network conductivities, across discrete and continuum models, including generalized and steady-state cases, enhancing understanding of biological transport systems.

Contribution

It generalizes Murray's law for discrete networks with multiple branches and metabolic coefficients, and proves its validity for continuum models derived from discrete limits and phenomenological assumptions.

Findings

Generalized 3/4-law for discrete networks with multiple branches

Proved Murray's law for continuum models from discrete limits

Validated Murray's law for stable steady states in phenomenological models

Abstract

We demonstrate the validity of Murray's law, which represents a scaling relation for branch conductivities in a transportation network, for discrete and continuum models of biological networks. We first consider discrete networks with general metabolic coefficient and multiple branching nodes and derive a generalization of the classical 3/4-law. Next we prove an analogue of the discrete Murray's law for the continuum system obtained in the continuum limit of the discrete model on a rectangular mesh. Finally, we consider a continuum model derived from phenomenological considerations and show the validity of the Murray's law for its linearly stable steady states.

Equations

$$E_p=Q\Delta P,$$

$$Q=\frac{\pi R^4}{8\mu L}\Delta P,$$

$$E_p=\frac{8\mu L}{\pi R^4}{Q}^2.$$

$$E_m=\nu \pi L{R}^2,$$

$$E:=E_p+E_m=\frac{8\mu L}{\pi R^4}{Q}^2+\nu \pi L{R}^2.$$

$$Q=\sqrt{\frac{{\pi }^2\nu }{16\mu }}{R}^3.$$

$$Q_0=\sum _{i=1}^n{Q}_i.$$

$$R_0^3=\sum _{i=1}^n{R}_i^3,$$

$$C_0^{3/4}=\sum _{i=1}^n{C}_i^{3/4}.$$

$$(\Delta P{)}_{ij}:=P_j-P_i.$$

$$Q_{ij}:=C_{ij}\frac{P_j-P_i}{L_{ij}}\text{for all }(i,j)\in E.$$

$$-\sum _{j\in \mathcal{N}(i)}{C}_{ij}\frac{P_j-P_i}{L_{ij}}=S_i\text{for all }i\in V.$$

$$\mathcal{N}(i):=\{j\in V;(i,j)\in E\},$$

$$\sum _{i\in V}{S}_i=0,$$

$$\mathcal{E}[C]:=\sum _{(i,j)\in E}\left(\frac{Q_{ij}[C{]}^2}{C_{ij}}+\frac{\nu }{\gamma }{C}_{ij}^{\gamma }\right)L_{ij}.$$

$$\frac{\partial }{\partial C_{kl}}\sum _{(i,j)\in E}\frac{Q_{ij}[C{]}^2}{C_{ij}}{L}_{ij}=-\frac{Q_{kl}[C{]}^2}{C_{kl}^2}{L}_{kl}.$$

$$\frac{\partial }{\partial C_{kl}}\sum _{(i,j)\in E}\frac{Q_{ij}[C{]}^2}{C_{ij}}{L}_{ij}=-\frac{Q_{kl}[C{]}^2}{C_{kl}^2}{L}_{kl}+2\sum _{(i,j)\in E}\frac{Q_{ij}[C]}{C_{ij}}\frac{\partial Q_{ij}[C]}{\partial C_{kl}}{L}_{ij},$$

$$\sum _{(i,j)\in E}\frac{Q_{ij}[C]}{C_{ij}}\frac{\partial Q_{ij}[C]}{\partial C_{kl}}{L}_{ij}=0.$$

$${\mathbb{A}}_{ij}=\{\begin{array}{ll}0& \text{mbox}if(i,j)\notin E,\\ 1& \text{mbox}if(i,j)\in E.\end{array}$$

$$2\sum _{(i,j)\in E}\frac{Q_{ij}[C]}{C_{ij}}\frac{\partial Q_{ij}[C]}{\partial C_{kl}}{L}_{ij}$$

$$=\sum _{j=1}^n{P}_j\sum _{i=1}^n{\mathbb{A}}_{ij}\frac{\partial Q_{ij}[C]}{\partial C_{kl}}-\sum _{i=1}^n{P}_i\sum _{j=1}^n{\mathbb{A}}_{ij}\frac{\partial Q_{ij}[C]}{\partial C_{kl}}$$

$$=-2\sum _{i=1}^n{P}_i\sum _{j=1}^n{\mathbb{A}}_{ij}\frac{\partial Q_{ij}[C]}{\partial C_{kl}}$$

$$=-2\sum _{i=1}^n{P}_i\frac{\partial }{\partial C_{kl}}\sum _{j\in N(i)}{Q}_{ij}.$$

$$-\sum _{j\in N(i)}{Q}_{ij}=\sum _{j\in N(i)}{C}_{ij}\frac{P_j-P_i}{L_{ij}}=S_i,$$

$$\sum _{(i,j)\in E}\frac{Q_{ij}[C]}{C_{ij}}\frac{\partial Q_{ij}[C]}{\partial C_{kl}}{L}_{ij}=0.$$

$$\frac{\partial }{\partial C_{ij}}\mathcal{E}[C]=-\left(\frac{Q_{ij}[C{]}^2}{C_{ij}^2}-\nu C_{ij}^{\gamma -1}\right)L_{ij}.$$

$$\frac{\partial }{\partial C_{ij}}\mathcal{E}[C]=0\text{mbox}forall(i,j)\in E,$$

$$Q_{ij}[C{]}^2=\nu C_{ij}^{\gamma +1}\text{mbox}forall(i,j)\in E.$$

$$\sum _{j\in \mathcal{N}(i)}{Q}_{ij}+S_i=0\text{for all }i\in V.$$

$${\mathcal{N}}^{+}(i):=\{j\in \mathcal{N}(i);Q_{ij}\ge 0\},{\mathcal{N}}^{-}(i):=\{j\in \mathcal{N}(i);Q_{ij}<0\}.$$

$$\sum _{j\in {\mathcal{N}}^{+}(i)}{Q}_{ij}+S_i=\sum _{j\in {\mathcal{N}}^{-}(i)}(-Q_{ij})\text{for all }i\in V,$$

$$\sum _{j\in {\mathcal{N}}^{+}(i)}|Q_{ij}|+S_i=\sum _{j\in {\mathcal{N}}^{-}(i)}|Q_{ij}|\text{for all }i\in V.$$

Full text

Murray’s law for discrete and continuum models of biological networks

Jan Haskovec111Mathematical and Computer Sciences and Engineering Division, King Abdullah University of Science and Technology, Thuwal 23955-6900, Kingdom of Saudi Arabia; jan.haskovec@kaust.edu.sa  Peter Markowich222Mathematical and Computer Sciences and Engineering Division, King Abdullah University of Science and Technology, Thuwal 23955-6900, Kingdom of Saudi Arabia; peter.markowich@kaust.edu.sa, and Faculty of Mathematics, University of Vienna, Oskar-Morgenstern-Platz 1, 1090 Vienna; peter.markowich@univie.ac.at  Giulia Pilli333Faculty of Mathematics, University of Vienna, Oskar-Morgenstern-Platz 1, 1090 Vienna; giulia.pilli@univie.ac.at

Abstract: We demonstrate the validity of Murray’s law, which represents a scaling relation for branch conductivities in a transportation network, for discrete and continuum models of biological networks. We first consider discrete networks with general metabolic coefficient and multiple branching nodes and derive a generalization of the classical 3/4-law. Next we prove an analogue of the discrete Murray’s law for the continuum system obtained in the continuum limit of the discrete model on a rectangular mesh. Finally, we consider a continuum model derived from phenomenological considerations and show the validity of the Murray’s law for its linearly stable steady states.

Keywords: Biological transportation networks; Murray’s law; Continuum limit.

AMS Subject Classification: 92C35, 05C21, 76S05

1. Introduction

Murray’s law is a basic physical principle for transportation networks which predicts the thickness or conductivity of branches, such that the cost for transport and maintenance of the transport medium is minimized. This law is observed in the vascular and respiratory systems of animals, xylem in plants, and the respiratory system of insects[18, 16]. It is also a powerful biomimetics design tool in engineering and has been applied in the design of microfluidic devices, self-healing materials, batteries, photocatalysts, and gas sensors; see, e.g., Refs. [3, 5, 19, 21].

Murray’s original analysis [14, 15] is based on the assumption that the radii inside a lumen-based system are such that the work for transport and upkeep is minimized. Larger vessels lower the expended energy for transport (pumping power), but increase the overall volume of blood in the system, which requires metabolic support. Murray’s law reflects the optimum balancing of the pumping power versus metabolic expenditure. Assuming laminar flow in a single tube, the pumping power $E_{p}$ is given by the Joule’s law,

[TABLE]

where $Q$ is the volumetric flux and $\Delta P$ is the pressure difference between the entry and exit of the tube. According to the Hagen-Poiseuille law for laminar flow in a three-dimensional tube,

[TABLE]

where $R>0$ is the radius of the tube, $L>0$ its length, and $\mu>0$ the dynamic viscosity of the fluid. Therefore, we have

[TABLE]

Furthermore, we assume that the metabolic cost $E_{m}$ is proportional to the number of blood cells in the tube[14], and that is proportional to the tube’s volume,

[TABLE]

where the proportionality constant $\nu>0$ is called the metabolic coefficient. This gives for the total energy

[TABLE]

A trivial calculation shows that the total energy is minimized if the following relation between the flow rate and the tube radius is verified,

[TABLE]

In a branching node with a parent branch with flow rate $Q_{0}>0$ and children branches with flow rates $Q_{1},Q_{2},\dots,Q_{n}>0$, we have the flux balance

[TABLE]

Using (3), we obtain for the respective radii $R_{0}>0$ and $R_{1},R_{2},\dots,R_{n}>0$,

[TABLE]

which is the classical Murray’s law for laminar blood flow[14, 15]. It can be also expressed in terms of the edge conductivities $C_{i}>0$, which for three-dimensional laminar flow are proportional to $R_{i}^{4}$, see Ref. [17], therefore

[TABLE]

The above formula represents the classical 3/4-law for optimal ratio of conductivities in 3D branching edges with laminar flow [14, 15]. It clearly demonstrates the nonlinear nature of optimal network conductivity models, in particular, the nonlinear (and nonconvex) dependence of the total energy (2) on the edge conductivity.

The first goal of this paper is to show the validity of a generalization of the Murray’s law (5) for general transportation networks with multiple branching nodes (Section 2). Then, in Section 3, we consider a continuum limit of the discrete model on a rectangular mesh, and derive an analogue of the discrete Murray’s law for the continuum PDE system. Finally, in Section 4 we consider another continuum model for biological transportation networks, derived from phenomenological considerations, and show that a continuum version of Murray’s law holds also for the linearly stable solutions of this model. The calculations made in Sections 3 and 4 are formal, however, they can be justified by assuming sufficient regularity of the steady state solutions of the respective continuum PDE systems. Details will be provided in the respective Sections.

2. Murray’s law for the discrete network model

We consider a given undirected connected graph $G=(V,E)$, consisting of a finite set of vertices (branching nodes) $V$ of size $N=|V|$ and a finite set of edges $E$. Any pair of vertices is connected by at most one edge, which represents a tube (channel), and no vertex is connected to itself. We denote the edge between vertices $i\in V$ and $j\in V$ by $(i,j)\in E$. Since the graph is undirected, $(i,j)$ and $(j,i)$ refer to the same edge. For each edge $(i,j)\in E$ of the graph $G$ we consider its length and its conductivity, denoted by $L_{ij}=L_{ji}>0$ and $C_{ij}=C_{ji}\geq 0$, respectively. The edge lengths $L_{ij}>0$ are given as a datum and fixed for all $(i,j)\in E$. With each vertex $i\in V$ there is associated the fluid pressure $P_{i}\in\mathbb{R}$. The pressure drop between vertices $i\in V$ and $j\in V$ connected by an edge $(i,j)\in E$ is denoted by

[TABLE]

Note that the pressure drop is antisymmetric, i.e., by definition, $(\Delta P)_{ij}=-(\Delta P)_{ji}$. The oriented flux (flow rate) from vertex $j\in V$ to $i\in V$ is denoted by $Q_{ij}$; again, we have $Q_{ij}=-Q_{ji}$. For the three-dimensional laminar flow, conductance (transportation capacity) $C_{ij}$ of an edge $(i,j)\in E$ is proportional to the fourth power of its diameter[17]. Then, (1) is reformulated in terms of the conductivities as

[TABLE]

The flux balance in each vertex is expressed in terms of the Kirchhoff law

[TABLE]

Here $\mathcal{N}(i)\subset V$ denotes the set of vertices connected to $i\in V$ through an edge, i.e.,

[TABLE]

and $S=(S_{i})_{i\in V}$ is the prescribed strength of the external flow source ($S_{i}>0$) or sink ($S_{i}<0$) in node $i\in V$. Clearly, a necessary condition for the solvability of (8) is the global mass conservation

[TABLE]

which we shall assume in the sequel. Given the vector of conductivities $C=(C_{ij})_{(i,j)\in E}$, the Kirchhoff law (8) is a linear system of equations for the vector of pressures $P=(P_{i})_{i\in V}$. With the global mass conservation (9), the linear system (8) is solvable if and only if the graph with edge weights $C=(C_{ij})_{(i,j)\in E}$ is connected[2], where only edges with positive conductivities $C_{ij}>0$ are taken into account (i.e., edges with zero conductivities are discarded). Note that the solution vector $P=(P_{i})_{i\in V}$ is unique up to an additive constant. The corresponding fluxes, which we denote by $Q_{ij}[C]$, are then calculated uniquely from (7).

The conductivities $C_{ij}$ are subject to an energy optimization and adaptation process, where, in analogy to (2), the cost functional is a sum of pumping and metabolic energies over all edges of the graph,

[TABLE]

Here we consider a generalization where the metabolic energy is proportional to a power $\gamma>0$ of the respective edge conductivity. We have $\gamma=1/2$ for blood vessel systems (recall that in a three-dimensional tube the conductivity is proportional to fourth power of the tube radius, while the metabolic cost scales quadratically with the radius, see Ref. [14]). For models of plant leaf venation the metabolic cost is proportional to the number of small tubes, which is proportional to $C_{ij}$, and the metabolic cost is due to the effective loss of the photosynthetic power at the area of the venation cells, which is proportional to $C_{ij}^{1/2}$. Consequently, the effective value of $\gamma$ typically used in models of leaf venation in plants lies between $1/2$ and $1$, see, e.g., Refs. [10, 11]. Again, $\nu>0$ is the so-called metabolic coefficient.

The derivation of the Murray’s law for critical points of the energy functional (10), constrained by the Kirchhoff law (8), is based on the following explicit formula for the derivative of the pumping term in (10) with respect to the conductivities, derived in Ref. [12]. We restate its proof here for the sake of completeness.

Lemma 1**.**

*Let $C=(C_{ij})_{(i,j)\in E}$ be a vector of nonnegative conductivities such that the underlying graph is connected if only edges with strictly positive conductivities are taken into account. Let $Q_{ij}[C]$ be given by (7), where $P=(P_{i})_{i\in V}$ is a solution of the linear system (8) with the vector of conductivities $C$. Then, for any fixed $(k,l)\in E$ we have *

[TABLE]

Proof.

Since

[TABLE]

it is sufficient to show that

[TABLE]

Let $\mathbb{A}=(\mathbb{A}_{ij})$ denote the adjacency matrix of the graph $G=(V,E)$, i.e. its coefficients are defined by

[TABLE]

Note that $G$ is an undirected graph, implying $\mathbb{A}_{ij}=\mathbb{A}_{ji}$. Due to the symmetry of $C_{ij}$ and $L_{ij}$ and antisymmetry of $Q_{ij}$ we have

[TABLE]

By the definition of the flow rate $Q_{ij}$ in (7) and Kirchhoff’s law (8) we have

[TABLE]

and since the sources/sinks $S_{i}$ are fixed, we conclude

[TABLE]

∎

Using the identity (11), it is straightforward to calculate the derivative of the energy (10), constrained by the Kirchhoff law (8), with respect to the edge conductivities,

[TABLE]

Since any critical point of ${\mathcal{E}}[C]$ satisfies

[TABLE]

we have

[TABLE]

Now, the Kirchhoff law (8) written in terms of the fluxes $Q_{ij}:=Q_{ij}[C]$ reads

[TABLE]

Let us write the set of neighbors $\mathcal{N}(i)$ of the node $i\in V$ as the disjoint union $\mathcal{N}(i)=\mathcal{N}^{+}(i)\cup\mathcal{N}^{-}(i)$, based on the flow directions along the respective edges, i.e.,

[TABLE]

Then, we can rewrite (14) as

[TABLE]

and, further,

[TABLE]

Using (13), we arrive at

[TABLE]

This is the (generalized) Murray’s law for the discrete model (8), (10). Recall that for the case of three-dimensional blood vessel systems we have $\gamma=1/2$, which leads to the power $\frac{\gamma+1}{2}=\frac{3}{4}$, so that (15) becomes a direct generalization of the classical Murray’s law (5).

3. Murray’s law for the continuum limit model on rectangular grids

In Ref. [12] the formal continuum limit of of the discrete model (8), (10) for $\gamma>1$ on rectangular equidistant grids was derived; the rigorous limit passage was studied in the consequent paper[13]. The resulting continuum energy functional is of the form

[TABLE]

with the metabolic coefficient $\nu>0$ and the exponent $\gamma>1$. The energy functional $\mathcal{E}[c]$ is defined on the set of nonnegative diagonal tensor fields $c=c(x)$ on a bounded domain $\Omega\subset\mathbb{R}^{d}$, $d\in\mathbb{N}$,

[TABLE]

The symbol $|c|^{\gamma}$ is defined as $\sum_{k=1}^{d}\left|c^{k}\right|^{\gamma}$. The scalar pressure $p=p(x)$ of the fluid occupying the domain $\Omega$ is subject to the Poisson equation

[TABLE]

equipped with no-flux boundary condition on $\partial\Omega$, and the datum $S=S(x)$ represents the intensity of sources and sinks in $\Omega$. Let us point out that the Poisson equation (18) is possibly strongly degenerate since the eigenvalues (i.e., diagonal elements) of the permeability tensor (17) may vanish. To overcome this problem, we follow Ref. [12] and introduce a regularization of (18) of the form

[TABLE]

with the permeability tensor

[TABLE]

where $r=r(x)\geq r_{0}>0$ is a prescribed function that models the isotropic background permeability of the medium, and $I\in\mathbb{R}^{d\times d}$ is the unit matrix. Again, (19) is equipped with the no-flux boundary condition

[TABLE]

where $n=n(x)$ denotes the outer unit normal vector to $\partial\Omega$. Then, the weak formulation of the boundary value problem (19)–(21) with a test function $\phi\in C^{\infty}(\Omega)$ reads

[TABLE]

Clearly, since we only admit nonnegative diagonal tensor fields $c=c(x)$, (22) is uniformly elliptic and it has solutions unique up to an additive constant. The energy functional (16) needs to be updated as

[TABLE]

We now calculate the Euler-Lagrange equations corresponding to critical points of the functional (23) constrained by (22).

Lemma 2**.**

*The Euler-Lagrange equations for the constrained energy minimization problem (22), (23) read *

[TABLE]

Proof.

We calculate the first variation of $\mathcal{E}$ in the direction $\phi$ where $\phi$ denotes a diagonal matrix with entries $\phi^{1},\ldots,\phi^{d}$, such that $c+\varepsilon\phi$ is nonnegative. Using the expansion

[TABLE]

we have

[TABLE]

Using $p_{0}$ as a test function in (22) with the permeability tensor $\mathbb{P}[c+\varepsilon\phi]=rI+c+\varepsilon\phi$ gives

[TABLE]

Subtracting the identity

[TABLE]

we obtain

[TABLE]

Plugging this into (26) gives

[TABLE]

and (24) follows. ∎

From (24) it follows that for $k=1,\dots,d$ there exist disjoint Lebesgue-measurable sets $\mathcal{A}_{k}^{+}$, $\mathcal{A}_{k}^{-}$ such that $\Omega={\mathcal{A}_{k}^{+}}\cup{\mathcal{A}_{k}^{-}}$ for all $k=1,\dots,d$ and

[TABLE]

Note that by assumption $c_{k}\geq 0$. Let us now choose an open subset $\Lambda\subset\Omega$ with boundary $\partial\Lambda\in C^{0,1}$ and use the characteristic function of $\Lambda$ as the test function $\phi$ in (22). This formal calculation can be justified by using a smoothened version of the characteristic function and passing to the limit, adopting the assumption that the flux $q:=-(rI+c)\nabla p$ is locally Lipschitz continuous on $\Lambda$. Then the generalized Green’s formula by De Giorgi-Federer[7] gives

[TABLE]

where $n=n(x)$ denotes the outward unit normal vector to $\partial\Lambda$, and $s=s(x)$ is the $(d-1)$-dimensional measure on $\partial\Lambda$. Using (28), we arrive at

[TABLE]

We define the inflow and outflow segments of $\partial\Lambda$ as

[TABLE]

with the flux $q=-(rI+c)\nabla p$. Then, (30) can be written in the form

[TABLE]

This is the Murray’s law for the continuum model (19)–(23).The two terms on the left-hand side represent the difference of in- and outgoing fluxes through the boundary $\partial\Lambda$, which is balanced by the volume integral of the external in- and outflows $S=S(x)$. Note the analogy with the Murray’s law for the discrete system (15).

Let us remark that the above formal calculation becomes rigorous if we assume sufficient regularity of the flux $q=-(rI+c)\nabla p$, in particular, we need $q$ to be locally Lipschitz continuous on $\Lambda$. In the next Section we construct solutions in the energy space which are less regular; the proof of higher regularity is beyond the scope of this paper and we postpone it to a future work.

3.1. Construction of solutions

From (24) it follows that

[TABLE]

Inserting this into (19), we obtain

[TABLE]

The weak formulation of this nonlinear elliptic equation, equipped with a no-flux boundary condition, reads

[TABLE]

for all test functions $\phi\in C^{\infty}(\Omega)$. For $\gamma>1$ solutions can be constructed by the direct method of calculus of variations.

Lemma 3**.**

Let $\gamma>1$. For every $S\in L^{2}(\Omega)$ there exists a unique solution $p\in H^{1}(\Omega)$ of (33) satisfying $\int_{\Omega}p(x)\,\mathrm{d}x=0$.

Proof.

We define the functional $\mathcal{F}:H^{1}(\Omega)\to\mathbb{R}$ by

[TABLE]

Since $\frac{2\gamma}{\gamma-1}>2$ for $\gamma>1$, the functional $\mathcal{F}$ is uniformly convex on $H^{1}(\Omega)$. Moreover, a straightforward application of the Poincaré inequality provides coercivity on the set $H^{1,0}(\Omega):=\left\{p\in H^{1}(\Omega);\;\int_{\Omega}p(x)\,\mathrm{d}x=0\right\}$. The classical theory (see, e.g., Ref. [6]) provides the existence of a unique minimizer $p\in H^{1,0}(\Omega)$. Since, as is easy to check, (34) is the Euler-Lagrange equation for the minimizer of $\mathcal{F}$, $p$ is the unique solution in $H^{1,0}(\Omega)$ of (33). ∎

Let us note that (33) is independent of the choice of the sets $\mathcal{A}_{k}^{+}$, $\mathcal{A}_{k}^{-}$ in (28), and so is the solution $p$ constructed in Lemma 3. We can therefore write, for $k=1,\dots,d$,

[TABLE]

3.2. Murray law for the continuum model with diffusion

Following Ref. [12], we equip the model with a linear diffusive term that accounts for random fluctuations in the medium. In particular, we update the energy functional (23) to

[TABLE]

where $D^{2}>0$ is the diffusivity constant and the symbol $|\nabla c|^{2}$ is defined as $\sum_{k=1}^{d}|\nabla c^{k}|^{2}$. The permeability tensor is given by (20), i.e., $\mathbb{P}[c]=rI+c$. As can be easily checked, the Euler-Lagrange equations corresponding to critical points of the functional (35) constrained by the Poisson equation (22) read

[TABLE]

Note that for sufficiently regular critical points of (35), (22) we have $\nu|c^{k}|^{\gamma-1}-D^{2}\Delta c^{k}\geq 0$ almost everywhere on $\Omega$. Then it follows that for $k=1,\dots,d$ there exist disjoint Lebesgue-measurable sets $\mathcal{A}_{k}^{+}$, $\mathcal{A}_{k}^{-}$ such that $\Omega={\mathcal{A}_{k}^{+}}\cup{\mathcal{A}_{k}^{-}}$ for all $k=1,\dots,d$ and

[TABLE]

As before we choose an open subset $\Lambda\subset\Omega$ with boundary $\partial\Lambda\in C^{0,1}$, and using (37) in (29), we obtain

[TABLE]

where $n=n(x)$ denotes the outward unit normal vector to $\partial\Lambda$, and $s=s(x)$ is the $(d-1)$-dimensional measure on $\partial\Lambda$. We define the inflow and outflow segments of $\partial\Lambda$ as

[TABLE]

where we denoted the flux $q:=-(rI+c)\nabla p$. Then, (38) can be written in the form

[TABLE]

This is the Murray’s law for the continuum model (35) constrained by the Poisson equation (22).

4. Murray’s law for the phenomenological continuum model

In this Section we consider the continuum model briefly proposed in Ref. [10] and analyzed in Refs. [8, 9, 1]. Its phenomenological derivation, carried out in Ref. [2], assumes that the network domain $\Omega\subset\mathbb{R}^{d}$ is occupied by a porous medium with the permeability tensor

[TABLE]

with $m=m(x)\in\mathbb{R}^{d}$ a smooth vector field on $\mathbb{R}^{d}$. Note that $\mathbb{P}[m]$ has the eigenvalue $|m|^{2}$ with eigenvector $m$, and [math] with eigenvectors orthogonal to $m$. Thus, it represents conduction along the direction $m$ with conductivity $|m|^{2}$, while there is no conduction in directions perpendicular to $m$. Consequently, we locally identify the network conductivity with the principal eigenvalue of $\mathbb{P}$, i.e., $C_{i}\simeq|m|^{2}$ in the neighborhood of a node $i\in V$. Assuming quasi-incompressibility of the fluid, i.e., constant fluid density along particle trajectories, we have the local mass conservation law

[TABLE]

where $q=q(x)\in\mathbb{R}^{d}$ is the flux and $S=S(x)$ is the density of external sources and sinks; see Ref. [2] for details. Assuming the validity of Darcy’s law for slow flow in porous media[20], we have

[TABLE]

where $p=p(x)$ is the fluid pressure. Combining with (40) and (39), we arrive at the Poisson equation

[TABLE]

Let us point out the possible strong degeneracy of (41), leading to solvability issues. That is why in Refs. [8, 9, 1] the following regularization was considered,

[TABLE]

where $I$ is the identity matrix and the scalar function $r=r(x)\geq r_{0}>0$ describes the isotropic background permeability of the medium. Then, (42) is uniformly elliptic and, equipped with suitable boundary conditions, solvable. We shall assume no flux through the boundary $\partial\Omega$, i.e., equip (42) with the homogeneous Neumann boundary condition

[TABLE]

where $n=n(x)$ denotes the outer unit normal vector to $\partial\Omega$. Then the solution $p=p(x)$ is unique up to an additive constant and the flux $q=-\mathbb{P}[m]\nabla p$ is defined uniquely.

The local conductivity of the network being identified with $|m|^{2}$, the continuum analogue of the discrete energy functional (10) is

[TABLE]

where $p[m]$ denotes a solution of the Poisson equation (42) and $\mathbb{P}[m]:=rI+m\otimes m$.

Lemma 4**.**

*The Euler-Lagrange equations for the constrained energy minimization problem (42), (44) read *

[TABLE]

Proof.

We calculate the first variation of $\mathcal{E}$ given by (44) in the direction $\phi=\phi(x)\in\mathbb{R}^{d}$. Using the expansion

[TABLE]

we have

[TABLE]

Multiplication of the Poisson equation (42) by $p_{0}$ and integration by parts, taking into account the no-flux boundary condition, gives

[TABLE]

Subtracting the identity

[TABLE]

we obtain

[TABLE]

Plugging this into (4) gives

[TABLE]

∎

Interpreting (45) as an eigenvalue-eigenvector problem for the matrix $\nabla p\otimes\nabla p$, and noting that the matrix has the eigenvalues $|\nabla p|^{2}$ (with unit eigenvector $\frac{\nabla p}{|\nabla p|}$) and zero, there exists a measurable set $\mathcal{A}\subset\Omega$ such that $m$ is a real multiple of $\nabla p$ on $\mathcal{A}$ and zero otherwise. Moreover, we have

[TABLE]

Therefore, there exist disjoint measurable sets $\mathcal{A}^{+}$, $\mathcal{A}^{-}$ such that ${\mathcal{A}^{+}}\cup{\mathcal{A}^{-}}=\mathcal{A}$ and

[TABLE]

where $\chi_{\mathcal{A}^{+}}$, resp., $\chi_{\mathcal{A}^{-}}$ are the characteristic functions of the sets $\mathcal{A}^{+}$, resp., $\mathcal{A}^{-}$. For $\gamma=1$ we define $\frac{m}{|m|}:=0$ for $m=0$. Note that the measurable sets $\mathcal{A}^{+}$, $\mathcal{A}^{-}$ can be chosen arbitrarily and each choice gives a solution of the Euler-Lagrange equations (45).

In Ref. [8], Theorem 7 and Remark 6, it was shown that if $\mathrm{meas}(\mathcal{A})<\mathrm{meas}(\Omega)$, i.e., if $m=0$ on a set of positive Lebesgue measure, then the solution $(m,p)$ is linearly unstable for the gradient flow system corresponding to (42)–(44). In other words, a necessary condition for asymptotic stability of the gradient flow system linearized at $(m,p)$ is that $\mathrm{meas}(\mathcal{A})=\mathrm{meas}(\Omega)$. We adopt this assumption for the sequel, since only the linearly stable steady states are likely to be observed in the nature or as results of numerical simulations.

Let us now choose an open subset $\Lambda\subset\Omega$ with boundary $\partial\Lambda\in C^{0,1}$ and integrate (42) over $\Lambda$, using the Green’s theorem,

[TABLE]

where $n=n(x)$ denotes the outward unit normal vector to $\partial\Lambda$, and $s=s(x)$ is the $(d-1)$-dimensional measure on $\partial\Lambda$. Here we assumed that the flux $q:=-(rI+m\otimes m)\nabla p$ is locally Lipschitz continuous on $\Lambda$, so that the generalized Green’s formula by De Giorgi-Federer[7] can be applied. We define the inflow and outflow segments of $\partial\Lambda$ as

[TABLE]

with the flux $:=-(rI+m\otimes m)\nabla p$. Note that due to (48) we have

[TABLE]

so that ${\rm sign}(n\cdot q)=-{\rm sign}(n\cdot\nabla p)$ and we can equivalently (and slightly more elegantly) define the inflow and outflow segments in terms of $\nabla p$ only,

[TABLE]

Inserting then (48) into (49) we obtain

[TABLE]

This is the Murray’s law for the linearly stable solutions of the phenomenological continuum model (42)–(44). Clearly, the two terms on the left-hand side represent the difference of in- and outgoing fluxes through the boundary $\partial\Lambda$, which is balanced by the volume integral of the external in- and outflows $S=S(x)$ on the right-hand side. Note the analogy with the Murray’s law for the discrete system (15).

Remark 1**.**

Should the restriction to only consider linearly stable states be relaxed, i.e., should it be admitted that $\mathrm{meas}(\mathcal{A})<\mathrm{meas}(\Omega)$, then (50) must be modified. Namely, since we have no information about $\nabla p$ on the complement $\mathcal{A}^{c}$ of $\mathcal{A}$ in $\Omega$, (50) needs to be rewritten as

[TABLE]

where we used the fact that $m\equiv 0$ on $\mathcal{A}^{c}$.

Let us again remark that the central assumption for the above calculations to be rigorously valid is that the flux $q=-(rI+m\otimes m)\nabla p$ is locally Lipschitz continuous on $\Lambda$.

4.1. Construction of solutions

We provide a short overview of rigorous existence results for the system (42), (45) that were obtained in Refs. [8] and [9]. Interpreting equation (45) as an eigenvalue problem for the matrix $\nabla p\otimes\nabla p$, we have the following expression for $m=m(x)$,

[TABLE]

where $\mathcal{A}_{+}\subseteq\Omega$ and $\mathcal{A}_{-}\subseteq\Omega$ are measurable disjoint sets. Inserting (51) into (42) gives

[TABLE]

subject to the no-flux boundary condition. Solutions of (52) for $\gamma>1$ are constructed by the direct method of calculus of variations.

Lemma 5**.**

Let $\gamma>1$. Then for any $S\in L^{2}(\Omega)$ and for any pair of measurable disjoint sets $\mathcal{A}_{+}$, $\mathcal{A}_{-}\subseteq\Omega$ there exists a unique (up to an additive constant) weak solution $p\in H^{1}(\Omega)\cap W^{1,2\gamma/(\gamma-1)}(\mathcal{A}_{+}\cup\mathcal{A}_{-})$ of (52).

Proof.

Direct method of calculus of variations, see Theorem 6 of Ref. [8]. ∎

In the case $\gamma=1$, (45) reads

[TABLE]

i.e., $m$ is either the zero vector or an eigenvector of the matrix $\nabla p\otimes\nabla p$ with eigenvalue $\nu>0$. The spectrum of $\nabla p\otimes\nabla p$ consists of zero and $|\nabla p|^{2}$, so that $m\neq 0$ is only possible if $|\nabla p|^{2}=\nu$. Therefore, for every stationary solution there exists a measurable function $\lambda=\lambda(x)$ such that

[TABLE]

and $p$ solves the highly nonlinear Poisson equation

[TABLE]

In Ref. [9], Section 4.2, it was shown that (53) is equivalent to the free boundary problem

[TABLE]

for some measurable function $a^{2}=a(x)^{2}$ on $\Omega$ which is the Lagrange multiplier for the condition (55). Solutions of (54)–(56) can be obtained as minimizers of the energy functional

[TABLE]

on the set $\{p\in H_{0}^{1}(\Omega),c^{2}|\nabla p|^{2}\leq 1\mbox{ a.e. on }\Omega\}$, or using the penalty method; see Section 4.2 of Ref. [9] for details.

Finally, for $1/2\leq\gamma<1$ the above approach fails due to the singularity of the term $|\nabla p(x)|^{\frac{2-\gamma}{\gamma-1}}\nabla p(x)$ at $|\nabla p|=0$ and the resulting non-boundedness from below of the associated functional. However, stationary solutions can be constructed by “cutting off” small values of $|\nabla p|$, see Section 4.3 of Ref. [9] for details.

4.2. Murray’s law for the phenomenological continuum model with diffusion

Following Refs. [8] and [9], we equip the model with a linear diffusive term that accounts for random fluctuations in the medium. In particular, we update the energy functional (44) to

[TABLE]

where $D^{2}>0$ is the diffusivity constant, $p[m]$ denotes a solution of the Poisson equation (42) and $\mathbb{P}[m]=rI+m\otimes m$. As can be easily checked, the Euler-Lagrange equations corresponding to critical points of the functional (57) constrained by the Poisson equation (42) read

[TABLE]

Note that for sufficiently regular critical points of (57), (42) we have $(-D^{2}\Delta m+\nu|m|^{2(\gamma-1)}m)\cdot m\geq 0$ almost everywhere in $\Omega$. Let us denote

[TABLE]

Then, there exist disjoint measurable sets $\mathcal{A}^{+}$, $\mathcal{A}^{-}$ such that ${\mathcal{A}^{+}}\cup{\mathcal{A}^{-}}=\mathcal{A}$ and

[TABLE]

where $\chi_{\mathcal{A}^{+}}$, resp., $\chi_{\mathcal{A}^{-}}$ are the characteristic functions of the sets $\mathcal{A}^{+}$, resp., $\mathcal{A}^{-}$. On the set $\mathcal{A}$ we have $(\nabla p\otimes\nabla p)m=0$, which is equivalent to $m\cdot\nabla p=0$.

Let us now choose an open subset $\Lambda\subset\Omega$ with boundary $\partial\Lambda\in C^{0,1}$ and define the inflow and outflow segments of $\partial\Lambda$ as

[TABLE]

where we denoted the flux $q:=-(rI+m\otimes m)\nabla p$. Using (58) in (49) leads to

[TABLE]

This is the Murray’s law for the phenomenological continuum model (44) constrained by the Poisson equation (42).

Acknowledgment

Giulia Pilli acknowleges support from the Austrian Science Fund (FWF) through the grants F 65 and W 1245.