Modeling Damage in Self-Regulated Transportation Networks with Tree-Like Structure

1. Introduction

Transportation network structures that convey energy, people, goods, or other matters from one point to another appear in a variety of different fields, including the vascular vein networks forming the ramified patterns in plant leaves [1,2], river networks [3,4], the circulatory system [5,6,7], cells, a network of chemicals linked by chemical reactions, and the internet, a network of routers and computers connected by wired or wireless links [8], and city streets [9,10,11]. Recent empirical studies [9,12,13] have shown that, at least at a coarse grained level, unexpected quantitative similarities exist.

Street networks allow people and goods to move through cities, but they can be strongly affected by natural and extreme events. Flooding caused by heavy rainfall or sea-level rise can make streets inaccessible, while earthquakes may block roads with debris, complicating evacuations and emergency responses. These disruptions highlight the importance of designing street networks that are both robust, meaning able to withstand disruptions, and resilient, meaning able to recover and adapt quickly in time. The factors that determine how street networks respond to damage are still not fully understood. Recent large-scale studies, such as [14], have measured the vulnerability of urban street networks worldwide and shown that disruptions affect central nodes. More generally, networks with higher connectivity, greater redundancy, fewer checkpoints, and lower circuity tend to be less vulnerable to a wide range of disruptions. For example, in [15], the authors describe the formation of rings, webs and circuits as pre-defined structural elements of the network, while in [16], the authors recognize the presence of complex structures such as central loops together with peripheral branches, such as subway and railway networks, and investigate the conditions for the emergence of these nontrivial topological structures in the context of human transportation in cities.

In the context of seismic events, damage to buildings and short-term emergency measures can directly affect road functionality by reducing capacity or blocking streets [17]. The post-event response can be broadly divided into two phases. Immediately after, the primary concern is the ability of emergency services to reach critical locations. In the longer-term phase, road capacity may still be reduced due to temporary countermeasures, displaced populations and unusable buildings. These modifications of the urban network leads to changes in road serviceability, and highlights the need for modeling approaches that account for these disruptions that affect the mobility patterns. This topic is also highly relevant in the broader context of climate change, where transport infrastructure is increasingly exposed to extreme events. In recent years, climate change has become a multidisciplinary research topic addressing the challenges faced by transport infrastructure in planning, construction, and operation. The study of transport system adaptation to evolving environmental conditions is commonly framed through the interconnected concepts of resilience, vulnerability, and criticality [18]. We emphasize that these examples are intended to provide a general motivation for the study of network disruptions, rather than to suggest that the present model is specifically designed to describe such complex real-world scenarios.

Over the past decades, the study of networks has shifted from idealized regular and random graphs to the recognition that many real-world systems exhibit structured and non-random connectivity patterns [8]. Advances in data acquisition and computational power have revealed that networks arising in contexts such as biology, technology, and computer science possess complex topologies that cannot be adequately captured by purely random models. While graph-theoretical approaches have been instrumental in identifying these topological properties, they often provide a discrete and static description that is insufficient to capture spatial organization, continuous transport processes, and the connection between structures and functionalities. In this context, continuous modeling frameworks offer a natural and powerful alternative, allowing network formation and adaptation to be described through variational principles and dynamical laws. Such approaches highlight that the networks are driven by robust self-organizing mechanisms, and they provide quantitative tools to relate network topology to functional properties such as robustness and resilience. In this framework, the aim of this work is to investigate mechanisms leading to loop formation in continuous models, which are typically associated with tree-like structures. For this reason, the focus is on qualitative features rather than on quantitative agreement with specific real-world systems.

In this context, we start from an elliptic-parabolic system of partial differential equations (PDE) describing the formation of biological network structures. The system is derived as gradient flows of an energy functional, consisting of a diffusive term, an activation term and a metabolic cost term. The energy functional can be seen as a continuum version of its discrete counterpart introduced in [19]. Here, the authors consider a total energy consumption function for a general class of biological transport networks (seen also in [20]), including a material cost function of the network. In the papers [19,21], the authors study the adaptation of biological transport networks by linking the dynamics of local information to an underlying optimization principle, a mechanism that is known to be widespread in nature. Within this framework, such networks are analyzed through an optimization perspective in which the energy consumption is minimized under the constraint of a constant total material cost [20].

Assuming the validity of Darcy’s law for slow flow in porous media (see, for instance,  [22,23]), the energy functional is constrained by a Poisson equation for the fluid pressure [1,2,24]. In the parabolic equation, a reaction term appears, representing the metabolic cost of the network, while the diffusion term describes the randomness in the material structure. The third term, called pumping term, describes the tendency of the network to align with the principal direction of the material flow. The elliptic Poisson equation describes local mass conservation and is equipped with a right-hand side describing the distribution of sources of the material. This distribution is supplied as a datum and is supposed to be time independent.

Transport networks obtained from optimization strategies that minimize a gradient-flow energy are known to be trees, i.e. loopless [25,26]. In contrast, leaf networks exhibit a large number of closed loops, which allow fluid transport even after local damage. If the leaf networks were tree structures, a damage to a single vein would interrupt transport to all downstream regions, leading to the death of entire leaf sectors. Instead, the presence of loops enables the flow to be rerouted around injuries, including damage to the main vein [27].

In this work, we investigate how damage can be introduced into optimization-derived networks and how such networks can be repaired. Since standard optimization models do not generate loops, we consider damage introduced artificially through two mechanisms: (i) by introducing holes in the physical domain, and (ii) by augmenting the optimization model with an additional potential term.

For the discretization of the space derivatives of the system, we refer to a ghost-nodal finite element method (ghost-FEM), a very recent numerical scheme first introduced in [28] in the context of Poisson’s equations, and further applied to biological network formation [29], and Multiscale Poisson-Nernst-Planck systems [30,31]. A rectangular region R is discretized by a uniform square mesh, and the physical domain $\mathsf{\Omega }\subset R$ intersects the mesh, giving rise to the to the computational domain ${\mathsf{\Omega }}_h$. The finite element space is then defined on ${\mathsf{\Omega }}_h$. Despite the abundance of unfitted finite element schemes available in the literature, such as CutFEM in [32] and AgFEM in [33], the ghost-FEM is of interest for the following reason: when dealing with arbitrary domains, it is common to face small-cut cells problem. While other schemes rely on stabilization techniques, the ghost-FEM employs a very simple algorithm based on snapping-back-to-grid technique to overcome the ill-conditioning issue, and it is straightforward to implement it numerically.

The paper is organized as follows: in Section 2 we describe the mathematical model. In Section 3 we propose two strategies to insert prescribed damages to the network, motivating our choices. Section 4 presents the numerical experiments based on the proposed approaches. At the end, we draw some Conclusions. Since these methods have already been introduced in previous papers, we defer their presentation to the end of the manuscript in order to emphasize the damage modeling aspect.

2. Mathematical Model

We introduce the energy functional $\mathcal{E}$ for the model,

$$\mathcal{E}\left[\mathbb{C}\right]={\int }_{\mathsf{\Omega }}{\displaystyle \frac{D}{2}}{\left|\right|\nabla \mathbb{C}\left|\right|}^2+c^2\nabla p\left[\mathbb{C}\right]·(r\mathbb{I}+\mathbb{C})\nabla p\left[\mathbb{C}\right]+M\left(\right|\left|\mathbb{C}\right|\left|\right)\mathrm{d}x,$$

(1)

where $\mathbb{C}=\mathbb{C}\left(x\right)\in {\mathbb{R}}^{2\times 2}$ is the conductivity tensor, $\nabla =({\displaystyle \frac{\partial }{\partial x}},{\displaystyle \frac{\partial }{\partial y}})$ the divergence operator, $D\in \mathbb{R}$ the diffusivity parameter that measures the effect of random fluctuations in the network, and the activation parameter $c^2>0$ controls the strength of the network formation feedback loop. The scalar $r>0$ describes the isotropic background permeability of the medium, and is assumed to be constant. The scalar pressure $p=p\left[\mathbb{C}\right]$ of the fluid transported within the network is the unique solution (up to an additive constant) of the Poisson equation

$$-\nabla ·\left((r\mathbb{I}+\mathbb{C})\nabla p\right)=S,$$

(2)

subject to homogeneous Neumann boundary condition on $\partial \mathsf{\Omega }$. The source/sink distribution $S=S\left(x\right)$, in the mass conservation Eq. (2), is considered as an input datum. The metabolic cost function $M:{\mathbb{R}}^{+}\to {\mathbb{R}}^{+}$ describes the dependence of the metabolic expenditure to maintain the network on its transport capacity, see [19]. The expression $\left|\right|\mathbb{C}\left|\right|$ denotes the Frobenius norm of $\mathbb{C}$.

Taking the $L^2$-gradient flow of the energy $\mathcal{E}\left[\mathbb{C}\right]$ in (1), constrained by (2), we obtain the parabolic-elliptic system

$$\begin{array}{cc}\hfill -\nabla ·\left(\left(r\mathbb{I}+\mathbb{C}\right)\nabla p\right)& =S\hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill {\partial }_t\mathbb{C}-\underset{\mathrm{randomness}}{\underbrace{D\Delta \mathbb{C}}}-\underset{\mathrm{pumping}\mathrm{term}}{\underbrace{c^2\nabla p\otimes \nabla p}}+\underset{\mathrm{metabolic}\mathrm{term}}{\underbrace{{\displaystyle \frac{M^{\prime }\left(\left|\right|\mathbb{C}\left|\right|\right)}{\left|\right|\mathbb{C}\left|\right|}}\mathbb{C}}}& =0,\hfill \end{array}$$

(3b)

see [2] for details of the derivation. The pumping/activation term $c^2\nabla p\otimes \nabla p$ in (3b) is quadratic and depends on $\mathbb{C}$ only through the solution p of the Poisson equation (3a). In particular, $\Delta ={\nabla }^2$ and ⊗ is the tensor product to create a matrix from two vectors. Moreover, the metabolic term becomes singular at $\mathbb{C}=0$ if $\gamma <1$. This obviously causes difficulties for both the analytical treatment of the system and its numerical resolution.

When considering biological networks, a common choice for the metabolic cost function (see, e.g., [1,24,34]) is the following,

$$M\left(s\right):={\displaystyle \frac{\nu }{\gamma }}{s}^{\gamma }\mathrm{for} s\ge 0,$$

(4)

where $\nu >0$ is the metabolic constant and $0<\gamma <1$ the metabolic exponent; see [19]. In this work, we adopt this choice; however, the proposed numerical scheme is flexible and can accommodate different metabolic cost functions. This will allow us, in future work, to consider alternative expressions that may be more appropriate for applications to urban networks.

In this work, we follow the simplified approach in [26], where the scalar case is used to highlight the intricate network structure of the results, or in [35], where the authors show that the scalar model has a gradient flow structure whose steady-state coincides with the solution of the p-Laplacian equation. Thus, with the assumption $\mathbb{C}\in \mathbb{R}$, system (3) becomes

$$\begin{array}{cc}\hfill \nabla ·\left((r+\mathbb{C})\nabla p\right)& =S\hfill \end{array}$$

(5a)

$$\begin{array}{cc}\hfill {\partial }_t\mathbb{C}-D\Delta \mathbb{C}-c^2\mathcal{P}+\nu {\left|\mathbb{C}\right|}^{\gamma -2}\mathbb{C}& =0\hfill \end{array}$$

(5b)

where $\mathcal{P}={\partial }_x{p}^2+{\partial }_y{p}^2.$ In our tests, the expression of the source term is the following

$$S=exp(-\omega ({(x-x_S)}^2+{(y-y_S)}^2)),\omega =1000,x_S,y_S\in \mathsf{\Omega }.$$

We define the equations on a bounded domain $\mathsf{\Omega }\subset {\mathbb{R}}^2$, and we define the boundary conditions on $\Gamma :=\partial \mathsf{\Omega }$ for the pressure and the conductivity. We choose homogeneous Neumann conditions for the two variables $\mathbb{C}$ and p, such that,

$$\nabla \mathbb{C}(t,\overrightarrow{x})·\widehat{n}=0,\nabla p(t,\overrightarrow{x})·\widehat{n}=0,\overrightarrow{x}\in \Gamma ,t\ge 0,$$

(6)

where $\widehat{n}$ is the outgoing normal vector to $\Gamma$. A direct consequence of the boundary condition defined for the pressure, that is also a necessary condition for the solvability of the Poisson equation, is that the source function has to be vanishing mean, i.e., ${\int }_{\mathsf{\Omega }}Sd\mathsf{\Omega }=0$.

To close the system, we prescribe an initial condition for $\mathbb{C}$

$$\mathbb{C}(t=0,\overrightarrow{x})={\mathbb{C}}^0\left(\overrightarrow{x}\right),\mathrm{in}\mathsf{\Omega }.$$

(7)

Figure 1. Plot of the solution $\mathbb{C}$ in a squared domain, for different choices of the diffusion coefficient D, and of the stabilization parameter ε. Other parameters: $c=15,\nu =\gamma =0.65,x_S=y_S=0.2,h=1/800,\Delta t=10h$ and final time $t_{\mathrm{fin}}=52$.

Figure 1. Plot of the solution $\mathbb{C}$ in a squared domain, for different choices of the diffusion coefficient D, and of the stabilization parameter ε. Other parameters: $c=15,\nu =\gamma =0.65,x_S=y_S=0.2,h=1/800,\Delta t=10h$ and final time $t_{\mathrm{fin}}=52$.

3. Modeling the Damage

In this section, we present two strategies to model known damages in a network. In the first approach, we introduce prescribed defects at specific, known locations within the computational domain, representing the response of an already formed network to a sudden localized disruption. In the second approach, the network develops in the presence of a pre-existing defect in the domain, thus forming while already accounting for the disruption.

Figure 2. Representation of the steps described in Algorithm 1. We solve system (5) until the energy decays, see panels (a) and (b). Panels (c) and (d) show the resulting network at time $t=50$. Holes are then introduced, as described in Algorithm 1, see panels (e) and (f). The system (5) is subsequently solved on the resulting perforated domain until time $t=t_{\mathrm{fin}}=52$, and the corresponding network is shown in panels (g) and (h). Parameters: $D={10}^{-5},c=15,\gamma =0.75,\nu =\gamma ,x_S=y_S=0.5,h=1/800$ and $\Delta t=10h$.

Figure 2. Representation of the steps described in Algorithm 1. We solve system (5) until the energy decays, see panels (a) and (b). Panels (c) and (d) show the resulting network at time $t=50$. Holes are then introduced, as described in Algorithm 1, see panels (e) and (f). The system (5) is subsequently solved on the resulting perforated domain until time $t=t_{\mathrm{fin}}=52$, and the corresponding network is shown in panels (g) and (h). Parameters: $D={10}^{-5},c=15,\gamma =0.75,\nu =\gamma ,x_S=y_S=0.5,h=1/800$ and $\Delta t=10h$.

3.1. Modified Level-Set (MLS) Approach

We assume the square region R is such that $\mathsf{\Omega }\subset R$. We discretize it by a uniform Cartesian mesh of size $h=1/N$, where $N\in \mathbb{N}$ is the number of cells in each space direction. Following the approach shown in [36,37], the open domain $\mathsf{\Omega }$ is implicitly defined by a level set function $\varphi (x,y)$ that is negative inside $\mathsf{\Omega }$, positive in $R\setminus \mathsf{\Omega }$ and zero on the boundary $\Gamma$:

$$\begin{array}{c}\hfill \mathsf{\Omega }=\left\{\right(x,y):\varphi (x,y)<0\},\Gamma =\left\{\right(x,y):\varphi (x,y)=0\}.\end{array}$$

(8)

The outer unit normal vector $\widehat{n}$ in (6) can be computed as $\widehat{n}={\displaystyle \frac{\nabla \varphi }{|\nabla \varphi |}}$. We observe that, for a given circle with center $(x_c,y_c)=(0,0)$ and radius R, there exist infinitely level-set functions. For instance, the functions ${\varphi }_1=R-\sqrt{{(x-x_c)}^2+{(y-y_c)}^2}$ and ${\varphi }_2=R^2-\left({(x-x_c)}^2+{(y-y_c)}^2\right)$ describe the same circle. Among all possible representations of a given bubble, the most convenient level-set function from a numerical stability standpoint is the signed distance function ${\varphi }_1$ defined as the distance between $(x,y)$ and the circumference (positive inside the circle and negative outside).

As already mentioned in the Introduction, the absence of loops in these networks is a direct consequence of the optimization framework described above, which selects tree-like structures by minimizing energy consumption subject to a fixed total material cost. For this reason, we resort to an alternative approach. In this subsection, we modify the computational domain by introducing holes to model damages in the network. To introduce other circular holes, we define other smaller circles in the same way, and then considering the maximum of the two resulting level-set functions. In Algorithm 1, we describe the procedure we follow. The idea is to solve the linear system (5) until the energy decays. In practice, we fix a tolerance and compute the solution, $\mathbb{C}$, until the variation of the energy, $\mathcal{E}\left[\mathbb{C}\right]$, between two consecutive time steps is smaller than this tolerance. At that point, we modify the level set function, introducing a certain number of holes or damages, and restart the computation of the solution of the linear system (5) until final time. Imposing homogeneous Neumann boundary conditions to the boundary of the new holes, do not modify the variational formulation. From a computational viewpoint, this imposition do not alter the numerical scheme, as the classification of grid points automatically sets the new boundary conditions at the corresponding ghost cells.

Remark 1. 

We observe that modifying the level-set function changes the computational domain and the corresponding boundary Γ, that now accounts for the presence of holes. We impose homogeneous Neumann boundary conditions also on the boundaries of the newly introduced holes, i.e.

$$\nabla \mathbb{C}(t,\overrightarrow{x})·\widehat{n}=0,\overrightarrow{x}\in \Gamma ,$$

where $\Gamma =\partial \mathsf{\Omega }\cup \partial {\mathcal{C}}_k,k=1,\dots ,\#holes$. With this choice, the variational formulation remains unchanged, and the numerical scheme can be applied without modification. Consequently, this approach allows us to incorporate the presence of holes into the simulations with only minimal adjustments to the code, and specifically to the level-set definition.

Following this strategy, we model prescribed damages occurring at specific, known locations of the computational domain. This setting can represent, for instance, the response of an already formed network to a sudden localized disruption, allowing us to investigate how the structure reorganizes in the presence of an externally imposed defect.

When holes are introduced at specific times and locations, part of the solution may be lost. To assess this effect, we compute the integral of the solution before and after the introduction of the holes. In Figure 4, we report the relative difference ${\mathrm{diff}}_{\mathrm{rel}}$, defined in Eq. (9), between the integral of the solution with and without holes at the same final time.

$${\mathrm{diff}}_{\mathrm{rel}}^k={\displaystyle \frac{\left|\right|{\int }_{\mathsf{\Omega }(k-\mathrm{holes})}\mathbb{C}(t=t_{\mathrm{fin}})d\mathsf{\Omega }-{\int }_{\mathsf{\Omega }(0-\mathrm{holes})}\mathbb{C}(t=t_{\mathrm{fin}})d\mathsf{\Omega }{\left|\right|}_2}{\left|\right|{\int }_{\mathsf{\Omega }(0-\mathrm{holes})}\mathbb{C}(t=t_{\mathrm{fin}})d\mathsf{\Omega }{\left|\right|}_2}},k=1,\dots ,\#\mathrm{holes}.$$

(9)

The results show that this difference increases as the total surface of the holes increases. Moreover, not all holes affect the solution in the same way. Disruptions introduced at locations where the value of the solution $\mathbb{C}$ is larger, e.g. point D, produce a larger difference than disruptions at locations where the solution is smaller, e.g. point C.

Figure 3. Panel (a) shows relative difference between the integrals of the solutions, defined in Eq. (9), computed with and without holes. Specifically, the horizontal axis reports the cumulative area (from 1 to 4), while the vertical axis shows this relative difference for configurations with 1, 2, 3, and 4 holes. The order of the introduction of holes is illustrated in Figure 2 (e),(f). Zoom-in in panel (b).

Figure 3. Panel (a) shows relative difference between the integrals of the solutions, defined in Eq. (9), computed with and without holes. Specifically, the horizontal axis reports the cumulative area (from 1 to 4), while the vertical axis shows this relative difference for configurations with 1, 2, 3, and 4 holes. The order of the introduction of holes is illustrated in Figure 2 (e),(f). Zoom-in in panel (b).

Algorithm 1 Computation of the modified level-set function.

$\varphi =R-\sqrt{{(x-x_c)}^2+{(y-y_c)}^2}$ while $\mathcal{E}\left[{\mathbb{C}}^n\right]-\mathcal{E}\left[{\mathbb{C}}^{n+1}\right]>{10}^{-5}$ do      solve system (5) in $\mathsf{\Omega }$(0-holes) end while for $k=1,\dots ,\#$ holes     do ${\varphi }_k:=\sqrt{{(x-x_1)}^2+{(y-y_1)}^2}-R_k$     $\varphi =max\{\varphi ,{\varphi }_k\}$ end for solve system (5) in $\mathsf{\Omega }$(k-holes) till final time

3.2. Potential Term

In [38], the authors say that the adaptation dynamics minimize the global energy consumption, leading to optimal networks that may develop hierarchical loop structures in the presence of strong fluctuations in the flow distribution. Motivated by this observation, an alternative way to model the creation of loops is to augment the optimization framework with an additional term, as follows

$$\begin{array}{cc}\hfill \nabla ·\left((r+\mathbb{C})\nabla p\right)& =S\hfill \end{array}$$

(10a)

$$\begin{array}{cc}\hfill {\partial }_t\mathbb{C}-D\Delta \mathbb{C}-c^2\mathcal{P}+\nu {\left|\mathbb{C}\right|}^{\gamma -2}\mathbb{C}-\lambda W& =0\hfill \end{array}$$

(10b)

where $\lambda =0.01$ and the function $W=W(x,y)$ has the following expression

$$W=exp(-\omega ({(x-x_W)}^2+{(y-y_W)}^2)),x_W=y_W=0.5.$$

(11)

Figure 4. Contour line plot of the external potential $W(x,y)$ introduce in Eq. (11).

Figure 4. Contour line plot of the external potential $W(x,y)$ introduce in Eq. (11).

Several numerical experiments show that the bifurcation of the network appears only with $x_W=y_W$. In other words, the loop forms when the potential W is centered at a point located in the main diagonal. This feature is consistent with the rest of the numerical results obtained with the MLS approach, where we are sure that the network creates loops around the holes (or circles) centered at the main diagonal.

4. Numerical Experiments

In this section, we present the numerical tests. To represent different possible city layouts, we consider domains with various geometries, namely circular, square, and rectangular shapes.

System (5) depends on several parameters, leading to different network configurations. In our numerical simulations, we illustrate the effect of these parameters. Some of these aspects have already been investigated in previous works (see, e.g., [25,26]). Figure 1 shows that smaller values of the diffusion coefficient D lead to thinner ramifications. We also investigate the role of the stabilization parameter $\epsilon$, introduced in Appendix B. As shown in [25], the parameter that most strongly influences the network formation is the background permeability r, which also affects the asymmetry of the network. The parameter $\epsilon$ is introduced to prevent division by zero in the metabolic term. However, it negatively affects the conditioning of the linear system that we solve at each time step. For this reason, it must be chosen carefully: in our simulations, we set $\epsilon <r$, trying to avoid excessively small values that would deteriorate the conditioning of the system. Figure 1 (b-c) show the results for $\epsilon ={10}^{-4},{10}^{-5}.$ In Figure 5, we observe that the metabolic exponent $\gamma$ influences the distribution of branches. For $\gamma =0.75$, the branches are more uniformly distributed across the domain.

For the circular domain, we take $R={[0,1]}^2$ with the source located at $x_S=y_S=0.5$. For the square domain, we again consider $R={[0,1]}^2$, with $x_S=y_S=0.2$. For the rectangular case with (the sum of) two sources, we set $R={[0,1.5]}^2$, with $x_{S1}=0.15$, $y_{S1}=0.4$, $x_{S2}=1.35$, and $y_{S2}=y_{S1}$. Finally, for the square domain with (the sum of) four sources, we define $R={[0,1]}^2$ with $x_{S1}=y_{S1}=0.1,x_{S2}=0.6,y_{S2}=0.1,x_{S3}=0.1,y_{S3}=0.6,x_{S4}=y_{S4}=0.6.$ In the following simulations, the initial condition in Eq. (7) is set to be constant, i.e., ${\mathbb{C}}^0\left(\overrightarrow{x}\right)={\mathbb{C}}^0=1.$

Figure 6. Comparison of the resulting network in a squared domain, with increasing number of damages/interruptions: three in panel (a), six in panel (b) and eight in panel (c). Parameters: $D={10}^{-5},c=15,\nu =\gamma =0.65,x_S=y_S=0.2,h=1/800,\Delta t=10h$ and final time $t_{\mathrm{fin}}=52$.

Figure 6. Comparison of the resulting network in a squared domain, with increasing number of damages/interruptions: three in panel (a), six in panel (b) and eight in panel (c). Parameters: $D={10}^{-5},c=15,\nu =\gamma =0.65,x_S=y_S=0.2,h=1/800,\Delta t=10h$ and final time $t_{\mathrm{fin}}=52$.

In Figure 2, we present results for a circular domain and illustrate the different steps of the Algorithm 1 used to generate holes in the domain. The idea is to introduce holes only after the solution has reached its steady state. To this end, we first wait until the energy reaches its minimum (panels (a)-(b)). We inspect the solution (panels (c)-(d)) and select ad-hoc points (A-D in panel (e), E-H in panel (f)) where holes are then introduced. The simulation is subsequently continued, and the resulting configurations are shown in panels (g)-(h). The results indicate that the solution adapts by flowing around the holes while preserving its original structure. The network is affected locally by the introduced damages, but remains essentially unchanged elsewhere. This suggests that, when the perturbations are ’sufficiently’ small, the network exhibits a self-repair mechanism, reorganizing itself to bypass the damaged regions while minimally impacting the overall structure.

In Figure 5, we present the results for a squared domain. As already observed in [26], plotting the logarithm of the solution (i.e., $log\left(\mathbb{C}\right)$) allows us to reveal finer ramification patterns that are not visible when displaying the solution itself (i.e., $\mathbb{C}$). Moreover, we consider larger holes in order to further assess the adaptability of the model and of the numerical scheme. In agreement with previous observations, the network adjusts by deviating as little as possible from its original configuration, even in the presence of more significant damage. This suggests that it is preferable, and less costly, to follow existing paths instead of creating new ones. However, the increased size of the holes induces more evident changes: in some regions, certain branches disappear (e.g., $x\approx 0.3,y\approx 0.6$, or $x\approx 0.5,y\approx 0.2$, from panel (c) to panel (e); also, $x\approx 0.6,y\approx 0.4$, from panel (d) to (g)), while new ones emerge to provide alternative pathways same region, panels (e,g). This reorganization reflects the same self-adaptive mechanism observed before, now acting to compensate for the stronger perturbations. Moreover, as expected, the newly formed branches do not reconnect with the existing ones, in agreement with the tree structure of the network. Another characteristic of the network is the presence of a preferred direction, namely the main diagonal when the source is placed in $x_S=y_S=0.2$. When holes are introduced along this direction, the network consistently adapts by flowing around them, even when we introduce significant damage. This feature is further highlighted in Figure 7, where we show the results of the second strategy introduced in Section 3.2. Among all the different numerical tests we performed, the only configuration in which a loop forms in the network occurs when the external potential $W(x,y)$ is placed in the main diagonal (see Figure 7). Figure 8 shows the same simulation for different values of $\lambda =0.1,0.01,0.001$. We see that for $\lambda =0.001$ the potential term is not strong enough to modify the network; on the contrary, with $\lambda =0.1$, the bifurcation is much larger than $\lambda =0.01$.

To investigate the impact of increasingly complex damage configurations, we consider a set of numerical experiments with a growing number of damaged regions. In Figure 6, we report a comparison between cases with three, six, and eight damages. Moreover, in Table 1, we compute the integral of the solution, showing that it increases as the number of damages grows. This behavior reflects the formation of additional branches and alternative pathways in the network. The extension to spatially varying damage intensity, representing a continuum of degradation levels, would require a more general formulation (e.g., through spatially distributed coefficients or potential terms) and is left for future investigation.

In Figure 9, we present the results in a rectangular domain. This configuration differs from the previous ones, as we introduce two sources. The results show that the domain appears to be split into two regions, each characterized by a distinct network, which are symmetric with respect to the axis $x=0.75$. It is important to emphasize that this partition of the domain is not imposed a priori, but it emerges naturally as a consequence of the presence of the two sources. Moreover, introducing holes, we are able to break the symmetry of the two networks, since they adapt according to the position and dimension of the damages. Similarly, in Figure 10, the domain appears to be divided into four regions, corresponding to the introduction of four sources. We show the location of the four sources in Figure 10 (b), and the corresponding network in panel (a). It is evident that the preferred direction is along the diagonal $y=x$, where the two leading networks develop. The other two networks (up-left and down-right) appear to be symmetric with respect to the same diagonal. This behaviour can be explained by the distribution of the sources: the two located on the right generate a network on their right side, since no other sources appear there. On the other hand, the two left sources spread on the remaining domain, resulting in a domain that is automatically divided into four parts.

In Figure 10, no holes are introduced. Nevertheless, this test shows an alternative strategy to handle possible disruptions that could affect the main paths coming from a single source (see Figure 5). In this test, we introduce multiple sources that cover different regions of the same domain, ensuring a more distributed transport mechanism.

5. Conclusions

In this work, we address the limitation of self-regulated transportation networks, which are typically tree-like structures and therefore cannot develop loops or alternative routes. Our goal is to introduce a general continuous framework capable to reorganize in the presence of damage in order to maintain functionality. For this reason, the focus is on qualitative features rather than on agreements with specific real-world systems. To overcome this limitation, we have proposed two different strategies to model damages in the domain. In the first approach, multiple holes are introduced after the solution has reached a steady state, allowing us to study the response of the system to sudden perturbations of the established network structure.

In the second approach, the location of the damaged region is prescribed before the beginning of the simulation. The numerical results show that the network is able to adapt effectively when damages occur at known and prescribed locations (A-D, E-H, see Figure 2 panels (e-f)). In particular, it reorganizes its structure and identifies alternative paths that preserve the transport of material. This suggests that, in applications such as energy or resource distribution, the system is robust against interruptions by being able to activate secondary routes. Moreover, we observe that the introduction of multiple sources in the domain leads to a natural partition of the domain into regions where transport occurs in parallel. This provides an efficient mechanism for distributing quantities across different parts of the domain. As future work, we aim to extend the model by coupling the system with an additional equation for the potential $W(x,y)$. This would enable a fully self-consistent description of both transport and damage formation. In addition, the observed tendency for loop formation along the diagonal direction is not yet fully understood and will require further investigation.

Moreover, future research will focus on extending the proposed framework to more realistic and application-driven scenarios. In particular, a promising direction is the integration of the model with real infrastructure systems, such as power grids, water distribution networks, and urban transportation systems, where resilience to localized failures is critical. This is particularly relevant for simulation environments in smart city contexts, such as the one presented in [39]. However, the proposed models have potential applications across multiple domains where resilient and adaptive transport systems are needed that will be investigated in future contributions: a) in biological contexts, such models can be used to better understand processes such as leaf patterns, vascular system, neural plasticity, and tissue regeneration; b) in engineering, they may contribute the design of robust urban transportation networks and communication infrastructures capable of maintaining functionality under failures or targeted attacks.