Flow rate of transport network controls uniform metabolite supply to tissue

Consider a cylindrical tube filled with fluid flowing at flow velocity $U(r,z)$ along the tube. Metabolites are advected with the flow and in addition disperse due to molecular diffusion with diffusivity $\kappa$. Considering the small scales of xylem vessels and xylem flow, we are in the regime of low Reynolds number, where flow is best described by laminar Poiseuille flow. Thus, we describe the flow along the longitudinal axis $z$ varying in radial direction $r$ in the circular tube with a radius $R$ and length $\ell$ as, $U(r,z)=2\left(1-\frac{r^{2}}{R^{2}}\right)\langle U(z)\rangle_{r}$. Here, $\langle U\rangle_{r}$ denotes the cross-sectionally averaged longitudinal velocity. The spread of metabolites of concentration $C$ is thus fully described by

Metabolites are absorbed into the surrounding tissue along the tube wall, given by the boundary condition

analogous to heat absorption or surface reactions. Here, the parameter $\nu$ denotes the metabolite absorption rate at the tube wall. Dividing $\nu$ by the molecular diffusivity we define the absorption parameter $\gamma$, where we consider $\nu$ as a constant tissue property. According to e._m._lungu_effect_1982 the advection-diffusion equation, Eq. (1), can be re-formulated employing both the boundary condition Eq. (2) and the Poiseuille profile as a single absorption-advection-diffusion equation in cylindrical coordinates.

In analogy to the derivation of Taylor Dispersion by G.I. Taylor taylor_dispersion_1953 , a simpler, though approximated, expression is possible where the concentration dynamics only depend on the longitudinal coordinate $z$. To this end the metabolite concentration is separated into the sum of a cross-sectional average concentration $\langle C\rangle_{r}$ and the radial variation $C^{\prime}$, $C(r,z)~=~\langle~C(z)~\rangle_{r}~+~C^{\prime}(r,z)$. The multidimensional diffusion-advection for $C(r,z)~=~\langle~C(z)~\rangle_{r}~+~C^{\prime}(r,z)$ can be simplified to an equation for the cross-sectionally averaged concentration $\langle~C(z)~\rangle_{r}$ if the cross-sectional variations of the concentration $C^{\prime}(r,z)$ are much smaller than the averaged concentration itself taylor_dispersion_1953 ; Cha_inprep , resulting in

This approach employs three approximations. First, the time-scale of diffusion across the tube’s cross-section has to be much smaller than the time-scale of advection, $\frac{\ell}{\langle U\rangle_{r}}\gg\frac{R^{2}}{\kappa}$. This sets an upper bond for the later choice of the fluid inflow rate. The second assumption states that the cross-sectional variations of the concentration have to be small $\langle C\rangle_{r}\gg C^{\prime}$. Since a high absorption parameter $\gamma$ implies a large concentration gradient across the cross-section, the second assumption implies $\gamma R\ll 1$. Third, the variation of $C^{\prime}$ has to be much greater in radial direction than in flow direction $\partial^{2}_{r}C^{\prime}\gg\partial^{2}_{z}C^{\prime}$. The last assumption implies the tube radius to be smaller than its length $R\ll\ell$. This is fulfilled by long slender tubes.

Employing these assumptions, the cross-sectional average metabolite concentration along a tube in steady state is given by an exponential decay from initial concentration $C_{0}$,

Here, we introduced two dimensionless variables Pe and $S$. $\text{Pe}=\frac{\langle U\rangle_{r}\ell}{\kappa}$ is the well known Peclet number describing the relation between diffusive and advective time scale. $S$ is the ratio of the time scale for absorption, given by the product of dimensionless absorption parameter and time to diffuse across the tube’s cross-section, and the time to be advected out of the tube, resulting in $S=\frac{\gamma\kappa\ell}{R\langle U\rangle_{r}}$.

Considering a constant influx $J_{0}$ by advection and diffusion of metabolites at the tube’s start, we find that the initial concentration $C_{0}$ is given by

The overall absorption $\phi$ along a tube is given by the integrated flux of metabolites across the tube wall $\mathcal{W}$, $\phi=2\pi R\int_{\mathcal{W}}\kappa\nabla C\text{d}z$, where $\text{d}z$ is integrating over the length of the tube. As in the derivation of the effective diffusion-advection-absorption equation Eq. (3), we use $R\gamma\ll 1$ to arrive at

where $\Lambda$ is an abbreviation $\Lambda=\sqrt{8S/\text{Pe}+3/4\gamma R+1}$. We identify two factors that control the absorption in a tube. The first is the total influx of metabolites over the cross-sectional area of the tube $\pi R^{2}J_{0}$. The total influx of metabolites is the upper limit for absorption in the tube. The second factor is the tube’s capacity to absorb metabolites as $\hat{\phi}=\phi/\pi R^{2}J_{0}$ with $\hat{\phi}\in[0,1]$. This absorption capacity is independent of the concentration of metabolites and only depends on the parameters of the tube and the flow velocity within the tube.

For the derivation of the optimal inflow rate, it is essential to approximate the absorption capacity as resulting from Eq. (6) above. We initially approximate the inverse of the absorption capacity $\hat{\phi}^{-1}$ by taking a finite $\text{Pe}>0$ and using $R\gamma\ll 1$ to find $\hat{\phi}^{-1}=\frac{1}{2S}+1$. Resubstituting the system’s parameters for $S$ we find for the absorption capacity of a tube

Note the simple dependence of the absorption capacity on the cross-sectionally average flow velocity in the tube. The approximation of the absorption capacity has been verified numerically to hold over the parameter space considered here, see supplemental information S4. Note, that this simplified expression is only used for the analytical derivation of the optimal inflow rate. For simulations the full expression Eq. 6 is used. From now on, we drop the brackets $\langle\rangle_{r}$ and only refer to cross-sectional averaged observables.

In a transport network, individual tubes are connected at nodes. Here, we aim to model the geometry of higher order xylem veins branching from a second order vein in dicotyledons, as shown in Fig. 1 b). We choose a planar transport network representing a rectangular excerpt of the leaf tissue. For the tissue excerpt we choose a general vascularisation using a slightly randomized tessellation of space, where the network is build with small triangles, known to bear least artifacts katifori_damage_2010 . A small Gaussian noise of a twentieth tube length $\ell$ is added to the positions of the tessellation nodes to avoid pattern artifacts arising from the underlying topology otherwise, see supplementary material S2. The tube length varies accordingly in a normal distribution around the mean tube length $\langle\ell\rangle$. In agreement with observations of diminishing hierarchy in the higher interwebbed xylem vessel radii ronellenfitsch_topological_2015 ; sack_leaf_2013 , we set the same radius $R$ for all tubes. Fluid and metabolites are flowing into the network at network nodes along one side of the rectangular region, representing the supply from secondary veins into the tissue. Xylem vessels are organized in vascular bundles in lower order veins that branch out into the interwebbed higher order xylem network thorne_structure_2006 presumingly supplying the same flow at every inflow node. Therefore, we approximate inflow rates $Q_{\text{in}}$ to be equal at all inflow nodes. To represent the effect of fluid evaporation at stomata, fluid, but not metabolite, is flowing out at every node in the network $Q_{\text{out}}$, see Fig. 1 c). Metabolites are absorbed across each tube wall. The absorption rate $\nu$ is constant throughout the network. Yet, as we have learned by studying a single tube in the previous section, the amount of metabolite absorbed depends on how much metabolite available in the fluid, and how much time the metabolite has to travel to the tube wall to get absorbed. Therefore, absorption despite a constant absorption rate varies largely within a network.

The flow of the metabolites is determined by the fluid flow in accordance with Eq. 1. The fluid flow throughout a network is fully defined by the network’s architecture, the inflow and outflow rates and Kirchhoff’s circuit law. The cross-sectionally averaged fluid flow velocity in a tube follows subsequently from pressure drops $\Delta P=R_{\text{hyd}}\cdot Q$ along the tube. Each tube is considered as straight cylinder with hydraulic conductance of $K_{\text{hyd}}^{-1}=\frac{8}{\pi}\ell\eta\frac{1}{R^{4}}$, where $\eta$ denotes the dynamic fluid viscosity. The pressure at every node is calculated by multiplying the inverse of the network conductivity matrix with the inflow or outflow rates at every node. The pressure drop along the tube is the difference between the pressure values at start and end node. The fluid flow is solved consistently through the whole network and takes network geometry and viscous and friction forces via the hydraulic resistance into account. Considering steady state solutions, the flow does not fluctuate over time.

The absorption of metabolite across a tube’s wall within the network depends on the metabolite available and the tube’s absorption capacity. The absorption capacity of tube follows directly from each tubes physical parameters and the fluid flow velocity within the tube. Next, we need to calculate the influx of metabolites $J_{0}$ in every tube which we solve for iteratively throughout the network starting with the influx nodes, see supplementary material S3A. As simplification we focus on stationary, steady state absorption patterns. We employ that metabolite flux is conserved at every network node. All metabolites flowing into an node are redistributed into tubes originating from this node. Redistribution is proportional to diffusion and flux into each tube. The metabolite outflux at the end of a tube is then given by the difference of metabolite influx and total absorption along the tube. Finally, at the lower end of the considered network excerpt, opposite the inflow nodes, remaining metabolites are flowing out of the network. Since the outflowing metabolites would lead to an accumulation of metabolites, we state the amount of metabolites not absorbed for every considered network excerpt.

Taking our initial motivation from plant leaves, we choose an average tube length $\ell=$0.1\text{\,}\mathrm{mm}$$ and tube radius $R=$3\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ in accordance with xylem vessels sack_leaf_2013 ; sack_hydraulic_2004 ; sack_developmentally_2012 . Note that there is a difference between leaf vein and xylem vessel radius, as leaf veins bundle both phloem and xylem vessels. We thus consider xylem vessels to be less than half the radius of leaf veins in our parameter choice. The order of magnitude of the total inflow rate is chosen to yield velocities observable in lower order xylem vessels $\langle U\rangle_{r}\approx$1\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{\,}{\mathrm{s}}^{-1}$$. We vary the fluid inflow rate between $Q_{\text{in}}~=~$0.8\text{\times}{10}^{-6}\text{\,}{\mathrm{mm}}^{3}\text{\,}{\mathrm{s}}^{-1}$$ to $Q_{\text{in}}~=~$6.4\text{\times}{10}^{-6}\text{\,}{\mathrm{mm}}^{3}\text{\,}{\mathrm{s}}^{-1}$$. The choice of the inflow is consistent with a average water evaporation of approximately $0.1\text{\,}\mathrm{mol}\text{\,}{\mathrm{m}}^{-2}\text{\,}{\mathrm{s}}^{-1}$ vico_effects_2011 ; zwieniecki_hydraulic_2002 and an average stomata density of $200\text{\,}{\mathrm{mm}}^{-2}$ woodward_influence_1995 . For the molecular diffusivity, we consider small molecules with $\kappa=$1\text{\times}{10}^{-10}\text{\,}{\mathrm{m}}^{2}\text{\,}{\mathrm{s}}^{-1}$$. For the network size, we choose a triangulation with $\mathcal{N}\approx 1000$ tubes. Since xylem vessels consist of highly pitted dead lignified tissue, no active absorption by chemical reactions but passive absorption by membrane permeation is expected. Values for membrane permeation are typically in the range of $\nu\approx$1\text{\times}{10}^{-9}\text{\,}\mathrm{m}$$ phillips_physical_2013 and depend on both membrane and metabolite properties. Alternatively an estimation for the absorption parameter $\gamma$ can be derived from concentration profiles in xylem veins vanBel:1979 . Translating the measured exponential decay for higher order veins, we find $\gamma\approx$10\text{\,}{\mathrm{m}}^{-1}$$ Horwitz:1958 ; vanBel:1979 . This estimate for $\gamma$ is in accordance with estimates using the membrane permeation. Thus, for the numerical calculations an absorption parameter of $\gamma=$10\text{\,}{\mathrm{m}}^{-1}$$ is chosen. Revisiting the three assumption made in Eq. 3 we verify that these assumptions hold for the chosen network topologies, see supplementary material S4.

We first study metabolite supply patterns in uniform transport networks, where all tubes have the same radius $R$. To compare different inflow rates we normalize the absorption by the total influx of metabolites. We find that the total fluid inflow rate is dominating the supply patterns, see Fig. 2. For small inflow rate, average flow velocities in the network are slow and the highest absorption is near the inflow nodes. Metabolites are not transported through to the end of the network limiting supply there. Calculating the mean absorption per row $\phi_{\text{row}}$ from inflow to opposing end, we can characterize this regime by $\phi_{\text{row}}>\phi_{\text{row}+1}$, see Fig. 2 right column. For high inflow rates, average flow velocities are fast and the absorption increases with the distance from the inflow nodes $\phi_{\text{row}}<\phi_{\text{row}+1}$. Metabolites arrive at the end of the network too quickly before getting absorbed limiting supply close to the inflow nodes. Between these limiting cases, we identify an inflow rate that gives rise to an optimally uniform supply pattern. We define the optimum by the lowest variance. The variance is calculated over the ensemble of all tubes. The overall variance in the optimal case is $0.07$ compared to $2.35$ and $0.75$ in the examples of low and high inflow rate shown in Fig. 2. In the optimally uniform supply pattern, absorption is the same constant rate in subsequent rows, $\phi_{\text{row}}=\phi_{\text{row}+1}$. On the basis of this simple relation stating uniform absorption, a scaling law for the optimal inflow rate is derived next.

To derive a scaling law for the optimal inflow rate that gives rise to the most uniform supply pattern, we consider a one-dimensional toy model of connected single tubes that captures the essential flow and transport characteristics along the rows of the two-dimensional transport networks investigated above. For this, we look at a straight pipeline of $N$ identical tubes. As for the networks considered in the previous section, see Fig. 2, all tubes are of the same radius $R$ and length $\ell$, in accordance with observations of diminishing radius hierarchy in higher order veins ronellenfitsch_topological_2015 ; sack_leaf_2013 . Metabolites and fluid are flowing into the first tube $Q_{\text{in}}$ and fluid is leaving at a constant rate $Q_{\text{out}}$ at every node between adjacent tubes. Metabolites cannot exit at nodes but remain in the fluid until the very end of the pipeline or are absorbed. Also the fluid inflow rate and total fluid outflow rate are equal, i.e. $Q_{\text{out}}=Q_{\text{in}}/N$. This results in a constant decrease in flow rate by $Q_{\text{in}}/N$ from one tube to the subsequent. To translate this to cross-sectionally averaged flow velocities, which are the flow properties determining absorption, we use $U=Q/\pi R^{2}$. Consequently, the flow velocity in segment $m+1$ is given by $U_{m+1}=U_{m}-Q_{\text{in}}/\pi R^{2}N$.

The outflux of metabolites from one tube is equal to the influx of metabolites in the subsequent tube $J_{\text{out},m}=J_{\text{in},m+1}$, since all tubes have the same radius $R$. If $\pi R^{2}J_{0}$ is the total amount of metabolites flowing into the first tube, then only the fraction $1-\hat{\phi}_{1}$ is flowing out while the fraction $\hat{\phi}_{1}$ is absorbed. Generalizing we determine the absorption in tube $m$ as

The state of optimally uniform absorption is now defined by absorption in subsequent tubes being equal. We use this constraint to determine the inflow rate that corresponds to the optimally uniform supply pattern. Using Eq. (8) to write the absorption in the $m+1$ tube as a function of the absorption in the previous tube and inserting the equality constraint, we arrive at an expression including absorption capacities only, $(1-\hat{\phi}_{m})\hat{\phi}_{m+1}=\hat{\phi}_{m}$. Inserting the simplified expression for the absorption capacity from Eq. (7), we find the scaling law determining the optimal inflow rate to yield uniform absorption

where $L=N\ell$ denotes the length of the sequence of pipes. Note, that this condition is independent of which tube segment $m$ is considered. The absorption is uniform along the entire sequence of tubes.

Our toy model is set up to capture the essential flow and transport characteristic along the rows of a two-dimensional network excerpt modulus a geometrical factor. To confirm that the same functional dependence of the optimal inflow rate holds for two-dimensional networks, we return to our simulations of rectangular two-dimensional networks. For a given parameter choice, we vary the inflow rate and determine its optimal value by the minimal variance in absorption. We independently vary the absorption parameter $\gamma$ and the tube radius $R$, equal for all tubes right now, as well as the overall size of the network $L$, see Fig. 3. While one parameter was varied, the other parameters were kept constant. To cover a large parameter range, base parameters values are chosen as $\ell=$0.1\text{\,}\mathrm{mm}$$, $R=$3\text{\,}\mathrm{\SIUnitSymbolMicro m}$$, and $10\text{\,}{\mathrm{m}}^{-1}$, see supplementary material S5. For each parameter combination the inflow rate was varied in step sizes of $\Delta Q=$1.5\text{\times}{10}^{-6}\text{\,}{\mathrm{mm}}^{3}\text{\,}{\mathrm{s}}^{-1}$$. Note, to increase the overall size of the network additional nodes where added. Therefore, $Q_{\text{out}}$ decreased, and thus the overall flow velocity gradient decreased as well. Each run over a parameter combination was repeated 15 times with different random Gaussian node perturbations of a rectangular two-dimensional network with tubes of the same radii $R$. We find a linear scaling between the optimal inflow rate and the absorption parameter, the radius, and the overall size of the network in agreement with the scaling law’s prediction. Even more, if we multiply the optimal inflow rate derived above for the one-dimensional tube network by a geometrical factor $\Gamma$ taking into account the two-dimensional network geometry the numerical results follow exactly the analytical prediction. The geometrical factor is product of three terms $\Gamma=\Gamma_{\text{L}}\cdot\Gamma_{\text{AR}}\cdot\Gamma_{\text{IF}}$, where the first term is correcting the length of the network and the later two are needed to correctly link the velocity profiles with the inflow in the network. The total length of the network in flow direction is effectively shortened as tubes of the two dimensional network are not connected with a angle of $180^{\circ}$ as in the toy model but the network is a tessellation of approximately equilateral triangles. The length of the network has thus to be shortened corresponding to the ratio of height and side length of such a triangle with $\Gamma_{\text{L}}=\sqrt{2}/3$. Considering the rows of nodes, the inflow of fluid in 16 nodes is distributed to 17 nodes in the next row. The flow in tubes connecting these two layers of nodes is thus reduced by the ratio $\Gamma_{\text{AR}}=16/17$. As the total inflow is the inflow over the complete width of the network, the optimal inflow has to flow into every of the 2 to 3 tubes connected to the 16 inflow node giving rise to $\Gamma_{\text{IF}}=47$.

We found that a global change in the total fluid inflow rate is the most important control mechanisms to generate uniform supply patterns in a tissue pervaded by a transport network. How does a network architecture need to change to compensate low or high inflow rates? How much more can we minimize the variance in absorption even if fluid flow rate is optimal? To answer these questions we now optimize our previously found supply patterns by allowing for local dilation or contraction of tubes starting with the randomized networks introduced above. In addition to tube dilation and contraction, changes to the network architecture by discarding entire tubes are allowed. A tube is regarded as cut, if it radius is reduced below a threshold of $0.05\text{\,}\mathrm{\SIUnitSymbolMicro m}$, compared to an average tube radius of R=$3\text{\,}\mathrm{\SIUnitSymbolMicro m}$. While locally changing the network architecture, we keep the total amount of material $M=\sum_{i}R_{i}\ell_{i}$ within the network constant as we redistribute changes in $M$ over all radii equally. For this we numerically optimize the network topology using Monte-Carlo methods, explained in detail in supplementary material S3B.

We optimize the network architecture regarding uniform tissue supply for the cases of low, high and optimal inflow rate, see Fig. 4. In all three cases overall variance in absorption was successfully decreased. For low inflow rates, we observe a contraction of tubes near the influx nodes and an expansion of tubes toward the opposing end. Contraction of tubes speeds up flow velocities thus reducing otherwise dominating absorption close to the inflow nodes and thereby making metabolites available for absorption further onwards. The increase in absorption follows spatially the rapid increase in radius. This indicates that shifts in the radius distribution impact the local absorption profile strongly. For high inflow rates, we observe the opposite optimization mechanism. Tubes dilate close to the inflow and contract toward the opposing end. Here, dilation decreases flow rate and increases the absorption early on, while at the same time reducing the amount of metabolite flushing through. For the optimal inflow rate, we observe slight dilation near the inflow as well as near the outflow nodes. These changes correct for network artifacts arising from the chosen rectangular form of the excerpt. In all optimized networks we find small fluctuations in the absorption pattern which result from the randomized node positions and random tube lengths.