Network Modeling of Murine Lymphatic System

: Animal models of diseases, particularly mice, are considered to be the cornerstone for translational research in immunology. The aim of the present study is to model the geometry and analyze the network structure of the murine lymphatic system (LS). The algorithm for building the graph model of the LS makes use of anatomical data. To identify the edge directions of the graph model, a mass balance approach to lymph dynamics based on the Hagen–Poiseuille equation is applied. It is the ﬁrst study in which a geometric model of the murine LS has been developed and characterized in terms of its structural organization and the lymph transfer function. Our study meets the demand for quantitative mechanistic approaches in the growing ﬁeld of immunoengineering to utilize or exploit the lymphatic system for immunotherapy.


Introduction
Networks of various natures, e.g., structural, functional, spatial, underlie the dynamics and mechanisms of regulation in live systems ranging from cells to physiological organs and to whole organisms [1].Network concepts are increasingly used to describe the structure and function of the immune system [2].The immune system represents an example of a highly complex network of interacting and migration cell populations embedded into the spatially distributed components of the lymphatic system (LS) [3].
The LS is a body-wide network of lymphatic vessels and lymphoid organs with two major functions: (1) fluid transport from tissues to the blood system to maintain fluid homeostasis and (2) trafficking of antigens and immune cells to lymph nodes where the immune responses take place [4,5].Lymphoid organs include a large number of lymph nodes as well as the spleen, thymus, tonsils, and bone marrow [5].The lymphatic vessels are the conduits that facilitate the directional lymph transport from peripheral tissues to secondary lymph nodes [4].Understanding of the lymphatic structural and functional organization is essential to discern how the LS interacts with different tissues and organs within the body [6].
Animal models of diseases, particularly mice, are considered to be the cornerstone for translational research in immunology [7,8].Research with laboratory mice enabled invaluable insight into mammalian immune systems [9].Despite numerous advances in understanding the immune system from mouse studies, there exist fundamental differences between mouse and human immune systems [8].The structural organization of the lymphatic system represents a straightforward example.However, a comprehensive mathematical network-based characterization of the LS in mice is still unavailable.
The primary function of the lymphatic system is the maintenance of the interstitial fluid homeostasis [10].The structure and topology of the LS network is heterogeneous and remains to be systematically explored [5].The existing mathematical models of the lymphatic system refer to specific parts of it, such as the lymphatic capillary network [11], collecting lymphatics [12], lymphangions [13], or branching networks of lymphatic vessels [14].Computational models of the whole lymphatic system network are still rare [15][16][17][18].Comprehensive reviews of existing approaches to modeling the lymphatic system structure and function can be found in [10,19].One of the major bottlenecks in developing the computational models of the lymphatic system is due to the lack of comprehensive anatomical and physiological data [10].The latest research activity has clearly stated, "Thus, gross lymphatic anatomy has not been updated for more than a century.... our knowledge of macro-lymphatic anatomy remains rudimentary" [6].The existing gaps [20], require further systematic research [21] including mathematical modeling [22], which serves to integrate available knowledge and to identify critical issues amenable to further biological testing.
The application of graph theory methods to describe the spatial organization of the human LS has been addressed in a number of recent studies (see for a review [18]).In [18], we developed a computational algorithm for representing the anatomy-based and rulebased graphs of the LS in humans.The graph models enabled the analysis of different metrics of complexity of the human LS such as spectral radius, clustering coefficient, average path length, and number of separators.A similar analysis for the mouse LS remains to be performed.
The aim of the present study is to model and analyze the network structure of the murine lymphatic system.The algorithm for building the graph model makes use of anatomical data.To identify the edge directions of the graph model, a mass balance approach to lymph dynamics based on the Hagen-Poiseuille equation is applied.Various matrix forms for graph representation are specified.The lymph transfer times between various nodes are estimated.We summarize the properties of the graph model of the murine LS using metrics similar to the human LS graph thus providing a quantitative basis for understanding essential structural differences of the LS between the mice and humans.

Anatomy and Physiology of Murine Lymphatic System
Available anatomical and physiological data provide the empirical basis for specifying the network structure of the lymphatic system in mice [23,24] using the notion of a simple graph.A simple graph G = (V, E) is a pair of sets V and E, with elements of V being vertices or nodes and E being edges [25].The simple graph with edges oriented in only one direction is called an oriented graph.There are some variations in the number of lymph nodes, i.e., ranging from 22 to 36 as indicated in Table 1.A generalized graph of the murine LS, consisting of 88 nodes and 87 edges is shown in Figure 1.It was developed using anatomical descriptions from [23,24].The vertices of the graph refer to either the lymph nodes, outlet vertices with out-degree deg + = 0 corresponding to the sink into jugular veins, the confluences of lymphatic vessels, or inlet vertices with in-degree deg − = 0 corresponding to the collecting lymphatics of various body tissues.
The geometric characteristics of the lymphatic vessels and the baseline parameters of the lymph flow through various parts of the LS network are detailed in Table 1.To set the pressure at the sink nodes p out , we used the estimate of the murine central venous pressure: 7.4 (5.9-8.9)cm H 2 O [26,27].Lymph viscosity is taken to be equal to 1.81 mPa•s.
The anatomy data enable specifying a simple graph of the murine LS.The adjacency matrix A of the LS graph is shown in Figure 2B.
As the LS functions to transport the lymph from the drained tissues to the venous part of cardiovascular system, additional analysis is required to transfer the simple graph representation to a physiologically meaningful oriented graph of the LS.To generate an oriented graph mode of the LS, we used a combination of experimental stud-ies on fluid dynamics in various parts of the LS in mice [28][29][30][31] and computation of the lymph flow through the system in accordance to an overall mass balance using the Hagen-Poiseuille equations.(BALB/cAnNCr) [24] (DD/NIH) [23] Diameter 1-2.3 mm 1-17.3 mm (C57Bl/6J, Nude, CB-17 SCID) [31,32] (DD/NIH) [23] Thoracic duct: Radius 300 µm [28] Flow

Oriented Graph Model of Murine LS
The graph of the lymphatic system G can be divided in two disconnected subgraphs: collecting the lymph to the left common jugular vein (g l ) and to the right common jugular vein (g r ) (Figures 1 and 3).The thoracic duct belongs to the subgraph g l .

Computing the Direction of Lymph Flows
As we aim to reproduce the target flow through thoracic duct q td = 10 mL/day (Table 1), we compute the flows in the left subgraph first.The pressure at the sink vertex is set to p out = 7.4 cm H 2 O, and the outflow from the left LS subgraph to the left jugular vein q (l) out is calibrated so that the computed flow through thoracic duct (the edge incident to jugular vein) is equal to q td .The inflows in the collecting lymphatics (inlet vertices with zero in-degree) are given to be the same and equal to q in is the number of inlet vertices.After obtaining the flows in the left subgraph g l , we compute the flows in the right subgraph g r by setting the same inflows q (r) in as in the left one, and the same output pressure p out .The outlet outflow is given by q (r) in .The distribution of the steady flows in the graph g(n, m) with n vertices and m edges is considered to be governed by the following:

•
The Hagen-Poiseuille equation links the flow q ij through the edge e ij with the drop of pressure from the tail i to the head j vertices (p i − p j ) by the conductances g ij that depend on lymph viscosity µ and the radii and the lengths of the edges; • The balance of flow through the vertices due to mass conservation: where A(i) is a set of vertices adjacent to i.
Using the oriented incidence matrix M ∈ R n×m and the diagonal conductance matrix G ∈ R m×m (with elements indexed by edges rather than nodes), one can rewrite Equations ( 1) and (2) as where p ∈ R n are the nodal pressures, q ∈ R n are the net flows through the vertices, and q ∈ R m are the flows through the edges.Hence, we get the linear system to solve for the nodal pressures: where L = MGM T is a symmetric weighted Laplacian matrix.
As at the outlet vertex the pressure is known (p out ), we substitute the vector p = [0, . . ., p out , . . ., 0] T of zeros with the known pressure at the corresponding index in (4).By shifting L p to the right-hand side, we obtain the system where L rect is the matrix L without the column corresponding to the index of the outlet vertex with known pressure.The pseudo-inverse of L rect provides the vector of unknown pressures: p unknown = L + rect ( q − L p).The graph was constructed and visualized using the R package igraph.The algorithm for computing the flows described in Section 3.1 was implemented in R using the ginv() function for pseudo-inverse calculation from the MASS package.In addition, we have verified the computation of the flows in Julia language using the singular value decomposition to obtain the pseudo-inverse.
An oriented graph of the murine LS resulting from the analysis of the global lymphatic flow balance is shown in Figure 1.It is derived assuming constant diameter of vessels in the LS.

Matrix Representation
To visualize the graph structure of the murine LS, we use the adjacency matrix, which indicates whether a pair of nodes are adjacent or not.For the oriented graph, the adjacency matrix is shown in Figure 2A.The adjacency matrix of the simple graph presented in Figure 2B is symmetric.A complementary representation of the graph is provided by the incidence matrix (see Figure 2C).The incidence matrix is different from an adjacency matrix, and it encodes the relation of node-vertex pairs.Finally, the graph Laplacian matrix is displayed in Figure 2D.It is related to the degree matrix D and the adjacency matrix A of the graph L = D − A, representing an edge-weighted graph.

Quantitative Characterization of Lymph Flows through the LS
The estimated values of the lymph flow through various vessels of the murine LS are specified in Figure 3.The upper panel shows the flow intensity in the major section of the LS, draining the left and lower parts of the body.The baseline values vary from about 20 to 420 µL/h.The lower panel characterizes the flow intensity in the minor section of the LS, draining the upper right part of the body.The baseline values vary from about 20 to only 100 µL/h.
The radius of the thoracic duct is known to be around 300 µm, while the radius of the lymphatic vessels afferent to popliteal nodes is around 20-40 µm (Table 1).For the human LS it is known that the largest lymphatic vessels have a diameter of about 2 mm and the diameter reduces to approximately 10-60 µm for initial lymphatic capillaries [10].Due to the lack of detailed anatomical and physiological data in mice on the diameters of all lymphatic vessels, we explored three complementary assumptions on the distribution of the radii of the edges of the lymphatic graph, as specified in the following three scenarios: • Scenario 1.All radii in the graph are assumed to be the same, equal to 150 µm (half of the radius on the thoracic duct).• Scenario 2. Edge radii decrease linearly with distance from the outlet vertices (jugular veins) to the inlet vertices.On the thoracic duct, the radius is assumed to be 300 µm, on the most distant edges (from the hindlimbs) it is assumed to be 41 µm.Therefore, on other edges from the inlet vertices, the value of the radius is equal to 41 µm and increases linearly when approaching the vein.On the subgraph collecting lymph into the right jugular vein, the radii are set symmetrically, equal to the radii in the left subgraph.• Scenario 3. Edge radii are distributed so that the cross-sectional area of incoming and outgoing vessels for each vertex of the graph is preserved.On all inlet edges, the radii are assumed to be the same and are estimated so that the radius on the thoracic duct would be equal to 300 µm.
The histograms of the vessel radii distribution for the above scenarios are presented in the left column of Figure 4.They clearly indicate that the median value of the vessel lumen decreases as we move from the first to the third scenario (from 150 µm (Scenario 1), to 115 µm (Scenario 2), to 60 µm (Scenario 3)).Note, that the estimated lymph flows shown in Figure 3 refer to Scenario 2.

Lymph Transfer Rates between Lymph Nodes
The key characteristic of the LS function is the rate of transfer of fluid through the system.Using the developed graph model, we estimated the lymph flow rates and the transition times between the lymph nodes.The computational results for three scenarios are summarized in Figure 4.The central column provides the estimates of the flow velocity.The median values turn out to be the smallest for the scenario of uniform vessel diameters and the largest under the assumption of conservation of the cross-sectional area at the vessel junctions.In particular, it increases from about 66 µm/s (Scenario 1), to 242 µm/s (Scenario 2), to 409 µm/s (Scenario 3).In addition, it is predicted to be practically homogeneous across the LS in the third case.The respective median transfer times between neighboring edges increase from about 19 s to 95 s with the longest transfer times from 2 min 25 s to 15 min 15 s when we compare the third and the first scenarios of the vessel geometry.

Sensitivity to Variations in Vessels Diameter
To analyze the sensitivity of the system to the diameter of the vessels, for all three scenarios, situations were simulated when the diameter of the vessel represented by the corresponding branch of the graph was varied by +10% and −10% from the initial one.The effect of change of the vessels' diameters was expressed in terms of the histograms of the relative change of pressure in the vessels of the LS as shown in Figure 5.In all three scenarios, the reduction of the radii led expectedly to a pressure increase, while the increase of the vessel radii had an opposite effect.The degree of variation was smallest for Scenario 1 and was largest for Scenario 2. The 10% diameter variation led to less than 1% change in the pressure for most of the vessels.

Topological Properties of the LS Graph
Following our previous study of the human LS [18], we characterized the topological properties of the murine LS graph model using some fundamental metrics.The quantified topological properties of the LS graph characterize the following structural and physiological features of the lymphatic network: the maximum number of lymphatics vessels entering and leaving a LN (maximum degree), the length of the shortest closed path (girth), the maximum distance between maximally separated nodes (diameter), the minimum distance between maximally separated nodes (radius), the typical distances between nodes (average path length), the characteristic of network regularity (energy and the spectral radius), characteristic of the network sparsity (edge density), the measure of how many nodes cluster together (clustering coefficient), the number of nodes critical for connectivity of the LS (number of separators), and the characteristics of lymph flow diversity (topological diversity of the vertices).
Let G = (n, m) be the graph with n nodes and m edges, respectively.Consider the following characteristics:

•
The number of input nodes N inp , i.e., the number of nodes with degree 1 and outdegree 0; • Maximum degree of graph ∆ G , i.e., the maximum degree of its vertices; • Girth of the graph g, which is the length of the shortest cycle in the simple graph; • Diameter, i.e., the longest geodesic distance (in other terms, maximum eccentricity of any vertex) • Average path length (mean geodesic distance)

•
The energy and the spectral radius of the graph are defined as follows, where λ j stand for the eigenvalues of the adjacency matrix A of the graph; • Edge density of the graph, i.e., the number of edges divided by the number of all possible edges,

•
The clustering coefficient C (transitivity) measures the probability that two neighbors of a vertex are connected.It can be computed as function of adjacency matrix A: • Number of separators n sep , i.e., the vertices removal of which disconnects the graph; • Topological diversity of the vertices as a function of the Shannon entropy associated with flow rates through the incident edges, where k is the number of v i 's incident edges and p ij is the proportion of the flow between the adjacent v i and v j to the total flow through the edges involving v i .The flow diversity is defined similar to the definition of network diversity in [33].
To analyze the robustness of the mouse LS to damage, we sequentially removed individual nodes of the graph and checked how many source vertices remained connected to the sink vertex into the circulatory system.The subgraphs of the LS and whole graph were analyzed.Accordingly, the robustness of the graph was estimated as the arithmetic mean of the ratio of the number of source vertices that retained the connection to the sink to their total number in the graph/subgraphs.The summary of topological properties of the murine LS graph model are presented in Table 2.The characteristics of the whole LS graph and the two subgraphs representing the LS parts draining the draining the left ∪ low and right parts of the body are specified.

Conclusions
In this work, we have developed a graph model of the LS network in mice.To define directions of edges in the original anatomy-based simple graph, we considered the mass balance of global lymph flow in the LS.The local interstitial pressure was not considered for the final definition of the oriented graph.It is the first study in which a geometric model of the murine LS has been developed and characterized in terms of its structural organization and the lymph transfer function.The study complements our previous analysis of the human lymphatic system [18].The developed graph model of the LS in mice provides a computational tool for studying the spatial aspects of the immune system functioning.It goes in line with recent advances in experimental techniques to characterize the whole-body dynamics of systemic infections in experimental mice [34].
The predicted orientation of the edges in the graph deserves further biological verification.To this end, various intravital imaging techniques could be used, e.g., near-infrared fluorescence imaging of lymphatic drainage patterns in mice [35], indocyanine green lymphangiography, or Doppler optical coherence tomography [4], to name just few of them [6].
We have previously implemented a similar approach to the analysis and modeling of the network structure of the human LS.The mouse and human LS differ fundamentally in terms of their cardinality, i.e., the number of elements comprising the system.The LS of mice consists of 28-36 LNs, whereas their number in the human LS ranges from about 500 to 1000.As a direct consequence, the graph model of the human LS is characterized by a much larger variability and more prominent randomness in its structure.The estimates of the probabilities of a new edge creation P e and the edge to connect nodes of different layers P o in random graph approximation of the LS quantify a larger uncertainty in the structure of the human LS graph compared to the mouse LS graph, i.e., 0.035 vs. 0.851 and 0.21 vs. 0.66, respectively, [18,36].In addition, the diameter, radius, average path length, and energy features of the respective LS graphs differ substantially.However, the topological properties of the graph models, such as the girth, spectral radius, edge density, and clustering coefficients are close for both human and mouse LSs.
The graph scheme of the mouse LS formulated to study the search time of antigen presenting cells by T cells in the LS has been recently presented in [36].However, the essential details of the LS, such as the adjacency matrix, the orientation of the edges, and the lymph flow estimates are not provided.In our study, we present a systematic analysis of the available anatomical and physiological data to develop the network model of the mouse LS, quantify the topological properties of the LS graph, and calculate flow through the network after making a series of assumptions about vessel diameters and terminal pressures.
The estimated parameters of the LS function in terms of the lymph flow rate and transfer time between various parts of the mouse body can be used in compartmental modeling for evaluation of the pharmacokinetic characteristics of drugs and adoptive cell therapies in advance of experiments.A remarkable example of the use of our recently developed graph model of human LS [18] is given in the study of how the topology of the lymphatic network affects the time required for an immune search through the lymphatic network to be completed [36].
To use the Hagen-Poiseuille equation for lymph flow, we assumed that the lymph is incompressible and Newtonian, the flow is laminar, and the vessels have constant circular cross-section with their length longer than the vessel diameter.However, the actual physiology of lymph flow through the LS is not considered, such as the lymphangion structure of the vessels, the pumping due to active contraction of the lymphatic muscle cells, adjacent tissue movement, and passive behavior properties of the vessels [5,14].
Further development of the presented graph model of the LS can be envisioned to proceed along three lines:

•
Considering the biomechanics of lymphatic pumping through a chain of lymphangions and lymph nodes; • Coupling the LS model with the cardiovascular system; • Integration with multi-physics models of the immune system.
The explored scenarios of lymph vessel radii reflect three different modalities of LS network construction.The computational results predict how the structural parameters impact the functional properties of the LS, such as lymph flow velocity and the transfer time between the nodes.These contribute to a better understanding of the LS in health and disease [6].Overall, our study meets the demand for quantitative rigorous approaches in the growing field of immunoengineering to utilize or exploit the lymphatic system for immunotherapy first in experimental animals and then to cure human immune-dependent diseases [4,37,38].

Figure 1 .
Figure 1.Oriented graph of the murine lymphatic system based on the anatomic data with 88 vertices and 87 edges.The vertices of the graph belong to four groups as detailed in the legend box: (1) lymph nodes (large blue), (2) outlet vertices with out-degree deg + = 0 corresponding to the sink into jugular veins (red), (3) connectors, i.e., the confluences of lymphatic vessels (light blue), (4) inlet vertices with in-degree deg − = 0 corresponding to the collecting lymphatics of various body tissues (orange).The vertex IDs and edge IDs are enumerated arbitrarily for correspondence with the matrix representation of the graph in Figure 2. The presented graph is specified in the CSV files containing the vertex and edge lists in Supplementary Materials.

Figure 2 .
Figure 2. Matrix representation of the anatomy-based graph of murine lymphatic system presented in Figure 1.The i-th and j-th vertices and edges are denoted as v i , v j and e i , e j , respectively.(A,B) Adjacency matrix A (for oriented and simple graph).(C) Incidence matrix, M. (D) Weighted Laplacian matrix, L = MGM T , normalized by the maximum matrix element.The conductance matrix G for constant vessel diameters is used for the illustration.The initial collecting lymphatics are likely to differ in the inflows as they absorb lymph from interstitial space characterized by volume and pressure varying across the body.However, the respective anatomical and physiological data to quantify the local impedance of the related edges are largely missing.Hence, we used a simplifying assumption that the flow velocities at all inlet vertices are the same and the pressure at the two sink nodes are equal.The following values of model parameters were used in computations of lymph flows: • number of inlet vertices n in = 52; • vessel radii range r ij = 40-300 µm; • vessel length l ij = 7-60 mm; • pressure at the sink nodes p out = 725 Pa; • lymph viscosity µ = 1.81 mPa•s; • lymph inflow q in = 0.4 mL/day.
lymph to the left jugular vein:

Figure 3 .
Figure 3.The distribution of the flows in the lymphatic vessels for Scenario 2, in which vessels radii decrease linearly with distance from the outlet vertices.Flows are shown in the subgraphs of the lymphatics collecting lymph into the left (A) and into the right (B) jugular vein.Colors and the arrows indicate the heads and the tails of the edges correspondingly.Edges are sorted by their distance from the sink.

Figure 4 .
Figure 4.The distribution of the velocities and transient transfer times of lymph flow in the vessels for the three scenarios of vessel radii distribution.At the bottom, the distribution of the edge lengths is shown, which is the same for all scenarios.

Figure 5 .
Figure 5.The distributions of the changes of pressures in the lymphatic vessels of the LS (vertices of the graph) after variation of the length diameters (±10%) for Scenario 1 (A,B), Scenario 2 (C,D), and Scenario 3 (E,F).
, v) is the geodesic distance (shortest oriented path connecting vertices u and v), (v) is the eccentricity of vertex v; • Radius of the graph (minimum eccentricity of any vertex),

Table 1 .
Physiological and anatomical properties of the murine lymphatic system.

Table 2 .
Summary statistics for the anatomy-based graph of murine lymphatic system.