# Self-Organized Homogenization of Flow Networks

Julien Bouvard, Swarnavo Basu, Charlott Leu, Onurcan Bektas, Joachim O Rädler, Gabriel Amselem, Karen Alim

► To cite this version:

Julien Bouvard, Swarnavo Basu, Charlott Leu, Onurcan Bektas, Joachim O Rädler, et al. Self-Organized Homogenization of Flow Networks. Physical Review X, 2025, 15 (4), pp.041038. 10.1103/j5ch-4vkh. hal-05469456

HAL Id: hal-05469456

https://hal.science/hal-05469456v1

Submitted on 21 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Self-organized homogenization of flow networks

Julien Bouvard,¹, * Swarnavo Basu,², * Charlott Leu,³ Onurcan

Bektas,²,³,⁴ Joachim O. Rädler,³ Gabriel Amselem,¹, † and Karen Alim², ‡

¹Laboratoire d'Hydrodynamique (LadHyX), CNRS, Ecole Polytechnique,
Institut Polytechnique de Paris, 91120 Palaiseau, France

²Center for Protein Assemblies and Department of Bioscience, School of
Natural Sciences, Technische Universität München, 85748 Garching, Germany

³Faculty of Physics and Center for NanoScience, Ludwig-Maximilians-Universität München, 80539 Munich, Germany

⁴Max Planck School Matter to Life, Max Planck Institute for Medical Research, Heidelberg, 69120, Germany

(Dated: January 21, 2026)

From the vasculature of animals to the porous media making up batteries, the core task of flow networks is to transport solutes and perfuse all cells or media equally with resources. Yet, living flow networks have a key advantage over porous media: they are adaptive and self-organize their geometry for homogeneous perfusion throughout the network. Here, we show that also artificial flow networks can self-organize toward homogeneous perfusion by the versatile adaption of controlled erosion. Flowing a pulse of cleaving enzyme through a network patterned into an erodible hydrogel, with initial channels disparate in width, we observe a homogenization in channel resistances. Experimental observations are matched with numerical simulations of the diffusion-advection-sorption dynamics of an eroding enzyme within a network. Analyzing transport dynamics theoretically, we show that homogenization only occurs if the pulse of the eroding enzyme lasts longer than the time it takes any channel to equilibrate to the pulse concentration. The equilibration time scale derived analytically is in agreement with simulations. Lastly, we show both numerically and experimentally that erosion leads to the homogenization of complex networks containing loops. Erosion being an omnipresent reaction, our results pave the way for a very versatile self-organized increase in the performance of porous media.

# I. INTRODUCTION

Whether in engineered or natural systems, fluid flows perfuse disordered porous media, from packed-bed reactors, filters, bio-engineered tissues and fuel cells, to rocks, geological sediments, and irrigated soils [1–5]. The disordered structure of their networks’ channels results in heterogeneous flows [6–14], diminishing perfusion in network channels of high resistance [15]. The consequences of the lack in perfusion are dramatic, reducing the efficiency of packed-bed reactors, filters and fuel cells and the survival of cells in bio-engineered tissues [16].

In contrast, in biological flow networks such as functional mycelia, slime mold veins, or our own microvasculature [17–19], network channels adapt such that network geometry is self-organized to transport solutes throughout the network, perfusing all parts of the organism with nutrients and signaling molecules [20–27]. Such perfusion is best achieved when the velocity of the fluid is homogeneous in all channels of a network like in intermeshed channels of equal resistance [20–22, 25, 26, 28, 29]. Such uniformity, though, is contrary to the natural tendency to generate disorder and thus inhomogeneities [13]. The key to living flow networks’ self-organization is the local adaptation of channels to local flow rate [30–32] or local solute concentrations [20, 28, 33, 34]. Such local

adaptation is sufficient to drive global self-organization of the flow networks as flows are globally coupled due to conservation of fluid volume: changing locally a channel width changes the local resistance to flow, which overall redistributes fluid throughout the entire network. Thus, even if inanimate flow networks cannot rise to the complexity of local adaptation as living networks, any local adaptation will hijack the global coupling of flows.

Engineered flow networks operating at dimensions similar to those of living networks are readily reproduced using microfluidic tools [35–37]. Integrating responsive materials renders microfluidic chips adaptive to local solutes controlled by transport within fluid flows [38, 39]. In particular, the swelling and contraction of hydrogels by local solute concentration, or temperature [40, 41], has so far only been employed to operate switches, closing or opening a channel [38, 42–46]. Yet, carving entire networks of channels into adaptive hydrogels opens the potential for self-organizing flow networks.

Here, we combine experiments, theory, and numerical simulations to develop adaptive inanimate flow networks capable of self-organization for homogeneous flow independent of their initial network geometry by controlled erosion. We build microfluidic channel networks out of an erodible hydrogel and let the device geometry progressively adapt as pulses of erosive chemicals flow into the chip. When the chemical diffuses into the hydrogel walls, these erode, which changes the hydraulic resistances of the channels and, thus, the overall flow pattern in the network. Experiments are successfully described by numerical solutions of the advection-diffusion-sorption

<footnote>* These authors contributed equally to this work.</footnote>

<footnote>† Contact author: gabriel.amselem@polytechnique.edu</footnote>

<footnote>‡ Contact author: k.alim@tum.de</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1. Erosion by the MMP-1 enzyme homogenizes the resistances of parallel channels whose walls are made of the hydrogel PEG-NB, cross-linked with a cleavable peptide sequence. (a) Experimental setup: glass slide with six microfluidic chambers. Each chamber contains a different microfluidic network, whose hydrogel walls are made of PEG-NB cross-linked with a cleavable peptide sequence. Inset: four-channel network. (b) Time-lapse of the erosion of a four-channel network, whose hydraulic resistances are initially imbalanced. The images, taken at t = 0 (top), t = 5 min (middle), and t = 8 min (bottom), show the hydrogel (light green) being eroded by the enzyme (red). Scale bar: 1 mm. (c) Time-evolution of the widths $w_i(t)$ of the four channels. The injection of the enzyme upstream of the network occurs from t = 2 min to t = 7 min (light red patch). Points: experimental data. Lines: numerical simulations. (d) Homogenization of the channels' hydraulic resistances $R_i(t)$ as the channel walls are getting eroded: the hydraulic resistances $R_i(t)/R_4(t)$ normalized by the narrowest channel's resistance homogenize as they increase over time. Channels are numbered from the largest to the smallest width. Points: experimental data. Lines: numerical simulations. Inset: normalized standard deviation $\sigma(t)/\sigma(0)$ of the normalized hydraulic resistances $R_i(t)/R_4(t)$. Similarly to the hydraulic resistances, the flow rates also homogenize, see Supp. Fig. 11 [47].</figcaption></figure>

dynamics of the flow-transported enzyme within the network. The theoretical description predicts homogenization by erosion to generally hold for any flow network as long as the pulse length is longer than the time for equilibration of enzyme concentration in any channel. We analytically derive the enzyme equilibration time scale as an easy-to-apply rule to achieve homogenization by controlled erosion. Successful homogenization of complex loopy microfluidic networks confirms the predicted general applicability. Our work paves the way for the design of a novel class of adaptive flow devices and versatile optimization of porous media performance by such a simple adaptation as erosion.

## II. RESULTS

### A. Erosion homogenizes resistances of parallel channels

We begin by considering an imbalanced network of four parallel channels of identical height ($H = 400$ µm) and length ($L = 4250$ µm), but whose widths $w_i$ vary from 210 to 545 µm, see Fig. 1a. The top and bottom of the channels are sealed with plastic, while their side walls consist of the hydrogel PEG-NB, cross-linked with a cleavable peptide sequence (see Methods). To probe the impact of controlled erosion, we inject a $t_e = 5$ min long pulse of the eroding enzyme MMP-1 at a concentration of $c_0 = 3 \times 10^{-5}$ mol L$^{-1}$ and a flow rate of $Q = 20$ µL min$^{-1}$, followed by the injection of Phosphate Buffer Saline (PBS) for 60 min. The boundaries of the

channel wall and the dynamics of the enzyme are monitored under a fluorescence microscope. The walls fluorescence in green and the enzyme solute in red due to the incorporation of a green fluorescent dye into the hydrogel mesh, and the mixing of the MMP-1 enzyme with a red fluorescent dye of the same molecular weight, Texasred dextran (see Methods). Snapshots of typical experiments are shown in Fig. 1b for a four-channel network and in Supp. Fig. 2a [47] for a single channel one (see also Movie 1 in the Supplemental Material [47]).

When the eroding solution is injected, the enzyme concentration increases both in the channels and, due to the enzyme's diffusivity, within the hydrogel walls. In the presence of the enzyme, the peptide sequence cross-linking the hydrogel is cleaved, leading to the erosion of the channel walls and, thus, to an increase in channel widths (see Fig. 1c). Flowing PBS into the network flushes away the eroding solution and gradually slows down the erosion.

To address how erosion affects the imbalance of the channels, we calculate the evolution of the hydraulic resistances of the four parallel channels from measurements of the channels' widths as a function of time. The hydraulic resistance for a rectangular channel of length $L$, width $w$ and height $H > w$, is given by

$$R = \frac{12\eta L}{Hw^3} \left(1 - 0.63\frac{w}{H}\right)^{-1}, \quad (1)$$

where $\eta$ is the fluid's viscosity [49]. Normalizing all channel resistances by the resistance of the narrowest channel in the network at any point in time reveals that
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Wall erosion follows Michaelis-Menten kinetics. (a) Spatio-temporal evolution of the channel width *w* of a single channel, see Movie 1 in the Supplemental Material [47], at a fixed position along the channel. The hydrogel (light green) is eroded by the enzyme (red, $c \ge 1.5 \times 10^{-6} \text{ mol L}^{-1}$) diffusing into the hydrogel and back into the channel after the pulse passed. The hydrogel walls, i.e. the boundaries between the channel and the hydrogel, are highlighted in blue. Three pulses of enzyme of duration $t_e = 5$ min are inter-spaced by 60 min of PBS flowing into the network. (b) Time-evolution of the channel width *w* (blue) as pulses of the enzyme are flowed in. The concentration of the enzyme in the channel (red) is averaged along the cross-section of the channel in *y* and along the channel in *x*, $\langle \hat{c} \rangle$. (c) Erosion rate $dw/dt$ of the hydrogel, as a function of the local enzyme concentration $c_w$ within 50 µm of the wall, follows Michaelis-Menten kinetics. Squares: raw data. Triangles: averaged data. Error bars show the standard deviation. Line: power-law fit of the data. Data over 10 experiments with 10 different channel geometries but keeping the PBS concentration and channel heights constant. The kinetic model is consistent over variations of height and buffer concentration in altered gels, see Supp. Fig. 2c [47].</figcaption></figure>

the network homogenizes in resistances with time. The smallest resistance in the network is initially ≈ 10 times smaller than the largest one; after one pulse of eroding enzyme, the smallest and the largest resistances are within a factor ≈ 3, see Fig. 1d. The normalized standard deviation $\sigma(t)/\sigma(0)$ of the normalized resistances $R_i(t)/R_4(t)$ decreases from 1 to 0.8 thanks to erosion, see inset of Fig. 1d. Likewise, the normalized flow rates homogenize, see Supp. Fig. 11 [47], with the normalized standard deviation halving from 1 to 0.54, see inset of Supp. Fig. 11 [47]. Thus, we find flow network adaptation by controlled erosion dynamics to homogenize both channel resistances and flow rates.

## B. Interplay of simulation and data yield quantitative prediction of erosion dynamics

To quantitatively understand the erosion dynamics, we follow the erosion in a device reduced in network complexity to a single PEG-NB channel, see Movie 1 in the Supplemental Material [47]. We measure the channel width and enzyme concentration across the channel and at the hydrogel wall over three subsequent erosion enzyme pulses of duration $t_e = 5$ min, each followed by 60 min of PBS solution to flush out the enzyme, see kymograph of the channel width in Fig. 2a. In the kymograph, diffusion of the enzyme from the channel into the hydrogel, and then from the hydrogel into the channel, leads to zones of high enzyme concentration in the hydrogel with a curved boundary (see red zones in Fig. 2a). The red zones correspond to regions of medium

to high enzyme concentrations, $c > 1.5 \times 10^{-6} \text{ g mol}^{-1}$, i.e. regions where the enzyme concentration is diluted less than 20 times compared to the injected concentration $c_0 = 3 \times 10^{-5} \text{ mol L}^{-1}$. For further analysis, both the enzyme concentration *c* and the channel width *w* are averaged over the entire channel length. Upon inflow of a pulse of the MMP-1 enzyme, the channel starts to erode, and its width increases at an average rate of $dw/dt \approx 5 \text{ µm min}^{-1}$, see Supp. Fig. 2b [47]. When, after a time $t_e = 5$ min, PBS is flowed in to replace the enzyme, the enzyme concentration in the channel decays to zero. During this time, enzyme molecules that had diffused into the hydrogel walls now diffuse back into the main channel and are advected away by the flow of PBS. The hydrogel walls keep eroding, albeit at a slower and slower rate, see Fig. 2a,b and Supp. Fig. 2b [47].

To understand how the enzyme controls the channel erosion rate $dw/dt$, we compute the wall erosion rate at all locations along the channel wall and measure the enzyme concentration locally within 50 µm of the channel wall, which we denote $c_w$ to distinguish it from the concentration within the channel. The higher the local enzyme concentration $c_w$, the larger the erosion rate, see Fig. 2c. For concentrations between $c_w = 1 \times 10^{-6} \text{ mol L}^{-1}$ and $c_w = 2 \times 10^{-5} \text{ mol L}^{-1}$, the erosion rate increases as $dw/dt \sim \sqrt{c_w}$, with typical values around $dw/dt \approx 5 \text{ µm min}^{-1}$ (see Supp. Fig. 2b [47]). We derive the erosion rates' scaling law from the enzyme's diffusion into the hydrogel, coupled to Michaelis-Menten kinetics for the enzymatic reaction. Call *A* the concentration of peptide crosslinker and $c_w$ the enzyme concentration. The erosion rate of the cross-linkers is given by
$dA/dt = -k_{cat}c_wA/(K_M + A)$, where $k_{cat}$ is the kinetic erosion constant in s⁻¹ and $K_M$ is the Michaelis constant in mol L⁻¹. In our experiments, the initial concentration of crosslinker $A_0 = 20$ mmol L⁻¹ is much larger than typical values reported for the Michaelis constant of the hydrolysis of similar peptide chains by MMP-1, $K_M \approx 0.5$ mmol L⁻¹ [50, 51]. Therefore, our characteristic time scale of the reaction is $T \sim A_0/(k_{cat}c_w)$. During this time, the enzyme has diffused a distance $Y \sim \sqrt{DT}$ into the hydrogel, with $D$ the diffusion coefficient of the enzyme in the hydrogel. The erosion rate $dw/dt$ is then proportional to

$$ \frac{dw}{dt} \sim \frac{Y}{T} \sim \sqrt{\frac{Dk_{cat}c_w}{A_0}} = \chi_0 \sqrt{c_w}, \quad (2) $$

recovering the square-root relationship. Note that the only unknown in Eq. (2) is the erosion constant $k_{cat}$. Fitting the scaling law to our data, we extract $k_{cat} \approx (2.1 \pm 0.7) \text{ s}^{-1}$, respectively $\chi_0 = (3.3 \pm 0.6) \times 10^3 \text{ µm min}^{-1} \text{ mol}^{-1/2} \text{ L}^{1/2}$, in good agreement with the literature [50–52], see Supp. Mat. 3 [47]. The model is consistent over variations of height and buffer concentration in altered gels, see Supp. Fig. 2c [47].

With a quantitative description of the erosion dynamics as a function of the enzyme concentration at hand, we next turn toward the dynamics of the enzyme within a channel of constant width $w$. We approximate the flow profile as parabolic along the width $w$. The cross-sectional flow velocity of the Poiseuille flow $\bar{U}$ under these assumptions is $\bar{U} = Q/wH$, where $Q$ denotes the flow rate. To optimize the numerical integration of enzyme dynamics, we reduce the dimension of the problem to the cross-sectional average concentration $\bar{c}(x, t)$ evolving only along the channel’s longitudinal axis $x$, commonly known as Taylor dispersion [53, 54]. We describe the enzyme diffusion into and out of the hydrogel with enzyme absorption and desorption kinetics into the wall, thereby neglecting the spatial variation in enzyme concentration within the channel wall, which is a good assumption for surface erosion-dominated kinetics of hydrogels [55]. For a solute getting absorbed in channel walls with rate $K_a$ and desorbed back into the channel stream with rate $K_d$, the series expansion of Zhang et al. [56] captures the cross-sectional average solute $\bar{c}(x, t)$ and absorbed solute dynamics $c_w(x, t)$ not limited by Péclet number,

$$ \frac{\partial \bar{c}}{\partial t} = -\frac{\alpha K_a}{w} \bar{c} + \frac{\alpha K_d}{w} C_w - U^{\text{eff}} \frac{\partial \bar{c}}{\partial x} + D^{\text{eff}} \frac{\partial^2 \bar{c}}{\partial x^2}, \quad (3a) $$

$$ \frac{\partial C_w}{\partial t} = K_a \bar{c} - K_d C_w, \quad (3b) $$

where we expand up to second order for the effective flow velocity $U^{\text{eff}}(t)$ and effective diffusivity $D^{\text{eff}}(t)$, which are functions of $K_a$, $K_d$, the enzyme’s molecular diffusivity $D$, and the cross-sectional flow velocity $\bar{U}$. For the absorption kinetics, we remain at first order to extract the surface concentration, $C_w = \bar{l}c_w$ directly with

$\bar{l} = 1$ mm, and $\alpha = 10^{-3}$ a non-dimensional parameter encoding the wall thickness, see Supp. Mat. 10 [47]. The advection-diffusion-sorption dynamics are therefore specified by three dimensionless parameters - the Péclet number defined as $Pe = \bar{U}w/D$, and two Damköhler numbers defined as $k_a = K_a w/D$ and $k_d = K_d w^2/D$. Note, that due to the desorption dynamics, effective flow velocity and effective diffusivity are time dependent variables in contrast to the static dimensionless parameters. We use $D = 30 \text{ µm}^2 \text{ s}^{-1}$ for the enzyme’s molecular diffusivity, based on experimental measurements (see Supp. Fig. 1 [47]), which are in agreement with the literature [57]. The absorption rate $K_a$ and desorption rate $K_d$ are extracted by fitting the experimental channel wall concentration to Eq. (3b). The rates are found to be $K_a = 8.4 \text{ µm s}^{-1}$ and $K_d = 11.4 \text{ s}^{-1}$ (see Supp. Fig. 5 [47]). Invoking our erosion law Eq. (2), the channel width $w(x, t)$ is predicted from $c_w(x, t)$ by

$$ w(x, t) = w_0 + \chi_0 \int_0^t \sqrt{c_w(x, \bar{t})} \, d\bar{t}, \quad (4) $$

where $w_0$ is the initial width of the channel. This numerical model is now our stepping stone in solving the erosion dynamics of channel networks and understanding the mechanism driving homogenization.

### C. Homogenization by erosion is robust for pulse concentration equilibration faster than pulse duration

To model erosion in channel networks, we map experimental network designs to a network skeleton consisting of network edges connected at network nodes. Each network edge represents a channel of constant width, whose hydraulic resistance is given according to channel geometry following Eq. (1). Imposing inflow and matching outflow at network inlet and outlet, respectively, flow rates in all channels and, thus, corresponding flow velocities follow from imposing conservation of fluid volume aka Kirchhoff’s law at all network nodes [15]. We then solve the advection-diffusion-sorption dynamics, Eqs. (3), in each channel of the network (see Methods). Numerically solving for the enzyme dynamics for our initial four-channel experimental parameters of $w_1 = 540 \text{ µm}$, $w_2 = 450 \text{ µm}$, $w_3 = 300 \text{ µm}$, $w_4 = 210 \text{ µm}$, $H = 400 \text{ µm}$ and $Q = 20 \text{ µL min}^{-1}$ we obtain the cross-sectional average enzyme concentration in the channel $\bar{c}(x, t)$, the absorbed enzyme concentration $c_w(x, t)$ and subsequently the channel width $w(t)$, see Fig. 1c,d, in excellent agreement with experimental observations. In particular, simulations successfully predict our hallmark of homogenization, the decrease in the normalized standard deviation of the normalized hydraulic resistances $\sigma(t)/\sigma(0)$, quantitatively. Notably, homogenization is achieved although all channels erode almost the same amount: The change in channel width over the course of a single pulse $\Delta w_i$ in
<figure><img src="image_4.png" /><figcaption>Figure 3. High absorption, low desorption and reducing Péclet number improve homogenization. (a) Normalized enzyme concentration averaged over the channel length $\langle c(x, t) \rangle$ for two parallel channels shows a slower rise in the narrow channel (green) relative to the wide channel (purple). Channels are of equal length $L = 4.5$ mm. The time scale to reach the upstream equilibrium concentration $t_{eq}$ is extracted as the time-point of $\langle \bar{c}(x, t) \rangle$ crossing 0.7. (b) Low values of $t_{eq}$ in the narrow channel, $t_{eq}^{\text{nar}}$, correspond to better homogenization - the shorter it takes for the concentration in the narrow channel to reach the upstream concentration, the more the channel erodes, leading to better homogenization of the hydraulic resistances reflected by a greater decrease in $\sigma(t = 60 \text{ s})/\sigma(0)$. Numerically determined $t_{eq}$ (blue dots) obtained from (a) for different values of $w_{\text{nar}}$ with a constant value for $w_{\text{wide}} = 300 \text{ µm}$ match the corresponding analytically calculated $t_{eq}$ from Eq. (6) (crosses). The simulation parameters are: diffusivity $D = 30 \text{ µm}^2 \text{ s}^{-1}$, absorption rate $K_a = 20 \text{ µm s}^{-1}$ and desorption rate $K_d = 30 \text{ s}^{-1}$. (c) Homogenization $\sigma(t = 60 \text{ s})/\sigma(0)$ (color coded) by erosion is best for high absorption, low desorption. Also, reducing the Péclet number increases homogenization for $Pe > 1$ considered here. Here, we only show the plane with $k_d = K_d w_{\text{nar}}^2 / D = 4$ for clarity; see Supp. Fig. 13a [47] and Movie 6 in the Supplementary Material for the full 3D phase plot [47]. $D$, $K_a$ and $K_d$ varied for $w_{\text{nar}} = 200 \text{ µm}$ and $w_{\text{wide}} = 300 \text{ µm}$. For all plots, homogenization dynamics employ a channel height $H = 1$ mm, enzyme concentration $c_0 = 3 \times 10^{-5} \text{ mol L}^{-1}$, erosion constant $\chi_0 = (3.3 \pm 0.6) \times 10^3 \text{ µm min}^{-1} \text{ mol}^{-1/2} \text{ L}^{1/2}$, and inflow rate $Q = 6 \text{ µL min}^{-1}$. Note that the pulse duration $t_e$ was set to $t_e = 1$ min for numerical efficiency. Increasing the pulse length increases the erosion, which results in a steeper drop in the homogenization metric $\sigma(t)/\sigma(0)$, indicating better homogenization, see Supp. Fig. 13c [47].</figcaption></figure>

the imbalanced four channels is within 6% for all channels.

To investigate analytically how the amount of wall erosion affects homogenization, we consider the resistance ratio of, for simplicity, two parallel channels of equal height and length but one narrow and one wide in width. As the hydraulic resistance $R$ of a channel depends nonlinearly on the channel width: $R \sim w^{-3}$, see Eq. (1), an initially small resistance ratio of a wide channel of low resistance $R_{\text{wide}}$ to a narrow channel of high resistance $R_{\text{nar}}$ is bound to increase if both channels erode equally $\Delta w_{\text{wide}} = \Delta w_{\text{nar}}$,

$$ \frac{R_{\text{wide}}(w_{\text{wide}} + \Delta w)}{R_{\text{nar}}(w_{\text{nar}} + \Delta w)} = \frac{R_{\text{wide}}(w_{\text{wide}})}{R_{\text{nar}}(w_{\text{nar}})} \left(1 + \frac{\Delta w_{\text{nar}}}{w_{\text{nar}}}\right)^3 \left(1 + \frac{\Delta w_{\text{wide}}}{w_{\text{wide}}}\right)^3 \quad (5) $$

In fact, homogenization will always take place as long as the fraction of the two cubic terms in Eq. (5) is bigger than one. Here, we remind ourselves that wall erosion is determined by Michaelis-Menten kinetics, captured by Eq. (4). Thus, homogenization by erosion is tied to the variation of enzyme absorbed into the channel walls $c_w$ and, consequently, enzyme distributed within channels across a channel network.

Therefore, we next address how the dynamics of enzyme distribution over a junction within a two-channel network affect homogenization. Mapping out the evolution of the normalized channel-averaged enzyme concentration $\langle \bar{c}(x, t) \rangle = \int_0^L \bar{c}(x, t) dx / (C_0 L)$ over the course of a 1 min pulse in two parallel channels disparate in width shows that both channels eventually reach the upstream enzyme concentration. Yet, the rise in enzyme concentration is much slower in the narrower channel compared to the wider channel, see Fig. 3a. We assess the impact of these disparate enzyme dynamics on homogenization by scoring the normalized standard deviation of the normalized hydraulic resistance $\sigma(t)/\sigma(0)$, while systematically varying the width of the narrower channel but keeping the width of the wide channel constant, see Fig. 3b. Over time $\sigma(t)/\sigma(0)$ decreases from 1 when the resistances homogenize. Thus, a lower final value indicates better homogenization. We observe that homogenization is less and less successful if the narrower channel is initially considerably narrower than the wide channel. We hypothesize that a narrower channel results in a longer time to equilibrate the inflowing enzyme concentration to the upstream enzyme concentration, such that there is less erosion of the narrower channel compared to the wide channel and, therefore, less homogenization. To test our hypothesis, we quantify the timescale to equilibrate to the upstream enzyme concentration $t_{eq}$ as the time at
<figure><img src="image_5.png" /><figcaption>Figure 4. Transition from heterogenization to homogenization at increasing erosion pulse length. (a) Two parallel channel design where parallel channels are now partitioned into three sections of varying width. The disparity in width in the entry zone of the parallel channels $w'_{nar} \ll w_{wide}$ creates a much longer time to equilibrate in the narrower channel than in the wider channel. The small difference in width in the center section $w_{nar} < w_{wide}$ reduces the required heterogeneity in channel width erosion for heterogenization, $\Delta w_{nar}/\Delta w_{wide} < w_{nar}/w_{wide}$. (b) The normalized standard deviation $\sigma(t)/\sigma(0)$ of the normalized hydraulic resistances $R_{nar}(t)/R_{wide}(t)$ of the center sections increases as the erosion pulse passing is shorter than 10 s but decreases at longer pulse length, much larger than the time to equilibrate in the narrow parallel channel $t_{eq} = 0.66$ s. The simulation parameters are: inflow rate $Q = 60 \mu L \min^{-1}$, diffusivity $D = 3 \mu m^2 s^{-1}$, absorption rate $K_a = 20 \mu m s^{-1}$ and desorption rate $K_d = 10^3 s^{-1}$ and channel height $H = 200 \mu m$.</figcaption></figure>

which the concentration reached 70% of the upstream enzyme concentration, see Fig. 3a. Indeed, long equilibration timescales accurately predict low homogenization success, see Fig. 3b. Mechanistically, the equilibration timescale is set by the time it takes for the influx of enzyme $J$ across the channel cross-section $w \cdot H$ to reach the enzyme concentration upstream of both channels $C_0$ within the entire channel volume considering a channel length $L$, i.e. $t_{eq} = C_0 L w H / J$. Neglecting diffusive influx at $Pe > 1$ considered here, the enzyme influx is given by $J = C_0 U^{\text{eff}} w H$. Therefore, the equilibration time results in

$$t_{eq} = \frac{L}{U^{\text{eff}}}, \qquad (6)$$

which is in agreement with the numerically determined equilibration timescales, see Fig. 3b. Note, that the closed expression in Eq. (6) hides the fact that the effective enzyme velocity $U^{\text{eff}}$ depends on all three non-dimensional parameters of the diffusion-advection-sorption dynamics, i.e. the Péclet number, $Pe = \bar{U}_{nar} w_{nar} / D$, absorption Damköhler number, $k_a = K_a w_{nar} / D$, and desorption Damköhler number, $k_d = K_d w_{nar}^2 / D$. To illustrate their role, we numerically assess homogenization by $\sigma(t)/\sigma(0)$ in two parallel channels for different sets of the three non-dimensional parameters,

see Fig. 3c. The channel widths are now fixed to $w_{nar} = 200 \mu m$ and $w_{wide} = 300 \mu m$, respectively, while varying kinetic constants $K_a$, $K_d$ and diffusivity $D$. We find that high absorption Damköhler number, $k_a \gg 1$, and low desorption Damköhler number, $k_d < 10$, increase homogenization. Also, reducing the Péclet number $Pe$ increases homogenization. The plane with $k_d = 4$ is depicted in Fig. 3c, see Movie 6 in Supplementary Material and Supp. Fig. 13a for the full phase space [47].

Again, homogenization is correctly predicted by the time to equilibrate $t_{eq}$ which is inversely proportional to the effective flow velocity, see Supp. Fig. 13b. Mechanistically, the effective flow velocity $U^{\text{eff}}$ is higher at high absorption kinetics, as high absorption reduces the number of slow enzymes close to the channel wall, thus effectively increasing the flow velocity of the cross-sectionally averaged enzyme concentration. Analogously, low desorption kinetics reduces the inflow of slow enzymes from the wall to the slow streamlines close to the channel wall, thus increasing the effective flow velocity of the cross-sectionally averaged enzyme concentration. We also observe that a reduced Péclet number, i.e. slower flow, increases homogenization, which is counterintuitive at first sight. However, this can be rationalized by noting that a smaller Péclet number means a smaller effective diffusivity $D^{\text{eff}}$. Reducing the Péclet number thus reduces the spread of the pulse front due to diffusion in the channels upstream; this, in turn, steepens the incoming pulse front such that, per inflow rate, more enzyme reaches a channel at the onset of the pulse, thereby reducing the time to equilibrate. Zooming from the enzyme dynamics within a single channel out to channel networks again, we find that short times to equilibrate to upstream enzyme concentration result in similar absorption of enzyme across channels and thus homogenization by equal erosion.

Yet, it is crucial to note that homogenization is not limited to fast equilibration. For homogenization to take place, in fact, narrow channels only need to erode more than the initial ratio of narrow to wide channel width, $\Delta w_{nar} / \Delta w_{wide} > w_{nar} / w_{wide}$, according to Eq. (5). Hence, channels may homogenize even if they do not erode equally but only roughly similar, which can be easily reached by adjusting the erosion pulse length. Channel wall erosion only depends on the integral over the square root of the absorbed enzyme, see Eq. (4). Thus, similar erosion of channels within a network is reachable if large times after concentration equilibration dominate wall erosion. Such large times are reached by increasing the pulse duration. We test the transition from heterogenization at low pulse length to homogenization at increased pulse length in parallel channels, where now the parallel channels are partitioned into three section of varying width, see Fig. 4. The very narrow entry section of the narrow channel ensures a much longer equilibration time in the narrow channel compared to the wider one. The small difference in width of the center section between the narrow and the wide channel, i.e. $w_{nar} / w_{wide} \approx 1$, ensures that only a little less ero-
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Erosion homogenizes complex, looped networks, and works even better with bubbles. Erosion of a loopy, hexagonal network with an initial bimodal distribution of channel widths perfused without (a,b,c,d) and with (e,f,g,h) air bubbles into wider channels before enzyme pulse (see Movies 3 & 4 in the Supplemental Material [47]). (a,e) Time-lapses of the erosion at *t* = 0 (top), *t* = 7 min (middle), and *t* = 60 min (bottom) show the hydrogel (light green) getting eroded by the enzyme (red). With bubbles blocking wide channels, see dark channels, the enzyme does not reach and thus does not erode wide channels hydrogel walls. Scale bar: 1 mm. (b,f) Increase of the channels’ width due to the hydrogel erosion. The initially wide (resp. narrow) channels are depicted in blue (resp. red). Black lines indicate numerical simulations. (c,d,g,h) Distribution of the channels’ hydraulic resistances *R*(*t*) normalized by the average hydraulic resistance of the narrow channels *R*<sub>nar</sub>(*t*), shown at the start of the experiment (*t* = 0, c & g) and after one enzyme pulse (*t* = 60 min, d & h). The value of the normalized standard deviation *σ*(*t*)/*σ*(0) of the normalized hydraulic resistance *R*(*t*)/*R*<sub>nar</sub>(*t*) of all channels (wide and narrow) is annotated for each distribution.</figcaption></figure>

sion in the narrow channel compared to the wide channel already hinders homogenization. Increasing the erosion pulse length relative to the equilibration time indeed promotes homogenization over heterogenization. Therefore, our theoretical framework predicts robust homogenization as long as the enzyme pulse length $t_e \gg t_{eq}$ is much longer than the time for enzyme concentration for any individual channel within a network to equilibrate with the enzyme concentration upstream - even independent of the network geometry.

D. Erosion homogenizes complex, looped networks, and works even better with bubbles

We put erosion-driven homogenization to the final test by eroding a loopy, hexagonal network of imbalanced channels, with initial widths either drawn from a bimodal distribution centered on *w* ≈ 65 µm and *w* ≈ 135 µm (see Fig. 5a and Supp. Fig. 6a,b [47]), or randomly distributed between 50 and 330 µm, see Supp. Fig. 6c [47].

We keep all the experimental parameters fixed except for the geometry of the network. Notably, we always set the duration of the enzyme pulse $t_e > t_{eq}$ to exceed the concentration equilibration time. We observe that all channels erode at the same rate, i.e. all channels’ widths increase at the same rate, see Fig. 5b and Movie 3 in the Supplemental Material [47]. The histogram of channel resistances shows a clear homogenization after only one enzyme pulse with the normalized standard deviation $\sigma(t)/\sigma(0)$ of normalized resistances shrinking from 1 to 0.36, see Fig. 5c,d and Supp. Fig. 7 [47]. Numerically solving the solute spread and the erosion dynamics results in excellent agreement with experimental observations in the early time regime, see Fig. 5b. Yet, neglecting the enzyme diffusion into the hydrogel in the theoretical descriptions accounts for the deviations at late-time erosion dynamics. Calculating the flow rates using the channel widths before and after erosion shows that the flow rate distribution homogenizes in the experiment shown in Fig. 5a-d, with its standard deviation shrinking from 0.31 to 0.22 over a single pulse, see Supp. Fig. 12a [47],
in perfect quantitative agreement with the decline of the standard deviation predicted by numerical simulations of the dynamics, see Supp. Fig. 12b [47]. The simulations reveal that a similar supply to all the channels, regardless of their widths, leads to a similar amount of erosion, which leads to the homogenization of the flow rates, see Supp. Fig. 12c [47]. Altogether, experimental and numerical data on loopy, hexagonal networks confirm the robustness of flow network homogenization via erosion: channels erode at the same rate, driving the homogenization of network-wide resistance and flow rates.

Although erosion homogenizes channel resistances in complex networks, full equality in resistance is never reached as both wide and narrow channels erode simultaneously. To selectively erode narrow channels only, we inject air bubbles into the network prior to eroding. Instead of carefully plugging the capillary tubing in the channel inlet while the inlet pool is already full of liquid, we plug it while it is mostly empty. This allows the formation of air bubbles, which follow the path of least resistance [58], and thus predominantly reach the low resistance large channels where they get stuck (see Supp. Fig. 10 [47] for a comparison of channel resistance and bubble occupation). The subsequently injected enzyme is blocked by the bubbles in the large, low resistance channels, and the enzyme predominantly reaches and erodes narrow channels, see Fig. 5e and Movie 4 in the Supplemental Material [47]. The erosion of the narrow channels drives them to catch up in width with wider channels, see Fig. 5f, and the resistance distribution homogenizes more than without bubbles, see Fig. 5g,h. The enhanced homogenization stands out as the initial modality of the resistances' distribution disappears in the presence of the bubbles, whereas it is still prominent without, see Fig. 5d,h. Once narrow channels reach a width comparable to that of wide channels, bubbles move freely toward their nearest low resistance channel, thus acting as equalizers of resistance by blocking the lowest resistance channels. Once the channel widths exceed $\approx 200~\mu\text{m}$, the bubbles are flushed out when the flow rate is $20~\mu\text{L min}^{-1}$. A lower flow rate allows the bubbles to remain longer in the network, increasing their efficiency to homogenize the channels' widths, see Movie 5 in the Supplemental Material [47]. Furthermore, we find that bubbles act robustly as equalizers independent of channel width randomness (see Supp. Mat. 6 [47]), as long as the bubbles are large enough to obstruct wide channels and the flow rates small enough not to flush out the bubbles.

### III. DISCUSSION

We showed that flowing pulses of an eroding solution into a network with initially imbalanced channel widths induces network self-organization, leading to the homogenization of channel resistances and flows throughout the network. Erosion works robustly as wall erosion only weakly depends, i.e. by square root, on the exact amount

of eroding solute concentration, and thus channels typically erode by a similar amount which leads to the homogenization of resistances simply by the non-linear scaling of the resistance with channel width. The key for similar wall erosion, and thereby homogenization, is the fast equilibration of solute concentration within individual channels compared to pulse length, whose dependence on Péclet and Damköhler numbers we determine analytically. To selectively erode the narrowest channels, we show the injection of bubbles prior to the pulse of the eroding solution as a promising method.

Previous work on adaptive microfluidics used the swelling and contraction of embedded responsive hydrogels to control flow in microchannels by switching channels open or closed [38–43, 46]. Yet, no adaptation of network geometry crucial for flow homogenization was attempted. Here, we embed changes in channel geometry and their associated fluid flows due to the progressive but controlled erosion of hydrogel walls, leading to continuous homogenization of the channel resistances, showing that even a simple adaptation, once controlled by timing, can lead to network homogenization. As flow network architecture plays a key role in biological and synthetic active fluids [59, 60], the feedback between hydrogel channels and flows introduced here has the potential to broaden the design rules not only for microfluidics but active fluid networks alike.

The homogenization mechanism introduced in this work is based on the local erosion of the channel wall by an eroding chemical transported by the fluid flow. The presence of bubbles usually considered a severe problem to circumvent in fluidic applications [61–63], actually plays to the advantage of homogenization by erosion here. As the hydrogel used in our experiments is bio-compatible [64], purposely tailoring the timing of erosive chemical pulsing emerges as a tool for designing adaptive, biomimetic self-controlled fluidic networks for biotechnological applications. In particular the challenge of tissue perfusion in bio-engineered organs could be overcome by the local channel adaptation dynamics introduced here [1, 16].

Nature's tendency is to create disorder. Flow networks in this regard resemble disordered systems in general, and in their statistical properties granular systems in particular [13]. Counteracting disorder by such simple means as erosion may, therefore, be applicable to other disordered systems. Note, that, as our theoretical description of the dynamics underlines, flow homogenization by erosion is general and not limited to the specific hydrogel and enzyme used in this work. Neither does erosion homogenize solely in two dimensions, used here for quantitative understanding: the non-linear dependency of channel resistance on channel width prevails in three dimensions, as do the transport dynamics, envisioning applications to three dimensional porous media such as packed-bed reactors. Self-homogenization by pulse-controlled erosion is a generally simple and attractive solution for applications where inhomogeneities cannot be avoided to arise
over time.

ACKNOWLEDGMENTS

This project has received funding from the French Agence Nationale de la Recherche (ANR), under grant ANR-21-CE30-0044, and the Deutsche Forschungsmeinschaft (DFG, German Research Foundation) under grant AL 1429/5-1 and the INST 95/1634-1 FUGG as well as funding by the IGSSE / TUM Graduate School for IPT SOFT. This research was conducted within the Max Planck School Matter to Life supported by the German Federal Ministry of Education and Research (BMBF) in collaboration with the Max Planck Society.

Appendix A: Experimental Methods

The experimental setup consists of microfluidic channels whose side walls are made of a norbornene polyethylene glycol (PEG-NB-8arm) backbone, cross-linked with a peptide chain (KCGPQGIWGQCK-OH, Iris Biotech, Germany). In the presence of the enzyme matrix metalloproteinase-1 (MMP-1, from Clostridium histolyticum, Sigma-Aldrich, Germany), the peptide chain is cleaved and the PEG hydrogel erodes.

Briefly, to fabricate erodible microfluidic channels, 30 µL of an aqueous polymer solution, containing 5 mmol L⁻¹ PEG-NB-8arm, 20 mmol L⁻¹ Crosslinker, 3 mmol L⁻¹ of the photoinitiator lithium phenyl-2,4,6-trimethylbenzoylphosphinate (LAP) and 2 mmol L⁻¹ of a green fluorescent dye (fluorescein-PEG-Thiol from NanoCS, USA) is flowed into a 400 µm-height plastic microfluidic chip (μ-Slide VI 0.4 ibiTreat, Ibidi, Germany). The chip is aligned on a shadow mask with the desired channel geometry and illuminated through a collimated light source at 365 nm; the exposed areas are cross-linked and form the walls of the microfluidic channels [64–66]. The channels are then washed with PBS (Phosphate Buffer Saline without Calcium nor Magnesium, VWR), and stored at 4°C filled with PBS in humid conditions for at least 5 days. This results in well-defined fluorescent hydrogel walls of 400 µm height and up to 10 mm in length, see Fig. 1a.

To erode the hydrogel walls, a solution containing $c_0 = 3 \times 10^{-5}$ mol L⁻¹ MMP-1 and Texas Red-labeled dextran with the same molecular weight, 70 kDa, is injected at a flow rate of $Q = 20$ µL min⁻¹ into the microfluidic networks. Beforehand, PBS is flowed into the channels for 30 min at the same flow rate $Q$ to flush out both the green dye, that had diffused from the hydrogel into the channels during the week-long storage, and bubbles that may have appeared upstream during the tubing connection. Glass syringes of 1 and 10 mL (CETONI, Germany) are set up on two syringe pumps (Nemesys S, CETONI, Germany), which are controlled by a computer through CETONI Elements software (CE-

TONI, Germany). This software allows for chronological automation of input flows via a script system, enabling perfect repeatability of the injection protocol. Images of the erosion process are taken by fluorescence and bright-field imaging every 20 s with a Hamamatsu ORCA-Flash 4.0 digital camera under an AxioZoom V.16 microscope (Zeiss, Germany), using a Zeiss PlanNeoFluar Z 1x objective. Zeiss Zen 3.2 (blue edition) software was used for imaging.

Appendix B: Numerical Methods

To numerically solve for the spread and absorption of a diffusive solute along the channel network, a Crank-Nicolson routine was employed to integrate the dynamics of the cross-sectionally averaged diffusion-advection-sorption dynamics Eqs. 3, see [56], implemented in MATLAB (MathWorks). The channel geometry used for the simulation in Figs. 1c, d and 5b followed the same design as the hydrogel channel networks: channels were reduced to their skeleton, decomposing every channel into segments which were further decomposed to boxes which is the smallest simulation unit. The connectivity of every channel was encoded into a connectivity matrix, thus defining the geometry of the network. Next, based on the channel width and length, each channel’s hydraulic resistance was computed. Then, setting the inflow rate at the inlet and the matching outflow rate at the outlet, Kirchhoff’s circuit law was employed, conserving fluid volume at every network node, to compute the flow profile in the channel network. Lastly, the network topology and the flow profile were used to numerically integrate the advection-diffusion-sorption equation along individual channels using the Crank-Nicolson integration routine. All model parameters were directly quantified from experimental data, except the absorption/desorption rate which were matched to fit experiments. At channel junctions, merging and splitting of solute concentration conditions were implemented following [15]. The solute concentration at the inlet is set to vary in time. To this end, we integrated the Taylor dispersion of the solute pulse along the tubing between the upstream switch and the chip inlet. In general, for a rectangular concentration pulse of concentration $\bar{c}_0$ extending from $x_0$ to $x_0 - \Delta x$ in one dimension, the concentration profile is given by the result of the spatial superposition i.e. integral overall solutions for individual $\delta$-peaks,

$$ \bar{c}(x,t) = \frac{\bar{c}_0 \sqrt{\pi}}{2} \left[ \erf\left( \frac{(x - U_{\text{tube}}t - x_0 + \Delta x)}{\sqrt{4D_{\text{tube}}t}} \right) - \erf\left( \frac{(x - U_{\text{tube}}t - x_0)}{\sqrt{4D_{\text{tube}}t}} \right) \right], \quad (\text{B1}) $$

where, $r_{\text{tube}}$ is the effective tube radius, $U_{\text{tube}} = Q/\pi r_{\text{tube}}^2$ and $D_{\text{tube}} = D \left(1 + \frac{r_{\text{tube}}^2 U_{\text{tube}}^2}{48D^2}\right)$. Both the flow
rate $Q$ and the distance between the switch and the chip inlet, $x_0$, are set by the experiments. The tube radius is variable throughout the length of the tube, such that $r_{\text{tube}}$, $D_{\text{tube}}$, and $\bar{c}_0$ are determined by fitting to the data. We fit Eq. (B1) to the experimentally obtained channel concentration at the channel inlet, normalized by the enzyme concentration injected in experiments, to obtain $\bar{c}_0 = 0.72$, which is less than one due to the cross-sectional averaging, and use these values to calculate the concen-

tration at the channel inlet $x=0$ using Eq. (B1). At the channel outlet, an open outflow condition was implemented, which estimated transport by advection and diffusion across the last spatial simulation point as in [15]. The effective coefficients $U^{\text{eff}}$ and $D^{\text{eff}}$ used for simulations were validated for $Pe = 10$, $k_a = 50$ and $k_d = 1$ against analytical solutions and plots [56] within a single channel.

[1] P. P. Stankey, K. T. Kroll, A. J. Ainscough, D. S. Reynolds, A. Elamine, B. T. Fichtenkort, S. G. Uzel, and J. A. Lewis, Embedding biomimetic vascular networks via coaxial sacrificial writing into functional tissue, *Advanced Materials*, e2401528 (2024).

[2] H. Sun, J. Zhu, D. Baumann, L. Peng, Y. Xu, I. Shakir, Y. Huang, and X. Duan, Hierarchical 3d electrodes for electrochemical energy storage, *Nat. Rev. Mater.* **4**, 45 (2019).

[3] C. Tien and A. C. Payatakes, Advances in deep bed filtration, *AIChE J.* **25**, 737 (1979).

[4] W. R. Bowen and F. Jenner, Theoretical descriptions of membrane filtration of colloids and fine particles: an assessment and review, *Adv. Colloid Interfac.* **56**, 141 (1995).

[5] K. Beven and P. Germann, Macropores and water flow in soils, *Water Resour. Res.* **18**, 1311 (1982).

[6] L. Lebon, L. Oger, J. Leblond, J. P. Hulin, N. S. Martys, and L. M. Schwartz, Pulsed gradient NMR measurements and numerical simulation of flow velocity distribution in sphere packings, *Phys. Fluids* **8**, 293 (1996).

[7] Y. E. Kutsovsky, L. E. Scriven, H. T. Davis, and B. E. Hammer, NMR imaging of velocity profiles and velocity distributions in bead packs, *Phys. Fluids* **8**, 863 (1996).

[8] M. D. Shattuck, R. P. Behringer, G. A. Johnson, and J. G. Georgiadis, Convection and flow in porous media. Part 1. Visualization by magnetic resonance imaging, *J. Fluid Mech.* **332**, 215 (1997).

[9] R. S. Maier, D. M. Kroll, Y. E. Kutsovsky, H. T. Davis, and R. S. Bernard, Simulation of flow through bead packs using the lattice Boltzmann method, *Phys. Fluids* **10**, 60 (1998).

[10] A. Sederman and L. Gladden, Magnetic resonance visualisation of single- and two-phase flow in porous media, *Magn. Reson. Imaging* **19**, 339 (2001).

[11] M. Moroni and J. H. Cushman, Statistical mechanics with three-dimensional particle tracking velocimetry experiments in the study of anomalous dispersion. II. Experiments, *Phys. Fluids* **13**, 81 (2001).

[12] S. S. Datta, H. Chiang, T. S. Ramakrishnan, and D. A. Weitz, Spatial Fluctuations of Fluid Velocities in Flow through a Three-Dimensional Porous Medium, *Phys. Rev. Lett.* **111**, 064501 (2013).

[13] K. Alim, S. Parsa, D. A. Weitz, and M. P. Brenner, Local pore size correlations determine flow distributions in porous media, *Phys. Rev. Lett.* **119**, 144501 (2017).

[14] M. Matyka, J. Gołembiewski, and Z. Koza, Power-exponential velocity distributions in disordered porous media, *Phys. Rev. E* **93**, 013110 (2016).

[15] F. J. Meigel, T. Darwent, L. Bastin, L. Goehring, and K. Alim, Dispersive transport dynamics in porous media emerge from local correlations, *Nat. Commun.* **13**, 5885 (2022).

[16] Z. A. Sexton, D. Rütsche, J. E. Herrmann, A. R. Hudson, S. Sinha, J. Du, D. J. Shiwarski, A. Masaltseva, F. S. Solberg, J. Pham, J. M. Szafron, S. M. Wu, A. W. Feinberg, M. A. Skylar-Scott, and A. L. Marsden, Rapid model-guided design of organ-scale synthetic vasculature for biomanufacturing, *Science* (New York, N.Y.) **388**, 1198 (2025).

[17] S. Isogai, M. Horiguchi, and B. M. Weinstein, The vascular anatomy of the developing zebrafish: an atlas of embryonic and early larval development, *Dev. Biol.* **230**, 278 (2001).

[18] L. Boddy, J. Hynes, D. P. Bebber, and M. D. Fricker, Saptrotrophic cord systems: dispersal mechanisms in space and time, *Mycoscience* **50**, 9 (2009).

[19] A. Tero, S. Takagi, T. Saigusa, K. Ito, D. P. Bebber, M. D. Fricker, K. Yumiki, R. Kobayashi, and T. Nakagaki, Rules for biologically inspired adaptive network design, *Science* **327**, 439 (2010).

[20] S.-S. Chang and M. Roper, Microvascular networks with uniform flow, *J. Theor. Biol.* **462**, 48 (2019).

[21] S.-S. Chang, S. Tu, K. I. Baek, A. Pietersen, Y.-H. Liu, V. M. Savage, S.-P. L. Hwang, T. K. Hsiai, and M. Roper, Optimal occlusion uniformly partitions red blood cells fluxes within a microvascular network, *PLoS Comput. Biol.* **13**, e1005892 (2017).

[22] J. B. Kirkegaard and K. Sneppen, Optimal Transport Flows for Distributed Production Networks, *Phys. Rev. Lett.* **124**, 208101 (2020).

[23] A. Simonin, J. Palma-Guerrero, M. Fricker, and N. L. Glass, Physiological significance of network organization in fungi, *Eukaryot. Cell* **11**, 1345 (2012).

[24] K. Alim, N. Andrew, A. Pringle, and M. P. Brenner, Mechanism of signal propagation in physarum polycephalum, *Proc. Natl. Acad. Sci. U. S. A.* **114**, 5136 (2017).

[25] S. Liese, L. Mahadevan, and A. Carlson, Balancing efficiency and homogeneity of biomaterial transport in networks, *Europhys. Lett.* **135**, 58001 (2021).

[26] F. Kramer and C. D. Modes, Biological flow networks: Antagonism between hydrodynamic and metabolic stimuli as driver of topological transitions, *Phys. Rev. Research* **5**, 023106 (2023).

[27] J. Karschau, A. Scholich, J. Wise, H. Morales-Naverrete, M. Zerial, and B. M. Friedrich, Resilience of three-dimensional sinusoidal networks in liver tissue, *PLoS
Comput. Biol. **16**, e1007965 (2020).

[28] F. J. Meigel and K. Alim, Flow rate of transport network controls uniform metabolite supply to tissue, *J. Roy. Soc. Interface* **15**, 20180075 (2018).

[29] Y. Qi and M. Roper, Control of low flow regions in the cortical vasculature determines optimal arterio-venous ratios, *Proc. Natl. Acad. Sci. U. S. A.* **118**, e2021840118 (2021).

[30] C. D. Murray, The Physiological Principle of Minimum Work: I. The Vascular System and the Cost of Blood Volume, *Proc. Natl. Acad. Sci. U. S. A.* **12**, 207 (1926).

[31] S. Marbach, N. Ziethen, L. Bastin, F. K. Bäuerle, and K. Alim, Vein fate determined by flow-based but time-delayed integration of network architecture, *Elife* **12**, e78100 (2023).

[32] A. Zareei, D. Pan, and A. Amir, Temporal Evolution of Erosion in Pore Networks: From Homogenization to Instability, *Phys. Rev. Lett.* **128**, 234501 (2022).

[33] T. Gavrilchenko and E. Katifori, Distribution networks achieve uniform perfusion through geometric self-organization, *Phys. Rev. Lett.* **127**, 078101 (2021).

[34] M. Kramar and K. Alim, Encoding memory in tube diameter hierarchy of living flow network, *Proc. Natl. Acad. Sci. U. S. A.* **118**, e2007815118 (2021).

[35] K. H. Wong, J. M. Chan, R. D. Kamm, and J. Tien, Microfluidic models of vascular functions, *Annu. Rev. Biomed. Eng.* **14**, 205 (2012).

[36] B. Sebastian and P. S. Dittrich, Microfluidics to mimic blood flow in health and disease, *Annu. Rev. Fluid Mech.* **50**, 483 (2018).

[37] M. Fenech, V. Girod, V. Claveria, S. Meance, M. Abkarian, and B. Charlot, Microfluidic blood vasculature replicas using backside lithography, *Lab Chip* **19**, 2096 (2019).

[38] D. J. Beebe, J. S. Moore, J. M. Bauer, Q. Yu, R. H. Liu, C. Devadoss, and B.-H. Jo, Functional hydrogel structures for autonomous flow control inside microfluidic channels, *Nature* **404**, 588 (2000).

[39] D. T. Eddington and D. J. Beebe, Flow control with hydrogels, *Adv. Drug Deliver. Rev.* **56**, 199 (2004).

[40] A. Richter, D. Kuckling, S. Howitz, T. Gehring, and K.-F. Arndt, Electronically controllable microvalves based on smart hydrogels: magnitudes and potential applications, *J. Microelectromech.* **S. 12**, 748 (2003).

[41] L. D'eramo, B. Chollet, M. Leman, E. Martwong, M. Li, H. Geisler, J. Dupire, M. Kerdraon, C. Vergne, F. Monti, et al., Microfluidic actuators based on temperature-responsive hydrogels, *Microsyst. & Nanoeng.* **4**, 1 (2018).

[42] A. Baldi, Y. Gu, P. Loftness, R. Siegel, and B. Ziaie, A hydrogel-actuated smart microvalve with a porous diffusion barrier back-plate for active flow control, in *Technical Digest. MEMS 2002 IEEE International Conference. Fifteenth IEEE International Conference on Micro Electro Mechanical Systems (Cat. No.02CH37266)*, MEMSYS-02, IEEE (IEEE, 2002) pp. 105-108.

[43] R. H. Liu, Q. Yu, and D. J. Beebe, Fabrication and characterization of hydrogel-based microvalves, *J. Microelectromech.* **S. 11**, 45 (2002).

[44] C. B. Goy, R. E. Chaile, and R. E. Madrid, Microfluidics and hydrogel: A powerful combination, *React. Funct. Polym.* **145**, 104314 (2019).

[45] F. Paratore, V. Bacheva, M. Bercovici, and G. V. Kaigala, Reconfigurable microfluidics, *Nat. Rev. Chem.* **6**, 70 (2022).

[46] H. Na, Y.-W. Kang, C. S. Park, S. Jung, H.-Y. Kim, and J.-Y. Sun, Hydrogel-based strong and fast actuators by electroosmotic turgor pressure, *Science* **376**, 301 (2022).

[47] See Supplemental Material at [URL will be inserted by publisher] for details on enzyme dynamics, more network architectures self-organizing, flow distributions of self-organizing networks, full parameters sweeps for homogenization, details on the simulation and derivations as well as movies of experimental observations. The experimental data, simulation code, and movies are available in the mediaTUM repository at doi:10.14459/2024mp1752696. The Supplemental Material also includes Ref. [48].

[48] L. M. Weber, C. G. Lopez, and K. S. Anseth, Effects of peg hydrogel crosslinking density on protein diffusion and encapsulated islet survival and function, *J. Biomed. Mater. Res. A* **90**, 720 (2009).

[49] H. Bruus, Acoustofluidics 1: Governing equations in microfluidics, *Lab Chip* **11**, 3742 (2011).

[50] M. P. Lutolf, J. L. Lauer-Fields, H. G. Schmoekel, A. T. Metters, F. E. Weber, G. B. Fields, and J. A. Hubbell, Synthetic matrix metalloproteinase-sensitive hydrogels for the conduction of tissue regeneration: engineering cell-invasion characteristics, *Proc. Natl. Acad. Sci. U. S. A.* **100**, 5413 (2003).

[51] H. Nagase and G. B. Fields, Human matrix metalloproteinase specificity studies using collagen sequence-based synthetic peptides, *Pept. Sci.* **40**, 399 (1996).

[52] J. Patterson and J. A. Hubbell, Enhanced proteolytic degradation of molecularly engineered peg hydrogels in response to mmp-1 and mmp-2, *Biomaterials* **31**, 7836 (2010).

[53] G. I. Taylor, Dispersion of soluble matter in solvent flowing slowly through a tube, *Proc. R. Soc. London, A* **219**, 186 (1953).

[54] R. Aris, On the dispersion of a solute in a fluid flowing through a tube, *Proc. R. Soc. London, A* **235**, 67 (1956).

[55] P. M. Kharkar, K. L. Kiick, and A. M. Kloxin, Designing degradable hydrogels for orthogonal control of cell microenvironments, *Chem. Soc. Rev.* **42**, 7335 (2013).

[56] L. Zhang, M. A. Hesse, and M. Wang, Transient solute transport with sorption in Poiseuille flow, *J. Fluid Mech.* **828**, 295 (2017), arXiv:1909.08720 [physics].

[57] P. B. Gaspers, C. R. Robertson, and A. P. Gast, Enzymes on immobilized substrate surfaces: diffusion, *Langmuir* **10**, 2699 (1994).

[58] H. Boukellal, Š. Selimović, Y. Jia, G. Cristobal, and S. Fraden, Simple, robust storage of drops and fluids in a microfluidic device, *Lab Chip* **9**, 331 (2009).

[59] C. Jorge, A. Chardac, A. Poncet, and D. Bartolo, Active hydraulics laws from frustration principles, *Nature Physics* **20**, 303 (2024).

[60] F. G. Woodhouse, A. Forrow, J. B. Fawcett, and J. Dunkel, Stochastic cycle selection in active flow networks, *Proceedings of the National Academy of Sciences* **113**, 8200 (2016), 1607.08015.

[61] P. Gravesen, J. Branebjerg, and O. S. Jensen, Microfluidics-a review, *J. Micromech. Microeng.* **3**, 168 (1993).

[62] J. H. Sung and M. L. Shuler, Prevention of air bubble formation in a microfluidic perfusion cell culture system using a microscale bubble trap, *Biomed. Microdevices* **11**, 731 (2009).

[63] C. Lochovsky, S. Yasotharan, and A. Günther, Bubbles no more: in-plane trapping and removal of bubbles in
microfluidic devices, Lab Chip **12**, 595 (2012).

[64] S. Stöberl, M. Balles, T. Kellerer, and J. O. Rädler, Photolithographic microfabrication of hydrogel clefts for cell invasion studies, Lab Chip **23**, 1886 (2023).

[65] J. G. Skinner, T. R. Groves, A. Novembre, H. Pfeiffer, and R. Singh, Photomask fabrication procedures and limitations,

in *Handbook of Microlithography, Micromachining, and Microfabrication. Volume 1: Microlithography*, Vol. 1 (SPIE PRESS, 1997) pp. 377–474.

[66] M. Dietrich, H. Le Roy, D. B. Brückner, H. Engelke, R. Zantl, J. O. Rädler, and C. P. Broedersz, Guiding 3d cell migration in deformed synthetic hydrogel microstructures, *Soft Matter* **14**, 2816 (2018).