A Matlab-based suite of tools for building models of microvascular transport and exchange from 3D structural information. The package is not a stand-alone model. Rather, these tools are generally adaptable to a wide range of problems involving advection, permeation, and diffusive transport in microvascular networks. To implement these tools, new users will need to become somewhat familiar with the methodology and the scripts, which are detailed below. Note that for large scale applications these scripts are best used with the MATLAB Compiler. The methodology is detailed in the published reference.
The methodology is outlined as an example problem simulating the washout of a permeant tracer from the rat cerebral microcirculation.
Components of Package
1. Description of geometry:Modeling microvascular transport starts with a description of the geometry of the network. In this package microvascular structure is modeled as a network of straight cylindrical vessels. A number of network structures are available at the Microvascular 3D Structural Information Page provided by Timothy Secomb at the University of Arizona. Here we use the example obtained from rat brain [Secomb, T.W., Hsu, R., Beamer, N.B. and Coull, B.M.Microcirculation 7:237-247, 2000], illustrated to the right.
For this example, the vessels sizes, connectivity, and flows are specified in the file net_data. Network data is loaded into MATLAB with the following commands:
load net_data.txt
S = net_data(:,2:3); % lists the 2 nodes for each segment
d = net_data(:,4); % vessel diameters (micrometers)
l = net_data(:,5); % vessel lengths (micrometers)
x1 = net_data(:,6:8); % start node positions
x2 = net_data(:,9:11); % end node positions
F = net_data(:,12); % flow in segment (micrometers^3 sec^{-1})
Npoint = max(max(S)); % number of nodes
The flows in the file “net_data” have been computed as described, (Ann. Biomed. Eng. 29:837-843, 2001). The function plot_network.m can be used to generate an image of the network.
2. The next step is to integrate the geometry of the vascular network within each volume element of the numerical grid. This is done using the scriptintegrate_network.m. The syntax for calling this script is illustrated below:
Lx = 150; % length of region in x direction (micrometers)
Ly = 160; % length of region in y direction (micrometers)
Lz = 140; % length of region in z direction (micrometers)
h = 10; % grid size (micrometers)
Nx = Lx/h;
Ny = Ly/h;
Nz = Lz/h;
N = Nx*Ny*Nz; % Number of volume elements
[V1,V2,Ax,Ay,Az] = integrate_network(x1,x2,d,h,Nx,Ny,Nz);
The function “integrate_network.m” returns the following vectors of length N: V1, V2, Ax, Ay, and Az, where N is the number of grid elements. V1 is the volume of the intravascular space; V2 is the volume of the extravascular space; Ax is the area of extravascular space on the x-face of the volume element; Ay is the area of extravascular space on the y-face of the volume element; and Az is the area of extravascular space on the z-face of the volume element. The geometry calculations are required for constructing the diffusion operator in Step 3. (For faster computing, it is recommended that “integrate network.m” be compiled.)
3. Next, the diffusion operator “Kdiff”, a sparse matrix of size (N,N), is constructed using the script initialize_Kdiff.m:
[Kdiff] = initialize_Kdiff(Nx,Ny,Nz,V2,Ax,Ay,Az,h);
Using this operator to simulate diffusion is described below (Step 5).
4. The final step in setting up the transport operators is to compute the avection operator for the intravascular space and the permeation operators which transport material between the intravascular and extravascular spaces. First, the following variables are input:
in_node = [1 2 3 8 11 49]; % Inflow nodes
out_node = [4 5 6 7 9 10]; % Outflow nodes
in_vessel = [1 8 9 10 13 14]; % Inflow vessel numbers
out_vessel = [6 15 28 29 31 35]; % Outflow vessel numbers
P = 10; % vessel permeability (um^{1} sec^{-1})
Ft = 0.5*Lx*Ly*Lz / 60; % flow into region um^{3} sec^{-1}
F = Ft .* F ./ sum(abs(F(in_vessel)));
Then the operators are calculated using the script discretize_network.m:
[Vs,flows,dt,Kadv,LQ,in_seg,out_seg] = ….
discretize_network(S,d,x1,x2,h,F,Npoint,Lx,Ly,Lz,P,in_node,out_node,V1,V2);
Where Vs is a vector listing the volume of vessel segments; flows lists the flows in vessel segments; dt is the numerical time step for advection operator; Kadv is the advection operator; LQ is the mass transfer operator for permeation between vessel segments and tissue volume elements; in_seg is the list of indices of inflow vessel segments; and out_seg is the list of indices of outflow vessel segments.
5. Example simulation: The impulse response washout is computed with a simple time-stepping algorithm:
D = 500; % Extravascular diffusion coefficient (um^{2} sec^{-1})
% Setting initial conditions:
Ns = length(Vs);
N = length(V1);
Cs = zeros(Ns,1); % Concentration in vessel segments
Ct = zeros(N,1); % Concentration in tissue volume elements
Cs(in_seg) = flows(in_seg)’ ./ Vs(in_seg); % Impulse at inflow segments
Cout = [];
i_sample = 10;
for i = 1:10/dt
% Advection
Cs = Kadv*Cs;
% Permeation
C = [Cs’ Ct’]';
C = C + dt.*(LQ*C);
Cs = C(1:Ns);
Ct = C(Ns+1:Ns+N);
% Diffusion
Ct = Ct + (dt*D).*Kdiff*Ct;
if mod(i,i_sample) == 0
Cout = [Cout sum(Cs(out_seg).*flows(out_seg)’)/ Ft];
end
end
t = i_sample.*dt.*(0:length(Cout)-1); plot(t,Cout); grid on;
You should get a plot that looks like this: