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

Jan Haskovec ∗

Peter Markowich †

Giulia Pilli ‡

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

∗ Mathematical 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

† Mathematical 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, OskarMorgenstern-Platz 1, 1090 Vienna; peter.markowich@univie.ac.at

‡ Faculty of Mathematics, University of Vienna, Oskar-Morgenstern-Platz 1, 1090 Vienna; giulia.pilli@univie.ac.at

law,

$$E _ { p } = Q \Delta P,$$

where Q is the volumetric flux and ∆ 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,

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

where R &gt; 0 is the radius of the tube, L &gt; 0 its length, and µ &gt; 0 the dynamic viscosity of the fluid. Therefore, we have

$$E _ { p } = \frac { 8 \mu L } { \pi R ^ { 4 } } Q ^ { 2 }. \\. \quad.$$

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,

$$E _ { m } = \nu \pi L R ^ { 2 },$$

where the proportionality constant ν &gt; 0 is called the metabolic coefficient. This gives for the total energy

$$E \coloneqq E _ { p } + E _ { m } = \frac { 8 \mu L } { \pi R ^ { 4 } } Q ^ { 2 } + \nu \pi L R ^ { 2 }. \quad \quad \quad ( 2 )$$

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

$$Q = \sqrt { \frac { \pi ^ { 2 } \nu } { 1 6 \$$

In a branching node with a parent branch with flow rate Q &gt; 0 0 and children branches with flow rates Q ,Q ,... , Q 1 2 n &gt; 0, we have the flux balance

$$\psi \psi ^ { \prime } \psi ^ { \prime } \psi$$

Using (3), we obtain for the respective radii R &gt; 0 0 and R , R , . . . , R 1 2 n &gt; 0,

$$\begin{array} {$$

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 &gt; i 0, which for threedimensional laminar flow are proportional to R 4 i , see Ref. [17], therefore

$$C _ { 0 } ^ { 3 / 4 } & = \sum _ { i = 1 } ^$$

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 ∈ V and j ∈ V by ( i, j ) ∈ E . Since the graph is undirected, ( i, j ) and ( j, i ) refer to the same edge. For each edge ( i, j ) ∈ E of the graph G we consider its length and its conductivity, denoted by L ij = L ji &gt; 0 and C ij = C ji ≥ 0, respectively. The edge lengths L ij &gt; 0 are given as a datum and fixed for all ( i, j ) ∈ E . With each vertex i ∈ V there is associated the fluid pressure P i ∈ R . The pressure drop between vertices i ∈ V and j ∈ V connected by an edge ( i, j ) ∈ E is denoted by

Note that the pressure drop is antisymmetric, i.e., by definition, (∆ P ) ij = -(∆ P ) ji . The oriented flux (flow rate) from vertex j ∈ V to i ∈ 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 ) ∈ E is proportional to the fourth power of its diameter[17]. Then, (1) is reformulated in terms of the conductivities as

$$( \Delta P ) _ { i j } \coloneqq P _ { j } - P _$$

$$Q _ { i j } \coloneqq C _ { i j } \frac { P _ { j } - P _ { i } } { L _ { i j } } \quad \text{ for all } ( i, j ) \in E.$$

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

$$\mu \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \colon \\ - \sum _ { j \in \mathcal { N } ( i ) } C _ { i j } \frac { P _ { j } - P _ { i } } { L _ { i j } } = S _ { i } \quad \text{ for all } i \in V. \quad \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{                                                                                                                                                                                                         }                                                                                                                                                                                                       
\int ( i \, \searrow \, V \, d o n o t o s \, \text{ the set of vertices connected to } i \, \in V \, \text{ through an } }$$

Here N ( ) i ⊂ V denotes the set of vertices connected to i ∈ V through an edge, i.e.,

$$\mathcal { N } ( i ) & \coloneqq \{ j \in V ; \, ( i, j ) \in E \}, \\ \colon = \colon = \max \max \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \max \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \max \max \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \max \min \max \min \min \min \min \min \min \min \min \min \min \min \max \min \min \min \min \min \min \min \max \max \min \min \min \min \min \min \min \max \max \max \min \min \min \min \min \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \sum \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \min \sum \min \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \min \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \text{ }$$

and S = ( S i ) i ∈ V is the prescribed strength of the external flow source ( S &gt; i 0) or sink ( S i &lt; 0) in node i ∈ V . Clearly, a necessary condition for the solvability of (8) is the global mass conservation

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

which we shall assume in the sequel. Given the vector of conductivities C = ( C ij ) ( i,j ) ∈ E , the Kirchhoff law (8) is a linear system of equations for the vector of pressures P = ( P i ) i ∈ 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 ) ∈ E is connected[2], where only edges with positive conductivities C ij &gt; 0 are taken into account (i.e., edges with zero conductivities are discarded). Note that the solution vector P = ( P i ) i ∈ 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,

$$\mathcal { E } [ C ] \coloneqq \sum _ { ( i, j ) \in E } ^ { \lambda } \left ( \frac { Q _ { i j } [ C ] ^ { 2 } } { C _ { i j } } + \frac { \nu } { \gamma } C _ { i j } ^ { \gamma } \right ) L _ { i j }. \quad \ ( 1 0 ) \\ \text{consider a generalization where the matabolic energy is non-functional}$$

Here we consider a generalization where the metabolic energy is proportional to a power γ &gt; 0 of the respective edge conductivity. We have γ = 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 1 2 / ij . Consequently, the effective value of γ typically used in models of leaf venation in plants lies between 1 / 2 and 1, see, e.g., Refs. [10, 11]. Again, ν &gt; 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 ) ∈ 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 ∈ V is a solution of the linear system (8) with the vector of conductivities C . Then, for any fixed ( k, l ) ∈ E we have

$$\begin{array} { c c c } \cdot \cdot & \cdot \cdot \\ & \frac { \partial } { \partial C _ { k l } } \sum _ { ( i, j ) \in E } \frac { Q _ { i j } [ C ] ^ { 2 } } { C _ { i j } } L _ { i j } = - \frac { Q _ { k l } [ C ] ^ { 2 } } { C _ { k l } ^ { 2 } } L _ { k l }. \end{array}$$

Proof. Since

$$\begin{array} { c } \frac { \ell } { \partial C _ { k l } } \sum _ { ( i, j ) \in E } \frac { Q _ { i j } [ C ] ^ { 2 } } { C _ { i j } } L _ { i j } = - \frac { Q _ { k l } [ C ] ^ { 2 } } { C _ { k l } ^ { 2 } } L _ { k l } + 2 \sum _ { ( i, j ) \in E } \frac { Q _ { i j } [ C ] } { C _ { i j } } \frac { \partial Q _ { i j } [ C ] } { \partial C _ { k l } } L _ { i j }, \\ \text{if is sufficient to show that} \end{array}$$

it is sufficient to show that

$$\sum _ { ( i, j ) \in E } \frac { Q _ { i j } [ C ] } { C _ { i j } } \frac { \partial Q _ { i j } [ C ] } { \partial C _ { k l } } L _ { i j } = 0.$$

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

$$\mathbb { A } _ { i j } = \begin{cases} 0 & \text{if $(i,j)\not\in E$}, \\ 1 & \text{if $(i,j)\in E$}. \end{cases} \\ \Big. \, \dots \, \dots \, \dots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \end{cases}$$

Note that G is an undirected graph, implying A ij = A ji . Due to the symmetry of C ij and L ij and antisymmetry of Q ij we have

$$\text{Note that $c$ is an unrefected graph, implying $A_{ij} = A_{ji}$. Due to the sym-} \\ \text{entry of $C_{ij}$ and $L_{ij}$ and antisymmetry of $Q_{ij}$ we have} \\ 2 \sum _ { ( i, j ) \in E } \frac { Q _ { i j } [ C ] } { C _ { i j } } \frac { \partial Q _ { i j } [ C ] } { \partial C _ { k l } } L _ { i j } & = \sum _ { i = 1 } ^ { n } \sum _ { j = 1 } ^ { n } \mathbb { A } _ { i j } \left ( \frac { P _ { j } - P _ { i } } { L _ { i j } } \frac { \partial Q _ { i j } [ C ] } { \partial C _ { k l } } \right ) L _ { i j } \\ & = \sum _ { j = 1 } ^ { n } P _ { j } \sum _ { i = 1 } ^ { n } \mathbb { A } _ { i j } \frac { \partial Q _ { i j } [ C ] } { \partial C _ { k l } } - \sum _ { i = 1 } ^ { n } P _ { i } \sum _ { j = 1 } ^ { n } \mathbb { A } _ { i j } \frac { \partial Q _ { i j } [ C ] } { \partial C _ { k l } } \\ & = - 2 \sum _ { i = 1 } ^ { n } P _ { i } \sum _ { j = 1 } ^ { n } \mathbb { A } _ { i j } \frac { \partial Q _ { i j } [ C ] } { \partial C _ { k l } } \\ & = - 2 \sum _ { i = 1 } ^ { n } P _ { i } \frac { \partial } { \partial C _ { k l } } \sum _ { j \in N ( i ) } Q _ { i j }. \\ \text{By the definition of the flow rate $Q_{ij}$ in $(7)$ and Kirchhoff's law $(8)$ we have} \\ - \quad - \quad P. - P.$$

By the definition of the flow rate Q ij in (7) and Kirchhoff's law (8) we have

$$- \sum _ { j \in N ( i ) } Q _ { i j } = \sum _ { j \in N ( i ) } C _ { i j } \frac { P _ { j } - P _ { i } } { L _ { i j } } = S _ { i }, \\ \dots \quad & \quad \dots \quad \alpha \quad \alpha \quad \dots \quad \dots$$

and since the sources/sinks S i are fixed, we conclude

$$\sum _ { ( i, j ) \in E } \frac { Q _ { i j } [ C ] } { C _ { i j } } \frac { \partial Q _ { i j } [ C ] } { \partial C _ { k l } } L _ { i j } = 0.$$

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,

$$\frac { \partial } { \partial C _ { i j } } \mathcal { E } [ C ] = - \left ( \frac { Q _ { i j } [ C ] ^ { 2 } } { C _ { i j } ^ { 2 } } - \nu C _ { i j } ^ { \gamma - 1 } \right ) L _ { i j }.$$

Since any critical point of E [ C ] satisfies

$$\frac { \partial } { \partial C _ { i j } } \mathcal { E } [ C ] = 0 \quad \text{for all } ( i, j ) \in E,$$

we have

$$Q _ { i j } [ C ] ^ { 2 } = \nu C _ { i j } ^ { \gamma + 1 } \quad \text{for all } ( i, j ) \in E. \quad \text{                                                                                                                                                                                                         (13) \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \ \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \end{cases}$$

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

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

Let us write the set of neighbors N ( ) of the node i i ∈ V as the disjoint union N ( ) i = N + ( ) i ∪ N -( ), i based on the flow directions along the respective edges, i.e.,

$$\mathcal { N } ^ { + } ( i ) \coloneqq \{ j \in \mathcal { N } ( i ) ; \, Q _ { i j } \geq 0 \}, \quad \mathcal { N } ^ { - } ( i ) \coloneqq \{ j \in \mathcal { N } ( i ) ; \, Q _ { i j } < 0 \}.$$

Then, we can rewrite (14) as

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

and, further,

$$\sum _ { j \in \mathcal { N } ^ { + } ( i ) } | Q _ { i j } | + S _ { i } = \sum _ { j \in \mathcal { N } ^ { - } ( i ) } | Q _ { i j } | \quad \text{for all } i \in V.$$

Using (13), we arrive at

$$\psi \cdot \cdot \cdot \cdot \cdot \\ \sqrt { \nu } \sum _ { j \in \mathcal { N } ^ { + } ( i ) } C _ { i j } ^ { \frac { \gamma + 1 } { 2 } } + S _ { i } = \sqrt { \nu } \sum _ { j \in \mathcal { N } ^ { - } ( i ) } C _ { i j } ^ { \frac { \gamma + 1 } { 2 } } \quad \text{ for all } i \in V \text{.} \quad \text{ } ( 1 5 ) \\ \pi \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \end{matrix}$$

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 γ = 1 2, / which leads to the power γ +1 2 = 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 γ &gt; 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

$$\mathcal { E } [ c ] & \coloneqq \int _ { \Omega } \nabla p \cdot c \nabla p + \frac { \nu } { \gamma } | c | ^ { \gamma } \, \mathrm d x, & ( 1 6 ) \\ \text{-aholo coefficient } \nu \succ \mathfrak { m } \text{ and } \mathfrak { t h o v n o n o n t } \sim \mathfrak { m } \text{ } 1 \text{ } \mathfrak { t h o n o r r }$$

$$c = \begin{pmatrix} c ^ { 1 } \\ \ddots \\ c ^ { d } \end{pmatrix}. \quad & ( 1 7 ) \\ \text{is defined as } \sum _ { k = 1 } ^ { d } | c ^ { k } | ^ { \gamma }. \text{ The scalar pressure } p = p ( x ) \text{ of}$$

with the metabolic coefficient ν &gt; 0 and the exponent γ &gt; 1. The energy functional E [ c ] is defined on the set of nonnegative diagonal tensor fields c = ( ) on a bounded domain Ω c x ⊂ R d , d ∈ N ,

The symbol | c | γ is defined as ∑ d k =1 ∣ ∣ c k ∣ ∣ γ . The scalar pressure p = p x ( ) of the fluid occupying the domain Ω is subject to the Poisson equation

$$- \nabla \cdot ( c \nabla p ) = S, \quad \ \ \ \ \ ( 1 8 ) \\. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad.$$

equipped with no-flux boundary condition on ∂ Ω, and the datum S = S x ( ) represents the intensity of sources and sinks in Ω. 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

$$- \nabla \cdot ( \mathbb { P } [ c ] \nabla p ) = S, \quad \quad \quad ( 1 9 )$$

with the permeability tensor

$$\mathbb { P } [ c ] \coloneqq r I + c,$$

where r = r x ( ) ≥ r 0 &gt; 0 is a prescribed function that models the isotropic background permeability of the medium, and I ∈ R d × d is the unit matrix. Again, (19) is equipped with the no-flux boundary condition

$$n \cdot \mathbb { P } [ c ] \nabla p = 0 \quad \text{ on } \partial \Omega, \quad \text{ (2)} \\. \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \colon \end{cases}$$

$$\cdot \cdot \cdot \\ \int _ { \Omega } \nabla \Phi \cdot ( r I + c ) \nabla p \, \mathrm d x = \int _ { \Omega } S \phi \, \mathrm d x. \\ \text{since we only admit nonnegative diagonal tensor fields $c=c(x)$}.$$

where n = n x ( ) denotes the outer unit normal vector to ∂ Ω. Then, the weak formulation of the boundary value problem (19)-(21) with a test function φ ∈ C ∞ (Ω) reads

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

$$\mathcal { E } [ c ] & \coloneqq \int _ { \Omega } \nabla p \cdot \mathbb { P } [ c ] \nabla p + \frac { \nu } { \gamma } | c | ^ { \gamma } \, \mathrm d x. \text{ \quad \ } ( 2 3 ) \\ \text{calculate the Fouler-Lagranger equations corresponding to critical}$$

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

$$( \partial _ { x _ { k } } p ) ^ { 2 } - \nu | c ^ { k } | ^ { \gamma - 1 } = 0, \quad \text{ for } k = 1, \dots, d. \quad \text{ } \quad \text{ } \quad ( 2 4 ) \\ \tilde { \ }. \ \stackrel { \bullet } { \bullet }, \quad \tilde { \ }. \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdots$$

Proof. We calculate the first variation of E in the direction φ where φ denotes a diagonal matrix with entries φ , . . . , φ 1 d , such that c + εφ is nonnegative. Using the expansion we have

$$p [ c + \varepsilon \phi ] = p _ { 0 } + \varepsilon p _ { 1 } + \mathcal { O } ( \varepsilon ^ { 2 } ),$$

$$& \text{we have} \\ & \frac { \text{d} } { \text{d} \varepsilon } \mathcal { E } [ c + \varepsilon \phi ] \Big | _ { \varepsilon = 0 } = \sum _ { k = 1 } ^ { d } \int _ { \Omega } ( \partial _ { x _ { k } } p _ { 0 } ) ^ { 2 } \, \phi ^ { k } + 2 ( r + c ^ { k } ) ( \partial _ { x _ { k } } p _ { 0 } ) ( \partial _ { x _ { k } } p _ { 1 } ) + \nu | c ^ { k } | ^ { \gamma - 2 } c ^ { k } \phi ^ { k } \, \text{d} x. \\ & \text{Using } p _ { 0 } \text{ as a test function in (22) with the permeability tensor } \mathbb { P } [ c + \varepsilon \phi ] =$$

Using p 0 as a test function in (22) with the permeability tensor P [ c + εφ ] = rI + c + εφ gives

$$& \sum _ { k = 1 } ^ { d } \int _ { \Omega } \left ( r + c ^ { k } + \varepsilon \phi ^ { k } \right ) ( \partial _ { x _ { k } } p _ { 0 } ) ^ { 2 } + \varepsilon ( r + c ^ { k } ) ( \partial _ { x _ { k } } p _ { 0 } ) ( \partial _ { x _ { k } } p _ { 1 } ) \, d x = \int _ { \Omega } S p _ { 0 } \, d x + \mathcal { O } ( \varepsilon ^ { 2 } ). \\ & \text{Subtracting the identity}$$

Subtracting the identity

$$\Phi \sum _ { k = 1 } ^ { d } \int _ { \Omega } ( r + c ^ { k } ) \, ( \partial _ { x _ { k } } p _ { 0 } ) ^ { 2 } \, \mathrm d x = \int _ { \Omega } S p _ { 0 } \, \mathrm d x,$$

we obtain

$$\sum _ { k = 1 } ^ { d } & \int _ { \Omega } ( \partial _ { x _ { k } } p _ { 0 } ) ^ { 2 } \, \phi ^ { k } + ( r + c ^ { k } ) ( \partial _ { x _ { k } } p _ { 0 } ) ( \partial _ { x _ { k } } p _ { 1 } ) \, \mathrm d x = 0. \\ \text{gging this into } ( 2 6 ) \text{ gives}$$

Plugging this into (26) gives and (24) follows.

$$\dashrightarrow \circ \dashrightarrow \circ \dashrightarrow \dashrightarrow \cup \dashrightarrow \\ \frac { d } { d \varepsilon } \mathcal { E } [ c + \varepsilon \phi ] \Big | _ { \varepsilon = 0 } & = \sum _ { k = 1 } ^ { d } \int _ { \Omega } \left [ - \left ( \partial _ { x _ { k } } p _ { 0 } \right ) ^ { 2 } + \nu | c ^ { k } | ^ { \gamma - 1 } \right ] \phi ^ { k } \, \mathrm d x, \\ \text{and} \left ( 24 \right ) follows. \quad & \square \\ \text{From } \mathcal { M } \, \text{it follow that for } k - 1 \quad \text{$d$} \text{or visit distinct to} \, \text{absome}.$$

/square

From (24) it follows that for k = 1 , . . . , d there exist disjoint Lebesguemeasurable sets A + k , A -k such that Ω = A + k ∪ A -k for all k = 1 , . . . , d and

$$\begin{array} {$$

Note that by assumption c k ≥ 0. Let us now choose an open subset Λ ⊂ Ω with boundary ∂ Λ ∈ C 0 1 , and use the characteristic function of Λ as the test function φ 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 ∇ p is locally Lipschitz continuous on Λ. Then the generalized Green's formula by De Giorgi-Federer[7] gives

$$& \cdot \cdot \\ & \quad - \int _ { \partial \Lambda$$

$$- \sqrt { \nu } \int _ { \partial \Lambda } \sum _$$

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

We define the inflow and outflow segments of ∂ Λ as

∂ Λ + := { x ∈ ∂ Λ; n x ( ) · q x ( ) &lt; 0 } , ∂ Λ -:= { x ∈ ∂ Λ; n x ( ) · q x ( ) &gt; 0 } , with the flux q = -( rI + ) c ∇ p . Then, (30) can be written in the form

∣ ∣ This is the Murray's law for the continuum model (19)-(23).The two terms on the left-hand side represent the difference of inand outgoing fluxes through the boundary ∂ Λ, 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).

$$\text{with the flux } q = - ( r I + c ) \nabla p.$$

Let us remark that the above formal calculation becomes rigorous if we assume sufficient regularity of the flux q = -( rI + ) c ∇ p , in particular, we need q to be locally Lipschitz continuous on Λ. 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

$$c ^ { k } = \nu ^ { - \frac { 1 } { \gamma - 1$$

Inserting this into (19), we obtain

$$- \sum _ { k = 1 } ^ { d } \partial _ { x _ { k } } \left ( \left ( r + \nu ^ { - \frac { 1 } { \gamma - 1 } } ( \partial _ { x _ { k } } p ) ^ { \frac { 2 } { \gamma - 1 } } \right ) \partial _ { x _ { k } } p \right ) = S. \\ \text{ne weak formulation of this nonlinear elliptic equation. examined with a}$$

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

$$\int _ { \Omega } & \sum _ { k = 1 } ^ { d } \left ( r + \nu ^ { - \frac { 1 } { \gamma - 1 } } ( \partial _ { x _ { k } } p ) ^ { \frac { 2 } { \gamma - 1 } } \right ) ( \partial _ { x _ { k } } p ) ( \partial _ { x _ { k } } \phi ) \, \mathrm d x = \int _ { \Omega } S \phi \, \mathrm d x \quad \quad ( 3 3 ) \\ \text{for all test functions $\phi \in \mathcal{C}^{\infty}(\Omega)$   For $\gamma > 1$ softions can be constructed } \text{hv}$$

for all test functions φ ∈ C ∞ (Ω). For γ &gt; 1 solutions can be constructed by the direct method of calculus of variations.

Lemma 3. Let γ &gt; 1 . For every S ∈ L 2 (Ω) there exists a unique solution p ∈ H 1 (Ω) of (33) satisfying ∫ Ω p x ( ) d x = 0 .

Proof. We define the functional F : H 1 (Ω) → R by

$$\mathcal { F } [ p ] \coloneqq & \int _ { \Omega } \left ( \sum _ { k = 1 } ^ { d } \frac { r } { 2 } | \partial _ { x _ { k } } p | ^ { 2 } + \nu ^ { - \frac { 1 } { \gamma - 1 } } | \partial _ { x _ { k } } p | ^ { \frac { 2 \gamma } { \gamma - 1 } } \right ) \, \mathrm d x - \int _ { \Omega } p S \, \mathrm d x. \quad \ \ ( 3 4 ) \\ \Im \infty \ \ 2 \gamma \ \ > \ \mathbb { m } { f o w } \sim 1 \ \ \text{the non-normal} \ \mathcal { T } \ \text{is uniformly common} \ \mathcal { m } { v } \ \in \mathcal { M } { 1 } ( \Omega )$$

Since 2 γ γ -1 &gt; 2 for γ &gt; 1, the functional F is uniformly convex on H 1 (Ω). Moreover, a straightforward application of the Poincar´ inequality provides e coercivity on the set H 1 0 , (Ω) := { p ∈ H 1 (Ω); ∫ Ω p x ( ) d x = 0 } . The classical theory (see, e.g., Ref. [6]) provides the existence of a unique minimizer p ∈ H 1 0 , (Ω). Since, as is easy to check, (34) is the Euler-Lagrange equation for the minimizer of F , p is the unique solution in H 1 0 , (Ω) of (33). /square

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

$$\mathcal { A } _ { k } ^ { + } = \{ x \in \Omega ; \partial _ { x _ { k } } p ( x ) > 0 \}, \quad \mathcal { A } _ { k } ^ { - } = \{ x \in \Omega ; \partial _ { x _ { k } } p ( x ) < 0 \}.$$

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

$$\mathcal { E } [ c ] \coloneqq & \int _ { \Omega } D ^ { 2 } | \nabla c | ^ { 2 } + \nabla p \cdot \mathbb { P } [ c ] \nabla p + \frac { \nu } { \gamma } | c | ^ { \gamma } \, \mathrm d x, \\ \intertext { a n d } D ^ { 2 } & > \, \partial \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \end{array}$$

where D &gt; 2 0 is the diffusivity constant and the symbol |∇ | c 2 is defined as ∑ d k =1 |∇ c k | 2 . The permeability tensor is given by (20), i.e., 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

$$( \partial _ { x _ { k } } p ) ^ { 2 } + D ^ { 2 } \Delta c ^ { k } - \nu | c ^ { k } | ^ { \gamma - 1 } = 0, \quad \text{ for } k = 1, \dots, d. \quad \text{ } \quad ( 3 6 ) \\ \intertext { s........................................................................................................................................................................................................ }$$

Note that for sufficiently regular critical points of (35), (22) we have ν c | k γ | -1 -D 2 ∆ c k ≥ 0 almost everywhere on Ω. Then it follows that for k = 1 , . . . , d there exist disjoint Lebesgue-measurable sets A + k , A -k such that Ω = A + k ∪ A -k for all k = 1 , . . . , d and

$$\partial _ { x _ { k } } p & = ( \chi _ { \mathcal { A } _ { k } ^ { + } } - \chi _ { \mathcal { A } _ { k } ^ { - } } ) \left ( \nu | c ^ { k } | ^ { \gamma - 1 } - D ^ { 2 } \Delta c ^ { k } \right ) ^ { 1 / 2 }, \quad \text{ for } k = 1, \dots, d. \quad ( 3 7 ) \\ & \text{$As before we choose an open subset $\Lambda \subset \Omega$ with boundary $\partial \Lambda \in C^{0,1}$, and}$$

As before we choose an open subset Λ ⊂ Ω with boundary ∂ Λ ∈ C 0 1 , , and using (37) in (29), we obtain

$$& - \int _ { \partial \Lambda } \sum _ { k = 1 } ^ { d } ( \chi _ { \mathcal { A } _ { k } ^ { + } } - \chi _ { \mathcal { A } _ { k } ^ { - } } ) ( r + c ^ { k } ) \left ( \nu | c ^ { k } | ^ { \gamma - 1 } - D ^ { 2 } \Delta c ^ { k } \right ) ^ { 1 / 2 } n _ { k } \, d s = \int _ { \Lambda } S \, \mathrm d x, ( 3 8 ) \\ & \text{where $n=n(x)$ denotes the outward unit normal vector to $\partial\Lambda$ and $s=s(x)$}$$

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

$$\partial \Lambda ^ { + } & \coloneqq \{ x \in \partial \Lambda ; \, n ( x ) \cdot q ( x ) < 0 \}, \quad \partial \Lambda ^ { - } \coloneqq \{ x \in \partial \Lambda ; \, n ( x ) \cdot q ( x ) > 0 \}, \\ \pi \pi & \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \nu \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \Pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \right \}.$$

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

$$& \text{the form} \\ & \quad - \int _ { \partial \Lambda ^ { + } } \left | \sum _ { k = 1 } ^ { d } ( \chi _ { \mathcal { A } _ { k } ^ { + } } - \chi _ { \mathcal { A } _ { k } ^ { - } } ) ( r + c ^ { k } ) \left ( \nu | c ^ { k } | ^ { \gamma - 1 } - D ^ { 2 } \Delta c ^ { k } \right ) ^ { 1 / 2 } n _ { k } \right | \, d s \\ & \quad + \int _ { \partial \Lambda ^ { - } } \left | \sum _ { k = 1 } ^ { d } ( \chi _ { \mathcal { A } _ { k } ^ { + } } - \chi _ { \mathcal { A } _ { k } ^ { - } } ) ( r + c ^ { k } ) \left ( \nu | c ^ { k } | ^ { \gamma - 1 } - D ^ { 2 } \Delta c ^ { k } \right ) ^ { 1 / 2 } n _ { k } \right | \, d s \\ & \quad = \int _ { \Lambda } S \, \text{d} x. \\ & \text{This is the Murray's law for the continuum model (35) constrained by the} \\ & \text{Poisson equation (22)}$$

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 Ω ⊂ R d is occupied by a porous medium with the permeability tensor

$$\mathbb { P } [ m ] \coloneqq m & \otimes m \\ \math$$

with m = m x ( ) ∈ R d a smooth vector field on R d . Note that P [ m ] has the eigenvalue | m | 2 with eigenvector m , and 0 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 P , i.e., C i /similarequal | m | 2 in the neighborhood of a node i ∈ V . Assuming

quasi-incompressibility of the fluid, i.e., constant fluid density along particle trajectories, we have the local mass conservation law

$$\nabla \cdot q = S \quad & \text{ in } \Omega,$$

where q = q x ( ) ∈ 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

$$q = - \mathbb { P } [ m ] \nabla p, \quad \begin$$

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

$$\begin{array} {$$

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,

$$\begin{array} {$$

where I is the identity matrix and the scalar function r = r x ( ) ≥ r 0 &gt; 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 ∂ Ω, i.e., equip (42) with the homogeneous Neumann boundary condition

$$\begin{array} {$$

where n = n x ( ) denotes the outer unit normal vector to ∂ Ω. Then the solution p = p x ( ) is unique up to an additive constant and the flux q = -P [ m ] ∇ 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

$$\mathcal { E } [ m ] = \int \nabla p [ m ] \cd$$

where p m [ ] denotes a solution of the Poisson equation (42) and P [ m ] := rI + m ⊗ m .

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

$$( \nabla p \otimes \nabla p ) m - \nu | m |$$

Proof. We calculate the first variation of E given by (44) in the direction φ = φ x ( ) ∈ R d . Using the expansion

$$p [ m + \varepsilon \phi ] = p _ { 0 } + \vare$$

we have

$$\frac { \text{d} } { \text{d} \varepsilon } \$$

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

$$\int _ { \Omega } \stackrel { \lambda } { | \nabla p _ { 0 } | ^ { 2 } + | m \cdot \nabla p _ { 0 } | ^ { 2 } + \varepsilon r \nabla p _ { 0 } \cdot \nabla p _ { 1 } + \varepsilon ( \nabla p _ { 0 } \cdot m ) ( \nabla p _ { 1 } \cdot m ) + 2 \varepsilon ( \nabla p _ { 0 } \cdot m ) ( \nabla p _ { 0 } \cdot \phi ) \, d x } \int _ { \Omega } \stackrel { \lambda } { | \nu } \quad \cdots \quad \varepsilon \quad \cdots$$

Subtracting the identity

$$= \int _ { \Omega } S p _ { 0 } \, \mathrm d x + \mathcal { O } ( \varepsilon ^ { 2 } ).$$

we obtain

$$\mu \quad \int _ { \Omega } r | \vec { \nabla p _ { 0 } } | ^ { 2 } + | m \cdot \nabla p _ { 0 } | ^ { 2 } \, \mathrm d x = \int _ { \Omega } S p _ { 0 } \, \mathrm d x,$$

$$\int _ { \Omega } r \nabla p _ { 0 } \cdot \nabla p _ { 1 } + ( \nabla p _ { 0 } \cdot m ) ( \nabla p _ { 1 } \cdot m ) \, \mathrm d x = - 2 \int _ { \Omega } ( \nabla p _ { 0 } \cdot m ) ( \nabla p _ { 0 } \cdot \phi ) \, \mathrm d x. \\ \text{Pln} \sigma \sin \sigma \text{ this into } ( 4 7 ) \text{ oives}$$

Plugging this into (47) gives

$$\cdots \circ \cdots \circ \cdots \circ \cdots \circ \cdots \circ \cdots \\ \frac { d } { d \varepsilon } \mathcal { E } [ m + \varepsilon \phi ] \Big | _ { \varepsilon = 0 } & = 2 \int _ { \Omega } \Big [ - ( \nabla p _ { 0 } \otimes \nabla p _ { 0 } ) m + \nu | m | ^ { 2 ( \gamma - 1 ) } m \Big ] \cdot \phi \, d x. \\ \square \quad \text{Interpreting } ( 4 5 ) \text{ as an eigenvalue-eigenvector problem for the matrix}$$

$$\nu | m | ^ { 2 ( \gamma - 1 ) } = | \nabla p | ^ { 2 } \quad \text{ on } \mathcal { A }. \\ \dots \quad \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \colon \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \ddots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots$$

Interpreting (45) as an eigenvalue-eigenvector problem for the matrix ∇ ⊗∇ p p , and noting that the matrix has the eigenvalues |∇ | p 2 (with unit eigenvector ∇ p |∇ | p ) and zero, there exists a measurable set A ⊂ Ω such that m is a real multiple of ∇ p on A and zero otherwise. Moreover, we have

Therefore, there exist disjoint measurable sets A + , A -such that A ∪A + -= A and where χ A + , resp., χ A -are the characteristic functions of the sets A + , resp., A -. For γ = 1 we define m | m | := 0 for m = 0. Note that the measurable sets A + , A -can be chosen arbitrarily and each choice gives a solution of the Euler-Lagrange equations (45).

$$\nabla p & = ( \chi _ { \mathcal { A } ^ { + } } - \chi _ { \mathcal { A } ^ { - } } ) \sqrt { \nu } | m | ^ { \gamma - 2 } m \quad \text{ on } \mathcal { A }, \quad \text{ (4.8)} \\ \mathfrak { r o w } \cdots \quad \text{ r o w to the constructor function of the code } \mathcal { A } ^ { + } \quad \text{ r o w }$$

In Ref. [8], Theorem 7 and Remark 6, it was shown that if meas( A ) &lt; meas(Ω), 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 meas( A ) = meas(Ω). 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 Λ ⊂ Ω with boundary ∂ Λ ∈ C 0 1 , and integrate (42) over Λ, using the Green's theorem,

$$\cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cd. & - \int _ { \partial \Lambda } n \cdot ( r I + m \otimes m ) \nabla p d s = \int _ { \Lambda } S \, \mathrm d x, & ( 4 9 ) \\ e \, n = n ( r ) \, \text{ denotes the outward unit normal vector to } \partial \Lambda \, \text{ and } s = s ( r )$$

where n = n x ( ) denotes the outward unit normal vector to ∂ Λ, and s = ( ) s x is the ( d -1)-dimensional measure on ∂ Λ. Here we assumed that the flux q := -( rI + m ⊗ m ) ∇ p is locally Lipschitz continuous on Λ, so that the

generalized Green's formula by De Giorgi-Federer[7] can be applied. We define the inflow and outflow segments of ∂ Λ as

$$& \partial \Lambda ^ { + } \coloneqq \{ x \in \partial \Lambda ; \, n ( x ) \cdot q ( x ) < 0 \}, \quad \partial \Lambda ^ { - } \coloneqq \{ x \in \partial \Lambda ; \, n ( x ) \cdot q ( x ) > 0 \}, \\ & \quad \text{with} \, \quad \text{$k$} \, \text{$two$} \, \text{$a$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$ } \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \quad \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \colon \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \tilde { $two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$.} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$} \, \text{$two$}$$

$$q = - ( r I + m \otimes m ) \nabla p = - ( r I + | m | ^ { 2 } ) \nabla p \quad \text{ on } \mathcal { A }, \\ \intertext { s t h e t i m } \lambda \lambda = - \intertext { s t h e t i m } \nabla \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \xi } }$$

with the flux := -( rI + m ⊗ m ) ∇ p . Note that due to (48) we have so that sign( n · q ) = -sign( n · ∇ p ) and we can equivalently (and slightly more elegantly) define the inflow and outflow segments in terms of ∇ p only,

Inserting then (48) into (49) we obtain

$$& \partial \Lambda ^ { + } \coloneqq \{ x \in \partial \Lambda ; \, n ( x ) \cdot \nabla p ( x ) > 0 \}, \quad \partial \Lambda ^ { - } \coloneqq \{ x \in \partial \Lambda ; \, n ( x ) \cdot \nabla p ( x ) < 0 \}. \\ & \tau _ { \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \ dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \$$

$$\text{inserting then } ( 4 8 ) \, & \text{into } ( 4 9 ) \, \text{ we obtain} \\ - \sqrt { \nu } \int _ { \partial \Lambda ^ { + } } \left ( r + | m | ^ { 2 } \right ) | m | ^ { \gamma - 2 } | m \cdot n | \, \text{d} s \quad & ( 5 0 ) \\ + \sqrt { \nu } \int _ { \partial \Lambda ^ { - } } \left ( r + | m | ^ { 2 } \right ) | m | ^ { \gamma - 2 } | m \cdot n | \, \text{d} s \ \ = \ \int _ { \Lambda } S \, \text{d} x. \\ \text{This is the Murray's law for the linearly stable solutions of the phenomena} \\ \text{ogical continuum model } ( 4 2 ) - ( 4 4 ). \text{ Clearlv. the two terms on the left-hand}$$

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 ∂ Λ, 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 meas( A ) &lt; meas(Ω) , then (50) must be modified. Namely, since we have no information about ∇ p on the complement A c of A in Ω (50) , needs to be rewritten as where we used the fact that m ≡ 0 on A c .

$$& \text{complement } \mathcal { A } \langle \text{ of } \mathcal { A } \ u n \, \underline { L }, \, ( \mathfrak { v } ) \text{ needs to be rewrttten as} \\ & - \sqrt { \nu } \int _ { \mathcal { A } \cap \partial \Lambda ^ { + } } \left ( r + | m | ^ { 2 } \right ) | m | ^ { \gamma - 2 } | m \cdot n | \, \mathrm d s - \int _ { \mathcal { A } ^ { c } \cap \partial \Lambda ^ { + } } r | \nabla p \cdot n | \, \mathrm d s \\ & + \sqrt { \nu } \int _ { \mathcal { A } \cap \partial \Lambda ^ { - } } \left ( r + | m | ^ { 2 } \right ) | m | ^ { \gamma - 2 } | m \cdot n | \, \mathrm d s + \int _ { \mathcal { A } ^ { c } \cap \partial \Lambda ^ { + } } r | \nabla p \cdot n | \, \mathrm d s \ \ = \ \int _ { \Lambda } S \, \mathrm d x, \\ & \text{where we used the fact that } m \equiv 0 \text{ on } \mathcal { A } ^ { c }. \\ & \quad \quad \quad$$

Let us again remark that the central assumption for the above calculations to be rigorously valid is that the flux q = -( rI + m m ⊗ ) ∇ p is locally Lipschitz continuous on Λ.

- 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 ∇ ⊗∇ p p , we have the following expression for m = m x ( ),

m x ( ) := ( χ A + ( x ) -χ A -( x ) ) ν 1 2(1 -γ ) |∇ p x ( ) | 2 -γ γ -1 ∇ p x , ( ) (51) where A + ⊆ Ω and A -⊆ Ω are measurable disjoint sets. Inserting (51) into (42) gives

$$\cdot \cdot \cdot \cdot \begin{array} { c } - \nabla \cdot \left [ \left ( r I + \nu ^ { \frac { 1 } { 1 - \gamma } } | \nabla p ( x ) | ^ { \frac { 2 } { \gamma - 1 } } \chi _ { \mathcal { A } _ { + } \cup \mathcal { A } _ { - } } ( x ) \right ) \nabla p ( x ) \right ] = S, \\ \text{subject to the no-flux boundary condition. Solutions of } ( 5 2 ) \text{ for } \gamma > 1 \text{ are} \end{array}$$

subject to the no-flux boundary condition. Solutions of (52) for γ &gt; 1 are constructed by the direct method of calculus of variations.

Lemma 5. Let γ &gt; 1 . Then for any S ∈ L 2 (Ω) and for any pair of measurable disjoint sets A + , A -⊆ Ω there exists a unique (up to an additive constant) weak solution p ∈ H 1 (Ω) ∩ W 1 2 , γ/ γ ( -1) ( A + ∪ A -) of (52) .

Proof. Direct method of calculus of variations, see Theorem 6 of Ref. [8]. /square

In the case γ = 1, (45) reads

$$( \nabla p \otimes \nabla p ) m = \nu m, \\ \dots \dots \dots \dots \dots \dots.$$

/negationslash i.e., m is either the zero vector or an eigenvector of the matrix ∇ ⊗ ∇ p p with eigenvalue ν &gt; 0. The spectrum of ∇ ⊗∇ p p consists of zero and |∇ | p 2 , so that m = 0 is only possible if |∇ | p 2 = ν . Therefore, for every stationary solution there exists a measurable function λ = λ x ( ) such that

$$m ( x ) = \lambda ( x ) \chi _ { \{ | \nabla p | ^ { 2 } = \nu \} } ( x ) \nabla p ( x ) \\ \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \vdots \end{cases}$$

and p solves the highly nonlinear Poisson equation

$$\begin{array} { c } \cdot & \cdot & \cdot \\ & - \nabla \cdot \left [ \left ( 1 + \lambda ( x ) ^ { 2 } \chi _ { \{ c ^ { 2 } | \nabla p | ^ { 2 } = \nu \} } ( x ) \right ) \nabla p \right ] = S. & & ( 5 3 ) \\ \text{Ref. } [ 9 ], \, \text{Section } 4. 2, \, \text{it was shown that } ( 5 3 ) \text{ is equivalent to the free} \\ \text{ndary roblam} \end{array}$$

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

$$\nu \text{.} \\ - \nabla \cdot \left [ ( 1 + a ( x ) ^ { 2 } ) \nabla p \right ] \ & = \ S, \quad \ p \in H _ { 0 } ^ { 1 } ( \Omega ), \quad \quad ( 5 4 ) \\ | \nabla p ( x ) | ^ { 2 } \ & \leq \ \nu, \quad \text{ a.e. on } \Omega, \quad \quad ( 5 5 ) \\ a ( x ) ^ { 2 } \left [ | \nabla p ( x ) | ^ { 2 } - \nu \right ] \ & = \ 0, \quad \text{ a.e. on } \Omega, \quad \quad ( 5 6 )$$

for some measurable function a 2 = a x ( ) 2 on Ω which is the Lagrange multiplier for the condition (55). Solutions of (54)-(56) can be obtained as minimizers of the energy functional

$$a ( x ) ^ { 2 } \left [ | \nabla p ( x ) | ^ { 2 } - \nu \right ] \ = \ 0, \quad \text{ a.e. on } \Omega, \quad \text{ (566)} \\ \text{me measurable function } a ^ { 2 } = a ( x ) ^ { 2 } \text{ on } \Omega \text{ which is the Lagrange mul-} \text{ for the condition } ( 5 \L ) \text{$\text{$Solution of $(5\L)/5\L)$ on $h$ obtained $e$}$$

$$\mathcal { F } [ p ] & \coloneqq \int _ { \Omega } \left ( \frac { | \nabla p | ^ { 2 } } { 2 } - S p \right ) \, \mathrm d x \\ H ^ { 1 } ( \mathcal { O } \ \mathcal { \ } ^ { 2 } | \nabla n | ^ { 2 } < 1 \ \mathcal { a } \ \mathcal { o } \ \mathcal { O } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \matrix \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathfrak { \ } \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathcal { o } \ \mathfrak { \ } \mathcal { o } \ \mathfrak { \ } \mathcal { o } \ \mathfrak { \ } \mathcal { o } \ \mathfrak { \ } \mathcal { o } \ \mathfrak { \ } \mathcal { o } \ \mathfrak { \ } \mathcal { o } \ \mathfrak { \ } \mathcal { o } \ \mathfrak { \ } \mathcal { o } \ \mathfrak { \ } \mathcal { o } \ \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfra \ \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathi \ \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathcfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak {\ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathi } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathcrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \ \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \ \mathi \ \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ \ } \mathi \ \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak{ \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \cdot \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathf ] \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathftrack { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ }\mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \mathfrak { \ } \$$

on the set { p ∈ H 1 0 (Ω) , c 2 |∇ | p 2 ≤ 1 a.e. on Ω , or using the penalty method; } see Section 4.2 of Ref. [9] for details.

Finally, for 1 / 2 ≤ γ &lt; 1 the above approach fails due to the singularity of the term |∇ p x ( ) | 2 -γ γ -1 ∇ p x ( ) at |∇ | 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 |∇ | 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

$$\mathcal { E } [ m ] & = \int D ^ { 2 } | \nabla m | ^ { 2 } + \nabla p [ m ] \cdot \mathbb { P } [ m ] \nabla p [ m ] + \frac { \nu } { \gamma } | m | ^ { 2 \gamma } \, \mathrm d x, \quad \quad ( 5 7 ) \\ \text{where $D^{2}$ > $0$ is the diffisivity constant $n[m]$ denotes a solution of the}$$

where D 2 &gt; 0 is the diffusivity constant, p m [ ] denotes a solution of the Poisson equation (42) and P [ m ] = rI + m ⊗ 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

$$( \nabla p \otimes \nabla p ) m = - D ^ { 2 } \Delta m + \nu | m | ^ { 2 ( \gamma - 1 ) } m.$$

Note that for sufficiently regular critical points of (57), (42) we have ( -D 2 ∆ m + ν m | | 2( γ -1) m ) · m ≥ 0 almost everywhere in Ω. Let us denote

Then, there exist disjoint measurable sets A + , A -such that A + ∪A -= A and

$$\mathcal { A } \coloneqq \left \{ x \in \Omega ; \, \left ( - D ^ { 2 } \Delta m ( x ) + \nu | m | ^ { 2 ( \gamma - 1 ) } m ( x ) \right ) \cdot m ( x ) > 0 \right \}. \\ \text{Then, there exist disjoint measurable sets } \mathcal { A } ^ { + }, \, \mathcal { A } ^ { - } \text{ such that } \mathcal { A } ^ { + } \cup \mathcal { A } ^ { - } = \mathcal { A }$$

$$& \text{and} \\ & \nabla p = ( \chi _ { \mathcal { A } ^ { + } } - \chi _ { \mathcal { A } ^ { - } } ) \left ( \nu | m | ^ { 2 \gamma } - D ^ { 2 } m \cdot \Delta m \right ) ^ { - 1 / 2 } \left ( \nu | m | ^ { 2 ( \gamma - 1 ) } m - D ^ { 2 } \Delta m \right ), ( 5 8 ) \\ & \text{where} \ \chi _ { \mathcal { A } ^ { + } }, \, \text{resp.,} \ \chi _ { \mathcal { A } ^ { - } } \text{ are the characteristic functions of the sets $\mathcal{A} ^ { + }, \, \text{resp.,} \\ & \mathcal { A } ^ { - }. \ \text{On the set $\mathcal{A}$ we have $(\nabla p\otimes \nabla p) m = 0$, which is equivalent to $m\cdot\nabla p = 0$.} \\ & \text{Let us now choose an open subset $\mathcal{A}\subset\Omega$ with boundary $\partial\in C^{0,1}$ and} \\ & \text{$c_{\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\dots.}$$

Let us now choose an open subset Λ ⊂ Ω with boundary ∂ Λ ∈ C 0 1 , and define the inflow and outflow segments of ∂ Λ as

$$& \partial \Lambda ^ { + } \coloneqq \{ x \in \partial \Lambda ; \, n ( x ) \cdot q ( x ) < 0 \}, \quad \partial \Lambda ^ { - } \coloneqq \{ x \in \partial \Lambda ; \, n ( x ) \cdot q ( x ) > 0 \}, \\ & \text{where we denoted the flux } q \coloneqq - ( r I + m \otimes m ) \nabla p. \, \text{ Using } ( 5 8 ) \text{ in } ( 4 9 ) \text{ leads}$$

where we denoted the flux q := -( rI + m ⊗ m ) ∇ p . to

$$& \text{where we denoted the flux } q \coloneqq - ( r I + m \otimes m ) \nabla p. \, \text{Using } ( 5 8 ) \text{ in } ( 4 9 ) \text{ leads} \\ & \text{to} \\ & - \int _ { \mathcal { A } \cap \partial \L ^ { + } } \left ( \nu | m | ^ { 2 \gamma } - D ^ { 2 } m \cdot \Delta m \right ) ^ { - 1 / 2 } \left | n \cdot ( r I + m \otimes m ) \left ( \nu | m | ^ { 2 ( \gamma - 1 ) } m - D ^ { 2 } \Delta m \right ) \right | \, d s \\ & \quad - \int _ { \mathcal { A } \cap \partial \L ^ { + } } \, r | \nabla p \cdot n | \, d s \\ & + \int _ { \mathcal { A } \cap \partial \L ^ { - } } \left ( \nu | m | ^ { 2 \gamma } - D ^ { 2 } m \cdot \Delta m \right ) ^ { - 1 / 2 } \left | n \cdot ( r I + m \otimes m ) \left ( \nu | m | ^ { 2 ( \gamma - 1 ) } m - D ^ { 2 } \Delta m \right ) \right | \, d s \\ & \quad + \int _ { \mathcal { A } ^ { c } \cap \partial \L ^ { - } } \, r | \nabla p \cdot n | \, d s \\ & \quad = \int _ { \Lambda } S \, d x. \\ & \text{This is the Murray's law for the phenomenonological continuum model } ( 4 4 ) \\ & \text{constrained by the Poisson equation } ( 4 2 ).$$

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.

## References

- [1] G. Albi, M. Artina, M. Fornasier and P. Markowich: Biological transportation networks: modeling and simulation. Analysis and Applications. Vol. 14, Issue 01 (2016).
- [2] G. Albi, M. Burger, J. Haskovec, P. Markowich and M. Schlottbom: Continuum Modelling of Biological Network Formation. In: N. Bellomo, P. Degond, and E. Tamdor (Eds.), Active Particles Vol.I - Theory, Models, Applications , Series: Modelling and Simulation in Science and Technology, Birkh¨user-Springer (Boston), 2017. a
- [3] M.E. Alston and E.R. Barber: Leaf venation, as a resistor, to optimize a switchable IR absorber. Sci Rep. 6:31611 (2016).
- [4] M. Burger, J. Haskovec, P. Markowich and H. Ranetbauer: A mesoscopic model of biological transportation networks. Preprint, arXiv:1806.00120v2 (2018).
- [5] D.R. Emerson and K. Cie´ slicki, X. Gu, R.W. Barber: Biomimetic design of microfluidic manifolds based on a generalised Murray's law. Lab Chip. 6(3), pp. 447-54 (2006).
- [6] L.C. Evans: Partial differential equations. Graduate Studies in Mathematics Vol. 19, American Mathematical Society (2010).

- [7] L.C. Evans and R.F. Gariepy: Measure Theory and Fine Properties of Functions. CRC press (2015).
- [8] J. Haskovec, P. Markowich and B. Perthame: Mathematical Analysis of a PDE System for Biological Network Formation. Comm. PDE 40:5, pp. 918-956 (2015).
- [9] J. Haskovec, P. Markowich, B. Perthame and M. Schlottbom: Notes on a PDE system for biological network formation. Nonlinear Analysis 138 (2016), pp. 127-155.
- [10] D. Hu: Optimization, Adaptation, and Initialization of Biological Transport Networks. Notes from lecture (2013).
- [11] D. Hu and D. Cai: Adaptation and Optimization of Biological Transport Networks. Phys. Rev. Lett. 111 (2013), 138701.
- [12] J. Haskovec, L. M. Kreusser and P. Markowich: ODE and PDE based modeling of biological transportation networks . arXiv:1805.08526 (2018).
- [13] J. Haskovec, L. M. Kreusser and P. Markowich: Rigorous continuum limit for the discrete network formation problem. arXiv: (2018).
- [14] C. Murray: The Physiological Principle of Minimum Work: I. The Vascular System and the Cost of Blood Volume . PNAS 12 (3): 207-214 (1926). doi:10.1073/pnas.12.3.207.
- [15] C. Murray: The Physiological Principle of Minimum Work: II. Oxygen Exchange in Capillaries . PNAS 12 (5): 299-304 (1926). doi:10.1073/pnas.12.5.299.
- [16] M. Razavi, E. Shirani and G. Kassab: Scaling Laws of Flow Rate, Vessel Blood Volume, Lengths, and Transit Times With Number of Capillaries . Frontiers in Physiology 9:581 (2018).
- [17] D. Rogers: Laminar flow analysis. Cambridge University Press, 1992.
- [18] T. Sherman: On connecting large vessels to small. The meaning of Murray's law . The Journal of General Physiology. 78 (4): 431-453 (1981). doi:10.1085/jgp.78.4.431.
- [19] D. Stephenson, A. Patronis, D. Holland and D. Lockerby: Generalizing Murray's law: An optimization principle for fluidic networks of arbitrary shape and scale . J. Appl. Phys. 118:174302 (2015).
- [20] S. Whitaker: Flow in porous media I: A theoretical derivation of Darcy's law. Transport in porous media, 1 (1986), pp. 3-25.
- [21] X. Zheng et al: Bio-inspired Murray materials for mass transfer and activity. Nature Communications 8:14921 (2017).
Jan Haskovec Peter Markowich Giulia Pilli
Murray's law for discrete and continuum models of biological networks
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.