Reproducibility study to simulate a chain of collapsible contracting lymphangions with progressive valve closure

Abstract

This paper presents a reproduced computational model of lymphatic collecting vessels. The model was originally introduced by Bertram et al. (2011), and comprises a series of contractile segments, known as lymphangions, interconnected by secondary lymphatic valves. The model focused on elucidating the pumping behavior of contracting lymphangions and deriving pump-characteristic curves, by incorporating pressure-dependent valve resistance, passive elasticity, and active contraction terms in multiple lymphangions. The original model was implemented and solved using MATLAB software. In our work, we reproduce Bertram et al. (2011)’s original model within the CellML framework, and present the results.

Keywords:CellMLOpenCORcomputational modellymphatic vesselslymphangionvalve

1Introduction

The primary paper of Bertram et al. (2011) introduced a model comprising a series of lymphangions connected via inter-lymphangion valves. The objective of their study was to take the first step toward developing a comprehensive model of the lymphatic system. The original contribution of Bertram et al. (2011)’s work lies in incorporating a pressure-dependent valve resistance and the consideration of the constitutive equation of the vessel wall, including passive elasticity and active contraction terms.

The model of Bertram et al. (2011) played a significant role in collecting lymphatic modelling, and subsequent studies Bertram et al., 2014Bertram et al., 2016Bertram et al., 2017 have improved and validated it by using ex-vivo experimental data published in later years Davis et al., 2011Davis et al., 2012. Over the last decade, there has been a surge of interest in lymphatic-related research Don et al., 2023. Despite this, the lymphatic system remains one of the least explored areas in biomedical research. Freely available and open-source reproducible lymphatic models aim to advance lymphatic research.

2Model description

A schematic diagram of the lymphangion-chain model displays the contractile vessel segments, each characterized by its diameter (D). These segments are separated by one-way flow valves, with the flow rate denoted as Q, and the flow resistance as R_{V}. At each connection point that links these components, effectively covering both the upstream and downstream ends of each segment, the pressure (P) is explicitly defined. P_{a} and P_{b} represent constant pressure values at the boundaries.

Figure 1:A schematic diagram of the lymphangion-chain model displays the contractile vessel segments, each characterized by its diameter (D). These segments are separated by one-way flow valves, with the flow rate denoted as Q, and the flow resistance as RVR_{V}. At each connection point that links these components, effectively covering both the upstream and downstream ends of each segment, the pressure (P) is explicitly defined. Pa and PbP_{b} represent constant pressure values at the boundaries.

Figure 1 shows a schematic diagram of a chain of contractile vessel segments linked via secondary valves. Both ends are connected to constant pressure boundary conditions, Pa and Pb. The chain is subjected to uniform external pressure, Pe. The pressure-dependent valve resistance is implemented using a sigmoidal function. Flow inside the vessels is modelled using the Poiseuille flow equation. Additionally, the force balance equation for the vessel walls includes both passive and active contractile components. The resultant ordinary differential equation (ODE) system is non-linear. The parameter values and boundary conditions are listed in Table 1 and Figures 4 and 5 in the captions of the primary paper.

The mathematical model was implemented within the open-source, extensible markup language (XML)-based CellML modeling environment. CellML is widely utilized for representing mathematical models in the field of biology, employing ODE and algebraic equations Cuellar et al., 2003. All simulations were executed using OpenCOR (Version 0.7) as the simulation platform Garny & Hunter, 2015. The original CellML file and all the associated code can be accessed by following this link within the PMR: https://models.physiomeproject.org/workspace/b44.

3Computational Simulation

The Bertram2011.cellml file contains the necessary code to execute the simulation. The simulation, conducted over 100 seconds with a time step of 0.0001 seconds, was completed in approximately 4.5 seconds of wall clock time. The CVODE solver was employed, and the simulation was executed on a Windows desktop computer equipped with an Intel Core i7 10th generation CPU. Additionally, the following solver settings were applied: a maximum step of 0.0001 and a total of 5 ×105 steps.

OpenCOR offers the capability to export results in CSV format. When exporting the CSV file for a specific set of parameters and boundary conditions, the following data must be selected from the drop-down menu: time, Q1 to Q5 (Valve_1 to Valve_5), D (Vessels_1 to Vessels_4), Mfunc_1 to Mfunc_4, P11, P12, P21, P22, P31, P32, P41, and P42 values. Subsequently, MATLAB scripts (PlotGenerator_figure_4.m and PlotGenerator_figure_5.m), available in the PMR, were utilized to generate plots after reading the CSV result files. Importantly, the column names of the output CSV file must be changed based on the defined variable names in the PlotGenerator.m files.

The column names should be adjusted as follows:

  • Rename "environment | time (second)" to "T"

  • Rename "Valve_1 | Q1 (ml_per_s)-Valve_5 | Q5 (ml_per_s)" to "Q1-Q5"

  • Rename "Vessel_1 | D (cm)-Vessel_4 | D (cm)" to "D1-D4"

  • Rename "Vessel_1 | Mfunc_1 (dyn_per_cm2)-Vessel_4 | Mfunc_4 (dyn_per_cm2)" to "M1-M4"

  • Remove vessel numbers and units from pressure values (e.g., "Vessel_1 | P11 (dyn_per_cm2)" becomes "P11").

Save the CSV file as "Bertram_2011_Data_figure4/5.csv". Then, running the MATLAB script will produce the correct figures.

4Model results

The simulated results for a chain of four lymphangions against two sets of prescribed pressures, Pa=2275 dyn cm-2 Pb=2875 dyn cm-2 and Pa=2275 dyn cm-2 Pb=2375 dyn cm-2, prescribed in the original work are demonstrated here. Figures 2 and 3 show active tension, diameter variation, pressure values at different locations, and flow rates for each lymphangions. Figures 2 and 3 correspond to Figures 4 and 5, respectively, in the primary paper. We employed the exact parameter values as specified in the primary paper. The parameter values used in the CellML model for Figures 2 and 3 are presented in Table 1.

Reproduced results corresponding to Figure 4 in the primary paper, demonstrate the cyclic pumping by a sequence of four lymphangions operating against a pressure differential of 600 dyn cm ^{-2}. The top panel illustrates the active tension waveform in each lymphangion. The second panel displays the dynamic variations in the diameters of the four lymphangions. The third panel presents the eight distinct time-varying pressures, with horizontal black lines denoting the inlet and outlet pressures P_{a} and P_{b}. The bottom panel illustrates the flow rate through each of the five valves.

Figure 2:Reproduced results corresponding to Figure 4 in the primary paper, demonstrate the cyclic pumping by a sequence of four lymphangions operating against a pressure differential of 600 dyn cm -2. The top panel illustrates the active tension waveform in each lymphangion. The second panel displays the dynamic variations in the diameters of the four lymphangions. The third panel presents the eight distinct time-varying pressures, with horizontal black lines denoting the inlet and outlet pressures Pa and Pb. The bottom panel illustrates the flow rate through each of the five valves.

Reproduced results corresponding to Figure 5 in the primary paper, demonstrate the cyclic pumping by a sequence of four lymphangions operating against a pressure differential of 100 dyn cm ^{-2}. The top panel illustrates the active tension waveform in each lymphangion. The second panel displays the dynamic variations in the diameters of the four lymphangions. The third panel presents the eight distinct time-varying pressures, with horizontal black lines denoting the pressures at points P_{a} and P_{b}. The bottom panel illustrates the flow rate through each of the five valves.

Figure 3:Reproduced results corresponding to Figure 5 in the primary paper, demonstrate the cyclic pumping by a sequence of four lymphangions operating against a pressure differential of 100 dyn cm -2. The top panel illustrates the active tension waveform in each lymphangion. The second panel displays the dynamic variations in the diameters of the four lymphangions. The third panel presents the eight distinct time-varying pressures, with horizontal black lines denoting the pressures at points Pa and Pb. The bottom panel illustrates the flow rate through each of the five valves.

Table 1:Parameter values used in the CellML model.

Parameter valuesFigure 4Figure 5
P_a2275 dyncm-22275 dyncm-2
P_b2875 dyncm-22375 dyncm-2
t_01_Vessel_10.50.5
P_d_01_Vessel_150 dyncm-250 dyncm-2
D_d_Vessel_10.025 cm0.025 cm
t_01_Vessel_211
P_d_01_Vessel_275 dyncm-275 dyncm-2
D_d_Vessel_20.022cm0.022cm
t_01_Vessel_31.51.5
P_d_01_Vessel_3100 dyncm-2100 dyncm-2
D_d_Vessel_30.019 cm0.019 cm
t_01_Vessel_422
P_d_01_Vessel_4125 dyncm-2125 dyncm-2
D_d_Vessel_40.016 cm0.016 cm

5Discussion

This lymphangion chain model is an efficient computational framework for comprehending the pumping performance and nonlinear dynamics of lymphatic vessel diameter fluctuations. There is no overlap between the authors of the primary paper and our Physiome paper. Our goal was to test the reproducibility of the model originally published by Bertram et al. (2011), independent of the original authors’ involvement. Utilizing the reported parameter values, we have effectively replicated and validated the results presented in Figures 4 and 5 of the primary paper. As previously mentioned, subsequent studies have enhanced the intricacy of this model. Consequently, we will be able to curate additional lymphatic models within the PMR, building upon this reproduced model of Bertram et al. (2011).

Acknowledgments

TDJD was supported in part by a New Zealand Ministry of Business, Innovation and Employment 12 Labours project grant. HMR was supported by a Sir Charles Health Research Fellowship from the Health Research Council of New Zealand grant 20/017. SS acknowledges support from the Aotearoa Foundation. AP and PR were supported in part by Health Research Council of New Zealand grants 20/035 and 21/714.

References
  1. Bertram, C. D., Macaskill, C., & Moore, J. E. (2011). Simulation of a chain of collapsible contracting lymphangions with progressive valve closure. Journal of Biomechanical Engineering, 133(1). 10.1115/1.4002799
  2. Bertram, C. D., Macaskill, C., & Jr., J. E. M. (2014). Incorporating measured valve properties into a numerical model of a lymphatic vessel. COMPUTER METHODS IN BIOMECHANICS AND BIOMEDICAL ENGINEERING, 17(14), 1519–1534. 10.1080/10255842.2012.753066
  3. Bertram, C. D., Macaskill, C., Davis, M. J., & Moore, J. E. (2016). Consequences of intravascular lymphatic valve properties: A study of contraction timing in a multi-lymphangion model. In American Journal of Physiology - Heart and Circulatory Physiology (Vol. 310, pp. H847–H860). 10.1152/ajpheart.00669.2015
  4. Bertram, C. D., Macaskill, C., Davis, M. J., & Moore, J. E. (2017). Valve-related modes of pump failure in collecting lymphatics: numerical and experimental investigation. In Biomechanics and Modeling in Mechanobiology (Vol. 16, pp. 1987–2003). 10.1007/s10237-017-0933-3
  5. Davis, M. J., Rahbar, E., Gashev, A. A., Zawieja, D. C., & Moore, J. E. (2011). Determinants of valve gating in collecting lymphatic vessels from rat mesentery. In American Journal of Physiology - Heart and Circulatory Physiology (Vol. 301). 10.1152/ajpheart.00133.2011
Funding and support by:
International Union of Physiological SciencesAuckland Bioengineering InstituteDigital ScienceVPH Institute