1.

Introduction

The lymphatic system in humans propels 2 to 4 liters of fluid from the interstitium back to the venous system daily. The forces behind this process can be divided in two categories: intrinsic and extrinsic forces. Intrinsic forces, also known as the intrinsic pump, are the forces generated within the lymphatic system. They are a major force behind lymph transport and are driven by phasic and tonic contractions. The extrinsic pump is the sum of all forces external to the lymphatic system that help drive lymph from peripheral organs back to the central circulatory system. These forces include pulsation of neighboring blood vessels, motion of adjacent skeletal, cardiac and gastrointestinal muscles, and the suction and compression effects of respiration.¹

A lymph vessel is formed by a series of lymphangions intermittently disrupted by valves to avoid back flow. These lymphangions are the essential functional units that consist of segment of the lymphatic vessel between two valves² and act as miniature hearts that phasically contract to pump lymph from one compartment to the next.³ The contractions within a lymphangion propagate along the lymphatic as waves with an average velocity of 4 to 8 mm/s⁴ and originated from pacemaker cells.^{5, 6} While the question of whether lymphatic pacemaker cells belong to a specific cell population, or if pacemaking is a functional role of the currently most depolarized lymphatic myocytes^{7, 8, 9, 10} has yet to be answered, the crucial role of intrinsic pacemaking in the initiation of the lymphatic contractile events is undoubtful. The understanding of the pacemaking mechanisms in lymphatics is necessary to develop the preventive measures and potential treatments of lymphedema that affects 140 to 250 million people worldwide.¹¹ However, the knowledge of the localization and propagation of electrical events in lymphatic vessels is far from complete. In this respect, the studies to determine the general principles and characteristics of the propagation of lymphatic contractile waves along the lymphatic vessels will provide important information of the role of lymphatic pacemaking in the coordination of effective pumping through collecting lymphatics. Such information is necessary to clarify the potential role of the disturbances of lymphatic electrical events in the pathogenesis of various diseases, which lead to diminished lymph flow and formation of local edema. In the past, the propagation velocity of the contractile waves has been measured using methods based on manual contraction tracking and therefore prone to error.^{4, 5, 12} In this study, we tested a system with a new automated algorithm to detect contraction propagation and to localize the pacemaker sites in the contracting rat mesenteric lymphatics. These techniques will be potentially helpful in developing a better understanding of the origin of the lymphatic contractile waves, characterization of their behavior, and their response to changes in the intramural pressure and lymph flow patterns, that are believed to have a substantial effect on the lymphatic contractile activity.^{1, 5, 6}

2.

Experimental Setup

3.

Methods

When imaging at 300 frames per s, a contractile wave propagating at 8-mm/s traverses less than 27 μm (54 pixels) in one frame period (1/300 of a s). The camera's field of view (FOV) in the direction perpendicular to the vessel was set to be slightly larger than the diameter of the vessel to minimize the camera memory needed. In the other direction, the FOV was set to the maximum capability of the optical setup (576 μm) to image the longest possible section of the vessel. The spatial sampling frequency of this system is 0.5 μm/pixel. This is sufficient to achieve the Nyquist sampling frequency and allows the placement of 21 interrogation windows separated by 27.4 μm with the ability to still detect the propagation of a wave moving at 8 mm/s, which is the high value reported by Zawieja ⁴ This gives a pacemaker localization resolution of approximately 28 μm. If better resolution is required, higher frame rates are possible but at the expense of recording time.

Diameter readings were simultaneously recorded from different interrogation windows along the vessel in the axial direction using a correlation based tracking algorithm.^{14, 15} The time delay between the contractile waveforms from two different axial positions was measured and used to calculate the propagation velocity as shown in Fig. 1.

Fig. 1

Schematic showing the principle of calculating the TD between the contractile waves at two different spots along the same lymphangion. This process was repeated for multiple points along the vessel wall and the TD between consecutive points was calculated leading to a spatial velocity map that we used to localize potential pacemaker sites, the origin of these waves.

Different frame rates ranging from 200 to 500 fps were used to acquire the images depending on the desired recording interval (intervals ranged from 25 to 58 s depending on the frame rate as well as the field of view needed, which varied with the diameter of the vessel under investigation). Multiple recordings of the same site were taken to sample multiple contraction cycles of a given vessel segment.

The algorithm started by searching the diameter tracings for time points of interest. Four types of time points were defined: systole start (SS), systole end (SE), diastole start (DS), and diastole end (DE), however, for the application discussed in this paper, we only used the two of them that are associated with the rapid contractile phase of the pumping activity, SS and SE.

The search criterion consisted of several steps. In the first stage, the derivative of the diameter was computed. This was compared to two thresholds; one negative signaling the start and end of contraction and the other positive indicating the start and end of relaxation. A square wave was created, from the comparison of the derivative to the thresholds, in which the rising and falling edges pointed to the beginning and end of a phase (systole or diastole), respectively. Squares closer to each other than 300 ms were combined together to neglect the effect of small quick oscillations in the diameter tracking due to noise artifacts, such as shaking of the vessel without the use of filtering methods that can result in a loss of details about the fast transition states (SS and DS).

Fine tuning was then used to minimize the error in locating the edges. In particular, a search for better matches in intervals around the edges of the square wave was performed. The algorithm was programmed to start from the edge of the square signal and search in both directions for points of the local derivative that change sign. The square wave was either extended or narrowed to the point found by this search criterion (Fig. 2). The derivative of the square signal leads to two impulse functions signaling the start and the end of the contraction or the dilation.

Fig. 2

Auto-tuning process applied to the diameter tracings corresponding to data shown in Fig. 5. (a) Diameter tracking (—) and square signal (—) delimiting the contraction period. (b) On the left, the first contraction before tuning. On the right, the same first contraction after tuning for better precision in determining the start and end of the systole phase.

Fig. 5

(a) The isolated lymphatic. (b) The contraction propagation velocity average over four consecutive contractions. Zero crossings that correspond to pacemaker locations are displayed as red dots where the central point is the average over the four contractions and the error bars show the range of change in pacemaker location. The zero crossings indicate the presence of two pacemaker sites around 328 and 525 μm. (Color online only.)

The output of the algorithm at this point is a display of the results of time point detection (SS, SE, DS, DE) and the user is allowed to intervene and change any of the automatically detected points to account for any errors in the detection.

To locate pacemaker sites, we started by measuring the time delay (TD) of the contractile waves from consecutive points along the vessel wall known distances apart (Fig. 1) and then calculated the propagation velocities. The sign of the velocity indicates the propagation direction, and the velocity axial profile can be described as in Fig. 3. Figure 3a depicts a simple scenario with one pacemaker site located on the far left-hand side, outside the FOV. In this case, the velocity profile is positive all along the vessel axis. Figure 3b depicts a similar situation with the single pacemaker located on the far right-hand side of the FOV, with the velocity profile negative along the whole length of the vessel. However, in more complex scenarios the sign of the velocity can change, as is seen in two cases in Figs. 3c and 3d. If a vessel FOV has a single pacemaker site near the middle, the contractile waves propagate away from the pacemaker in opposite directions and the propagation velocity changes from negative, left of the pacemaker, to positive right of the pacemaker [type I zero crossing, Fig. 3c]. If a FOV has two pacemakers near opposite sides of the FOV, then two contractile waves are propagating in opposite directions toward each other resulting in a collision site, where the velocity changes from positive to negative [type II zero crossing, Fig. 3d].

Fig. 3

The sign of the contractile waveform propagation velocity carries information on the location of pacemaker sites. Four cases are shown in the schematic above. The black dots denote the pacemaker sites; the FOV is represented by the dashed gray lines. (a) and (b) represent cases where the only pacemaker sites in the vessel are situated outside the FOV. Panel (c) shows the change in contraction propagation velocity sign at a pacemaker site within the FOV (Type I zero crossing). Panel (d) corresponds to the case of a collision site within the FOV (Type II zero crossing).

We also used this algorithm to output the diameter shortening velocity and the contraction frequency. The derivative of the diameter within the limit of the square signal gives the diameter shortening velocity. The separation between two consecutive SS peaks is the contraction period, which is the inverse of the contraction frequency.

3.1.

Simulation of Lymphatic Vessel Contraction

The algorithm was first tested on user generated images simulating the lymphatic system contractile behavior. To the authors’ knowledge, no one has modeled the propagation of the contractile wave in lymphatic vessels. From our observations, a propagating saturated sine wave is a good model to test the performance of our algorithm, its ability to measure contraction propagation, and subsequently locate pacemaker sites. A wave propagating along the x direction can be written as a function of (Kx – ωt) where the propagation velocity is given by V = ω/K.¹⁶ The contractile wave was modeled as a function varying in time and space, propagating outward from a central point, imitating the pacemaker location x ₀ as shown in Eq. 1. The parameters K, ω, and A were varied to change the contraction frequency, propagation velocity, and amplitude as shown in the equations below.

Eq. 1

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \psi (x,t) = \left\{ {\begin{array}{l@{\quad}c} {DD} & {A.{\rm Sin}(K\left| {x - x_0 } \right| - \omega.t + \varphi) > DD} \\ {SD} & {A.{\rm Sin}(K\left| {x - x_0 } \right| - \omega.t + \varphi) < SD} \\ {A.{\rm Sin}(K\left| {x - x_0 } \right| - \omega.t + \varphi)} & {\rm otherwise} \\ \end{array}} \right.\vspace*{-3pt} \end{equation}\end{document} $$\psi (x,t)=\left\{\begin{array}{cc}\hfill DD& \hfill A.\mathrm{Sin}(K\left|x-x_0\right|-\omega .t+\varphi )>DD\\ \hfill SD& \hfill A.\mathrm{Sin}(K\left|x-x_0\right|-\omega .t+\varphi )<SD\\ \hfill A.\mathrm{Sin}(K\left|x-x_0\right|-\omega .t+\varphi )& \hfill \mathrm{otherwise}\end{array}\right.$$

x ₀: pacemaker location; ω = 2πf: contractions' angular frequency; f: contractions' frequency; K: wave number K = ω/V; V: contractions’ propagation velocity; A: contractions’ amplitude; SD: systole minimum diameter; and DD: diastole maximum diameter.

1 shows the simulated images at different times illustrating the propagation of the contractile wave.

Video 1

Simulated lymphatic in different states with an end diastolic diameter of 128.5 μm and end systolic diameter of 91 μm. (a) Vessel in diastole, (b) start of the contraction at the pacemaker site at t = t ₀, (c) propagation of the contraction away from the pacemaker site at t = t ₀+8 ms, and (d) end of systole. (QuickTime, 2 MB).

10.1117/1.3544512.1

Imaging noise in a CCD has a Poisson distribution, which in brightfield microscopy due to the relatively high levels of light, can be approximated as additive Gaussian noise.¹⁷ To imitate that potential noise source, Gaussian noise was added to the images which mimics the challenges presented in the real case to the tracking algorithm. The diameter was recorded at 10 different spots along the simulated vessel wall and the propagation velocity was calculated using our automated algorithm manually, as shown in 1. The manual contractile wave tracking from the diameter tracings showed no difference from the automated tracking shown in Fig. 4.

Fig. 4

Results from the simulation images show the change in the contractile wave propagation direction at the pacemaker site. The simulated contractile wave propagation velocity spatial map is shown in (-). These simulated images were processed by the developed algorithm and the output is shown in (*). The propagation velocity was also calculated manually from the diameter tracings using the algorithm discussed by Zawieja (Ref. 4) and the results showed no systematic differences.

To verify the results and determine the source of the small deviations, we calculated the contraction propagation velocity manually from the diameter tracings and we received the same results as the automated algorithm, proving that the source of the error is the inherent errors of the diameter tracking algorithms due to pixelation of the vessel wall in the simulation. The mean difference between the contraction propagation velocity calculated by the algorithm and that implemented in the simulation images was 0.26 mm/s with a mean squared error (MSE) = 0.659 mm/s.

4.

Results and Discussion

4.1.

Isolated Lymphatic Experiments

The developed system was used to determine the pacemaker cells location(s), contraction propagation velocity, contraction diameter shortening velocity, and contraction frequency. Multiple pacemaker locations were located in a section of an isolated mesenteric microlymphatic. Figure 5 shows the average contraction propagation velocity of four consecutive contractions. The contraction propagation velocity changes sign at four different locations, two of which are changing from propagation in the negative direction to the positive direction (type I) indicating the presence of pacemaker cells located around 328 and 525 μm. The propagation velocity changes from positive to negative (type II) in the other two zero crossings indicating that these points are collision sites located in between two pacemakers. The first type II zero crossing is located at 415, 110 μm to the left of the pacemaker site at 525 and 87 μm to the right of the other pacemaker site. The second type II zero crossing is located at 248.5 μm which indicates the presence of a pacemaker on its left side.

To verify the repeatability of the method under various physiologic conditions, we investigated the presence of pacemakers in a different section of the same vessel under different pressure and flow conditions increasing the contractions frequency from 7 contractions per min in a previous data set to 21 contractions per min covering most of the physiologic range for microlymphatic vessels. In addition, this section also showed the presence of one pacemaker site as shown in Fig. 6. This data set includes 15 consecutive contractions compared to only 4 in the previous one (Table 1, set 2). Figure 5 averages the contraction propagation velocity over four contraction cycles. However, different contractions might be initiated from different pacemaker sites. We conducted a temporal study of the same vessel, analyzing all four contractions independently and calculating the propagation velocity for each contraction period. The pacemaker site at 328 μm seems to be active during all four contractions. The second pacemaker site at 525 μm changes its location to 397, 525, 450, and 493 μm during contractions 1, 2, 3, and 4, respectively. Figure 5 shows this variation as red error bars with the central point being the average of all four pacemaker positions during the four different contractions.

Fig. 6

(a) 2 of the image of the lymph vessel under investigation. (b) The average propagation velocity over 15 consecutive contractions. On average, there is only one zero crossing around 250 m indicating the presence of a possible pacemaker site. The absolute values of propagation velocities changed between 0 and 27.6 mm/s and averaged 8.1 mm/s. (QuickTime, 2.9 MB)

10.1117/1.3544512.2

Table 1

Summary of results from vessel 1 and a comparison showing the effect of user intervention on the results.

Data set		Average frequency (contractions/min)	Contraction velocity (μm/s)	Maximum diameter (μm)	Minimum diameter (μm)	Average diameter (μm)
Set 1	Before user intervention	7.794	20.1140	179	53	129.46
	After user intervention	7.428	20.4363	179	53	129.46
Set 2	Before user intervention	21.4335	28.583	167	82.5	126.85
	After user intervention	21.6790	27.1844	167	82.5	126.85

We did the same temporal analysis for the data set shown in Fig. 6. Six of the 15 contractions showed the presence of two pacemaker sites separated on average by 167 μm and centered at 290 μm. The rest of the contractions showed only one pacemaker site at the axial position 251.5 μm with a standard deviation of 68.27 μm.

To show the repeatability of our method, data was collected from two additional vessels and the results are shown in Figs. 7 and 8. The vessel in Fig. 7 showed a pacemaker site around 540 μm. This pacemaker site was stable through all three contractions recorded. The propagation velocities recorded on this vessel ranged from 0 to 10.37 mm/s. This falls in the physiologic range reported by Zawieja ⁴

Fig. 7

Average propagation velocity over three consecutive contractions. A pacemaker site is detected around 540 μm.

Fig. 8

Data from the third vessel showed the presence of two neighboring pacemakers around 450 and 577 μm. The propagation velocity shown in (b) was measured from four consecutive contractions.

Figure 8 shows the data from the third vessel investigated. Two pacemaker sites were detected around 450 and 577 μm. These pacemakers are close to each other and they both were present around the same area in all four contractions recorded.

5.

Discussion

Our system is based on tracking the mechanical behavior of lymphatics to deduce the origin of the electrical activity. This is validated by the fact that phasic contractions in lymphatics are driven by Action Potentials^{18, 19, 20, 21} and collecting lymphatics contain one to three layers of smooth muscle cells, providing a complete coverage of the vessel which helps in the propagation of this action potential from a pacemaking site to all neighboring muscle cells^{21, 22} resulting in the propagation of the phasic contractions. Direct readings of electrical activity using microelectrodes cannot discriminate between the pacemaker generated potential and the propagated action potential from neighboring muscle cells,²¹ thus making it difficult to localize pacemaker sites. The measurement of phasic contractions propagation leads to the action potential propagation. The contractions starting point should correspond to the pacemaking site.

Although the algorithm has a self-correction mechanism as shown in Fig. 2, we left the option of adjusting the location of the detected peaks (SS, SE, DS, DE) to the user. Table 1 shows the effect of user adjustments on the final results from the first vessel and this effect was minimal due to the fine tuning performed automatically by the algorithm minimizing errors. These correction techniques, autotuning and possibly user intervention, will have more of an effect in our future in vivo studies because of the higher background noise and lower contrast in these experiments.

The results we showed in this paper are taken from isolated mesenteric microlymphatic vessels. Using an animal preparation protocol similar to the one described by Dixon,²³ the same technique can be applied in vivo with the added difficulty of finding clear sections of lymphatic vessels to detect contraction propagation. In that regard, the use of high speed imaging gave us the advantage of being able to detect contraction propagation within relatively short sections compared to the previous manual method.⁴ In the future, this system will allow us to conduct multiple in vitro and in vivo experiments, such as studying pacemaker activity under different pressure gradients and monitor flow and shear stress that are believed to be the main modulators of the contractile activity in the lymphatic system.^{1, 24} Coupling this system with flow measurements will help us separate the contractile wave propagation due to the electrical activity of pacemaker cells and to the propagation of pressure waves.

Compared to high speed video microscopy, there are many other powerful imaging techniques. However, we demonstrate that this modality has enough resolution, contrast, and speed to locate pacemaker regions in microlymphatics. In addition, the simplicity of this system allows its use with almost any commercial microscope. Furthermore, this technique has proven to be the most successful in measuring flow and shear stress in lymphatics, as discussed by Dixon ^{23, 24} in conjunction with the determination of contraction propagation characteristics, which as mentioned earlier, is one of our future goals for this study. Our system outputs a variety of parameters other than pacemaker sites such as the contraction frequency and the contraction shortening velocity. Table 1 summarizes these results for the two data sets shown previously.

The results we presented from all three vessels showed the presence of at least one pacemaker site in each vessel. Beside the data in Fig. 5, all the data sets had a pacemaker site that was very stable throughout the imaging period (3 to 4 contractions). The data in Fig. 5 represents a special case, where one of the zero crossings was changing its spatial location from one contraction to the other. One explanation is that this crossing is not an actual pacemaking site but a discontinuity in muscle cell coverage. This is possible since this data was taken next to the valve region (the valve can be seen in the left part of the vessel image in the top panel of Fig. 5) where the muscle coverage is believed to be incomplete. This discontinuity in muscle cell coverage causes an interruption of electrical conductivity and the contractions on both sides of the discontinuity are not completely synchronized, causing the delay between the points on both sides to fluctuate.

6.

Conclusions

It is a well-accepted paradigm that lymph vessels possess pacemaker cells that initiate contractile waves. However, no one was able to locate the sites of pacemaking precisely and study the propagation of electrical activity, generated by them, their distribution along the lymphatics or their response to various imposed conditions (transmural pressure gradient, shear stress, etc.). The system we presented in this manuscript is the first semi-automated tool (with minimal subjective user modifications) with the ability to detect pacemaker sites in lymphatics. The presented results proved the high potential of our new tool to investigate self-initiated contractions in lymph vessels and pacemaker sites were located in three different isolated rat mesenteric lymphatics. Also, two sections of the same lymphatic vessel were investigated under different physiologic conditions and pacemaker sites were located. The algorithm was first tested on simulated images and the results showed no difference from manual tracking of the contractile wave. Ultimately, this system will be used to better understand the nature of contractile activity in lymphatics and study its regional and local variability, as well as principles of its regulation.