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.


Model Description
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][29][30] in response to increased shear stress, and contractions of the collecting lymphatic are determined by the concentration of + Ca 2 in the lymphatic muscle cells. Briefly, + Ca 2 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 + Ca 2 systems 15,16 . The details of the model are given in the Methods section.

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 as where τ is the relaxation time. f t x ( , ) i is the distribution function at time t and location x. e i is the discrete velocity. f i eq 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 ν τ = − . (2 1) 6 (2) For the D2Q9 model, one form of the equilibrium distribution is 27 www.nature.com/scientificreports www.nature.com/scientificreports/ i eq i i i 2 2 where, = W 4/9 0 , = ∼ W 1/9 1 4 , and = ∼ W 1/36 5 8 . ρ 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 e u e u f 1 6 1 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 f i on S in Eq. (5) must be known, and f i 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][34][35][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 ξ = . 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: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][41][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][29][30] , which is allowed to convect and diffuse in both the fluid and tissue spaces 15 : can be approximated using a second order finite difference scheme: where δ = 1 is the discretized distance. ⋅ ∇C u NO 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. − K NO and + K NO are the decay rate and production coefficient of NO. Here we consider that the production of NO is proportional to the stress component ∂ ∂ v x l n 16,44,45 . v l indicates fluid velocity along the tangential direction of the membrane surface and x n is the normal direction vector. Within the www.nature.com/scientificreports www.nature.com/scientificreports/ lymphangion, the valve structures can also generate NO. In fact, as shown in Fig. 2, we can indirectly calculate ∂ ∂ v x l n using the hydrodynamic force dF l we calculate in Eq. (6).
Ca 2+ enters and leaves the cytoplasm of the lymphatic muscle cells and can also pass through junctions to neighboring cells [46][47][48][49] ; therefore, we use the reaction diffusion equation, with diffusion restricted to the 1D space of the vessel wall 15 : Here C Ca is the concentration of Ca 2+ . D Ca is the effective diffusion coefficient of + Ca 2 spreading from one cell to the neighboring cells along the vessel wall. Note that the speed of the + Ca 2 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 . ∇ C 2 Ca can also be approximated through Eq. (8), but here, δ should be changed to a cell length of 1.51 on the lattice. − K Ca is the decay rate of + Ca 2 . The total decay rate of + Ca 2 is multiplied by + K C (1 ) Ca,NO NO based on the fact that NO acts on the lymphatic muscle cells through myosin light chain phosphatase to reduce force generation [51][52][53] . + K Ca is the production rate of l l Ca 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 R Ca are the local radius of the vessel and the reference radius of the nonlinear term respectively. R l is the limiting minimum radius of the vessel. δ↑ is an asymmetric Kronecker δ-function. If C Ca increases from below to exceed the threshold C th , this function is set to 1; otherwise, it is set to 0. This is based on the fact that when C Ca reaches the threshold C th , it triggers action potentials mediated by voltage gated and calcium-induced calcium channels 16,54 . δ + K defines a steep increase of C Ca while > C C Ca th . There are five forces that act on the wall of the vessel: 1. 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. 2. The lymphatic muscle force depends on the concentrations of + Ca 2 and NO according to 15 where K M is the coefficient determining the strength of action. 3. Elastic force from the tissue: In order to limit the contraction amplitude and reproduce tissue mechanical resistance, if < + Δ R R l , we simply increase F E by multiplying an 11th power item. Eq. (11) changes to: 4. 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. y m and y n are the y-direction positions of the neighboring points of the vessel. www.nature.com/scientificreports www.nature.com/scientificreports/ 5. Considering the vessel is viscoelastic, a viscous resistance force is introduced as: The minus sign means that F r 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][58][59][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.
The membrane also has elastic force f E v , bending force f B v and viscous force f r v shown in Eqs (11)(12)(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 where i indicates a discrete segment of the valve. When the two membranes of a valve come into close proximity, we also increase f E v by multiplying an 11th power term: where y c is the center line of the vessel (dash line in Fig. 3(b)). When the segment arrives at the limit position y l shown as in Fig. 3 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 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. k v 0 is the maximum of k B v . k R v is an approximation of the minimum of k B v , the coefficient A adjusts the soft range of the valve. Figure 4 shows how k B v changes over the length of the valve leaflet. www.nature.com/scientificreports www.nature.com/scientificreports/ 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 : 2 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.
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 F N lub (which acts in the y direction, because valve segments move only in the y direction) to determine the force on the valves.
According to lubrication theory 63 , the Stokes equation can be represented by  www.nature.com/scientificreports www.nature.com/scientificreports/ where: The flux is: The fluid volume between two planes in Fig. 5 During virtual variation of δt, the variation of h 0 is: then the variation of the volume is: The flux caused by the relative movement in the normal direction is: is caused by flow from upstream. Because Eq. (22) is linear and the relationship between Q and u x is also linear, the general flux can be written as: x , the Rayleigh's dissipation function is 64 : www.nature.com/scientificreports www.nature.com/scientificreports/ 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 ( ≥ v 0 y ) and others move down ( < v 0 y ), for example, n 1 segments move up and n 2 segments move down, the lubrication force exerted on segment s i can be estimated as: L c m 0 0004 . 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 + Ca 2 concentration within the vessel wall is set to 0.0149, which is close to, but less than, C th . 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  . 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 . K E 1 51 B and . k E 1 56 v 0 , 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.

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 www.nature.com/scientificreports www.nature.com/scientificreports/ 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 L 0 6 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 =  www.nature.com/scientificreports www.nature.com/scientificreports/ The changes in + Ca 2 concentration agree with previous reports using one-dimensional models 16 . In Eq. (9), the equilibrium of C Ca can be estimated as 

Dynamics of Ca 2+ concentration.
Ca l l Ca , thus we choose < + −

K K
Ca Ca and the equilibrium of C Ca is less than 1. Meanwhile + K Ca is large enough to ensure that the equilibrium C Ca (when = C 0 NO ) is larger than the threshold C th . Initially, = C 0 NO , and C Ca increases from the initial value to the equilibrium value. When C Ca becomes larger than C th (dotted line in Fig. 7(a)), in equation (9), δ↑ will be 1, and C Ca increases steeply. The + Ca 2 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 C NO 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 C NO cause C Ca 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 C Ca level recedes below the threshold C th , the NO concentration also decreases to the baseline. After C Ca reaches the minimum, decreased NO allows C Ca to increase again, starting A 0 5, increasing B of Eq. (14)). The inset figure shows how B affects the rest position of the valve.  Figure 7(b) shows how the diameter at the midpoint of the lymphangion changes with time. When > p p in out , 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 R 0 . In Fig. 7(b), from t 1 to t 2 , after the previous cycle of relaxation, NO reaches the baseline, and C Ca begins to increase (dashed curve in Fig. 7(a)). When = t t 2 , C Ca just passes through the threshold, the + Ca 2 concentration in the wall suddenly increases and a contraction begins. When = t t 3 , the radius of the vessel reaches its minimum; the duration of the contraction is

Vessel contractions and fluid flow.
. After the contraction, the relaxation stage begins and the + Ca 2 concentration returns to its minimum value at = t t 4 . Because of the inertia of the system the vessel relaxes to its maximum radius at = t t 5 . The relaxation time = − = .
T t t s 1 509 r 5 3 , and the total period of the contraction cycle is = .
T s 1 725 . T c can be adjusted by the rate constant K Ca,NO , while T r can be adjusted by + K Ca , − K NO and D NO . 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 < < t t t 3 5 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 t 4 ). As a result, the higher + Ca 2 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( .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), = t t A ). Because of the backflow, some NO flows out from the left valve ( Fig. 8(a), =  t t B ). 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 t B ) (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 t t , C D ), 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 t D ), the fluid velocity becomes small, and little NO is generated. NO rapidly degrades or diffuses, and its concentration drops. Meanwhile, the + Ca 2 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( . 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 t C , 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 t D , the fluid almost fills the vessel and the fluid velocity becomes small everywhere, resulting in low NO generation. Meanwhile the + Ca 2 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 t E ). 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 t A ). 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( . As the intraluminal pressure reaches a minimum ( Fig. 8(b), = t t C ), flow enters from the left. The decreasing pressure difference across the valve leaflet in Fig. 8, = t t D , 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( . 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 + Ca 2 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 t C , 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 ρ = . ⋅ − 1 0 g cm 3 , we varied the outlet density as follows: 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 ′ y 2 0 as: where y′ is the distance between the leaflet and the vessel center line. 2y′ is the distance between two leaflets. Here we choose ′ = y 2 0 lattice units and = . k 0 0004 l . As shown in Fig. 9, when the left side valve tries to close, this force will maintain a small gap. The + Ca 2 dynamics are shown in Fig. 10(a). As expected, in the range of ≤ < t T 0 1 , Δ = p 0 1 , the periodic contractions are self-sustaining. In the range of ≤ < T t T 1 2 , the forward pressure difference is Δ = .
⋅ ⋅ − − p 6 032 g cm s 2 1 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 + Ca 2 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 R 0 ), implying that the pressure inside the vessel is lower than outside. In the range of ≤ < T t T 2 3 , the forward pressure difference changes to Δ = − . ⋅ ⋅ − − p 6 032 g cm s 3 1 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 R 0 , the nonlinear term in Eq. (7) begins to activate, creating high levels of + Ca 2 and making > C C Ca th . 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 C Ca from reaching the threshold. Thus the contraction phase is extended as in Fig. 12(a). Figure 10(c) shows how flux is affected by the pressure difference. When Δ = p 0 1 , the vessel always starts to contract from the right side, as in Fig. 7(c), resulting in a negative spike of flux. When Δ = .

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 + Ca 2 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.