Authors: Huabing Li, Yumeng Mei, Nir Maimon, Timothy P. Padera, James W. Baish, Lance L. Munn
Categories: Article, Computer modelling, Lymphatic vessels
Source: Scientific Reports
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 Ca^2+^ -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 Ca^2+^ fluxes. The model reproduces NO and Ca^2+^ 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.
**Subject ** Computer modelling, Lymphatic vessels
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 Ca2+ 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.
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 Ca2+ in the lymphatic muscle cells. Briefly, Ca2+ 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 Ca2+ systems^15,16^. The details of the model are given in the Methods section.
Figure 1 Diagram of the lymphatic vessel. R
0= 100 μm is the relaxed radius of the vessel.
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 as
where τ is the relaxation time. fi(x,t) is the distribution function at time t and location x. ei is the discrete velocity. fieq 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, W0=4/9, W1∼4=1/9, and W5∼8=1/36. ρ 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 δij is the Kronecker delta function and i,j=x,y. 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 fi on S in Eq. (5) must be known, and fi 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 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 ξ=0.85, 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. 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 CNO is the concentration of NO. DNO is the diffusion coefficient of NO. ∇2CNO(x,t) can be approximated using a second order finite difference
where δ=1 is the discretized distance. u⋅∇CNO can also be approximated using an upwind scheme
In Eq. (7), λ is the chemical reaction rate constant. λΔt determines the chemical reaction time scale. Δt is the lattice time scale. As λ increases, the chemical reaction rate increases. KNO− and KNO+ are the decay rate and production coefficient of NO. Here we consider that the production of NO is proportional to the stress component |∂vl∂xn|^16,44,45^. vl indicates fluid velocity along the tangential direction of the membrane surface and xn 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 ∂vl∂xn using the hydrodynamic force dFl we calculate in Eq. (6). Thus ∂vl∂xn=dFl/(ρυds).
Figure 2 Hydrodynamic force dF acting on surface element ds.
Ca^2+^ 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 Ca^2+^. DCa is the effective diffusion coefficient of Ca2+ spreading from one cell to the neighboring cells along the vessel wall. Note that the speed of the Ca2+ 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^. ∇2CCa can also be approximated through Eq. (8), but here, δ should be changed to a cell length of 1.51 on the lattice. KCa− is the decay rate of Ca2+. The total decay rate of Ca2+ is multiplied by (1+KCa,NOCNO) based on the fact that NO acts on the lymphatic muscle cells through myosin light chain phosphatase to reduce force generation^51–53^. KCa+ is the production rate of Ca2+. KCa+((R−Rl)(RCa−Rl))11 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. Rl 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^. Kδ+ defines a steep increase of CCa while CCa>Cth.
There are five forces that act on the wall of the
The minus sign means that Fr always acts in the direction opposing wall velocity v.
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 fEv, bending force fBv and viscous force frv 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 is
where y0 indicates the rest or limit position of the valve. H is the rest position of the vessel, x0 is the anchor point of the valve. ‘−’ and ‘+’ indicate upper and lower valve leaflets. B can adjust the rest position to be biased to stay open or close. For the limit position, B is maximum. In the open configuration, the fEv and fBv of the valve are designated
where i indicates a discrete segment of the valve. When the two membranes of a valve come into close proximity, we also increase fEv by multiplying an 11th power
where yc is the center line of the vessel (dash line in Fig. 3(b)). When the segment arrives at the limit position yl shown as in Fig. 3(b), we also increase fEv
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 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
where i indicates segment number and n is the total number of segments. A membrane of the valve is discretized from left to right into n segments. k0v is the maximum of kBv. kRv is an approximation of the minimum of kBv, the coefficient A adjusts the soft range of the valve. Figure 4 shows how kBv 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
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
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 δt=1/100 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.
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 FNlub (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 by
where p is the pressure, ηs is the dynamic viscosity, and ux is the fluid velocity along the x-direction. The boundary conditions are
Integrating Eq. (22) of y
where:
From Eq. (23)
The flux
The fluid volume between two planes in Fig. 5 is
During virtual variation of δt, the variation of h0
then the variation of the volume
The flux caused by the relative movement in the normal direction
The flux Q′0(uf) is caused by flow from upstream. Because Eq. (22) is linear and the relationship between Q and ux is also linear, the general flux can be written
Because ux≫uy, ∂ux∂y≫∂ux∂x, the Rayleigh’s dissipation function is^64^:
then the normal lubrication force
Let vb=−va (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
If some segments move up (vy≥0) and others move down (vy<0), for example, n1 segments move up and n2 segments move down, the lubrication force exerted on segment si can be estimated
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++.
The system time scale depends on the fluid kinematic viscosity and diameter of the vessel T=νν′(D′D)2, where ν=(2τ−1)/6, D=2R0 and ν′, D′ are the actual fluid kinematic viscosity and diameter, respectively. T is the duration of one lattice step. The space scale depends on L=D′D, where L is the length of one lattice unit. Here we choose τ=0.75, ν′=0.01cm2/s, D=25 and D′=0.01cm. Thus, real time and space scale are T=1.33×10−6s and L=0.0004cm. The chemical reaction speed depends on λT. In our simulation, we choose the simulation parameters as given in Table 1.
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 Ca2+ concentration within the vessel wall is set to 0.0149, which is close to, but less than, Cth. The density at the boundary (red area in Fig. 1) is kept constant as ρin by imposing an equilibrium distribution with a velocity of 0, and the density at the outlet is kept constant as ρout 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 Ca2+ 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 ρin=1. Since the fluid is treated as water, the pressure unit is P=(L/T)2g⋅cm−3=9.045×104g⋅cm−1⋅s−2, thus the inlet pressure is held at 3.015×104g⋅cm−1⋅s−2. 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 E=L2/T2g⋅cm−3. 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 1.51KBE and 1.56k0vE, giving 13658 and 28266 dynes/cm^2^ 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 10^6^ dynes/cm^2^, but are reasonable, given the much smaller size of our valves. Table 1 lists the model parameter values and their origin.
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 B=0.6L−1, 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 A=5, 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 A=0.5, increasing B of Eq. (14)). The inset figure shows how B affects the rest position of the valve.
Figure 7 shows how Ca^2+^ concentration, vessel diameter and fluid flux change with time when the vessel undergoes sustained spontaneous contractions with ρout=0.99999; in this case, the pressure difference between the tissue boundary and lymphangion outlet is Δp=0.3015g⋅cm−1⋅s−2. The changes in Ca2+ concentration agree with previous reports using one-dimensional models^16^. In Eq. (9), the equilibrium of CCa can be estimated as KCa+(1+((R−Rl)(RCa−Rl))11)/KCa−, thus we choose KCa+<KCa− and the equilibrium of CCa is less than 1. Meanwhile KCa+ is large enough to ensure that the equilibrium CCa (when CNO=0) is larger than the threshold Cth. Initially, CNO=0, and CCa increases from the initial value to the equilibrium value. When CCa becomes larger than Cth (dotted line in Fig. 7(a)), in equation (9), δ↑ will be 1, and CCa increases steeply. The Ca2+ 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 CNO 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 CNO cause CCa 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 CCa level recedes below the threshold Cth, the NO concentration also decreases to the baseline. After CCa reaches the minimum, decreased NO allows CCa to increase again, starting another cycle. Because of the rapid diffusion of Ca2+ between neighboring cells, the movements of adjacent vessel segments are coordinated.
Figure 7 The Ca2+ concentration, diameter of the vessel and flux during periodic contractions, Δp=0.3015g⋅cm−1⋅s−2. (a) Ca2+ concentration, (b) diameter, and (c) flux. (d and e) are the enlarged curve of Ca2+ concentration in (a). The horizontal dotted line indicates the threshold of Ca2+ concentration.
Figure 8 (a) NO concentration, (b) pressure and velocity field dynamics during a contraction cycle, Δp=0.3015g⋅cm−1⋅s−2. 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.
Figure 7(b) shows how the diameter at the midpoint of the lymphangion changes with time. When pin>pout, 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 t1 to t2, after the previous cycle of relaxation, NO reaches the baseline, and CCa begins to increase (dashed curve in Fig. 7(a)). When t=t2, CCa just passes through the threshold, the Ca2+ concentration in the wall suddenly increases and a contraction begins. When t=t3, the radius of the vessel reaches its minimum; the duration of the contraction is Tc=t3−t2=0.216s. After the contraction, the relaxation stage begins and the Ca2+ concentration returns to its minimum value at t=t4. Because of the inertia of the system the vessel relaxes to its maximum radius at t=t5. The relaxation time Tr=t5−t3=1.509s, and the total period of the contraction cycle is T=1.725s. Tc can be adjusted by the rate constant KCa,NO, while Tr can be adjusted by KCa+, KNO− and DNO. 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 because the fluid flows in from the left (during the relaxation phase t3<t<t5 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), t=t4). As a result, the higher Ca2+ 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), t=t1).This behavior has been reported by other researchers^22,67^.
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), t=tA). Because of the backflow, some NO flows out from the left valve (Fig. 8(a), t=tB). 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), t=tB) (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), t=tC,tD), 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), t=tD), the fluid velocity becomes small, and little NO is generated. NO rapidly degrades or diffuses, and its concentration drops. Meanwhile, the Ca2+ 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), t=tE).
Because of the smaller effective vessel diameter between the valve leaflets, the shear stresses are largest at the valves (Fig. 8(b), t=tA). 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), t=tC, 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), t=tD, the fluid almost fills the vessel and the fluid velocity becomes small everywhere, resulting in low NO generation. Meanwhile the Ca2+ 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), t=tE). 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), t=tA). 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), t=tB). As the intraluminal pressure reaches a minimum (Fig. 8(b), t=tC), flow enters from the left. The decreasing pressure difference across the valve leaflet in Fig. 8, t=tD, 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), t=tE). 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 Ca2+ 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), t=tC, between the right valve and 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 ρ=1.0g⋅cm−3, we varied the outlet density as
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 2y′0
where y′ is the distance between the leaflet and the vessel center line. 2y′ is the distance between two leaflets. Here we choose y′0=2 lattice units and kl=0.0004. 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 Ca2+ dynamics are shown in Fig. 10(a). As expected, in the range of 0≤t<T1, Δp1=0, the periodic contractions are self-sustaining. In the range of T1≤t<T2, the forward pressure difference is Δp2=6.032g⋅cm−1⋅s−2. 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 Ca2+ levels below the threshold, and the lymphangion stops contracting (see Fig. 10(b)). The radius of the vessel is smaller than the rest radius (R<R0), implying that the pressure inside the vessel is lower than outside. In the range of T2≤t<T3, the forward pressure difference changes to Δp3=−6.032g⋅cm−1⋅s−2, 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 Ca2+ and making CCa>Cth. 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 CCa from reaching the threshold. Thus the contraction phase is extended as in Fig. 12(a).
Figure 10 Pumping state changes with pressure difference Δp1=0, Δp2=6.032g⋅cm−1⋅s−2, Δp3=−6.032g⋅cm−1⋅s−2. (a) Concentration of Ca2+, (b) diameter near the outlet of the lymphangion, (c) flux changing with time, and (d) shows the plots in (b) for t<T1 and T2≤t<T3, 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. Q¯1, Q¯2, Q¯3 are the average flux for Δp=Δp1, Δp
2, Δp3, respectively. See also Supplemental Video S3.
Figure 11 NO distribution in the relaxed state. (Δp2=6.032g⋅cm−1⋅s−2).
Figure 12 Pumping cycle and average fluxes are affected by Δp.
Figure 10(c) shows how flux is affected by the pressure difference. When Δp1=0, the vessel always starts to contract from the right side, as in Fig. 7(c), resulting in a negative spike of flux. When Δp2=6.032g⋅cm−1⋅s−2, the vessel does not contract, resulting in a constant flux. When Δp3=−6.032g⋅cm−1⋅s−2, 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 Q¯ 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), T2≤t<T3). 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):
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 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 Ca2+ 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.
This work was supported by NIH HL128168 (LLM, TPP) and NSFC grant No. 11362005 (HL).
Huabing Li created the computational codes, produced the simulations and wrote the manuscript; Yumeng Mei and Nir Maimon provided useful discussions; Timothy P. Padera provided experimental data, helped with project conception and wrote the manuscript; James W. Baish was involved with project development, analysis of results and writing the manuscript; Lance Munn initiated and guided the project, analyzed results and wrote the manuscript. All authors reviewed the manuscript before submission.
The authors declare no competing interests.
Huabing Li, Email: hbli@guet.edu.cn.
Lance L. Munn, Email: munn@steele.mgh.harvard.edu
Supplementary information accompanies this paper at 10.1038/s41598-019-46669-9.