Network Modelling of Murine Lymphatic System

1. Introduction

Networks of various nature, 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 a spatially distributed components of the lymphatic system (LS) [3].

Animal models of diseases in particular, the mice, are considered to be the cornerstone for translational research in immunology [4,5]. The research with laboratory mice enabled invaluable insight into mammalian immune systems [6]. Despite numerous advances in understanding the immune system from mouse studies, there exist fundamental differences between mouse and human immune systems [5]. The structural organization of the lymphatic system represents a straightforward example. However, a mathematical network-based characterization of the LS in mice is still unavailable.

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 [7]). In [7], we developed a computational algorithm for representing the anatomy-based and rule-based graphs of the LS in humans. The graph models enabled the analysis of different metrics of complexity of human LS such as spectral radius, clustering coefficient, average path length, number of separators. A similar analysis for the mouse LS remains to be done.

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 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.

2. 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 [8,9]. 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 in shown in Figure 1. It is developed using anatomical descriptions from [8,9]. The vertices of the graph refer to either the lymph nodes, outlet vertices with out-degree $de{g}^{+}=0$ corresponding to the sink into jugular veins, the confluences of lymphatic vessels, or inlet vertices with in-degree $de{g}^{-}=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₂O [15,16]. Lymph viscosity is taken to be equal to 1.81 mPa·s.

The anatomy data enable to specify a simple undirected graph of the murine LS. The adjacency matrix A of the LS graph is shown in Figure 2 (B).

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 undirected 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 studies on fluid dynamics in various parts of the LS in mice ([10,12,13,14]) and computation of the lymph flow through the system in accordance to an overall mass balance using the Hagen-Poiseuille equations.

3. 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$) (Figure 1,Figure 3). The thoracic duct belongs to the subgraph $g_l$.

3.1. 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₂O, and the outflow from the left LS subgraph to the left jugular vein $q_{out}^{\left(l\right)}$ 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}^{\left(l\right)}=q_{out}^{\left(l\right)}/n_{in}^{\left(l\right)}$, where $n_{in}^{\left(l\right)}$ 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_{in}^{\left(r\right)}=q_{in}^{\left(l\right)}$ as in the left one, and the same output pressure $p_{out}$. The outlet outflow is given by $q_{out}^{\left(r\right)}=q_{in}^{\left(r\right)}{n}_{in}^{\left(r\right)}$.

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 Hagen-Poiseuille equation:

  $$q_{ij}=g_{ij}(p_i-p_j)=\frac{\pi r_{ij}^4}{8\mu l_{ij}}(p_i-p_j),$$

  (1)

  which 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 $\mu$ and the radii and the lengths of the edges,
• and by the balance of flow through the vertices due to mass conservation:

  $$\sum _{j\in \mathcal{A}\left(i\right)}{q}_{ij}=\left\{\begin{array}{cc}{q}_{in},\hfill & \mathrm{if}i\mathrm{is}\mathrm{inlet}\mathrm{vertex}\hfill \\ -q_{out},\hfill & \mathrm{if}i\mathrm{is}\mathrm{outlet}\mathrm{vertex}\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

  (2)

  where $\mathcal{A}\left(i\right)$ is a set of vertices adjacent to i.

Using the oriented incidence matrix $M\in {\mathbb{R}}^{n\times m}$ and the diagonal conductance matrix $G\in {\mathbb{R}}^{m\times m}$ (with elements indexed by edges rather than nodes), one can rewrite the equations (1-2) as

$$q=-G{M}^{T}p,Mq=-\widehat{q},$$

(3)

where $p\in {\mathbb{R}}^n$ are the nodal pressures, $\widehat{q}\in {\mathbb{R}}^n$ are the net flows through the vertices and $q\in {\mathbb{R}}^m$ are the flows through the edges. Hence, we get the linear system to solve for the nodal pressures:

$$MG{M}^{T}p=Lp=\widehat{q},$$

(4)

where $L=MG{M}^T$ is a symmetric weighted Laplacian matrix.

As at the outlet vertex the pressure is known ($p_{out}$), we substitute the vector $\widehat{p}={[0,\cdots ,p_{out},\cdots ,0]}^T$ of zeros with the known pressure at the corresponding index in (4). By shifting $L\widehat{p}$ to the right hand side, we obtain the system

$$L_{rect}{p}_{unknown}=\widehat{q}-L\widehat{p},$$

(5)

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}^{+}\left(\widehat{q}-L\widehat{p}\right)$.

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.

3.2. 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 directed graph, the adjacency matrix is shown in Figure 2 (A). The adjacency matrix of the undirected graph presented in Figure 2 (B) is symmetric. A complementary representation of the graph is provided by the incidence matrix (see Figure 2 (C)). 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 2 (D). It is related to the degree matrix D and the adjacency matrix A of the graph $L=D-A$, representing an edge-weighted graph.

4. 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 $\mu L/hr$. 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 $\mu L/hr$.

The radius of the thoracic duct is known to be around 300 $\mu$m, while the radius of the lymphatic vessels afferent to popliteal nodes is around 20-40 $\mu$m (Table 1). Due to the lack of detailed data on the diameters of all lymphatic vessels, we made 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 $\mu$m (half of the radius on the thoracic duct).
• Scenario 2. Edge radii vary linearly with distance from the outlet vertices (jugular veins) to the inlet vertices. On the thoracic duct, the radius is assumed to be 300 $\mu$m, on the most distant edges (from the hind legs)– 41 $\mu$m. Therefore, on other edges from the inlet vertices, the value of the radius is equal to 41 $\mu$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 $\mu$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 1st to the third scenario (from 150 $\mu m$ (Scenario 1), to 115 $\mu m$ (Scenario 2), to 60 $\mu m$ (Scenario 3)).. Note, that the estimated lymph flows shown in Figure 3 refer to Scenario 1.

4.1. 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 one for the scenario of a uniform vessel diameters and the largest one under assumption of conservation of the cross-sectional area at the vessel junctions. In particular, it increases from about 66 $\mu m/s$ (Scenario 1), to 242 $\mu m/s$ (Scenario 2), to 409 $\mu m/s$ (Scenario 3). In addition, it is predicted to by practically homogeneous across the LS in the third case. The respective median transfer times between neighboring edges increase from about 19 seconds to 95 second with the longest transfer times from 2 minutes 25 seconds to 15 minutes 15 seconds when we compare the third and the first scenarios of the vessel geometry.

4.2. 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 les than 1% change in the pressure for most of the vessels.

5. Topological properties of the LS graph

Following our previous study of the human LS [7], we characterize the topological properties of the murine LS graph model using some fundamental metrics. 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 out-degree 0;
• maximum degree of graph ${\Delta }_G$, i.e., the maximum degree of its vertices;
• girth of the graph g, which is the length of the shortest (undirected) cycle in the graph;
• diameter, i.e., the longest geodesic distance (in other terms, maximum eccentricity of any vertex),

  $$D=\underset{v\in V}{max}ϵ\left(v\right)=\underset{v\in V}{max}\underset{u\in V}{max}d(u,v),$$

  (6)

  where $d(u,v)$ is the geodesic distance (shortest directed path connecting vertices u and v), $ϵ\left(v\right)$ is the eccentricity of vertex v;
• radius of the graph (minimum eccentricity of any vertex),

  $$r=\underset{v\in V}{min}ϵ\left(v\right)=\underset{v\in V}{min}\underset{u\in V}{max}d(u,v).$$

  (7)
• average path length (mean geodesic distance),

  $$l_G=\frac{1}{n(n-1)}\sum _{u,v\in V,u\ne v}d(u,v).$$

  (8)
• the energy and the spectral radius of the graph are defined as follows,

  $$En\left(A\right)=\sum _{j=1}^n\mid {\lambda }_j\mid ,\rho \left(A\right)=max\left\{\right|{\lambda }_j\left|\right\},$$

  (9)

  where ${\lambda }_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,

  $${\rho }_d=\frac{m}{n(n-1)}.$$

  (10)
• 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:

  $$C\left(A\right)=\frac{{\sum }_{i=1,j=1,k=1}^{n,n,n}{a}_{ij}·a_{jk}·a_{ki}}{{\sum }_{i=1}^n(\left({\sum }_{j=1}^n{a}_{ij}\right)·(\left({\sum }_{j=1}^n{a}_{ij}\right)-1))}.$$

  (11)
• 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,

  $$D_{flow}\left(v_i\right)=\frac{H\left(v_i\right)}{log\left(k\right)}=\frac{-{\sum }_{j=1}^k{p}_{ij}log\left(p_{ij}\right)}{log\left(k\right)},p_{ij}=\frac{|Q_{ij}|}{{\sum }_{j=1}^k|Q_{ij|}},$$

  (12)

  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 [?].

To analyze the robustness of the mouse LS to damage, we sequentially removed individual nodes of the graph and checked how many source vertices remain 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 topological properties of the murine LS graph model are presented in Table 1. The characteristics of the whole LS graph and the two subgraphs representing the LS parts draining the right∪low and left parts of the body are specified.

6. Conclusions

In this work, we have developed a graph model of the LS network in mice. 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 [7]. 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 [17].

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 on compartmental modelling 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 [7] is given in the study of how the topology of the lymphatic network affects the time for an immune search through the lymphatic network to be completed [18].

Further development of the presented graph model of the LS can be visioned 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.

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 [19,20,21].