The effects of valve leaflet mechanics on lymphatic pumping assessed using numerical simulations.

The lymphatic system contains intraluminal leaflet valves that function to bias lymph flow back towards the heart. These valves are present in the collecting lymphatic vessels, which generally have lymphatic muscle cells and can spontaneously pump fluid. Recent studies have shown that the valves are open at rest, can allow some backflow, and are a source of nitric oxide (NO). To investigate how these valves function as a mechanical valve and source of vasoactive species to optimize throughput, we developed a mathematical model that explicitly includes Ca2+ -modulated contractions, NO production and valve structures. The 2D lattice Boltzmann model includes an initial lymphatic vessel and a collecting lymphangion embedded in a porous tissue. The lymphangion segment has mechanically-active vessel walls and is flanked by deformable valves. Vessel wall motion is passively affected by fluid pressure, while active contractions are driven by intracellular Ca2+ fluxes. The model reproduces NO and Ca2+ dynamics, valve motion and fluid drainage from tissue. We find that valve structural properties have dramatic effects on performance, and that valves with a stiffer base and flexible tips produce more stable cycling. In agreement with experimental observations, the valves are a major source of NO. Once initiated, the contractions are spontaneous and self-sustained, and the system exhibits interesting non-linear dynamics. For example, increased fluid pressure in the tissue or decreased lymph pressure at the outlet of the system produces high shear stress and high levels of NO, which inhibits contractions. On the other hand, a high outlet pressure opposes the flow, increasing the luminal pressure and the radius of the vessel, which results in strong contractions in response to mechanical stretch of the wall. We also find that the location of contraction initiation is affected by the extent of backflow through the valves.

Introduction

Lymphatic function is critical for fluid homeostasis, and local failure leads to edema with significant morbidity and local immune dysfunction[1,2]. The collecting lymphatic vessels are responsible for pumping fluid from the tissue and delivering it to lymph nodes in order to prevent edema and maintain immune function. Although potentially beneficial for a wide variety of patients, there are currently no pharmacological modulators of lymphatic function[3-5]. The prerequisite for creating such drugs is a fundamental understanding of the mechanisms that influence lymph clearance.

The lymphatic system consists of lymphatic capillaries - termed initial lymphatics - that absorb interstitial fluid from the tissue and have unique overlapping endothelial cells[6], which function as one -way valves to allow fluid to enter the vessels and prevent backflow into the tissue[7,8]. Initial lymphatics pass the fluid to collecting lymphatics, which contain distinct compartments defined by luminal one-way valves, known as lymphangions, that contract spontaneously to drive flow[9-12]. Thus, a collecting lymphatic vessel behaves as a group of distributed pumps aligned in series and regulated by local physiological conditions[13]. Depending on the fluid pressure conditions in the tissue and lymphatic system, collecting lymphatic vessels can adjust their behavior to maintain fluid homeostasis. Lymphatic vessels minimize contractions if existing fluid pressure gradients can drive flow, but actively pump to drive fluid otherwise[10,14,15].

Our previous work has elucidated how these behaviors are coordinated throughout the lymphatic system by mechano-sensitive feedback[15,16]. In summary, our model shows that and NO concentrations establish complementary and oscillatory feedback loops that are self-regulating, maintaining normal lymphatic function, in agreement with experimental observations[12,17-19]. However, in our previous work, the intraluminal valves were modeled by mathematically inserting a semi-permeable wall in the vessel when the flow reversed direction. Although this approach was able to reproduce the correct behaviors seen in vivo, it simplified valve performance and neglected valve leaflet structure and mechanics.

The intraluminal valves that bias the flow in the direction away from the peripheral tissues are extremely important for lymphatic function, but little is known about their mechanical properties or how they influence lymph clearance or contraction efficiency. Previous studies have determined that lymphatic valves are biased to the open position, especially when the vessel is partially distended by a trans-wall pressure gradient[20]. Experimental studies and mathematical models show that this property can increase flow efficiency by reducing flow resistance when the vessel is not actively contracting, even though it results in significant backflow as the downstream lymphangion contracts[21-24]. The valve leaflets are also important as sources of biomolecules. Experiments have shown that the leaflets are a significant source of nitric oxide, presumably produced by dynamic shear stress on the leaflet structures themselves[25].

Here, we extend our previous model of collecting lymphatic vessels by introducing realistic intraluminal valves. Furthermore, we embed the pumping vessel within a poro-elastic tissue that has fixed pressure boundaries in order to simulate fluid transport and edema at the tissue level.

Model Description

Figure 1 shows the model domain and dimensions. Fluid can enter the tissue from anywhere on the boundary except at the vessel outlet. Fluid percolates through the tissue to enter the initial lymphatic capillary at left (represented by the dashed line). The solid regions in this segment are impermeable, while the open regions are semipermeable, simulating the primary valves in the lymphatic capillaries. The collecting lymphatic vessel is downstream from the initial lymphatic, and fluid passes through it to exit at right. This vessel has two intraluminal valves that bias the flow toward the exit, and the walls can move to actively pump fluid. We use the lattice Boltzmann method (LBM)[26,27] to calculate the fluid flow, shear stresses and pressures in the tissue and lymphatic vessel. The endothelium on the inner wall of the vessel can generate nitric oxide (NO)[28-30] in response to increased shear stress, and contractions of the collecting lymphatic are determined by the concentration of in the lymphatic muscle cells. Briefly, can be depleted over time due to recharging of the cytoplasm ion concentrations, and can increase due to mechanical or chemical triggers. The self-regulating contraction dynamics that effect lymph drainage are governed by mutual mechanical feedback of the NO and systems[15,16]. The details of the model are given in the Methods section.

Figure 1

Diagram of the lymphatic vessel. R0 = 100 μm is the relaxed radius of the vessel.

Methods

The lattice Boltzmann method

We use the lattice Boltzmann method (LBM)[26,27] to simulate fluid flow through the tissue and lymphatic vessel. The curved boundary condition[31] is used to keep track of the location of the vessel wall and the valves with respect to the lattice grid. The single relaxation time approximation to the Boltzmann equation can be discretized in both space and time aswhere is the relaxation time. is the distribution function at time t and location x. e is the discrete velocity. is an equilibrium distribution function that depends on fluid density and velocity. Through a Chapman-Enskog expansion, the Navier-Stokes equations can be recovered and the kinematic viscosity depends on as

For the D2Q9 model, one form of the equilibrium distribution is[27]where, , , and . and u are the fluid density and velocity defined by

To smooth out the hydrodynamic forces on the solid surfaces (vessel wall and valves), we use the stress -integration method. The stress tensor in the lattice Boltzmann method can be calculated as[32]where is the Kronecker delta function and . The hydrodynamic force acting on area S can be calculated by integrating the stress tensor and momentum flux over S[32]:where n is the unit normal vector on S oriented towards the fluid. V is the velocity of S. In order to calculate the hydrodynamic force exerted on S, the distribution function f on S in Eq. (5) must be known, and f can generally be calculated by extrapolating from the nearby fluid nodes. Here, to simplify the calculations, we use the nearest fluid node values to estimate the distribution function on S. Under conditions of low Reynolds number, the resulting errors are negligible. We have previously confirmed the robustness of this model in simulations of cells and blood flow[33-36].

The lymphatic vessel model

The schematic diagram of a lymphatic vessel with single lymphangion defined by two valves and embedded in tissue is shown in Fig. 1. The diameter of the vessel is 100 μm. Fluid from the tissue can enter the porous vessel with fixed geometry at left, which represents the initial lymphatic capillaries[7,37,38]. The solid regions in this segment are impermeable, while the open regions are semipermeable, simulating the primary valves in the lymphatic capillaries. Lymphatic capillary valves allow fluid entry but resist leakage back into the tissue, so we institute a partial bounce back boundary condition in the “open” regions of this segment[39]. We use a constant bounce back ratio , which means that 15 percent of the fluid is allowed to leak back to the tissue when the pressure conditions favor flow in that direction. The area ratio between the semipermeable and impermeable parts is 1:1. The right end of the lymphatic vessel represents the outlet, which would be connected with additional downstream lymphangions or a lymph node. In this region, we introduce a fixed segment 234 μm long to reduce the effect of the right side boundary and provide structural support.

Although many vasoactive substances can alter vessel contractility[40-42], NO appears to play a major role in lymphatic muscle cells. NO is rapidly produced by endothelium in response to changes in shear stress and is quickly degraded[43]. For simplicity, our analysis focuses on the effects of NO on vessel contractions, but could be generalized to other flow-mediated relaxation factors. The endothelium on the inner wall of the vessel can generate nitric oxide (NO)[28-30], which is allowed to convect and diffuse in both the fluid and tissue spaces[15]:where is the concentration of NO. is the diffusion coefficient of NO. can be approximated using a second order finite difference scheme:where is the discretized distance. can also be approximated using an upwind scheme

In Eq. (7), λ is the chemical reaction rate constant. determines the chemical reaction time scale. Δt is the lattice time scale. As λ increases, the chemical reaction rate increases. and are the decay rate and production coefficient of NO. Here we consider that the production of NO is proportional to the stress component [16,44,45]. v indicates fluid velocity along the tangential direction of the membrane surface and x is the normal direction vector. Within the lymphangion, the valve structures can also generate NO. In fact, as shown in Fig. 2, we can indirectly calculate using the hydrodynamic force dF we calculate in Eq. (6). Thus .

Figure 2

Hydrodynamic force d acting on surface element ds.

Ca2+ enters and leaves the cytoplasm of the lymphatic muscle cells and can also pass through junctions to neighboring cells[46-49]; therefore, we use the reaction diffusion equation, with diffusion restricted to the 1D space of the vessel wall[15]:

Here CCa is the concentration of Ca2+. DCa is the effective diffusion coefficient of spreading from one cell to the neighboring cells along the vessel wall. Note that the speed of the signal in the confined space of the vessel wall is faster than 3D diffusivity, and is related to the speed of an action potential[50]. can also be approximated through Eq. (8), but here, δ should be changed to a cell length of 1.51 on the lattice. is the decay rate of . The total decay rate of is multiplied by based on the fact that NO acts on the lymphatic muscle cells through myosin light chain phosphatase to reduce force generation[51-53]. is the production rate of . determines the strength of the nonlinear elasticity of the membrane[16]. The 11th power term was arbitrarily chosen to provide a rapidly increasing force that prevents the solid elements from colliding, which would cause a logic error. Other similar functions could be substituted without significant consequence. R and RCa are the local radius of the vessel and the reference radius of the nonlinear term respectively. R is the limiting minimum radius of the vessel. is an asymmetric Kronecker δ–function. If CCa increases from below to exceed the threshold Cth, this function is set to 1; otherwise, it is set to 0. This is based on the fact that when CCa reaches the threshold Cth, it triggers action potentials mediated by voltage gated and calcium-induced calcium channels[16,54]. defines a steep increase of CCa while .

There are five forces that act on the wall of the vessel:

The hydrodynamic force F, which is calculated by stress integration using the LBM. In order to calculate F, the walls of the vessel are discretized into N segments, and the segments can only move along the y direction.

The lymphatic muscle force depends on the concentrations of and NO according to[15]:where is the coefficient determining the strength of action.

Elastic force from the tissue:

In order to limit the contraction amplitude and reproduce tissue mechanical resistance, if , we simply increase F by multiplying an 11th power item. Eq. (11) changes to:

Bending force:

Here we consider that the left and right sides of the vessel wall are fixed; in between, the wall can move along the y-direction. and are the y-direction positions of the neighboring points of the vessel.

Considering the vessel is viscoelastic, a viscous resistance force is introduced as:

The minus sign means that always acts in the direction opposing wall velocity .

As shown in Fig. 3, the bileaflet valve (illustrated in Fig. 3(b)) has structure similar to that of a real lymphatic valve (Fig. 3(a))[55,56]. The valve is treated as two opposing deformable parabolic - shaped membranes. Compared with previous models of valves[57-60], our 2D valve model has been simplified to provide computational efficiency for the multi-lymphangion simulations while, at the same time, recapitulating the relevant physics, biology, physiology and biochemistry of the lymphatic vessel.

Figure 3

Simulating intraluminal lymphatic valves. (a) Intravital image of a lymphatic valve in the popliteal collecting lymphatic of a mouse. (b) Parabolic valve shape assumed for the model.

The membrane also has elastic force , bending force and viscous force shown in Eqs (11–13). The rest position of lymphatic valves has been shown to be biased in the open position[22,23]. Thus, we set our valve position to a parabolic line (solid line in Fig. 3(b)), the equation of which iswhere indicates the rest or limit position of the valve. is the rest position of the vessel, is the anchor point of the valve. ‘−’ and ‘+’ indicate upper and lower valve leaflets. can adjust the rest position to be biased to stay open or close. For the limit position, is maximum. In the open configuration, the and of the valve are designatedwhere i indicates a discrete segment of the valve. When the two membranes of a valve come into close proximity, we also increase by multiplying an 11th power term:where is the center line of the vessel (dash line in Fig. 3(b)). When the segment arrives at the limit position y shown as in Fig. 3(b), we also increase as:

The valve mechanical properties have evolved to ensure proper valve operation. The leaflets must be flexible, so that they can form an effective seal, but rigid enough to resist the hydrodynamic forces as the valve closes during backflow. Therefore, we assume the valve mechanical properties vary along the length. Although this assumption has not yet been verified for these small lymphatic valves, mathematical models and experiments have shown that large valves in veins have such spatially- variable stiffness[60] (see Fig. 3(b)). To satisfy both these requirements, we specify that the leaflets have variable rigidity by dividing them into segments: near the attachment region at the wall, the leaflets are stiff, but they become more flexible at the tip. Specifically, the varying bending coefficient is set according to:where indicates segment number and is the total number of segments. A membrane of the valve is discretized from left to right into segments. is the maximum of . is an approximation of the minimum of , the coefficient adjusts the soft range of the valve. Figure 4 shows how changes over the length of the valve leaflet.

Figure 4

Bending modulus of the valve decreases from the anchor point to the free end.

To ensure that the total inner force is zero, Eq. (12) changes to:

There are special considerations as the valve leaflets close. If we specify that the valve can close completely, then there will be no unoccupied lattice points on the center line of the vessel, and the two membranes could approach without being separated by any fluid nodes. We avoid the consequent computational difficulties and unrealistic fluid dynamics by introducing a lubrication force that establishes the hydrodynamic force in this region (see “Lubrication force” below). When the membranes of the valve approach the center line of the vessel, we change the rest position of the valve to a straight line and the bending force is changed to:

For motion of the valves and vessel wall, we assume the structures have the same viscosity. The vessel and the valves move according Newton’s laws of motion. At each Newtonian dynamics time step, the center of mass of each segment is updated by a so-called half-step ‘leap-frog’ scheme[61]:where V is the velocity of the center of mass of each segment, δt is the Newtonian time step, here chosen to be lattice time step, and a is the acceleration of each segment. When the segments move, new fluid nodes are revealed; we extrapolate the status of each new fluid node based on the known quantities of those fluid nodes located on the same side of the moving boundary.

Lubrication force

When two membranes approach each other, they can displace all the fluid, and the hydrodynamic force cannot be calculated using LBM. In this case, we use lubrication theory to calculate the hydrodynamic force[62]. As shown in Fig. 5, we only consider the normal lubrication force (which acts in the y direction, because valve segments move only in the y direction) to determine the force on the valves.

Figure 5

Schematic of lubrication force between to planes.

According to lubrication theory[63], the Stokes equation can be represented bywhere p is the pressure, is the dynamic viscosity, and is the fluid velocity along the x-direction. The boundary conditions are

Integrating Eq. (22) of ywhere:

From Eq. (23)

The flux is:

The fluid volume between two planes in Fig. 5 is

During virtual variation of δt, the variation of h0 is:then the variation of the volume is:

The flux caused by the relative movement in the normal direction is:

The flux is caused by flow from upstream. Because Eq. (22) is linear and the relationship between Q and u is also linear, the general flux can be written as:

Because , , the Rayleigh’s dissipation function is[64]:then the normal lubrication force is:

Let (upper and lower leaflets of the valve are symmetric):

For two planes constructed with n connected segments, each having length l, if each segment moves with the same velocity, the lubrication force exerted on an individual segment should be:

If some segments move up () and others move down (), for example, segments move up and segments move down, the lubrication force exerted on segment can be estimated as:

The simulation domain is 433 × 66 lattice points. The calculations are performed using an NVIDIA Quadro M5000 graphics card with 2048 Cuda cores. The number of parallel threads depends on the calculations being performed. For example, for the calculations of fluid node density with LBM, we send 28578 threads to Cuda, which distributes them amongst 2048 cores. Host memory is allocated to the GPU to avoid excessive data transfer between the CPU and GPU. Codes were written in C++.

Simulation parameters

The system time scale depends on the fluid kinematic viscosity and diameter of the vessel as: , where , and , D′ are the actual fluid kinematic viscosity and diameter, respectively. T is the duration of one lattice step. The space scale depends on , where L is the length of one lattice unit. Here we choose , , and . Thus, real time and space scale are and . The chemical reaction speed depends on λT. In our simulation, we choose the simulation parameters as given in Table 1.

Table 1

Chemical parameters of NO and Ca2+; Mechanical parameters of the fluid, vessel wall and valves.

Parameter	Definition	Value	Units	Source
FLUID
L	Lattice space unit	4	μm
T	Lattice time step	1.33 × 10−6	s
ρ_in	Fluid density at inlet	1	g/cm3	As water
ν′	Fluid viscosity	0.01	cm2/s	As water
τ	Relaxation time	0.75	T
Chemical properties of NO and Ca 2+
D NO	NO diffusivity	1.2 × 10−4	cm2/s	3.0×10^−5(cm^2/s) [75]
KNO−	NO degradation rate constant	3.7594	s^−1	Estimated
KNO+	NO production rate constant	400	Dimensionless	Estimated
D_Ca	Ca^2+ signal propagation rate	6.5×10^−6	cm^2/s	Estimated
KCa−	Ca^2+ degradation rate constant	37.6	s^−1	Estimated
KCa+	Ca^2+ production rate constant	1.2	s^−1	Estimated
Kδ+	Ca^2+ production rate constant	15038	s^−1	Estimated
C_th	Ca^2+ threshold	0.015	Dimensionless	Estimated
R_Ca	Threshold radius for Ca^2+ channel sensitization	R_0		Estimated
K_Ca,NO	Rate constant for NO inhibition of Ca^2+	5.3	Dimensionless	Estimated
λ	Chemical reaction rate constant	0.03	Dimensionless	Estimated
VESSEL
K_M	Force constant for Ca^2+	7.6×10^−5	dynes	Estimated
K_E	Elastic modulus the vessel	4.52	dynes/cm^2	Estimated
K_B	Bending modulus the vessel	9045	dynes/cm^2	10^6(dynes/cm^2) [66]
K_r	Viscosity coefficient of vessel	4.8×10^−9	dynes⋅s/cm	Estimated
K_NO	NO inhibition of force	1	Dimensionless	Estimated
R_l	Limit radius	0.003	cm	40% contraction[20]
R_0	Rest radius of the vessel	0.005	cm
VALVE
A	How soft the valve is	5	Dimensionless	Estimated
B	How much the valve biased to open	1500	cm^−1	Estimated
kEv	Elastic modulus of valves	9.0 × 10−4	dynes/cm^2	Estimated
k0v	Bending modulus of the base of valves	18090	dynes/cm^2	10^6(dynes/cm^2) [66]
kRv	Bending modulus of the tip of valves	0.0091	dynes/cm^2	0.1k0v [60]
VESSEL & VALVE
Δ		2 × 10−4	cm	Estimated

Parameter values were estimated in order to yield physiologically reasonable results when independently measured values were unavailable.

To simulate the different masses of the vessel wall and valve leaflets, we specify different densities for these structures on the lattice. Specifically, the density of the vessel wall is assumed to be 80 times that of the valve leaflet. This incorporates both the increased cellularity of the wall vs. leaflet, but also the structural anchoring of the wall to the surrounding tissue. Initially, all the fluid node densities are set to 1, the velocities are set to 0, and the NO concentration is set to 0. The initial concentration within the vessel wall is set to 0.0149, which is close to, but less than, . The density at the boundary (red area in Fig. 1) is kept constant as by imposing an equilibrium distribution with a velocity of 0, and the density at the outlet is kept constant as through a constant pressure boundary condition[65]. NO concentration on the boundary is held constant at 0. NO can diffuse through the fluid, valve structures and vessel wall, but distributes only within the vessel wall. At the left and right sides of the domain, the calcium is treated similar to a bounce-back boundary condition. In our simulation, we keep . Since the fluid is treated as water, the pressure unit is , thus the inlet pressure is held at . The upper and lower vessel walls are discretized into 200 segments; the length of each is 6.03 μm. Each valve is discretized into 32 segments; the length of each is 6.25 μm. The Young’s modulus unit is . For the bending rigidity, the shear modulus is Kdl, where dl is length of each segment of vessel or valve. Thus, the shear modulus of the vessel or valve can be calculated as and , giving 13658 and 28266 dynes/cm2 respectively. Here 1.51 and 1.56 are the dimensionless lengths of each segment of vessel and valve on the lattice. These values for the moduli are less than those estimated for valves in large veins[66], which are on the order of 106 dynes/cm2, but are reasonable, given the much smaller size of our valves. Table 1 lists the model parameter values and their origin.

Results

The effect of valve mechanical parameters A and B on pumping

We first investigated how valve mechanical properties affect lymph clearance and pumping behavior under equal pressure at inlet and outlet. To do this, we set the bounce back condition of the porous part of the initial lymphatic region (Fig. 1, dashed line) to allow more leakage. Setting the bounceback ratio to zero means that fluid can leave the vessel as easily as it enters. This puts more burden on the intralumenal valves to control backflow. We then varied the mechanical properties of the valve leaflets and measured performance. We first varied the parameter A, which describes the spatial dependence on bending modulus along the length of the leaflet. When A is zero, the tip section of the leaflet is as rigid as the base; increasing A makes the tip more flexible than the base (Eq. 17 and Fig. 4). When the rest position condition is set at , increasing A results in a continuous increase in output (Fig. 6(a)). When A is small, the valve is too rigid, and does not close easily, resulting in excessive backflow that limits the flux through the vessel. As A increases, the free end of the valve becomes more flexible, and the valve can deform to resist the backflow. In Fig. 6(b), we keep , and increase B from 0.35L−1 (valve closed at rest) to 0.7L−1 (valve open at rest, see inset figure). Larger values of B mean that the valve is more open at rest. Increasing B causes the cycle-averaged flux to increase before 0.5. This implies that, before 0.5, valves biased to stay open may allow some leakage, hurting efficiency, but this is more than offset by the decrease in overall flow resistance presented by the open valves. Above 0.6, the leakage significantly decreases the efficiency, resulting in decreased flux. This is consistent with previous models and experimental observations, which show that lymphatic valves are biased toward the open position[20,22-24]. In the remaining simulations, we set A = 5, B = 0.6L−1.

Figure 6

Pumping flux is affected by the mechanical properties and rest state of the valves. (a) Output flux is affected by the rigidity parameter A (defined in Fig. 4, increasing A means increased flexibility of the leaflet tip). Flux increases as the tip of the leaflet becomes more flexible, relative to the base. (b) Output flux is increased when the valve rest position is more open (Keeping , increasing B of Eq. (14)). The inset figure shows how B affects the rest position of the valve.

Dynamics of Ca2+ concentration

Figure 7 shows how Ca2+ concentration, vessel diameter and fluid flux change with time when the vessel undergoes sustained spontaneous contractions with ; in this case, the pressure difference between the tissue boundary and lymphangion outlet is . The changes in concentration agree with previous reports using one-dimensional models[16]. In Eq. (9), the equilibrium of can be estimated as , thus we choose and the equilibrium of is less than 1. Meanwhile is large enough to ensure that the equilibrium (when ) is larger than the threshold . Initially, , and  increases from the initial value to the equilibrium value. When becomes larger than (dotted line in Fig. 7(a)), in equation (9), will be 1, and increases steeply. The concentration reaches saturation and stops increasing; consequently, becomes 0. Meanwhile, the vessel begins to contract and drives fluid downstream. This increases the shear stress on the vessel downstream. In response, the endothelial cells on the inner vessel wall and on the both surfaces of the valve leaflet generate NO that appears in the nearby fluid nodes. The NO concentration at the inner fluid node near the vessel and the valve downstream increases (see Fig. 8), and can reach saturation (value = 1). Because of the second term of Eq. (9), large values of cause to decrease and the vessel begins to relax and pull in fluid, generating fluid flow upstream. As a result, the upstream fluid node near the inner surface of the vessel and both surfaces of the valve leaflet increase their generation of NO. When the level recedes below the threshold , the NO concentration also decreases to the baseline. After reaches the minimum, decreased NO allows to increase again, starting another cycle. Because of the rapid diffusion of between neighboring cells, the movements of adjacent vessel segments are coordinated.

Figure 7

The concentration, diameter of the vessel and flux during periodic contractions, . (a) concentration, (b) diameter, and (c) flux. (d and e) are the enlarged curve of concentration in (a). The horizontal dotted line indicates the threshold of concentration.

Figure 8

(a) NO concentration, (b) pressure and velocity field dynamics during a contraction cycle, . The minimum is blue and maximum is red. The arrows lengths are proportional to the local fluid velocity. See also Supplemental Videos S1 and S2.

Vessel contractions and fluid flow

Figure 7(b) shows how the diameter at the midpoint of the lymphangion changes with time. When , fluid flow is driven by the hydrodynamic pressure; in this condition, because of the Venturi effect, the radius of the vessel can be less than R0. In Fig. 7(b), from to , after the previous cycle of relaxation, NO reaches the baseline, and begins to increase (dashed curve in Fig. 7(a)). When , just passes through the threshold, the concentration in the wall suddenly increases and a contraction begins. When , the radius of the vessel reaches its minimum; the duration of the contraction is . After the contraction, the relaxation stage begins and the concentration returns to its minimum value at . Because of the inertia of the system the vessel relaxes to its maximum radius at . The relaxation time , and the total period of the contraction cycle is . T can be adjusted by the rate constant , while T can be adjusted by , and . h controls the chemical reaction speed – that is, it can adjust T.

Figure 7(c) shows the fluid flux averaged across the vessel lumen. The large positive flux peak is preceded by a brief negative flux at the onset of contraction. This is because the right side of the vessel contracts first due to the lower NO concentration in this region: because the fluid flows in from the left (during the relaxation phase in Fig. 7(b)), more NO is produced in that part of the vessel, and this delays the contraction in the left hand side relative to the right (see Fig. 8(a), ). As a result, the higher on the right side initiates the contraction there. The fluid is then pushed to the left briefly (at the reference point of the vessel midpoint, which is plotted in Fig. 7(c)). The backflow is also evident in the velocity vectors, which are directed to the left when a valve is closing (Fig. 8(b), ).This behavior has been reported by other researchers[22,67].

Distribution of NO, pressure and deformation of the vessel and valve

As the vessel begins to contract, the left valve closes and the right valve opens, and the resulting changes in wall shear stresses affect NO production (Fig. 8(a)). The movement of the valve further increases the shear stress on both surfaces of the leaflet, generating high levels of NO near the valves. (see Fig. 8(a), ). Because of the backflow, some NO flows out from the left valve (Fig. 8(a), ). The high concentration of NO exits from the opened right valve, as is evident from the convex shape of the NO gradient surrounding the right valve (Fig. 8(a), ) (see also Fig. 8(b) for flow velocities). The higher concentrations of NO near the valves predicted by the simulations have been observed experimentally in lymphatic vessels in rats[25].

During lymphangion relaxation, (Fig. 8(a), ), the left valve opens and the right valve is closed. Fluid flows into the vessel from the left, and more NO is generated in this region compared with the right side. After relaxation, the concentration of NO remains higher. As the vessel reaches peak diastole, (Fig. 8(a), ), the fluid velocity becomes small, and little NO is generated. NO rapidly degrades or diffuses, and its concentration drops. Meanwhile, the concentration is lower than the threshold, and the vessel is primed for another contraction. The contraction produces high shear stresses and generates NO, which floods the vessel (Fig. 8(a), ).

Because of the smaller effective vessel diameter between the valve leaflets, the shear stresses are largest at the valves (Fig. 8(b), ). This, and the fact that both surfaces of the leaflet can produce NO, is responsible for the higher generation of NO shown in Fig. 8(a). In Fig. 8(a), , the vessel is relaxing and the left valve is open, resulting in fluid flowing in from the left side. Thus, the left side generates more NO than the right side. In Fig. 8(a), , the fluid almost fills the vessel and the fluid velocity becomes small everywhere, resulting in low NO generation. Meanwhile the concentration is lower than the threshold and is primed for another contraction. During the contraction, NO is generated, and it diffuses and convects, filling the vessel (Fig. 8(a), ). In the initial lymphatic segment at left, the porous wall with simulated check valves allow some lymph leakage back to the interstitium (in the junction regions, 85 percent of leakage is reflected, while 15 percent passes).

Figure 8(b) also shows the pressure color map corresponding to Fig. 8(a). Unlike NO, fluid can not pass though the vessel wall or the valve leaflet. Thus, there are pressure discontinuities across these structures. As the vessel contracts, the pressure inside is higher than outside (Fig. 8(b), ). In the initial lymphatic segment at left, the biased valves in the wall maintain the pressure higher than the surrounding tissue. As the contraction ends, the pressure inside the vessel decreases (Fig. 8(b), ). As the intraluminal pressure reaches a minimum (Fig. 8(b), ), flow enters from the left. The decreasing pressure difference across the valve leaflet in Fig. 8, , causes the valve to return to its rest position. Because of the small pressure difference between the inlet and outlet, the fluid velocity is very small, and the pressure is almost homogeneous, implying that the vessel almost completely relaxes. As another contraction starts, the left valve closes and the lymphangion pressure increases again (Fig. 8(b), ). Retrograde flow at a valve can occur either from a contraction in the lymphangion(s) downstream, or by relaxation of the lymphangion(s) upstream from the valve. In Fig. 8(c) we can see a rapid backflow due to the contraction starting at the downstream region of the lymphangion. In this case, the spreads rapidly from downstream to the upstream. This is in contrast to the situation when the backflow is caused by dilation of the upstream vessel[67], which results in slower backflow, and is also seen in our simulations when the outlet valve is closing (see Fig. 8(b), , between the right valve and outlet).

Pumping cycle is affected by the pressure difference between the tissue and the lymphangion outlet

Our simulations predict that lymphatic contractions depend on NO concentration dynamics, which depend on the changing fluid stresses. Thus, adjusting the pressure drop from the tissue to the lymphatic can affect  the fluid velocity and the pumping behavior. To see how the lymphangion responds to changes in tissue pressure that might be encountered in vivo, we performed simulations with step changes in pressure drop. Holding the inlet density , we varied the outlet density as follows:

Thus the pressure difference is

To reproduce valve performance characteristics reported for actual lymphatic vessels[68,69], we allow some backflow through the valves. We introduce an elastic force when the distance between the two valve leaflets is smaller than as:where y′ is the distance between the leaflet and the vessel center line. 2y′ is the distance between two leaflets. Here we choose lattice units and . As shown in Fig. 9, when the left side valve tries to close, this force will maintain a small gap.

Figure 9

Backflow through the gap between the leaflets influences NO distribution.

The dynamics are shown in Fig. 10(a). As expected, in the range of , , the periodic contractions are self-sustaining. In the range of , the forward pressure difference is . In this case, the pressure drives flow and increases the shear stress on the endothelial surfaces to generate high levels of NO. As in Fig. 11, high concentrations of NO depress levels below the threshold, and the lymphangion stops contracting (see Fig. 10(b)). The radius of the vessel is smaller than the rest radius (), implying that the pressure inside the vessel is lower than outside. In the range of , the forward pressure difference changes to , meaning that the outlet pressure is higher than inlet. Because of backflow, The high pressure forces some fluid back into the vessel and expands the vessel. When the diameter of the vessel exceeds R0, the nonlinear term in Eq. (7) begins to activate, creating high levels of and making . This triggers the next contraction (see Fig. 10(a,b)) which strongly forces fluid out of the lymphangion (see Fig. 10(d)). Meanwhile backflow generates NO, which delays from reaching the threshold. Thus the contraction phase is extended as in Fig. 12(a).

Figure 10

Pumping state changes with pressure difference , , . (a) Concentration of , (b) diameter near the outlet of the lymphangion, (c) flux changing with time, and (d) shows the plots in (b) for and , shifted to the same baseline time and diameter to show the differences in period and amplitude. (e) Expanded view of the plots in (c), showing details of the waveforms. , , are the average flux for , Δp2, Δp3, respectively. See also Supplemental Video S3.

Figure 11

NO distribution in the relaxed state. ().

Figure 12

Pumping cycle and average fluxes are affected by Δp.

Figure 10(c) shows how flux is affected by the pressure difference. When , the vessel always starts to contract from the right side, as in Fig. 7(c), resulting in a negative spike of flux. When , the vessel does not contract, resulting in a constant flux. When , the outlet pressure is higher than inlet. The vessel expands beyond R0 and contracts strongly. In this case, the backflow through the right valve during relaxation is sufficient to increase the concentration of NO in this region; therefore, the vessel starts to contract from the left side, resulting in a positive peak of flux at the beginning of the contraction (see Fig. 10(e)). The average flow rate values in Fig. 10 show that an opposing pressure gradient decreases overall flow rate because of the backflow caused by the high pressure at the outlet (see Fig. 10(c), ). Perhaps the most important parameter to reproduce in the simulations is the wall shear rate, which depends on the vessel diameter as well as the fluid velocity. Examining the data in Fig. 10, our shear rate values are comparable to measurements made non-invasively in vivo[67,70]: the flow velocity range in vivo was approximately −100–600 μm/s, in a ~34 μm vessel, which translates to wall shear rates of −23.5 to 141.2 s−1. Our simulations in Fig. 10 fall within this same range, with the wall shear rate varying in the range ~−1.4 to 140.0 s−1.

Analyzing the total output from the vessel, averaged over time, we identify a number of regimes where the vessel exhibits interesting behaviors (Fig. 12):

: When the opposing pressure is this high, the vessel tries to pump, but the flow is retrograde, driven by pressure through the leaky valves.

: At more moderate opposing pressures, the vessel contracts to drive flow, and the output, in general, increases as the opposing pressure decreases. The local maximum in flow at , is due to the generation of NO from backflow at the outlet valve; this NO inhibits contractions near the outlet, and encourages them to start near the inlet. These pseudo-peristaltic contractions tend to increase efficiency slightly. The peak in the period at , is where NO has a large effect on the calcium dynamics, depressing Ca2+ far below the threshold level and thus delaying the next contraction.

: When the pressure gradient is favorable for flow (i.e., the inlet pressure is higher than the outlet), a combination of pressure-driven flow and contractions can move the fluid. Increasing pressure causes a linear increase in flow, roughly according to Poiseuille’s law.

: When the favorable pressure gradient is large enough, the vessel no longer contracts because large amounts of NO are produced. Flow follows Poiseuille’s law in this regime.

Discussion

In summary, we have developed a mathematical model based on tissue mechanics, mecho-sensitive feedback and lymph fluid dynamics that reproduces self sustained cyclic lymphatic contractions and fluid clearance from tissue. To optimize the model, we introduced valves with variable rigidity that are stiffer at the base than the tip. This is consistent with electron microscopic observations that the density of fibers is increased at the valve base[71]. For the fluid dynamics between the leaflets, it was necessary to incorporate lubrication forces, a simple and effective method to estimate forces when two moving fluid boundaries come into close proximity and there is no fluid node between them.

Our simulations reveal a number of interesting implications of the interplay between lymphatic valve structure and mechanobiological feedback control of contractions. Because of the inhibitory effect of NO on contractions, and the fact that its release is dependent on shear stresses, its spatiotemporal dynamics along the vessel wall –especially around valves and during backflow – can affect the propagation of contractions in interesting ways. In actual lymphatic vessels, the diameter is smallest at the valve[56,72], and the shear stress is correspondingly higher in this region. Two other considerations tend to increase NO production at the valve leaflets: first, valve motion results in constantly changing shear stresses, which can increase NO release[73,74] and second, fluid circulation between the valve leaflet and vessel wall means that NO can be released from both sides of the leaflet structure and become concentrated in this area. Thus, the valve performance (e.g., rest position and amount of leakage allowed) can determine lymph flow rates and the apparent direction of contraction propagation. Although the one-way valves ensure positive flow for any sequence of contractions, experimental observations conclude that the contractions can propagate either forward or retrograde with respect to the flow direction[18,19]. The reasons for this were previously unknown. We previously reported that if NO dominates, the contraction propagates retrograde, but if dominates, the contraction propagates forward[15]. Our current results with explicit mechanical valves extend this result, showing that the onset of contraction depends on the NO distribution during diastole. Many factors can shift the NO distribution along the lymphangion during diastole, including valve leakage, the double-sided nature of the valve leaflets and the pressure at the outlet. As these vary, the NO distribution can change, and the contractions initiate wherever NO is lowest. More importantly, these considerations can change the efficiency of fluid or cell transport through the lymphatic system.

A limitation of the model is that we have only a single lymphangion. Nevertheless, it is possible to observe at which end of the vessel the contractions initiate and to determine overall direction of propagation. Another limitation is that the lymphangion has to be anchored to fixed walls at the left and right to establish the boundary conditions. This limits motion of the valve structures somewhat compared to what is seen in vivo. In addition, we set the NO concentration at the boundaries to zero. However, in vivo we would expect there to be dynamic changes in the NO levels in nearby tissue contributed by other nearby lymphatic elements, including upstream and downstream. Future model development will include multiple lymphangions in series, an improvement that will address all these current limitations.

Supplementary video 1

Supplementary video 2

Supplementary video 3