Abstract

1. Introduction

The lymphatic vascular system contains endothelial primary valves in the wall of initial lymphatics, and intravascular secondary valves (Bohlen et al. 2009) which divide collecting lymphatics into lymphangions. The valves are essential for lymph transport, since there is no equivalent of the heart at the upstream end of the system to propel flow. Instead lymphangions, as their name implies, act as pumping chambers in series, utilising either the intermittent passive squeezing of lymphatic vessels resulting from the relative motion of surrounding tissues or the contractions of the unique muscle type in their walls. Like venous valves, lymphatic valves also break up the hydrostatic column that would otherwise create high transmural pressure in lower leg lymphatics when standing.

Unlike cardiac valves, most of the secondary valves are located in very small vessels; even the largest lymphatic vessel in the body is only some 2.2 mm in diameter (Telinius et al. 2010). Working with rat spinotrapezius muscle Mazzoni et al. (1987) found valves in lymphatic vessels as small as 13 μm wide. Consequently the valves must operate at low Reynolds number and in an almost completely quasi-steady flow field. These viscous-flow conditions are thought to be the reason why lymphatic valves have a complex downstream-extending funnel shape. Unlike cardiac valves, where the leaflets insert into the wall at essentially a single axial location called the valve ring, the leaflet insertions of a secondary lymphatic valve extend over an axial distance greater than one vessel diameter (Zawieja 2009). The shape also differs from that of a typical valve in large veins of the leg, which more resembles two pouches butting up to one another, but has much in common with that of valves in micro-veins (Phillips et al. 2004). For more details of valve configuration, see Moore and Bertram (2018). Whereas venous valves are typically more diaphanous than the vessels they inhabit, the walls of small lymphatic vessels are scarcely more substantial than the valve leaflets themselves; see fig. 1 of Zawieja (2009). Mazzoni et al. (1987) proposed a simple mechanistic model to explain how the valve operated, identifying the essential features as being the valve morphology and the flexibility of the leaflets.

Open in a separate window

Figure 1

Views of the three components of the fluid/structure model with 10 μm half-gap. (A) Grey: the closed surface mesh of triangles bounding the valve wall, including the sinus. Red: the fluid boundaries of the valve leaflet. (B) The closed surface mesh of triangles bounding the fluid space in the valve. The curved near surface is rendered with a degree of translucency, so that the interior bounding surfaces where the valve leaflet excludes fluid, and the planes bounding the model, can be perceived. (C): Cyan: the outer surface of the valve leaflet. Yellow: the insertion surface that is bonded to the valve wall.

The sole major study of secondary lymphatic valve function thus far is that by Davis et al. (2011). Their measurements of secondary lymphatic valve function showed that the valves are biased to the open position, to an extent depending on transmural pressure, and that they display opening/closing hysteresis. However, the closure measurements, which involved backflow through the valve, did not account for pressure drops due to micropipette resistance, and so are of uncertain accuracy. Thus the extent of the hysteresis remains unknown at this point. Davis has gone on to evolve his valve tests into an effective tool for assessing valve incompetence in genetically altered mice (Lapinski et al. 2017).

A numerical model can provide a concise description of valve properties and fill in gaps in knowledge from experiments. The simplest possible model of a secondary lymphatic valve simply prohibits backward flow and assigns constant resistance for forward flow (Reddy et al. 1975; Venugopal et al. 2007; Kunert et al. 2015). Macdonald (2008) added finite opening and closing time, but this had more to do with stabilizing the numerical scheme than emulating reality more closely. Bertram et al. (2011a) introduced a model in which valve resistance varied sigmoidally with the trans-valvular pressure difference. This utilized four constants to define minimum and maximum valve resistance, the pressure difference offset from zero on which the transition was centred, and the steepness of the transition. Later, this model was elaborated (Bertram et al. 2014a) to include both the dependence of the offset or bias on the valve transmural pressure and the hysteresis which had been observed experimentally (Davis et al. 2011). The performance of the resulting model subsequently threw doubt on the closure measurements, and motivated attenuation of the valve-closure offset (Bertram et al. 2014b). Adapting a pre-existing circulatory model to lymphatic vessels, Contarino & Toro (2018) used the Mynard cardiac valve model, thereby reintroducing parameters for opening and closing time to the lymphatic valve field. However, lacking data for these parameters, they resorted to guessed values¹.

All the above valve models were lumped-parameter descriptions of function only; they did not attempt to explain how the valve architecture led to a given behaviour. The task of understanding why the valves behave as they do requires a higher-dimensional numerical model. Macdonald (2008) conducted limited 2D finite-volume simulations of steady flow through idealised rigid lymphatic valve geometries, showing how small valves and creeping flow inhibit fluid recirculation behind the leaflets.

Recently, Li et al. (2019) simulated lymphatic valves in 2D again, using the lattice-Boltzmann method. The valves had leaflet stiffness decreasing from the root to the free edge by a factor of almost 2 million, and were incorporated in a 2D model of a lymphangion with flexible walls. It is difficult to picture leaflet stiffness varying spatially in reality by this amount, and it is likely that extreme property variation was more a response to the artefactual difficulties of obtaining realistic lymphatic valve behaviour from a 2D flexible model. The model lymphangion did not conform to what has been measured (Davis et al. 2011); it ceased to produce forward flow when the adverse pressure difference was only 0.025 cmH₂O.

Wilson et al. (2013) modelled the space inside a lymphatic valve in 3D, using data from a stack of confocal images of a rat mesenteric lymphatic. Both steady and time-dependent flows through this space were computed using a finite-volume code and assuming rigid boundaries. Streamlines and wall shear stress (WSS) distributions were computed, and the focus was on determining the concentration distribution of nitric oxide produced in a WSS-dependent rate at the flow domain boundaries, i.e. the vessel wall and valve leaflets.

Wilson et al. (2015) defined an idealised geometry for a 3D lymphatic valve based on measurement of critical dimensions in confocal images from 74 isolated vessels. Three computations were performed successively. In the first, fully developed flow was simulated with the leaflets 2 μm apart on the plane of symmetry cutting them, to determine a pressure field over each leaflet surface and from those data evaluate an average trans-leaflet pressure difference. Then this pressure difference was applied to leaflets with defined thickness and shear modulus, rigidly attached at the rigid vessel boundary, to determine leaflet deflection. Then flow as in step 1 was applied to the rigid valve with this leaflet geometry, to determine valve resistance defined as the trans-leaflet pressure difference divided by the volume flow-rate. Results were obtained for various values of the ratio of maximum sinus diameter to inlet vessel diameter. They found that orifice area was a poor predictor of resistance, and that the combined resistance to forward flow of valve and sinus was less than that of a uniform duct having the diameter of the inlet when the sinus-to-vessel diameter ratio was 1.39 or more.

Wilson et al. (2018) went on to perform two-way coupled fluid-structure interaction computations on the same geometry, still with a rigid vessel boundary. Despite having specified a neo-Hookean elastic material for the leaflet, with no viscoelasticity, they found less deflection of the leaflet for a given trans-valvular pressure during valve opening than during valve closure, for waveforms which varied trans-valvular pressure from zero to peak and back to zero in 1.5 s, with the bulk of each pressure change occurring in 0.2 s. The source of this effect, which is qualitatively different from the static hysteresis observed by Davis et al. (2011), was not specified, but is likely to have been the combined inertia of the leaflet and fluid.

A similar response to a cyclic pressure gradient was found by Ballard et al. (2018), who applied the lattice-Boltzmann method to create a 3D fluid/structure-interaction simulation of a valve consisting of flexible leaflets that were planar when undeformed, intersecting with and fixed to a rigid cylinder representing the vessel wall. They applied steady pressure differences in both forward and backward directions, and varied the leaflet length and stiffness. They found that resistance to forward flow increased with valve length, but that very short valves, of axial length scarcely more than diameter, failed to close fully when flow reversed.

Thus existing models have not yet explained some of the observed behaviour of real lymphatic valves, in particular what structural features set the extent of the bias, and why the bias to the open position is transmural-pressure-dependent (Davis et al. 2011; Bertram et al. 2014b). Accordingly the present paper reports the first 3D computational fluid-dynamic study of lymphatic valves involving fluid/structure interaction with flexibility of both leaflets and vessel wall. By varying the stiffness of the valve wall material and of the leaflet material separately, the contribution of each to the overall valve properties is defined.

2. Methods

2.1. Description of the model

The fluid/structure-interaction models were realized in the finite-element software ADINA v9.4 (ADINA R&D Inc., Watertown, MA, USA), utilizing dimensions based on confocal images of real valves as defined by van Loon (Wilson et al. 2015). Two different geometries were modelled, with the difference being the distance between the leaflets on the plane of symmetry cutting them in the unstressed state. For each model, surface mesh-point grids defining the boundaries of the valve wall, and the fluid space within, excluding the valve leaflet, were constructed in Matlab (see Figure 1) and imported as STL files into ADINA-Structures and ADINA-Fluids, where volume meshes were generated for the wall and the fluid. In the case of the leaflet, point coordinates only were imported into ADINA, and all the surface and volume meshing was performed therein, before marrying the resulting leaflet with the tube wall. Symmetry was assumed about two perpendicular cut-planes intersecting on the axis of the valve; thus only a quarter valve was computed. For each geometry, a rigid-wall model was also defined, in which the wall was merely a fixed no-slip boundary to the fluid space.

Starting from the unstressed or resting configuration with no flow, simulated flow through the valve occurred when the trans-valvular pressure difference was gradually ramped up. To open the valve, the inlet pressure was raised; to close it, the outlet pressure. The coordinate position at each end of the vessel was kept constant (corresponding to securing to pipettes experimentally), but elsewhere the entire structure was free to deform. Over 2s, the pressure was linearly ramped up to 1000 dyn cm⁻² in time-steps of 5 ms; given the minute dimensions of the valve, this was sufficiently slow that the flow was quasi-steady, although it was solved with no simplifying assumptions. For runs which involved pressurisation of the whole structure, the inlet and outlet pressures were initially ramped up simultaneously before one of the two pressures was held constant and the other continued to climb. In every case, the initial fluid mesh was used only for one or a few time-steps, then a re-meshing was carried out which concentrated cells in regions of high shear. Re-meshing was repeated regularly, so that the mesh remained well-formed in the face of the deforming valve and well adapted to the changing flow (Figure 2). An initial list of times for re-meshing was prescribed, but additional re-meshings were organized by restarting the run as needed when the rate of geometry change was high. Whereas the initial fluid grid consisted of 84,892 tetrahedral elements in the first fluid model, and 59,462 in the second, the re-meshed grids consisted of in the order of half a million elements. In many instances, particularly when computing the highest-compliance part of the valve opening process, the time-step had to be reduced to achieve convergence or avoid fluid/structure element overlap; time-steps down to 0.5 ms were employed as necessary.

Open in a separate window

Figure 2

The initial volume mesh (upper panel) underwent successive revisions during the simulation, using refinement criteria based on the local flow. These views of the model with 10 μm half-gap show the exterior of the volume meshes at the plane of symmetry x = 0 which bisects each leaflet.

The wall and leaflet were both modelled as Mooney-Rivlin materials, via the neo-Hookean formulation

$$W=\frac{G}{2}(I_1-3)+\frac{\kappa }{2}{(J_3-1)}^2$$

(1)

where W is the strain energy density, G is the shear modulus, l₁ is the volume-preserving first invariant of the Cauchy-Green deformation tensor C, κ is the bulk modulus and the reduced invariant J₃ is the ratio of the deformed volume to the undeformed volume. The first and second terms on the right-hand side of the equation are the deviatoric and the volumetric strain energy density respectively. The shear modulus was 20 or 40 kPa for the wall, and 10, 20 or 40 kPa for the leaflet. In addition models were run which replaced the wall with a rigid no-slip boundary to which the leaflet insertion adhered, corresponding to the work of Wilson et al. (2018), and an extremely high value of G_wall (2000 kPa) was computed, to check that the flexible-wall model converged to the rigid-wall model in the stiff limit. Where the wall was flexible, both materials were given density 10³ kg m⁻³; the bulk modulus was calculated as

$$\kappa =\frac{2G(1+v_{\mathrm{P}})}{3(1-2v_{\mathrm{P}})}$$

(2)

with Poisson ratio v_P = 0.49. As noted by Wilson et al. (2018), a shear modulus of 20 kPa is slightly less than that of arterial elastin; this is commensurate with the imposition in these simulations of relatively small deformations, such that the lymphatic solids are confined to the part of their passive stress-strain characteristic where the properties of elastin dominate. A Poisson ratio close to 0.5 assumes that lymphatic vascular wall and valve materials are essentially incompressible as a result of being water-saturated, like the material of blood vessel walls (Carew et al. 1968); for the same reason their density is very close to that of water.

The fluid was taken as Newtonian, with density (10³ kg m⁻³) and bulk modulus (2×10⁹ Pa) similar to water, but with viscosity of 1.5 mPa s for intestinal lymph (Moore and Bertram 2018). All dimensional properties were entered into the simulations in μm-gm-s units, then internally made dimensionless on the basis of L₀ = 1 μm, U₀ = 10⁴ μm s⁻¹, ρ₀ = 10⁻¹² gm μm⁻³.

The modelled valve occurred in a vessel of inside diameter 100 μm, and the associated sinus had maximum diameter 160 μm and total length 417 μm. In the unstressed configuration, the valve leaflets were separated on the plane of symmetry x = 0 by 20 μm in the first model. This compares with a real valve (Figure 3). In a second model, this distance was halved. A small but important difference between the implementations of the idealised geometry here and by Wilson et al. (2015) is that they assumed the orifice of the valve projected onto a constant-z plane to be elliptical; here, circular arcs define the inboard trailing edges of the two leaflets and thus the initial valve orifice in z-projection (see the left panel of Fig. 3).

Open in a separate window

Figure 3

Left: diagram from the construction of the first model (half-gap 10 μm), seen in projection onto an axial plane, i.e. z-location. Magenta: the outer and inner curves are respectively the insertion and the trailing edge of the inner surface of each leaflet. Straight lines join equi-spaced points on the insertion and the trailing edge. Black: curves at constant azimuth (which therefore appear straight) outlining the vessel inside wall and sinus. Broken blue: the initial separation of the trailing edges of the inner surface of each leaflet in the second model (half-gap 5 μm) is superimposed. Right: a rat mesenteric lymphatic valve imaged by confocal microscopy, looking downstream. From Zawieja (2009).

Both the valve wall and the leaflet were assumed to have uniform thickness of 5 μm. The maximum pressure difference applied during opening was 1000 dyn/cm², or approximately 1 cmH₂O; during closure, the leaflets reached apposition before this value was attained, except in one instance. Closure was computed up to the last possible apposition configuration before actual leaflet contact; thus contact was not explicitly modelled.

2.2. Numerical considerations

All simulations were computed as simultaneous (rather than iterative) solutions of the fluid and structure component models. In the structure model with 5 μm half-gap, the tube wall and valve leaflet consisted of 9,613 and 3,908 11-node tetrahedral elements respectively. In the equivalent 10 μm model, the wall and leaflet had 10,261 and 10,271 11-node elements respectively.

The fluid space was initially meshed with 84,892 tetrahedral cells in the 10 μm model, and 59,462 tetrahedral cells in the 5 μm model. The greatly increased number of cells in the re-meshed fluid grids was controlled by an adaptive-meshing parameter which nominally set the maximum number of cells. In practice the re-meshed grids typically contained some 50% more cells than this setting. In early runs the parameter was set to 800,000. To check for grid independence, a run with this number set to 400,000 was compared, using the four characterising plot parameters that will be shown below in Figs. 6–9. The curves from the two runs merged everywhere to within the thickness of the (thin) plotted line, i.e. the results with the lower value fully overlaid those from the run with the higher value. With the result thus confirming grid independence, the lower value was used for all subsequent runs, including those giving rise to Figs. 6–9.

Open in a separate window

Figure 6

Trailing-edge position, maximum sinus radius, flow-rate, and valve resistance, all as functions of the pressure difference Δp, while the shear modulus G of the valve wall is varied, from the model with maximum valve-leaflet half gap = 5 μm. G_leaflet = 20 kPa for all curves.

Upper left: the maximum inter-leaflet gap occurs at the plane x = 0; half this gap is shown here.

Upper right: the maximum sinus radius at the plane x = 0.

Lower left: the volume flow-rate Q through the valve.

Lower right: the valve resistance Δp/Q.

Open in a separate window

Figure 9

Trailing-edge position, maximum sinus radius, flow-rate, and valve resistance, all vs. Δp, while the initial y-position of the inner trailing edge of the leaflet at the plane x = 0 is varied. Cyan: half-gap = 10μm. Red: half-gap = 5μm. G_leaflet = G_wall = 20 kPa for all curves.

Volume flow-rate was calculated from the computed axial velocity on the centre-line x = y = 0 at the extremities z = 150 and 667 μm of the conduit, assuming a fully developed Poiseuille profile. The assumption of fully developed flow was checked both by comparison of the inlet and outlet flow-rates so calculated, and by comparing computed velocities with assumed velocities on the four radial lines x = 0, z = 150 μm, y = 0, z = 150 μm, x = 0, z = 667 μm, and y = 0, z = 667 μm. These comparisons were performed for the maximum favourable pressure drop of 1000 dyn/cm², giving rise to the maximum flow-rate through the valve. The Reynolds number based on diameter for this flow was 2.8. The departures from the velocities expected in a parabola having the same centre-line velocity reached a maximum discrepancy of 0.73%, but most were very much less. The maximum outlet-minus-inlet difference in computed velocity was 0.77% of the outlet centre-line velocity; on the centre-line it reached 0.75%.

3. Results

3.1. Pressure and velocity fields

The results from a 3D simulation can be interrogated in many different ways; examples of pressure and velocity distributions when the opening pressure difference is 100 dyn cm⁻² are shown in Figures 4 and ​and5.5. The jet that emerges from the leaflet orifice gives rise to recirculation behind the leaflet, but the flow is 3-dimensional, with complex secondary velocities, so fluid is not trapped indefinitely in closed eddies. Velocities in Fig. 4 on the axis of symmetry x = y = 0 reach 6.25 mm/s at the undilated inlet/outlet regions, and a maximum of 11.6 mm/s through the valve aperture, reflecting the imposition of the pressure difference over an axial distance of just 0.517 mm. The velocity vectors in Fig. 4 show a single recirculation eddy in the lee of the valve leaflet; Margaris et al. (2016) observed a similar phenomenon in a rat mesentery lymphatic valve subjected to an axial pressure difference of 1 cmH₂O (981 dyn cm⁻²), of which most would have been dissipated in their perfusion system. The centre-line velocity in their undilated vessel reached 4.5 mm/s where the diameter was about 92 μm; see their fig. 9.

Open in a separate window

Figure 4

Pressure contours (above) and velocity vectors (below) at the plane x = 0. Units: pressure = 10⁴ dyn/cm², velocity = μm/s. The vector arrow length is fixed; speed is indicated by the colour. The imposed pressure difference creating forward flow is 100 dyn cm⁻². The valve aperture at x = 0 has increased from the initial half-gap of 10 μm to 20 μm.

Open in a separate window

Figure 5

Velocity magnitude contours (left) and secondary-flow velocity vectors (right) at the plane z = 336 μm, 72.5 μm before the z-location of the maximum sinus radius. This axial plane is tangent to the inner trailing edge of the undeflected leaflet, and transects the deflected leaflet upstream of the trailing edge.

3.3. Variation of valve leaflet stiffness

Whereas the shear modulus of the vessel wall was varied over a wide range, including up to infinity for the sake of comparison with a previous model, the shear modulus of the leaflet was varied only over a modest range relating to that of arterial elastin (Wilson et al. 2018). Consequently the dramatic qualitative changes in curve shape seen in the upper left panel of Fig. 6 were not seen when the leaflet stiffness was varied. However, quantitatively the valve is more sensitive to variation of the stiffness of the leaflet than to variation of that of the wall; see Figure 7, where G_wall = 20 kPa, and compare the extent of the curve changes between G_wall = 20 kPa and G_wall = 40 kPa in Fig. 6. The maximum valve opening here (trailing edge y-position 35.4 μm when G_leaflet = 10 kPa, purple curve) is greater than any attained in Fig. 6, but the leaflet here starts from a more open position (10 μm half-gap). Across all simulations reported here, the maximum change in trailing-edge position occurred when the half-gap was 5 μm initially and became 32.4 μm at maximum opening (Fig. 6, red curve).

Open in a separate window

Figure 7

Trailing-edge position, maximum sinus radius, flow-rate, and valve resistance, all vs. Δp, while the shear modulus of the leaflet is varied (values as in legend), from the model with maximum leaflet half gap = 10 μm. G_wall = 20 kPa for all curves. The inset to the lower right panel shows the valve resistance during opening on an enlarged scale, compared with the Poiseuille resistance of an equal-length cylindrical conduit of diameter 100 μm. Figs. 4 and ​and55 illustrate the flow at Δp = 100 dyn/cm² from the curves here for G_leaflet = 10 kPa (purple).

Because the valve leaflets span the lymphatic vessel conduit, their stiffness has a considerable effect on the deformation of the sinus when pressure is applied at valve inlet or outlet to open or close the valve respectively (Fig. 7, upper right). The maximum sinus radius varies through both overall distension of the flexible conduit and variation in the degree of sinus cross-sectional circularity.

Again here there is a retrograde flow-rate maximum during closure (Fig 7, lower left). For leaflet stiffness 10, 20 and 40 kPa, the maximum is −16.4, −23.9 and −33.7 μL/hr, at pressure difference −80, −124 and −178 dyn/cm², when the trailing-edge half-gap at x = 0 is 4.64, 4.57 and 4.72 μm respectively.

The maximum valve resistance Δp/Q reached here (Fig. 7, lower right) was 9.53 × 10⁷ dyn s/cm⁵, at G_leaflet = 40 kPa (grey curve), or about 12.5 times the interpolated resistance at Δp = 0. However, the orifice flow is not properly resolved where the inter-leaflet gap approaches zero, so the computed resistance is not reliable in this limit. Note that there is still backward flow at this extremity, a point that will be returned to below. Fig. 7 also shows (lower right, dashed line) the Poiseuille resistance that a straight cylindrical conduit having the diameter of the inlet and outlet of the modeled valve but no leaflets and no sinus would have. Because of the sinus, the additional valve resistance falls to zero at favourable pressure drops, even becoming slightly negative (see inset). The minimum open-valve resistance reached in this study was 2.59 × 10⁶ dyn s/cm⁵, for the model with 10 μm half-gap and G_leaflet = 10 kPa, G_wall = 20 kPa, at the maximum forward pressure difference of 10³ dyn/cm².

3.4. Variation of valve Inflating pressure

The model with 5 μm half gap and G_leaflet = G_wall = 20 kPa was made the subject of an experiment in which the two boundary pressures were simultaneously ramped up to 500 dyn/cm² (≈ 0.51 cmH₂O) or 1000 dyn/cm² (≈ 1.02 cmH₂O) before any differential pressure was applied. The same procedure as usual was then applied, i.e. in one run the inlet pressure was ramped up a further 1000 dyn/cm² while outlet pressure remained unchanged, and in the other run the outlet pressure was ramped up while the inlet pressure stayed fixed. In the case of the 1 cmH₂O initial inflation, the maximum sinus radius was thereby increased from 80 μm to 81.8 μm, and the inter-leaflet gap on the plane of symmetry also increased, from 10 μm to 16.8 μm; see Figure 8. Because of these enlarged dimensions, the valve admitted slightly more flow-rate for a given pressure difference, in both forward and backward directions². This in turn translated into a slightly lower resistance at a given pressure difference, which can alternatively be viewed as shifting the resistance characteristic to the left, i.e. to greater adverse pressure differences. Thus, relative to the uninflated valve (Fig. 8, lower right, red curve), the inflated valve (beige and light blue curves) opens at a more negative trans-valvular pressure, i.e. its bias to stay open is increased.

Open in a separate window

Figure 8

Trailing-edge position, maximum sinus radius, flow-rate, and valve resistance, all vs. Δp, while the initial inflation pressure is varied. Red: as Fig. 6. Beige: Δp_inflation = +500 dyn/cm² (≈0.5 cmH₂O). Light blue: Δp_inflation = +1000 dyn/cm² (≈1 cmH₂O). G_leaflet = G_wall = 20 kPa for all curves.

3.6. Comparison with lumped-parameter modeling

Using a lumped-parameter model, Bertram et al. (2017) showed that realistically low rates of leakback flow are achieved with a closed-valve resistance much higher than can be computed here. The maximum resistance of a competent valve is a parameter for which the lower bound is around 10¹⁰ dyn s/cm⁵, established by simulations involving slow leak-back through a closed valve over many cycles of lymphatic contraction; see, e.g., figs. 9 and ​and1010 of Bertram et al. (2017). On the other hand, the upper limit to closed-valve resistance in the simulations here is set by the topological change that occurs when the inter-leaflet gap on the symmetry plane x = 0 reaches zero. At this point there is still space between the leaflets off the symmetry plane; see Figure 10. As with the crossing-over of curves for the y-position of the trailing edge of the leaflet during valve opening noted in Fig. 6, this is another manifestation of the three-dimensional deformation of the valve leaflet, but here seen during valve closure.

Open in a separate window

Figure 10

Sections of the model with (left) unstressed half-gap 10 μm and (right) at valve closure, showing that a significant gap between the leaflets continues to exist off-centre even while the gap on the symmetry plane x = 0 becomes vanishingly small. G_leaflet = G_wall = 20 kPa. The two inner arcs show the position of the wall in the cylindrical entrance (slight axial taper owing to distending pressure causes the apparent thickening at valve closure); the two outer arcs show the wall position at the middle of the sinus. Only the inner trailing edge of the leaflet is shown, i.e. that edge which contacts the opposite leaflet (not shown) across the symmetry plane y = 0.

Thus the 3D modelling here does not supersede the values for maximum and minimum valve resistance arrived at previously in the lumped-parameter modelling through comparison with experiments (Bertram et al. 2014b, 2017). Nor does it supersede the lumped-parameter modelling of the valve bias to the open state, since the real extent to which lymphatic valves are anatomically slightly open in the unstressed state is unknown. The comparison of numerical and real orifices shown in Fig. 3 would suggest that even a half-gap of 10 μm may underestimate the true value, but further experimental data are needed before the extent of actual bias can be made more precise. However, one parameter in the lumped-parameter modelling can be compared rather well with the present 3D modelling: the constant s_o, which controls the steepness of the open/close transition in resistance with trans-valvular pressure difference via the equation

$$R_{\mathrm{V}}(\Delta p_{\mathrm{V}})=R_{\mathrm{V}\mathrm{n}}+\frac{R_{{\mathrm{V}}_{\mathrm{x}}}}{1+\exp [-s_{\mathrm{o}}(\Delta p_{\mathrm{V}}-\Delta p_{\mathrm{o}})]}$$

(4)

where R_V = valve resistance, Δp_V is the trans-valvular pressure difference, R_Vn = minimum valve resistance, R_Vx + R_Vn is the maximum valve resistance, s_o is the steepness parameter, and Δp_o is the opening pressure-difference offset. Recent lumped-parameter modelling (Bertram et al. 2017, 2018) used s_o = 0.2 cm²/dyn; this value (chosen to limit closing regurgitation and closed-valve leakage) produces very steeply rising transitions. A resistance curve which closely matches that in Fig. 7 for 10 μm half-gap and G_leaflet = G_wall = 20 kPa (cyan curve) is produced with the parameter values R_Vn = 2.7 × 10⁶ dyn s/cm⁵, R_Vx = 10¹⁰ − R_Vn dyn s/cm⁵, s_o = 0.01 cm²/dyn, and Δp_o = −755 dyn/cm²; see Figure 11³.

Open in a separate window

Figure 11

Comparison of valve resistance vs. Δp, as computed via 3D simulation, with a lumped-parameter model of a lymphatic valve. Thick-line cyan curve: 3D simulation with 10 μm half-gap and G_leaflet = G_wall = 20 kPa. Thin-line green and red curves: lumped-parameter model (parameter values as shown in Table 1). Thin dashed black line: Poiseuille resistance for an equivalent-length cylinder of 100 μm diameter.

The result of changing s_o from 0.2 to 0.01 cm²/dyn was tested in the lumped-parameter model of a single lymphangion (Bertram et al. 2017); see Figure 12. Under control conditions as stipulated in the figure [see parameter values above left panels—explanation of all parameters is given in the Appendix], including pumping against an adverse pressure difference of 2.5 cmH₂O, the cycle-average flow-rate was 27.4 μL/hr. At the reduced s_o value, the cycle-average flow-rate increased to 36.3 μL/hr. The intra-lymphangion pressures varied with time during the contraction over a larger range than in the control run, because the more gradually opening valves provided more resistance to flow when nominally open. By itself, this factor would have reduced the overall ejection volume during the cycle, but its effect was more than compensated by the absence of the regurgitation prior to valve closure that was prominent in the control situation. With the reduced s_o value, Δp_o varied between −0.24 and −0.07 cmH₂O for valve V1, with cycle-average −0.13 cmH₂O, and between −0.34 and +0.08 cmH₂O for V2, with cycle-average −0.10 cmH₂O. To match more closely the value of −0.77 cmH₂O, the model was again run with s_o = 0.01 cm²/dyn, but the two curves⁴ of Δp_o(Δp_tmv) which define the thresholds at which a valve switches from open to closed and vice versa were both shifted by −0.6 cmH₂O, i.e. to more negative values of Δp_o, for both valves. Under these conditions (right panels), the cycle-average flow-rate was 36.4 μL/hr, i.e. almost unchanged; a small amount of regurgitation is compensated by reduced valve resistance when open. But Δp_o varied between −0.90 and −0.67 cmH₂O for V1, with cycle-average −0.76 cmH₂O, and between −0.94 and −0.58 cmH₂O for V2, with cycle-average −0.74 cmH₂O⁵. Other small differences attributable to the changed valve-switching characteristics are noticeable; for instance the isovolumic periods at the beginning and end of systole when both valves are closed are slightly lengthened in the middle panels of Fig. 12, and considerably shortened in the right-hand panels. The leftward shift (in terms of Fig. 11) of the valve characteristic created by the offsets in the right-hand panels of Fig. 12 restores the prompt diastolic filling of the lymphangion that was present in the control run but was hindered under the conditions of the middle panels of Fig. 12. This is seen both in the peak value of flow-rate Q₁ through the inlet valve and in terms of the speed with which the intra-lymphangion pressures p₁, p_m and p₂ all reach the inlet reservoir value of 2.5 cmH₂O after the inlet valve opens.

Open in a separate window

Figure 12

Comparison of one complete cycle from three runs of a lumped-parameter model of a single contracting lymphangion. Left panels: control run with s_o = 0.2 cm²/dyn; see panel surtitles for values of other parameters (CGS units unless stated, resistance in millions). Centre panels: s_o = 0.01 cm²/dyn. Right panels: s_o = 0.01 cm²/dyn and valve open/close thresholds shifted by −0.6 cmH₂O. Each column: top panel shows pressures, second panel shows lymphangion diameter D (blue) and activation time-course M (cyan), third panel shows inlet (blue) and outlet (green) flow-rates, bottom panel shows timing of valve opening and closure (blue = V1, green = V2). The plotted pressures are the upstream and downstream boundary pressures p_a and p_b respectively (thin black lines: Pa,b in legend), and the intravascular pressures p₀ (magenta, upstream of valve V1), p₁ (blue, downstream of V1), p_m (red, mid-lymphangion), p₂ (green, upstream of V2), and p_p (beige, downstream of V2). Dashed line shows external pressure p_e. Traces are plotted in the order shown in the legend; earlier-plotted traces may be largely overwritten by nearly identical ones plotted later.

4. Discussion

The simulations predict that, during that part of valve closure which takes the valve from its rest configuration to first leaflet apposition, there is a retrograde (leakage) flow-rate which goes through a maximum before decaying towards zero. Because it is slightly open at rest, the valve requires an adverse trans-valvular pressure difference to cause the leaflets to deflect enough to achieve closure. This necessarily invokes retrograde flow, which is eventually gradually suppressed as increasing degrees of closure are achieved. Transient backflow in lymphatics with valves was observed by Dixon et al. (2006) and Blatter et al. (2016). The specific prediction will be difficult to verify experimentally, where it is hard to attain perfectly smooth control of the tiny pressures involved, and closure tends to become an unstable process, occurring suddenly⁶. The leaflets are then immediately pushed together in a transient which may momentarily bring unsteady dynamics to the situation.

With rigid wall, the model with 5 μm half-gap and G_leaflet = 20 kPa had minimum open-valve resistance of 3.26 × 10⁶ dyn s/cm⁵. This compares with the value of 2.68 × 10⁶ dyn s/cm⁵, computed by Wilson et al. (2018) in their rigid conduit, with otherwise only small geometry differences. But there is a difference in definition: Δp as used here is the pressure difference from end to end of the conduit, whereas Wilson et al. (2018) used the pressure difference from end to end of the sinus alone. An equivalent value can be approximated here by subtracting the Poiseuille resistance of the cylindrical inlet and outlet tubes, giving 2.65 × 10⁶ dyn s/cm⁵.

In agreement with Wilson et al. (2015), the computations here showed that the structure simulated, which included both the valve proper and the sinus, reached an overall resistance in the open state which was as low or lower than that of a uniform cylindrical conduit having the diameter of the inlet and outlet tubes (3.16 × 10⁶ dyn s/cm⁵). However, Wilson et al. used a rigid tubular geometry, whereas all the structure here is flexible. The comparison here ignores the distension of the structure when one of the boundary pressures is increased to create a pressure difference. An initially cylindrical but flexible tube would also distend; thus overall the resistance of the open valve may not become less than that of the equivalent lymphatic vessel with no valve.

The only measurement of open-valve resistance is the value 0.6 × 10⁶ dyn s/cm⁵, measured by Bertram et al. (2014b). However, the experimental value was for the added resistance of the valve mechanism itself, over that of an empty equivalent conduit. It does not include any component of viscous resistance along the conduit at all. To make approximate comparison with the computations here, it is necessary to subtract the Poiseuille resistance of a valve-less tube having the dimensions of the sinus plus the inlet and outlet tubes. Using the rigid-wall figure from above, the result is 1.48 × 10⁶ dyn s/cm⁵. Using instead the minimum open-valve resistance for the flexible-wall model with 10 μm half-gap (section 4.3), the result is 0.81 × 10⁶ dyn s/cm⁵. That the experimental value is some 25% lower may relate to the numerical valve here, based on the dimensions defined by Wilson et al. (2015), not being able to open as fully as the real valves measured by Bertram et al. (2014b), which were not part of the Wilson dataset.

Fig. 10 showed that, despite the application of significant adverse pressure drops (reaching 425 dyn/cm² in the case of the simulation with 10pm half-gap and G_leaflet = G_wall = 20 kPa⁷), the valve orifice never closed entirely, with a narrow gap remaining on each side of the valve centre-line. It is not known whether such configurations ever occur transiently in vivo. It is likely that in real valves any remaining off-centre spaces are progressively squeezed to ever-smaller dimensions by more adverse pressure difference than can be applied here. Many other factors may operate to allow the real valve to achieve full leak-free closure; these may include subtle small-scale spatial material-property changes where the leaflet trailing edge attaches to the vessel wall in the sinus. A lymphatic valve leaflet consists of a monolayer of endothelial cells on either side of an extracellular matrix which was reported by Rahbar et al. (2012) to be composed predominantly of elastin, with collagen bands at the insertion points. Clearly a homogeneous and isotropic hyper-elastic solid describing simple rubber-like materials, as used here, cannot be expected to emulate the fine detail of the in vivo contact conditions when two such leaflets are pressed against each other by a substantial adverse pressure drop across the valve.

The simulations here suggest that the value of one of the parameters in the lumped-parameter model of a secondary lymphatic valve (of comparable dimensions pertaining to rat mesentery collecting lymphatics) should be modified. Specifically s_o, the parameter which controls the steepness of the open/closed transition in resistance with changing trans-valvular pressure, produces a characteristic which matches that simulated here with 10 μm half-gap and G_leaflet = G_wall = 20 kPa if its value is decreased by a factor of 20 and Δp_o is set to an appropriately changed value. Note that what is approximately fitted here is only a small part of the foot of a valve resistance characteristic that in the lumped-parameter model extends up to a maximum closed-valve resistance of 10¹⁰ dyn s/cm⁵; this is why the fitted values of s_o and Δp_o are so interdependent. When the lumped-parameter model is run with the reduced s_o value (and the valve characteristic offsets set to produce the changed Δp_o value on average), the nature of the pumping cycle changes in small but important ways. Slightly increased open-valve resistance at modest transvalvular pressure drops in the open-valve part of the characteristic does not reduce the mean flow-rate in the trials shown here, involving an adverse pressure difference of 2.5 cmH₂O across a single lymphangion, because reduced regurgitation before valve closure more than compensates. However, it should be emphasized that this result was achieved with the lumped-parameter valve open/close threshold curves shifted by −0.6 cmH₂O relative to those indicated by experiment (Davis et al. 2011; Bertram et al. 2014b). The comparison shows that the 3D simulations are compatible with lumped-parameter valve characteristics that produce a workable lymphangion pump. However, even the detailed 3D valve model is still greatly simplified in terms of material properties, etc., relative to real lymphatic valves in vivo, so this agreement does not guarantee that these shared characteristics are then necessarily closely aligned with the properties of real valves.

To elaborate on some of the simplified aspects, all structures here were assumed neo-Hookean, whereas the passive mechanics of collecting lymphatic vessels is highly non-linear (Rahbar et al. 2012). Partly for that reason, the simulated changes in transmural pressure here (0.5 and 1 cmH₂O) were considerably smaller than the range investigated experimentally by Davis et al. (2011). Furthermore, their evidence suggested that valve function can be influenced by vessel tone, here effectively assumed to be zero or a small constant value.

Lymphatic muscle can generate both tone and phasic contractions. In a lymphangion, the active deformation of the vessel wall generates the pressure within the vessel. In the model here the wall passively deforms in response to externally applied inlet and outlet pressures; this is analogous to the experiments performed by Davis et al. (2011) on lymphatic valves under calcium-free conditions to inhibit vessel contraction. However, muscle cells are scarce (Bridenbaugh et al. 2013a; Zawieja et al. 2018) or absent (Sabine and Petrova 2014) in the immediate vicinity of lymphatic valves. The model here was confined to the vicinity of one valve, with total length 517 μm, or just over five non-sinus vessel diameters.

Another factor that can vary considerably is the valve geometry. The idealised geometry used here is closely based on that defined by Wilson et al. (2015) using dimensions measured in 74 rat mesenteric lymphatic vessels. However, their fig. 2 showed that each dimensional parameter varied quite widely, even in the one tissue and species, while Ballard et al. (2018) show images comparing valve morphology across four locations in the rat and one in the sheep.

This paper has shown how the characteristics of a lymphatic valve are affected by both the stiffness of the leaflets and that of the vessel wall into which the leaflets insert. For the first time, they confirm experiments (Davis et al. 2011) which showed that the valve opens and closes at greater adverse pressure difference if distended by transmural pressure. They also show, again for the first time, how the extent of the valve bias to remain open until significant adverse pressure difference is a function of the extent of the gap between the leaflets in the unstressed condition, with all other geometric parameters unchanged. This work extends our knowledge of how the naturally occurring geometry creates a valve structure which functions competently even at the microscopic spatial scale of typical lymphatic collecting vessels.

Acknowledgements

Appendix - equations of the model

Valve characteristics:

the twin curves that determine Δp_o (eqs. 4 and A4) for opening and closure of each valve are shown in Figure A1.

Figure A1

Open in a separate window

Variation of Δp_o with valve transmural pressure. Solid black and magenta: curves fitted to data measured by Davis et al. (2011) for opening and closing respectively. Broken magenta: closure curve, times constant cfact (Bertram et al. 2014b). Solid and broken green: opening and closing characteristics respectively, offset by −0.6 cm H₂O as for right-hand panels of Fig. 12. Pressure p_int = p_in for opening, p_out for closure.

Footnotes

References