Advection and diffusion of substances in biological tissues with complex vascular networks

Abstract

For highly diffusive solutes the kinetics of blood-tissue exchange is only poorly represented by a model consisting of sets of independent parallel capillary-tissue units. We constructed a more realistic multicapillary network model conforming statistically to morphometric data. Flows through the tortuous paths in the network were calculated based on constant resistance per unit length throughout the network and the resulting advective intracapillary velocity field was used as a framework for describing the extravascular diffusion of a substance for which there is no barrier or permeability limitation. Simulated impulse responses from the system, analogous to tracer water outflow dilution curves, showed flow-limited behavior over a range of flows from about 2 to 5 ml min(-1) g(-1), as is observed for water in the heart in vivo. The present model serves as a reference standard against which to evaluate computationally simpler, less physically realistic models. The simulated outflow curves from the network model, like experimental water curves, were matched to outflow curves from the commonly used axially distributed models only by setting the capillary wall permeability-surface area (PS) to a value so artifactually low that it is incompatible with the experimental observations that transport is flow limited. However, simple axially distributed models with appropriately high PSs will fit water outflow dilution curves if axial diffusion coefficients are set at high enough values to account for enhanced dispersion due to the complex geometry of the capillary network. Without incorporating this enhanced dispersion, when applied to experimental curves over a range of flows, the simpler models give a false inference that there is recruitment of capillary surface area with increasing flow. Thus distributed models must account for diffusional as well as permeation processes to provide physiologically appropriate parameter estimates.

FIGURE 1

Hexagonal distribution of axially aligned segments. The network is modeled as a periodic arrangement of 16 axial segment positions in the x–y plane.

FIGURE 2

The network model is arranged periodically in the z (axial) direction. An arteriolar source occurs at a z position chosen from a Gaussian distribution with mean of 500 μm and standard deviation of 50 μm. Venular sinks occur at 125±25 and 875±25 μm.

FIGURE 3

Topology of a network generated by randomly placing 160 cross-connecting segments between nearest-neighbor axial vessels. This is a two-dimensional representation of the three-dimensional model. The axial vessel number (see Fig. 1) is plotted vs the axial position. Nodes are indicated as small circles; the large filled circle indicates the arteriolar source; the large open circles indicate venular sinks. The minimum distance along a capillary between cross connections d_min is 15 μm.

FIGURE 4

Hexagonal grid in the x–y plane used for numerical computation. The positions of axial vessels are indicated by closed circles.

FIGURE 5

Concentration at grid point is denoted by C_i. Concentrations in the six surrounding points in the x–y plane are denoted by $C_i^1$ through $C_i^6$.

FIGURE 6

Intravascular tracer washout simulation from the network model for various flows. The shapes of the curves are similar. But as flow increases, the effects of diffusion are diminished. At low flows the outflow curve is relatively smooth, while at the highest flow shown (10 ml min⁻¹ g⁻¹), the outflow curve is not unim odal.

FIGURE 7

Simulated concentration profiles in a slice through the x-z plane are shown following injection of an impulse bolus at the arterial inflow. Each profile corresponds to a different time in the simulation. The injection point is indicated in the top slice (t=0) by an arrow. The color bar to the right indicates concentration in arbitrary units. Tracer spreads from the injection point into the venular zones with concentration remaining relatively constant in the x (and y) direction. This simulation was performed using a flow of 2 ml min⁻¹ g⁻¹.

FIGURE 8

Tracer washout simulation from the network model for various flows. Lower flow washouts show an earlier peak due to diffusional shunting. The diffusion coefficient D is 10⁻⁵cm²s⁻¹.

FIGURE 9

Tracer washout simulation from the network model for various flows. A flow-limited regime exists for flows between 2 and 5 ml min⁻¹ g⁻¹ in which the normalized outflow curves appear similar. At higher flow (10 ml min⁻¹ g⁻¹) an early peak occurs due to a bolus of tracer that remains intravascular during its passage through the network. The diffusion coefficient D is 10⁻⁵cm²s⁻¹.

FIGURE 10

Estimated permeability–surface area product PS_e obtained from fitting simulated outflow data to the axially distributed two-region model described in Bassingthwaighte et al. (Ref. 10). Diamonds correspond to model fits using a heterogeneous flow model (20 parallel units) with a relative dispersion of flows equal to 40%. The circles correspond to model fits obtained from a single-tube axial distributed model.

FIGURE 11

Estimated extravascular dispersion coefficient D_t obtained from fitting simulated outflow data to an axially distributed two-region single-capillary model (see Refs. and 10) with the capillary PS_e set to the value of 500 ml g⁻¹ min⁻¹, corresponding to the expected effective PS based on Eq. (1). Plasma diffusion D_p is 10⁻⁵cm²s⁻¹.

FIGURE 12

Various axially distributed model fits are shown for the simulated outflow curve from the 3D network model. The solid line indicates output from the network model calculated at a flow of 5 ml min⁻¹ g⁻¹. The axial model approximates the network model behavior only for unrealistically low PS or high D_t.

FIGURE 13

Comparison of upwind method to minmod flux-limiter method. Outflow concentration curves following an impulse injection are plotted for the axial-distributed model with one advecting region and no permeation. Theoretical output curve (solid arrow) is a delayed impulse. The upwind method introduces substantially more dispersion than the minmod method. Grid size g=5 μm and a time step k = 0.0025 s.

FIGURE 14

Numerical and analytic solutions to Eqs. (B2) and (B3) are compared. The numerical solution (dashed line) is virtually indistinguishable from the analytic solution (solid line). The numerical method is a splitting scheme which uses the upwind method for the advection step. Grid size g=5 μm and a time step k= 0.0025 s.