Wet-tip versus dry-tip regimes of osmotically driven fluid flow

The secretion of osmolytes into a lumen and thereby caused osmotic water inflow can drive fluid flows in organs without a mechanical pump. Such fluids include saliva, sweat, pancreatic juice and bile. The effects of elevated fluid pressure and the associated mechanical limitations of organ function remain largely unknown since fluid pressure is difficult to measure inside tiny secretory channels in vivo. We consider the pressure profile of the coupled osmolyte-flow problem in a secretory channel with a closed tip and an open outlet. Importantly, the entire lateral boundary acts as a dynamic fluid source, the strength of which self-organizes through feedback from the emergent pressure solution itself. We derive analytical solutions and compare them to numerical simulations of the problem in three-dimensional space. The theoretical results reveal a phase boundary in a four-dimensional parameter space separating the commonly considered regime with steady flow all along the channel, here termed “wet-tip” regime, from a “dry-tip” regime suffering ceased flow downstream from the closed tip. We propose a relation between the predicted phase boundary and the onset of cholestasis, a pathological liver condition with reduced bile outflow. The phase boundary also sets an intrinsic length scale for the channel which could act as a length sensor during organ growth.

www.nature.com/scientificreports www.nature.com/scientificreports/ some considered the bi-directionally coupled osmolyte-flow system but had to apply simplifying assumptions for solving the model, i.e., either negligible hydrodynamic pressure compared to the osmotic pressure, thereby uncoupling the problem 21-23 , or considering a limit case of small solute secretion rate 4 and, additionally, diffusion dominating over advection 19 , or replacing the mixed set of boundary conditions by pure velocity boundary conditions and not solving for the pressure profile explicitly [12][13][14][15][16][17][18] .
In this work we explore the effects of the mechanical feedback loops, see Fig. 1A, by taking a combined analytical and numerical approach to the bi-directionally coupled system and evaluate its implications for bile secretion in the liver. In the following sections we derive the model equations for the distributed secretion-advection problem of osmotically driven fluid flow, then map the model equations at steady-state to a single ordinary differential equation (ODE) of Abel type and obtain relevant particular solutions which predict the fluid pressure profile as a function of geometrical and material parameters. When parameter values are changed, the analytic solution reveals a transition between continuous flow from the closed tip on versus stalled flow with osmolyte accumulation near the closed tip, which we will term "wet"-tip and "dry"-tip parameter regimes, respectively. These analytical results are confirmed and complemented by numerical simulations of the full model in three-dimensional space. We finally apply the generic model to distributed bile secretion and flow in the liver and propose a relation to biliary pathology. only covered with water pores but also with active transporters for osmolytes with flux density g. For the example of bile canaliculi in the liver, the lateral cylinder surface is constituted by the apical membrane of many hepatocytes 7 . Tight junctions seal the gaps between the apical membranes of neighboring hepatocytes. Energyconsuming transmembrane pumps, including BSEP, MDR2, MRP2, actively transport osmolytes even against any osmolyte concentration gradient across the apical membrane 24 . We regard any process in the surrounding tissue at steady state and factored into the osmolyte secretion flux density g such that g becomes a fixed parameter, independent of osmolyte concentrations in lumen and surrounding cells. A cylindrical coordinate system φ = r z x ( , , ) is introduced and the z-axis is aligned along the centerline of the channel with its origin at the closed tip. This closed tip represents the dead end of secretory channels like the tips of bile canaliculi in the pericentral area of a liver lobule where no direct contact exists between the channel and any circulatory system 7 . Moreover, for the cleft system spanned by the lateral membranes of epithelial cells, this closed tip corresponds to a sealed tight junction preventing any leakage 23 . Hence, the closed tip ( = z 0) is modeled as a rigid non-permeable wall while the other end ( = z L) is an open outlet, representing e.g. a large bile duct or a reservoir.

Mathematical
This model system can be studied by three-dimensional computational fluid dynamics simulations of the Navier-Stokes equations, see Fig. 1B,B' and Methods. A typical resulting flow pattern is illustrated in Fig. 1B with streamlines on the center plane of the computed three-dimensional velocity field colored by the magnitude of the local velocity vector. Note the beginning of streamlines at the lateral cylinder surface, reflecting the distribution of osmotic water sources, and the steep velocity increase towards the open outlet. As flow is observed continuously all the way from the very tip of the channel, we will below introduce the term "wet-tip regime" for the underlying parameter choices. When repeating the simulation for different parameter values, then also a different flow regime is encountered in which flow ceases in a self-organized manner near the tip of the channel, see Fig. 1B' . Below, we will term this regime "dry-tip" as no fluid enters the channel at steady state for some portion downstream of the tip. The data from our numerical analysis of the three-dimensional problem will be discussed in a corresponding section below. We here aim to understand the transition between these flow patterns as well as their properties and parameter dependencies.
Towards analytical analysis, we first derive a one-dimensional model of coupled fluid and osmolyte transport by considering (i) flow in a narrow channel (blue arrows in Fig. 1A), (ii) osmotic driving and pressure feedback (upper red arrows in Fig. 1A) and (iii) the mass-balance of osmolytes (green and lower red arrows in Fig. 1A). A one-dimensional model here captures all important aspects since the cross-sectional Péclet number is of the order −  10 1 4 for relevant systems, see next Section and Methods. The derivation is given in a Methods section and yields coupled equations for the cross-section averaged velocity profile w z t ( , ) following Darcy's law, the pressure profile p z t ( , ) in the continuity equation and the osmolyte concentration profile c z t ( , ) governed by mass-balance: and  Table 1.
A special case of this general model has recently been applied to bile flow in the mouse liver and was validated using live intra-vital microscopy data 25 . Therein, specific parameter values inferred from untreated mice together with measured solute secretion upon treatment allowed to correctly predict bile flow in treated mice. Complementing these studies, we here aim at understanding the model's behavior, notably different flow regimes, for any parameter combination and also including heterogeneous osmolyte concentration profiles.
Mapping of the general steady state problem to an Abel equation and analytical solution for limit cases. To gain analytical insight into parameter dependencies, we aim at a closed-form solution for the pressure profile and concomitant velocity and water influx profiles. In Eq. 3, a potential time dependency in c z t ( , ) can only result from time-dependent solute secretion g z t ( , ), which then renders all other variables of the fluid dynamics problem time-dependent. However, for temporally constant solute secretion g(z) and as a consequence of the saturating role of the negative feedback loops without time delays (see Fig. 1A), both the fluid flow and the concentration profile converge to a stable steady state, which is reproduced with high accuracy by our numerical simulations (see Fig. 1B,C). We are restricting our analytical analysis to this steady state. Additionally, one can account for slow time dependencies in parameters or solute secretion g z t ( , ) by adiabatic approximation and substitution of slow time dependencies into the solutions derived here.
Moreover, there are two regimes of the axial Péclet number (Pé = Lw /D) in which either diffusion (Pé  1) or advection (Pé  1) is the dominant transport mechanism of osmolytes such that ongoing secretion with positive g(z) gets balanced. The range of different secretory systems with channel lengths between 10 μm (lateral cleft in www.nature.com/scientificreports www.nature.com/scientificreports/ epithelia) and 1mm (intrahepatic bile canaliculi or pancreatic ducts), with flow velocities around the μm/s range 25 and with diffusion constants between 10 μm 2 /s (large molecules like bile salts) and 10 3 μm 2 /s (K + ions in water) results in Péclet numbers that span across both regimes and we therefore consider both regimes. The regime of dominating diffusion (Pé  1), leading to uniform = c const and allowing to derive a general solution for any choice of g(z), will be analyzed in a subsequent section below while mixed cases of comparable contributions by diffusion and advection (Pé~1) can be studied numerically. Here, we first consider the more difficult regime of dominating advection (Pé  1) over diffusion and analyze Eq. 3 without the diffusion term and at steady state: Taking the definite integral with respect to the axial coordinate from 0 to z on both sides of Eq. 4, we can solve for solute concentration c(z) at steady state: a dp z dz 8 represents the cumulative upstream solute secretion flux and Eq. 1 has been inserted to relate c and p. Substituting c(z) from Eq. 5 into Eq. 2, yields a single nonlinear ordinary differential equation for p(z): a dp z dz To appreciate the fundamental scaling relations between the physical quantities, we partially non-dimensionalize Eq. 6 by introducing the relative coordinate  z and the non-dimensional parameter grouping M: a L 8 2 3 4 that represents the ratio of osmotic conductivity for water (membrane surface area times membrane permeability) to hydrodynamic conductivity and is known as the Münch number 26 . Eq. 6 then reads: We first consider two limit cases of small and large Münch number: , then integration yields 2. In the case when  M 1 and for physiological pressure profiles with bounded values and derivatives, the value of the bracket in Eq. 8 has to vanish. Reverting the substitution of c(z) from Eq. 5 for this bracket, we obtain Hence, the right hand side in Eq. 2 becomes zero and the gradient of the velocity vanishes, keeping the velocity at its boundary value = w 0. So, all parameter combinations corresponding to a large Münch number prevent normal outflow in agreement with the definition of M.
For arbitrary M, we can integrate Eq. 8 with respect to  z and obtain a special nonlinear ordinary differential equation of Abel type (Eq. 26), see Methods, for which particular solutions are known. However, Eq. 26 still poses difficulties for finding the general solution given an arbitrary solute secretion profile  g z ( ). Hence, we consider particular profiles  g z ( ) and then identify commonalities and differences between the solutions of relevant particular cases.
www.nature.com/scientificreports www.nature.com/scientificreports/ First, we choose the case that yields uniform concentration =  c z c ( ) , which in the following section permits us to look for the pressure solution  p z ( ) and corresponding solute secretion profile  g z ( ) independently. Then we extend this approach to non-homogeneous osmolyte profiles  c z ( ) in the subsequent section.
Analytical solution for uniform osmolyte concentration. We let = = RTc p const osm in Eqs 2 and 8, and apply the non-dimensionalization from Eq. 7, to get are fulfilled, yields the solution: This pressure profile is monotonously decreasing (solid curves in Fig. 1C,C') with its maximum p max located at the closed tip ( =  z 0): Hence the maximum hydrostatic pressure, which could potentially damage cells, is intrinsically limited by the osmotic pressure. This intrinsic limitation may safeguard the surrounding cells, that establish the osmotic pressure by actively pumping osmolytes into the lumen, as higher hydrostatic pressures interfering with cell geometry and cell functions would reduce the cells' pump efficiency and cause both pressures to relax. Another characteristic of the solution is the spatially averaged pressure as a function of Münch number M: max 0 1 Figure 2A,B show these dependencies of the fluid pressure on the non-dimensional parameter grouping M. Moreover, Taylor expansion of Eq. 12 reveals linear dependencies , both in agreement with the limit cases studied for the general case above. Notably, these qualitatively different limits, one with strong versus the other with vanishing fluid secretion, indicate the existence of two different flow regimes.
Calculating the osmolyte concentration from given maximum osmolyte secretion, see Methods, and inserting this relation Eq. 30 into Eq. 13 gives the emergent pressure range: max 0 where the pressure drop between the system boundaries (p max ) is a model result and depends on the solute secretion rate as g 0 as well as all other parameters, whereas in previously studied fluid flow models a pressure drop between the system boundaries had to be pre-defined. This model prediction reveals a non-trivial pressure dependence on the membrane permeability κ, exhibiting non-monotonic behavior, as κ affects the pressure explicitly and implicitly through M. Figure 2B shows this increasing pressure at low permeability values, when high water influx steepens the pressure gradient, versus decreasing pressure at high permeability values, when larger water inflow washes away the osmolyte and prevents its accumulation and concomitant pressure build-up.
As a consequence of conservation of fluid volume, we consistently observe that the cumulative inflow π 2 equals the terminal fluid outflow, linking Eq. 32 with Eq. 27, see Methods. In Fig. 2C we observe that the average water influx 〈 〉 j vanishes for large Münch number and hence the overall outflow of the channel stalls ( = →  w z ( 1) 0). The transition between the unobstructed flow regime at low Münch number and this stalled-flow regime is investigated in a subsequent section.
General solution for dominating osmolyte diffusion. Complementary to the above case of dominating advective transport over osmolyte diffusion, we next analyze the case (Pé  1) of osmolyte diffusion dominating over advective transport and heterogeneities in osmolyte secretion. Eq. 3 then describes the relaxation towards a stable uniform concentration state with = c z t c ( , ) . Hence all derivations and solutions from above, i.e. Eqs 11 to 15, apply here as well and, moreover, they apply for arbitrary forms of g(z) as long as diffusion dominates. The uniform concentration = c z t c ( , ) is then established as balance between cumulative secretion and terminal outflow. Correspondingly, integrating Eq. 3 over the channel volume and inserting =  w z ( 1) from Eq. 27 yields: Retrograde flow for heterogeneous osmolyte concentration profiles. In hyper-osmotic secretory processes, an osmolyte concentration gradient along the channel is observed that combines water secretion with re-absorption 21 . We therefore consider the case of an exponentially distributed osmolyte concentration profile at steady state: where c 1 , z 2 and  z 0 are parameters (see Fig. 3A). This case is solved in Eqs 34 to 36, see Methods. To study this solution, we consider the pressure profile concavity at the closed end (see Fig. 3B) and determine the Münch number ranges with characteristic behaviors of the system. We study several qualitatively different spatially-inhomogeneous osmolyte concentration profiles, but all with the same spatial average to ease comparison, see Fig. 3A. Figure 3B shows that the behavior of the solution changes with the Münch number. We choose specific M values that show different trends of pressure profile concavity and, thereby, cause qualitatively different profiles  Table 2. www.nature.com/scientificreports www.nature.com/scientificreports/ for pressure, axial velocity and osmotic water influx, see Fig. 3C,D". All shown spatial profiles are normalized by the respective maximal values for a uniform osmolyte concentration at the same Münch number, see previous section. This comparison shows that by varying fluid properties and channel geometry, one may achieve qualitatively different behaviors of the system, ranging from fluid secretion to its re-absorption.
In particular, a large decay length (black curves in Fig. 3A,C,D) yields spatial profiles similar to those observed in the uniform concentration case (see Fig. 1C). A more pronounced spatial inhomogeneity in the osmolyte profile causes a significant change in the flow profile compared to the uniform case. A rising osmolyte concentration profile with a short length scale generates retrograde flows to the closed end and a steep rise in the velocity profile at the outlet, peaking much higher than the uniform case (red curves in Fig. 3A,C,D). This fluid velocity profile results from the steep rise in the osmotic water influx close to the outlet. Monotonously decreasing osmolyte concentration profiles never cause retrograde flows but may create a water re-absorption region in some part of the channel (blue curves in Fig. 3A,C,D). Importantly, this result is independent of other system parameters including the osmotic pressure p osm . The particular value ≈ M 4 cr is a result of the considered feedback structure (see Fig. 1A) and governing equations (Eqs 2 and 3). From the definition of M (Eq. 7), we predict that the osmotic conductivity of the channel surface must remain smaller than 4 times the hydrodynamic conductivity of the channel to be able to completely drain the channel.
We now introduce the term "wet-tip regime" for an unobstructed steady state outflow from the very tip of the channel. This wet-tip regime occurs for < M M cr . Conversely, the term "dry-tip regime" will be used for stalled flow at the very tip of the channel and some part of the channel downstream of the tip. The term "dry" serves to indicate the lack of continuous water influx despite ongoing osmolyte secretion. This dry-tip regime occurs for > M M cr . The surface κ μ = M a L M ( , , , ) cr separates the wet-tip from the dry-tip regime in the 4-dimensional parameter space. Using the definition of M in Eq. 7, it is then straight forward to decide, for a given set of geometrical and material parameters, whether a channel can be drained all the way from its tip on.
The analytical results obtained above were derived for the case of uniform osmolyte concentration, corresponding to spatially increasing osmolyte secretion rate  g z ( ). An analogous result for a stagnant zone was derived by Rademaker et al. 20 . To verify the broader applicability of the criterion = M M cr for the regime boundary, we numerically solved the complementary case of uniform osmolyte secretion rate =  g z const ( ) as presented in Fig. 4 and found a qualitatively similar transition. Details of the numerical analysis are given in the Methods section. These analytical and numerical results together with the limiting solutions (Eqs 9 and 10) establish the wet-tip versus dry-tip transition as a general phenomenon.
Applications of the theory. As a quantitative test of our theory, it has recently been applied in an intra vital microscopy study of fluorescent tracer transport with bile flow in murine liver 25 . Parameter estimation for normal control mice and model predictions had been performed analytically and numerically using the simulation software Morpheus 27 . Comparison of model predictions for treatment with the commonly used analgesic acetaminophen (APAP) to experimental data had confirmed the model and identified an additional contractile driving mechanism especially in a zone near the closed tip. Potentially, this additional mechanism safeguards the mouse liver against the dry-tip regime as the Münch number calculated from their data (parameters κμ μ = = .
− p m 16 6 9710 www.nature.com/scientificreports www.nature.com/scientificreports/ therewith slightly exceeds M cr calculated here. This inferred Münch number is consistent with the strong curvature of the experimentally measured velocity profile.
On the other hand, typically wider channel radii in other systems will, with the third power of the radius, lead to smaller Münch numbers and we evaluate Eq. 7 for two such systems. For the classical problem of phloem flow in a plant leaf, we consider μ = a m 10 , = L mm 10 , μ = . ⋅ mPa s 1 5 , κ μ = ⋅ ⋅ − m Pa s 2 10 /( ) 5 following Table 3 in 16 . The resulting Münch number = .
M 0 048 is below M cr , confirming free flow 16 . For the kidney, models of osmotic fluid transport have addressed the question of urine concentration in the nephron due to water uptake from the nephron into the interstitium [28][29][30] . While the geometry of the interstitium is different from the tube considered here, a similar flow problem and dry-tip regime may in principle arise inside the interstitium. To estimate the Münch number for such a geometry around a nephron in rat, we consider μ ∼ a m 10 , ∼ L mm 1 , μ = . ⋅ mPa s 0 7 (viscosity of water at 37 °C) and membrane permeability κ for water at different positions in the range μ . … . ⋅ ⋅ − m Pa s (0 6 516 8) 10 /( ) 5 derived from experimental data (last column of Table 2 in 31 and assuming a typical osmotic pressure from 16 ). The resulting range of Münch numbers spans ⋅ − 7 10 5 to 0.06 below M cr , predicting free flow in the kidney.
Model implications for organ-size. Continuous drainage of the whole channel and thereby avoiding hyper-concentration and aggregation of solutes is of high importance for maintaining proper function of any secretory organ. Hence, it is further useful to derive the largest channel length that can still be drained completely by osmotic driving, so we invert Eq. 7 and obtain Should the channel length increase beyond L max then the dry-tip regime is reached for the channel region ] max , while the part of length L max near the open outlet sustains unobstructed outflow. Secretory organs may rely on this transition to control their size during growth in development and regeneration. Analogously, the maximum observed size for plant needles of around 5 cm has been explained 20 .

Discussion
Here we studied the biophysical effects of osmotic driving during fluid secretion using a mathematical modeling approach. Opposite to many fluid dynamics problems with given pressure gap across the system boundaries, here the pressure gap is an emergent property of the feedback loops in the system. Our numerical and analytical results have revealed the relation between tissue geometry and rheological properties of the secreted fluid and the function of secretory organs. The results identify two fundamental flow regimes: a free outflow from all the length of the channel, termed wet-tip regime, and a combination of free outflow with ceased flow within some portion of the channel near the closed end, the dry-tip regime.
The transition was found to depend on the value of a non-dimensional parameter grouping, the Münch number κμ = M L a 16 / 2 3 . In particular for  M 1, the total outflow π μ ∼ * a w L a M L M ( ) ( /( )) tanh 2 4 from the channel outlet yields κ μ a / 5 independent of L. This refutes intuition stemming from approximated uniform osmotic inflow that would cumulatively grow with L as ~κaL. The above counterintuitive parameter dependency originates from the competition of osmotic pressure against hydrostatic pressure (see pressure difference in Fig. 1A) for water flux through transmembrane water channels. This competition limits the functional part of the channel to the downstream portion of length L max while the upstream portion suffers the dry-tip condition. For comparison, an oversimplified model neglecting the pressure feedback (no pressure difference in Fig. 1A) would yield the predictions of the wet-tip regime for any parameter combination (Münch number) and miss the dry-tip regime.
By considering heterogeneous osmolyte concentration profiles, we have explored how osmosis sustains steady and spatially stratified patterns of water secretion and re-absorption. Our results to lowest order also extend to systems with weakly heterogeneous parameter profiles, a(z), κ(z) or μ(z), where these functions can be substituted into our pressure and velocity solutions and their spatial average enters the Münch number. From a methodological perspective, the deeper understanding of physiological behavior was obtained by the close integration of numerical and analytical methods as opposed to just applying numerical simulations alone. This integrated approach has been fruitfully applied before to questions of signaling, spatio-temporal patterning and function in liver cells 32 , pancreas 33 and tissue regeneration 34 .
We next explore implications of our theory for bile flow in the liver. Pathologies of bile flow range from solute aggregation into gallstones via stalled flow in cholestasis to bile infarct at too high bile fluid pressure 24 . Although recent theoretical work has considered bile flow in larger non-secreting ducts downstream of the gallbladder 35-37 , feedback-controlled bile formation as studied here was not part of that work. Using our Eq. 19, we can now predict a critical canalicular path length beyond which the dry-tip regime sets in. Then, the flow in part of the chan- ] max is stalled, which is reminiscent of cholestasis 24 . Our model suggests that the dry-tip regime may arise either due to reduction in the width of the channel (e.g. as in biliary constriction) thereby reducing ∼ L a max 3/2 or by increasing the viscosity of bile (e.g. as in Byler disease) thereby reducing μ ∼ − L max 1/2 . A gradual parameter shift over the course of disease progression corresponds to a continuous increase of M and Fig. 4 shows the resulting flow re-distribution.
As a subject of future work on more realistic channel geometries, it will be interesting to consider bile canaliculi in the vertebrate liver as an interconnected network including closed loops of branches. Such closed loops will result in many local pressure maxima and corresponding flow stagnation points, thereby reducing the potential output of the system while gaining robustness against perturbation of individual branches. Potentially, the www.nature.com/scientificreports www.nature.com/scientificreports/ stable fluid flow solutions with osmotic driving in looped network architectures are no longer restricted to stationary states as for the unbranched channel considered herein. These extensions may also connect to results for flow routing in capillary networks 38,39 . Whether channel geometry and fluid properties are actively adjusted to the capacity of osmotic driving in order to sustain an outflow from the channel's closed tip remains an interesting open question which can now be considered in the light of the wet-to-dry tip transition. Finally, it will be interesting to survey a broader set of secretory systems to compare their physiological range of Münch numbers and relate these to the threshold at ≈ M 4 cr .

Conclusions
We developed a mathematical model of the coupled osmolyte transport and fluid secretion problem and considered its steady-state dynamics. The obtained solutions for both homogeneous and inhomogeneous osmolyte distributions reproduce physiologically relevant flows. Our analytic results revealed a wet-to-dry tip transition between physiological outflow and pathological stalled-flow conditions based on a ratio between four parameters of the system, known as the Münch number gives a theoretical prediction for the onset of pathological conditions in vertebrate glands, like cholestasis in the liver, and may set a physical limit on organ size during growth and regeneration. Numerical simulations of the full-scale model confirmed these analytical predictions.

Methods
Model derivation. Fluid flow in long and narrow cylindrical channels allows for approximations to the Navier-Stokes equations. Here, = u v w u x ( ) ( , , ) is the fluid velocity in cylinder coordinates, p(z) is the fluid pressure and μ is the dynamic viscosity. We consider the steady laminar flow that preserves the axial symmetry of the channel geometry since the Reynolds numbers for the considered biological systems are low, e.g. for bile flow it was estimated as Re~10 −6 which lies 9 orders of magnitude below the threshold to turbulent flow 25 . Also, for  a L we can neglect the radial components of the fluid velocity, thereby reducing the problem to one spatial dimension. As it then suffices to analyze the cross-section average w of the axial velocity component, we set it proportional to the negative local pressure gradient, following Darcy's law Numerical solution of the three-dimensional full-scale problem. To study the system for arbitrary osmolyte secretion profiles g(z) and to test the applicability of the one-dimensional approximation for long and narrow channels with osmotic driving, numerical simulations of the coupled three-dimensional transport equations for water and solute were performed using the COMSOL software package (COMSOL Multiphysics ® v. 5.2, www.comsol.com, COMSOL AB, Stockholm, Sweden), see Fig. 1B,B' . The computational domain was defined as a cylinder with the following boundary conditions: no-flux at the closed tip; open-boundary at the outlet; solute secretion and solute-dependent water mass-inflow at the remaining boundaries. For the case of uniform osmolyte secretion rate (see Fig. 4B), the osmolyte influx g is defined as fixed uniform source in the simulation. For the uniform concentration case (see Figs 1C,C' and 4A), solute secretion was defined via a flexible source at the boundary that maintains a constant and homogeneous osmolyte concentration profile. The osmotic water influx density at the lateral boundaries is given by Eq. 21 and is used as input for the fluid dynamics solver. To avoid numeric convergence problems, ≥ j 0 was set explicitly, corresponding to fluid secretion and outflow from the channel. The Stokes equations were coupled to the diffusion-advection equation for osmolyte transport and solved iteratively at steady-state on the whole computational domain using a direct solver. We have verified the numerical accuracy of the simulations by showing convergence of a solution metric for increasingly refined meshes (data not shown). Table 2 with their respective references when appropriate. Regarding the membrane water permeability per unit osmotic pressure, we estimated the order of magnitude as follows. A flux density = ⋅ ⋅ − − j 151 10 m s 6 1 was measured 45 for a fixed osmolyte gradient in shrinking apical vesicles obtained from primary hepatocytes of male Fisher rats. The permeability κ is obtained by dividing this flux density by the pressure difference. In this experiment, sucrose has been used as artificial osmolyte and we approximate the resulting osmotic pressure difference using the osmotic potential of a related and well-studied sugar (see Appendix A.3 in 16  . Regarding fluid viscosity, we considered bile as an example and assume that it is slightly more viscous than water, adopting the experimentally measured value for human hepatic bile 36 . Also note that under disease conditions, the viscosity of gallbladder bile 36 may increase to a group median of ⋅ 3 mPa s 46 . However, only gallbladder