Modelling Lymph Propulsion in a 3D Model of Murine Collecting Vessel with Three Lymphangions in Series

Abstract

The lymphatic system maintains tissue homeostasis by transporting the excess fluid from the interstitium and ultimately returning it to the venous circulation against an adverse pressure gradient and gravitational force. The spontaneous contractions of lymphangions, the building blocks of collecting vessels, and the secondary lymphatic valves play key roles in lymph propulsion. The aim of this study was to investigate lymph propulsion in a series of three contracting lymphangions in a 3D reconstructed model segmented from micro-CT scans of the collecting lymphatics in the hind limb of mice. Computational Fluid Dynamics and Fluid–Structure Interaction were used to study the behavior of flow within the collecting vessel, as well as the behavior and deformations of the vessel wall and the poroelastic interstitium. The secondary valves were modelled as porous membranes with closed or open states depending on their permeability. A sensitivity study revealed that the parameters having the most impact on the total volume of lymph propelled by active contraction of the lymphangions were the elastic modulus of the interstitium and the permeability of the secondary valves during the open states.

1 Introduction

The lymphatic system is a tissue drainage network of converging vessels that maintains the tissue fluid balance by transporting the excess fluid, proteins, and metabolic waste products from the interstitium and ultimately returning it to the venous circulation [1]. The lymphatic system also acts as a conduit for immune cells and facilitates the immune response. The main components of the lymphatic route are initial lymphatics, the collecting lymphatics, and the lymph nodes [2]. The initial lymphatics are thin-walled blunt vessels composed of a single layer of endothelial cells serving as the site where the interstitial fluid (IF) enters the lymphatic system (a process referred to as lymph formation or lymph uptake) [3]. To return to the venous circulation, lymph travels against an adverse pressure gradient and gravitational forces. Hence, lymphatic transport is accomplished as a result of the extrinsic and intrinsic pumping mechanisms that passively and actively cause the compression of the lymphatic vessel walls, as well as the primary and secondary lymphatic valves that inhibit backflow [4]. The deformation of tissues surrounding lymphatic vessels, due to muscle contractions, arterial pulsations, etc., contributes to the extrinsic (passive) pumping mechanism, while the intrinsic spontaneous contractions in the collecting lymphatic vessels actively propel lymph forward.

Collecting lymphatics are comprised of easily identifiable segments, called lymphangions, which are separated by secondary lymphatic valves (Fig. 1). Lymphangions are the contractile components of the lymphatic system that actively transport the lymph along the lymphatic system (i.e., lymphatic propulsion).

Fig. 1

The anatomy of collecting vessels. A The schematic anatomy of a collecting vessel with lymphangions being the segments in between two secondary valves. B A snapshot of the collecting vessel from the hind limb of a Prox¹-GFP mouse under a fluorescence microscope, courtesy of Dr. Mike Davis

Impairment of lymph propulsion leads to a variety of pathologies including lymphedema. Lymphedema is the general term used to refer to a set of pathologic conditions causing excessive accumulation of protein-rich fluid in the interstitium [5]. Lymphedema can be a disfiguring condition, usually affecting the extremities, causing pain and discomfort, tissue fibrosis, and high susceptibility to infections which can require hospitalization. Potential causes include the removal of lymph nodes during a surgical procedure for breast cancer treatment, infections, or congenital problems during the development of the lymphatic system [6]. Breast cancer survivors face up to 40% risk of lymphedema occurrence [7]. Breast cancer related lymphedema may show immediately or even years after the treatment [8]. To our knowledge, there is no proven drug treatment for lymphedema. Currently, its management aims at reducing or delaying progression of the swelling and preventing its associated infections. In severe stages of lymphedema, the composition and organization of the interstitium may drastically change due to the degradation of the extracellular matrix leading to changes in its mechanical properties, such as its hydraulic permeability and Young’s modulus [9]. Lack of a proper understanding of the fluid transport mechanisms during lymphedema has partly led to ineffective treatments.

Mathematical models may provide a flexible and powerful gateway to understanding lymphatic fluid transport and complement ex vivo and in vivo experimental research in animal models and in humans. Drake et al. [10] were the first to apply equivalent electric circuit techniques to the lymphatic system with each equivalent circuit consisting of a single resistor in series with a single pressure source, while the valves were modelled as diodes. However, the model did not account for the active contractions of the lymphangions. Quick and Venugopal published several papers [11,12,13,14,15] modeling the lymphatic system with a more elaborated model in which the parameters were time-dependent. Bertram et al. [16] used a more sophisticated model for the secondary valves in which a pressure-dependent resistance was incorporated instead of a simple diode. Bertram and colleagues’ multi-lymphangion model consists of multiple actively contracting segments in between secondary valves. The same group then incorporated realistic and measured parameter values in collaboration with Davis [17]. Jamalian et al. [18] performed a parametric study on a lumped parameter model of a chain of lymphangions in series. They showed that the model was highly sensitive to the minimum valve resistance, which could cause the time averaged flow rate to change by several orders of magnitude. Later, Jamalian et al. [19] and Bertram et al. [20] used a more realistic formulation of secondary valve properties based on experimental measurements and concluded that sequential contractions with an appropriate time delay in between would produce a more effective pumping than synchronized contractions. In Ikhimwin et al. [21], that lumped parameter model was further extended to be used for a network of initial lymphatics and pre-collectors with sparse secondary valves with the network being surrounded by a permeable interstitium, allowing to study the interplay between primary and secondary valve function.

Lumped parameter models as described above are computationally efficient and can be easily applied to large and branched networks of collecting vessels. However, they fail to include the spatial variations of different parameters. That limitation is partially resolved when solving the 1D formulation of the Navier–Stokes equations, yielding 1D models of lymphangions that provide more fundamental results [22]. Reddy et al. [23] created a 1D model of the collecting lymphatics in which a set of equations was discretized using differencing techniques with each lymphangion being considered as one grid point. Ultimately, a lumped parameter model was obtained here, but based on the Navier–Stokes equations. MacDonald developed a more refined 1D model starting from the Navier–Stokes equations, modelling either a single lymphangion [24], or a short series of lymphangions [25]. Contarino and Toro [26] developed a 1D model for collecting vessels where the electrical activity of the lymphatic wall was based on the FitzHugh–Nagumo model, coupled to the vessel wall mechanics, and a lumped parameter model for valve dynamics. In this model, the contractions were initiated depending on environmental and stretch-activated calcium influx and inhibited depending on the wall shear stresses. However, the adjacent lymphangions were not electrically coupled.

MacDonald et al. [24] modelled a single 2D bi-leaflet secondary valve where the valve’s motion was imposed rather than being solved by fluid–structure interaction (FSI). Rahbar and Moore [27] published a model describing the flow in 3D in a straight cylindrical tube with a contracting segment, modelling a single lymphangion with static inlet and outlet segments. The secondary valves were not included in the model and the FSI was implemented by imposing the wall motion as a boundary condition. Wilson et al. [28] modelled a 3D idealized secondary valve based on dimensions from confocal images using an uncoupled approach. These authors [29] then proceeded with a two-way coupled FSI approach on the same model of the 3D idealized secondary valve. They were able to show the eddy structure in the sinus of the valve as well as the influence of the valve leaflets on the resistance to forward flow, albeit the vessel wall was assumed rigid. Ballard et al. [30] also developed a fully coupled 3D computational model of a secondary valve with idealized geometry using a lattice Boltzmann model to simulate the lymph and a lattice spring model for the secondary valve. Similar to Wilson et al. [29], the vessel wall was assumed rigid in Ballard et al.’s steady-state simulation. Bertram [31] also modelled a 3D secondary valve using FSI with both the vessel wall and the valve being deformable. While contractions were not incorporated in this 3D model, the impacts of both the vessel wall and the leaflets’ stiffness on the valves’ behavior were shown.

Higher dimensional models of not only the lymphatic collecting vessels but the lymphatic system in general are limited. The models of the collecting vessels have often relied on simplified geometries that do not fully capture the anatomical characteristics observed in medical imaging [24, 27,28,29,30,31]. Additionally, the previous 3D models of collecting vessels have primarily focused on an isolated lymphangion without considering its surrounding tissue [24, 27,28,29,30,31]. Other simplifications noticed in some of the 3D computational models include assumptions on the collecting vessel wall to be rigid [29, 30], not incorporating the contractions of the lymphangions, and not studying the collective behavior of multiple lymphangions [27,28,29,30,31].

With that in mind, the aim of this study was to investigate lymph propulsion in a series of three contracting lymphangions in a 3D model reconstructed from micro-CT scans of the collecting lymphatics in the hind limb of mice. We used Computational Fluid Dynamics (CFD) and FSI to study the behavior of flow within the collecting vessel, as well as the behavior and the deformations of the vessel wall and the interstitium modelled as a poroelastic medium. After developing the baseline model, we also performed a parametric study to evaluate their effect on the transport of lymph.

2 Materials and Methods

2.1 Animal Experiments

The ex vivo animal experiments were performed to collect morphological information from the collecting lymphatic vessels in mice limbs. The morphology of the lymphatic system is difficult to capture due to its complex nature. Lymph is transparent and due to the presence of the secondary valves, contrast agents can only be introduced into the lymphatics via antegrade interstitial injections, which makes imaging of the lymphatic system even more difficult. The animal experiments were approved by the Ethical Committee of the Faculty of Medicine and Health Sciences of Ghent University Hospital (ECD 19/25) and were conducted with respect to animal welfare. A total of 4 C57BL/6 male mice aged between 15 and 17 weeks (young adults) were sacrificed for these experiments, where the mice were first euthanized via ${\text{CO}}_{2}$ overdose, before being prepared for micro-CT imaging. To visualize the lymphatic vessels in the micro-CT scanner, a contrast agent was injected into the lymphatic vessels via antegrade interstitial injections. An appropriate contrast agent for reaching the lymphatic vessels through interstitial injections should only be taken up by the lymphatic capillaries and should not leak. The contrast agent used for this study was Iopromide (Ultravist-300®, Bayer Inc.), which is a commercially available, water-soluble iodinated contrast medium. By means of a stereomicroscope, approximately 10 µl of the contrast agent was manually injected into the toe skin using a 29-gauge needle and a 500 µl syringe. Care was taken to ensure that the injections were performed superficially, into the dermis layer, since that is where the initial lymphatics are more dominantly present [32]. The possibility that the extra-interstitial pressure due to the injection could lead to more lymph production and expansion of the vessel should also be noted. Immediately after injection, micro-CT scans were taken using a High-Energy CT system Optimized for Research (HECTOR, UGent). This high energy CT scanner is equipped with a 240 kV X-ray tube from X-RAY WorX, a PerkinElmer 1620 flat-panel detector, and a rotation stage on which the sample is mounted [33]. This system provided a resolution of 40 µm for this specific sample.

2.2 Segmentation and Mesh Generation

The images acquired from the aforementioned micro-CT scanning were imported into the commercially available software Mimics Research® (Materialise, Leuven, Belgium) where they were segmented, 3D reconstructed, and smoothened. A section of the collecting vessel including three lymphangions was located at the top of the foot region of the mouse’s right hind limb, as shown in Fig. 2. After improving the surface mesh of the segment of the vessel in 3-Matic (®Materialise, Leuven, Belgium), the .stl file was imported into ICEM (Ansys Inc., Canonsburg, PA, USA) to create the volume mesh resulting in 63,843 volume elements. The volume mesh was then imported into COMSOL Multiphysics 5.6 where the interstitium with a radius of 500 µm was also added to the computational domain (Fig. 2F). The lengths of the entire collecting vessel, lymphangion 1, lymphangion 2, and lymphangion 3 were 1.83 mm, 0.41 mm, 0.67 mm, and 0.31 mm, respectively. Cross-sections are non-circular, but the equivalent diameter of the collecting vessel ranged from 130 to 50 µm. The locations of the valves and lymphangions were estimated based on the vessel contours and the assumption that usually an enlargement is seen at the location of the secondary valves. Such enlargement can also be seen in Fig. 1B at the location of the secondary valve.

Fig. 2

From micro-CT images to the 3D reconstructed collecting vessel. The raw micro-CT images resulting from imaging the mouse (A), the segmented and 3D reconstructed collecting vessel that leads to the popliteal lymph node in the right hind limb of the mouse (B), the region of interest at the top of the foot (C), the smoothened and chosen segment of the vessel used to create the 3D reconstructed model (D), a cross-section of the vessel in the valve region (E), and finally, the dimensions and geometry of the final 3D reconstructed model (F)

For the mesh sensitivity analysis, four different mesh densities were created with around 76,000, 137,000, 178,000, and 350,000 volume elements and around 283,000, 528,000, 684,000, and 1,356,000 degrees of freedom, respectively. To verify the accuracy of the model, these four meshes along with coarser to finer timesteps were computed. The accuracy was evaluated by checking the conservation of mass first and then comparing the conveyed volume at the outlet for the different meshes. Ultimately, the mesh with 178,000 elements and a relative tolerance of ${10}^{-4}$ was selected (a smaller tolerance would mean smaller timesteps and higher accuracy).

3 Computational Domain

The model was built based on a realistic morphology acquired from in vitro experiments and was created using COMSOL Multiphysics® simulation software v.5.6. using the CFD, Microelectromechanical Systems (MEMS), and FSI modules. The model describes a series of three pumping lymphangions immersed in the interstitium. It includes four different domains being the inner volume of the collecting vessel, the lymphatic vessel wall, the interstitium, and the secondary valves. Lymph and the interstitial fluid were modelled as a free fluid phase with laminar flow, a density of 998 ${\text{kg/m}}^{{3}}$ and a viscosity of 1.5 ${\text{mPa}}\;{\text{s}}$. The interstitium was modelled as a poroelastic medium with a Young’s modulus of 2 kPa (for more information on poroelasticity see Supplementary file1). The interstitium can have a wide range of values for its elastic modulus depending on where it is situated in the body and in which species [34]. We chose the value of 2 kPa but also included higher values of Young’s modulus in the parametric study later on. The vessel wall was modelled as a non-permeable linear elastic material with an elastic modulus of 1.2 kPa [25]. The secondary valves were modelled as thin porous disks that followed the contours of the vessel and had an arbitrary thickness of 5 µm (Fig. 2E and F). Situated in between the lymphangions, the valves were implemented as porous media with a pressure-dependent permeability ranging between $10^{ - 11}$ m² (open state) and $10^{ - 30}$ m² (closed state) to prevent lymph backflow whenever the transvalvular pressure is higher than a threshold (for more detailed information, see Sect. 4). An overview of the mechanical properties of these domains deduced from literature is presented in Table 1.

Table 1 Material properties of the different domains of the baseline model

4 Secondary Lymphatic Valve Mechanism

To implement the one-way secondary lymphatic valves in between lymphangions, we used a sigmoidal function as proposed by Bertram et al. [16] (where the permeability κ of the secondary valves fluctuated in response to transvalvular pressure variations according to the equation below) to create the opened and closed status of the valves

$$\kappa = \frac{{10^{ - a} }}{{1 + e^{{ - s\left( {\Delta p - P_{{{\text{th}}}} } \right)}} }} + 10^{ - b} .$$

(1)

In this equation, the pressure drop Δp, in Pa, is the transvalvular pressure of the secondary valves, $P_{{{\text{th}}}}$, also in Pa, indicates the threshold at which point the valve transitions between the open and closed states, and s is the slope of the sigmoidal function. The terms $10^{ - a}$ and $10^{ - b}$ represent the maximum and minimum values of permeability serving as the open and closed states of the valves, respectively. An example of the permeability function, as well as the flow rate for a secondary valve with circular cross-section and a diameter of 100 µm, is presented in Fig. 3A.

Fig. 3

The boundary conditions. The permeability function of the secondary valve and the flow rate of a valve with a diameter of 100 µm, as a function of transvalvular pressure gradient (A). An overview of the boundary conditions with the yellow hachures representing the fixed surfaces (B), a cross-section of the collecting vessel showing valve 1 and the fixed surrounding vessel wall (C), and the smoothened piecewise boundary load applied on each lymphangion’s vessel wall for the baseline model, with a maximum magnitude of 150 Pa (D)

The slope and the range of permeability were defined in such a way to have the least amount of lymph backflow at the valves, based on the initial results, avoiding valve leakage while maintaining computational stability. Hence, the values for maximum and minimum secondary valve permeabilities of the baseline model were specified as $10^{ - a} = { }10^{ - 11}$ ${\text{m}}^{{2}}$ and $10^{ - b} = { }10^{ - 30}$ ${\text{m}}^{{2}}$, respectively. The slope was set to 20 and a pressure threshold of 0.2 Pa.

5 Boundary Conditions and Materials

The overview of all the boundary conditions is shown in Fig. 3. At the inlet and outlet of the lymph vessel, open boundaries with zero pressure were imposed, so that the lymph propulsion would only be due to the contractions of the lymphangions. To emulate contractions, a smoothened piecewise function was applied as the boundary load at the vessel wall of each lymphangion with a maximum magnitude of 150 ${\text{Pa}}$ for the baseline model. The load is applied in such a manner that the diameter of the lymphangions would reduce synchronously everywhere [5]. The lymphangions contract with a duration of 2 s; they contract one after another with a time delay (Δt) of 0.4 s in between their contractions, as shown in Fig. 3D [41]. The boundary load applied at the interstitium was set to zero for the baseline model. At the inner side of the vessel wall, no-slip boundary conditions were specified. Aside from the lymphangion walls, the rest of the collecting vessel walls were fixed. The lymphangion walls can freely move, except at the valve locations that are fixed in space, as well as the two side boundaries of the interstitium. The interstitial fluid, however, could freely pass all outer boundaries of the interstitium.

6 Parameter Study

A parametric study was performed to determine the impact of the elastic modulus of the interstitium $\left( {E_{{{\text{interstitium}}}} } \right)$, the elastic modulus of the lymphatic collecting vessel wall $\left( {E_{{{\text{vessel}}}} } \right)$, the magnitude of the lymphangions’ contraction load (${\text{BL}}_{{{\text{vessel}}}}$), the refractory period in between multiple contraction cycles (RP), the adverse pressure at the outlet ($P_{{{\text{outlet}}}}$), the load applied at the interstitial boundaries (${\text{BL}}_{{{\text{interstitium}}}}$), and the maximum permeability of the secondary valves, $10^{ - a}$ in Eq. (1) ($k_{{{\text{max}}}}$). The corresponding values of each parameter are presented in Table 2. For parameters having a wide range of values reported in the literature (such as the elastic modulus of the collecting vessel wall with values from 1200 ± 700 Pa to 10 kPa [25, 36, 40]), this range was incorporated into the parametric study. For parameters for which we could not find specific values in literature, a physiologically relevant range was chosen.

Table 2 List of parameters and their corresponding values studied

7 Results

7.1 Baseline Model

The volume averaged speed of lymph flow in the lymphangions, their volume averaged pressure, the lymphangions’ surface averaged vessel wall displacements, the variations of the valve permeabilities over time, and the conveyed volume are shown in Fig. 4.

Fig. 4

Results from the baseline model. The graphs show the volume averaged speed of lymph in lymphangions (A), the volume averaged pressure of lymphangions (B), the surface averaged vessel wall displacement of each lymphangion (C), the permeabilities of each valve (D), the volumetric flow rates passing each valve (E), and the conveyed volume (F)

Taking both the lymph velocity (Fig. 4A) and the valves’ permeability magnitudes (Fig. 4D) into consideration, it is observed that as soon as the permeability of a valve starts decreasing, the lymph speed at the valves decreases towards zero with the valve transitioning into the closed state ($10^{ - b}$ in Eq. 1, permeability k = $10^{ - 30} \; {\text{m}}\left. {^{2} } \right)$. Correlating the contraction causing boundary loads (Fig. 3D), and the valves’ permeability graphs (Fig. 4D), it can be deduced that with each contraction and as the lymphangion starts compressing, its upstream valve closes to prevent backflow, while the downstream valve remains open until the deformation of the downstream lymphangion starts.

Not only the surface averaged wall displacements of the three lymphangions, shown in Fig. 4C, were non-uniform despite having the same boundary load, but also the displacement along each single lymphangion’s wall was non-homogeneous due to the shape irregularities and varying diameters (see supplementary materials, Fig. S1).

Another parameter to take into consideration was the volume conveyed, the accumulated volume of fluid conveyed at the outlet ($V_{{{\text{outlet}}}}$), which was calculated by integrating the volumetric outflow over time. The volume conveyed for this 3D reconstructed model by the end of a single contraction cycle is shown in Fig. 4E.

By correlating the pressure and volume variations of a lymphangion throughout its contraction cycle, a pressure–volume loop (P–V loop) is obtained, as shown in Fig. 5, for lymphangion 1. Point A represents the start of systole, when the contraction of lymphangion 1 starts. Point B represents peak intralymphatic pressure. Point C represents the minimum diameter and the end of the systole. Point D represents the start of rapid filling of lymphangion 1 and the end of isovolumetric relaxation. Point E represents the onset of downstream lymphangion (lymphangion 2) relaxation. Point F represents the end of rapid filling phase and the beginning of slow filling. Ultimately, point G marks the end of the contraction cycle. Hence, points A–C represent the systole, while the rest the diastole. The ejection fraction (EF) of the lymphangions is defined as the ratio of the end-diastolic volume of lymph ejected by each lymphangion over its initial volume. The calculated EF for lymphangion 1, 2, and 3 were equal to 0.17, 0.18, and 0.17, respectively.

Fig. 5

Contraction cycle characteristics of lymphangion 1. Representative graphs of the average intralymphatic pressure, diameter, and P–V loop of lymphangion 1. Some significant points throughout this contraction cycle were chosen to showcase some significant points throughout its contraction cycle. These points (A–G) represent the start of systole (point A), the peak intralymphatic pressure (point B), the end of systole (point C), the start of rapid filling of the lymphangion 1 (point D), the onset of downstream lymphangion (lymphangion 2) relaxation (point E), the end of the rapid filling phase (point F), and finally the end of the contraction cycle (point G)

8 Parametric Study

8.1 Adverse Pressure at the Outlet ($P_{{{\text{outlet}}}}$)

As the pressure in the outlet was increased from 0 Pa in the baseline model to 1 Pa, 5 Pa, and 20 Pa, the volume conveyed decreased linearly by 0.61, 3.24, and 13.21%, compared to the conveyed volume for the baseline model (Fig. S2). The reason behind this reduction of $V_{{{\text{outlet}}}}$ can be found by looking at the behavior of the last valve, valve 4, which is in the immediate vicinity of the vessel outlet and consequently the $P_{{{\text{outlet}}}}$ (Fig. 6). As a result and to prevent backflow, whenever there is an adverse pressure at the outlet, valve 4 would stay closed until the pressure behind it equals or surmounts the pressure after it. Correspondingly, with an adverse pressure at the outlet, valve 4 stayed closed for the majority of the contraction cycle. In contrast, without the adverse pressure, valve 4 remained mainly open and the only instant it was closed was during the relaxation of lymphangion 3.

Fig. 6

The effect of adverse pressure at the outlet on the behavior of valve 4. The volume conveyed and the behavior of the last valve, valve 4, throughout the contraction cycle with 20 Pa adverse pressure at the outlet (upper row) and without any adverse pressure (bottom row)

With valve 4 staying closed for longer due to $P_{{{\text{outlet}}}}$, the volume of lymph that was pumped as a result of contractions starts to accumulate causing an increase in intralymphatic pressure (Fig. S3) and swelling in the lymphangions upstream valve 4. The swelling of lymphangion 2 and 3 starts with the initiation of the first lymphangion’s contraction and lasts until valve 4 is opened. Due to the elevated intralymphatic pressure and extra swelling in the presence of adverse pressure, a smaller wall displacement was also observed (Fig. 7).

Fig. 7

The displacements of the collecting vessel wall along its length at specific time points. The displacements along the vessel wall of the collecting vessel (magnified by a factor of 5) without and with adverse pressure at each time point. The last column shows the different time points during the contraction cycle

Figure 8A shows the pressure–volume loop of a representative lymphangion, lymphangion 1, for increasing $P_{{{\text{outlet}}}}$ showing an increase of the end-systolic volume, as well as of the end-systolic intralymphatic pressure.

Fig. 8

The P–V loops of lymphangion 1. The effect of changing pressure at the outlet (A), lymphangions’ contraction load (B), interstitial boundary load (C), refractory period (D), and maximum permeability of the secondary valve (E) on the P–V loops of a representative lymphangion, lymphangion 1. The black arrow in the legends points to the baseline model

8.2 Boundary Loads at the Lymphangions’ Vessel Walls (${\text{BL}}_{{{\text{vessel}}}} {)}$

Increased ${\text{BL}}_{{{\text{vessel}}}}$ means more deformation of the lymphangions’ wall, and hence, an increased ejection fraction. The linear increase in the ejection fraction is directly correlated to the linear increase in volume conveyed. The percentage variations of the volume conveyed and EF of each lymphangion as a result of changing ${\text{BL}}_{{{\text{vessel}}}}$ and relative to the baseline model are shown in Fig. S4.

The P–V loops of the lymphangions for varying ${\text{BL}}_{{{\text{vessel}}}} { }$ were compared (Fig. 8B). As the magnitude of the contraction load increases, the minimum volume at the end of systole decreases due to increased deformation leading to a bigger variation in intralymphatic pressure during a contraction cycle, more prominently for the systolic pressures than the diastolic ones.

8.3 Boundary Load at the Interstitium (${\text{BL}}_{{{\text{interstitium}}}} {)}$

A constant load with a range of magnitudes from 5 to 50 Pa was applied at the boundaries of the interstitium (see Fig. 3B) to recreate the extrinsic pumping, while ${\text{BL}}_{{{\text{vessel}}}}$, representing the intrinsic pumping, stayed the same as in the baseline model. Applying different boundary loads of 5 Pa, 10 Pa, 25 Pa, and 50 Pa at the outer surfaces of the interstitium, increased the volume conveyed by 1.31, 2.54, 6.39, and 12.78%, compared to the baseline model where no external forces (0 Pa) were applied to the interstitium (Fig. S5 Top). As expected, applying higher loads to the interstitium generated more $V_{{{\text{outlet}}}}$ as the direct result of enhanced deformation of the collecting vessel both during the contraction of lymphangions and afterwards. The same reasoning can explain the percentage variation of lymphangions’ EFs relative to the baseline model versus ${\text{BL}}_{{{\text{interstitium}}}}$, presented in Fig. S5 (Bottom).

As the magnitude of the interstitial boundary loads increases, the P–V loops widen due to higher intralymphatic systolic pressure and smaller end-systolic volumes (Fig. 8C). The starting point of contraction (point A in Fig. 5), representing the initial pressure and volume, remains almost the same. Since the boundary load applied at the interstitium is constant and continuous, the compression of the collecting vessel continues beyond the contraction cycle. That is why, in Fig. 8C, the ultimate volume of the lymphangion by the end of the contraction cycle (point G in Fig. 5) did not fall on the starting point of the contraction (point A Fig. 5), creating a seemingly open loop.

8.4 Refractory Period (RP)

The refractory period is the recovery time in between consecutive contractions of a lymphangion. An RP of zero means that the lymphangion’s second contraction starts immediately as the first one finishes.

During diastole, relaxation and dilatation of the lymphangions cause a pressure drop, resulting in intralymphatic pressures lower than the pressure at the inlet. This pressure drop creates a suction effect that causes the first valve to open and drags the fluid from upstream (the fast and slow filling phase, refer to Fig. 5). However, an RP that is too small would cut the duration of the suction effect short with the second contraction cycle starting before the intralymphatic pressures are fully recovered to the initial value (0 Pa for this case). In our model, the impact of reducing RP on the amount of conveyed volume was only prominent for RPs that were smaller than 1 s with the RP of 0.5 s leading to 0.9% and the RP of 0 s to 9.36% less conveyed volume than the case with a RP of 1 s, 1.5 s, and 2 s after a given number of cycles (Fig. S6).

The P–V loops for different RPs are presented in Fig. 8D. It is shown that with a refractory period of 0.5 s, the slow filling phase (points F–G in Fig. 5), where the pressure is to rise back to its initial value, is not completely fulfilled and even less for an RP of 0 s. Hence, contraction cycles with very short refractory periods are not efficient as they propel less fluid forward despite having the same ejection fraction. On the other hand, for a given amount of time, a collecting vessel with high contraction frequency would propel more lymph while having less conveyed volume per contraction. It is a matter of efficiency versus frequency.

8.5 Maximum Permeability in the Secondary Valves ($k_{{{\text{max}}}}$)

As the value of the maximum permeability $k_{{{\text{max}}}}$ increases ($10^{ - a}$in Eq. 1), the valve’s resistance to forward flow decreases allowing the lymph to pass more easily. However, increased permeability also means less resistance to backflow and hence more leakage from the valves. With the balance between the forward flow resistance and the resistance backflow, in our model, the conveyed volume was found to be maximum for the $k_{{{\text{max}}}}$ of $10^{ - a} = 10^{ - 11} \; {\text{m}}^{{2}}$, the value used for the baseline model. However, the behavior is not linear (Fig. S7). Reducing $k_{{{\text{max}}}}$ to $10^{ - 12} \;{\text{m}}^{2}$ and $5 \; \times \; 10^{ - 13} \; {\text{m}}^{2}$ reduced the conveyed volume for 0.89 and 2.05%, respectively, relative to the baseline model. While increasing $k_{{{\text{max}}}}$ to $10^{ - 10} \;{\text{m}}^{2}$ and $5 \times 10^{ - 9} \;{\text{m}}^{2}$ reduced the conveyed volume for 0.45 and 96.64%, respectively, compared to the baseline model. Compared to the highest studied permeability ($k_{{{\text{max}}}}$ of $5 \times 10^{ - 9} \;{\text{m}}^{2}$), the effect of other values on conveyed volume seems almost negligible.

By reducing the $k_{{{\text{max}}}}$ of the secondary valves, the P–V loops stretch out in a vertical direction which corresponds to a bigger variation range for intralymphatic pressures (Fig. 8E). There is also a slight increase in the minimum volume of the lymphangion and hence a slightly smaller deformation which is the result of the increased intralymphatic pressure.

In Fig. 9, the flow rate passing each valve is shown for multiple values of $k_{{{\text{max}}}}$. For the $k_{{{\text{max}}}}$ of $5\; \times \; 10^{ - 9} \;{\text{m}}^{2}$, the permeability is so high that the results were as if there were no valve resistance present, with the backflow being comparable to the forward flow, since the transvalvular pressure was just not high enough to close the valve (see valve function in Fig. 3A).

Fig. 9

Effect of maximum valve permeability on flow rate. The flow rate of lymph passing valves 1–4 over time for different values of maximum valve permeabilities

8.6 Elastic Modulus of the Vessel Wall ($E_{{{\text{vessel}}}}$)

The results presented in Fig. 10A show that the higher $E_{{{\text{vessel}}}}$, the less the volume conveyed. Increasing $E_{{{\text{vessel}}}}$ means that the wall becomes less prone to deformation in response to the boundary load applied to the lymphangion walls. Hence, increasing $E_{{{\text{vessel}}}}$ also leads to a smaller ejection fraction (Fig. 10B).

Fig. 10

Impact of elastic modulus on volume conveyed, EF, and the P–V loops. The variations of the volume conveyed, in percentage, relative to the baseline model, versus the elastic modulus of the vessel wall (A), the variations of the lymphangions’ EF in percentage, from the baseline condition versus the elastic modulus of the vessel wall (B), and the P–V loops of lymphangion 1 (C) are shown for various elastic modulus values of the collecting vessel wall. Results for the volume conveyed, EF, and P–V loops of lymphangion 1 were also shown for various values of the elastic modulus of the interstitium in D–F, respectively

Figure 10C shows the P–V loops of lymphangion 1 for varying $E_{{{\text{vessel}}}}$ presenting that as the $E_{{{\text{vessel}}}}$ increases the P–V loop gets smaller which would translate into less deformation of the vessel wall during the contraction cycle and, as a result, a smaller intralymphatic pressure variation between systolic and diastolic phases. With increasing $E_{{{\text{vessel}}}}$, the reduction of intralymphatic systolic pressure is more significant.

8.7 Elastic Modulus of the Interstitium ($E_{{{\text{interstitium}}}}$)

Increasing $E_{{{\text{interstitium}}}}$ leads to reduced conveyed volume at the outlet (Fig. 10D) by diminishing the deformation of the vessel wall during contraction, which would mean a smaller ejection fraction for the lymphangions, as shown in Fig. 10E. Increasing $E_{{{\text{interstitium}}}}$ also leads to smaller P–V loops meaning a smaller diameter reduction during contraction and a narrower pressure variation of the lymphangion (Fig. 10F).

9 Discussion

In this study, we simulated the transport of lymph in collecting lymphatics consisting of three lymphangions in a 3D model reconstructed from micro-CT scans of the collecting vessels at the top of the foot of a mouse hind limb. 3D morphological data of lymphatic vessels are rare: they are either very small segments of a collecting vessel captured through microscopy or via invasive injections into the vessel itself [42, 43]. We managed to find a commercially available contrast agent, Ultravist, that is commonly used for cardiovascular imaging and administer it via superficial injection into the interstitium.

In Fig. S8a and S8b, two graphs by Benoit et al. [44] are presented that were obtained from in vivo pressure and diameter measurements of collecting lymphatics in male Sprague–Dawley rats using in vivo microscopic techniques, while Fig. S8c and S8d presents results from our Baseline model. In both sets of graphs, significant instances throughout a lymphangion’s contraction cycle are highlighted by correlating the variations seen in pressure and diameter of the lymphangion (points A–F). Comparing these results shows that the general trends are similar for both figures from Benoit et al. [44] and ours, as represented by the following instances: the start of systole (A), the peak intralymphatic pressure (B), the minimum diameter and the end of the systole (C), the start of rapid filling of the lymphangion (D), the end of the rapid filling phase and the beginning of slow filling (E), and the end of the slow filling phase (F).

There are also clear differences. For instance, in Fig. S8a, both base and peak pressures are higher than our model. It can also be seen in Fig. S8b that the diameter of the collecting vessel reduces for 12 µm, at its peak deformation; however, at the same instance (point C), the deformation of the lymphangion in our model was equal to 8.34 µm (Fig. S8d). Higher deformation in our 3D model could be achieved by increasing the lymphangion’s boundary load which would also result in higher intralymphatic pressure values. A pressurized vessel (non-zero inlet and outlet) and an adverse pressure at the outlet can also contribute to much higher intralymphatic pressures.

Based on the parametric study, the most impactful parameter overall was found to be the maximum permeability of the secondary valves for the case of $k_{{{\text{max}}}} \; = \;5 \; \times \; 10^{ - 9} \;{\text{m}}^{2}$, as it could lead to faulty valves for which the volume conveyed dropped drastically. The permeability of the valve in the open state is an important factor, since it provides the necessary pressure drop that leads to the closure of the valve. However, for the extreme case of $k_{{{\text{max}}}} \; = \;5\; \times \;10^{ - 9} \;{\text{m}}^{2}$, the transvalvular pressure varied only between − 0.02 and 0.01 Pa throughout the contraction cycle, which is smaller than $P_{{{\text{th}}}}$. In other words, the pressure downstream of the valve was just not high enough, relative to the pressure upstream of the valve, to close it. Amongst the different values of $k_{{{\text{max}}}}$ in the parametric study, our results also showed that the maximum $V_{{{\text{outlet}}}}$ was achieved for $k_{{{\text{max}}}} = 10^{ - 11} \; {\text{m}}^{2}$ (Fig. S7). For this value of $k_{{{\text{max}}}}$, allowing higher permeability (high enough but not too high) for the forward flow outweighs the backflow (Fig. 9), which is in line with the findings of Jamalian et al. [18]. In their paper, a parametric study was performed on the improved lumped parameter model of Bertram et al. [16]. Their work emphasized that besides the importance of preventing backflow for the overall function of the lymphatic system, the lymph flow is highly sensitive to any source of resistance to forward flow. In our model, for a valve such as valve 1, having the maximum permeability of $k_{{{\text{max}}}} = 10^{ - 11} \; {\text{m}}^{2}$, is the equivalent of having the forward resistance (valve resistance in the open state) of $6.2\; \times \;10^{10}$ $\frac{{{\text{Pa}}\;{\text{ s}}}}{{{\text{m}}^{{3}} }}$. Bertram et al. [17] reported an open-valve resistance of approximately $0.6 \times$ ($6 \; \times \; 10^{10}$ $\frac{{{\text{Pa}}\;{\text{ s}}}}{{{\text{m}}^{{3}} }}$) for segments of rat mesenteric lymphatic vessels including one valve which makes the better performed value of $k_{{{\text{max}}}}$ a good approximation.

In the parametric study, the properties of all valves were changed simultaneously. It may, however, also be relevant to assess the impact of one or multiple dysfunctional valves ($k_{{{\text{max}}}} \; = \;5 \; \times \;10^{ - 9} \;{\text{m}}^{2}$). The results for different scenarios of faulty valves and the impact on $V_{{{\text{outlet}}}}$ are presented in Fig. S9 which shows that having a very high permeability (in the open state) for either valve 1 or valve 4 has more impact on $V_{{{\text{outlet}}}}$ than a case where both valves 2 and 3 were not able to close properly. When both intermediate valves were leaking, the conveyed volume would be reduced to 1.36 nl, compared to 1.43 nl with one intermediate leaking valve. Nonetheless, this value was still more than the conveyed volume when either valve 1 or 4 was dysfunctional (i.e., 0.91 nl).

The parameters influencing the maximum deformation of the lymphangion by the end of systole and, consequently, the volume conveyed, from the most impactful to the least, were the elastic modulus of the interstitium, the elastic modulus of the vessel, the magnitude of the lymphangion’s boundary load, and the magnitude of the boundary load at the interstitium.

The smaller the elastic modulus of either the interstitium or the vessel wall, the easier the vessel wall deforms, creating a bigger compression of the lymphangions. Interestingly, decreasing the elastic modulus of the interstitium had a sharper effect on the volume conveyed than decreasing the elastic modulus of the vessel wall (Fig. 10), which is due to the fact that, along with the thin collecting vessel wall deforming, the poroelastic interstitium is also stretched.

The higher the lymphangions’ contraction load, the higher the compression, ejection fraction, and volume conveyed (Fig. S4). Within the range of values studied, the volume conveyed changed linearly with the boundary load. However, we expect that by applying higher boundary loads, a plateau would be reached for the volume conveyed, since above a certain threshold, the lymphangion would collapse. On the other hand, once the lymphangion collapses, it would hinder further lymph propulsion. This was shown by Jamalian et al. [18] describing that increased contraction load initially increased the pumping, but then it impeded downstream lymphangions to incoming forward flow.

The final parameter directly contributing to wall deformation is the load at the interstitium. This parameter is in a sense simulating the extrinsic pumping or the force that is applied to the extremities in the form of a massage, as part of lymphedema treatment. The constant and continuous interstitial load kept compressing the collecting vessels and propelling the lymph even after their contraction cycle had ended, albeit with much lower velocities (25 µm/s). Figure S5 (top) shows improved lymph propulsion as a result of this interstitial boundary load. Other configurations of this load could also be studied (for instance peristaltic interstitial loads in the direction of lymph flow) to find the most effective configuration in terms of the conveyed volume.

Under physiological conditions, the collecting vessels propel lymph against adverse pressure gradients. In our baseline model, we decided to create a non-pressurized vessel with zero pressure gradient ($\Delta P = P_{{{\text{outlet}}}} - P_{{{\text{inlet}}}}$) as a neutral state to focus on the impact of intrinsic pumping on lymph propulsion. However, our parametric study with $P_{{{\text{outlet}}}}$ showed that the higher the adverse pressure at the outlet, the less volume was conveyed. The last valve in immediate contact with the outlet, valve 4, remained closed for the majority of the computation time except during the systole of the last lymphangion. We also noticed fluid accumulation and swelling in the lymphangions due to the long closure of valve 4. With increased intralymphatic pressures, diminished deformation and compression of lymphangions were observed. Correspondingly, Bertram and colleagues’ lumped parameter model of a chain of lymphangions [16] also showed a generally reduced mean flow rate, distension of the lymphangions, and smaller lymphangion diameter variations during the contraction cycle. In Fig. S3, we showed how $P_{{{\text{outlet}}}}$ increased the intralymphatic pressures. Similarly, Bertram et al. [16] reported that the intralymphatic pressure of the lymphangion upstream of the final valve was capped at values barely exceeding $P_{{{\text{outlet}}}}$, while the intralymphatic pressures of intermediate lymphangions along the chain reached peak values much higher than that of $P_{{{\text{outlet}}}}$. It is also important to note that Bertram and colleagues reported these results for a pressure gradient ($P_{{{\text{outlet}}}} - P_{{{\text{inlet}}}}$) of 10 Pa, close to the adverse pressure ranges in this study. On the other hand, they stated that for adverse pressures in the order of 60 Pa and higher, the nature of the pumping changed with flow patterns becoming more sinusoidal and regular, as well as showing more backflow. Davis et al. [45] studied the responses of an isolated single lymphangion to elevated $P_{{{\text{outlet}}}}$ experimentally and documented the changes in the lymphangion’s P–V loops. In agreement with our results (Fig. 8A), Davis and coworkers reported that as $P_{{{\text{outlet}}}}$ increased, the P–V loops became more vertical and rectangular in shape with a leftward shift in the end-systolic P–V relationship.

A pressurized collecting vessel could have been replicated by applying pressure at the inlet and the outlet of the collecting vessel. However, our research question was mainly focused on the volume of lymph that was transported and specifically the impact of varying different parameters. On the other hand, the pressurized vessel would introduce extra challenges for the model’s convergence, which would have been more critical to incorporate if we were looking at the stresses in both the vessel wall and the interstitium. We also believe that the conclusions that we made based on our parametric study would ultimately still be very similar.

Our parametric study on the refractory period showed that in order for the lymphangions to have an efficient pumping contraction, enough time is needed in between a lymphangion’s contractions for it to recover (Fig. S6). In other words, enough time is needed for the lymphangion to fill up again, both fast and slow filling as the suction effect manifest itself during diastole (Fig. 5). For our contraction cycle, the optimum RP was 1 s. An RP of less than 1 s did not provide enough time for the lymphangion to recover from the previous contraction; however, increasing the RP beyond 1 s did not lead to more conveyed volume either. In line with our findings, Bertram et al. [20] also showed in their improved lumped parameter model that when the RP was long enough to allow diastolic lymphangion filling to complete before the start of the next round of contractions, the volume pumped per cycle reached a maximum. On the other hand, they stated that eliminating the refractory period (RP = 0 s) would be highly inefficient. Jamalian et al. [18] also suggested that contracting less often than the maximum rate (RP = 0 s) would be a more efficient use of the metabolic energy, as the short RP basically cuts the suction effect short. Jamalian et al. [46] emphasized the importance of the suction effect for lymph uptake of the initial lymphatics at subatmospheric pressures.

To assess the importance of considering the 3D geometry and the morphology of a collecting vessel, we also created a 2D axisymmetric model where the volume of each lymphangion corresponded to the ones of the 3D model of the collecting vessel. The comparison of the adapted 2D axisymmetric model and the 3D model is shown in Fig. S10. The only result correctly predicted by the 2D axisymmetric model with adjusted volumes was found to be the amount of propelled lymph, the conveyed volume. The deformation of the vessel wall of the 2D model also provided a good approximation, except for lymphangion 2 where the collecting vessel has a clear tightening. The other results regarding the characteristics of lymph flow such as its speed and pressure, as well as the permeability of the secondary valves are quite different, which showcases the added value of the 3D model despite it being more computationally expensive. The change in morphological parameters along the 3D reconstructed collecting vessel leads to more (local) variations in terms of intralymphatic pressures and lymph velocities. As Bertram et al. [16] showed, downstream lymphangions with smaller diameters contribute to the flow resistance against the upstream lymphangion to an increasing extent. With higher transvalvular pressure drops, the permeability of the valves would also decrease more and contribute to the extra resistance.

10 Limitations

To avoid computational instabilities, we had to make some compromises, for instance in the level of deformation, which leads to smaller intralymphatic pressure gradients and ejection fractions relative to the values reported in the literature. We also modelled the collecting vessel wall as linear elastic which may not completely represent its true physiological behavior. The secondary valves were simplified as a thin porous membrane; including their true morphology in the model would also impact the flow patterns and boost the rationale for the 3D modelling. At the same time, it would also further drastically complicate the modelling, including challenges both in terms of achieving stable models given the contact of the secondary valve leaflets during closure, and significant increases in computational costs.

By testing different RPs over multiple cycles, we realized that for a small chain of lymphangions, such as ours, the periodicity was easily obtained, as long as the next contraction cycle did not interrupt the previous one. However, in a much longer chain of lymphangions, it would not have been possible to avoid running many contraction cycles before one could have achieved periodicity.

11 Conclusion

In conclusion, the present study evaluated a 3D reconstructed model of a multi-lymphangion collecting vessel in a poroelastic interstitium. Despite a number of simplifications, our results and behavioral patterns were in harmony with the literature, encouraging using these models as a stepping stone to gradually build up on their complexities. While previous lumped parameter models could provide valuable information on the average values or the total volume of lymph transported, a 3D model with a realistic morphology offers the advantage of providing more information on a local scale and its effect on lymph propulsion. Furthermore, the 3D model represents a building block for future work on an integrated multi-scale model of the lymphatic network that we will use to support the testing and optimization of the in silico design of an implantable lymphatic pump.