# A Graph-Based Framework for Multiscale Modeling of Physiological Transport M. Deepa Maheshvare, et al.

#### ByM. Deepa Maheshvare

Jan 12, 2022

GRAPHICAL ABSTRACT

## 1 Introduction

The physiological system is a complex network in which each organ forms a subsystem and different subsystems interact to maintain overall homeostasis of the body. Within each subsystem, functional networks exist at different levels of complexity. Metabolic and signaling networks within a cell, cell-to-cell communication networks in the extravascular region of tissue, cell-to-vessel communication networks, and the vascular network which couples the local dynamics to the global dynamics determine the functional behavior of all tissues. Bottom-up and top-down modeling approaches emulate the cellular dynamics and organ-level physiology. The ability to simultaneously capture the local and global dynamics by hierarchically bridging the communication networks existing across diverse scales is the key challenge in the holistic modeling of physiology.

Microscale cellular models use a bottom-up approach in which multi-omics data assimilated from high-throughput sources are employed to formulate and validate parameter-intensive kinetic models. These models capture the dependency of intracellular dynamics on metabolic steady states and flux distributions (König et al., 2012; Berndt et al., 2018a; Masid et al., 2020). Owing to cellular heterogeneity and the existence of intercellular communication, autocrine and paracrine signaling (Koh et al., 2012; Watts et al., 2016; Rao and Rizzo, 2020), the response elicited by a single cell, cannot be scaled to a cell community. Cell population models, therefore, employ discrete modeling approaches for examining cell-to-cell interactions such as intra- and inter-islet synchronization established by gap junctional coupling (Pedersen et al., 2005; Barua and Goel, 2016).

Macroscopic organ scale compartment models (Sorensen, 1985) employ a top-down modeling approach for predicting the bulk flow and elimination kinetics of biomolecules. These organ scale models rely on single-tube, parallel-tube, or tank-in-series approximations for idealizing distribution volume of blood into compartments (Gray and Tam, 1987). For improving the mechanistic understanding of tissue–vessel interaction, multiphase porous media–based models representing the tissue volume as intravascular and multi-region extravascular compartments (e.g., capillary–interstitial–parenchymal exchange unit) (Deussen and Bassingthwaighte, 1996; Chalhoub et al., 2007) emerged. However, these frameworks do not offer the possibility to fuse macroscale and microscale models. Consequently, the effect of network architecture on microperfusion patterns (Dolenšek et al., 2015) and its influence on the nutrient exchange cannot be investigated by these compartment models.

To overcome the limitations of compartment-based models, continuum approaches have been put forth for understanding the implications of morphological changes on the functional response of an organ. In the extravascular domain, continuum approaches are helpful in estimating the collective response of a tissue mass where the bulk of the tissue is smeared and treated as a homogenous domain. Although homogenization simplifies the complexity of the computational domain, the approach is limited in its ability to probe aspects such as the influence of heterogenous arrangement of cells on nutrient release patterns. In the intravascular domain, continuum approaches are suitable for analyzing the effects of dilation of blood vessels, deformations that occur as a result of fluid–structure interactions. With the advancements in imaging studies, the availability of microvascular datasets offers the possibility to model large-scale networks. However, discretizing the tortuous microvasculature vasculature for 3D modeling of the advection–dispersion physics gives rise to extensive computational overhead while employing continuum approaches. This has led to the development of graph- and hybrid-based approaches in which the vasculature infiltrating the tissue volume is represented as a one-dimensional network of pipes for modeling the flow and delivering resources over networks (Beard and Bassingthwaighte, 2000; Fang et al., 2008; Heaton et al., 2012; Kojic et al., 2017; Safaei et al., 2018; Erlich et al., 2019). In summary, while efforts have been made to resolve: 1) spatial heterogeneity at subcellular scale (Ii et al., 2019), 2) short-range communication in the microenvironment of cell communities, and 3) metabolic zonation in single sinusoid models (Berndt et al., 2018b), explicit models of long-range communications mediated by the vascular system remain underdeveloped at both intra-organ and inter-organ scale.

Toward this end, we need a scalable hierarchical framework that allows us to bridge diverse scales for modeling production, consumption, and distribution of biochemicals in a tissue microenvironment. We introduce a discrete modeling framework for simulating gradient-driven advection–dispersion-reaction physics of multispecies transport. Graph-theoretic approaches that have been proven successful in examining flow of information through large-scale real-world networks are applied (Kumar et al., 2019; Bellocchi and Geroliminis, 2020; Besse and Faye, 2021) in this study. We resort to discrete-vector calculus and use the operators defined on a finite-graph to spatially discretize and formulate the transport dynamics in the vascular domain as a “tank-in-series” model. Further, the computational domain for establishing the vessel-to-cell exchange and cellular dynamics within the cell are set up by combining ideas from other multiscale and Krogh cylinder models (Berndt et al., 2018b; Frost et al., 2019). Dynamics of nutrient exchange from the blood vessel to the layer of cells that lie in close proximity to the vessel surface is modeled; cellular reactions are explicitly modeled by representing cells as discrete volume nodes. Differential equations defining the interactions over nodal volumes embedded in the graph are solved by translating the physical domain into a metamodel in which the biophysical attributes are subsumed. This framework is suitable for the following key applications: 1) to reduce the computational cost involved in the spatial discretization of large tissue volumes (Section 3.1.2); our discrete approach is geared toward obtaining fast solution by reducing the system dimension, and the metamodel is scalable into any domain; 2) to probe the effect of flow topology on scalar transport and the sensitivity of concentration dynamics to network parameters and variations in physiological set points. (Section 3.2, 3.4); and 3) to assimilate multi-omics data from in vitro and in vivo studies and vascular topology from imaging studies (Section 3.3) for examining the influence of structural changes on the functional response of a tissue. Our graph-based discrete modeling framework differs from the existing approaches in the following aspects: conventional finite-difference, finite volume– or finite-element–based formulations operate on a continuous domain, and the equations are discretized, and approximate solutions are obtained. In our approach, we discretize the physical space and solve the equations on the graph which forms the discrete domain. This makes it possible to scale our framework to large networks and offers the flexibility to fuse multiscale models by encoding imaging data of vascular topology and omics data of cellular reactions to enhance systems-level understanding.

The outline of this study is as follows: the procedure followed in translating the capillary vasculature into a weighted graph and the preliminary assumptions considered for setting up the computational domain are discussed in Sections 2.1 and 2.2. The governing equations of the flow distribution and the mathematical formulation of the discrete model of advection–diffusion-reaction physics are presented (Sections 2.3, 2.4, 2.5). We use our framework to model two physiologically relevant test systems: 1) advection–dispersion dynamics of glucose transport in the microvasculature and 2) advection–dispersion-reaction dynamics of glucose–lactate exchange in the functionally coupled tissue–vascular domains (glucose–lactate dynamics is relevant in tumor metabolism where metabolic activity alters in tumor microenvironments (Yang et al., 2021) and in modeling fuel-stimulated insulin secretion (Jiang et al., 2007; Prentki et al., 2013)). By applying our method, we predict glucose drift in the islet vasculature and cross-validate the flow and concentration fields of the multiphysics simulation with COMSOL simulations (Section 3.1.1). We establish the cell–vessel link and predict the spatio-temporal evolution of glucose–lactate exchange in the extravascular and intravascular domains (Section 3.3). We test the model behavior for various flow topologies (Section 3.2) and different pressure drops and glucose doses (Section 3.4). The network configurations illustrated in the applications presented in this work are the capillary blood vessels.

## 2 Methodology

For setting up the discrete modeling framework to study the multiphysics coupling in multiscale systems, we start by introducing the steps involved in constructing the computational domain which is a metamodel of the physical space, dissection of the metamodel into subdomains which form the functional units, and formulation of the mathematical operators. The three main steps involved in our workflow are illustrated in Figure 1. 1) Create skeleton: The topological organization of the capillary network and biophysical characteristics such as length and diameter of the vessels in the network constitute the structural and anatomical characteristics relevant for setting up the computational domain. These characteristics are extracted in this step by skeletonizing the reconstructed vasculature (Section 2.1). 2) Solve flow distribution in the network: The physical space is translated into a weighted graph representing a hydraulic circuit. The pressure and flow fields are computed over the network by establishing the relationship between the node and edge entities of the graph using the Hagen–Poiseuille equation (Section 2.3). 3) Solve advection–dispersion-reaction dynamics: The metamodel is subdivided into functional units composed of cell and vessel subdomains. We combine multiple scales by coupling the uptake and release flux of the cell domain ωt with the carrier-mediated exchange flux occurring at cell–vessel interface. The metabolic reactions occurring at the cellular scale are modeled by biochemical rate laws, and the mass transport in the capillary domain Ωbv is described by coupling the cell-to-vessel influx or outflux with the gradient-driven advection–dispersion transport in the capillary domain Ωbv (Section 2.4, 2.5). Henceforth, the superscripts bv and t denote the parameters and variables defined in the capillary blood vessel and tissue domains, respectively.

FIGURE 1. Workflow involved in setting up the system for simulating the spatio-temporal evolution of biochemical species is illustrated: 1. Create skeleton: top, z-stack of microscopic images; middle, three-dimensional volume of the reconstructed vascular geometry; bottom, the length of blood vessel branches are color-coded. 2. Solve flow distribution in the network: bottom, the vascular network represented as a hydraulic circuit; middle, from the estimated nodal pressures, flow through the vessel branches are computed, and checked whether the mass conservation holds at the junctions; top: 3D visualization of the pressure gradient observed across the vasculature. 3. Solve advection–dispersion-reaction dynamics: top, creation of the computational domain which includes the blood vessel and the layer of the cells located in the vicinity of the outer surface of the capillary; middle, functional unit of the computational domain consisting of the vascular and cell sub-domains; bottom, concentration gradient observed across the vasculature after solving the advection–dispersion-reaction dynamics.

### 2.1 Construction of Capillary Networks

The topology of the blood vessels that exchange nutrients with the surrounding tissue is represented by the one-dimensional network of pipes. The biophysical attributes, such as diameter and length of each segment and the blood vessel fragment connecting two branching points, are quantified as edge weights (

$W$

) of the graph network (

$G$

). We investigate four network configurations in this work. The 3D vasculature of pancreatic islet displayed in Figure 2A and both the configurations of tumor network shown in Figures 2B,C were reconstructed from binary images of vessels which were examined in Chen et al. (Chen et al., 2020). The tiff stacks containing binary images were generated in their study by imaging the blood vessels, labeled with a fluorescent dye, using light-sheet microscopy followed by segmentation of the vessels in the ilastik (Sommer et al., 2011) toolkit which leverages a machine learning–based classification algorithm. We follow the workflow detailed below for generating a weighted graph from the binary images: 1) the multipage .tiff image was rendered into a 3D volume using the ray-casting technique available in 3D Slicer (Fedorov et al., 2012) by specifying the spacing of the image stack in the xyz directions. The dimensions of the input stack and domain size of the reconstructed volume are given in Table 1. 2) The largest connected region of the segmented volume was filtered using the Island effect available in Slicer, and the 3D object was exported in a stereolithography file (.stl) for skeletonization in the Vascular Modeling Toolkit library (VMTK) (Antiga et al., 2008). 3) vmtksurfaceclipper was employed to open the surface at the network inlet, and the vmtknetworkextraction algorithm was utilized to skeletonize the geometry. This yielded a network with nodes (vertices) (

$V$

) positioned at the N-furcation points or terminal ends and edges (

$E$

) formed by the vessel segment linking two nodes. 4) Segment length lbv and diameter dbv of the blood vessels were computed by tracing the shortest path between two nodes and extracting the maximum inscribed sphere radius (thickness), respectively. The 2D structure of the mesentery vasculature displayed in Figure 2D was generated by parsing the diameter, length of vessels, and topology information available in Amira mesh file provided in Esposito et al. (d’Esposito et al., 2018). Statistics of diameter and length distributions are shown in Figures 2E,F. We represent the skeletonized geometry of the capillaries as a weighted graph,

$G(V∈Rm,E∈Rn,W∈Rn)$

, for investigating the advection–diffusion-reaction dynamics. The cardinality of all the properties and the operators defined on

$G$

are summarized in Supplementary Table S1.

FIGURE 2. Microvascular network configurations: (A) pancreatic islet, (B) tumor design 1, (C) tumor design 2, and (D) mesentery generated after skeletonization of 3D volumes reconstructed from the image geometries studied in Chen et al. (2020) and d’Esposito et al. (2018). Diameter (E) and length (F) distributions of the vessel segments present in the four geometries.

TABLE 1. Specifications of the computational domain, values of flow boundary conditions, and the values of transport parameters used in the model.

### 2.2 Preliminary Assumptions

The time-dependent uptake and release of biochemicals by a tissue is determined by the gradient-driven transport across the capillary vessels which facilitate the transvascular exchange of solute with the tissue interstitium. Subsequently, the carrier-mediated sites present on the plasma membrane of cells aid the uptake of nutrient resources from the interstitial space for metabolism. Here, we introduce a graph-based mathematical framework for capturing blood-tissue exchange. The following simplifying assumptions are made in our model similar to the assumptions considered in the other multiscale studies: 1) assuming that the scalar concentration in the interstitial fluid attains rapid equilibrium with the concentration in the blood (Chalhoub et al., 2007), the interstitial compartment is not modeled; 2) the endothelial layer of the capillary surface is lined with metabolite transporters, which promote facilitated diffusion of biochemicals across the capillary wall; a similar approach has been presented in Heaton et al. (Heaton, 2012). We consider this as a reasonable assumption since a large fraction of the endocrine cells of the pancreas lie in close contact with the surface of the capillaries (Cohrs et al., 2017); 3) due to the deficiency of lymphatic drainage in the islets of Langerhans (Korsgren and Korsgren, 2016), fluid exchange with lymph vessels is not modeled (Thurber and Weissleder, 2011). Based on these considerations, we subdivide each blood vessel branch and the layer of tissue surrounding the capillary into discrete functional units, diagrammed in Figure 3D. The molar transport occurring through the blood vessel compartment of these functional units is modeled by the one-dimensional advection–dispersion equation (Taylor, 1954) that accounts for convection flux and the axial and radial diffusive flux of the solute molecules. The transcapillary exchange flux and the cellular processes occurring in the tissue compartment of each functional unit are modeled by rate expressions that capture the kinetics of metabolite-specific transporters and enzyme-catalyzed reactions, respectively. The governing equations that model the inter-compartment dynamics of these well-mixed volumes embedded in the finite connected network representation of the capillary bed is detailed in Sections 2.4 and 2.5.

FIGURE 3. Schematic representation of functional units in the model. Computational domain for studying (A) flow distribution, (B) and (C) advective–dispersive transport of biochemicals in the discrete volumes of the blood vessel Ωbv, and (D) advective–dispersive-reactive transport of biochemicals in the discrete volumes of the blood vessel Ωbv and tissue domains ωt.

### 2.3 Mathematical Formulation of Flow Distribution in the Network

To simulate the spatio-temporal evolution of chemical species, we first solve for the flow field in the vascular branches using the approach generally applied in studies that focus on modeling flow distribution in branching networks (Kirkegaard and Kim, 2020; Poelma, 2017; d’Esposito et al., 2018; Erlich et al., 2019).

##### 2.3.1 Domain

To perform fluid flow simulations, the oriented graph depicted in Figure 3A was used as the computational domain for setting up the linear system of equations. As illustrated, each blood vessel branch was treated as an axially symmetric cylinder of axial length lbv and circular radius rbv derived from the diameter dbv extracted from the skeletonized geometry (Section 2.1).

##### 2.3.2 Equations

We consider blood as an incompressible, viscous Newtonian fluid and apply the Hagen–Poiseuille equation (Eq. 2) for modeling the conductance of an edge as a function of viscosity μ, radius rbv, and length lbv. The linear analysis of the flow distribution presented here is applicable for laminar flow, Re

$<$

1 in all segments, and this can be extended further to study the nonlinear rheology of blood using the procedure illustrated in Pries et al. (Pries et al., 1990). Therefore, the distribution of flow in a microvascular network is determined by the pressure gradient and the resistance offered to flow:

Here, pi and pj are the pressures at tail i and head j nodes of the oriented edge eij, and Gij is the conductance associated with eij.

$Qij=Gijpi−pj=πrbv48μlbvΔp.(2)$

The vector

$Q∈Rn$

of volumetric flow rates in n edges is defined in Eq. 3. The negative of the oriented incidence matrix

$M(G)∈Rn×m$

, denoted as M henceforth, is the gradient operator that acts on the vector

$P∈Rm$

of nodal pressures to result in the vector

$ΔP∈Rn$

of pressure gradients. We obtain Q by premultiplying ΔP with the diagonal conductance-matrix

$G∈Rn×n$

; this scales the pressure gradient across each edge by the corresponding edge conductance.

In our model, we consider blood as an incompressible fluid and determine the unknown nodal pressures by imposing mass conservation at all nodes. Consequently, the net flow at any given node i is zero (Eq. 4). Here, Qij is positive when flow enters node i

$Qi←j$

, negative when flow leaves node i

$Qi→j$

, and

$A(i)$

denotes the set of nodes that are adjacent to i.

Over m nodes in the network, the vector MTQ (

$MT∈Rm×n$

is a divergence operator which is given by the transpose of M) defines the flow conservation at all nodes excluding the terminal nodes, where the flow boundary conditions are specified. The non-zero entries of vector

$q̃∈Rm$

(Eq. 5) contain the values of inflow or outflow rates at the boundary nodes (Erlich et al., 2019). In addition to the flow rate, we specify one value of known pressure at the inlet pin or outlet pout. The values of these boundary conditions were specified based on experimental measurements of blood flow velocities reported in Table 1.

$q̃i=−Qin,ifi=inlet nodeQout,ifi=outlet node0,otherwise(5)$

Here, MTGM is a square matrix. After substitution of known pressures in the vector P, the columns of MTGM scaled by the values of known pressure are shifted to the RHS of Eq. 6. This operation results in a non-square matrix on the LHS of Eq. 6. The resulting system of linear equations is solved for the unknown nodal pressures by finding the pseudoinverse (Golub and Pereyra, 1973), which is the generalization of inverse for rank-deficient matrices. Pseudoinverse was computed in MATLAB using singular-value decomposition. From the estimated nodal pressures, the centerline velocity uij and the volumetric flow rates are computed using Eq. 2.

### 2.4 Advection–Dispersion of Chemical Species in the Blood Vessel

##### 2.4.1 Domain Discretization

The axial lines (edges) of the pipe network in Figure 3A were spatially discretized to set up the computational domain and study the transport of biochemical species in the microvasculature. For discretizing the edges into 1D elements, in Gmsh (Geuzaine and Remacle, 2009), the vasculature was represented as a geometry with point and line entities. The length of the mesh elements, denoted by h in Figure 3B, was constrained by specifying the maximum and minimum characteristic lengths, i.e., h ∈ (cl + δ, clδ). We derive the characteristic length cl based on the average diameter of a biological cell, approximately 11.5 μm, calculated from the volume ranges reported in the pancreas (Pisania et al., 2010; Parween et al., 2016), and we consider a deviation δ of 2.5 μm from cl.

##### 2.4.2 Domain Volume Elements

Each node in the discretized domain forms the center of the volume surrounding it. For instance, the volume of ith node

$Vibv$

(Equation 7) located at the bifurcation point, depicted in Figure 3B, is the sum of half-cylinder volumes formed between nodes i–j, where j

$∈A(i)$

(Reichold et al., 2009). Here,

$A(i)$

is the set {i-1, i+1, i+2} of nodes adjacent to i, and Aij and lij are the cross-sectional area and length of the cylindrical volume between nodes i–j, respectively. At the junction nodes, we merge the half-cylinder volumes of the adjoining edges and equate the sum to a spherical volume (Figure 3C). Therefore, each branch in the network is dissected into cylindrical elements, and the branches are assembled together by the spherical elements at the junctions. Since the rate of momentum transport is three orders of magnitude greater than the rate of mass transport, the accumulation term at any N-furcation junction is zero while solving for flow field and non-zero while solving the mass transport problem. As an example, for glucose species, the ratio of momentum to mass transport defined by the Schmidt number ν/D is around 4,000.

$Vibv=12∑j∈AiAijlij.(7)$
##### 2.4.3 Equations

The continuous formulation of the advection–dispersion physics is given by the following partial differential equation, which describes the solute mole balance:

$A∂Cbv∂t=∇.−QCbv+DA∇Cbv.(8)$

Here, A denotes the area of flow, D is the dispersion coefficient of species, Q denotes the flow field, and Cbv is the scalar concentration of a chemical species. To shift from the continuous to discrete counterpart, first we assign the scalar concentration field C0 to all vertices, and the flow field and dispersion coefficients form the edge weights. The net change in the molar concentration of a species

$Cibv$

in the control volume

$Vibv$

at the ith node of the blood vessel depends on the contributions from advective

$JiA$

and dispersive fluxes

$JiD$

presented in Equation 9. The influx and outflux of

$JiA$

are the summation of molar fluxes through the edges feeding into i (ji) and leaving i (il), respectively, and the molar flux through an edge is the flow-weighted concentration of the compartment from which the oriented edge originates (Chapman, 2015; Hošek and Volek, 2019). The third term contributes to the dispersive transport. Furthermore,

$A(i)$

is the set that contains nodes that are adjacent to node i,

$A+$

(i) is the set of in-neighbor nodes of i,

$A−$

(i) is the set of out-neighbor nodes of i, and Dik, Aik, and lik are the dispersion coefficient, cross-sectional area, and length of edge eik, respectively.

$VibvdCibvdt=∑j∈A+iCjbvQj→iin︸influx−∑l∈A−iCibvQi→lout︸outflux︸advection fluxJiA+∑k∈AiDikAiklik︸DikQCkbv−Cibv︸dispersion fluxJiD(9)$
$Dik=D̃ik+uik2rik248D̃ik.(10)$

The dispersion coefficient is determined from Aris–Taylor’s relation in Eq. 10. Here,

$D̃ik$

, uik, and rik are the diffusion coefficient of a species, centerline velocity, and radius of eik, respectively.

$VbvdCbvdt=MTQMoCbv−MTDQMCbv.(11)$

For setting up the discrete advection–diffusion equation on the network, we express the first and second order partial derivatives in Eq. 8 in terms of weighted advection Laplacian proposed by Rak et al. (Rak, 2017) and weighted diffusion Laplacian, respectively. This yields a system of ordinary differential equations shown in Eqs 11 and 12. In Eq. 11,

$DQ∈Rn×n$

is the diagonal matrix with diagonal entries the volumetric dispersion coefficient DA/l of each edge specified in Eq. 9,

$Q∈Rn×n$

is the diagonal matrix with diagonal entries the volumetric flow rate of each edge,

$Mo∈Rn×m$

is the modified incidence matrix (Rak, 2017; Erlich et al., 2019), MTDQM is the weighted dispersion Laplacian matrix

$LDG∈Rm×m$

, and MTQMo is the weighted advection Laplacian matrix

$LQG∈Rm×m$

(Eq. 12).

$VbvdCbvdt=LQCbv−LDCbv.(12)$

##### 2.5.1 Domain

We adapt the Krogh cylinder approach presented in multiscale models (Thurber and Weissleder, 2011; Chalhoub et al., 2007; Berndt et al., 2018b) and approximate the layer of tissue surrounding the vessels as a hollow cylindrical volume element ωt illustrated in Figure 3D. The outer diameter of the region ωt is given by the summation of dbv and dt. A value of 12.4 µm was used in our model for dt. This value is derived from the physical volume of an insulin-secreting beta cell (1020 μm3) (Finegood et al., 1995) located in the islets of pancreas.

##### 2.5.2 Equations
$VibvdCibvdt=∑j∈A+iCjbvQjiin−∑l∈A−iCibvQilout︸advection fluxJiA+∑k∈AiDikAiklikCkbv−Cibv︸dispersion fluxJiD−VicelljiE︸exchange fluxJiE.(13)$

In the last term of Eq. 13,

$jiE$

(Eq. 15) is the net exchange rate that governs the bidirectional transport between tissue compartment

$ωit$

and blood vessel compartment

$Ωibv$

of the functional unit. The uptake or release flux

$JiE$

is computed by multiplying

$jiE$

with

$Vicell$

since the maximal rates are often reported in per unit volume of a biological cell. This ensures mole balance when a species moves from

$ωit$

to

$Ωibv$

which differs in volume. Furthermore,

$Cicell$

is the concentration of species in the tissue cell that interacts with ith node of Ωbv, the half-saturation constant Km quantifies the affinity of a transporter protein or an enzyme for a metabolite, and Vm is the maximal rate of metabolite transport. When

$Ωibv$

is encompassed by

$ωit$

, interaction between the node associated with both the compartments exists, and

$JiE$

appears as a source or sink term in Eq. 13. At the junction nodes of the blood vessel where no interaction exists,

$JiE$

is zero.

$VbvdCbvdt=MTQMoCbv−MTDQMCbv−JE(14)$
$jiE=VmKmCibv−Cicell1+CicellKm+CibvKm.(15)$

As an example, we consider a minimal model of glucose metabolism in β-cells in the islets. Under basal conditions, the concentration of glucose and lactate are at the basal level in the bloodstream. In the fed state, glucose transporters sense the high blood glucose level and export glucose from Ωbv to Ωt. The glucose-to-lactate conversion in Ωt is presented as a lumped reaction in our model for simplicity. Ωt acts as a source of lactate, and the lactate transporter facilitates its release into the bloodstream. The direction of the exchange flux is dictated by the concentration gradient CbvCcell across the vessel wall.

### 2.6 Cell

The mole balance of each species in the cell domain ωt is modeled by

where N is the stoichiometric matrix, and v is the reaction flux vector. The rate expressions of the kinetic reactions and metabolite transporters and the values of the kinetic parameters used in the model are presented in Table 2. The values of half-saturation constants, Km and Vm, were chosen from panmin (König and Deepa Maheshvare, 2021), a minimal model of glucose metabolism and insulin secretion in the pancreatic β-cell. The Vm values from the minimal model were scaled in this study so that the reaction fluxes are comparable to the diffusive and convective fluxes. To apply our framework for predictive modeling in clinical applications, the Vm values can be determined by model calibration to achieve good agreement between model predictions and experimental measurements. Steady-state and transient values of the metabolite and flux distributions measured from biochemical assays can be used as the inputs for calibrating the tunable parameters in the model (Berndt et al., 2018b). Since the focus of this study is in introducing and illustrating the applicability of our mathematical framework for bridging multiple scales, the parameters were not calibrated in the glucose–lactate test system presented here.

TABLE 2. Reactions, rate laws, and parameter values of exchange and cellular reactions.

## 3 Results and Discussion

The discrete model framework developed allows bridging the cell-to-vessel exchange and explicitly modeling the cellular dynamics. In this section, the mathematical formulations presented in Section 2 are verified using several test systems, and the results are compared with finite element analysis in COMSOL. The analysis includes solving two steps that are coupled: first, a flow field analysis which involves computation of pressure and velocity distribution in the blood vessel branches; and second, analysis of the spatio-temporal evolution of the concentration fields.

We first present the results of the pressure gradient and flow distribution observed across the islet vasculature. After validating the results of nodal pressures and edge velocities with the results from COMSOL, we proceed with the simulations of advection–dispersion dynamics of glucose species in the blood vessel. Here, we compare the transient change in the blood glucose concentration obtained from our discrete model versus COMSOL simulations for islet and mesentery networks. Next, we investigate the effect of change in the perfusion pattern on glucose distribution by varying the inlet and outlet locations in the tumor vasculature. Furthermore, we examine the effect of different pressure gradients applied across the network and the effect of glucose doses supplied at the inlet on metabolite rise times observed in the islet vasculature.

All the 3D visualizations of flow and concentration fields presented in this article are rendered using vedo (Musy et al., 2021), a python-based module for analyzing and visualizing multidimensional point-cloud, mesh, and volume data.

### 3.1 Comparison of Flow and Concentration Fields

To illustrate how the results from our discrete formulation compare with the finite element implementation available in COMSOL, we first solve the static flow problem (Eq. 6) and use the flow profile for simulating the advection–dispersion dynamics (Eq. 12) in the islet vasculature.

The pressure and velocity fields computed across the islet vasculature are shown in Figures 4A,B. We observe a net pressure drop of 34.83 Pa for an inlet pressure and flow rate of 60 Pa and 3.76 nl/min, respectively. The velocity distribution lies in the range reported by Diez et al. (Diez et al., 2017). The conservation of flow at each node is cross-verified by computing the divergence of the flow field MTQ, and the visualization is presented in Figure 4C.

FIGURE 4. Distribution of the static properties of islet vasculature. (A) Pressure gradient from inlet to outlet and (B) directional flow from inlet to outlet. (C) Flow conservation at all nodes. The summation of flows through the edges entering and leaving each node, except the boundary nodes, adds to zero. Validation of static results of the islet vasculature. (D): Comparison of nodal pressures between the current work and COMSOL simulation. (E): Comparison of edge velocities between the current work and COMSOL simulation. Comparison of the concentration profiles at various locations in the islet and mesentery vasculatures. (F) Scalar concentrations observed at positions 12.07 μm (red), 54.85 μm (green), and 110.27 μm (blue) from the inlet node. (G) Scalar concentrations observed at positions 1,267.78 μm (blue), 1,692.3 μm (cyan), 3,476.96 μm (green), 5,085.33 μm (magenta), and 6,319.57 μm (blue) from the inlet node. Solid and dashed lines indicate the transient change in glucose concentration from our model and COMSOL, respectively.

##### 3.1.1 Validation in COMSOL Multiphysics

Here, we describe the simulation implementation of the coupled multiphysics problem in COMSOL to validate the flow fields and concentration fields of advection–dispersion simulations from our discrete model. Geometry: The generation of the geometrical model of the islet and mesentery vasculature (Figure 2A,D) was automated using an AutoLISP script (Mac, 2020). The coordinates of points and the connectivity information of the lines were specified as inputs for creating the CAD geometry. The DXF file containing the geometry data was imported into COMSOL, and the normal mesh size was used to generate the mesh elements. Parameters: Values of the diffusion coefficient of glucose species, viscosity of blood, inlet pressure and outlet flow rate specified in Table 1, and diameter of all branches in the vasculature were defined as input parameters. Computation: Fluid flow was studied as a static problem in the Pipe Flow Module considering blood as a Newtonian fluid. This stationary problem was solved as a linear system using the direct solver by specifying the pressure and flow boundary conditions. Then, one-way coupling of the flow physics was established with the Transport of Dilute Species in Pipes module to solve for the time-dependent advection–dispersion physics of glucose species in the finite element solver. For this transient simulation, a value of 5 mM was used for the Dirichlet boundary condition defined at the inlet, and the mass outflow was modeled by setting the diffusive flux to zero. The initial concentration was set to zero in the volume elements of our discrete framework. In COMSOL, the concentration was initialized to zero, and a smoothed step function was applied to avoid discontinuity with the boundary condition.

Figures 4D,E show that the nodal pressures and edge velocities computed from our model are consistent with the results from COMSOL for the islet vasculature. The time-varying concentration profile of glucose species obtained from our discrete model is compared with COMSOL simulations at nodes highlighted in the inset of Figure 4F. After comparing the results of the islet vasculature containing 52 edges and 125 discretized elements in our discrete model, we further extend this analysis to compare the results (Figure 4G) of the large mesentery network with 489 edges and 9,033 discretized elements. Supplementary videos S1 and S2 show the comparison of evolution of concentration profiles from our model versus the results from COMSOL for the entire region of the islet and mesentery vasculature, respectively.

##### 3.1.2 Scalability

The procedure involved in solving the partial differential equations can be split into two steps: assembly step and solve step. In the assembly step, discretization is performed, and matrices are generated for the second step which solves the system. For the islet vasculature, we report the time taken for solving the stationary fluid flow and the transient advection–dispersion dynamics for a time span of 15 s. Our discrete method implementation in MATLAB takes 7 s, and the finite element solver in COMSOL takes 10 s. For the mesentery network, running the advection–dispersion dynamics for a time span of 300 s takes 267 and 168 s in MATLAB and COMSOL implementations, respectively. We provide the sparse Jacobian pattern as an additional input to ode15 s to speed up the compute time. Solving the same system in Julia using the QNDF method, which is a translation of MATLAB’s ode15 s, gives a speed up of 127x when compared to that of COMSOL. A relative error tolerance of 1e-3 is set for carrying out all simulations, and an absolute error tolerance of 1e-6 is set for both ode15 s and QNDF.

Some of the challenges involved in expanding this approach to a vascular network composed of arterioles, arteries, and veins would be in scaling flow parameters and resistances to flow (such as molecular interaction in capillary versus viscosity in larger vessels) to fit the formulation so that the model still retains its physical fidelity. Computationally, there can be stiffness in the differential equations and condition numbers of matrices may be affected, especially when both capillary and larger vessels are present in the model. This may need smaller time steps in simulation and regularization of affected matrices.

### 3.2 Influence of Flow Topology on Scalar Transport

Engineered perfusable vasculatures have been useful for investigating the structure–function relationship of complex vasculatures (Kinstlinger et al., 2020). Computational models that capture the influence of flow topology on the scalar transport will enable experimental scientists to design and test the efficacy of optimized drug delivery systems.

Motivated by the experiments carried out by Chen et al. (Chen et al., 2020) on microfluidic devices imprinted with blood vessel vasculatures, we study the effect of variation in perfusion patterns on the metabolite rise times tr in two different geometric configurations (Figures 2B,C) of the tumor vasculature examined in their study. Flow distribution was computed for an inlet pressure of 60 Pa and fluid flow rate of 23.8 nl/min. Spatial distribution of glucose species was obtained, and the rise time in each volume element

$Ωibv$

is computed by finding the time taken for the concentration to rise from 10 to 90% of the steady-state value. Figure 5 illustrates the time taken for the distribution of glucose in tumor design 1 (Figure 5B) is much greater than that in tumor design 2 (Figure 5A). These findings support the dye distribution patterns observed in these two configurations reported in Chen et al. Supplementary Video S3 shows the distribution of glucose computed by our discrete model for both the configurations of tumor vasculature.

FIGURE 5. Comparison of rise time tr for two different configurations of tumor vasculature. (A) Tumor design 2: the outlet is positioned away from the inlet, and the rise time is shorter in this configuration. (B) Tumor design 1: the outlet is positioned closer to the inlet, and the rise time is longer in this configuration.

### 3.3 Functional Coupling of Blood Vessel–Cell Exchange

In the mathematical framework presented in Section 2, it is practical to model heterogenous cell types exhibiting heterogenous enzyme activity at different vascular and cell densities observed under various pathophysiological conditions.

In this study, for the ease of demonstration of our method, we present the results of a minimal model of cellular glucose metabolism with uniform enzyme activity and homogenous cell type exchanging biochemicals with bloodstream. In Figure 6, we present the volumetric spatial analysis of the islet vasculature enveloped by a layer of homogenous cell mass, which forms the annular region of the tissue domain ωt. The extracellular glucose and lactate concentrations in Ωbv are initialized to 5 and 1.2 mM at the network inlet, respectively. The intracellular concentrations of glucose and lactate in ωt are initialized to zero. Due to the high concentration of glucose in the bloodstream, the gradient established between Ωbv and ωt in response to the advective–dispersive transport promotes the uptake of glucose by the annular region ωt. The cellular enzymes metabolize glucose-to-lactate, and the concentration of lactate in ωt rises above the basal value in the bloodstream. This gradient drives the export of lactate into the bloodstream until the system equilibrates.

FIGURE 6. Spatial concentrations of metabolites in the blood vessel and tissue domains of the islet vasculature observed at three different instants in time (0.6, 6, and 30 s). (A) (eglucose)/ (B) (elactate) and (C) (glucose)/ (D) (lactate) denote the concentration of glucose/lactate species observed in blood vessel Ωbv and tissue ωt compartments, respectively. In panel (A) and (D), the concentration of cylindrical and spherical volume elements embedded in the blood vessel domain is simulated, and the gradient displayed along the length of the blood vessel is generated by interpolation. In panels (B) and (C), the color-coded spheres represent the concentration evaluated in the annular region of ωt.

It is known from clinical observations that vascular phenotypes alter in cohorts with disease conditions such as diabetes and tumor. For example, in a tumor condition, the glycolytic enzymes are upregulated to produce more energy molecules that aid in the rapid proliferation of tumor cells (Rojas et al., 2018; St Clair et al., 2018). As a result, glucose–lactate dynamics is perturbed when compared to normal cells. In case of diabetes, histopathological and image reconstruction studies reveal that a decrease in the cell mass and reduction in vascular density can alter insulin release patterns (Richards et al., 2010; Cohrs et al., 2017). The framework proposed here can be useful in such clinical applications for carrying out systematic analysis in which 1) the change in vascular density can be induced in the form of mutations in the network (e.g., by deleting or inserting blood vessels) to examine the effect of anatomical changes on functional response of a tissue; 2) the alterations in the expression levels of enzymes quantified as fold changes in proteomics studies (Haythorne et al., 2019; Malinowski et al., 2020) can be mapped to the reaction velocities (i.e., Vm, a function of enzyme abundance, can be scaled using fold change) to examine the effect of genetic or environmental perturbations on cellular dynamics.

### 3.4 Sensitivity of Concentration Dynamics to ΔP and Glucose Dose

The pressure conditions observed in a vascular tissue may vary due to several physiological factors. To study the influence of the pressure gradient on the time taken to reach the steady-state concentration, we vary ΔP across the vasculature by specifying the inlet pressure and zero outlet pressure. Figure 7A illustrates the concentration profiles generated by simulating the advection–dispersion dynamics of glucose species by varying the pressure drop from 20 Pa to 200 Pa. When ΔP is high, the velocity of fluid is high in each branch. Consequently, convective flow dominates over dispersion, and this results in short rise times. The effect of change in the pressure gradient results in change in transit time of the fluid from 90.76 s at ΔP = 20 Pa to 9.08 s at ΔP = 200Pa. In Figure 7B, we show the sensitivity of the net glucose uptake flux and net lactate release flux to different glucose doses set at the inlet. It is observed that with an increase in the glucose dose the tissue units uptake more glucose from blood, and this drives the formation of lactate in the tissue subdomains. The excess lactate is then transported to the vessel subdomains until equilibrium is attained and the driving potential is zero.

FIGURE 7. Sensitivity of concentration dynamics to ΔP and glucose dose. (A) Left: influence of the varying pressure gradient across the network on metabolite rise times. Glucose concentration observed at positions 12.07 μm (red), 54.85 μm (green), and 110.27 μm (blue) from the inlet node. Inset displays the distribution of the Peclet number obtained by computing ul/

$D̃$

for each branch. Dashed and solid lines indicate the dynamic change in concentration observed at 20Pa and 200Pa, respectively. (B) Right: variation in the net glucose uptake (red) and net lactate release (blue) from the cells in response to change in glucose dose is displayed.

In conclusion, here we have presented a mathematical framework for understanding the multiscale connectivity existing in the functional networks of tissues. The test cases presented previously demonstrate how experimental data from different sources (e.g., kinetic data available in databases such as SABIO-RK (Wittig et al., 2012) and BRENDA (Chang et al., 2021), proteomics and metabolomics data, imaging data of vascular phenotypes, and flux measurements from perfusion experiments) can be encoded in our framework to build explicit models of cellular and vessel-to-cell interaction dynamics to better predict pathomechanisms. At the intra-organ scale, this work can be further extended to include cell-to-cell communication networks and decipher the order of communication occurring in the microenvironments with different vasculature architectures (e.g., periphery to the center, center to periphery, and one pole to other patterns in islets (El-Gohary and Gittes, 2018)) and cytoarchitectures (e.g., mantle-core and heterogenous distribution of β and α cells in islets (Dolenšek et al., 2015)). The metamodel is easily mutable; it is possible to induce different vascular phenotypes and disease states by altering vascular density and incorporating fold changes of metabolic and enzyme concentrations in the cellular units. To further research efforts involved in carrying out virtual experiments of inter-organ communication, our framework can be utilized for piecing together the top-down and bottom-up modeling approaches and accommodates each organ at the desired spatial resolution. For example, a comprehensive understanding of multi-organ disease states such as diabetes can be developed by interpreting the intra- and inter-organ interaction as communication occurring within “network of networks”; in whole-body models, the inter-organ communication can be modeled by abstracting organs as compartments forming nodes of the global network, and the detailed local dynamics occurring in the functional networks of an organ can be modeled by including subnetworks in the global network.

## Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

## Author Contributions

DM, SR, and DP conceived and designed the study. DM implemented the code, performed the analysis, and drafted the manuscript. All authors read, discussed the results, revised, and approved the manuscript.

## Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## Acknowledgments

DM thanks Matthias König, Institute for Theoretical Biology, Humboldt-University Berlin, for the valuable discussions, and Pradip Dutta, Department of Mechanical Engineering, Indian Institute of Science, for offering the computational facility to carry out COMSOL simulations.