Next Article in Journal
Method to Determine the Constitutive Permeability Parameters of Non-Linear Consolidation Models by Means of the Oedometer Test
Next Article in Special Issue
Polynomial Representations of High-Dimensional Observations of Random Processes
Previous Article in Journal
Promoters versus Adversaries of Change: Agent-Based Modeling of Organizational Conflict in Co-Evolving Networks
Previous Article in Special Issue
Update-Based Machine Learning Classification of Hierarchical Symbols in a Slowly Varying Two-Way Relay Channel
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Graph Theory for Modeling and Analysis of the Human Lymphatic System

by
Rostislav Savinkov
1,2,3,*,†,
Dmitry Grebennikov
1,2,4,*,†,
Darya Puchkova
5,
Valery Chereshnev
6,
Igor Sazonov
7 and
Gennady Bocharov
1,2,4,*
1
Marchuk Institute of Numerical Mathematics, Russian Academy of Sciences (INM RAS), 119333 Moscow, Russia
2
Moscow Center for Fundamental and Applied Mathematics at INM RAS, 119333 Moscow, Russia
3
Nikolsky Mathematical Institute, Peoples’ Friendship University of Russia (RUDN University), 117198 Moscow, Russia
4
Institute for Personalized Medicine, Sechenov First Moscow State Medical University, 119991 Moscow, Russia
5
Moscow Institute of Physics and Technology, National Research University, 141700 Dolgoprudny, Moscow Region, Russia
6
Institute of Immunology and Physiology, Ural Branch of Russian Academy of Sciences, 620000 Yekaterinburg, Russia
7
College of Engineering, Swansea University, Swansea SA1 8EN, UK
*
Authors to whom correspondence should be addressed.
These authors contributed equally to this work.
Mathematics 2020, 8(12), 2236; https://doi.org/10.3390/math8122236
Submission received: 27 November 2020 / Revised: 11 December 2020 / Accepted: 14 December 2020 / Published: 17 December 2020
(This article belongs to the Special Issue Random Processes on Graphs)

Abstract

:
The human lymphatic system (HLS) is a complex network of lymphatic organs linked through the lymphatic vessels. We present a graph theory-based approach to model and analyze the human lymphatic network. Two different methods of building a graph are considered: the method using anatomical data directly and the method based on a system of rules derived from structural analysis of HLS. A simple anatomical data-based graph is converted to an oriented graph by quantifying the steady-state fluid balance in the lymphatic network with the use of the Poiseuille equation in vessels and the mass conservation at vessel junctions. A computational algorithm for the generation of the rule-based random graph is developed and implemented. Some fundamental characteristics of the two types of HLS graph models are analyzed using different metrics such as graph energy, clustering, robustness, etc.

1. Introduction

The human lymphatic system (HLS) is a complex network of lymphatic organs linked through the lymphatic vessels. Its structure and functioning are critical to immunity by transporting the immune cells and antigens [1] and maintaining the fluid balance in tissues [2]. Currently, there are only a few studies considering the spatial structure of the human lymphatic system in terms of the network models [3,4,5]. The mathematical properties of the HLS network remain poorly investigated. The graph theory provides a powerful tool to implement, visualize and analyze the characteristics of the network models. Recently, it has been successfully applied for analysis of the topological properties and robustness of conduit networks in lymphatic nodes [6].
Generally, two complementary approaches to modeling the network structure of the lymphatic system. One is based on employment of available anatomical data [7,8]. The second approach uses some general rules of lymphatic vessels network organization as discussed in [5,9]. In this study, we develop oriented graph-type models of the HLS following both approaches. To specify the links direction in the anatomically derived HLS graph model, the analysis of the homeostatic lymph flow through the network is performed using the Poiseuille law. The topological characteristics of the graph models are analyzed using some general metrics such as graph energy, clustering, robustness, etc. The aim of the present study is to explore the fundamental mathematical properties of the lymphatic vessels network, for which the guiding organizational principles remain to be identified [2].
The paper is structured as follows. In Section 2, we present two anatomical data-based graph models, i.e., the Reddy’s [3] and the Plastic boy model [5] of the HLS. The lymph flow balance is analyzed to specify the respective directed graph of the HLS. In Section 3, computational algorithm for generating a random rule-based directed graph of the HLS is formulated and applied to derive the rule-based model. In Section 4, the parameters of the algorithm for the rule-based graph model are estimated to best-fit the topological characteristics of the anatomy-based model. In Section 5, a systematics analysis of the anatomy-based and the rule-based graph models of the HLS is performed. Conclusions and future work are discussed in the last section.

2. Graph Models Based on Anatomical Data

The physiological data on the structural organization of the HLS provide a foundation for developing anatomically based models of the HLS network. These days, the input information can be obtained from the Plasticboy project (http://www.plasticboy.co.uk/store/Human_Lymphatic_System_no_textures.html) or the CGTrader (2011–2020) at https://www.cgtrader.com/3d-models/character/anatomy/lymphatic-system-inhuman-body.

2.1. Reddy’s Directed Graph Model of the HLS

The first network model of the HLS has been systematically developed by Reddy [3]. Due to the scarcity of the available data at that time, the network model considers only the major lymphatic vessels, many organs are lumped into single vessels and the right side of the head, neck, thorax and upper right extremity are nor considered. The simplified Reddy’s network model of the HLS can be represented by oriented graph G ( V , E ) with 29 vertices and 28 edges as shown in Figure 1 generated using the network analysis package i g r a p h (https://igraph.org/).

2.2. Plastic Boy-Derived Graph Model of the HLS

We have previously developed a 3D computational geometry model of the HLS based on available anatomical data from the P l a s t i c b o y project [5]. This network is a graph G ( V , E ) containing 996 vertices and 1117 edges as shown in Figure 2. Unlike the Reddy’s graph model of the HLS shown in the Figure 1, the above model is an undirected graph. To transform it to an oriented graph, we need to determine the directions of the lymph flow through the graph. To this end, the simulation of a steady lymph flow based on Poiseuille equation is employed.

2.2.1. Undirected Graph of the HLS

The anatomical data-based 3D graph model contains several vertices of degree 2 corresponding to the lymphatic vessels which are composed of several lymphangions. We remove from the graph all vertices of degree 2 corresponding to individual lymphangions to represent the multi-lymphangion vessels as single graph edges. Moreover, we add one edge to correct the connection of right arm to the body (which was missing in the original graph) and three edges to get a single output vertex for the whole HLS graph (modeling the excretory vessels emptying the lymph to the venous system, i.e., the jugular lymphatic trunk, thoracic duct and jugular vein). Overall, the following specifications are made to obtain the HLS presented in Figure 2:
  • coordinates of all vertices in the graph are scaled to correspond to the basal 1750 mm height of a man;
  • the right hand and the head right side vertices are attached to the major part of HLS graph by adding the edge from vertex 987 to the vertex 993;
  • single output node (sink): (1) two vertices are added (numbered as 995 and 996) with vertex 996 being the output of the system, (2) three edges, i.e., between nodes 993–995, 994–995 and 995–996 were added;
  • we iterate over graph vertices and search for irrelevant vertices that meet all the following conditions:
    (a)
    their degree is equal to two;
    (b)
    they are not lymph nodes;
    (c)
    their neighbor vertices v a , v b are not connected with each other by an edge;
    at each iteration, the irrelevant vertex was removed, and the edge v a v b was created.
The topological properties of the HLS data-based graph models, i.e., the vertices degree- and edge-length distributions for the Plastic boy data-based graph and Reddy’s graph of the HLS are summarized in Figure 3A,B, respectively. For the data-based HLS graph, we additionally show the edge-length distribution and degree distribution of the vertices that represent lymph nodes (LNs) in Figure 3A.

2.2.2. Lymph Flow Analysis for Specifying the Directed Graph

The lymphatic system collects excess of intercellular fluid in the body tissues, drains it through the lymph nodes (LN) and redirects the lymph further along the lymphatic vessels, delivering the collected fluid into the central venous duct. In the LNs, about 90% of the afferent lymph flow takes a peripheral path through the subcapsular and medullary sinuses [10]. The remaining 10% enters interstitium to finally reach the efferent lymphatic vessel or to be absorbed by parenchymal blood vessels of the LN. Earlier experimental data show that about 10% of the lymph is absorbed into the blood vessels that penetrate the lymph nodes [11,12], while the remaining 90% of lymph entering LN is transported further through the lymphatic vessels to higher-level LN.
We study the balance of lymph flow in the undirected anatomical-based graph presented in Figure 2 to obtain the direction of lymph flow through the graph. We use the similar approach as described in [13]. We apply the Poiseuille equation to every internal vessel represented by edge ( i , j ) ,
Q i j = π R i j 4 ( p i p j ) 8 η L i j
where Q i j is the lymph flow rate along the vessel, R i j = R is the vessel lumen radius, η is the dynamic viscosity of the lymph, p i , p j are the pressures at the ends of the vessel, L i j is the length of the vessel. In the internal nodes, the mass conservation law is imposed:
j Q i j = 0 , i = 1 , , N N i n p N o u t ,
where N i n p = 357 is the number of input nodes with zero in-degree that provide the inflow of the lymph Q i n = 0.081 mm 3 /s, N o u t = 1 is the only output node with outflow Q o u t = N i n p Q i n = 28.93 mm 3 / s. To close the system (1) and (2) we specify the pressure p o u t = 1127.73 Pa on the output node and Q i n flow rate on the input nodes.
The above set of equations on HLS network in the matrix form
K p = b ,
to calculate the lymph flows where p = [ p 1 , , p N ] T is an unknown pressure vector for all the vertices of the graph G ( V , E ) . Matrix K and vector b are defined as
K i j = A i j π R 4 8 η L i j j i i V i n p V i n t k = 1 N A i k π R 4 8 η L i k j = i i V i n p V i n t δ i j i V o u t b i = Q i n i V i n p 0 i V i n t P o u t i V o u t
where A is an adjacency matrix of the LS graph; L i j is a distance between vertices i and j; V i n p , V i n t , V o u t are sets of input, internal and output vertices; δ i j is Kronecker’s delta. For the inflow vertices, we assume that all input vertices obtain equal amount of lymph ( Q i n ). For internal vertices, the conservation of mass requires the sum of incoming and outgoing flows to be zero. For the outflow vertices, we fix the pressure values P o u t .
To calculate the stationary flows, the following assumptions were made:
  • there are no other sinks of the lymph except for the exit vertex number 996: V o u t = { 996 } ;
  • all vertices of the graph that have only one connection with other vertices (except for the vertex 996) are the points of lymph entry (collecting lymphatics) into the lymphatic network;
  • the pressure at all inflow vertices is adjusted to provide the lymph flow rate on the input edges to be equal to the output flow rate of the system divided by number of vertices with zero input edges;
  • the radii of all the vessels are set to be same, equal to 1 mm;
  • the constant viscosity is set in the lymphatic network.
The value of the dynamic viscosity of lymph at a temperature of 37 ° C ( η = 1.8 × 10 3 Pa·s) is used. The pressure in the central venous duct ranges from 8 to 15 mmHg [2,14,15,16], and the value of 11.5 mmHg is considered in the calculations. According to existing data, from 2 to 3 L of lymph per day gets into the central venous duct [17,18], hence we take this parameter to be 2.5 L per day. The resulting stationary lymph flow distribution in the HLS is shown in Figure 4.
The above analysis enabled us to specify the lymph flow directions in the graph edges, and ultimately, to generate an oriented graph of the human lymphatic system, which is shown in Figure 5A. As the dimension of the graph HLS model is high, it can be visualized in 2D being subdivided into six segments representing specific parts of the human organism: arms, legs, head and body. The resulting oriented subgraphs are presented in the Figure 5B.
The adjacency matrix of this graph is presented in Figure 6.

3. Computational Algorithm for Generating a Random Rule-Based Directed Graph of the HLS

A quantitative modeling of the HLS presents a challenge due to enormous morphological complexity and variability in its appearance [5]. Scarcity of available anatomical data calls for the development of the so-called rule-based models of the HLS network structure. General rules of the lymphatic vessels network organization can be specified as follows: ( i ) no long-distance edges are allowed, ( i i ) nodes can have multiple inflow connections, ( i i i ) the nodes connections are locally acyclic, ( i v ) the lymphatic system network is a circular graph [9], and additionally, the links between nodes of the same layer are allowed. To build a random graph of the lymphatic system, we developed a rule-based algorithm that generates a random oriented graph with a predetermined number of layers N l , number of vertices N v , number of sources N i n p (input nodes, which should only have out-degree 1) and with one sink (output node). The layers in the algorithm are numerated from the output to the input vertices, from ground layer with index 1 to top layer with index N l . Iterative process on vertices is executed layer by layer starting from the next to the top layer. The rule-based graph HLS model generation consists of three major steps as specified in Algorithm 1 (C++ language realization of the algorithm is provided in Supplementary Materials, for any questions please contact authors of the article).
Algorithm 1: Generation of a rule-based directed graph of the HLS.
Mathematics 08 02236 i001
The optimal estimation of the algorithm parameters is performed in Section 4. To this end a topological fitness function specifying a mismatch between the anatomy-based graph and the rule-based graph is used. One realization of a random graph constructed using Algorithm 1 with the optimal parameters is presented in Figure 7, its degree distribution is shown in Figure 8, its adjacency matrix—in Figure 9.

4. Parameters of the Rule-Based Algorithm Providing Best Match to the Anatomy Data Graph

The basic topological properties of the graph are the node (vertex) degree- and the edge (arc) length distributions. These were presented for two types of graph in the previous sections. As the coordinates of nodes (and edge lengths) of the rule-based graph are not specified, we can only obtain statistics on the degree distribution of the vertices, and this indicator varies significantly depending on the parameter settings of the algorithm. We consider some fundamental characteristics of the HLS graph properties to make a deeper comparison of the graph models such as the energy, clustering, robustness, etc. Let G = ( V , E ) be an undirected version of the graph with n = V ( G ) and m = E ( G ) nodes and edges, respectively. Let A denote the n × n adjacency matrix of G.
  • The number of input nodes N i n p , i.e., the number of nodes with degree 1 and out-degree 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 (undirected) cycle in the graph.
  • Diameter, i.e., the longest geodesic distance (in other terms, maximum eccentricity of any vertex),
    D = max v V ϵ ( v ) = max v V max u V d ( u , v ) ,
    where d ( u , v ) is the geodesic distance (shortest directed path connecting vertices u and v), ϵ ( v ) is the eccentricity of vertex v.
  • Radius of the graph (minimum eccentricity of any vertex),
    r = min v V ϵ ( v ) = min v V max u V d ( u , v ) .
  • Average path length (mean geodesic distance),
    l G = 1 n ( n 1 ) u , v V , u v d ( u , v ) .
  • The energy and the spectral radius of the graph are defined as follows,
    E n ( A ) = j = 1 n λ j , ρ ( A ) = max { | λ j | } ,
    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,
    ρ d = m n ( n 1 ) .
  • 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 ( A ) = i = 1 , j = 1 , k = 1 n , n , n a i j · a j k · a k i i = 1 n ( ( j = 1 n a i j ) · ( ( j = 1 n a i j ) 1 ) ) .
  • Number of separators n s e p , i.e., the vertices removal of which disconnects the graph.
  • The robustness R of the graph can be defined as the fraction of peripheral vertices that retained the connection with the output vertex after removing 5% of the vertices selected randomly, averaged over k = 1 , , N r removal realizations. Given adjacency matrix A ˜ k of the k-th realization of the graph, indices i 1 , , i n i n p of its input vertices and index o of the output vertex, the robustness is computed as
    R = 1 N r k = 1 N r R k , R k = i = 1 n i n p F ( A ˜ k , i i , o , n ) n i n p , F ( A ˜ k , i i , o , n ) = 1 , i f j = 1 n 1 | ( A ˜ k j ) i i , o | > 0 , 0 , i f j = 1 n 1 | ( A ˜ k j ) i i , o | = 0 ,
  • Topological diversity of the vertices as a function of the Shannon entropy associated with flow rates through the incident edges,
    D f l o w ( v i ) = H ( v i ) log ( k ) = j = 1 k p i j log ( p i j ) log ( k ) , p i j = | Q i j | j = 1 k | Q i j | ,
    where k is the number of v i ’s incident edges and p i j 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 analogous to the social and spatial diversity of networks from [20].
The rule-based graphs generated with Algorithm 1 have five tuning parameters: (1) number of vertices N v , (2) number of input vertices N i n p , (3) number of layers N l , (4) probability of new edge creation P e at each step of the algorithm, probability that the created edge connects vertices from different layers P o . The first two parameters are set to match explicitly the properties of the anatomy-based graph: N v = 996 , N i n p = 357 . The number of layers can be estimated as N l = D + 1 = 41 , where D is the diameter of the directed target graph, which is increased by one because the output vertex is represented as a separate (ground) layer in Algorithm 1. Given the metrics characterizing the topology of an anatomy-based graph, we search for values of the remaining parameters P e and P o that can produce the graphs with similar topological structure. Specifically, we introduce the following state-vector
s ( G ) = m , g , D , r , l G , E n , ρ , ρ d , C , n s e p , n d e g 1 , n d e g 2 , n d e g 3 , n d e g 4 , n d e g 5 , n d e g 7 , n d e g 8 , Δ G T ,
which describes the topological properties of graph G ( V , E ) , and the objective function
Φ ( G ) = i = 1 18 s i ( G ) s i ( G ) s i ( G ) 2 ,
which penalizes the topological discrepancies of graph G from the target graph G and weighs them with s i ( G ) 2 to bring discrepancies of different components of vector s to a single scale. Here, we denote the number of vertices with degree i as n d e g i = v V , d e g ( v ) = i 1 , and exclude from the state-vector n d e g 6 , which equals zero for an anatomy-based graph, to avoid singularity in (14).
First, we investigate the landscape of topological fitness of algorithmically generated graphs over uniform ranges of admissible values of parameters P e and P o . The dependence of the objective function (14) on P e and P o is presented in Figure 10A. The colors correspond to logarithmically transformed mean values of objective function Φ m e a n over 640 algorithm realizations for each set of parameter values. One can see that there is a narrow region in which the objective function reaches its minimum on the grid cells where P e = 0.035 . Next, we examine more closely the dependence of the objective function distribution on parameter P o for fixed P e = 0.035 , which is presented in Figure 10B. One can see that for a long range of parameter P o , the objective function values remain low, until they increase for large values P o > 0.8 , where the maximum degree of graphs become too large (the dependence of topological metrics on P o is presented in Figure S1). From Figure 10, one can estimate the following ranges of viable parameter values for P e and N l : P e ( 0.001 , 0.07 ) , P o ( 0.01 , 0.8 ) . Two examples of rule-based graphs with satisfactory topological fitness generated by the algorithm with the use of optimal parameters P e = 0.035 , P o = 0.21 are shown in Figure 7.

5. Comparative Analysis of the HLS Graph Models

The differences between the anatomy data-based and the rule-based graphs in terms of the metrics introduced in Section 4 are summarized in Table 1 and Table 2. The histograms of robustness distributions ( R k values) for each type of graph calculated over N r = 100,000 ( N r = 435 for Reddy’s graph is the number of random deletions of 5 % of nodes corresponding to 29 one node deletions and 406 deletions of randomly chosen pairs of nodes) realizations of random deletion of 5% of vertices are presented in Figure 11. The distributions of flow diversity of anatomy-based graphs are presented in Figure 12. Please note that for data-based graph we use the flows computed in Section 2.2.2. For Reddy’s graph we obtain the flows assuming the equidistant inflow along 16 input edges equal to the 1/16-th of the outflow to the blue output node and the conservation law.
Given the algorithmically generated rule-based graph with optimal parameters presented in Figure 7A, we used another method to estimate the level of its resemblance to the anatomy-based graph. Specifically, to evaluate the level of isomorphism ( L I ) between the data-based and the rule-based graphs we searched for the row and column permutations to find the best match of them. Please note that examination of the graph isomorphism level belongs to the NP complexity calls. For this, we formed an expanded connectivity matrix M 1 = A O O O , M 2 = O O O B , D m a x = i , j ( M 1 i , j + M 2 i , j ) , where A is the adjacency matrix of the data-based graph, and B is the adjacency matrix of the rule-based graph, O is the zero matrix with equal sizes, D m a x is the initial Hamming distance for the given adjacency matrices. We used the simulated annealing method to minimize the Hamming distance between matrices M 1 , M 2 by performing permutations in the M 2 matrix (based on graph isomorphism):
i , j | M 1 i , j M 2 i , j | m i n
The permutations in the M 2 matrix are performed as follows: a pair of numbers i , j is chosen at random. Then, in the matrix M 2 , columns number i , j are swapped. After that, the same operation is performed with rows i , j . From the point of view of the graph, such transformations mean a permutation of the designations of the vertices i , j . Using this, we carry out a stochastic optimization of the M 2 matrix by the simulated annealing method to minimize (15). After that, we calculate the relative distance between graphs, presented by the adjacency matrices M 1 , M 2 as D r e l = i , j | M 1 i , j M 2 i , j | D m a x . Here, D r e l = 1 means a complete mismatch, D r e l = 0 means the complete equivalence of the adjacency matrices. By solving the minimization problem for the entire graphs, we achieved the upper bound estimate of the mismatch level of the range L I 0.49 .

6. Conclusions

The human lymphatic system is a complex system consisting of many components and performing several important functions related to the metabolism, fluid tissue homeostasis and the immune system [21,22]. Currently, there are only a few models addressing the network structure of the human lymphatic system in a systematic way [3,4,5]. In this work, we propose an approach based on the graph theory to modeling the human lymphatic system and analyzing the fundamental mathematical properties of the resulting graph models. Two different types of graph models are developed, the anatomical data- and the rule-based ones. The transformation of the anatomical data-based simple graph into the oriented graph is implemented on the base of the steady-state lymph flow balance in the lymphatic vessel network. This homeostatic balance is quantified by applying the Poiseuille equation to every vessel and the mass conservation law to every vessel junction. The available empirical data on lymph flow through the system are rather scarce and range from 0.004 mm 3 / s in rats to 46.3 mm 3 / s in humans [2]. Our computations consistently suggest the range of lymph flows from 0.001 to 10 mm 3 / s, depending on the location of the respective vessel.
Our study is aimed to understand the fundamental topological properties of the human lymphatic system and to critically analyze the underlying organizational principles (rules) of the LS that have been proposed for modeling the systemic spreading of HIV infection. The combination of the developed HLS data-based graph and the publicly available software SimVascular (https://simvascular.github.io/index.html) designed for blood flow simulations provide a solid basis to further move on to the 1D, 2D- and 3D-simulations of the lymphatic system flows. The knowledge gained on the lymphatic network topology can be useful for systems pharmacokinetics studies and drug therapies design.
Recent progress in tissue bioengineering and 3D bio-printing enables the development of artificial organs including those of the lymphatics system, e.g., the LN-on-a-chip [23]. The artificial LNs can be used as a functional cure for certain pathological conditions like lymphoedema characterized by disruption of the interstitial fluid transport. The graph model of the human lymphatic system provides a platform for an optimal positioning of the artificial LNs and lymphatic vessels to restore the functioning and performance of the HLS.
The lymphatic system (LS) is responsible for maintaining the fluid balance in tissues and the functioning of the immune system. As the complexity of the organization and regulation of the immune system is extremely high, the move towards a mechanistic understanding of the immune system functioning takes the form of specific rules (see [24,25,26,27]). The respective set of organizational principles (rules) of the LS structure was recently proposed and used for modeling the systemic spreading of HIV infection [9]. In our study we analyzed the consistency of the rule-based graph of the human lymphatic system with that based on the anatomical data. To this end, a novel computational algorithm for generation of the rule-based random graph has been developed. Some fundamental characteristics of the two types of HLS graph models are analyzed using metrics such as graph energy, clustering, robustness, etc.
Additionally, we estimated the optimal parameters of the proposed algorithm for generating the rule-based graph model which provide the best-fit to the anatomy data-based graph model with respect to basic topological metrics. Our analysis revealed that the currently available set of rules specifying universal properties of HLS for the rule–based graph modeling requires further maturation to be applied in clinical studies.
Overall, it is the first study in which the full complexity of the HLS network is addressed by methods of the graph theory methods. We probe to elucidate the efficacy of the rule-based approach to generate the HLS models. Further specification of the basic rules underlying the structure and function of the HLS is needed in collaboration with immunologists and physiologists. The graph models which we make publicly available, provide an important step towards a quantitative analysis of transport phenomena in the whole lymphatic system. This is a problem of a paramount importance for systems immunology, pharmacokinetics and medicine [1,2,28].

Supplementary Materials

The following are available online at https://www.mdpi.com/2227-7390/8/12/2236/s1, Figure S1: extended rule-based graphs metrics analysis; graph_edges.txt, graph_vertices.txt: data-based graph of the HLS; graph_generator.cpp, Makefile: C++ code to build rule-based graphs.

Author Contributions

Conceptualization, G.B., R.S. and D.G.; methodology, R.S., D.G., I.S., V.C. and G.B.; software, R.S., D.G., D.P.; validation, R.S., D.G. and I.S.; formal analysis, R.S., D.G.; investigation, R.S., D.G., D.P.; data curation, R.S., D.G.; writing—original draft preparation, R.S., D.G., G.B.; writing—review and editing, R.S., D.G., G.B., V.C., I.S.; supervision, G.B.; project administration, G.B.; funding acquisition, G.B. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Russian Science Foundation grant number 18-11-00171. R.S., D.G. and G.B. were partly supported by Moscow Center for Fundamental and Applied Mathematics (agreement with the Ministry of Education and Science of the Russian Federation No. 075-15-2019-1624) and by RFBR (project number 20-01-00352). R.S. Savinkov was partly supported by the «RUDN University Program 5-100».

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations

The following abbreviations are used in this manuscript:
HLSHuman Lymphatic System
LNslymph nodes

References

  1. Randolph, G.J.; Ivanov, S.; Zinselmeyer, B.H.; Scallan, J.P. The Lymphatic System: Integral Roles in Immunity. Annu. Rev. Immunol. 2017, 35, 31–52. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  2. Moore, J.E., Jr.; Bertram, C.D. Lymphatic System Flows. Annu. Rev. Fluid Mech. 2018, 50, 459–482. [Google Scholar] [CrossRef] [PubMed]
  3. Reddy, N.P.; Krouskop, T.A.; Newell, P.H., Jr. A Computer Model of the Lymphatic System. Comput. Biol. Med. 1977, 7, 181–197. [Google Scholar] [CrossRef]
  4. Mozokhina, A.S.; Mukhin, S.I. Pressure Gradient Influence on Global Lymph Flow. In Trends in Biomathematics: Modeling, Optimization and Computational Problems; Springer International Publishing: Berlin/Heidelberg, Germany, 2018; pp. 325–334. [Google Scholar] [CrossRef]
  5. Tretyakova, R.; Savinkov, R.; Lobov, G.; Bocharov, G. Developing Computational Geometry and Network Graph Models of Human Lymphatic System. Computation 2017, 6, 1. [Google Scholar] [CrossRef] [Green Version]
  6. Novkovic, M.; Onder, L.; Bocharov, G.; Ludewig, B. Topological Structure and Robustness of the Lymph Node Conduit System. Cell Rep. 2020, 30, 893–904.e6. [Google Scholar] [CrossRef] [Green Version]
  7. Qatarneh, S.M.; Kiricuta, I.C.; Brahme, A.; Tiede, U.; Lind, B.K. Three-dimensional atlas of lymph node topography based on the visible human data set. Anat. Rec. B New Anat. 2006, 289, 98–111. [Google Scholar] [CrossRef]
  8. Plasticboy. Plasticboy Pictures 2009 CC. Available online: http://www.plasticboy.co.uk/store/Human_Lymphatic_System_no_textures.html (accessed on 21 December 2017).
  9. Nakaoka, S.; Iwami, S.; Sato, K. Dynamics of HIV infection in lymphoid tissue network. J. Math. Biol. 2016, 72, 909–938. [Google Scholar] [CrossRef]
  10. Jafarnejad, M.; Woodruff, M.C.; Zawieja, D.C.; Carroll, M.C.; Moore, J.E., Jr. Modeling Lymph Flow and Fluid Exchange with Blood Vessels in Lymph Nodes. Lymphat. Res. Biol. 2015, 13, 234–247. [Google Scholar] [CrossRef] [Green Version]
  11. Adair, T.H.; Guyton, A.C. Modification of Lymph by Lymph Nodes. II. Effect of Increased Lymph Node Venous Blood Pressure. Am. J. Physiol. Heart Circ. Physiol. 1983, 245, H616–H622. [Google Scholar] [CrossRef]
  12. Adair, T.H.; Guyton, A.C. Modification of Lymph by Lymph Nodes. III. Effect of Increased Lymph Hydrostatic Pressure. Am. J. Physiol. Heart Circ. Physiol. 1985, 249, H777–H782. [Google Scholar] [CrossRef]
  13. Grebennikov, D.; Van Loon, R.; Novkovic, M.; Onder, L.; Savinkov, R.; Sazonov, I.; Tretyakova, R.; Watson, D.J.; Bocharov, G. Critical Issues in Modelling Lymph Node Physiology. Computation 2017, 5, 3. [Google Scholar] [CrossRef] [Green Version]
  14. Russell, P.S.; Hong, J.; Windsor, J.A.; Itkin, M.; Phillips, A.R.J. Renal Lymphatics: Anatomy, Physiology, and Clinical Implications. Front. Physiol. 2019, 10, 251. [Google Scholar] [CrossRef] [PubMed]
  15. Hariri, G.; Joffre, J.; Leblanc, G.; Bonsey, M.; Lavillegrand, J.-R.; Urbina, T.; Guidet, B.; Maury, E.; Bakker, J.; Ait-Oufella, H. Narrative Review: Clinical Assessment of Peripheral Tissue Perfusion in Septic Shock. Ann. Intensive Care 2019, 9, 37. [Google Scholar] [CrossRef] [PubMed]
  16. Martin, G.S.; Bassett, P. Crystalloids vs. Colloids for Fluid Resuscitation in the Intensive Care Unit: A Systematic Review and Meta-Analysis. J. Crit. Care 2019, 50, 144–154. [Google Scholar] [CrossRef]
  17. Sherwood, L. Human Physiology: From Cells to Systems; Cengage Learning: Belmont, CA, USA, 2012. [Google Scholar]
  18. Hall, J.E. Guyton and Hall Textbook of Medical Physiology, 13th ed.; Elsevier: Philadelphia, PA, USA, 2016; ISBN 9781455770052. [Google Scholar]
  19. Kamada, T.; Kawai, S. An Algorithm for Drawing General Undirected Graphs. Inf. Process. Lett. 1989, 31, 7–15. [Google Scholar] [CrossRef]
  20. Eagle, N.; Macy, M.; Claxton, R. Network Diversity and Economic Development. Science 2010, 328, 1029–1031. [Google Scholar] [CrossRef]
  21. Swartz, M. The Physiology of the Lymphatic System. Adv. Drug Deliv. Rev. 2001, 50, 3–20. [Google Scholar] [CrossRef]
  22. Hsu, M.C.; Itkin, M. Lymphatic Anatomy. Tech. Vasc. Interv. Radiol. 2016, 19, 247–254. [Google Scholar] [CrossRef]
  23. Shanti, A.; Samara, B.; Abdullah, A.; Hallfors, N.; Accoto, D.; Sapudom, J.; Alatoom, A.; Teo, J.; Danti, S.; Stefanini, C. Multi-Compartment 3D-Cultured Organ-on-a-Chip: Towards a Biomimetic Lymph Node for Drug Development. Pharmaceutics 2020, 12, 464. [Google Scholar] [CrossRef]
  24. Zinkernagel, R.M. Immunology and immunity against infection: General rules. J. Comput. Appl. Math. 2005, 184, 4–9. [Google Scholar] [CrossRef]
  25. Farber, D.L.; Netea, M.G.; Radbruch, A.; Rajewsky, K.; Zinkernagel, R.M. Immunological memory: Lessons from the past and a look to the future. Nat. Rev. Immunol. 2016, 16, 124–128. [Google Scholar] [CrossRef] [PubMed]
  26. Grossman, Z.; Paul, W.E. Dynamic Tuning of Lymphocytes: Physiological Basis, Mechanisms, and Function. Annu. Rev. Immunol. 2015, 33, 677–713. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  27. Grossman, Z.; Min, B.; Meier-Schellersheim, M.; Paul, W.E. Concomitant regulation of T-cell activation and homeostasis. Nat. Rev. Immunol. 2004, 4, 387–395. [Google Scholar] [CrossRef] [PubMed]
  28. O’Melia, M.J.; Lund, A.W.; Thomas, S.N. The Biophysics of Lymphatic Transport: Engineering Tools and Immunological Consequences. iScience 2019, 22, 28–43. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Figure 1. Oriented graph of a simple network model of the human lymphatic system suggested by Reddy et al. [3] with 29 nodes and 28 edges. The black node with out-degree d e g + = 0 corresponds to jugular vein (sink).
Figure 1. Oriented graph of a simple network model of the human lymphatic system suggested by Reddy et al. [3] with 29 nodes and 28 edges. The black node with out-degree d e g + = 0 corresponds to jugular vein (sink).
Mathematics 08 02236 g001
Figure 2. Graph of a human lymphatic system based on anatomical P l a s t i c b o y data [5] with 996 vertices and 1117 edges (graph files are provided in Supplementary Materials). (A) Face- and (B) side views. Blue colored vertices represent 272 lymph nodes. The black vertex corresponds to the output node.
Figure 2. Graph of a human lymphatic system based on anatomical P l a s t i c b o y data [5] with 996 vertices and 1117 edges (graph files are provided in Supplementary Materials). (A) Face- and (B) side views. Blue colored vertices represent 272 lymph nodes. The black vertex corresponds to the output node.
Mathematics 08 02236 g002
Figure 3. Statistical distributions of the anatomy data-based HLS graphs. (A) Statistics of the data-based graph, (B) of the Reddy’s graph. Statistics include the distributions of total degrees of the vertices, distributions of in- and out-degrees, as well as the edge-length distribution and degree distribution of the vertices that represent the lymph nodes (LNs) for data-based graph.
Figure 3. Statistical distributions of the anatomy data-based HLS graphs. (A) Statistics of the data-based graph, (B) of the Reddy’s graph. Statistics include the distributions of total degrees of the vertices, distributions of in- and out-degrees, as well as the edge-length distribution and degree distribution of the vertices that represent the lymph nodes (LNs) for data-based graph.
Mathematics 08 02236 g003
Figure 4. Steady-statedistribution of lymph flows in the HLS graph model. The output sink node is presented as black circle.
Figure 4. Steady-statedistribution of lymph flows in the HLS graph model. The output sink node is presented as black circle.
Mathematics 08 02236 g004
Figure 5. Oriented graph model of the HLS based on the analysis of the lymph flow directions in the network. (A) Spatial view of the graph with numbered vertices (face-view projection of 3D coordinates). (B) View of the graph on a plane with no-overlapping layout. The 2D layout was obtained using spring-based model algorithm [19] implemented in package igraph, followed by some manual tweaking to space out the vertices. The colored subgraphs correspond to different parts of the body, i.e., arms (dark and light blue), head/neck (magenta), torso (light green), legs (red and gray). LNs are represented by larger circles.
Figure 5. Oriented graph model of the HLS based on the analysis of the lymph flow directions in the network. (A) Spatial view of the graph with numbered vertices (face-view projection of 3D coordinates). (B) View of the graph on a plane with no-overlapping layout. The 2D layout was obtained using spring-based model algorithm [19] implemented in package igraph, followed by some manual tweaking to space out the vertices. The colored subgraphs correspond to different parts of the body, i.e., arms (dark and light blue), head/neck (magenta), torso (light green), legs (red and gray). LNs are represented by larger circles.
Mathematics 08 02236 g005
Figure 6. Adjacency matrix of the oriented data-based graph of human lymphatic system shown in Figure 5. Colors represent the lengths L i j of the edges e i j .
Figure 6. Adjacency matrix of the oriented data-based graph of human lymphatic system shown in Figure 5. Colors represent the lengths L i j of the edges e i j .
Mathematics 08 02236 g006
Figure 7. Two examples (A,B) of algorithmically generated realizations of rule-based directed random graphs of the HLS. The output nodes are marked with black color. Algorithm parameters: N v = 996 , N i n p = 357 , N l = 41 , P e = 0.035 , P o = 0.21 . The metrics of topological fitness to anatomy-based graphs presented in figure are defined in Section 4.
Figure 7. Two examples (A,B) of algorithmically generated realizations of rule-based directed random graphs of the HLS. The output nodes are marked with black color. Algorithm parameters: N v = 996 , N i n p = 357 , N l = 41 , P e = 0.035 , P o = 0.21 . The metrics of topological fitness to anatomy-based graphs presented in figure are defined in Section 4.
Mathematics 08 02236 g007
Figure 8. Degree distribution of the rule-based graph presented in Figure 7A.
Figure 8. Degree distribution of the rule-based graph presented in Figure 7A.
Mathematics 08 02236 g008
Figure 9. Directed graph adjacency matrix of the rule-based random graph presented in Figure 7A.
Figure 9. Directed graph adjacency matrix of the rule-based random graph presented in Figure 7A.
Mathematics 08 02236 g009
Figure 10. The landscape of topological fitness of rule-based graphs over different values of algorithm parameters. (A) Image of log 10 ( Φ m e a n ) values as function of varying parameters P e and P o and fixed N l = 41 . The mean values are calculated in each image cell over 640 graphs realizations of rule-based graphs generated with corresponding parameter sets. (B) Dependence of objective function on P o with fixed P e = 0.035 , N l = 41 . The statistics on values of objective function (14) are calculated over 10,000 realizations of the algorithm for each value of P o .
Figure 10. The landscape of topological fitness of rule-based graphs over different values of algorithm parameters. (A) Image of log 10 ( Φ m e a n ) values as function of varying parameters P e and P o and fixed N l = 41 . The mean values are calculated in each image cell over 640 graphs realizations of rule-based graphs generated with corresponding parameter sets. (B) Dependence of objective function on P o with fixed P e = 0.035 , N l = 41 . The statistics on values of objective function (14) are calculated over 10,000 realizations of the algorithm for each value of P o .
Mathematics 08 02236 g010
Figure 11. Histograms of the robustness of (A) Reddy’s graph, (B) anatomical data-based graph and (C) rule-based random graph presented in Figure 7A.
Figure 11. Histograms of the robustness of (A) Reddy’s graph, (B) anatomical data-based graph and (C) rule-based random graph presented in Figure 7A.
Mathematics 08 02236 g011
Figure 12. Histograms of the flow diversity of the inner vertices of anatomy-based graphs: (A) for data-based graph, (B) for Reddy’s graph.
Figure 12. Histograms of the flow diversity of the inner vertices of anatomy-based graphs: (A) for data-based graph, (B) for Reddy’s graph.
Mathematics 08 02236 g012
Table 1. Summary statistics for Reddy’s, anatomy-based and rule-based graphs characterizing their topological properties. For algorithmically generated graphs, we present the statistics obtained over 10,000 graphs generated with algorithm parameters P e = 0.035 , P o = 0.21 , N l = 41 . Robustness was calculated over 1000 algorithm realizations, with 1000 removal attempts for each graph.
Table 1. Summary statistics for Reddy’s, anatomy-based and rule-based graphs characterizing their topological properties. For algorithmically generated graphs, we present the statistics obtained over 10,000 graphs generated with algorithm parameters P e = 0.035 , P o = 0.21 , N l = 41 . Robustness was calculated over 1000 algorithm realizations, with 1000 removal attempts for each graph.
Lymphatic Vascular
System Graph Model
Reddy’s
Model
Anatomy-Based
Model
Rule-Based
Model (Mean)
Rule-Based
Model (SD)
Rule-Based m.
(Min-Max Range)
G ( n , m ) ( 29 , 28 ) ( 996 , 1117 ) ( 996 , 1029 ) ( 0 , 6.4 ) (996, 1009–1056)
N i n p 163573570(357–357)
Maximum degree, Δ G 48161.58(8–21)
Girth, g0340.9(3–14)
Diameter, D540 39.96 0.22 (37–40)
Radius, r43028 2.3 (23–38)
Average path length, l G 2.46 12.79 15.3 0.86 (13.6–18)
Energy, E n 32.1 1224.5 1190 4.9 (1173–1203)
Spectral radius, ρ 2.58 3.51 4.18 0.19 (3.28–4.72)
Edge density, ρ d 0.034 0.001127 0.001038 6.4 × 10 6 (0.00102–0.00107)
Clustering coefficient, C0 0.027 0.0004 0.0008 (0–0.0036)
Number of separators, n s e p 1340149620(437–548)
Robustness, R 0.7 0.6 0.66 0.05 (0.45–0.77)
Average flow diversity, D f l o w ¯ 0.9080.996
Table 2. Summary statistics for the subgraphs of the anatomy-based graph which correspond to different parts of the body.
Table 2. Summary statistics for the subgraphs of the anatomy-based graph which correspond to different parts of the body.
SubgraphLeft ArmRight ArmHead & NeckTorsoLeft LegRight Leg
G ( n , m ) ( 149 , 181 ) ( 141 , 170 ) ( 198 , 208 ) ( 168 , 183 ) ( 163 , 176 ) ( 177 , 192 )
Number of input nodes464564686671
Number of output nodes112111
Maximum degree, Δ G 884854
Girth, g334333
Diameter, D232221182222
Radius, r1091391313
Average path length, l G 8.77 8.37 7.24 7.02 8.95 9.55
Energy, E n 192.22 180.6 236.24 164.4 196.57 214.1
Spectral radius, ρ 3.3 3.32 2.95 3.51 2.92 3.08
Edge density, ρ d 0.0082 0.0086 0.0053 0.0065 0.0067 0.0062
Clustering coefficient, C 0.049 0.035 0 0.055 0.01 0.009
Number of separators, n s e p 414146 (left), 46 (right)816672
Number of LNs363651 (left), 51 (right)592118
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Savinkov, R.; Grebennikov, D.; Puchkova, D.; Chereshnev, V.; Sazonov, I.; Bocharov, G. Graph Theory for Modeling and Analysis of the Human Lymphatic System. Mathematics 2020, 8, 2236. https://doi.org/10.3390/math8122236

AMA Style

Savinkov R, Grebennikov D, Puchkova D, Chereshnev V, Sazonov I, Bocharov G. Graph Theory for Modeling and Analysis of the Human Lymphatic System. Mathematics. 2020; 8(12):2236. https://doi.org/10.3390/math8122236

Chicago/Turabian Style

Savinkov, Rostislav, Dmitry Grebennikov, Darya Puchkova, Valery Chereshnev, Igor Sazonov, and Gennady Bocharov. 2020. "Graph Theory for Modeling and Analysis of the Human Lymphatic System" Mathematics 8, no. 12: 2236. https://doi.org/10.3390/math8122236

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop