A one-dimensional mathematical model of collecting lymphatics coupled with an electro-fluid-mechanical contraction model and valve dynamics

TL;DR

This paper introduces a comprehensive one-dimensional mathematical model of collecting lymphatics coupled with an innovative electro-fluid-mechanical contraction model and valve dynamics, enabling detailed simulation of lymphatic function and pathologies.

Contribution

The study presents a novel coupled PDE-ODE model for lymphatic contractions and valve behavior, including stability analysis and simulation of stenosis and regurgitation effects.

Findings

The EFMC model replicates pressure and shear stress effects on contraction frequency.

Stenotic valves significantly reduce pump flow at high contraction frequencies.

Incompetent valves diminish net lymph flow as incompetence severity increases.

Abstract

We propose a one-dimensional model for collecting lymphatics coupled with a novel Electro-Fluid-Mechanical Contraction (EFMC) model for dynamical contractions, based on a modified FitzHugh-Nagumo model for action potentials. The one-dimensional model for a compliant lymphatic vessel is a set of hyperbolic Partial Differential Equations (PDEs). The EFMC model combines the electrical activity of lymphangions (action potentials) with fluid-mechanical feedback (stretch of the lymphatic wall and wall shear stress) and the mechanical variation of the lymphatic wall properties (contractions). The EFMC model is governed by four Ordinary Differential Equations (ODEs) and phenomenologically relies on: (1) environmental calcium influx, (2) stretch-activated calcium influx, and (3) contraction inhibitions induced by wall shear stresses. We carried out a complete mathematical analysis of the stability of the stationary state of the EFMC model. Overall, the EFMC model allows imitating the influence of pressure and wall shear stress on the frequency of contractions observed experimentally. Lymphatic valves are modelled using a well-established lumped-parameter model which allows simulating stenotic and regurgitant valves. We analysed several lymphodynamical indexes of a single lymphangion for a wide range of upstream and downstream pressure combinations. Stenotic and regurgitant valves were modelled, and their effects are here quantified. Results for stenotic valves showed in the downstream lymphangion that for low frequencies of contractions the Calculated Pump Flow (CPF) index remained almost unaltered, while for high frequencies the CPF dramatically decreased depending on the severity of the stenosis (up to 93% for a severe stenosis). Results for incompetent valves showed that the net flow during a lymphatic cycle tends to zero as the degree of incompetence increases.

Tables

Table 1. Table 1: Parameters used for the one-dimensional EFMC model for lymph flow. We adopted the geometrical structure of collecting lymphatics from rat mesentery.

Parameter	Description	Value	Units
	Unknowns
$A$	Lymphatic cross-sectional area	$A(x,t)$	$\mu$m2
$q$	Lymph flow	$q(x,t)$	$\mu$L min-1
$v$	Excitable variable	$v\left(t\right)$	-
$w$	Recovery variable	$w\left(t\right)$	-
$s$	State of contraction ($0\le s\le 1$)	$s\left(t\right)$	-
$I$	Stimulus	$I\left(t\right)$	-
$\xi$	State of the lymphatic valve ($0\le \xi \le 1$)	$\xi \left(t\right)$	-
$q_v$	Flow across the lymphatic valve	$q_v\left(t\right)$	$\mu$L min-1
	Material parameters
$\gamma$	Parameter for velocity profile	2	$-$
$\mu$	Lymph dynamic viscosity	1	cP
$\rho$	Lymph density	998	kg m-3
$\nu$	Poisson ratio	0.5	-
$E_{max}$	Maximum Young modulus	135000	Pa
$E_{min}$	Minimum Young modulus	35000	Pa
$r_0$	Radius at zero transmural pressure	47.7	$\mu$m
$h_0$	Wall-thickness at zero transmural pressure	$0.3r_0$	$\mu$m
$A_0$	Cross-sectional area at zero transmural pressure	$\pi r_0^2$	$\mu$m2
$L$	Lymphangion length	1.5	mm
$p_e$	External pressure	0	cmH2O
	Tube law
$m$	Parameter	0.5	-
$n$	Parameter	-5.0	-
$z$	Parameter	19.0	-
$C$	Parameter	1.0e-16	-
	Electro-Fluid-Mechanical Contraction (EFMC) model
$t_{excited}$	Required time to perform an action potential	$\approx$ 2 $s$	s
$t_{activation}$	Time required to activate an action potential	Eq. (38)	s
$f_{min}$	Minimum frequency without stretch	3.0	min-1
$f_{Ca}$	Maximum frequency at $\tilde{A}=A_{Ca}$	20.0	min-1
$R_I$	Radius of the activation region	0.1	-
$n_{Ca}$	Stretch-activation parameter	10	-
$A_{Ca}$	Reference stretch-activation cross-sectional area	7.75$A_0$	$\mu$m2
$k_{Ca}^{\left(1\right)}$	Baseline increasing rate of stimulus $I$	Eq. (41)	s-1
$k_{Ca}^{\left(2\right)}$	Stretch-activated increasing rate of stimulus $I$	Eq. (41)	s-1
$k_{rel}$	Decreasing rate of the stimulus $I$	10	s-1
$a_1$	Parameter	100	s-1
$a_2$	Parameter	0.5	-
$a_3$	Parameter	25.0	-
$b_1$	Parameter	3.0	s-1
$b_2$	Parameter	0.0	s-1
$c_1$	Increasing rate of contraction state $s$	110	s-1
$c_2$	Decreasing rate of contraction state $s$	3	s-1
$\tilde{I}$	Approximated stimulus required to trigger an action potential	Eq. (43)	-
$k_{NO}$	Contraction inhibition parameter ($0\le k_{NO}\le 1$)	0.8	-
${\tau }_{NO}$	Reference wall shear stress	6.0	dyne cm-2
$n_{NO}$	Wall shear stress inhibition parameter	1.2	-
	Valve model
$\Delta p_{open}$	Valve opening threshold pressure difference	0	cmH2O
$\Delta p_{close}$	Valve closure threshold pressure difference	0	cmH2O
$K_{vo}$	Rate coefficient valve opening	10.0	Pa-1 s-1
$K_{vo}$	Rate coefficient valve closure	10.0	Pa-1 s-1
$B$	Bernoulli resistance	Eq. (50)	cmH2O s2 $\mu$L-2
$L$	Lymphatic inertia	Eq. (50)	cmH2O s2 $\mu$L-1
$R$	Viscous resistance to flow	Eq. (50)	cmH2O s $\mu$L-1
$M_{st}$	Maximum valve opening ($0\le M_{st}\le 1$)	1.0	-
$M_{rg}$	Minimum valve closure ($0\le M_{rg}\le 1$)	0.0	-
$L_{eff}$	Effective length	0.1	mm

Table 2. Table 2: Lymphodynamical indexes. t 1 subscript 𝑡 1 t_{1} and t 2 subscript 𝑡 2 t_{2} corresponds to the initial and ending time of the lymphatic contractions and L 𝐿 L is the lymphangion length.

	Index	Formula	Description	Units
ESD	End-Systolic Diameter	-	Diameter at the end of lymphatic contraction	$\mu$m
EDD	End-Diastolic Diameter	-	Diameter at the beginning of filling	$\mu$m
ESP	End-Systolic Pressure	-	Pressure at the end of lymphatic contraction	cmH2O
EDP	End-Diastolic Pressure	-	Pressure at the beginning of filling	cmH2O
EF	Ejection Fraction	$1-\frac{\text{EDV}}{\text{ESV}}$	Fractional amount of ejected lymph	-
SV	Stroke Volume	$\text{EDV}-\text{ESV}$	Ejected volume amount	nL
FPF	Fractional Pump Flow	$\text{EF}\times \text{FREQ}$	Fractional change in lymphatic volume per minute	min-1
CPF	Calculated Pump Flow	$\text{SV}\times \text{FREQ}$	Flow produced by lymphatic contraction	$\mu$L h-1
AMP	Amplitude	$\text{EDD}-\text{ESD}$	Difference between diastolic and systolic diameter	$\mu$m
SW	Stroke Volume	$\int PdV$	Area inside the pressure-volume loop	nL cmH2O
CPFI	Calculated Pump Flow Index	$\left|\frac{\text{CPF}}{q_{\text{mean}}}\right|$	Ratio between propelled and passive flows	-
${\tau }_{mean}$	Time-averaged Wall Shear Stress	$\frac{1}{t_2-t_1}{\int }_{t_1}^{t_2}\tau (\frac{L}{2},t)dt$	Averaged WSS during a lymphatic cycle	dyne cm-2
$q_{mean}$	Time-averaged Flow	$\frac{1}{t_2-t_1}{\int }_{t_1}^{t_2}q(\frac{L}{2},t)dt$	Averaged flow during a lymphatic cycle	$\mu$L h-1

Table 3. Table 3: Sensitivity analysis of the one-dimensional lymph flow equations coupled to the EFMC model and valve dynamics. Negative pressure gradient . Symbol " − " " " "-" indicates that the absolute value of the sensitivity index S ¯ i , j subscript ¯ 𝑆 𝑖 𝑗 \bar{S}_{i,j} was less than 0.5 and therefore parameter x i subscript 𝑥 𝑖 x_{i} did not influence index P j subscript 𝑃 𝑗 P_{j} . Likewise, green-coloured parameters show an influence between 25 and 50, blue-coloured parameters show an influence between 50 and 100, and red-coloured parameters show an influence greater than 100. Results are shown using as mean ± plus-or-minus \pm SD.

$P_{in}$ = 3 - $P_{out}$ = 4	$𝐏\left(𝐗\right)$		7.80 $\pm$ 5.43	0.74 $\pm$ 0.05	47.19 $\pm$ 14.78	5.75 $\pm$ 4.03	22.87 $\pm$ 16.60	-0.12 $\pm$ 0.09	118.93 $\pm$ 18.05	232.78 $\pm$ 34.08	4.49 $\pm$ 0.15	3.00 $\pm$ 0.00	3.07 $\pm$ 0.05	21.73 $\pm$ 15.79
$𝐗$	$\overline{𝐒}$ [%]		Frequency	EF	SV	FPF	CPF	WSS	ESD	EDD	ESP	EDP	Mean Pressure	Mean Flow
			[min-1]	[-]	[nL]	[min-1]	[$\mu$L h-1]	[dyne cm-2]	[$\mu$m]	[$\mu$m]	[cmH2O]	[cmH2O]	[cmH2O]	[$\mu$L h-1]
47.52 $\pm$ 8.21	$r_0$	[$\mu$m]	-	-	207.7 $\pm$ 61.9	0.8 $\pm$ 1.2	249.1 $\pm$ 170.3	121.6 $\pm$ 82.9	119.2 $\pm$ 18.2	105.6 $\pm$ 15.1	-	2.0 $\pm$ 1.3	-	235.7 $\pm$ 161.6
100.23 $\pm$ 16.92	$a_1$	[s-1]	5.7 $\pm$ 6.6	-0.9 $\pm$ 1.3	-1.1 $\pm$ 1.6	3.2 $\pm$ 5.5	3.4 $\pm$ 6.2	-1.4 $\pm$ 6.5	1.5 $\pm$ 2.2	-	-	-	-	3.0 $\pm$ 7.1
0.49 $\pm$ 0.09	$a_2$	[-]	-76.9 $\pm$ 108.6	22.7 $\pm$ 29.5	26.3 $\pm$ 34.0	-30.6 $\pm$ 73.6	-33.4 $\pm$ 80.8	-4.3 $\pm$ 83.0	-40.1 $\pm$ 50.7	-	-	1.5 $\pm$ 0.7	3.6 $\pm$ 2.3	-33.4 $\pm$ 81.0
25.09 $\pm$ 4.31	$a_3$	[-]	-30.5 $\pm$ 46.4	6.5 $\pm$ 8.2	7.5 $\pm$ 9.4	-16.1 $\pm$ 34.7	-17.0 $\pm$ 36.2	4.6 $\pm$ 32.7	-11.5 $\pm$ 14.2	-	-	0.6 $\pm$ 0.3	1.5 $\pm$ 0.9	-17.1 $\pm$ 36.4
3.00 $\pm$ 0.52	$b_1$	[s-1]	28.5 $\pm$ 49.1	-2.5 $\pm$ 3.0	-3.5 $\pm$ 4.3	19.8 $\pm$ 39.1	20.4 $\pm$ 40.1	-14.4 $\pm$ 33.1	4.0 $\pm$ 4.9	-	-	-	-0.6 $\pm$ 0.5	20.6 $\pm$ 40.4
109.79 $\pm$ 19.28	$c_1$	[s-1]	-0.9 $\pm$ 1.7	4.6 $\pm$ 6.4	5.4 $\pm$ 7.4	4.0 $\pm$ 5.8	4.9 $\pm$ 7.4	-10.8 $\pm$ 13.7	-8.2 $\pm$ 11.0	-	-	-	-	5.5 $\pm$ 8.4
3.01 $\pm$ 0.53	$c_2$	[s-1]	15.3 $\pm$ 20.5	-	-	14.0 $\pm$ 18.7	17.5 $\pm$ 23.8	-13.5 $\pm$ 19.0	-	-	-	-1.3 $\pm$ 0.4	-	17.3 $\pm$ 23.3
0.10 $\pm$ 0.02	$R_I$	[s-1]	4.6 $\pm$ 4.4	-	-	4.6 $\pm$ 4.5	5.6 $\pm$ 5.8	-4.8 $\pm$ 5.4	-	-	-	-	-	5.3 $\pm$ 6.6
10.10 $\pm$ 1.70	$k_{rel}$	[s-1]	2.0 $\pm$ 2.9	-1.4 $\pm$ 2.2	-1.6 $\pm$ 2.5	-	-	1.8 $\pm$ 4.0	2.5 $\pm$ 3.8	-	-	-	-	-0.5 $\pm$ 3.8
9.92 $\pm$ 1.72	$n_{Ca}$	[-]	-52.0 $\pm$ 34.5	-	-	-47.1 $\pm$ 32.4	-51.4 $\pm$ 40.1	45.6 $\pm$ 33.2	-	-	-	-	-1.1 $\pm$ 0.8	-52.3 $\pm$ 40.0
7.78 $\pm$ 1.31	$A_{Ca}/A_0$	[-]	-377.7 $\pm$ 281.3	-	-	-349.8 $\pm$ 257.3	-405.2 $\pm$ 289.7	339.3 $\pm$ 255.3	-	-	-	-	-8.2 $\pm$ 5.8	-406.0 $\pm$ 287.1
0.50 $\pm$ 0.09	$k_{No}$	[-]	-	-	-	-	-	-	-	-	-	-	-	-
5.96 $\pm$ 1.03	$t_{No}$	[dyne cm-2]	-	-	-	-	-	-	-	-	-	-	-	-
1.20 $\pm$ 0.20	$n_{No}$	[-]	-	-	-	-	-	-	-	-	-	-	-	-
9.89 $\pm$ 1.76	$K_{vo}$	[Pa-1 s-1]	-	-	-	-	-	-	-	-	-	-	-	-
10.01 $\pm$ 1.74	$K_{vc}$	[Pa-1 s-1]	-	-	-	-	-	-7.2 $\pm$ 4.5	-	-	-	-	-	5.5 $\pm$ 4.2
99.49 $\pm$ 17.43	$L_{eff}$	[$\mu$m]	-	-	-	-	-	-6.4 $\pm$ 4.4	-	-	-	-0.6 $\pm$ 0.4	-	5.0 $\pm$ 4.9
1007.99 $\pm$ 176.63	$\rho$	[kg m-3]	-	-	-	-	-	-	-	-	-	-	-	-
1.00 $\pm$ 0.17	$\mu$	[cP]	-	-	-	-	-	-114.7 $\pm$ 79.6	-	-	-	-1.0 $\pm$ 0.6	-	4.9 $\pm$ 4.8
2.01 $\pm$ 0.35	$\gamma$	[-]	-	-	-	-	-	-56.7 $\pm$ 39.8	-	-	-	-0.5 $\pm$ 0.4	-	2.6 $\pm$ 3.5
134280.92 $\pm$ 22699.60	$E_{max}$	[Pa]	-5.7 $\pm$ 7.9	29.1 $\pm$ 15.7	33.1 $\pm$ 19.2	24.6 $\pm$ 20.4	28.5 $\pm$ 24.9	-78.6 $\pm$ 55.8	-54.3 $\pm$ 25.8	-	-	-	-	29.7 $\pm$ 26.5
34591.78 $\pm$ 5844.81	$E_{min}$	[Pa]	-195.7 $\pm$ 171.4	-19.6 $\pm$ 17.3	-80.6 $\pm$ 52.3	-193.2 $\pm$ 159.6	-270.3 $\pm$ 191.3	199.8 $\pm$ 166.2	-1.9 $\pm$ 3.5	-32.9 $\pm$ 22.6	-	1.5 $\pm$ 0.8	-3.8 $\pm$ 3.5	-275.6 $\pm$ 192.3

Table 4. Table 4: Sensitivity analysis of the one-dimensional lymph flow equations coupled to the EFMC model and valve dynamics. Positive pressure gradient . Symbol " − " " " "-" indicates that the absolute value of the sensitivity index S ¯ i , j subscript ¯ 𝑆 𝑖 𝑗 \bar{S}_{i,j} was less than 0.5 and therefore parameter x i subscript 𝑥 𝑖 x_{i} did not influence index P j subscript 𝑃 𝑗 P_{j} . Likewise, green-coloured parameters show an influence between 25 and 50, blue-coloured parameters show an influence between 50 and 100, and red-coloured parameters show an influence greater than 100. Results are shown using as mean ± plus-or-minus \pm SD.

$P_{in}$ = 4 - $P_{out}$ = 2	$𝐏\left(𝐗\right)$		5.75 $\pm$ 4.62	0.78 $\pm$ 0.03	50.36 $\pm$ 13.91	4.54 $\pm$ 3.72	17.44 $\pm$ 14.75	-12.98 $\pm$ 2.05	108.40 $\pm$ 15.53	231.99 $\pm$ 31.88	3.37 $\pm$ 0.08	3.01 $\pm$ 0.01	3.02 $\pm$ 0.01	4821.88 $\pm$ 2221.97
$𝐗$	$\overline{𝐒}$ [%]		Frequency	EF	SV	FPF	CPF	WSS	ESD	EDD	ESP	EDP	Mean Pressure	Mean Flow
			[min-1]	[-]	[nL]	[min-1]	[$\mu$L h-1]	[dyne cm-2]	[$\mu$m]	[$\mu$m]	[cmH2O]	[cmH2O]	[cmH2O]	[$\mu$L h-1]
47.59 $\pm$ 8.44	$r_0$	[$\mu$m]	-18.5 $\pm$ 18.4	-	301.1 $\pm$ 82.2	-19.1 $\pm$ 19.3	965.1 $\pm$ 786.6	-47.5 $\pm$ 18.8	108.6 $\pm$ 15.5	118.5 $\pm$ 16.3	-	-	-	328.7 $\pm$ 134.1
99.68 $\pm$ 17.26	$a_1$	[s-1]	14.1 $\pm$ 19.9	-	-0.6 $\pm$ 0.8	13.1 $\pm$ 20.2	17.1 $\pm$ 28.3	-	0.6 $\pm$ 0.8	-	-	-	-	-
0.50 $\pm$ 0.09	$a_2$	[-]	-182.1 $\pm$ 287.9	8.5 $\pm$ 10.4	12.6 $\pm$ 15.6	-152.6 $\pm$ 273.5	-200.6 $\pm$ 342.7	-27.7 $\pm$ 17.0	-15.8 $\pm$ 19.5	-	-	-	-	-13.7 $\pm$ 8.8
24.87 $\pm$ 4.24	$a_3$	[-]	-67.3 $\pm$ 115.0	2.4 $\pm$ 2.9	3.6 $\pm$ 4.4	-60.6 $\pm$ 116.9	-78.7 $\pm$ 147.9	-9.8 $\pm$ 6.0	-4.5 $\pm$ 5.4	-	-	-	-	-4.9 $\pm$ 3.2
2.99 $\pm$ 0.52	$b_1$	[s-1]	55.8 $\pm$ 109.9	-0.9 $\pm$ 1.1	-1.7 $\pm$ 2.1	54.3 $\pm$ 111.5	69.0 $\pm$ 132.9	2.8 $\pm$ 3.5	1.5 $\pm$ 1.9	-	-	-	-	1.4 $\pm$ 1.6
110.24 $\pm$ 18.67	$c_1$	[s-1]	-2.4 $\pm$ 5.1	1.7 $\pm$ 2.2	2.5 $\pm$ 3.3	2.8 $\pm$ 5.5	4.4 $\pm$ 8.3	-1.5 $\pm$ 1.5	-3.2 $\pm$ 4.1	-	-	-	-	-0.8 $\pm$ 0.8
2.99 $\pm$ 0.52	$c_2$	[s-1]	36.6 $\pm$ 59.2	-	-	39.4 $\pm$ 64.4	57.8 $\pm$ 90.5	7.6 $\pm$ 4.5	-	-	-	-	-	3.0 $\pm$ 1.9
0.10 $\pm$ 0.02	$R_I$	[s-1]	10.6 $\pm$ 13.1	-	-	12.3 $\pm$ 15.0	17.4 $\pm$ 21.0	-0.7 $\pm$ 1.0	-	-	-	-	-	-
9.92 $\pm$ 1.66	$k_{rel}$	[s-1]	3.8 $\pm$ 6.4	-	-0.8 $\pm$ 1.1	1.9 $\pm$ 5.3	2.6 $\pm$ 7.4	0.7 $\pm$ 0.9	0.9 $\pm$ 1.4	-	-	-	-	-
10.08 $\pm$ 1.69	$n_{Ca}$	[-]	-132.8 $\pm$ 89.8	-	-	-139.2 $\pm$ 95.4	-194.1 $\pm$ 148.8	7.2 $\pm$ 5.4	-	-	-	-	-	3.0 $\pm$ 2.5
7.75 $\pm$ 1.31	$A_{Ca}/A_0$	[-]	-1272.7 $\pm$ 1070.8	-	-	-1350.4 $\pm$ 1136.1	-1938.7 $\pm$ 1626.2	68.0 $\pm$ 54.9	-	-	-	-	-2.1 $\pm$ 1.7	26.5 $\pm$ 22.1
0.51 $\pm$ 0.09	$k_{No}$	[-]	-94.3 $\pm$ 53.1	-	-	-99.0 $\pm$ 55.8	-145.3 $\pm$ 92.1	5.1 $\pm$ 3.1	-	-	-	-	-	2.3 $\pm$ 1.6
6.00 $\pm$ 1.04	$t_{No}$	[dyne cm-2]	56.5 $\pm$ 34.2	-	-	59.6 $\pm$ 36.3	83.8 $\pm$ 52.1	-3.1 $\pm$ 2.0	-	-	-	-	-	-1.3 $\pm$ 1.0
1.22 $\pm$ 0.20	$n_{No}$	[-]	-25.6 $\pm$ 16.8	-	-	-26.9 $\pm$ 17.6	-39.7 $\pm$ 29.8	1.5 $\pm$ 1.3	-	-	-	-	-	0.6 $\pm$ 0.6
9.98 $\pm$ 1.68	$K_{vo}$	[Pa-1 s-1]	-	-	-	-	-	-	-	-	-	-	-	-
10.06 $\pm$ 1.72	$K_{vc}$	[Pa-1 s-1]
99.81 $\pm$ 17.06	$L_{eff}$	[$\mu$m]	31.0 $\pm$ 20.5	-	-	33.4 $\pm$ 22.5	46.5 $\pm$ 31.9	43.0 $\pm$ 6.1	-	-	-	-0.8 $\pm$ 0.3	-	-50.3 $\pm$ 16.4
999.59 $\pm$ 172.23	$\rho$	[kg m-3]	8.5 $\pm$ 7.4	-	-	9.0 $\pm$ 7.8	13.7 $\pm$ 13.5	11.9 $\pm$ 5.0	-	-	-	-	-	-17.7 $\pm$ 12.7
0.99 $\pm$ 0.17	$\mu$	[cP]	-17.1 $\pm$ 13.4	-	-	-18.2 $\pm$ 14.3	-28.3 $\pm$ 26.0	-23.6 $\pm$ 10.0	-	-	-	-	-	-63.0 $\pm$ 20.6
1.99 $\pm$ 0.34	$\gamma$	[-]	-8.4 $\pm$ 7.4	-	-	-8.9 $\pm$ 8.0	-13.5 $\pm$ 13.5	-11.7 $\pm$ 4.8	-	-	-	-	-	-30.8 $\pm$ 10.1
136201.26 $\pm$ 22028.42	$E_{max}$	[Pa]	-14.2 $\pm$ 21.7	12.5 $\pm$ 5.8	17.9 $\pm$ 8.9	21.2 $\pm$ 19.6	30.8 $\pm$ 30.5	-9.3 $\pm$ 6.7	-23.9 $\pm$ 10.2	-	-	-	-	-5.4 $\pm$ 3.9
35051.20 $\pm$ 5948.93	$E_{min}$	[Pa]	-716.3 $\pm$ 643.2	-19.2 $\pm$ 15.9	-118.9 $\pm$ 82.6	-796.9 $\pm$ 692.8	-1368.9 $\pm$ 1072.0	-28.1 $\pm$ 48.5	-0.5 $\pm$ 1.0	-38.7 $\pm$ 26.7	1.3 $\pm$ 1.4	3.1 $\pm$ 1.0	-	-1.7 $\pm$ 13.4

Table 5. Table 5: Analysis of the effect of lymphatic valve deficits on lymphodynamical indexes . Here we compare indexes for healthy and defective valves. In the first column, we show lymphodynamical indexes, while in the second and third we show results for the stenotic and regurgitant valve, respectively. In the last column results for healthy valves are shown. A green-coloured result indicates a normalised, percentage change in absolute compared to the healthy case value between 15 and 50 % percent \% . Likewise, a blue-coloured result indicates a change between 50 and 200 % percent \% and a red-coloured result indicates a change above 200 % percent \% . The arrows indicate a positive or a negative change.

	Stenotic valve				Regurgitant valve				Healthy valves
	Left		Right		Left		Right
	$M_{st}^L=$ 0.5	$M_{st}^L$= 0.1	$M_{st}^R=$ 0.5	$M_{st}^R=$ 0.1	$M_{rg}^L=$ 0.1	$M_{rg}^L=$ 0.8	$M_{rg}^R=$ 0.1	$M_{rg}^R=$ 0.8	$M_{st}=1$ $M_{rg}=0$
Frequency [min-1]	7.51 $\approx$	5.87 $\downarrow$	7.54 $\approx$	7.54 $\approx$	7.54 $\approx$	7.54 $\approx$	17.19 $\uparrow$	17.54 $\uparrow$	7.54
SW [nL cmH2O]	160.74 $\approx$	174.27 $\uparrow$	167.63 $\approx$	122.43 $\downarrow$	162.01 $\approx$	54.15 $\downarrow$	126.10 $\approx$	24.65 $\downarrow$	148.10
EF [-]	0.66 $\approx$	0.65 $\approx$	0.65 $\approx$	0.32 $\downarrow$	0.77 $\uparrow$	0.81 $\uparrow$	0.70 $\approx$	0.69 $\approx$	0.66
SV [nL]	46.81 $\approx$	45.30 $\approx$	45.58 $\approx$	22.07 $\downarrow$	53.78 $\uparrow$	56.66 $\uparrow$	54.97 $\uparrow$	54.76 $\uparrow$	45.71
FPF [min-1]	4.97 $\approx$	3.83 $\downarrow$	4.94 $\approx$	2.39 $\downarrow$	5.83 $\uparrow$	6.14 $\uparrow$	11.98 $\uparrow$	12.19 $\uparrow$	4.95
CPF [$\mu$L h-1]	21.09 $\approx$	15.97 $\downarrow$	20.62 $\approx$	9.99 $\downarrow$	24.34 $\uparrow$	25.64 $\uparrow$	56.70 $\uparrow$	57.65 $\uparrow$	20.68
Mean Flow [$\mu$L h-1]	19.09 $\approx$	12.90 $\downarrow$	18.04 $\approx$	9.45 $\downarrow$	11.38 $\downarrow$	1.42 $\downarrow$	4.36 $\downarrow$	-0.55 $\downarrow$	17.70
CPFI [-]	1.10 $\approx$	1.24 $\approx$	1.14 $\approx$	1.06 $\approx$	2.14 $\uparrow$	18.06 $\uparrow$	13.01 $\uparrow$	103.90 $\uparrow$	1.17
WSS [mdyne cm-2]	-75.35 $\approx$	-51.90 $\uparrow$	-72.99 $\approx$	-24.75 $\uparrow$	-35.04 $\uparrow$	-20.86 $\uparrow$	-19.41 $\uparrow$	1.23 $\uparrow$	-70.81
Peak Velocity [mm s-1]	6.56 $\approx$	6.88 $\approx$	5.74 $\downarrow$	0.94 $\downarrow$	5.29 $\downarrow$	9.49 $\uparrow$	7.39 $\approx$	6.66 $\approx$	6.77
ESD [$\mu$m]	142.36 $\approx$	143.04 $\approx$	142.86 $\approx$	200.90 $\uparrow$	115.98 $\downarrow$	104.91 $\downarrow$	142.40 $\approx$	142.84 $\approx$	142.45
EDD [$\mu$m]	244.95 $\approx$	242.72 $\approx$	243.09 $\approx$	243.10 $\approx$	243.10 $\approx$	243.11 $\approx$	258.72 $\approx$	258.63 $\approx$	243.08
ESP [cmH2O]	6.56 $\approx$	6.41 $\approx$	6.95 $\approx$	9.48 $\uparrow$	6.37 $\approx$	4.21 $\downarrow$	6.43 $\approx$	6.52 $\approx$	6.40
EDP [cmH2O]	3.00 $\approx$	2.93 $\approx$	3.00 $\approx$	3.00 $\approx$	3.00 $\approx$	3.00 $\approx$	5.93 $\uparrow$	6.00 $\uparrow$	3.00

Equations

$$\{\begin{array}{l}{\partial }_{t}A+{\partial }_{x}q=0,\\ {\partial }_{t}q+{\partial }_x\left(\alpha \frac{q^2}{A}\right)+\frac{A}{\rho }{\partial }_{x}p=-\frac{f}{\rho },\end{array}$$

$$p(x,t)=K(x,t)\psi \left(A(x,t);A_0(x)\right)+p_e(x,t),$$

$$\psi (A(x,t);A_0(x))$$

$$+C\left[{\left(\frac{A(x,t)}{A_0(x)}\right)}^z-1\right],$$

$$p_{transm}(x,t):=p(x,t)-p_e(x,t).$$

$$\displaystyle K(x,t)=\dfrac{E(x,t)}{12(1-\nu^{2})}\Bigg{(}\dfrac{h_{0}(x)}{r_{0}(x)}\Bigg{)}^{3}\;,$$

$$E(x,t)=$$

$$s(x,t)\in [0,1],$$

$$\tau (x,t)=-u(x,t)\mu \frac{\gamma +2}{r(x,t)},$$

$${\partial }_t\mathbf{Q}(x,t)+{\partial }_x\mathbf{F}(\mathbf{Q}(x,t),x,t)=\mathbf{S}(\mathbf{Q}(x,t),x,t),$$

$$\mathbf{Q}\left(x,t\right)$$

$$\mathbf{F}\left(\mathbf{Q},x,t\right)$$

$$\mathbf{S}\left(\mathbf{Q},x,t\right)$$

$$\Psi$$

$$\displaystyle\begin{aligned} =A_{0}\Bigg{(}&\dfrac{1}{m+1}\bigg{(}\dfrac{A}{A_{0}}\bigg{)}^{m+1}-\dfrac{1}{n+1}\bigg{(}\dfrac{A}{A_{0}}\bigg{)}^{n+1}\\ +C&\frac{1}{z+1}\bigg{(}\dfrac{A}{A_{0}}\bigg{)}^{z+1}\Bigg{)}\;,\end{aligned}$$

$${\partial }_{A_0}\Psi$$

$$\displaystyle\begin{aligned} =-\Bigg{(}&\dfrac{m}{m+1}\bigg{(}\dfrac{A}{A_{0}}\bigg{)}^{m+1}-\dfrac{n}{n+1}\bigg{(}\dfrac{A}{A_{0}}\bigg{)}^{n+1}\\ +C&\frac{z}{z+1}\bigg{(}\dfrac{A}{A_{0}}\bigg{)}^{z+1}\Bigg{)}\;.\end{aligned}$$

$$\mathbf{S}(\mathbf{Q},x,t)=\mathbf{S}(\mathbf{Q})=\left[\begin{array}{c}0\\ -2\left(\gamma +2\right)\pi \frac{\mu }{\rho }u\end{array}\right].$$

$${\partial }_t\mathbf{Q}+\mathbf{A}(\mathbf{Q},t){\partial }_x\mathbf{Q}=\mathbf{S}(\mathbf{Q}),$$

$$\mathbf{A}(\mathbf{Q},t)=\left[\begin{array}{cc}0& 1\\ \frac{A}{\rho }{\partial }_A\psi -u^2& 2u\end{array}\right],\mathbf{S}(\mathbf{Q})=\left[\begin{array}{c}0\\ -\frac{f}{\rho }\end{array}\right].$$

$${\lambda }_1=u-c,{\lambda }_2=u+c,$$

$$c=\sqrt{\frac{A}{\rho }{\partial }_A\psi }=\sqrt{\frac{K}{\rho }\left[m{\left(\frac{A}{A_0}\right)}^m-n{\left(\frac{A}{A_0}\right)}^n+Cz{\left(\frac{A}{A_0}\right)}^z\right]}.$$

$${\mathbf{R}}_1={\gamma }_1\left[\begin{array}{c}1\\ u-c\end{array}\right],{\mathbf{R}}_2={\gamma }_2\left[\begin{array}{c}1\\ u+c\end{array}\right],$$

$$G\left(\frac{A}{A_0}\right)$$

$$\displaystyle+Cz\left(z+2\right)\bigg{(}\dfrac{A}{A_{0}}\bigg{)}^{z}\;.$$

$$\begin{array}{l}{\lambda }_1-\text{GRI}:u+\int \frac{c\left(A\right)}{A}\mathrm{d}A=constant,\\ {\lambda }_2-\text{GRI}:u-\int \frac{c\left(A\right)}{A}\mathrm{d}A=constant.\end{array}\}$$

$$\frac{d}{dt}\mathbf{Y}=\mathbf{L}\left(\mathbf{Y}\right),$$

$$\mathbf{Y}\left(t\right)=\begin{bmatrix}v\left(t\right)\\ w\left(t\right)\\ I\left(t\right)\\ s\left(t\right)\end{bmatrix}\;,\quad\mathbf{L}\left(\mathbf{Y}\right)=\begin{bmatrix}a_{1}\big{[}v\big{(}v-a_{2}\big{)}\big{(}1-a_{3}v\big{)}-w+vI\big{]}\\ b_{1}v-b_{2}w\\ f_{I}(\bar{A},\bar{\tau},v,w,I)\\ f_{s}(v,w,s)\end{bmatrix}\;,$$

$$f_{I}(\bar{A},\bar{\tau},v,w,I)=\begin{dcases}\left(k_{Ca}^{(1)}+k_{Ca}^{(2)}\bigg{(}\dfrac{\bar{A}}{A_{Ca}}\bigg{)}^{n_{Ca}}\right)f_{NO}\big{(}\bar{\tau}\big{)}\;,&\sqrt{v^{2}+w^{2}}\leq R_{I}\;,\\ -Ik_{rel}\;,&\sqrt{v^{2}+w^{2}}>R_{I}\;,\end{dcases}\;$$

$$f_{NO}\big{(}\bar{\tau}\big{)}=1-k_{NO}\left(\dfrac{2}{1+exp\left({-\left|\dfrac{\bar{\tau}}{\tau_{NO}}\right|^{n_{NO}}}\right)}-1\right)\;,$$

$$f_s(v,w,s)=\{\begin{array}{ll}+c_{1}vw\left(1-s\right),& v>0,\\ -c_{2}s{\left|1-w\right|}_{+},& v\le 0,\end{array}$$

$${\left|x\right|}_{+}=\{\begin{array}{ll}x,& x>0,\\ 0,& x\le 0.\end{array}$$

$$\underset{|\bar{\tau }|\to +\infty }{\lim }{f}_{NO}(\tau )=1-k_{NO}\le f_{NO}\le 1=f_{NO}(0).$$

$$\bar{A}(t)=\frac{1}{L}{\int }_0^{L}A(x,t)\mathrm{d}x,\bar{\tau }(t)=\frac{1}{L}{\int }_0^L\tau (x,t)\mathrm{d}x,$$

Full text

1. 1 Introduction

2. 2 Mathematical models

3. 2.1 A one-dimensional model for lymph flow

4. 2.1.1 Conservative formulation of the one-dimensional lymph flow equations

5. 2.1.2 Mathematical analysis of the one-dimensional lymph flow equations

6. 2.2 A dynamical Electro-Fluid-Mechanical Contraction (EFMC) model

7. 2.2.1 Mathematical analysis of the modified FitzHugh-Nagumo model

8. 2.2.2 Qualitatively analysis of the EFMC model

9. 2.2.3 Analysis of the EFMC frequency

10. 2.3 A lumped-parameter model for lymphatic valves

11. 3 Numerical methods

12. 3.1 A finite volume method for the one-dimensional model

13. 3.2 Coupling between valves and lymphangions

14. 3.3 Imposed pressure at boundaries

15. 3.4 Numerical methods for the systems of ODEs

16. 3.5 Complete algorithm

17. 4 Results and discussion

18. 4.1 Test problem with piecewise initial condition: a Riemann problem

19. 4.2 Four representative test cases of a collecting lymphatic

20. 4.2.1 Test case 1: representative case of a single lymphangion

21. 4.2.2 Test case 2: contraction frequency increases as the intraluminal pressure increases

22. 4.2.3 Test case 3: contraction frequency decreases with increasing WSS

23. 4.2.4 Test case 4: representative case of a collecting lymphatic composed of ten lymphangions

24. 4.3 Pressure versus normalised cross-sectional area (PA) plots for a single lymphangion

25. 4.4 Analysis of lymphodynamical indexes by varying $P_{in}$ and $P_{out}$

26. 4.5 Sensitivity analyses

27. 4.6 A quantitative study on the lymphodynamical effect of stenotic and regurgitant lymphatic valves

28. 5 Concluding Remarks

1 Introduction

The lymphatic system is an intricate network of vessels and nodes which connect tissues to the bloodstream. The main functions of the lymphatic system comprise the hydrodynamic homoeostasis of the fluid tissues through the drainage of the excess interstitial fluid, the transport of proteins and waste products, as well as the transport of immune cells. From the pumping action of the heart, blood flows through arteries, arrives at capillaries, and is filtered into interstitial fluids through Starling forces. A portion of it, known as lymph, enters into the lymphatic system, is pushed uphill through several mechanisms, is filtered at lymph nodes, and is drained into the venous system at the junction between the internal jugular veins and the subclavian veins. In addition, according to the classical theory, cerebrospinal fluid drains into nasal lymphatics through the cribriform plate, flows into deep cervical lymph nodes, and enters into the subclavian veins. Moreover, interstitial fluid from the central nervous system parenchyma drains into lymph nodes within the basement membrane of cerebral capillaries and arteries and then enters into the subclavian veins [1, 2, 3, 4, 5]. From there on, lymph mixes with blood again, enters into the heart and the cycle restarts. During these regular cycles, leukocytes travel throughout the lymphatic system and provide an immune response for the body.

Several pathologies are associated with the impairment of the lymphatic system. Very recently, Louveau et al. [6] and Aspelund et al. [7] have just discovered the lymphatic system in the brain, known as meningeal or dural lymphatics. This discovery has significant effects on the immune privilege concept of brain, on brain-waste clearance of parenchyma, on the genesis of neurological disorders [6, 8, 9, 10, 11, 2], and on possibly new therapeutic treatments for removal of macromolecules in the central nervous system [12].

The building block of the lymphatic system is the lymphangion: a mini-heart like, compliant vessel, which contracts and pushes lymph centripetally, and has several mechanobiological auto-regulatory systems to provide optimal flows in various scenarios [13, 14]. It is enclosed between valves which prevent reverse flow and allow for unidirectional flows. The contraction of lymphangions has been described for several animals in [15, 16, 17]; it is regulated by Ca2+ fluxes from extracellular and intracellular stores [18], is divided in systolic and diastolic phases [16], and can be described using cardiac indexes [19, 20, 21, 22]. Several works have been devoted to studying the membrane potential in lymphangions [23, 24, 25].

The mechanobiological systems which underlie the auto-regulatory homeostatic functions of lymphatic contractions, are still being investigated [26, 27, 13]. The frequency of the lymphatic contractions depends on the rate of stretch of the vessel wall. Indeed, stretch-activated calcium channels in the vessel wall change conformation in response to membrane tension [13, 23]. McHale and Roddie [28] described the dynamics of frequencies in bovine mesenteric vessels varying intraluminal pressures and showed that the more the lymphangions are stretched, the faster they contract. Lymphatic muscle also exhibits rate-sensitive contractile responses to stretch, as described for rat mesenteric lymphatics by Davis et al. [29]. The authors analysed the responses in amplitude and frequency to time-varying preload and pressure. Bursts of contraction occurred when positive ramps were imposed. Lymphatic contractions are also inhibited by flows. Gashev et al. [30] studied rat mesenteric lymphatics in response to imposed flows and showed that the frequency dropped from $9.0\pm 1.6$ min*-1* to $3.1\pm 1.4$ min*-1* when flow changed from zero to a transaxial-pressure-gradient induced flow of $7$ cmH2O. In addition, the activity of the lymphatic contractions differs from region to region. Gashev et al. [31] showed that the flow-induced inhibition of contraction in the thoracic duct of rats is more evident when compared to that of femoral lymphatics. As a matter of fact, the authors showed that the frequencies of thoracic duct and femoral lymphatics are $4.6\pm 0.6$ min*-1* and $15.2\pm 2.6$ min*-1*, respectively, at zero flow conditions, and $0.1\pm 0.1$ min*-1* and $5.0\pm 2.4$ min*-1*, respectively, at transaxial-pressure-gradient induced flow of $5$ cmH2O. An underlying mechanism of the above-mentioned flow-induced contraction inhibition is the local production of NO induced by wall shear forces [14, 32, 33], which results in blunting the Ca2+-dependent contraction. A lymphatic vessel composed of two or more lymphangions is called a collecting lymphatic. Several studies have analysed the coordination of adjacent lymphangions [34, 15, 35], wave-propagations throughout the collector [36, 35, 37] and the microstructure of the vessel wall [36, 38, 39]. Overall, the above-mentioned processes which underlie and permit optimal lymph flows, are more complex and difficult to quantify, when compared to those of arteries and veins. For complete reviews of the mechanics of lymphangions and collectors, see [13, 26, 27, 40, 41, 42, 43].

The lymphatic system has two different types of valves called primary and secondary valves. The former is located at the initial lymphatics at the level of the endothelium, while the latter is located between lymphangions in collectors [44, 45]. Primary and secondary lymphoedema, a lymphatic disease that leads to tissue swelling, is linked to lymphatic valve deficits [46, 47, 48, 49, 50]. For instance, the lack of valves in the lymphoedema distichiasis impairs lymphatic flow due to the inability to properly pump lymph forward [47, 51, 52, 45]. Also, the chronic venous insufficiency leads to fibrotic lymph vessels due to hypertension, then it compromises the functionality of lymphatic valves, and finally result in accumulation of fluid in tissues [53, 54]. In addition, in secondary lymphoedema after lymph node dissection, lymph retention and lymphatic hypertension occur, and valvular dysfunction induces retrograde lymph flow [50]. Despite the connections of lymphatic valve deficits and lymphoedema, to the authors’ knowledge, the effect of stenotic and regurgitant valves in the lymphatic system has not properly been investigated and quantified. This is probably due to the difficulties in performing experiments on animal lymphatic valves, though the effects of genes mutations in engineered mice can be studied.

A huge gap is currently present in the literature between mathematical models for the circulatory [55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66], and lymphatic systems. The first mathematical attempt to model for the lymphatic system is attributed to Reddy et al. [67], and was based on one-dimensional flow equations. The model was then refined by MacDonald et al. [68] including a tension and a damping term. However, the convection term, tube law and valve model were relatively simple. Extensive work has been done in lumped-parameter models [69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81]. Jamalian et al. [80] constructed a lumped-parameter model able to simulate lymph transport in a network of rat lymphangions based on experimental measurements for vessel and valve parameters [75]. The authors analysed the effect of pumping coordination in branched network structures. Caulk et al. [39] described a detailed microstructurally four-fibre family constitutive law for rat thoracic duct based on experimental measurements of collagen, elastin and solid structures of the lymphatic wall. Then, the authors combined the lumped-parameter model described by Bertram et al. [77] with their four-fibre family constitutive law, and studied the variation of muscle contractility in response to a sustained elevation in afterload [81]. Rahbar et al. [82] investigated the validity of assuming Poiseuille flow to estimate the wall shear stress. A mechanobiological oscillatory model for the lymphatic contraction has been proposed by Kunert et al. [14]. Compared to the other existing mathematical models for lymphangions, the authors provided a biologically-based, dynamical model for the contractibility of the vessel wall. They proposed 1) evolutionary equations for Ca2+ and NO, and 2) a coupling between wall mechanical forces, Ca2+ and NO. The resulting model was able to control lymphatic transport via mechanobiological feedback loops, given by stretch-activated contractions and flow-induced relaxations. Very recently, Baish et al. [83] proposed a model of a vascular oscillator and studied the interaction of Ca2+ and NO in the context of lymphatic vessel contractions.

With the aim of constructing a mathematical model of the entire lymphatic system, the above-mentioned mathematical models for lymphangions, except for [14, 83], are based on a relatively simple contraction dynamic. As a matter of fact, these models 1) prescribe muscle contractility dynamics by using trigonometric functions, and 2) artificially prescribe time delays between adjacent lymphangions. This result in non-dynamical models insofar as contractions are triggered by artificial parameters rather than by local dynamical factors, and contraction frequencies in those models do not depend on the rate of stretch of the lymphatic wall, nor on local shear forces.

In the study of cardiac contraction, there exist an extensive literature. From the pioneer, detailed model of Hodgkin and Huxley in 1952 [84] of how action potential in neurons are initiated and propagated, several extended and even simplified models have been proposed in the literature for neurons and heart contractions [85, 86, 87, 88]. The FitzHugh-Nagumo model [86] is an example of a significantly simplified, but elegant two-parameter formulation of the original Hodgkin-Huxley model [89]. It is a set of two Ordinary Differential Equations (ODEs) with a fast and a slow variable. The former represents the action potential, while the latter phenomenologically summarises all the effects of all ionic currents. The FitzHugh-Nagumo model has been adapted to simulate rhythmically discharging cells in the heart membrane, such as the atrioventricular and sinoatrial node [90]. Thanks to its simplicity, it also permits a rigorous analysis of the stability. Action potential arises to be triggered by an unstable equilibrium state of the ODEs. Many studies have been done to couple modified versions of the FitzHugh-Nagumo model to the heart contractions [90]. However, to date no studies have been performed to model contractions of lymphangions with the over-mentioned dynamical and phenomenological set of ODEs for action potentials.

In the present paper, we propose a one-dimensional model for collecting lymphatics coupled with an Electro-Fluid-Mechanical Contraction (EFMC) model for dynamical lymphatic contractions based on a modified FitzHugh-Nagumo model. The one-dimensional model for a compliant lymphatic vessel adopt a tube law which fits the experimental measurements shown in Bertram et al. [75] and performed by Davis et al. [19]. The resulting system is a set of hyperbolic Partial Differential Equations (PDEs), is written in conservative formulation, allows for space-variable geometrical parameters and permits simulating contraction-wave propagations based on a prescribed space-time variation of the Young modulus. The EFMC model is a set of four ODEs. We carried out a complete mathematical analysis of the EFMC model on the stability of its equilibrium state. Contractions turn out to be triggered by the instability of the stationary state and phenomenologically depend on: (1) environmental calcium influx, (2) stretch-activated calcium influx, and (3) contraction inhibitions induced by wall shear stresses. Lymphatic valves are modelled using a well-established, versatile lumped-parameter model proposed by Mynard et al. [91] and analysed by Toro et al. [92]. This valve model also allows simulating stenotic and regurgitant valves and quantify the lymphodynamical effects. Here, we simulate lymphatic vessel from the mesentery of rats. We use the SLIC method [93] to numerically solve the one-dimensional equations and an implicit Runge-Kutta solver for the systems of ODEs. We show some numerical tests of our mathematical model inspired by experiments [19, 20, 21, 22]. Then we analyse several lymphodynamical indexes of a single lymphangion for different combinations of upstream and downstream pressures and considering both positive and negative pressure gradients. We performed two sensitivity analyses based on [94, 59] for both positive and negative pressure gradients and investigate the most influential parameters affecting the lymphodynamical indexes. Finally, we quantified the effects of stenotic and regurgitant valves.

The rest of this paper is structured as follows: in Section 2 we propose a one-dimensional model for lymph flow that allows for variable parameters, we describe the EFMC model, analyse it and then we briefly review the valve model. In Section 3 we briefly review the finite volume schemes used for the one-dimensional lymph flow equations, explain how to couple valves with lymphangions, how to impose boundary pressure at the interfaces of the lymphangions, and describe the numerical methods used for the valve and EFMC models. In Section 4 we show some numerical results our mathematical model, analyse the lymphodynamical indexes for a wide range of boundary pressures, perform two sensitivity analyses for positive and negative pressure gradients and finally investigate the effects of defective lymphatic valves, namely stenotic and regurgitant valves. Section 5 gives a summary and concluding remarks.

2 Mathematical models

A lymphangion is a lymphatic vessel delimited by two valves, also called secondary lymphatic valves, located on the upstream and downstream boundaries of the lymphangion. A chain of lymphangions is called collecting lymphatic. Here we aim to model the dynamics of flowing lymph inside a collecting lymphatic propelled by lymphatic contractions and pressure gradients, and the dynamic of lymphatic valves. First, we present a one-dimensional mathematical model for lymph flow and a set of four ODEs to model the lymphatic contraction, and then we describe a mathematical model for lymphatic valves. Fig. 1 illustrates a collecting lymphatic, a single lymphangion and two lymphatic valves. It also shows a generic cross-sectional area at a given position $x$ and time $t$.

2.1 A one-dimensional model for lymph flow

Here we assume the lymph to be an incompressible Newtonian fluid. To derive the one-dimensional flow equation for a compliant lymphatic vessel, one can follow the procedure done for arteries and veins, where Reynold’s transport theorem is used to obtain the conservation of mass and momentum in a flexible tube, see [95, 96, 97]. The one-dimensional flow equations for a compliant lymphatic vessel are

[TABLE]

where $x$ is the space variable, $t$ is time, $\alpha$ is the Coriolis coefficient assumed to be $\alpha=1$, $A(x,t)$ is cross-sectional area of the vessel, $q(x,t)=A(x,t)u(x,t)$ is flow, $u(x,t)$ is velocity, $p(x,t)$ is pressure, $\rho$ is lymph density, $f(x,t)=2(\gamma+2)\pi\mu u(x,t)$ is friction force per unit length of the tube with parameter $\gamma$ chosen from the velocity profile [98], and $\mu$ is the dynamic viscosity. There are two governing partial differential equations and three unknowns, namely $A(x,t)$, $q(x,t)$ and $p(x,t)$. For this reason, an extra relation is required to close the system, the tube law, which relates pressure $p(x,t)$ and cross-sectional area $A(x,t)$. Equations (1) have been widely described in the literature for blood flow [97, 99, 95, 96], where different tube laws were proposed depending on the characteristics of the wall material. The lymphatic wall is characterized by elastin, collagen, smooth muscle cells and other extracellular matrix constituents [39, 38]. Elastin fibres give to lymphatic vessels a compliant, elastic behaviour, while collagen prevents vessels from stretching beyond their physiological limits. The overall dynamics of elastin and collagen is reflected in highly non-linear tube laws. Several works have been proposed to describe tube laws for rats [39, 38, 77], bovine [68] and human thoracic ducts [100]. Here we propose the following general and purely elastic tube law:

[TABLE]

with

[TABLE]

where $p_{e}(x,t)$ is the external pressure, $A_{0}(x)$ is vessel cross-sectional area at equilibrium, $K(x,t)$ is the stiffness of the vessel wall, $m\geq 0$, $n\leq 0$, $z\geq 0$, and $C\geq 0$ are real numbers to be specified. The transmural pressure is defined as

[TABLE]

To the authors’ knowledge, a mechanical study of the stiffness of collecting lymphatic is not present in the literature. It is clear that lymphatics are highly deformable and their physiological pressure is comparable or smaller than that of veins. As a matter of fact, the lymphatic pressure approximately ranges from the interstitial fluid pressure to the subclavian vein pressure. Here we assume parameter $K$ of lymphatics to be equal to that of veins [60], namely

[TABLE]

where $E(x,t)$, $\nu$, $h_{0}(x)$ and $r_{0}(x)$ are respectively the Young modulus of elasticity, the Poisson ratio, the wall-thickness and cross-sectional radius at zero transmural pressure (equilibrium), respectively. The assumption to use parameter $K$ equal to that of veins, might need further investigations. Following [68, 101], lymphatic contractions are modelled by varying the Young modulus $E(x,t)$ from a minimal value $E_{min}(x)$ to a maximum value $E_{max}(x)$ as follows

[TABLE]

where $s(x,t)$ is the state of contraction. The lymphangion can be relaxed and contracted and this is described by variable $s(x,t)$. The lymphangion is contracted when $s(x,t)=1$ and the Young modulus reaches its maximum $E(x,t)=E_{max}$ while the lymphangion is relaxed when $s(x,t)=0$ and the Young modulus reaches its minimum $E(x,t)=E_{min}$. In Section 2.2 we provide the governing equation for $s$ with a set of four ODEs.

In the present paper, we simulated lymphatic vessel from the mesentery of rats, whose parameters are found in Table 1. The tube law used in Eq. (3) is a relationship between normalised cross-sectional area and pressure. The parameters of the tube law and the minimum and maximum Young’s modulus were tuned to fit the experimental measurements shown in Bertram et al. [75] and performed by Davis et al. [19]. The experimental measurements and the tube law with and without contraction are shown in Fig. 2.

The Wall Shear Stress (WSS) is fundamental in connecting the fluid mechanics and the Nitric Oxide (NO) production. The fluid movement increases the shear at the lymphatic wall and increases the production of NO through the activation of local endothelium NOS. The production of NO results in vasodilating the lymphatic vessel and in decreasing the contraction frequency by blunting the $Ca^{2+}$-dependent contraction. For references, see [14, 32, 33, 30, 26]. According to [98], the WSS in our formulation is

[TABLE]

where $r(x,t)$ is cross-sectional radius.

2.1.1 Conservative formulation of the one-dimensional lymph flow equations

It is possible to write the lymph flow equations in conservative form as follows:

[TABLE]

where

[TABLE]

with

[TABLE]

and

[TABLE]

The constants arising from the integrals (12) and (13) are set to zero for consistency with (1) and (2), see [102, 103, 97].

The present formulation allows for space-time variable Young modulus. This let us simulate travelling contraction-waves through the lymphatic wall by prescribing a space-time varying trigonometric contraction state $s(x,t)$. However, in the present work we consider the simpler case in which the contraction state is constant throughout the lymphangion, namely $s=s(t)$, and we also neglect the interaction between adjacent lymphangions. Then, instead of prescribing a trigonometric function for $s$, here we propose a set of governing ODEs given in Section 2.2.

Here, parameters $h_{0}(x)$, $r_{0}(x)$, $E_{min}(x)$, $E_{max}(x)$ and $p_{e}(x)$ are constant in space. As a result, the source term simplifies in

[TABLE]

The general case of variable material properties poses numerical and mathematical [104, 105] challenges and require the use of well-balanced schemes, see [106, 107].

2.1.2 Mathematical analysis of the one-dimensional lymph flow equations

Here we study the mathematical properties of (8) assuming constant parameters along the lymphatic vessel. The equations in (8) are a generalization of the one-dimensional blood flow equations presented in [95]. As a matter of fact, the main difference is an additional term in the tube law (3). For this reason, here we summarize the main mathematical structure of the lymph flow equations without proofs. System (8) can be written in quasi-linear form as

[TABLE]

where

[TABLE]

The eigenvalues of matrix $\mathbf{A}$ are

[TABLE]

where $c$ is the wave speed

[TABLE]

We assume parameters $m\geq 0$, $n\leq 0$, $z\geq 0$, and $C\geq 0$ for the tube law. Thus, the wave speed $c$ is always real. The eigenvectors of $\mathbf{A}$ are

[TABLE]

where $\gamma_{1}$ and $\gamma_{2}$ are arbitrary scaling factors. It can be shown that system (8) is hyperbolic, as the eigenvalues are real and distinct and the eigenvectors $\mathbf{R}_{1}$ and $\mathbf{R}_{2}$ are linearly independent. Following proofs in [95] and [97], the $\lambda_{1}$ and $\lambda_{2}$ characteristic fields are genuinely non-linear outside the locus of the following function

[TABLE]

With the choice of parameters $m$, $n$ $z$ and $C$ in Table 1, there exist at least one solution of $G\left(\frac{A}{A_{0}}\right)=0$. This means that the two characteristic fields are neither genuinely non-linear nor linearly degenerate. The consequences of this are unclear to the authors, and might require further investigations. See [108] and [109] for details. The Generalized Riemann Invariants (GRIs) for $\lambda_{1}$ and $\lambda_{2}$ characteristic fields are respectively

[TABLE]

In the present work, the generalized Riemann invariants will be used to couple valves with lymphangions and to impose the pressure at the terminal interfaces of the collector.

2.2 A dynamical Electro-Fluid-Mechanical Contraction

(EFMC) model

Lymphatic contractions are preceded by action potentials [41, 24, 26]. Action potentials occur at smooth muscle cells of the collecting lymphatic wall. These are triggered by the rapid influx of $Ca^{2+}$ and other ions through different type of channels (low voltage T-channels, high-voltage L-channels, stretch-activated channels, gap junctions) [13]. From the pioneering work of Hodgkin and Huxley in 1952 [84] on action potentials in neurons, mathematical models for action potentials have been widely used in the literature and applied in several contexts [85, 86, 87, 88]. Here we propose an Electro-Fluid-Mechanical Contraction (EFMC) model for lymphatics, based on the FitzHugh-Nagumo model for action potentials. This model combines the electrical activity of lymphangions (action potential) with fluid-mechanical feedback (stretch of the lymphatic wall and wall shear stress), and the mechanical variation of the lymphatic wall properties (contractions).

The modelling system of ODEs is

[TABLE]

where

[TABLE]

and

[TABLE]

[TABLE]

and

[TABLE]

where $\left|x\right|_{+}$ is the positive value of $x$

[TABLE]

The unknowns of the above system are: the excitable variable $v(t)$ (membrane potential), the recovery variable $w(t)$, the stimulus $I(t)$, and the contraction state $s(t)$ introduced in Eq. (6). The first two equations of (23) are the classical FitzHugh-Nagumo (FHN) model with a simple modification. In the classical formulation of the FHN model, the stimulus $I$ has a constant value, while here it varies in time and multiplies the excitable variable $v$. As described in the next section, this modification allows us to control the equilibrium of the stationary solution and to trigger action potentials under certain conditions on the stimulus $I$.

The evolution in time of $I$ is controlled by the function $f_{I}$ and phenomenologically depends on three mechanisms: (1) environmental calcium influx, (2) stretch-activated calcium influx, and (3) contraction inhibitions induced by WSS. The environmental baseline influx is regulated by the parameter $k_{Ca}^{(1)}$. The stretch-activated calcium influx is regulated by the parameters $k_{Ca}^{(2)}$, $A_{Ca}$ and $n_{Ca}$. The contraction inhibitions induced by WSS are regulated by the function $f_{NO}$, which depends on parameters $k_{NO}$, $\tau_{NO}$ and $n_{NO}$. The function $f_{NO}$ is bounded by $1-k_{NO}$ and $1$, namely

[TABLE]

The contraction state is controlled by the function $f_{s}$, which depends on parameters $c_{1}$ and $c_{2}$. Functions $f_{I}$ and $f_{NO}$ are evaluated at the space-averaged cross-sectional area and at the space-averaged WSS, respectively, at the current time

[TABLE]

where $L$ is the length of the lymphangion. Parameters of the EFMC model are found in Table 1.

2.2.1 Mathematical analysis of the modified FitzHugh-Nagumo model

Here we analyse the modified FitzHugh-Nagumo model on which the EFMC model is based. First we find the stationary state solution, and then we study its nature depending on the stimulus $I$.

The stationary points are found by solving the following system

[TABLE]

From now on, we assume $b_{2}=0$. Therefore, the only stationary point of (30) is $(v,w)=(0,0)$ and is independent of $I$. To study the nature of the stationary point, one has to evaluate the Jacobian of $F_{FHN}$ at the stationary state and study the sign of its eigenvalues. The resulting eigenvalues are

[TABLE]

Then we conclude that the stationary solution $(0,0)$ is a

[TABLE]

When $I<a_{2}$, the nature of the stationary solution is stable, and therefore action potentials are not automatically triggered. When $I>a_{2}$ the nature of the stationary solution is unstable, and therefore action potentials can be triggered. Fig. 3 summarizes the study of the nature of the stationary solution.

We stress that the mathematical analysis reported here is based on well-known mathematical criteria and is valid only near the stationary solution. If $I>a_{2}$ and the excitable and recovery variables are not near enough to the stationary solution $(0,0)$, then we might not have an action potential. In other words, when $I$ exceeds the threshold level $a_{2}$, the stationary solution becomes unstable, and an action potential develops as soon as the variables $v$ and $w$ are near enough to the stationary solution. For a time-varying stimulus $I$, we assume that the needed stimulus $\tilde{I}$ to trigger an action potential is between a minimum and a maximum value $\tilde{I}_{min}$ and $\tilde{I}_{max}$, defined as follows:

[TABLE]

These two values will be useful to estimate the frequency of contractions of the EFMC model.

Note that the initial condition for the excitable variable $v$ or the recovery variable $w$ needs to be different from zero. Otherwise, action potentials do not occur in time.

2.2.2 Qualitatively analysis of the EFMC model

Here we qualitatively explain the EFMC model. Fig. 4 depicts a representative numerical result of the EFMC model. Excitable variable $v$, recovery variable $w$, contraction state $s$ and stimulus $I$ are shown for two lymphatic cycles. When the excitable variable $v$ and the recovery variable $w$ are near the stationary state ($\sqrt{v^{2}+w^{2}}<R_{I}$), the stimulus $I$ increases as given by Eq. (24). When a certain value $\tilde{I}$ is reached (in this case $\tilde{I}_{min}<\tilde{I}<\tilde{I}_{max}$), an action potential is triggered: variables $v$ and $w$ perform a cycle, increase in absolute value and move far from the stationary state ($\sqrt{v^{2}+w^{2}}>R_{I}$) and consequently the stimulus $I$ exponentially decreases to zero. The state of contraction $s$ increases until $v>0$, and then decreases to zero, see Eq. (26). When the action potential ends, variables $v$ and $w$ return to the equilibrium point. From there on, the stimulus $I$ restarts to increase and possibly triggers a new contraction.

Fig. 5 shows the numerical result in the phase space. Three nullclines are represented. They are: the nullcline for the recovery variable $\partial_{t}w=0:v=0$, and two representative nullclines with $I=0$ and $I=\tilde{I}$ for the excitable variable \partial_{t}v=0:w=v\big{(}v-a_{2}\big{)}\big{(}1-a_{3}v\big{)} and \partial_{t}v=0:w=v\big{(}v-a_{2}\big{)}\big{(}1-a_{3}v\big{)}+v\tilde{I}. As soon as a contraction occurs, the stimulus $I$ decreases and the third nullcline tends on the second one. The sphere of radius $R_{I}$ centred at $v=w=0$, intersection of all the nullclines, divides the phase space into two regions. We call the region outside the sphere the excited region, while the latter the activation region. In the excited region, the solution quickly performs a cycle, while the stimulus $I$ exponentially decreases to zero. In the activation region, the numerical solution tends to the equilibrium state, while the stimulus $I$ increases.

The time required for the numerical solution to perform a cycle (from the excited region into the activation region) can be numerically evaluated and is denoted with $t_{excited}$. The time required to activate an action potential is denoted with $t_{activation}$.

2.2.3 Analysis of the EFMC frequency

Lymphangions contract differently according to the location of the vessel, the stretch of the lymphatic wall and wall shear stress feedback [31]. Here we aim to analyse the frequency of contraction of the EFMC model and to estimate parameters $k_{Ca}^{(1)}$ and $k_{Ca}^{(2)}$ in order to prescribe a baseline frequency $f_{min}$ and a frequency $f_{Ca}$ at $\bar{A}=A_{Ca}$. The time $t_{total}$ between two cycles can be written as follows:

[TABLE]

and the relative frequency is

[TABLE]

The excited time $t_{excited}$ can be assumed to be constant and can be evaluated numerically from the FHN model. The activation time, on the other hand, depends strongly on the rate of increase of $I$ given by (24). We now estimate the activation time, namely the time required for the stimulus $I$ to attain a certain triggering value $\tilde{I}$. Near the stationary solution $v=w=0$, it is reasonable to assume $\sqrt{v^{2}+w^{2}}<R_{I}$. We solve the following initial value problem

[TABLE]

whose solution is

[TABLE]

We assume that during the activation time, $\bar{A}$ and $\bar{\tau}$ are constant in time because the lymphangion is already at the end of the diastolic phase. Thus, the above integral can be exactly computed and is

[TABLE]

Consequently, the activation time $t_{activation}$ required to attain a triggering value $\tilde{I}$ is

[TABLE]

The maximum activation time (when $\textstyle\left({\bar{A}}/{A_{Ca}}\right)^{n_{Ca}}\approx 0$) and the activation time at $\bar{A}=A_{Ca}$ both at zero WSS ($\bar{\tau}=0$) are

[TABLE]

The maximum activation time, corresponding to the minimum frequency, depends on parameter $k_{Ca}^{(1)}$. In our model, parameter $k_{Ca}^{(1)}$ phenomenologically represents the environmental calcium influxes. Indeed, in the limiting case in which there is no environmental calcium influxes ($k_{Ca}^{(1)}=0$), the activation time becomes infinite, which means that the lymphangion does not autonomously contract. Parameter $k_{Ca}^{(2)}$, on the other hand, phenomenologically regulates the stretch-induced calcium influxes. The greater the parameter $k_{Ca}^{(2)}$, the less the activation time, and thus the greater the influence of the stretch of the lymphatic wall on the frequency contractions.

If we assume the frequencies $f_{min}$ and $f_{Ca}$, corresponding to $t^{max}_{activation}$ and $t^{Ca}_{activation}$ respectively, to be known, then we have

[TABLE]

Using (39) and (40), we can explicitly find parameters $k_{Ca}^{(1)}$ and $k_{Ca}^{(2)}$:

[TABLE]

To assure positive activation times $t^{max}_{activation}>0$, $t^{Ca}_{activation}>0$, the following conditions need to be satisfied:

[TABLE]

Here we assume the triggering value $\tilde{I}$ to be the mean value of $I_{max}$ and $I_{min}$ defined in Eq. (32), namely

[TABLE]

Numerical results confirmed that this is a good choice, even though $\tilde{I}_{min}$ and $\tilde{I}_{max}$ can be used as triggering values too. Substituting $k_{Ca}^{(1)}$ and $k_{Ca}^{(2)}$ and the activation time $t_{activation}$ defined in (38) into (34), one obtain a frequency function as

[TABLE]

Then, one can easily prove the following inequalities

[TABLE]

and

[TABLE]

The first property (45) says that the frequency decreases as the absolute value of the WSS increases, and maximum contractions are attained at zero WSS. The second property (46) says that the frequency increases as the cross-sectional area increases, that is the higher the intraluminal pressure, the higher the frequency of the contraction.

Fig. 6 shows theoretical results of the EFMC model. The parameters were taken from Table 1. The top row shows pressure against frequency, for different values of $f_{Ca}$, $f_{min}$ and $n_{Ca}$, while the bottom row shows WSS against frequency for a given pressure, for different values of $n_{NO}$, $\tau_{NO}$ and $k_{NO}$. In the top row, we considered a linear variation for pressure between 0 and 12 cmH2O and obtained the corresponding cross-sectional area $A$ using the inverse of the tube law (2). Then we depicted a coloured shaded area for the frequency in Eq. (44) using $\tilde{I}_{min}$ and $\tilde{I}_{max}$ as triggering values. Indeed, note that

[TABLE]

Numerical results of the complete model composed of the one-dimensional lymph flow equations in Eq. (8) and of the EFMC model in Eq. (22) shows that the resulting frequency follows the theoretical trend shown in Fig. 6. These results show that it is possible to imitate the experimental measurements of a specific lymphangion by fitting parameters $f_{min}$, $f_{Ca}$ and $n_{Ca}$. In the bottom row, we considered a linear variation of flow $q$ between 0 and 240 $\mu$L min*-1* and then obtained the WSS from Eq. (7) assuming $\bar{A}=7A_{0}$. These results show that given the experimental trend of the frequency of a specific lymphangion for different flows, it is possible to adjust parameters $n_{NO}$, $\tau_{NO}$ and $k_{NO}$ to imitate that trend. For instance, from experimental measurements we know that the WSS has a greater negative chronotropic effect on the thoracic duct than on femoral lymphatics [31], and this can be modelled by using $k_{NO}\approx 1$ for the thoracic duct and, for example, $k_{NO}\approx 0.5$ for femoral lymphatics.

2.3 A lumped-parameter model for lymphatic valves

Here we follow the mathematical model proposed by Mynard et al. [91], which has been already adapted to simulate human venous valves by [61, 110, 92]. The time variation of the flow across the valve $q_{v}\left(t\right)$ is approximated as

[TABLE]

where

[TABLE]

Here $p_{u}$ and $p_{d}$ are the upstream and downstream pressures, respectively. Coefficients $B$, $L$ and $R$ are the Bernoulli resistance, the lymphatic inertia and the viscous resistance to flow, given respectively as

[TABLE]

where $L_{eff}$ is the effective length and $A_{eff}$ is the effective area, which varies from a minimum value to a maximum value as

[TABLE]

Compared to the work of Mynard et al. [91], we have added the Poiseuille-type viscous losses insofar as numerical experiments showed that this term plays an important role in the lymphatic context. See [111] for the derivation of $L$ and $R$. The minimum and the maximum effective areas are evaluated as follow

[TABLE]

where $M_{rg}$ is a parameter that controls the minimum closure. A normal closure is modelled by setting $M_{rg}=0$, while to model a regurgitant valve we use $0<M_{rg}\leq 1$. The parameter $M_{st}$ controls the maximum opening. A normal opening is modelled by setting $M_{st}=1$, while to model a stenotic valve we use $0\leq M_{st}<0$. Then $A_{0}$ is here taken as the mean value between the cross-sectional areas at equilibrium of the adjacent lymphangions. The valve state $\xi\left(t\right)$ is governed by the following ODE

[TABLE]

where $K_{vo}$ and $K_{vc}$ are the valve opening/closure rates, and $\Delta p_{vo}$ and $\Delta p_{vc}$ are the opening/closure threshold pressures. Davis et al. [19] noticed that lymphatic valves are biased to stay open and display hysteresis. The authors observed that the threshold pressure to open and close the valve are strictly related to the transmural pressure of the lymphatic vessel. Bertram et al. [77] proposed formulas for $\Delta p_{open}$ and $\Delta p_{close}$ to imitate those experimental measurements. In the present work, we assume both the opening and closure threshold pressures to be zero. Hence, the system of ODEs to be solved is

[TABLE]

where

[TABLE]

3 Numerical methods

Here we briefly describe the finite volume schemes used for the one-dimensional lymph flow equations, explain how lymphatic valves and lymphangions are coupled, and illustrate the treatment of the boundary conditions at the terminal interfaces of the lymphangion. Then, we summarize the numerical methods used for the valves and the EFMC models.

3.1 A finite volume method for the one-dimensional model

Consider the system of $m$ hyperbolic balance laws

[TABLE]

By integrating (56) over the control volume $V=[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]\times[t^{n},t^{n+1}]$ we obtain the exact formula

[TABLE]

with definitions

[TABLE]

[TABLE]

Eq. (58) gives the spatial-integral average at time $t=t^{n}$ of the conserved variable $\mathbf{Q}$ while Eqs. (59) give the time-integral average at interface $x=x_{i+\frac{1}{2}}$ of the physical flux $\mathbf{F}$ and the volume-integral average in $V$ of the source term $\mathbf{S}$. Spatial mesh size and time step are $\Delta x=x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}}$ and $\Delta t=t^{n+1}-t^{n}$ respectively. Finite volume methods for (56) depart from (57) to (59), where integrals are approximated, and then formula (57) becomes a finite volume method, where the approximated integrals in (59) are called numerical flux and numerical source, respectively. Here index $i$ runs from $1$ to $M$, where the cell $i=1$ is the leftmost cell with $x_{\frac{1}{2}}$ being the first interface, and the cell $i=M$ is the rightmost cell with $x_{M+\frac{1}{2}}$ being the last interface. See Fig. 7 for an illustration of the finite volume framework. To compute the time step $\Delta t$, the Courant-Friedrichs-Lewy condition is applied for each lymphangion

[TABLE]

with $CFL=0.9$. Superindex $j$ indicates the $j$-th lymphangion. Then, the time step $\Delta t$ to be used is the minimum of all the time steps, namely $\Delta t=\min\limits_{j}\left(\Delta t^{j}\right)$.

In the present paper we used the SLIC method to evaluate the numerical fluxes within the domain ($\mathbf{F}_{\frac{3}{2}},\dots,\mathbf{F}_{M-\frac{1}{2}}$) [93]. This method is second-order accurate in space and time and is based on the MUSCL-Hancock scheme where the Godunov upwind flux is replaced by the FORCE flux, see Section 14.5.3 of [112] and references therein. The numerical source was approximated using a second order in space and time method, see Chapter 19 of [112]. For the numerical fluxes at the boundaries ($\mathbf{F}_{\frac{1}{2}}$ and $\mathbf{F}_{M+\frac{1}{2}}$) we used a first-order Godunov-type method based on the solution of a classical Riemann problem at the interface.

3.2 Coupling between valves and lymphangions

Here we aim to couple valves and lymphangions. For each lymphangion, we need to calculate the numerical flux at the interface in which the valve is located, which can be either $\mathbf{F}_{\frac{1}{2}}$ or $\mathbf{F}_{M+\frac{1}{2}}$ according to Fig. 7. There are three possible configurations for a lymphatic valve. It can be the leftmost or rightmost valve of a collector, or it can be interposed between two lymphangions. In every case, the flow across the lymphatic valve is calculated from (54), where the pressure gradient $\Delta p$ in (48) is evaluated at the current time $t^{n}$ using either the two lymphangions, or one of the lymphangions and a prescribed time-varying pressure. Specifically, at $t=t^{n}$ the pressure gradient $\Delta p(t^{n})$ is

[TABLE]

where values $p_{u}(t^{n})$ and $p_{d}(t^{n})$ are

[TABLE]

and

[TABLE]

where pressures $p_{M}^{n}$ and $p_{1}^{n}$ refers to the upstream and downstream lymphangions, respectively, and $P_{in}$ and $P_{out}$ are prescribed functions of time. The three possible configurations are summarized here

[TABLE]

Once we numerically solve system (54), the flow across the valve at the future time $q_{v}^{n+1}$ is determined.

In the present paper, to find $A^{*}$ and calculate the numerical flux at the boundary we follow the numerical methodology proposed by Alastruey et al. [111]. This method has already been used in [60, 61, 113, 55]. To impose the valve flow at the boundaries of a lymphatic vessel, we use the Riemann invariant

[TABLE]

where $A^{n}$ and $u^{n}$ are the cross-sectional area and the velocity at the cell adjacent to the boundary at current time $t=t^{n}$, $q_{v}^{n+1}$ is the known flow rate across the valve and

[TABLE]

To find $A^{*}$, we solve the following non-linear algebraic equation

[TABLE]

using the Newton-Raphson iterative method. Then the numerical flux at the boundary is

[TABLE]

where

[TABLE]

When a valve is interposed between two lymphangions, then the non-linear problem (67) has to be solved twice: one for the upstream lymphangion ($\beta=1$) and one for the downstream lymphangion ($\beta=-1$). The above numerical flux is a Godunov-type flux. Indeed, the numerical flux is evaluated at the solution along the $t$-axis of a classical Riemann problem where one of the two unknowns is given by an ODE, which in this case is the valve model, see [114] for similar works.

Fig. 8 illustrates the two non-linear problems to be solved to couple two lymphangions connected by a valve. First, we need to solve system (54), and then the resulting flow rate $q_{v}$ is imposed at both lymphangions. Once we solve the two numerical problems, then the numerical fluxes at the boundaries are provided by (68).

3.3 Imposed pressure at boundaries

In various experiments reported in the literature [19, 20, 21, 22], lymphatic collectors containing two or more valves were isolated from the animal, cannulated at each end with a glass micropipette and pressurised. Here we present a numerical method to simulate this kind of experiments. The procedure is similar to the coupling method for valves and lymphangions: both cross-sectional area and flow rate need to be found at the ending interface where the pressure needs to be imposed.

Consider a time-varying pressure $p_{I}\left(t\right)$ at a terminal interface. From pressure $p_{I}\left(t\right)$, cross-sectional area $A_{I}\left(t\right)$ can be calculated by using the inverse of the tube law (2). The flow rate $q^{*}$ can be found by applying the Riemann invariants as described in 3.2, and in this case it can be explicitly calculated as

[TABLE]

where $A^{n}$, $u^{n}$ and $\beta$ are the cross-sectional area and the velocity at the cell adjacent to the boundary at current time $t=t^{n}$ and $\beta$ is given by Eq. (66). As before, the numerical flux at the boundary is

[TABLE]

where

[TABLE]

3.4 Numerical methods for the systems of ODEs

The systems of ODEs (54) and (22) were numerically solved with a second-order implicit Runge-Kutta method using the Lobatto IIIC method. The Butcher tableau is

[TABLE]

In the next section, we present the coupling of the systems of PDEs and ODEs, through an algorithm description.

3.5 Complete algorithm

Here we provide the complete algorithm to update the solution from time $t^{n}$ to time $t^{n+1}=t^{n}+\Delta t$. When not specified, the initial conditions are: $p(x,0)=P_{in}(0)$, $u(x,0)=0$, $v(0)=0.1$, $w(0)=s(0)=I(0)=0$ and $q_{v}(0)=\xi(0)=0$.

Assume data for all variables at $t=t^{n}$. 2. 2.

Calculate the time step $\Delta t$ as explained in Section 3.1. 3. 3.

Evolve the valve flow $q$ and valve state $\xi$ of each lymphatic valve from time $t^{n}$ to $t^{n+1}$ by applying a second-order implicit Runge-Kutta method to the system of ODEs (54) and assuming the pressure difference $\Delta p$ at time $t^{n}$. 4. 4.

Calculate the numerical fluxes at the boundaries $\mathbf{F}_{\frac{1}{2}}$ and $\mathbf{F}_{M+\frac{1}{2}}$ of each lymphangion, as described in Sections 3.2 and 3.3, using the lymphatic valve flow rates at time $t^{n+1}$. 5. 5.

Using the contraction state $s$ at the current time $t^{n}$, calculate the numerical fluxes $\mathbf{F}_{i+\frac{1}{2}}$ within each domain of the lymphangions using the SLIC method (Section 14.5.3 of [112]). 6. 6.

Using the contraction state $s$ at the current time $t^{n}$, calculate the numerical sources $\mathbf{S}_{i}$ within each domain of the lymphangions using a second-order method in space and time (Chapter 19 of [112]). 7. 7.

Update the conserved variables $\mathbf{Q}$ of the PDEs of each lymphangion from time $t^{n}$ to $t^{n+1}$ according with finite volume formula (57). 8. 8.

Evolve the variables of the EFMC model of each lymphangion from time $t^{n}$ to $t^{n+1}$ by applying a second-order implicit Runge-Kutta method to the system of ODEs (22) and using the space-time averaged cross-sectional area and WSS at time $t^{n+1}$.

The EFMC model and the system of PDEs are coupled through the contraction state $s$. The variable $s$ gives the actual value of the Young modulus in Eq. (6) to be used to calculate the physical flux in Eq. (10). Observe that even though we use second-order methods for every model, the accuracy of the global algorithm is formally of only first order. This is caused by the coupling methods. As a matter of fact, we couple the set of ODEs and PDEs using only use a first-order method. There are more sophisticated high-order coupling methods in the literature, see for instance [114].

4 Results and discussion

In this section, we assemble the components of the model, implement the explained numerical methods and show numerical results. First, we numerically solve a Riemann problem for the PDEs and compare the numerical solution with the exact solution; then we show results for a single, three and ten lymphangions. We analyse lymphodynamical indexes for a wide range of boundary pressures, perform two sensitivity analyses for positive and negative pressure gradients and finally investigate the effect of defective lymphatic valves. Table 1 gives parameters used in the numerical simulations.

4.1 Test problem with piecewise initial condition: a Riemann problem

Here we solve numerically a Riemann problem for the one-dimensional model for a single lymphangion without valves. The Riemann problem is a particular Cauchy problem where the initial conditions are piecewise constant with a single discontinuity. The exact solution in subcritical regime is available for this PDE system but is not reported here. For the exact solution of the Riemann problem for arteries and veins see [99, 95, 115]. Here, the chosen initial condition for $A$ and $u$ is

[TABLE]

where $L$ is the lymphangion length and $A_{0}$ is the cross-sectional area at equilibrium, and $u(x,0)=0$. We assumed no contractions, dynamic viscosity $\mu=0$ and transmissive boundary conditions. Fig. 9 shows the numerical results using $M=40$ and $M=1000$ cells and the exact solution at the output time time $t_{output}=0.0015$ $s$. The numerical solution with $M=40$ is comparable with the exact solution, which is composed of a left rarefaction and a right shock. The numerical solution with $M=1000$ confirms that the numerical solution converges to the exact solution. This result is typical of a $2\times 2$ non-linear system of hyperbolic differential equations and is comparable with the Riemann problem for the Euler equations, shallow water equations and the blood flow equations [112].

4.2 Four representative test cases of a collecting lymphatic

Here we show four representative test cases where we simulated from one to ten lymphangions. In the first test case, we highlight the diastolic and systolic phases of the lymphatic cycle. The second test clearly shows the frequency increases as the intraluminal pressure increases. The third case shows the negative chronotropic effect given by the WSS. The four case shows an example of a bigger collecting lymphatic composed of ten lymphangions and eleven valves. Figs. 11, 12 and 13 show the numerical results of cases from one to three. From top to bottom frames we show: an illustration of lymphangions and lymphatic valves, time-varying valve states (open $\xi=1$ and closed $\xi=0$), flow rates across the valves, and pressures, velocities, diameters, WSS, Young’s modulus at the centre of the lymphangions, excitable variable and stimulus. The line colours shown from the second to the last panels refer the colour configuration shown in collector illustrated in the first panel. In the last panel, the unstable spiral-node region is represented by a blue-shaded area, while the unstable node region is represented by a red-shaded area.

4.2.1 Test case 1: representative case of a single lymphangion

Here we show a representative test where the EFMC model, the valve model and the one-dimensional model for lymph flow are coupled together. We show results for a single lymphangion with two terminal valves. Upstream (left) and downstream (right) boundary pressures are $P_{in}=3$ cmH2O and $P_{out}=4$ cmH2O, respectively. We applied the numerical method at the boundaries explained in Section 3.2. Fig. 11 shows the numerical results. As soon as the stimulus goes beyond the unstable spiral-node region and falls into the unstable node region, an action potential occurs followed by a contraction. The lymphatic pressure increases and the downstream valve opens: both flow rate across the downstream valve and the lymphatic velocity increase, the diameter decreases and the absolute value of the WSS increases. When the lymphatic pressure decreases below the downstream boundary pressure $P_{out}$, the downstream valve closes, and then the flow rate and the velocity return to zero. Then the lymphatic pressure decreases below the upstream boundary pressure $P_{in}$ and the upstream valve opens. Subsequently, the lymphangion is filled with lymph.

Fig. 10 shows the space-time numerical solution of diameter, flow, pressure and WSS. The blue and red line represent the numerical solutions close to the upstream and downstream valve, respectively. The green line represents the numerical solution at the centre of the lymphangion. The diameter decreases almost homogeneously throughout the lymphangion. The same happens to the pressure, but it can be noticed that the systolic pressure varies in space, as the maximum is reached at the upstream side, while the minimum is reached at the downstream one. The flow rate and the WSS share a similar behaviour with opposite signs. During the systolic phase, the flow rate (WSS) reaches its maximum (minimum) at the downstream side, while it reaches its minimum (maximum) at the upstream one. The red and blue lines in the flow rate are similar to the valve flow rates shown in Fig. 11. This result highlights that the mathematical model can give quantitative information throughout the length of the lymphangion.

4.2.2 Test case 2: contraction frequency increases as the intraluminal pressure increases

The numerical test shown here was inspired by the experiments performed in several works [19, 20, 21, 22] where time-varying pressures were imposed at the boundaries of the collector. More specifically, this test imitates the ramp-wise $P_{out}$ elevation shown in [20]. We simulated a collector

of three lymphangions and we imposed the following time-varying pressures

[TABLE]

where $p_{1}=p_{3}=3$ cmH2O, $p_{2}=13$ cmH2O, $t_{1}=60$ s, $t_{2}=t_{1}+120$ s, $t_{3}=t_{2}+30$ s and $t_{output}=t_{3}+60$ s. Applying the numerical methods explained in 3.3, the inlet pressure $P_{in}$ was imposed at the leftmost interface of the upstream lymphangion, while the output pressure was imposed at the rightmost interface of the downstream lymphangion. Only a negative transaxial-pressure gradient is taken into account.

Fig. 12 shows the results of the numerical simulation. Even though the upstream and downstream lymphangions contract, their pressures are controlled. The downstream pressure (red line) follows the behaviour of the imposed output pressure $P_{out}$, while the upstream pressure (blue line) is almost constantly $P_{in}$. The lymphangion at the centre responds to these changes of boundary pressures (green line). Initially, both valves close and open, but when the downstream pressure $P_{out}$ reaches a certain value ($\approx$ 11 cmH2O), the centred lymphangion cannot open the downstream valve anymore, and this can be seen by the valve state and the valve flow of the downstream lymphangion. Since the pressure of the downstream lymphangion increases during the numerical simulation, its frequency of contraction increases as well. As soon as $P_{out}$ falls back to 3 cmH2O, then the contraction frequency suddenly decreases to the initial value. The frequency of centred lymphangion is not affected by the increase in the frequency of the downstream lymphangion. This comes from the fact that the current model does not take into account the interaction between adjacent lymphangions. From this numerical results, we can see that: 1) the frequency of contractions depends on the intraluminal pressure and 2) the centred lymphangion tries to overcome the downstream time-varying pressure by increasing the end-systolic pressure, but this is possible up to a certain threshold.

4.2.3 Test case 3: contraction frequency decreases with increasing WSS

The test proposed here simulates a collector of three lymphangions with two valves and highlights the effect of the WSS on the frequency of contractions. As done for the test case 2, we imposed the following time-varying pressures at the terminal interfaces of the collector

[TABLE]

where $p_{1}=p_{4}=3.5$ cmH2O, $p_{2}=8$ cmH2O, $p_{3}=1$ cmH2O, $t_{1}=20$ s, $t_{2}=t_{1}+12$ s, $t_{3}=t_{2}+30$ s, $t_{4}=t_{3}+12$ s, $t_{5}=t_{4}+60$ s, $t_{6}=t_{5}+12$ s and $t_{output}=t_{6}+20$ s. In this test, both negative and positive transaxial pressure gradients are taken into account. The inlet pressure $P_{in}$ is fixed to $3.5$ cmH2O, while the output pressure $P_{out}$ is initially equal to the inlet pressure, it increases up to $8$ cmH2O, it decreases to $1$ cmH2O, and finally returns back to the initial value.

Fig. 13 shows the numerical results. The downstream lymphangion initially contracts at the same frequency of the other lymphangions, then as the output pressure $P_{out}$ increases, it contracts faster. The end-systolic pressure of the centre lymphangion follows the behaviour of the downstream pressure $P_{out}$. When a favourable pressure gradient occurs ($P_{in}>P_{out}$), then the absolute value of the WSS of each lymphangion increases, and thus the frequencies drastically decrease. Interestingly, when one of the lymphangions contracts, then the flow at the downstream valve decreases from $\approx$ 50 $\mu$L min*-1* to $\approx$ 18 $\mu$L min*-1*. This means that the contractions of the lymphangions decrease the outflow, and therefore the negative chronotropic effect given by the increment of the WSS gives an overall positive effect on the averaged outflow.

4.2.4 Test case 4: representative case of a collecting lymphatic composed of ten lymphangions

Here we illustrate a representative example of a collecting lymphatic composed of ten lymphangions and eleven valves. We applied the numerical method at the boundaries explained in Section 3.2 with prescribed upstream and downstream pressures $P_{in}=P_{out}=4$ cmH2O. We used $M=20$ computational cells to discretize each one-dimensional lymph vessels and we set the output time $t_{output}=93$ $s$.

Fig. 14 shows the numerical results of the simulation. Figure shows an illustration of the collector at time $t=89.4$ (shown by the black vertical line), and time-variation of flow rate, pressure and diameter at the centre of each lymphangion. Lymphatic contractions occur without prescribed delays between lymphangions. When the leftmost lymphangion contracts, the adjacent lymphangion becomes stimulated as its diameter increases and will thus contract as well. This leads to a chain of consecutive contractions, which occur in a non-linear manner.

4.3 Pressure versus normalised cross-sectional area (PA) plots for a single lymphangion

In this section we show plots of pressure against normalised cross-sectional area (PA) of a single lymphangion with fixed upstream pressure $P_{in}=3$ cmH2O and different downstream pressures $P_{out}$, from $3$ to $10$ cmH2O. The aim of this exercise is to show that the numerical results of the mathematical model imitate the experimental measurements of the pressure-volume relationship [20]. Fig. 15 shows the numerical results and also the contracted and relaxed tube laws (green lines). Results show a qualitatively good agreement with [20, 21]. As the downstream pressure $P_{out}$ increases, the PA plots shrink and the systolic pressure increases. As a consequence, the stroke work and the ejection fraction decrease towards zero. The systolic pressure can increase up to a certain level, depending on the baseline pressure. As a matter of fact, in the current case $P_{in}=3$ cmH2O, the maximum systolic pressure is $\approx 11$ cmH2O and decreases as $P_{in}$ decreases. For instance, for $A/A_{0}=4$ e $P_{in}\approx 2$ cmH2O, the maximum systolic reachable pressure is $\approx$ 8 cmH2O. From the literature, we know that lymphangions can increase the strength of contraction under stresses [20], which means that the slope of the Ending Systolic Pressure-Volume Relationship (ESPVR) increases [21] or, in analogy to our terminology, that the maximum Young’s modulus $E_{max}$ increases somewhat. For an example of a mathematical model with adaptation of contractility, see Caulk et al. [81]. This represents a limitation of our current mathematical model as the tube laws do not undergo modification under continuous stresses.

4.4 Analysis of lymphodynamical indexes by varying $P_{in}$ and $P_{out}$

The aim of this study is to show lymphodynamical indexes for several combinations of $P_{in}$ and $P_{out}$ from 0 to 9 cmH2O. We constructed a grid of points $\left(P_{out},P_{in}\right)$ with all the possible combinations of $P_{in}=\left(0,0.1,\dots,8.9,9\right)$ and $P_{out}=\left(0,0.1,\dots,8.9,9\right)$. For each combination of $P_{in}$ and $P_{out}$, we simulated a single lymphangion with two terminal valves, $t_{output}=160$ $s$ and $M=20$ computational cells. The total number of simulations was $91\times 91=8281$. We applied the numerical method at the boundaries as explained in Section 3.2. We analyzed the indexes reported in Table 2.

CPF is the flow produced only by the contraction and does not take into account the passive flow induced by a positive transaxial-pressure gradient. To understand the influence of the lymphatic pump on the total flow during a lymphatic cycle, we introduce the Calculated Pump Flow Index, which indicates if the flow is produced by contractions (CPFI $\approx$ 1) or by a positive transaxial-pressure gradient (CPFI $\approx$ 0). We observe that assuming constant frequency of contraction, the CPFI can also be calculated as

[TABLE]

which is the ratio between the ejected lymph volume due to the lymphatic contraction and the total lymph volume difference. Davis et al. [20] performed an analysis of the lymphodynamical indexes by keeping fixed the upstream pressure $P_{in}$ and varying the downstream pressure $P_{out}$ from 0 to 18 cmH2O. Similar experiments were performed by Scallan et al. [21] where the effect of increasing the upstream pressure $P_{in}$ while keeping fixed the downstream pressure $P_{out}$ was studied. In the present paper, both experimental setups are comprised in the numerical simulations.

Fig. 16 shows the numerical results of the analysis through contour plots and isolines. To describe the figure, we divide the $P_{in}$-$P_{out}$ space in two regions: 1) the negative transaxial-pressure gradient $\Delta P=P_{in}-P_{out}<0$ region (lower triangle) and 2) the positive transaxial-pressure gradient $\Delta P=P_{in}-P_{out}>0$ region (upper triangle). The two regions are therefore divided by the bisector $P_{in}=P_{out}$.

Negative pressure gradient $\Delta P=P_{in}-P_{out}<0$. Here the upstream and downstream valves open non-linearly and close during the lymphatic cycle (results not shown). The frequency increases as $P_{in}$ rises, and this is in agreement with Fig. 6. This comes from muscle-stretch feedback from the EFMC model in Eq. (24). The frequency does not increase when $P_{out}$ rises. This might be surprising because it is well-known that contraction-waves propagate between lymphangions through gap-junctional communications [35, 116]. Therefore, if the downstream pressure rises, then it would be natural to expect a rise in frequency among the neighbouring lymphangions. However, we modelled only a single lymphangion with variable downstream and upstream pressures, and we did not take into account the interaction between adjacent lymphangions. Moreover, the gap-junctional communications between lymphangions have not been modelled in this paper for the sake of simplicity, though we speculate that it is mathematically possible. EF tends to decrease as the $P_{out}$ increases, while increases when $P_{in}$ increases. FPF combines both frequency and EF: it increases when $P_{in}$ rises, and it decreases when $P_{out}$ increases. Intriguingly, the maximum of FPF is when $P_{in}\approx P_{out}$ and near 6 cmH2O. For higher pressures, FPF decreases. The results for FPF are not comparable with Davis et al. [20], as the frequency remain constant when $P_{out}$ rises. The SV and AMP follow the same behaviour of EF. ESD increases only when $P_{out}$ rises, and it remains constant when $P_{in}$ increases. On the contrary, EDD remains constant when $P_{out}$ increases, and this is in agreement with [20]. SW is maximum for $P_{out}\approx 7$ cmH2O and $P_{in}\approx 3$ cmH2O, and it tends to decrease elsewhere. The mean flow and CPF are comparable almost everywhere. This means that lymph flow is produced only by contractions, and not by a passive flow induced by pressure gradients. This is summarised in CPFI, which is almost 1 everywhere in this region. Finally, note that WSS is generally low in this region.

Positive pressure gradient $\Delta P=P_{in}-P_{out}>0$. Here the valves remain open for most of the time during the lymphatic cycle (results not shown). The mean flow increases as $P_{in}$ rises, and it decreases as $P_{out}$ increases. Subsequently, the absolute value of WSS rises when $P_{in}$ increases, and this gives a negative chronotropic effect on the frequency. In other words, the frequency decreases as the absolute value of WSS increases. This comes from function $f_{NO}$ in Eq. (24): the greater the absolute value of the WSS, the greater the negative chronotropic effect on the lymphatic contractions. The decrease in frequency is controlled by parameters $k_{NO}$, which in these simulations is set to $k_{NO}=0.8$. Lower values of $k_{NO}$ would have weakened the negative chronotropic effects, leading to higher rates of lymphatic contractions. Observe that the CPF differ from the mean flow insofar as the CPF only takes into account the flow given by contractions. This is summarised in CPFI which is almost zero in this region. The ESD and EDD increase when $P_{in}$ and $P_{out}$ rise. EF, SV and AMP share a similar behaviour and reach their maximum value at $P_{in}=P_{out}\approx 4$ cmH2O.

4.5 Sensitivity analyses

[FIGURE:]

The mathematical model for lymphatic collectors proposed here depends on several parameters shown in Table 1, and the lymphodynamical indexes studied in Section 4.4 are strongly related to these parameters. To investigate the influence of each parameter on the indexes, we performed a sensitivity analysis following the idea of [94]. This procedure was already applied, for instance, to find the most relevant model parameters of a zero-dimensional model of the cardiovascular system [59].

The method is divided into a local and global analysis. In the local analysis we calculated $N$ local sensitivity matrixes $S_{i,j}^{k}$, for $k=1,\dots,N$, as follows: starting from the reference value in Table 1, we randomly varied each parameter from $70$ to $130\%$ and obtained a new set of parameters. Here, the baseline value for $k_{NO}$ was set to 0.5. With this varied set of parameters, we calculated the local sensitivity matrix as follows

[TABLE]

where $\mathbf{X}=\left(x_{1},x_{2},\dots,x_{m}\right)$ is the vector of the varied model parameters, $\mathbf{P}=\left(P_{1},P_{2},\dots,P_{n}\right)$ is the vector of the lymphodynamical indexes and $\mathbf{S}^{k}=\left(S_{i,j}^{k}\right)_{i,j}$ is local sensitivity matrix. The value $S_{i,j}^{k}$ represents the non-dimensional percentage change in index $P_{j}$ induced by a small change in parameter $x_{i}$. For instance, if the model parameter $x_{i}$ does not influence index $P_{j}$, then $S_{i,j}^{k}$ will be almost zero. Viceversa, if there is a significant influence of $x_{i}$ on $P_{j}$, then the absolute value of $S_{i,j}^{k}$ will be greater than zero. A positive sign of $S_{i,j}^{k}$ indicates that an increase of parameter $x_{i}$ induces an increase of index $P_{j}$. Viceversa, a negative sign of $S_{i,j}^{k}$ indicates that an increase of parameter $x_{i}$ induces a decrease of index $P_{j}$.

Subsequently, in the global analysis we performed a statistical analysis of $S_{i,j}^{k}$ by calculating its mean $\bar{S}_{i,j}$ and its standard deviation $\sigma_{i,j}$, namely

[TABLE]

[TABLE]

A large standard deviation $\sigma_{i,j}$ indicates a strong correlation of the studied parameter with the remaining parameters in determining the sensitivity index. To calculate (80) and (81), we removed possible outliers by discarding the data below the $3$rd percentile and above the $97$th percentile.

The partial derivative in Eq. (79) were approximated using a second-order finite difference method based on a percentage change of the parameter as follows:

[TABLE]

where

[TABLE]

The parameter $\epsilon$ was chosen as $\epsilon=0.05$. Compared to [94] and [59], we did not constructed a stratified sampling space, but rather a simple random variation in the considered range. Here, we studied parameter ${A_{Ca}}/{A_{0}}$ with baseline value $7.75$.

Based on the results of the lymphodynamical indexes shown in Section 4.4, we performed two sensitivity analyses: one for a positive pressure gradient $\Delta P=P_{in}-P_{out}>0$ and one for a negative pressure gradient $\Delta P=P_{in}-P_{out}<0$. For each analysis, we calculated $N=500$ local sensitivity matrixes. We used $M=20$ computational cells to discretize the one-dimensional lymph vessel and output time $t_{output}=60$ $s$. Results are shown in Tables 3 and 4. In both tables, results are shown in the following manner: the second column from the left shows parameters $\mathbf{X}$ of the model while the first column shows their means $\pm$ SDs. The second row from the top shows the studied lymphodynamical indexes $\mathbf{P}$ while the first row shows the resulting means $\pm$ SDs. Symbol $"-"$ indicates that the absolute value of the sensitivity index $\bar{S}_{i,j}$ was less than 0.5 and therefore parameter $x_{i}$ did not influence index $P_{j}$. Likewise, green-coloured parameters show an influence between 25 and 50, blue-coloured parameters show an influence between 50 and 100, and red-coloured parameters show an influence greater than 100.

[FIGURE:]

Negative pressure gradient $\Delta P=P_{in}-P_{out}<0$. Here we used $P_{in}$ = 3 cmH2O and $P_{out}$ = 4 cmH2O. The cross-sectional radius at equilibrium $r_{0}$ positively influence indexes SV, CPF, ESD, EDD, WSS and mean flow. Among the parameters of the EFMC model, $a_{2}$ is the most influential parameter. Indeed, it is the threshold to change the nature of the stationary point described in Section 2.2.1 from stable to unstable. Parameter $R_{I}$ does not have a remarkable influence on the indexes. Likewise for parameters $k_{rel}$. The frequency, and thus FPF and CPF, are strongly influenced by $A_{Ca}/A_{0}$ and $n_{Ca}$. Indeed, the greater these parameters, the lower the frequency. Results for $n_{Ca}$ are in agreement with Frame 6c. Moreover, these parameters also influence the mean flow. Then parameters $k_{NO}$, $\tau_{NO}$, $n_{NO}$, which are related to WSS and passive flows, do not affect the lymphodynamical indexes since $\Delta P<0$. The parameter of the valve model $K_{vo}$, $K_{vc}$ and $L_{eff}$ do not affect the indexes. The fluid properties of lymph $\mu$ and $\gamma$ only affect the WSS, while the density $\rho$ has no effects on the indexes. The maximum and minimum Young modulus $E_{max}$ and $E_{min}$ affect the ESD and EDD, respectively, and also influence most of the parameters, such as the frequency, EF, SV, FPF, CPF and mean flow. To conclude, among the analysed parameters, nothing appears to influence the ESP.

Positive pressure gradient $\Delta P=P_{in}-P_{out}>0$. Here we used $P_{in}$ = 4 cmH2O and $P_{out}$ = 2 cmH2O. The cross-sectional radius at equilibrium $r_{0}$ influences the CPF, SV, WSS, mean flow, ESD, EDD and WSS. The most influential parameter of the EFMC model is $a_{2}$, followed by $a_{3}$ and $b_{1}$, which influences the frequency, CPF, FPF and WSS. Radius $R_{I}$ and parameter $k_{rel}$ do not affect the lymphodynamical indexes. As before, an increase of parameters $n_{Ca}$ and $A_{Ca}/A_{0}$ influence the WSS and leads to a negative chronotropic effect, which in turn influences the FPF and the CPF. An increase of parameter $k_{NO}$ decreases the frequency, indeed the greater this parameter, the greater the influence of the contraction inhibition given by the WSS. On the contrary, the increase of parameter $\tau_{NO}$ increases the frequency. Results for $k_{NO}$ and $\tau_{NO}$ are in agreement with Frames 6f and 6e, respectively, of Fig. 6. To a lesser degree, parameters $L_{eff}$, $\rho$, $\mu$ and $\gamma$ influence the frequency, FPF, CPF, WSS and the mean flow. Among these parameters, the greatest effect is seen on the mean flow given by the dynamic viscosity $\mu$. Between the minimum and maximum Young’s modulus, the most influential one is $E_{min}$, as it affects the frequency, SV, FPF and CPF. It also affects WSS and EDD, though to a lesser extent.

4.6 A quantitative study on the lymphodynamical effect of stenotic and regurgitant lymphatic valves

The mathematical model for collectors proposed in the present paper includes a well-established model for valves proposed by Mynard et al. [91]. It has already been used for the heart valve modelling [63], as well as for the venous valves [92]. More interestingly, the model allows for a quantitative study of the effect of stenotic and regurgitant valves. For instance, the model was already used to study the impact on brain haemodynamics of bilateral stenotic and regurgitant valves of the internal jugular veins [92, 117].

Here we aim to analyse the lymphodynamical effect of stenotic and regurgitant lymphatic valves. We modelled a single collector with three lymphangions and two valves. Here we used $M=20$ computational cells to discretize the one-dimensional lymph vessel. We simulated a collector cannulated at each end, that is, we imposed a fixed pressure at the leftmost and rightmost interfaces of the collector, as described in 3.3. The imposed pressures were $P_{in}=3$ cmH2O and $P_{out}=4$ cmH${}_{2}O$. A healthy valve is characterized by parameters $M_{st}=1$ and $M_{rg}=0$. A stenotic valve can be simulated by reducing the maximum effective area $A_{eff,max}=M_{st}A_{0}$, and this can be done by decreasing parameter $M_{st}$. Thus, a severe stenosis can be modelled by setting $M_{st}\approx 0$. Viceversa, a regurgitant valve can be simulated by increasing the minimal effective area $A_{eff,min}=M_{rg}A_{0}$, and this can be done by increasing parameter $M_{rg}$. A severe regurgitant valve can be modelled by setting $M_{rg}\approx 1$. For further details, see the Eqs. reported in 2.3. See also for [91, 63, 92].

We consider four possible situations: a left stenotic valve, a right stenotic valve, a left regurgitant valve and a right regurgitant valve. The numerical results of the centred lymphangion are shown in Fig. LABEL:fig:incompetentValves. From top to bottom we show in order: the PA loops, the time-varying lymphatic pressure, diameter and flow at the centre of the lymphangion. The boundary pressures $P_{in}$ and $P_{out}$ are shown in the PA loops and in time-varying pressure plots. The first two columns show results for the left and right stenotic valves, while the remaining two columns show results for the left and right regurgitant valves. Parameters $M_{st}^{L/R}$ and $M_{rg}^{L/R}$ were varied from 0 to 1.

A left stenotic valve diminishes the inflow from the upstream valve. This results in the following: the greater the severity of the left stenosis, the greater the time required to fill the centre lymphangion after a contraction. For the tests considered here, contractions occur at a frequency of $\approx$ 7 min*-1*, which means approximately every $\approx$ 8.6 s. For a severe left stenosis ($M_{st}^{L}$ = 0.025), the time required to fill the lymphangion is $\approx$ 16.5 s. A severe reduction of the EDD can happen when the lymphangion does not have enough time to fill itself, and this may happen when the contraction period is less than 8.6 s. To verify this hypothesis, we performed additional simulations with a left stenotic valve, varying the frequency and for different severities of the stenosis. We set $f_{min}=f_{Ca}$ from 4 to 24 min*-1* and calculated the resulting frequency, CPF and EDD. The numerical results are shown in Fig. 17. For a mild stenosis ($M_{st}^{L}>0.1$) and low frequencies, the CPF does not suffer any changes, but as soon as the frequency increases (e.g. above approximately 10 min*-1* for $M_{st}^{L}=0.1$), the CPF decreases depending on the severity of the stenoses. At the frequency of $f=25$ min*-1* and $M_{st}^{L}=0.1$, the CPF reduces from $\approx$ 67.4 $\mu$L h*-1* to $\approx$ 33.0 $\mu$L h*-1*, that is it reduces of the 51 $\%$. For even more severe left stenoses ($M_{st}^{L}<0.1$), the CPF drastically decreases and the lymphangion becomes unable to push the lymph forward. At the frequency of $f=25$ min*-1* and $M_{st}^{L}=0.025$, the CPF reduces of the 93 $\%$, namely it reduces to $\approx$ 4.1 $\mu$L h*-1*. This comes from a decrease of the EDD for high frequencies. The PA loops for different frequencies and a severe left stenosis are also shown. The higher the frequencies, the greater the shrinkage of the PA loops. Overall, a left stenosis causes a decrease of the CPF for high frequencies of contractions.

A right stenotic valve drastically increases the ESP and the ESD. This comes from the difficulties for the lymph to be pushed downstream through a stenotic passage. As a matter of fact, the outflow greatly decreases.

From simulating right and left stenotic valves we speculate the following: a stenotic valve causes an increase of the systolic peaks in the upstream lymphangions and maintains almost unchanged the downstream pressures. Moreover, it causes a reduction of the CPF for high frequencies of contractions in the downstream lymphangions. It would be interesting to experimentally quantify and compare the effect of stenotic lymphatic valves in regions where the frequency of contraction is high and low.

A left regurgitant valve has a great impact on the effective pump flow, namely the real amount of flow ejected from the lymphangion. As the severity of the left regurgitant valve increases, backflows increase and this can be seen with the negative values of flows. This means that during contractions, the lymph is ejected backwards into the upstream lymphangion, and not forward into the downstream one. Moreover, the ESD diameter decreases and for a severe left regurgitant valve the downstream valve stays closed most of the time (result not shown) insofar as the ESP decreases.

To conclude, a right regurgitant valve increases the leakage from the downstream valve, even for small values of $M_{st}^{R}$. This results in increasing the EDP from 3 cmH2O to 6 cmH2O, which corresponds to the downstream boundary pressures $P_{out}$. For severe right regurgitant valves, the upstream valve does not open during the lymphatic cycle (results not shown).

The effects of regurgitant and stenotic valves are summarised in the lymphodynamical indexes shown in Table 5. The table shows from the left to the right column: the indexes, the stenotic valve case, the regurgitant valve case and the healthy control. For each valve dysfunction, results for a moderate and a severe case are shown. For a left stenotic case, there are almost no variations in any of the indexes. As discussed before, problems arise for high frequencies of contractions. For a right stenotic valve, EF, SV, FPF and CPF halve, while the ESP increases. The cases of regurgitant valves show interesting properties on some of the lymphodynamical indexes. As a matter of fact, EF, SV, FPF and CPF indexes do not indicate any reduction in the pumping performance. Instead, based on the results in Table 5, it seems that the pumping action of the lymphangion has undergone improvements with the incompetence of the valves. For instance, with a left regurgitant valve case, EF, SV, FPF and CPF increase. The same happens for the right regurgitant valve case for indexes SV, FPF and CPF, though FPF and CPF might have increased because the frequency was increased. This is obviously misleading: since a significant amount of lymph is flowing retrograde due to the deficit, the effective time-averaged flow is approximately zero. Thus, we would expect CPF to be zero. As it was pointed out by Scallan et al. [118], indexes EF, SV, FPF and CPF are usually assumed to represent forward lymph flow and therefore they can give inaccurate results for dysfunctional valves. Interestingly, CPFI gives unrealistic results and shows the inaccuracy of CPF. This index should be bounded between 0 and 1, but since the mean flow reduces to zero while CPF even increases depending on the severity and on the side of the regurgitant valve, the CPFI results to be much greater than 1 (CPFI $\approx$ 103 for a severe right regurgitant valve). Thus, we speculate that index CPFI can indicate a dysfunctional regurgitant valve.

5 Concluding Remarks

In this paper, we have proposed a one-dimensional model for collecting lymphatics coupled to a novel Electro-Fluid-Mechanical Contraction (EFMC) model for dynamical contractions based on a modified FitzHugh-Nagumo model for action potentials and to a modern lumped model for valve dynamics. We performed an analysis of lymphodynamical indexes for a wide range of upstream and downstream pressure combinations. Numerical experiments showed that the contraction frequency resulting from the mathematical model strongly depends on the baseline frequency contraction, on the intramural pressure and on the wall-shear stress. We also performed numerical tests inspired by experiments where terminal lymphangions were cannulated, and the numerical results showed a good agreement with the experimental trend. The most influential model parameters were found by performing two sensitivity analyses for positive and negative pressure gradients. We then quantified the effect of stenotic and regurgitant valves. A stenotic valve caused an increase of the systolic peaks in the upstream lymphangions and maintained almost unchanged the downstream pressures. Moreover, it caused a reduction of the CPF in the downstream lymphangions for high frequencies of contractions (up to 93$\%$ for a severe stenosis), while the CPF remained unaltered for low frequencies. A regurgitant valve was unable to prevent backflows, and this resulted in zero net flows during the lymphatic cycle. The lymphodynamical indexes EF, SV, FPF and CPF were misleading for a regurgitant valve, insofar as from these indexes, it seemed that the pumping action of the lymphangion had undergone improvements with the incompetence of the valves. As a matter of fact, these indexes are usually assumed to represent forward lymph flow and therefore they can give inaccurate results for dysfunctional valves. To overcome this problem, we introduced the Calculated Pump Flow Index (CPFI), ideally bounded between 0 and 1, that indicates if the lymph is driven by contractions (CPFI = 1) or by a passive flow induced by a positive transaxial-pressure gradient (CPFI = 0). For a regurgitant valve, CPFI gave unrealistic results (CPFI $\approx$ 103 for a severe regurgitant valve).

Several improvements can be made to the mathematical model. The EFMC model can be generalised to include a diffusion term. The resulting model would be able to simulate contraction waves travelling in the lymphatic wall and induced by pacemaker cells. In this way, gap-junctional communications would be included in the mathematical model. Then, the tube law can be improved in several ways. For instance, it can include a viscous term for the vessel wall. High-order methods can be used to improve the accuracy of the numerical solution, and these methods should properly couple the systems of ODEs with the one-dimensional model. Moreover, variable geometric parameters have been neglected, but they can be included in the simulation, as the mathematical formulation proposed in the present work allows for their presence. Finally, we also assumed an ex vivo setting, in which there is no interaction with the environment, such as skeleton muscle contraction or lymphatic contractions due to neuro-activities.

The present EFMC model has been coupled to a one-dimensional model, but we speculate that it can be coupled to lumped parameter lymphatic models. Networks of collecting lymphatics can be simulated through the proposed model, but the lack of quantitative experimental measurements represents a great problem in validating the numerical results. We believe that the current mathematical model of collecting lymphatic can be coupled to multi-scale, closed-loop mathematical model of the cardiovascular system and can give quantitative lymphodynamical information in healthy and pathological cases.

Acknowledgment

The authors gratefully acknowledge the suggestions given by Prof. Christian Vergara from the Department of Mathematics in Milan.