The lymphatic system is responsible for transporting interstitial fluid back to the bloodstream, but unlike the cardiovascular system, lacks a centralized pump-the heart–to drive flow. Instead, each collecting lymphatic vessel can individually contract and dilate producing unidirectional flow enforced by intraluminal check valves. Due to the large number and spatial distribution of such pumps, high-level coordination would be unwieldy. This leads to the question of how each segment of lymphatic vessel responds to local signals that can contribute to the coordination of pumping on a network basis. Beginning with elementary fluid mechanics and known cellular behaviors, we show that two complementary oscillators emerge from i) mechanical stretch with calcium ion transport and ii) fluid shear stress induced nitric oxide production (NO). Using numerical simulation and linear stability analysis we show that the newly identified shear-NO oscillator shares similarities with the well-known Van der Pol oscillator, but has unique characteristics. Depending on the operating conditions, the shear-NO process may i) be inherently stable, ii) oscillate spontaneously in response to random disturbances or iii) synchronize with weak periodic stimuli. When the complementary shear-driven and stretch-driven oscillators interact, either may dominate, producing a rich family of behaviors similar to those observed in vivo.
To maintain fluid homeostasis, interstitial fluid drains into the lymphatic system through initial lymphatic vessels that carry it to the collecting lymphatic vessels. The collecting lymphatic vessels transport the fluid (known as lymph) both passively and actively to lymph nodes and back to the systemic blood circulation. Collecting lymphatic vessels are surrounded by specialized lymphatic muscle cells (LMCs) [1] and sub-divided by valve structures that define individual segments called lymphangions (Fig 1) [2]. Lymphangions serve as both pumps and conduits. In contrast to the blood circulation, where a single pump drives flow through relatively passive conduits, each lymphangion has the ability to pump lymph through the converging network to lymph nodes and eventually to the thoracic duct. Pumping occurs when expansions in radius draw fluid into the upstream end of the lymphangion and then expel it downstream during a contraction. Directional flow is enforced by intraluminal valves, which favor flow toward the thoracic duct. Lymphatic vessel contractions are triggered when cytosolic Ca2+ entering from intravascular stores and outside the cell surpasses a threshold concentration in the cytoplasm of the LMC, resulting in actin and myosin cross-bridging within the LMCs [3]. The contraction phase ends as transmembrane pumps restore cytoplasmic Ca2+ concentration to equilibrium allowing actin-myosin binding to relax, and the trans-wall pressure and the passive elastic properties of the wall to reopen the vessel. The effects of Ca2+ on LMC contraction are moderated by endothelial-derived relaxation factors (EDRFs) that act as potent dilators of lymphatic and blood vessels when produced by the vessel-lining lymphatic endothelial cells (LECs) in response to dynamic fluid shear stresses. The best known EDRF is nitric oxide although others such as histamine have been shown to be important [4,5]. For notational simplicity we represent the entire class of EDRFs herein as NO. The NO and Ca2+ levels are both subject to mechanical regulation; Ca2+ can enter the cell through stretch-activated ion channels [6,7], and NO is produced by LECs when they are exposed to increased fluid shear stress [8]. Although rhythmic contractions can be produced by purely chemical oscillations in Ca2+ within the LMCs [9,10], it is likely that feedback regulation is necessary for robust homeostasis. Indeed, we previously used a relatively complex numerical simulation of lymphatic pumping to demonstrate that a wide spectrum of oscillatory behaviors is possible, and that the behavior is very sensitive to local levels of stretch and stress [11]. Our present aim is to reduce the complexity of our previous model to examine how the observed oscillations arise from the integration of simple mechanical and chemical processes within a physiological system. Once simplified, we employ tools such as linear stability analysis to identify key parameter groups that determine the qualitative dynamic behaviors of such systems. Linear stability analysis seeks to determine which parameter values cause a small disturbance to grow rapidly from an initial state, or alternatively decay back to equilibrium. In addition, our approach allows us to show how oscillations can arise from the interactions between mechanical and chemical processes that lack intrinsic oscillators when considered separately. Here we develop generic formulations of the mechano-chemical processes in the muscular lymphatic vessel wall based on Ca2+ and NO signaling, then explore the dynamics of each chemical species while holding the effects of the other constant. Finally we examine the behavior of the fully coupled system. Linear stability analysis reveals a new class of oscillator arising from the dynamics of shear and NO that can act alone or in concert with the better recognized Ca2+ dynamics. The most remarkable feature of the shear-NO mechanism is its ability to offer distributed control of the pumping process, which is essential for managing a decentralized network of pumps and conduits.
Our model is based on a single lymphangion (Fig 1) bounded at each end by one-way valves. The radius of the lymphangion is governed by the radial forces which are determined by the contractile (Ca2+) and dilatory (NO) signaling molecules. Neglecting inertial effects, the radial forces on the vessel wall balance as DdRdt=−E (R) −F (R, CCa, CNO) +AR (t) (1) where the left hand side represents rate dependent effects with D incorporating the visco-elastic material properties of the vessel wall as well as viscous losses in the flow and lags in the transduction of concentrations into force, E is a restoring force-including elastic forces and vessel “tone” -imparted by the material properties of the vessel wall and F is a dynamic inward acting contractile force produced by muscle cells that surround the vessel. To a first approximation, we assume that the concentration of Ca2+ is transduced into a contractile force as F=FCaCCa/ (1+αCNO) where α scales the possible desensitizing effect of NO [12]. The activation term AR can include a steady component from the mean transmural pressure difference pm that influences the baseline radius as well as extrinsic disturbances to the radius from the surrounding tissue and adjacent lymphangions. The restoring force is typically highly nonlinear[13–15]. Here we adopt the form E (R) =AeaR−poffset with stiffening coefficient a, scaling coefficient A, and offset pressure poffset selected to give a good fit of Shirasawa and Benoit (see the third figure of reference [15]) at typical operating pressures. In our numerical simulations we retain the full nonlinear form of E (R), but for the stability analysis that follows we linearize the elastic force near an equilibrium radius R1 in the absence of dynamic increments to Ca or NO as E (R) ≈E0+E1 (R−R1) where E0 is the elastic force at equilibrium and a Taylor series expansion near equilibrium yields E1=AaeaR1. We find the equilibrium radius by solving 0=DdRdt=− (AeaR1−poffset) −FCaSCa0/KCa+pm for R1 where pm is the mean transmural pressure. Given the stiffening behavior of the wall (E1∝eaR1), we expect that appropriate values of E1 will be larger at higher mean transmural pressures where the equilibrium radius will be somewhat larger. The concentrations of the signaling molecules (i ∈ {Ca, NO}) are governed by the generic conservation law dCidt=−Ki (CCa, CNO) +Si (R, R˙, CCa, CNO) +Ai (t) (2) where all concentrations are taken to be dimensionless ratios relative to a suitable reference, Ki is clearance of the signaling species through chemical reaction, transmembrane ion pumps and advective-diffusive transport, Si is a dynamic source term for the signaling molecule and Ai is an additional source term that can include the effects of imposed flow from upstream fluid pressure, inflammation, pace-making signals from adjacent cells, neural signaling, random disturbances, etc. Since our focus is on the interactions between Ca2+ and NO, we introduce a minimal representation of the Ca2+ dynamics rather than a fully detailed model of Ca2+ oscillations as may be found in the literature [10,16]. We retain the following features: i) at rest, Ca2+ is at a low concentration in the cytoplasm of LMCs; ii) a contraction is initiated when Ca2+ is rapidly admitted to the cytoplasm through ion-selective channels thereby triggering cross bridge formation between actin and myosin chains creating a contractile force [17]; iii) relaxation of LMCs coincides with a drop in cytoplasmic Ca2+ concentration due to a drop in the rate of influx and the restoration of baseline conditions by ion pumps in the cell membrane and sarcoplasmic reticulum; and iv) the LMC is refractory to a new contraction cycle until Ca2+ levels have returned to near equilibrium. As the Ca2+ levels approach their threshold level, we hypothesize that the membrane acquires sensitivity to small perturbations. Furthermore, the sensitivity is enhanced when the membrane is stretched to a larger radius. This models stretch-sensitive ion channels found in LMCs [6,18]. Each step in the process has the potential for modulation by NO. Alternatively, each form of modulation by NO can be disabled to demonstrate behaviors that have been observed in experimental preparations, for example after removal of LECs (which produce NO) or the genetic or pharmacological suppressions of NO [19–21]. We mathematically express the release of Ca2+ into the cytoplasm from intracellular stores and the extracellular fluid as the sum of a steady source Sca0 needed to maintain the baseline Ca2+ concentration and a transient component that is sufficiently rapid to be modeled as an impulse function δ (t) where t is the time since the Ca2+ concentration most recently passed the threshold necessary to trigger another contraction CCaThresh. This is expressed as, SCa=SCa0+SCa1δ (t) (1+γCNO) (3) where SCa1/ (1 + γCNO) is the magnitude of a bolus of Ca2+ with a possible reduction due to NO that is scaled by γ. We model the clearance of Ca2+ from the cytoplasm with KCa=KCa (1+βCNO) (CCa−CCaThresh) (4) where KCa is a rate constant and β scales the possible enhancement of Ca2+ clearance attributed to NO [12,16]. At high concentrations of Ca2+ the clearance rate may be limited by the membrane pump capacity, but near the threshold required to trigger a contraction we assume clearance rates proportional to the concentration increment. The threshold itself may include a random component that we incorporate into the activation term. The linearized form of the model for a constant level of NO can now be written as DdRdt=−E1R−FCaCCa+AR (t) (5) and dCCadt=−KCaCCa+SCa1δ (t) +ACa (t) (6) where the constants have been absorbed into the activation terms so that the radius and concentration now represent the increments from baseline values. The parameters SCa1 and KCa now include the adjustments due to NO introduced in Eqs 3 and 4. In the previous section we developed the model so that it reproduces Ca2+ induced contractions. We next considered how NO, created when LECs experience increased shear stress, can modulate contractions when it diffuses rapidly into adjacent LMCs. A suitable form for the NO source term can be obtained by considering steady laminar flow in a circular tube with negligible inertia [22]. Conservation of mass in a tube of time-varying radius requires that ∂Q∂z=−2πRdRdt which when integrated with respect to z along a vessel yields Q (z) =−2πRdRdtz+C1 where the constant of integration depends on the end conditions for the lymphangion. When the segment is contracting dRdt<0 and the upstream valve is closed [Q (0) = 0] we have Q (z) =−2πRdRdtz. Alternatively, if the segment is expanding dRdt>0 and the downstream valve is closed (Q (L) = 0), we obtain Q (z) =2πRdRdt (L−z). We can express the mean flow along the length more compactly as Q¯=πR|dRdt|L. When additional flow Q0 is imposed on the segment by an axial pressure gradient we have Q¯=πR|dRdt|L+Q0. Approximating the velocity profile with that of steady laminar flow with negligible inertia [22] we relate the mean shear stress τ to the flow rate by τ=4μQ¯πR3. The mean shear stress along the segment due to dynamic changes during contraction or expansion is therefore approximately τ=4μLR2|dRdt|+4μQ0πR3. Recent studies show that valves in collecting lymphatic vessels are biased toward the open condition [23], but the simplification employed here allows us to study the basic stability of the system, at the possible expense of some accuracy in the predictions of pumping efficiency. There may be levels of shear stress below which NO production is negligible and above which NO production saturates at a maximum, but here we linearize the transduction of shear stress into the production of NO in an intermediate range to yield SNO=SNOR2|dRdt|+SNO0 (7) where SNO has absorbed the remaining constants in the shear stress expression and SNO0 represents NO released due to the through-flow term Q0 or chronic sources of NO such as might arise during inflammation. The source term SNO0 can be time varying, but arises from the local environment of the lymphangion and mathematically acts as an input to our model of a single lymphangion rather than as an interaction within the system itself. SNO0 therefore can serve as an external trigger to the system or as a steady offset, but does not directly impact the dynamics of an individual lymphangion, except by parametrically (rather than dynamically) changing the equilibrium radius. We can examine the effects of fluid viscosity on the pressure by using the same set of assumptions. The pressure will vary due to viscous flow effects according to ∂p∂z=−8μQπR4. When contracting we have ∂p∂z=16μR3dRdtz, which when integrated along the length gives p (z) =16μR3dRdtz2+p (0). Averaging over the length of the segment yields p¯=16μL23R3dRdt+p (0) where the first term gives the magnitude of the pressure decrement (or increment for vessel expansion) due to flow induced by the contraction of a single lymphangion. We see that the pressure increment due to flow induced by the single lymphangion also multiplies dRdt, so it can be absorbed into the overall damping term D. For typical vessel sizes, we find that the lag due to viscosity is orders of magnitude smaller than that from chemical and mechanical lags which are on the order of one second. NO does not produce a true outward force. However, it is conceptually equivalent to consider an effective force produced by NO that has the effect of countering FCa and the elastic effects. Mathematically, for small αCNO, we can write this as F=FCaCCa01+αCNO≈FCaCCa0 (1−αCNO) (8) By defining FNO ≡ FCaCCa0α, we can write FNO=−FNOCNO (9) And the net force from Eq 8 becomes F=FCaCCa0−FNOCNO (10) where the net contractile force is decomposed into a positive term set by the baseline Ca2+ levels and a negative term that represents how the Ca2+ levels are modulated by NO. Thus, the NO-dependent term is not a true outward force, but arises mathematically from a reduction in the Ca2+-dependent contractile forces. The parameter values used in the simulations that follow are given in Table 1. The parameter values were based on experimental data where possible, but were chosen to demonstrate a wide range of mathematical behaviors for the system rather than to mimic a particular experimental data set in detail. As a representative example, we show simulations based on measurements in rats [24] which offer data relating Ca2+ concentrations to lymph vessel diameter and contractile tension. Our own experiments discussed later [20] were done on mice which have smaller collecting lymphatic vessels than rats. Specific parameters governing the effects of NO are difficult to estimate, but fortunately may not be necessary here. As will be shown in the results section, we require only estimates of combinations of parameters such as SNO and FCaCCa0α, rather than values for each parameter individually. Unlike the geometrically-detailed continuum model of lymphatic NO transport in Wilson et al [25] that includes shear-induced production and clearance by diffusion, convection and reaction, our present model employs averages over a single lymphangion and combines the sensitivity to shear stress with the rate of production of NO. To that end, we employ parameter values that yield diameter changes due to NO on the order of 10% as observed in [20,21]. Moreover, we expect the effects of NO to be rapid. NO is released by endothelial cells about 2 seconds or less after increases in shear stress as observed previously [26–29]. Lymphatic vessels can be expected to dilate faster than blood vessels [30] because their muscle cells contain more rapid-acting contractile proteins than those of blood vessels [31]. We take the clearance of NO to be similar to, but somewhat faster than, that of Ca2+ [32]. Parameters such as the contractility, NO production and the mechanical stiffness appear to depend on anatomical location, species and age [13,28,33–36], suggesting that the full range of possibilities realizable in vivo awaits further investigation.
In the absence of dynamic activation, Eq 6 implies the Ca2+ concentration during each contraction cycle will decay as CCa (t) =SCa1e−KCat. Using this as an input to Eq 5, we find that the radius varies from its baseline value during each contraction cycle as ΔR (t) =FCaSCa1E1KCa (e−t/tCa−e−t/tmech) (tmech−tCa) (11) where we see that the return to equilibrium depends on two characteristic times, one set by the rate of Ca2+ clearance tCa = 1/KCa and the other by the mechanical lag tmech = D/E1 which can include the lag between the concentration increase and force production. As a reference time scale we have selected tCa = 1/KCa = 1 s which is in the range observed by Shirasawa and Benoit [15]. They observed a similar lag between the rise in Ca2+ concentration and the peak force generation. Here we use this lag as an estimate of tmech which we take to incorporate the visco-elastic and chemical-mechanical transduction lags. While both time constants contribute to the overall response, the slower of the two characteristic times gives the dominant time constant tc that determines the return to equilibrium. Experimental observations of the magnitude of radius change as a function of pressure show that it decreases with increased internal pressure [21]. This phenomenon is reproduced by our model as the vessel wall stiffens (larger E1 in the denominator) at greater radius. The amplitude may be further modified if the myosin cross bridging is length dependent as seen in skeletal muscle [36]. The frequency of contractions at constant NO levels is set by the interplay between the characteristic time for calcium tCa and the magnitude of the random activation term. In addition to the steady source of Ca2+ that establishes the vascular tone, we include a random component ACaRand (R, t) with zero mean and a standard deviation σ that can be applied to either the concentration itself or to the threshold level at which a new release of Ca2+ is triggered. A higher radius leads to more stretch in the LMC membrane and therefore greater sensitivity of ion channels. This can be modeled by increasing the noise level (for example let σ ∝ pm). This will lead to a higher frequency at larger radii as found experimentally [21]. Exponential decay of Ca2+ near equilibrium leads to a latency period between contractions that varies as T = tc log (Cmax/σ) where Cmax is the magnitude of the Ca2+ increment from baseline. We further investigated random activation with the aid of numerical simulations implemented with the Euler-Maruyama method [37], which properly scales the computational time step with the standard deviation of the noise (Fig 2). The simulation presented in Fig 2d–2f results from a higher mean pressure than Fig 2a–2c. Thus, the baseline radius is larger in Fig 2d–2f, which in turn yields a stiffer wall (larger E1) leading to a smaller mechanical time constant (smaller tmech) and reduced amplitude for change in the radius. We note that our simple model of stretch-Ca2+ dynamics replicates important features of the contraction cycles observed in vivo [20] where contractions are generally similar to one another in magnitude and duration, but may be separated by inconsistent periods of latency. In Fig 2, we see that the simulated contractions are nearly identical to each other (that is, the trajectories nearly retrace one another in the phase portraits shown Fig 2c and 2f), but occur on inconsistent intervals. Even though the intervals between contractions are not perfectly uniform, they are well estimated by tc log (Cmax/σ). We also see that increased transmural pressure can reduce the interval between contractions by stiffening the wall (pm ↑⇒ R1 ↑⇒ E1 ↑⇒ tmech ↓⇒ tc ↓⇒ T↓) and also by increasing the sensitivity of the Ca2+ channels by stretching the vessel wall (pm ↑⇒ σ ↑⇒ T ↓) yielding higher frequency contractions (compare Fig 2b and 2e). We also find that our model of the stretch-Ca2+ process readily synchronizes when we impose extrinsic rhythmic pace-making since only small variations in the Ca2+ concentration relative to the threshold level are needed to initiate the next contraction cycle (Fig 3). Such small variations in Ca2+ concentration can be readily introduced by diffusion or voltage signals from adjacent LMCs. Alternatively, the vessel may be locally stretched by lymph arriving from upstream, which can also trigger a local contraction. In this way, neighboring LMCs can synchronize contractions to coordinate flow along a series of lymphangions throughout a connected network of collecting lymphatic vessels [38,39]. The conditions for oscillations in radius to arise near baseline Ca2+ levels in the absence of sharp spikes in Ca2+ as considered in the previous section are available from linear stability analysis of the shear-NO process near a point R1 which yields [R˙CNO. ]=[−E1DFNOD−SNOE1DR12 sgn (R˙) SNOFNODR12−KNO][RCNO]+[inputs] (12) where the inputs include all extrinsic disturbances from the adjacent lymphangions and surrounding tissue. We treat small variations in R parametrically so that the dynamics of the system may be characterized by the eigenvalues of the Jacobian matrix [40] which are roots of the characteristic polynomial: λ2+ (E1D+KNO−sgn (R˙) SNOFNODR12) λ+E1KNOD=0 (13) Since all of the coefficients are positive, stability requires only that the second term be positive. Thus the system is always stable during contraction (R˙<0). However during dilation (R˙>0) the second term can be positive or negative which allows the system to switch between stability and instability (Fig 4). Herein is the key feature from which NO can induce spontaneous oscillations in the radius without the sharp spikes in Ca2+ concentration described in Eq 6. If the radius is large enough so that E1/D + KNO > SNOFNO/DR2, then the fixed point is inherently stable. If instead, the radius is small enough that E1/D + KNO < SNOFNO/DR2 then the radius will unstably increase when perturbed. The instability arises because a slight increase in radius (from point 0 on Fig 5) pulls fluid into the lymphangion, increasing shear and temporarily creating a runaway effect wherein more NO is released from the LEC further increasing the radius and drawing in still more fluid (upper branch from point 0 to point 2 on Fig 5b). The instability persists until the radius becomes large enough at point 2 that the shear stresses begin to drop because more cross-sectional area is available for lymph flow. Thereafter the release of NO occurs more slowly than its degradation so that the system can return stably to equilibrium along the lower branch of the trajectory from point 2 to 0. Mathematically, the unstable increase in radius persists until the sign of R˙ changes at point 2. A change in the sign of R˙ does not require that the eigenvalues move to the left half of the complex plane at point B on Fig 4 as required for inherent stability, but rather requires only that the radius increase sufficiently to move the eigenvalues off of the real axis beyond point A, thus permitting at least a partial cycle of oscillation that includes a time at which R˙=0. As the vessel begins to contract, the sign of SNOFNO/DR2 changes at the R-nullcline where R˙=0, leading to an unconditionally stable return to the original radius. The time scale for contraction is approximated by tc ≈ tmech + tNO + tFNO where the three contributions arise from mechanical lag tmech as before, the clearance of NO tNO = 1/KNO, and the rate of force modulation by NO tFNO = SNOFNO/E1KNOR2. To a similar degree of approximation, the instability of the NO-shear dynamics requires tmech + tNO ≤ tFNO. In other words, the change in force elicited by shear stress must persist longer than processes that tend to dissipate its effects. Exact algebraic expressions for the eigenvalues may be employed if desired, but this approximation captures the key dependencies. See Table 2. The NO cycle can be generalized into a controllable and synchronizable oscillator. Fig 6 shows the behavior of the NO cycle in response to small random disturbances. During the stable contraction process, the vessel remains refractory to disturbances until close enough to equilibrium for a random disturbance to trigger another cycle, much as we found with the stretch-Ca process. Here the period of the NO-induced oscillations varies as tc log (Cmax/σ) as before but tc and Cmax now refer to NO rather than Ca2+. As with Ca2+, a relatively quiet environment or reduced sensitivity to disturbance will elicit longer latency periods between cycles but will not significantly change the shape of the shear-NO cycle. The NO dynamics can also readily synchronize with externally-imposed, small-amplitude sinusoids (Figs 7 and 8). Fig 7b and 7c show how the radius oscillates at precisely the input frequency for frequencies reasonably close to the response when noise triggered (Fig 6). However, when the input frequency is too high (Fig 7a) or too low (Fig 7d) synchronization occurs, but at half or double the input frequency, respectively. Fig 8 shows a parametric study of synchronization over a wide range of input frequencies and amplitudes where we find that synchronization can include a variety of integer ratios between input and output frequencies as explained further in the Discussion. Having explored the system dynamics when Ca2+ and NO are taken to be constant relative to each other we now consider their combined, dynamic effects. Our model includes three possible interactions reported in the literature: NO may i) desensitize the LMCs to Ca2+ as modeled by Eq 10, ii) modify the availability of Ca2+ by 1/ (1 + γCNO) or iii) speed clearance of Ca2+ by 1 + βCNO [41,42]. Our simulations (Fig 9) show that the dynamic effects of NO are most pronounced when the shear-NO dynamics are unstable. When the shear-NO dynamics are unstable, the radius can overshoot the nominal radius before or after a Ca2+-induced contraction, yielding oscillations in radius that are more symmetrical about equilibrium than when shear-NO is stable. At marginal stability (Fig 9d), the NO concentration rings at a frequency determined by the point where the eigenvalues of the shear-NO oscillator cross the imaginary axis f = (E1KNO/D) 1/2/2π. At larger radii, the shear-NO mechanism is inherently stable, but can still reduce the magnitude of the oscillations driven by the stretch-Ca2+ process. This process is important in the presence of an assisting pressure gradient because the dilation induced by the forced flow can put the vessel into the range of radii where the shear-NO mechanism can inhibit contractions that would otherwise tend to restrict free flow through the vessel. The overall frequency is set by a complex interplay of stretch-Ca2+ and shear-NO mechanisms, but will typically be dominated by the faster of the two processes. Long latency intervals between Ca2+-induced contractions can permit NO to produce an unstable dilation, whereas, short intervals due to Ca2+ can suppress the autonomous oscillations possible through the NO mechanism. Interestingly, the published clearance rates for Ca2+ and NO cover a wide enough range that either possibility exists in vivo [15,32]. Experimental observations of diameter in vivo show cycles consistent with the model predictions (Fig 10) (data from [20]). In the absence of direct measurements of concentrations, we employ an alternative phase portrait of diameter plotted against the rate of change of diameter. Fig 10a and 10b are from a wild-type mouse in which Ca2+ and NO effects can operate normally. Here, we observe complex oscillations that include both rapid contractions and occasional strong dilations above the baseline diameter as expected from the shear-NO mechanism. In contrast, when the NO effects have been genetically deleted in eNOS-/- mice in Fig 10c and 10d, we see wave forms that are nearly identical to each other but dominated by contraction with the dilatory effects of NO appearing to be substantially weakened. In all cases, we see cycles occurring on irregular intervals as we expect from noise-triggered oscillators.
While the correspondence between the experiments and the model is encouraging, we should not expect a model of an isolated lymphangion to reproduce all features of a vessel in an intact, in vivo vascular network. For example, the effects of flow introduced from upstream or disturbances from surrounding tissue are inputs present in the animal, but are not included in the modeled dynamics. We have also not yet included effects due to nonlinear valve efficiencies or the bias of the check valves toward the open position [43]. Nonetheless, phase portraits, such as those newly employed here, promise to assist further study of the nonlinear dynamics that govern vascular oscillations. While our results await further experimental validation and improved estimates of key parameter values, it is interesting to consider the newly identified shear-NO oscillator more generically from a nonlinear dynamics perspective. The shear-NO oscillator has some important similarities and differences from the well-known Van der Pol oscillator [44] which has the form x¨−A (1−x2) x˙+x=g (t) (14) The form of the characteristic polynomial in Eq 13 implies that the shear-NO oscillator may be written as x¨+A (1−sgn (x˙) B/x2) x˙+ (x−x1) =g (t) (15) where the generic variable x fills the role of the radius in the shear-NO model, x1 is the nominal operating point and g (t) is a forcing function that can include steady, random or periodic components. The stability-determining second term in both oscillators can change sign based on the magnitude of the variable with a positive second term implying stability. The Van der Pol oscillator is known to self-sustain oscillations about the origin in phase space (x, x˙) = (0,0) when g (t) = 0, as its second term changes sign during different phases of each cycle. In contrast, the shear-NO oscillator operates near (x, x˙) = (x1,0), but with x > x1. Therefore, the second term in Eq 15 can be either i) always be positive (x12>B) regardless of the sign on x˙ implying inherent stability or ii) can be conditionally positive depending on both the magnitude of B and the sign of x˙. As a result, the shear-NO oscillator cannot produce self-sustained oscillations for large radius (x12>B). Furthermore, even when the radius is sufficiently small (x12<B), the radius will return unequivocally to equilibrium as long as x˙<0 unless a non-zero forcing function is present to change the sign of x˙. However, we find that when x12<B the magnitude of the forcing needed to start a new cycle can be arbitrarily small and in the form of either random noise or a periodic stimulus provided that enough time has passed for the system to approach its equilibrium point. In the context of the shear-NO dynamics, the key to oscillations is the inverse dependence on radius for the NO source due to shear stress in Eq 7. As long as the exponent on R−2 remains negative (increasing radius leading to lower shear stress and less NO production), then the NO-shear mechanism will be capable of a mathematical transition from unstable to stable as seen above in the generic oscillator in Eq 15. The physiological impact of this result then depends on the relative magnitude of the time scales identified herein, not on any single parameter value. For example, the stability of the NO-shear mechanism depends on groups of parameters such as tFNO = SNOFNO/E1KNOR2, which combines the sensitivity of the vessel to shear stress, the contractile force, the wall stiffness and the NO clearance rates. Inputs of constant magnitude have the effect of adjusting the equilibrium point. Using the shear NO oscillator as an example, an increase in transmural pressure will dilate the vessel, as will a pressure gradient that assists flow by inducing NO production via a steady shear stress. Likewise, a steady source of NO from local inflammation will chronically dilate the vessel [20]. If the vessel becomes sufficiently large, the stability criterion found above suggests that the shear-NO process will not support self-sustaining oscillations, in part due to the direct effect of radius on the stability criterion, but also due to greater stiffness of the wall at larger radius (higher E1). Nonlinearities in the force production and chemical source/elimination terms may also alter the stability in similar ways. Numerical simulation and examination in the phase plane reveal that the stretch-Ca2+ and shear-NO processes possess numerous symmetries that offer intriguing possibilities when the processes act together (Figs 2 and 6, Table 3). Most notably, we see that the shear-NO process produces rapid and unstable dilation toward a larger radius, followed by stable contraction, while the stretch-Ca2+ process causes the vessel to contract rapidly and unstably toward a smaller radius and then to dilate stably. An essential feature of both the Ca2+ and NO mechanisms is that taken separately they do not produce traditional, self-sustaining limit cycles, but instead have a one-sided stability near equilibrium from which a new cycle begins only with a perturbation from the local environment. Interestingly, a suitable trigger for the stretch-Ca2+ oscillator can be an increase in radius produced by the shear-NO mechanism. And conversely, the shear-NO oscillatory can be triggered by a contraction arising from the stretch-Ca2+ mechanism. Balanov et al [45] reviews a variety of similar, so-called “noise-induced” oscillators in contexts outside of lymphatic physiology such as neurons and electrical monovibrators, but to our knowledge, the coupling of symmetric, noise-induced oscillators described in the present study has not been previously investigated. Balanov et al [45] also review how nonlinear oscillators can synchronize with small-amplitude sinusoidal inputs. Here we found that synchronization of either oscillator can occur over a wide range of frequencies (shear-NO shown in Fig 8, similar behaviors for stretch-Ca2+ acting alone and in combination with shear-NO can be observed). The synchronization behavior seen here is similar to that for the forced Van der Pol oscillator in its ability to produce so-called Arnold tongues which are broad domains within which the input and output frequencies are locked in ratios of m: n where m and n are small integers [44,46]. Kornuta et al [47] recently showed that lymphatic vessels studied ex vivo synchronize their contractions in a 1: 1 fashion with imposed oscillatory variations in shear stress when the amplitude of the stimulus is sufficient large and the frequency of the input is relatively close to the autonomous frequency. Interestingly, they also observed that small amplitude variation in transmural pressure did not yield 1: 1 frequency locking. However, our examination of their results (Fig 8 in [47]) suggests that 2: 3 locking may have occurred. In the absence of imposed flow, they also found that the vessel continued to contract, but with a lower and more erratic frequency consistent with our simulated noise-triggered oscillator (Figs 2 and 3) in the absence of the shear-NO mechanism. Ohhashi et al [48] also examined sinusoidal variations in transmural pressure at frequencies well away from the spontaneous frequency. Here too, 1: 1 frequency locking did not arise, but the frequency of the contractions responded strongly to the input waveform. Given the subtlety of identifying non-1: 1 synchronization, further examination of the experimental record may be warranted. In conclusion, we have presented a model of a vascular oscillator. The present analysis is sufficiently general to point toward several features that are likely found in other systems. The linear stability analysis shows: (i) complementary mechanisms for dilation and contraction of collecting lymphatic vessels, (ii) a fast, unstable process that recovers slowly and stably to a one-sided equilibrium, (iii) disturbance-based triggering that facilitates either synchronization with a cyclic pacemaker or spontaneous oscillations from random disturbances and (iv) the capability for reciprocal modulation between contractile and relaxation effects. Those features are not only limited to the presented example of Ca2+ and EDRFs but can be extended into other fields. The ability of the Ca2+ and NO based oscillators to respond to each other and external stimuli explains how lymphatic pumping can be coordinated along extended lengths of collecting lymphatic vessels without the need for higher order coordination. This new class of coupled, noise-driven oscillator can help to explain the diverse pumping behavior of lymphatic vessels.
For decades, cardiovascular physiology has been an area of intense research, and we have a fundamental understanding of the mechanisms the heart uses to drive blood flow through the distributed network of vessels in the body. The lymphatic system is now receiving similar attention as more is learned about its functional role in disease processes. The importance of the lymphatic system in collecting excess fluid from tissues and returning it to the blood is well known, but how the lymph flow is regulated without a central pump is poorly understood. Each segment of collecting lymphatic vessel can independently contract yielding a network of distributed pump/conduits. This paper shows how the lymphatic muscle cells that squeeze fluid along the lymphatic vessels can be effectively regulated using only chemical and mechanical signals that they receive from their immediate microenvironment. Using stability theory and the tools of nonlinear dynamics we identify two complementary oscillators that respond to stretch of the vessel wall and shear of fluid flowing over the vessel wall. Numerical simulations of the combined oscillators show that they have characteristics well suited to the regulation of distributed systems in general and may have application in other biological and physical contexts.
acoustics
medicine and health sciences
neurochemistry
body fluids
classical mechanics
fluid mechanics
radii
mechanical stress
geometry
neuroscience
muscle contraction
mathematics
algebra
sound pressure
lymph
muscle physiology
neurochemicals
nitric oxide
fluid dynamics
continuum mechanics
fluid flow
physics
shear stresses
biochemistry
linear algebra
anatomy
physiology
biology and life sciences
physical sciences
eigenvalues
Synchronization and Random Triggering of Lymphatic Vessel Contractions