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

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 deformable lymphatic vessel is a nonlinear system of hyperbolic Partial Differential Equations (PDEs). The EFMC model combines the electrical activity of lymphangions (action potentials) with fluid-mechanical feedback (circumferential stretch of the lymphatic wall and wall shear stress) and lymphatic vessel wall 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 stability analysis of the stationary state of the EFMC model. Contractions turn out to be triggered by the instability of the stationary state. Overall, the EFMC model allows emulating the influence of pressure and wall shear stress on the frequency of contractions observed experimentally. Lymphatic valves are modelled by extending an existing lumped-parameter model for blood vessels. Modern numerical methods are employed for the one-dimensional model (PDEs), for the EFMC model and valve dynamics (ODEs). Adopting the geometrical structure of collecting lymphatics from rat mesentery, we apply the full mathematical model to a carefully selected suite of test problems inspired by experiments. We analysed several indices of a single lymphangion for a wide range of upstream and downstream pressure combinations which included both favourable and adverse pressure gradients. The most influential model parameters were identified by performing two sensitivity analyses for favourable and adverse pressure gradients.

Appendices

Appendices

1.1 A Mathematical analysis of the one-dimensional lymph flow equations

Fig. 13

Framework for a finite volume scheme. Top: illustration of a computational volume for a lymphangion. Bottom: illustration of the space–time control volume

Here, we study the mathematical properties of (10) assuming constant parameters along the lymphatic vessel. The equations in (10) are a generalization of the one-dimensional blood flow equations (Toro 2016). 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 (10) can be written in quasi-linear form as

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

(51)

where

$$\begin{aligned} \mathbf {A}(\mathbf {Q},t)=\begin{bmatrix} 0&1 \\ \frac{A}{\rho }K\partial _{A}\psi -u^2&2u \\ \end{bmatrix}\;, \quad \mathbf {S}(\mathbf {Q})= \begin{bmatrix}0\\-\frac{f}{\rho }\end{bmatrix} \;. \end{aligned}$$

(52)

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

$$\begin{aligned} \lambda _1=u-c\;, \quad \lambda _2=u+c\;, \end{aligned}$$

(53)

where c is the wave speed

$$\begin{aligned} c= & {} \sqrt{\frac{A}{\rho }K\partial _{A}\psi }\nonumber \\= & {} \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] }\;. \end{aligned}$$

(54)

We assume parameters $m\ge 0$, $n\le 0$, $z\ge 0$, and $C\ge 0$ for the tube law. Thus, the wave speed c is always real. The wave speed increases actively during contraction as it depends on the coefficient K. This means that during lymphatic contraction, the lymphatic wall becomes stiffer and waves propagate at a faster rate. The eigenvectors of $\mathbf {A}$ are

$$\begin{aligned} \mathbf {R}_1=\gamma _1 \begin{bmatrix}1\\u-c\end{bmatrix} \;, \quad \mathbf {R}_2=\gamma _2 \begin{bmatrix}1\\u+c\end{bmatrix} \;, \end{aligned}$$

(55)

where $\gamma _1$ and $\gamma _2$ are arbitrary scaling factors. It can be shown that system (10) 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 Toro (2016) and references therein, the $\lambda _1$ and $\lambda _2$ characteristic fields are genuinely nonlinear outside the locus of the following function

$$\begin{aligned} G\left( \frac{A}{A_0}\right)= & {} m\left( m+2\right) \bigg (\frac{A}{A_0}\bigg )^m-n \left( n+2\right) \bigg (\frac{A}{A_0}\bigg )^n\nonumber \\&+\,Cz\left( z+2\right) \bigg (\frac{A}{A_0}\bigg )^z\;. \end{aligned}$$

(56)

With the choice of parameters m, n, z, and C in Table 1, there exists at least one solution of $G\left( \frac{A}{A_0}\right) =0$. This means that the two characteristic fields are neither genuinely nonlinear nor linearly degenerate. The consequences of this are unclear to the authors and might require further investigations. See LeFloch (2002) for details. The Generalized Riemann Invariants (GRIs) for $\lambda _1$ and $\lambda _2$ characteristic fields are, respectively,

(57)

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.

B 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.

1.1 B.1 A finite volume method for the one-dimensional model

Consider the system of m hyperbolic balance laws

$$\begin{aligned} \partial _t\mathbf {Q}+\partial _x\mathbf {F}(\mathbf {Q})=\mathbf {S}(\mathbf {Q})\;. \end{aligned}$$

(58)

By integrating (58) 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

$$\begin{aligned} \mathbf {Q}^{n+1}_i=\mathbf {Q}^n_i-\frac{\varDelta t}{\varDelta x}\big (\mathbf {F}_{i+\frac{1}{2}}-\mathbf {F}_{i-\frac{1}{2}}\big )+\varDelta t \mathbf {S}_i\;, \end{aligned}$$

(59)

with definitions

(60)

(61)

Eq. (60) gives the spatial–integral average at time $t=t^n$ of the conserved variable $\mathbf {Q}$ while Eqs. (60) 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 $\varDelta x=x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}}$ and $\varDelta t=t^{n+1}-t^n$, respectively. Finite volume methods for (58) depart from (59) to (60), where integrals are approximated, and then formula (59) becomes a finite volume method, where the approximated integrals in (60) 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. 13 for an illustration of the finite volume framework. To compute the time step $\varDelta t$, the Courant–Friedrichs–Lewy condition is applied for each lymphangion

$$\begin{aligned} \displaystyle \varDelta t^j = CFL \frac{\varDelta x^j}{\max \limits _{i=1,\dots , M^j} \left( |u_i^j|+c_i^j \right) }\;, \end{aligned}$$

(62)

with $CFL=0.9$. Superindex j indicates the j-th lymphangion. Then, the time step $\varDelta t$ to be used is the minimum of all the time steps, namely $\varDelta t=\min \limits _{j}\left( \varDelta 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}}$) (Toro and Billett 2000). 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 Toro (2009) and references therein. The numerical source was approximated using a second order in space and time method; see Chapter 19 of Toro (2009). 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.

1.2 B.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. 13. There are three possible modelling 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 (45), where the pressure difference $\varDelta p$ in (39) is evaluated at the current time $t^n$ using either of the two lymphangions, or one of the lymphangions and a prescribed time-varying pressure. Specifically, at $t=t^n$ the pressure difference $\varDelta p(t^n)$ is

$$\begin{aligned} \varDelta p(t^n) = p_u(t^n) - p_d(t^n)\;, \end{aligned}$$

(63)

where values $p_u(t^n)$ and $p_d(t^n)$ are

(64)

and

(65)

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

(66)

Once we numerically solve system (45), 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. (2008). This method has already been used in Müller and Toro (2014), Contarino et al. (2016). To find $A^*$, we solve the following nonlinear algebraic equation based on the Riemann invariant

$$\begin{aligned} \mathcal {F}(A^*)\,{:=}\,q_v^{n+1}+A^*\left( -u^n+\beta \int _{A^n}^{A^*} \frac{c\left( \tau \right) }{\tau }\mathrm {d}\tau \right) =0, \end{aligned}$$

(67)

using the Newton-Raphson iterative method. $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

(68)

Then, the numerical flux at the boundary is

$$\begin{aligned} \mathbf {F}_{\frac{1}{2} \text { or } M+\frac{1}{2}}=\mathbf {F}\left( \mathbf {Q}^*\right) , \end{aligned}$$

(69)

where

$$\begin{aligned} \mathbf {Q}^*= \begin{bmatrix}A^*\\q_v^{n+1}\end{bmatrix} . \end{aligned}$$

(70)

As shown in Fig. 14, when a valve is interposed between two lymphangions, then the nonlinear problem (67) has to be solved twice: one for the upstream lymphangion ($\beta =1$) and one for the downstream lymphangion ($\beta =-1$).

Fig. 14

Illustration of the coupling method between two lymphangions and one valve. Riemann invariants are used to couple the flow through the valve and the two lymphangions

1.3 B.3 Imposed pressure at boundaries

The numerical procedure to impose a pressure in one of the extremities of a lymphatic vessel is similar to the coupling method for valves and lymphangions. 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 B.2, and in this case it can be explicitly calculated as

$$\begin{aligned} q^*=A_{I}\left( t^n\right) \left( u^n-\beta \int _{A^n}^{A_I\left( t^n\right) } \frac{c\left( \tau \right) }{\tau }\mathrm {d}\tau \right) \;, \end{aligned}$$

(71)

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. (68). As before, the numerical flux at the boundary is given by (69) with

$$\begin{aligned} \mathbf {Q}^*= \begin{bmatrix}A_I\left( t^n\right) \\q^*\end{bmatrix} \;. \end{aligned}$$

(72)

1.4 B.4 Numerical method for the systems of ODEs

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

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

1.5 B.5 Complete algorithm

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

• 1. Assume data for all variables at $t=t^n$.

• 2. Calculate the time step $\varDelta t$ as explained in Section B.1.

• 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 (45) and assuming the pressure difference $\varDelta p$ at time $t^n$.

• 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 B.2 and B.3, using the lymphatic valve flow rates at time $t^{n+1}$.

• 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 Toro (2009)).

• 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 Toro (2009)).

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

• 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 (17) 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 coefficient K(t) in Eq. (5) to be used to calculate the physical flux in Eq. (12). Observe that even though we use second-order methods for every model, the accuracy of the global algorithm is formally of order one. This is caused by the coupling methods. As a matter of fact, we couple the set of ODEs and PDEs using only a first-order method. There are more sophisticated high-order coupling methods in the literature; see, for instance, Borsche and Kall (2016).

C Sensitivity Analysis

The method is divided into a local and global analysis. In the local analysis, we calculated N local sensitivity matrices $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:

$$\begin{aligned} S_{i,j}^k=\left| \frac{x_i}{P_j(\mathbf {X})}\right| \frac{\partial P_j\left( \mathbf {X}\right) }{\partial x_i}\times 100, \end{aligned}$$

(73)

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 lymphatic indices, 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 relative change in $P_j$ to the relative change in parameter $x_j$, expressed as a percentage. If the model parameter $x_i$ does not influence index $P_j$, then $S_{i,j}^k$ will be zero. Vice versa, 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. For instance, if $1\%$ change in $x_i$ leads to $1\%$ change in $P_j$, then $S_{i,j}^k$ is $100\%$. A positive sign of $S_{i,j}^k$ indicates that an increase of parameter $x_i$ induces an increase of index $P_j$. Vice versa, 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}$. 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 $\bar{S}_{i,j}$ and $\sigma _{i,j}$, we removed possible outliers by discarding the data below the third percentile and above the 97th percentile.

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

$$\begin{aligned} S_{i,j}^k\approx \frac{\text {sgn}(x_i)}{\left| P_j(\mathbf {X})\right| } \frac{P_j\left( \mathbf {X}^{i,\epsilon _+}\right) -P_j\left( \mathbf {X}^{i, \epsilon _-}\right) }{2\epsilon }\times 100, \end{aligned}$$

(74)

where $\mathbf {X}^{i,\epsilon _{\pm }}=\left( x_1,\dots ,x_i\left( 1\pm \epsilon \right) ,\dots ,x_m\right) .$ The parameter $\epsilon$ was chosen as $\epsilon =0.05$. Compared to van Griensven et al. (2006), we did not construct a stratified sampling space, but rather a simple random variation in the considered range.