A Graph-Based Framework for Multiscale Modeling of Physiological Transport

Trillions of chemical reactions occur in the human body every second, where the generated products are not only consumed locally but also transported to various locations in a systematic manner to sustain homeostasis. Current solutions to model these biological phenomena are restricted in computability and scalability due to the use of continuum approaches in which it is practically impossible to encapsulate the complexity of the physiological processes occurring at diverse scales. Here, we present a discrete modeling framework defined on an interacting graph that offers the flexibility to model multiscale systems by translating the physical space into a metamodel. We discretize the graph-based metamodel into functional units composed of well-mixed volumes with vascular and cellular subdomains; the operators defined over these volumes define the transport dynamics. We predict glucose drift governed by advective–dispersive transport in the vascular subdomains of an islet vasculature and cross-validate the flow and concentration fields with finite-element–based COMSOL simulations. Vascular and cellular subdomains are coupled to model the nutrient exchange occurring in response to the gradient arising out of reaction and perfusion dynamics. The application of our framework for modeling biologically relevant test systems shows how our approach can assimilate both multi-omics data from in vitro–in vivo studies and vascular topology from imaging studies for examining the structure–function relationship of complex vasculatures. The framework can advance simulation of whole-body networks at user-defined levels and is expected to find major use in personalized medicine and drug discovery.

GRAPHICAL ABSTRACT

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

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 ($\mathbb{W}$) of the graph network ($\mathcal{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) ($\mathbb{V}$) positioned at the N-furcation points or terminal ends and edges ($\mathbb{E}$) formed by the vessel segment linking two nodes. 4) Segment length l^bv and diameter d^bv 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, $\mathcal{G}(\mathbb{V}\in {\mathbb{R}}^m,\mathbb{E}\in {\mathbb{R}}^n,\mathbb{W}\in {\mathbb{R}}^n)$, for investigating the advection–diffusion-reaction dynamics. The cardinality of all the properties and the operators defined on $\mathcal{G}$ are summarized in Supplementary Table S1.

FIGURE 2

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

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

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 l^bv and circular radius r^bv derived from the diameter d^bv 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 r^bv, and length l^bv. 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:

$Q_{ij}=\frac{1}{R_{ij}}\left(p_i-p_j\right).(1)$

Here, p_i and p_j are the pressures at tail i and head j nodes of the oriented edge e_ij, and G_ij is the conductance associated with e_ij.

$Q_{ij}=G_{ij}\left(p_i-p_j\right)=\frac{\pi {\left(r^{bv}\right)}^4}{8\mu l^{bv}}\mathrm{\Delta }p.(2)$

The vector $Q\in {\mathbb{R}}^n$ of volumetric flow rates in n edges is defined in Eq. 3. The negative of the oriented incidence matrix $M(\mathcal{G})\in {\mathbb{R}}^{n\times m}$, denoted as M henceforth, is the gradient operator that acts on the vector $P\in {\mathbb{R}}^m$ of nodal pressures to result in the vector $\mathrm{\Delta }P\in {\mathbb{R}}^n$ of pressure gradients. We obtain Q by premultiplying ΔP with the diagonal conductance-matrix $G\in {\mathbb{R}}^{n\times n}$; this scales the pressure gradient across each edge by the corresponding edge conductance.

$Q=-GMP.(3)$

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, Q_ij is positive when flow enters node i $\left(Q_{i\leftarrow j}\right)$, negative when flow leaves node i $\left(Q_{i\to j}\right)$, and $\mathcal{A}(i)$ denotes the set of nodes that are adjacent to i.

$\sum _{j\in \mathcal{A}\left(i\right)}{Q}_{ij}=0.(4)$

Over m nodes in the network, the vector M^TQ ($M^T\in {\mathbb{R}}^{m\times 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 $\tilde{q}\in {\mathbb{R}}^m$ (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 p_in or outlet p_out. The values of these boundary conditions were specified based on experimental measurements of blood flow velocities reported in Table 1.

${\tilde{q}}_i=\left\{\begin{array}{cc}-Q_{in},\hfill & \text{if}i=\text{inlet\hspace{0.17em}node}\hfill \\ Q_{out},\hfill & \text{if}i=\text{outlet\hspace{0.17em}node}\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array}\right.(5)$

$M^{T}GMP=\tilde{q}.(6)$

Here, M^TGM is a square matrix. After substitution of known pressures in the vector P, the columns of M^TGM 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 u_ij 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 ∈ (c_l + δ, c_l—δ). We derive the characteristic length c_l 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 c_l.

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 $V_i^{bv}$ (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 $\in \mathcal{A}(i)$ (Reichold et al., 2009). Here, $\mathcal{A}(i)$ is the set {i-1, i+1, i+2} of nodes adjacent to i, and A_ij and l_ij 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.

$V_i^{bv}=\frac{1}{2}\sum _{j\in \mathcal{A}\left(i\right)}{A}_{ij}{l}_{ij}.(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\frac{\partial C^{bv}}{\partial t}=\nabla .\left(-Q{C}^{bv}+DA\nabla C^{bv}\right).(8)$

Here, A denotes the area of flow, D is the dispersion coefficient of species, Q denotes the flow field, and C^bv is the scalar concentration of a chemical species. To shift from the continuous to discrete counterpart, first we assign the scalar concentration field C₀ to all vertices, and the flow field and dispersion coefficients form the edge weights. The net change in the molar concentration of a species $C_i^{bv}$ in the control volume $V_i^{bv}$ at the ith node of the blood vessel depends on the contributions from advective $J_i^A$ and dispersive fluxes $J_i^D$ presented in Equation 9. The influx and outflux of $J_i^A$ are the summation of molar fluxes through the edges feeding into i (j → i) and leaving i (i → l), 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, $\mathcal{A}(i)$ is the set that contains nodes that are adjacent to node i, ${\mathcal{A}}_{+}$(i) is the set of in-neighbor nodes of i, ${\mathcal{A}}_{-}$(i) is the set of out-neighbor nodes of i, and D_ik, A_ik, and l_ik are the dispersion coefficient, cross-sectional area, and length of edge e_ik, respectively.

$\begin{array}{c}\hfill V_i^{bv}\frac{d{C}_i^{bv}}{dt}=\underset{\text{advection\hspace{0.17em}flux}{J}_i^A}{\underbrace{\underset{\mathit{\text{influx}}}{\underbrace{\sum _{j\in {\mathcal{A}}_{+}\left(i\right)}{C}_j^{bv}{Q}_{j\to i}^{in}}}-\underset{\mathit{\text{outflux}}}{\underbrace{\sum _{l\in {\mathcal{A}}_{-}\left(i\right)}{C}_i^{bv}{Q}_{i\to l}^{out}}}}}\\ \hfill +\underset{\text{dispersion\hspace{0.17em}flux}{J}_i^D}{\underbrace{\sum _{k\in \mathcal{A}\left(i\right)}\underset{D_{ik}^Q}{\underbrace{\frac{D_{ik}{A}_{ik}}{l_{ik}}}}\left(C_k^{bv}-C_i^{bv}\right)}}\end{array}(9)$

$D_{ik}={\tilde{D}}_{ik}+\frac{u_{ik}^2{r}_{ik}^2}{48{\tilde{D}}_{ik}}.(10)$

The dispersion coefficient is determined from Aris–Taylor’s relation in Eq. 10. Here, ${\tilde{D}}_{ik}$, u_ik, and r_ik are the diffusion coefficient of a species, centerline velocity, and radius of e_ik, respectively.

$V^{bv}\frac{d{C}^{bv}}{dt}=M^{T}Q{M}^o{C}^{bv}-M^T{D}^{Q}M{C}^{bv}.(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, $D^Q\in {\mathbb{R}}^{n\times n}$is the diagonal matrix with diagonal entries the volumetric dispersion coefficient DA/l of each edge specified in Eq. 9, $Q\in {\mathbb{R}}^{n\times n}$ is the diagonal matrix with diagonal entries the volumetric flow rate of each edge, $M^o\in {\mathbb{R}}^{n\times m}$ is the modified incidence matrix (Rak, 2017; Erlich et al., 2019), M^TD^QM is the weighted dispersion Laplacian matrix $L^D\left(\mathcal{G}\right)\in {\mathbb{R}}^{m\times m}$, and M^TQM^o is the weighted advection Laplacian matrix $L^Q\left(\mathcal{G}\right)\in {\mathbb{R}}^{m\times m}$ (Eq. 12).

$V^{bv}\frac{d{C}^{bv}}{dt}=L^Q{C}^{bv}-L^D{C}^{bv}.(12)$

2.5 Advection–Dispersion-Reaction

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 d_bv and d_t. A value of 12.4 µm was used in our model for d_t. This value is derived from the physical volume of an insulin-secreting beta cell (1020 μm³) (Finegood et al., 1995) located in the islets of pancreas.

2.5.2 Equations

$\begin{array}{cc}\hfill V_i^{bv}\frac{d{C}_i^{bv}}{dt}& =\underset{\text{advection\hspace{0.17em}flux}{J}_i^A}{\underbrace{\sum _{j\in {\mathcal{A}}_{+}\left(i\right)}{C}_j^{bv}{Q}_{ji}^{in}-\sum _{l\in {\mathcal{A}}_{-}\left(i\right)}{C}_i^{bv}{Q}_{il}^{out}}}\hfill \\ \hfill & +\underset{\text{dispersion\hspace{0.17em}flux}{J}_i^D}{\underbrace{\sum _{k\in \mathcal{A}\left(i\right)}\frac{D_{ik}{A}_{ik}}{l_{ik}}\left(C_k^{bv}-C_i^{bv}\right)}}\hfill \\ \hfill & -\underset{\text{exchange\hspace{0.17em}flux}{J}_i^E}{\underbrace{V_i^{cell}{j}_i^E}}.\hfill \end{array}(13)$

In the last term of Eq. 13, $j_i^E$ (Eq. 15) is the net exchange rate that governs the bidirectional transport between tissue compartment ${\omega }_i^t$ and blood vessel compartment ${\mathrm{\Omega }}_i^{bv}$ of the functional unit. The uptake or release flux $J_i^E$ is computed by multiplying $j_i^E$ with $V_i^{cell}$ since the maximal rates are often reported in per unit volume of a biological cell. This ensures mole balance when a species moves from ${\omega }_i^t$ to ${\mathrm{\Omega }}_i^{bv}$ which differs in volume. Furthermore, $C_i^{cell}$ is the concentration of species in the tissue cell that interacts with ith node of Ω^bv, the half-saturation constant K_m quantifies the affinity of a transporter protein or an enzyme for a metabolite, and V_m is the maximal rate of metabolite transport. When ${\mathrm{\Omega }}_i^{bv}$ is encompassed by ${\omega }_i^t$, interaction between the node associated with both the compartments exists, and $J_i^E$ appears as a source or sink term in Eq. 13. At the junction nodes of the blood vessel where no interaction exists, $J_i^E$ is zero.

$V^{bv}\frac{d{C}^{bv}}{dt}=M^{T}Q{M}^o{C}^{bv}-M^T{D}^{Q}M{C}^{bv}-J^E(14)$

$j_i^E=\frac{V_m}{K_m}\left(\frac{C_i^{bv}-C_i^{cell}}{1+\frac{C_i^{cell}}{K_m}+\frac{C_i^{bv}}{K_m}}\right).(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 C^bv—C^cell across the vessel wall.

2.6 Cell

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

$V^{cell}\frac{d{C}_i^{cell}}{dt}=Nv,(16)$

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, K_m and V_m, were chosen from panmin (König and Deepa Maheshvare, 2021), a minimal model of glucose metabolism and insulin secretion in the pancreatic β-cell. The V_m 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 V_m 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

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 M^TQ, and the visualization is presented in Figure 4C.

FIGURE 4

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 t_r 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 ${\mathrm{\Omega }}_i^{bv}$ 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

FIGURE 5. Comparison of rise time t_r 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

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

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/$\tilde{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.

Funding

The study is supported by intramural funds provided by the Ministry of Education, Government of India. Research of DM is supported by the Senior Research Fellowship from the Ministry of Human Resource Development (MHRD), Government of India. The preprint of this article has appeared online on bioRxiv (Deepa Maheshvare et al., 2021).

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.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnetp.2021.802881/full#supplementary-material

References

Antiga, L., Piccinelli, M., Botti, L., Ene-Iordache, B., Remuzzi, A., and Steinman, D. A. (2008). An Image-Based Modeling Framework for Patient-specific Computational Hemodynamics. Med. Biol. Eng. Comput. 46 (11), 1097–1112. doi:10.1007/s11517-008-0420-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Barua, A. K., and Goel, P. (2016). Isles within Islets: the Lattice Origin of Small-World Networks in Pancreatic Tissues. Physica D: Nonlinear Phenomena 315, 49–57. doi:10.1016/j.physd.2015.07.009

CrossRef Full Text | Google Scholar

Beard, D. A., and Bassingthwaighte, J. B. (2000). Advection and Diffusion of Substances in Biological Tissues with Complex Vascular Networks. Ann. Biomed. Eng. 28 (3), 253–268. doi:10.1114/1.273

PubMed Abstract | CrossRef Full Text | Google Scholar

Bellocchi, L., and Geroliminis, N. (2020). Unraveling Reaction-diffusion-like Dynamics in Urban Congestion Propagation: Insights from a Large-Scale Road Network. Sci. Rep. 10 (1), 4876–4911. doi:10.1038/s41598-020-61486-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Berndt, N., Bulik, S., Wallach, I., Wünsch, T., König, M., Stockmann, M., et al. (2018a). Hepatokin1 Is a Biochemistry-Based Model of Liver Metabolism for Applications in Medicine and Pharmacology. Nat. Commun. 9 (1), 2386–2412. doi:10.1038/s41467-018-04720-9

CrossRef Full Text | Google Scholar

Berndt, N., Horger, M. S., Bulik, S., and Holzhütter, H.-G. (2018b). A Multiscale Modelling Approach to Assess the Impact of Metabolic Zonation and Microperfusion on the Hepatic Carbohydrate Metabolism. Plos Comput. Biol. 14 (2), e1006005. doi:10.1371/journal.pcbi.1006005

PubMed Abstract | CrossRef Full Text | Google Scholar

Besse, Christophe., and Faye, Grégory. (2021). Dynamics of Epidemic Spreading on Connected Graphs. J. Math. Biol. 82 (6), 1–52. doi:10.1007/s00285-021-01602-5

CrossRef Full Text | Google Scholar

Chalhoub, E., Xie, L., Balasubramanian, V., Kim, J., and Belovich, J. (2007). A Distributed Model of Carbohydrate Transport and Metabolism in the Liver during Rest and High-Intensity Exercise. Ann. Biomed. Eng. 35 (3), 474–491. doi:10.1007/s10439-006-9217-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, A., Jeske, L., Ulbrich, S., Hofmann, J., Koblitz, J., Schomburg, I., et al. (2021). Brenda, the Elixir Core Data Resource in 2021: New Developments and Updates. Nucleic Acids Res. 49 (D1), D498–D508. doi:10.1093/nar/gkaa1025

PubMed Abstract | CrossRef Full Text | Google Scholar

Chapman, Airlie. (2015), Advection on Graphs. Semi-Autonomous Networks, 3–16. Springer. doi:10.1007/978-3-319-15010-9_1

PubMed Abstract | CrossRef Full Text | Google Scholar

Cohrs, C. M., Chen, C., Jahn, S. R., Stertmann, J., Chmelova, H., Weitz, J., et al. (2017). Vessel Network Architecture of Adult Human Islets Promotes Distinct Cell-Cell Interactions In Situ and Is Altered after Transplantation. Endocrinology 158 (5), 1373–1385. doi:10.1210/en.2016-1184

PubMed Abstract | CrossRef Full Text | Google Scholar

Deepa Maheshvare, M., Raha, S., and Pal, D. (2021). A Graph-Based Framework for Multi-Scale Modeling of Physiological Transport. bioRxiv. doi:10.1101/2021.09.14.460337

CrossRef Full Text | Google Scholar

d’Esposito, Angela., Sweeney, Paul. W., Ali, Morium., Saleh, Magdy., Ramasawmy, Rajiv., Roberts, Thomas. A., et al. (2018). Computational Fluid Dynamics with Imaging of Cleared Tissue and of In Vivo Perfusion Predicts Drug Uptake and Treatment Responses in Tumours. Nat. Biomed. Eng. 2 (10), 773–787.

PubMed Abstract | Google Scholar

Deussen, A., and Bassingthwaighte, J. B. (1996). Modeling [15O]oxygen Tracer Data for Estimating Oxygen Consumption. Am. J. Physiology-Heart Circulatory Physiol. 270 (3), H1115–H1130. doi:10.1152/ajpheart.1996.270.3.h1115

PubMed Abstract | CrossRef Full Text | Google Scholar

Diez, J. A., Arrojo e Drigo, R., Zheng, X., Stelmashenko, O. V., Chua, M., Rodriguez-Diaz, R., et al. (2017). Pancreatic Islet Blood Flow Dynamics in Primates. Cel. Rep. 20 (6), 1490–1501. doi:10.1016/j.celrep.2017.07.039

PubMed Abstract | CrossRef Full Text | Google Scholar

Dolenšek, J., Rupnik, M. S., and Stožer, A. (2015). Structural Similarities and Differences between the Human and the Mouse Pancreas. Islets 7 (1), e1024405. doi:10.1080/19382014.2015.1024405

PubMed Abstract | CrossRef Full Text | Google Scholar

El-Gohary, Yousef., and George, Gittes. (2018). Structure of Islets and Vascular Relationship to the Exocrine Pancreas. Pancreapedia: The Exocrine Pancreas Knowledge Base.

Google Scholar

Erlich, A., Pearce, P., Mayo, R. P., Jensen, O. E., and Chernyavsky, I. L. (2019). Physical and Geometric Determinants of Transport in Fetoplacental Microvascular Networks. Sci. Adv. 5 (4), eaav6326. doi:10.1126/sciadv.aav6326

PubMed Abstract | CrossRef Full Text | Google Scholar

Fang, Q., Sakadžić, S., Ruvinskaya, L., Devor, A., Dale, A. M., and Boas, D. A. (2008). Oxygen Advection and Diffusion in a Three-Dimensional Vascular Anatomical Network. Opt. Express 16 (22), 17530. doi:10.1364/oe.16.017530

PubMed Abstract | CrossRef Full Text | Google Scholar

Fedorov, A., Beichel, R., Kalpathy-Cramer, J., Finet, J., Fillion-Robin, J.-C., Pujol, S., et al. (2012). 3d Slicer as an Image Computing Platform for the Quantitative Imaging Network. Magn. Reson. Imaging 30 (9), 1323–1341. doi:10.1016/j.mri.2012.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Finegood, D. T., Scaglia, L., and Bonner-Weir, S. (1995). Dynamics of -cell Mass in the Growing Rat Pancreas: Estimation with a Simple Mathematical Model. Diabetes 44 (3), 249–256. doi:10.2337/diab.44.3.249

PubMed Abstract | CrossRef Full Text | Google Scholar

Frost, Timothy. S., Estrada, Victor., Jiang, Linan., and Zohar, Yitshak. (2019). Convection–diffusion Molecular Transport in a Microfluidic Bilayer Device with a Porous Membrane. Microfluidics and Nanofluidics 23 (10), 1–13. doi:10.1007/s10404-019-2283-1

CrossRef Full Text | Google Scholar

Gabriel, E. M., Kim, M., Fisher, D. T., Powers, C., Attwood, K., Bagaria, S. P., et al. (2020). Dynamic Control of Tumor Vasculature Improves Antitumor Responses in a Regional Model of Melanoma. Sci. Rep. 10 (1), 13245–13313. doi:10.1038/s41598-020-70233-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Geuzaine, C., and Remacle, J.-F. (2009). Gmsh: A 3-D Finite Element Mesh Generator with Built-In Pre- and post-processing Facilities. Int. J. Numer. Meth. Engng. 79 (11), 1309–1331. doi:10.1002/nme.2579

CrossRef Full Text | Google Scholar

Golub, G. H., and Pereyra, V. (1973). The Differentiation of Pseudo-inverses and Nonlinear Least Squares Problems Whose Variables Separate. SIAM J. Numer. Anal. 10 (2), 413–432. doi:10.1137/0710036

CrossRef Full Text | Google Scholar

Gray, M. R., and Tam, Y. K. (1987). The Series-Compartment Model for Hepatic Elimination. Drug Metab. Dispos. 15 (1), 27–31.

PubMed Abstract | Google Scholar

Haythorne, E., Rohm, M., van de Bunt, M., Brereton, M. F., Tarasov, A. I., Blacker, T. S., et al. (2019). Diabetes Causes Marked Inhibition of Mitochondrial Metabolism in Pancreatic β-cells. Nat. Commun. 10 (1), 2474–2517. doi:10.1038/s41467-019-10189-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Heaton, L. L., López, E., Maini, P. K., Fricker, M. D., and Jones, N. S. (2012). Advection, Diffusion, and Delivery over a Network. Phys. Rev. E Stat. Nonlin Soft Matter Phys. 86 (2), 021905. doi:10.1103/PhysRevE.86.021905

PubMed Abstract | CrossRef Full Text | Google Scholar

Heaton, L. M. (2012). Biological Transport Networks. PhD thesis (Oxford: University of Oxford).

Hošek, R., and Volek, J. (2019). Discrete Advection–Diffusion Equations on Graphs: Maximum Principle and Finite Volumes. Appl. Math. Comput. 361, 630–644.

Google Scholar

Ii, K., Mashimo, K., Ozeki, M., Yamada, T. G., Hiroi, N., and Funahashi, A. (2019). Xitosbml: A Modeling Tool for Creating Spatial Systems Biology Markup Language Models from Microscopic Images. Front. Genet. 10, 1027. doi:10.3389/fgene.2019.01027

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, N., Cox, R. D., and Hancock, J. M. (2007). A Kinetic Core Model of the Glucose-Stimulated Insulin Secretion Network of Pancreatic Beta Cells. Mamm. Genome 18 (6), 508–520. doi:10.1007/s00335-007-9011-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Kinstlinger, I. S., Saxton, S. H., Calderon, G. A., Ruiz, K. V., Yalacki, D. R., Deme, P. R., et al. (2020). Generation of Model Tissues with Dendritic Vascular Networks via Sacrificial Laser-Sintered Carbohydrate Templates. Nat. Biomed. Eng. 4 (9), 916–932. doi:10.1038/s41551-020-0566-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Kirkegaard, J. B., and Sneppen, K. (2020). Optimal Transport Flows for Distributed Production Networks. Phys. Rev. Lett. 124 (20), 208101. doi:10.1103/physrevlett.124.208101

PubMed Abstract | CrossRef Full Text | Google Scholar

Koh, D.-S., Cho, J.-H., and Chen, L. (2012). Paracrine Interactions within Islets of Langerhans. J. Mol. Neurosci. 48 (2), 429–440. doi:10.1007/s12031-012-9752-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Kojic, M., Milosevic, M., Simic, V., Koay, E. J., Fleming, J. B., Nizzero, S., et al. (2017). A Composite Smeared Finite Element for Mass Transport in Capillary Systems and Biological Tissue. Comp. Methods Appl. Mech. Eng. 324, 413–437. doi:10.1016/j.cma.2017.06.019

CrossRef Full Text | Google Scholar

König, Matthias., Bulik, Sascha., and Holzhütter, Hermann-Georg. (2012). Quantifying the Contribution of the Liver to Glucose Homeostasis: a Detailed Kinetic Model of Human Hepatic Glucose Metabolism. Plos Comput. Biol. 8 (6), e1002577.

PubMed Abstract | Google Scholar

König, Matthias., and Deepa Maheshvare, M. (2021). pancreas_minimal: Minimal Model of Pancreas Glucose Metabolism and Insulin Secretion. November. doi:10.5281/zenodo.5729745

CrossRef Full Text | Google Scholar

Korsgren, E., and Korsgren, O. (2016). An Apparent Deficiency of Lymphatic Capillaries in the Islets of Langerhans in the Human Pancreas. Diabetes 65 (4), 1004–1008. doi:10.2337/db15-1285

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, D., Gupta, J., and Raha, S. (2019). Partitioning a Reaction-Diffusion Ecological Network for Dynamic Stability. Proc. R. Soc. A. 475 (2223), 20180524. doi:10.1098/rspa.2018.0524

PubMed Abstract | CrossRef Full Text | Google Scholar

Mac, L. (2020). Converting a Graph to a 2d Diagram.

Google Scholar

Malinowski, R. M., Ghiasi, S. M., Mandrup-Poulsen, T., Meier, S., Lerche, M. H., Ardenkjær-Larsen, J. H., et al. (2020). Pancreatic β-cells Respond to Fuel Pressure with an Early Metabolic Switch. Sci. Rep. 10 (1), 15413–15511. doi:10.1038/s41598-020-72348-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Masid, M., Ataman, M., and Hatzimanikatis, V. (2020). Analysis of Human Metabolism by Reducing the Complexity of the Genome-Scale Models Using Redhuman. Nat. Commun. 11 (1), 2821–2912. doi:10.1038/s41467-020-16549-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Musy, M., Jacquenot, G., Dalmasso, G., neoglez,, , de Bruin, R., Pollack, A., et al. (2021). marcomusy/vedo: 2021.0.2. Zenodo. doi:10.5281/zenodo.4609336

CrossRef Full Text | Google Scholar

Parween, S., Kostromina, E., Nord, C., Eriksson, M., Lindström, P., and Ahlgren, U. (2016). Intra-islet Lesions and Lobular Variations in β-cell Mass Expansion in Ob/ob Mice Revealed by 3D Imaging of Intact Pancreas. Sci. Rep. 6 (1), 34885–34911. doi:10.1038/srep34885

PubMed Abstract | CrossRef Full Text | Google Scholar

Pedersen, M. G., Bertram, R., and Sherman, A. (2005). Intra- and Inter-islet Synchronization of Metabolically Driven Insulin Secretion. Biophysical J. 89 (1), 107–119. doi:10.1529/biophysj.104.055681

PubMed Abstract | CrossRef Full Text | Google Scholar

Pisania, A., Weir, G. C., O'Neil, J. J., Omer, A., Tchipashvili, V., Lei, J., et al. (2010). Quantitative Analysis of Cell Composition and Purity of Human Pancreatic Islet Preparations. Lab. Invest. 90 (11), 1661–1675. doi:10.1038/labinvest.2010.124

PubMed Abstract | CrossRef Full Text | Google Scholar

Poelma, C. (2017). Exploring the Potential of Blood Flow Network Data. Meccanica 52 (3), 489–502. doi:10.1007/s11012-015-0255-4

CrossRef Full Text | Google Scholar

Prentki, M., Matschinsky, F. M., and Madiraju, S. R. M. (2013). Metabolic Signaling in Fuel-Induced Insulin Secretion. Cel. Metab. 18 (2), 162–185. doi:10.1016/j.cmet.2013.05.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Pries, A. R., Secomb, T. W., Gaehtgens, P., and Gross, J. F. (1990). Blood Flow in Microvascular Networks. Experiments and Simulation. Circ. Res. 67 (4), 826–834. doi:10.1161/01.res.67.4.826

PubMed Abstract | CrossRef Full Text | Google Scholar

Rak, Annie. (2017). Advection on Graphs. PhD thesis.

Reichold, J., Stampanoni, M., Keller, A. L., Buck, A., Jenny, P., and Weber, B. (2009). Vascular Graph Model to Simulate the Cerebral Blood Flow in Realistic Vascular Networks. J. Cereb. Blood Flow Metab. 29 (8), 1429–1443. doi:10.1038/jcbfm.2009.58

CrossRef Full Text | Google Scholar

Richards, O. C., Raines, S. M., and Attie, A. D. (2010). The Role of Blood Vessels, Endothelial Cells, and Vascular Pericytes in Insulin Secretion and Peripheral Insulin Action. Endocr. Rev. 31 (3), 343–363. doi:10.1210/er.2009-0035

PubMed Abstract | CrossRef Full Text | Google Scholar

Rojas, J. D., Papadopoulou, V., Czernuszewicz, T. J., Rajamahendiran, R. M., Chytil, A., Chiang, Y. C., et al. (2018). Ultrasound Measurement of Vascular Density to Evaluate Response to Anti-angiogenic Therapy in Renal Cell Carcinoma. IEEE Trans. Biomed. Eng. 66 (3), 873–880. doi:10.1109/TBME.2018.2860932

PubMed Abstract | CrossRef Full Text | Google Scholar

Safaei, S., Blanco, P. J., Müller, L. O., Hellevik, L. R., and Hunter, P. J. (2018). Bond Graph Model of Cerebral Circulation: toward Clinically Feasible Systemic Blood Flow Simulations. Front. Physiol. 9, 148. doi:10.3389/fphys.2018.00148

PubMed Abstract | CrossRef Full Text | Google Scholar

Sommer, C., Straehle, C., Kothe, U., and Hamprecht, F. A. (2011). “Ilastik: Interactive Learning and Segmentation Toolkit,” in 2011 IEEE International Symposium on Biomedical Imaging, From Nano to Macro, 230–233. doi:10.1109/ISBI.2011.5872394

CrossRef Full Text | Google Scholar

Sorensen, J. T. (1985). A Physiologic Model of Glucose Metabolism in Man and its Use to Design and Assess Improved Insulin Therapies for Diabetes. PhD thesis (Cambridge: Massachusetts Institute of Technology).

St Clair, J. R., Ramirez, D., Passman, S., and Benninger, R. K. P. (2018). Contrast-enhanced Ultrasound Measurement of Pancreatic Blood Flow Dynamics Predicts Type 1 Diabetes Progression in Preclinical Models. Nat. Commun. 9 (1), 1742–1812. doi:10.1038/s41467-018-03953-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Taylor, Geoffrey. Ingram. (1954). The Dispersion of Matter in Turbulent Flow through a Pipe. Proc. R. Soc. Lond. Ser. A. Math. Phys. Sci. 223 (1155), 446–468. doi:10.1098/rspa.1954.0130

CrossRef Full Text | Google Scholar

Thurber, G. M., and Weissleder, R. (2011). A Systems Approach for Tumor Pharmacokinetics. PloS One 6 (9), e24696. doi:10.1371/journal.pone.0024696

PubMed Abstract | CrossRef Full Text | Google Scholar

Rao, V. P., and Rizzo, M. A. Diffusion of Metabolites across gap Junctions Mediates Metabolic Coordination of β-islet Cells. bioRxiv, 2020.

Google Scholar

Watts, M., Ha, J., Kimchi, O., and Sherman, A. (2016). Paracrine Regulation of Glucagon Secretion: the β/α/δ Model. Am. J. Physiology-Endocrinology Metab. 310 (8), E597–E611. doi:10.1152/ajpendo.00415.2015

CrossRef Full Text | Google Scholar

Wittig, U., Kania, R., Golebiewski, M., Rey, M., Shi, L., Jong, L., et al. (2012). SABIO-RK--database for Biochemical Reaction Kinetics. Nucleic Acids Res. 40 (D1), D790–D796. doi:10.1093/nar/gkr1046

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, J., Davis, T., Kazerouni, A. S., Chen, Y-I., Bloom, M. J., Yeh, H-C., et al. (2021). Longitudinal Fret Imaging of Glucose and Lactate Dynamics and Response to Therapy in Breast Cancer Cells. Mol. Imaging Biol., 1–12. doi:10.1007/s11307-021-01639-4

CrossRef Full Text | Google Scholar

Chen, Y. Y., and Kingston, B. R., and Chan, W. C. W. (2020). Transcribing In Vivo Blood Vessel Networks into In Vitro Perfusable Microfluidic Devices. Adv. Mater. Technol., 5 (6), 2000103. doi:10.1002/admt.202000103

CrossRef Full Text | Google Scholar