Mathematical Modelling of the Structure and Function of the Lymphatic System

Abstract

This paper presents current knowledge about the structure and function of the lymphatic system. Mathematical models of lymph flow in the single lymphangion, the series of lymphangions, the lymph nodes, and the whole lymphatic system are considered. The main results and further perspectives are discussed.

1. Introduction

The lymphatic system (LS) in the human body complements the cardiovascular system and provides drainage of interstitial fluid through a complex system of lymphatic vessels and lymph nodes. Various molecules, pathogens, and immune cells are transported in the organism along the LS. It performs drainage, transport, immune, and deposit functions, namely, lymph accumulation in the lymphatic capillary network, lymph nodes, and spleen.

For a long time, the role of the LS was underestimated and it was not sufficiently well studied. Only in the 20th century, its importance for normal functioning in humans and its role in the pathogenesis of various diseases were realized [1]. One of the most common disorders related to the LS is lymphedema, which can develop as a result of surgical treatment and because of other reasons. It is also reported that solid tumors can create metastases due to malignant cells travelling along the LS [2]. Recent observations correlate plastic bronchitis with dysfunction of the LS [3]. Other examples (lymphangioleiomyomatosis, Kaposi’s sarcoma, lymphoceles, lymphatic malformations, and chylothorax) can be found in [1,4].

The presence of lymphatic structures—meningeal lymphatics or dural lymphatics—in the brain was recently discovered [5,6], opposing the previous belief that there are no lymphatic vessels in the brain. This discovery can be important for understanding the pathogenesis of neural diseases, including Alzheimer’s disease, and for the development of efficient treatments of neurological disorders [7].

Despite progress in investigations of the LS, the parameters of lymph flow in the human body remain basically unknown because of their difficulty to be investigated: lymphatic vessels are small, so it is hard to place sensors; lymphs are low scattering and its signals are near the noise threshold, so the Doppler-type technique has failed until recently; the LS starts from the interstitial space, so contrast cannot be injected globally; and the LS contains lymph nodes filtrating lymphs, so regional injections allow contrasting lymphs only to the first node. Moreover, injections of contrast alter the interstitial pressure and, therefore, the measured parameters. Modern methods of lymph investigation in vivo [8] are more suitable for lymphangiography than for determination of the flow parameters. Data about lymph flow velocity are usually calculated by dividing the length of the imaged vessel by the time of observation [9,10,11,12,13,14], which allows getting mean lymph velocity but not the the characterization of lymph pulsations. Moreover, these techniques are invasive and can alter the flow [15]. Recently, a Doppler optical coherence tomography platform with Doppler algorithm operating on low signal-to-noise measurements has been developed [16] to get lymph flow characteristics. It allows getting parameters of lymph flow in vivo including visualization of its pulsating nature. This approach has certain limitations which do not allow using it widely. However, it can open new perspectives in the determination of lymph flow parameters.

Some data about lymph flow are available from in vitro experiments. In the common setup, a part of the cattle mesenteric lymphatic vessel is fixed in the bath with a special solution and the pressure at both ends, the diameter, and the output flux are measured [17,18,19,20]. These data can differ from those under physiological conditions [15], but they are used in mathematical models as model parameters and for checking the model output. Mathematical models are used for investigation of lymph flow in the series of contracting parts of lymphatic vessels called lymphangions. The underlying assumptions in these models are essential, and depending on them, the models give different, sometimes opposite results (e.g., [21,22], see Section 3.2). In what follows, we consider lymph transport in different parts of the LS, the influence of valves and contractions on the flow, their modeling, and results. We will illustrate the main approaches to lymph flow modeling on the example of the models described in [23,24] including models developed in the last years and will consider some possible applications.

2. Transport Function of the LS

The lymphatic system (LS) includes lymphatic vessels and lymph nodes. It complements the cardiovascular system and performs drainage, transport, immune, and deposit functions. About 10% of blood goes to the LS in the process of capillary filtration [25,26]. The total lymph flux returning to the venous system is 2–4 L/day [25,26,27,28]. The LS absorbs excess fluid, waste products, and large molecules from the interstitial space and transports it to the veins of the cardiovascular system—to the left and right venous angles between neck veins—so the LS starts from the interstitial space and opens into the upper vena cava (Figure 1A). The initial part of the LS is called initial lymphatics. They absorb lymph (tissue fluids, cells, and large extracellular molecules) from the interstitial space. Initial lymphatics have one layer of endothelial cells and no or almost no basal membrane. They connect with anchoring filaments to the surrounding extracellular matrices, and this positioning prevents them from collapsing even under increased pressure in the interstitial fluid. Endothelial cells of initial lymphatics form junctions, where big cells can go into, and prevent lymphs from going out from the initial lymphatics into the interstitial space (Figure 1B). These overlapping cells are called primary (flap) valves. Here and below, anatomical descriptions are provided according to [1,27,29,30,31].

From the initial lymphatics, the lymph goes to the collectors, which have three layers of cells similar to blood vessels. There are secondary valves in these collectors which prevent the backward flow of lymphs. These valves consist of two or three leaflets located along with the flow (Figure 2A). Valves divide lymphatic vessels into functional parts, and the parts between adjacent valves are called lymphangions (Figure 2B). These lymphangions can produce active contractions. Contractions in the vessels with valves provide the unidirectional flow of lymph (Figure 2B), which is a characteristic property of lymph flow in the LS.

The network of lymphatic vessels is highly unstructured and contains many lymph nodes, and its hierarchical structure is shown in Figure 3. Numerous lymph nodes contained in the LS are lymphoid organs actively participating in the transport function of the LS [32].

It is believed now that lymph flows in the LS are because of the pressure gradient between the interstitial space and the upper vena cava and because of passive (extrinsic) contractions which are caused by the movements of muscle tissues, the diaphragm, and big blood vessels nearby and of active (intrinsic) contractions which are spontaneous contractions of lymphangions. The precise mechanism of active contractions is not yet well known [33]. One of the open questions is whether active contractions can produce enough force to drive the lymphs throughout the whole LS. The lymph flow cycle includes the entrance of lymphs into the initial lymphatics, the flow into lymphangions, and the flow into lymph nodes. The biomechanics of lymph flow in all these steps is not yet elucidated.

The majority of models of lymph flow in the lymphangion or series of lymphangions are based on the zero-dimensional (including lumped ones) [17,22,37,38] or one-dimensional (including quasi-one-dimensional ones) [18,39,40,41] approaches. “Lumped” models here are models based on electrical circuit theory, where the network of vessels is substituted with the circuit and where the voltage, current, and charge are equivalent to pressure, lymph flux, and mass, respectively [38,42,43]. Zero-dimensional models are models with no space distribution of velocity (flux) and pressure functions. Quasi-one-dimensional models are one-dimensional models in which the cross section may change. A three-dimensional formulation for lymph flow in the single lymphangion is considered in [44]. In the works [40,45], models of lymph flow in the whole LS are developed on the basis of zero-dimensional and quasi-one-dimensional approaches, respectively. The influence of valves and contractions on the flow distinguishes lymph flows from other physiological flows.

3. Flow in Lymphangions

We will describe lymph flow in this section using the quasi-one-dimensional approach [40] and other models when they are available. This approach is widely used to model systemic blood circulation [46,47] and, more recently, to model lymph flow in the lymphatic vessels [18,40,41]. In this approach, fluid flow is described by a system of hemodynamic equations which are derived from mass and momentum conservation (or from the Navier–Stocks equations [48]):

$$\frac{\partial s}{\partial t}+\frac{\partial us}{\partial x}=0,$$

(1)

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+\frac{1}{\rho }\frac{\partial p}{\partial x}=f,$$

(2)

where $u(x,t)$ is fluid velocity, $p(x,t)$ is pressure, $s(x,t)$ is cross sections area, $\rho =const$ is density, x is spatial coordinate, and t is time. The external force $f(x,t)$ corresponds to viscous friction. Systems (1) and (2) are completed by the relation between the cross-sectional area s and pressure p. While this system is suitable to describe blood flow in arteries and veins, it cannot be used to describe lymph flow in lymphatic vessels without modifications because of the presence of valves and vessel wall contractions.

3.1. Influence of Valves

Secondary valves in lymphatic vessels prevent backward flow of lymphs and provide its unidirectional flow. Some investigations report that valves are included in the contraction of lymphangion [31]. Leaflets of valves can influence the flow. The models of secondary valves are usually limited by the absence of backward flow and sometimes by the resistance to the incoming flow.

Secondary valves can prolapse in a high adverse pressure gradient, and the presence of backward flow is one of the diagnostic criteria for lower limb lymphedema [49]. However, it is not known whether it is the reason for or a consequence of the disease [17]. Some models of valves include the possibility of prolapse.

There are several ways to represent the influence of valves on the lymph flow. One of them is to consider adjacent lymphangions as two segments connected with the bifurcation point. In this representation, the flow in each vessel is described by System (1, 2) and a closing relation. At the point of bifurcation, conditions on the pressure and flux are stated. Valves are taken into account in the flux condition at the bifurcation point, and the simplest model for their behavior is as follows [17,37]:

$$q_v\ge 0,$$

(3)

where $q_v$ is the lymph flux through the valve. Condition (3) restricts the backward lymph flow through the valve. Valve resistance to the incoming flow is considered in [37] on the basis of the data from [50].

A valve model in blood vessels for lumped and one-dimensional approximations is proposed in [51]. It is sufficiently simple and detailed to describe normal functioning of valves and their pathologies. This model is used in [39] to describe the valve influence on lymph flow. The equation for time-varying flux through valve $q_v$ has the following form:

$$\frac{\mathrm{d}}{\mathrm{d}t}{q}_v=\frac{1}{L\left(\xi \right)}\left(\Delta p\left(t\right)-R\left(\xi \right)q_v-B\left(\xi \right)q_v\left|q_v\right|\right),$$

(4)

where $\Delta p$ is pressure drop on the valve and where $L\left(\xi \right)$, $B\left(\xi \right)$, $R\left(\xi \right)$ are functions of lymphatic inertia (proportional to $l{s}^{-1}$), Bernoulli resistance (proportional to $s^{-2}$), and viscose resistance to flow (proportional to $l{s}^{-2}$), respectively.

In [22], the valve resistance $R_v$$\left(R=\frac{p}{q}\right)$ changes in time, depending on pressure gradient in the lymphangion:

$$R_v=R_{V_{min}}+R_{V_{max}}\left(\frac{1}{1+exp\left(s_{open}\left(\Delta p-p_{open}\right)\right)}+\frac{1}{1+exp\left(-s_{fail}\left(\Delta p-p_{fail}\right)\right)}-1\right).$$

(5)

The transition from the maximal resistance $R_{V_{max}}$ to the minimal $R_{V_{min}}$ occurs when pressure difference $\Delta p$ before and after the valve reaches some threshold value $p_{open}$ is sufficient to open a valve. The second sigmoid describes the prolapse effect at the large adverse pressure gradient $p_{fail}$. In this formulation, resistance to backward flow has a large but limited value.

Another way to model the valve influence on lymph flow is to take it into account in the right-hand part of Equation (2). The valve provides some force influencing the flow in one direction but not in the opposite direction. This approach is used in [52] to model valve functioning in blood vessels. In [40], the valves in trunks and ducts are modelled by Equation (3). The valve influence in the smaller vessels is modelled by the following force in the right-hand side of the momentum equation:

$$f=-8\pi \nu \left(u\right)\frac{u}{s},$$

(6)

where $\nu \left(u\right)$ is a viscosity coefficient, considered as a function of velocity u having a sigmoidal form. This function characterizes friction force acting on flow in the backward direction. The sigmoidal form of this coefficient describes valve resistance to the upcoming flow. A similar approach was used in [41].

The choice of the valve model depends on the modeling task. If valve resistance to the flow is not important in the investigation, the simple boundary Condition (3) is sufficient. Formulas (4) and (5) allow for consideration of changing the valve resistance. Formula (6) is suitable for modeling multiple valves in big networks of lymphatic vessels [40]. It also allows for application of analytical methods to study lymph flow [53].

3.2. Influence of Contractions

Lymph flow in lymphatic vessels is influenced by active and passive contractions. Active contractions include a tonic part induced by intravascular pressure, and phase contractions are spontaneous and periodic (Figure 4). Active contractions are determined by lymphatic vessel wall properties and lymph chemical composition, while passive contractions occur due to some external reasons (e.g., muscle contractions or diaphragm movement). Biophysics of active contractions is not yet known. Many investigations consider the influence of different humoral factors on the amplitude, phase, and period of contractions, but a consistent theory does not exist. One more important question is about the coordination of contractions of adjacent lymphangions [17,22], since these data cannot be obtained from physiological experiments.

In early works, contractions of lymphangions were modeled using the same approach as in the models of heart pumping with time-varying elasticity [17]:

$$E\left(t\right)=\frac{P_t\left(t\right)}{V\left(t\right)-V_0},$$

(7)

where $V_0$ is dead volume (theoretical volume at zero pressure), $V\left(t\right)$ is measured volume, and $P_t$ is transmural pressure. Therefore, passive and active contractions are taken into account in this model. This model was extended in [21], where contractions of two and more (up to four) lymphangions and their coordination were investigated. In this series of works, it was demonstrated that coordination of lymphangions minimally affects lymph flow. However, later investigations showed opposite results [22]. Both works are supported by physiological experiments, and, as discussed in [23], the reason for such a difference can be related to their design, in particular, the presence of valves can be essential.

Lymphangion contractions depend on the pressure gradient. Lymphangions contract if the pressure gradient is negative, and they act as passive conduits if the pressure gradient is positive [56]. Active contractions decrease flow under a positive pressure gradient [17,18]. Therefore, under these conditions, more effective drainage occurs without contractions. This result can bring a new understanding of lymphedema pathogenesis and treatment [17]. In particular, the reduction of pumping activity in some edemas can facilitate drainage of excessive fluid [17], while it is conventionally considered a pathogenesis of the disease.

Both components of active contractions are considered together in (7) [17,21]. In other models, they can be singled out in the pressure–area relationship closing Systems (1) and (2). They are considered below in Section 3.2.1 and Section 3.2.2, and passive contractions are considered in Section 3.2.3.

3.2.1. Modelling of Tonic Response

Tonic response determines how the vessel wall reacts to changing pressure inside the lumen. It is reported in [31] that elastic properties of lymphatic vessels are close to elastic properties of veins. More precise characteristic of lymph vessels is not yet determined.

In the quasi-one-dimensional approach, elastic properties of the vessel wall are usually taken into account in the third equation (also called state equation, tube law, or pressure–area relationship), which closes Systems (1) and (2) and describes dependence of the cross-sectional area on the intraluminal pressure:

$$s=s\left(p\right).$$

(8)

Despite the existence of models concerning the cell structure of blood vessels, in the majority of investigations, empirical curves for the s–p relationship are used [57]. Similar representations are also used for all existing models of lymph flow on the basis of physiological experiments [18,22,37,39]. Despite the similarity of s–p curves for different pressure–area dependencies, mathematical solutions (in particular, the velocity of small perturbations) can differ significantly [57]. Let us consider implementations of tonic reactions in different models of lymph flow in lymphatic vessels.

The simplest pressure–area relationship with a linear dependence is used in [40]:

$$s=s\left(p\right)=s_0+\theta (p-p_0),$$

(9)

where $s_0$ is the minimum cross-sectional area at pressure $p_0$ and where $\theta$ is the coefficient characterising vessel elasticity. A nonlinear pressure–area relationship

$$p=p_{ext}+\frac{h}{r}\left(\frac{4E}{3}\frac{(r-r_0)}{r_0}+{\sigma }_{act}\right),$$

(10)

is considered in [37]. Here, h is the lymphangion wall thickness, r is the time-dependent radius, $r_0$ is the resting radius, and E is the elasticity coefficient (Young’s module). Let us note that E in (10) has an opposite meaning compared to $\theta$ in (9): for larger E, the wall is more rigid. The tonic part is represented by the first term in brackets, while ${\sigma }_{act}$ characterises phase contractions and $p_{ext}$ is external passive pressure considered as a given function. In this formulation of the pressure–area relationship, the assumption of the thin vessel wall is imposed, i.e., wall thickness is supposed to be much smaller than lymphangion radius. A thick wall model is considered in [18]:

$$\Delta p=E\Delta r_{out}\frac{r_{out}^2-r_{in}^2}{2(1-{\sigma }^2)r_{in}^2{r}_{out}}-2\pi rT\frac{{\partial }^{2}r}{\partial x^2}+2\pi r\gamma \frac{\partial r}{\partial t},$$

(11)

where $\Delta p=p-p_e$ is the transmural pressure (gradient between pressure inside of the vessel p and outside of the vessel $p_e$), $r_{in}$ and $r_{out}$ are the internal and external vessel radius respectively, $\Delta r_{out}$ is the radius variation caused by the pressure difference, and $\sigma$ is the Poisson ratio. This formulation takes into account tube tension T resulting from longitudinal bending of the vessel and the damping term with coefficient $\gamma$. The damping term $\gamma$ serves as a natural regularization for the explicit numerical scheme [18], helping to avoid the numerical instability observed in [45].

Another dependence was used in [22,55]:

$$\Delta p=p_d\left(exp\left(d/d_d\right)-{\left(\frac{d_d}{d}\right)}^3\right),$$

(12)

where d is diameter, and $p_d$ and $d_d$ are some parameters. The exponent describes wall-stiffening with strain at positive transmural pressure, and the cubic term describes a progressively diminishing vessel compliance at negative transmural pressure. Two more refined forms of this relation are used by the same authors in [58]. A similar approach with more parameters is used in [39]. These parameters are determined by fitting the data from physiological experiments.

These models describe the elastic properties of lymphatic vessels. The parameters of Equation (9) can be directly measured from the physiological experiments (for s linearly depending on p). In (10) and (11), the additional model for Young’s module E is required [59]. Formula (12) closely represents the considered physiological experiment. However, as discussed in [57], the same s–p dependence can give different solutions depending on the exact s–p function. Since the level of uncertainty in parameter determination for lymph flow is very high, a more complex model does not necessarily work better. Therefore, the final choice of the model depends on the modelling task.

3.2.2. Phase Contractions

Phase contractions represent another type of active contraction. They provide lymph propagation in the low or inverse pressure gradient between the interstitial space and the upper vena cava. Lymphedema is associated with failure of these contractions. The precise mechanism of phase contractions is not known. It is generally recognized that phase contractions are periodic and have relaxation and contraction times and that the amplitude of contractions seems to depend on the lymph volume inside the lymphangion just before contraction (diastolic volume). During the contraction phase, the lymphangion diameter can decrease up to 40% compared to its relaxed state [60]. The lymphangion length also decreases in the contraction, but existing mathematical models do not take it into account because of its small contribution compared to the diameter decrease. Furthermore, in some investigations, valves are reported to actively participate in contraction activity [31]. In the existing mathematical models, the valves are considered as a passive mechanism restricting backward flow.

It is not known whether the contractions can produce sufficient force to provide lymph flow against the gravity force in low or inverse pressure gradients between the interstitial space and upper vena cava. The models of lymph flow in the series of contracting lymphangions and in the whole LS are developed to answer this question [17,22,37,38,61]. In [61], and this question is answered positively with the developed mathematical model (Section 4). However, the evaluated parameters exceed the physiological values reported in [31,62].

A 3D model for different types of phase contractions shows that lymph flow velocity can be approximated by Poiseuille profile [44].

Since the precise mechanism of phase contractions is not known, they are modelled by some periodic, usually trigonometric, function. For example, active contractions in [40] have the following form:

$$p_{ph}=Asin\left(\frac{2\pi }{\lambda }\left(x-at\right)\right),$$

(13)

where $p_{ph}$ is the phase part of active contractions to be added in Equation (9) and where A, $\lambda$, and a are the parameters of a contraction wave: amplitude, wave length, and wave velocity, respectively.

An analytical estimate of pumping efficiency in the vessel with valves described by (6) is presented in [53]. It is shown that increased frequency of contractions leads to an increase in output flux. The same is true for increased amplitude. These results are extended in [63], where lymph flow is considered in a detailed graph of the LS (Figure 5c).

The trigonometric function is also used in [37], taking into account the current volume of lymphs in lymphangions and relax time: if these thresholds are not reached, contraction does not happen; otherwise, contractions are modelled by sinus function with an amplitude depending on lymphangion volume. A threshold value for the pressure to activate a fixed impulse is also used in [38], and this is a rare example of the work where lymphangion contraction duality is taken into account. Namely, contractions under negative pressure gradient are triggered if internal pressure in the lymphangion rises above a pressure threshold value, and it acts as a simple conduit under positive pressure gradient. This model reproduces a pulsating flow, and it can be used for the construction of big vessel networks.

Both sinus waves, simple like in (13) and with time delay like in [37], are used in [18], showing good correlation with the results of physiological experiments. It is supposed that contraction waves in lymphatic vessel networks propagate in the backward direction, from downstream to upstream, and that this behavior should be more efficient. Backward propagation (against the valves) of contraction waves was recorded in these experiments [18]. Cosine waves are used in [22] for every segment in the chain of lymphangions, and their propagation is investigated. All models reproduce the pulsating nature of lymph flow.

Phase contractions with switching from the minimal to the maximal values

$$p_{ph}\left(\frac{s}{s_0},\kappa \right)=\kappa \left(K_{max}-K_{min}\right)\psi \left(s,s_0\right),$$

(14)

in tube law $p=p\left(s\right)=p_p\left(s\right)+p_{ph}\left(s\right)+p_e$ are modelled in [39]. Tonic $p_p$, phase $p_{ph}$, and passive $p_e$ components are taken into account. Parameter $\kappa \in [0,1]$ is a state of contractions (designated by s in original work), which is a solution of the Electro-Fluid-Mechanical contraction (EFMC) model based on the FitzHugh–Nagumo model for action potentials [64], and $\psi$ is a function of a cross-sectional area similar to (12).

In this model, lymphatic vessel contraction activation is similar to heart pacemakers. According to this assumption, the FitzHugh–Nagumo model is modified and used to describe calcium-dependent contractions. Modeling results are in agreement with the experimental trends. To our knowledge, this is the first model of lymph flow dynamics where humoral factors influencing contractions are taken into account. The EFMC model opens new perspectives in the study of humoral influence on the lymph flow contractions in models of lymph flow dynamics.

3.2.3. Passive Contractions

Passive contractions are determined by external factors, such as muscle contractions, massage, diaphragm movement, etc. They are independent of the characteristics of lymphatic vessels and of the lymph chemical composition. Passive contractions are usually modeled as some function of pressure outside the vessel wall, and they are taken into account in $p_e$ ($\Delta p=p-p_e$, see (11)). Consideration of passive contractions can be essential in models of systemic lymph flow, e.g., in [40], because big muscles, especially of the lower limbs, and diaphragm excursions can have a significant effect on lymph flow. In models of contractions in one lymphangion or a series of lymphangions, the function of external pressure corresponds to the pressure inside the bath containing lymph vessel segment in in vitro experiments.

4. Systemic Lymph Flow and Graphs of the LS

The LS is highly unstructured, and it has almost no tree structure (Figure 3). As a consequence, the construction of graphs of the LS is quite difficult. However, there are several works where this attempt was implemented.

A graph of the LS was constructed in [45] (Figure 5a), and to our knowledge, it is the first reported attempt to consider systemic lymph flow. This graph starts from the interstitial space in several organs and contains feeders of the thoracic duct and the thoracic duct itself, while the right lymphatic duct is not included. It has 296 lymphangions, some of which represent lymph nodes. Lymph flow in the lymphangions is described with the zero-dimensional model from [37], and the pressure equilibrium and flux conservation are imposed in the bifurcation points. This model gives a description of fluid exchange in the interstitial space with equations similar to those describing flow in lymphangions but including flux from the interstitial space instead of flux from the previous lymphangion. Flux in the interstitial space is given by a model of Starling law for fluid exchange, and it takes into account hydraulic and osmotic pressure. The model for tissue fluid exchange was further developed in [62] with the one-dimensional approach, including the influence of anchoring filaments and primary valves in the initial lymphatics. This systemic model shows the pulsating nature of lymph flow from the thoracic duct that is consistent with physiological observations. The parameters of the flow were in agreement with available physiological data.

The algorithm of automatic construction of a graph of the LS is suggested in [65]. The data from the Plastic Boy project [66] is used as a source of anatomical structure. These data are not suitable for calculations as-is, and the constructed algorithm turns this image into a consistent graph of the LS (Figure 5b). Structural analysis of the LS is provided, including the average degree of lymph nodes (how many vessels connected to the node). This algorithm can be applied to other anatomical models. The presence of a graph of the lymphatic system makes it possible to analyze the redistribution of lymph flow in the human lymphatic system in case of damage and/or difficulty in the outflow of lymphs from individual lymph nodes, which will allow for prediction of the effects of various injuries and changes in the organs of the human lymphatic system as well as for prediction of the best methods for correcting these problems.

One more graph of the LS is constructed in [40] (Figure 5c). It is anatomically adequate, spatially oriented, and topology consistent with a similar graph of the cardiovascular system [47]. The graph contains 271 collecting lymphatic vessels, thousands of lymphangions, and 161 lymph nodes. Groups of lymph nodes form the basis of the graph, and it contains 46 main groups [54]. Flow in the vessels is described by the quasi-one-dimensional equations including valve and active contraction influence, and the conditions of pressure equilibrium and mass conservation are imposed in the bifurcation points. The model of active pumping in this graph and its influence on the total flow are investigated without the influence of gravity force [63] and with influence of gravity force [61]. The increase in amplitude and frequency of phase contractions is shown to increase the output flux. Also, the elasticity and contractions of lymph nodes have a considerable effect on the output flux [63]. Therefore, the lymph nodes affect the transport function of the LS, and it is consistent with physiological observations [32]. The frequency of contractions, required to get output flux in physiology adequate intervals, is higher than reported in the physiological experiments [31,62]; however, in vivo parameters can also differ from them [15]. The consideration of passive contractions, which are not included in the model now, may allow for a reduction in the frequency. The results show that active contractions can produce enough force to make lymphs flow though the whole LS against the pressure gradient [61].

The approach used to construct a spatially oriented graph in [65] requires source images to contain information about both lymph nodes and lymphatic vessels. Another approach is used to create a graph in [40]. The basis of this graph is the regional groups of lymph nodes which connect according to the communication matrix filled in with data from [29]. This approach can be used to construct personalized graphs of the LS by extracting information about lymph nodes from MRI images and by connecting them with vessels using the communication matrix according to some anatomical source (i.e., [29]). The described concept could be the first step towards the creation of personalized graphs of the LS.

5. Lymph Drainage and Flow in Initial Lymphatics

The LS starts from the initial lymphatics (lymphatic capillaries, terminal lymphatics), where lymphs arrive from the interstitial space. Characteristic features of lymph drainage from the interstitial space are the influence of primary valves, which allow lymphs to enter the LS and which stop it from flowing out of the lymphatic capillaries, and the presence of anchoring filaments, which prevent initial lymphatics from collapse even in excess pressure in the interstitial fluid (Figure 1B).

In contrast to models of secondary valves, models of the primary valves often take into account the geometry of the valve—the curvature of the cell [67,68].

Flow into the initial lymphatics can be described by a system of ordinary differential equations, including flux and pressure balance and taking into account osmotic pressure according to Starling’s hypothesis [45]. Another model takes into account the primary valves influence given by Condition (3) and anchoring filaments with additional force in the pressure–area relationship similar to (10) [62]. This model describes pumping of income lymph flow in the LS due to fluctuations in the interstitial fluid pressure and to the suction effect of adjacent pumping lymphatic vessels.

Primary valve dynamics can also be described by the Stockes equation [67]:

$$\frac{\mathrm{d}p}{\mathrm{d}x}=-\frac{12Q\mu }{\omega {\left(x\right)}^{3}T},$$

(15)

where T is a unit tissue thickness, $p\left(x\right)$ is pressure inside the fluid gap between the cells, Q is flux through the gap, and $\mu$ is dynamic viscosity. This equation must be solved simultaneously with the equation for cell deflection $\omega \left(x\right)$:

$$\frac{{\mathrm{d}}^4\omega \left(x\right)}{\mathrm{d}{x}^4}=\frac{1}{EI}\Delta p\left(x\right),$$

(16)

where $\Delta p=p\left(x\right)-p_L$ is transendothelial pressure drop, $p_L$ is pressure inside the lumen, and E is Young’s module. Flux through the gap is given by the following expression:

$$Q={\int }_0^{\omega \left(x\right)}\mathbf{u}\left(\mathbf{x}\right)\mathrm{d}y,$$

(17)

where $\mathbf{u}\left(\mathbf{x}\right)$ is flow velocity. Simulation results demonstrate primary valve dynamics for this model: the flux nonlinearly depends on pressure p with positive pressure gradient $\Delta p$, and it equals 0 for negative pressure gradient.

This model was extended in [69], where a more precise description of flow around the overlapping endothelial cells is considered in the same geometry. The computational results correspond to the available physiological data.

The geometrical structure of the interstitial space with initial lymphatics and blood capillaries is considered in [68], and a mechanical model for primary valve dynamics concerning their curvatures is used to describe flow through the junctions between endothelial cells. As a result of this two-dimensional analysis, the formula for lymph flux per unit length q in terms of the pressure difference between blood $\widehat{p}$ and lymphatic capillaries $p_0$ was obtained:

$$q=\left\{\begin{array}{c}\begin{array}{cccc}\hfill & 0,\hfill & & \widehat{p}<p_0+p_{crit}\\ \hfill & \frac{\xi \alpha \pi r_b}{\alpha \pi r_b+2\xi }\frac{k}{\mu }\left(\widehat{p}-\left(p_0+p_{crit}\right)\right),\hfill & & \widehat{p}>p_0+p_{crit}\end{array}\hfill \end{array}\right.,$$

(18)

where $p_{crit}$ is the critical pressure needed to open a valve, $\xi$ is a dimensionless parameter that describes the effect of geometry on the fluid flux, k is permeability of intestitium, $\mu$ is viscosity of interstitial fluid, and $r_b$ is radius of the blood capillary. This model describes fluid exchange between blood and lymphatic capillaries through the interstitial space.

All considered models describe periodical entering of lymphs into the initial lymphatics. They provide more precise boundary conditions for models of flow in the lymphangions, and they can be used to describe lymph drainage in the investigations of lymphedema pathology.

6. The Structure and Function of the Lymph Node

Lymph nodes (LNs), which represent an important part of the lymphatic system, in addition to merging lymph flows, perform a number of functions related to the human immune system. This is due to the fact that they contain localized populations of cells that perform protective functions in the body. Thus, populations of T-lymphocytes are located in the so-called T-cell zone and B-lymphocytes are grouped in B-cell follicles, which are involved in the production of antibodies specific for the antigens entering the body. Dendritic and fibroblastic reticular cells are also present in the T-cell zone and are involved in all aspects of the functioning of lymph nodes. Figure 6 shows the key elements of the structure of the lymph node. These data are provided according to [70,71,72].

In the lymph node, T and B lymphocytes are in an inactive, naive state. The populations of these cells are regulated by fibroblastic reticular cells (FRC) secreting interleukin-7, which is a survival factor for naive lymphocytes [73]. The FRC network is also involved in maintaining the structural integrity of the entire organ through threads formed by polymerized collagen (Figure 6) secreted by fibroblastic reticular cells. The process of immune response in LNs begins with the helper T lymphocyte that has detected the antigen or with another activated antigen-presenting cell [74]. In this case, a cascade process occurs: T helper lymphocytes start activation of resting antigen-presenting cells, T effector lymphocytes and B-lymphocytes, and interleukin-2 is produced to increase the proliferative activity of T and B-lymphocytes [75]. Activated antigen-presenting cells (APC) fix onto the FRC network for a period of about 6 hours with the possibility of subsequent reactivation. After that, APCs begin to express interferon-$\alpha$. Under the influence of interferon-$\alpha$, naive T effector lymphocytes specialize according to the components of the antigen, presented by the antigen-presenting cells’ major histocompatibility complex class I, II (MHC class I, II) [76]. At the same time, the FRC network is actively involved in stimulating immune response, producing chemokines CCL19 and CCL21. According to existing data, the lymphocytes interacting with these chemokines exhibit increased motor activity. This effect does not improve immune response alone, but in synergy with the chemotaxis effect used by cells to search for foci of infection, chemokines shortened the time from activation to interaction with the target. Activated B-lymphocytes can also activate resting cells, but their main function is the production of antibodies specific for the target antigen [74]. At the same time, activated T effector lymphocytes search for and destroy cells expressing signals characteristic of the desired antigen. T effector lymphocytes mainly destroy body cells infected with intracellular parasites (viruses or bacteria). Macrophages are also involved in the process of the immune response, destroying free viruses and bacteria. After activation in the lymph nodes, cells of the immune system move into tissues in which antigens are present to fulfill functions of the immune response. In connection with such a heterogeneous structure of the lymph nodes, it becomes extremely important to understand how the lymph flows into the node are distributed ([77,78]). The size of particles that can enter different zones of the lymph node is limited [77]. Thus, it turns out that huge viruses and bacteria cannot penetrate into the T-cell zone through the walls of the FRC network due to too large sizes. Therefore, they get there either through blood vessels or inside the cells that came from outside the lymph node. It should be noted that the detection of antigens is difficult not only because of the limited methods of their entry into the T-cell zone. About 93% of lymphs generally passes through the subcapsular sinus of the LN. This is shown, inter alia, in [78]. It should also be taken into account that part of the lymph from the lymph node enters the blood vessels system inside LNs.

Due to the complexity of the processes occurring in the lymph nodes and the heterogeneity of the lymph nodes, it becomes necessary to be able to form various configurations of the lymph nodes to simulate the processes of cell interactions, fluid flow, and collagen coagulation. Since it is expensive and time-consuming to restore structures using layer-by-layer scanning, an important area of research is the development of algorithmic methods for constructing lymph nodes. Algorithms were developed, and 3D models of the components of the lymph node were built in [79,80]. Three-dimensional models of key elements of LNs in a voxel format were created with further conversion to a solid-state 3D model: the shell of the lymph node, trabecular and medullary sinuses, B-cell follicles, and a network of fibroblastic reticular cells.

Separately, algorithms were developed to build a 3D model of a network of fibroblastic reticular cells [80]. It appears that a highly damaged network (up to 50%) can still support lymph flow. This is important for the case of HIV infection where the structure of the fibroblastic reticular cells as well as the ducts formed by them are gradually destroyed. Damage of the network significantly affects not only the ability of the T-cell zone of the lymph node to pass lymphs but also the homeostasis of cell populations and the immune response to HIV in general [80,81].

7. Conclusions and Perspectives

The main goal of investigations of lymph flow is to understand the underlying mechanisms in order to propose effective therapeutic techniques. In the 90s, the specific lymphatic factor (VEGF–C) was discovered, and it led to the explosive investigation of LS [82]. The mechanisms of lymph flow are more clear now, but there are still many open questions.

Mathematical modeling is used with in vitro experiments to get more data, to check the hypothesis about lymph flow, and to construct new experiments. In [17,21], the authors established conditions for pumping and conduit acting. In negative pressure difference, the lymphangion acts as an active pump, and in positive pressure, it acts as a simple conduit. Moreover, in large positive pressure gradients, active contractions are shown to inhibit lymph flow. This result is supported by a physiological experiment [18,56]. The meaning of this result in the application to lymphedema pathogenesis and treatment is discussed in [17]. For some edemas, the inhibition of lymphangion pumping may have a positive effect.

Pathogenesis and treatment of lymphedema are one of the most important applications of models of lymph flow in lymphatic vessels. Direct models of lymph flow in lymphatic vessels are not the only way to address edema pathology: in [83], lymph filtration through the extremity was considered with diffusion equations instead of considering flow in separate vessels, and some clinical recommendations (for using of bandages) were proposed.

Another important goal of lymph flow modelling is the investigation of drug distribution in the LS and various tissues and organs. Insufficient lymphatic circulation is one of possible reasons for ineffective drug therapy, and such a model for systemic circulation can find its applications not only for lymphedema treatment but also for other drug-delivery tasks. Lymph drainage from the interstitial fluid and sorption effect of the secondary lymphatic circulation (vessels with contractions and secondary valves) are essential for understanding the mechanisms of more efficient drainage in the cases of different pathologies. The proposed concept for construction of personalized graphs of the LS could be the first step towards achieving a personalized model of human LS. It opens new perspectives in these investigations.

Lymph flow in the lymph node may influence the immune response occurring in them. Modeling the processes of cellular interactions in the lymph nodes as well as the processes of transfer between nodes through the lymphatic system are critical for understanding the processes that occur in HIV-infected people. It gives an opportunity to develop the best strategies for drug therapy, to analyze the effectiveness of existing methods, and to find key relationships in complex immune response processes. The human immune system depends on many populations of cells represented in the T-cell zone of the lymph nodes: T-lymphocytes (helpers and effectors); B-lymphocytes, which produce antibodies when pathogens enter the body; macrophages; and antigen-presenting cells. Certain cell populations, such as fibroblast reticular cells, are involved in maintaining the homeostasis of naive lymphocyte populations, and the whole lymphatic system provides a reserve of bandwidth for the fluid circulation in the body and acts as another filter for intercellular fluid collected by lymphatic capillaries from tissues.

There are still many open questions in lymph flow mechanics, but close cooperation between physiologists and mathematicians will help understand it better.