Synchronization and Random Triggering of Lymphatic Vessel Contractions

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.


Author Summary
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

Introduction
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 Ca 2+ 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 Ca 2+ 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 Ca 2+ 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 Ca 2+ levels are both subject to mechanical regulation; Ca 2+ 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 Ca 2+ 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 Ca 2+ 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 Ca 2+ 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.

Model Formulation
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 (Ca 2+ ) and dilatory (NO) signaling molecules. Neglecting inertial effects, the radial forces on the vessel wall balance as where the left hand side represents rate dependent effects with D incorporating the viscoelastic 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 Ca 2+ is transduced into a contractile force as F ¼ F Ca C Ca =ð1 þ aC NO Þ where α scales the possible desensitizing effect of NO [12]. The activation term A R can include a steady component from the mean transmural pressure difference p m 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][14][15]. Here we adopt the form EðRÞ ¼ Ae aR À p offset with stiffening coefficient a, scaling coefficient A, and offset pressure p offset 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 R 1 in the absence of dynamic increments to Ca or NO as EðRÞ % E 0 þ E 1 ðR À R 1 Þ where E 0 is the elastic force at equilibrium and a Taylor series expansion near equilibrium yields E 1 ¼ Aae aR 1 . We find the equilibrium radius by solving 0 ¼ D dR dt ¼ À ðAe aR 1 À p offset Þ À F Ca S Ca 0 =K Ca þ p m for R 1 where p m is the mean transmural pressure. Given the stiffening behavior of the wall (E 1 / e aR 1 ), we expect that appropriate values of E 1 will be larger at higher mean transmural pressures where the equilibrium radius will be somewhat larger.
The concentrations of the signaling molecules (i 2 {Ca, NO}) are governed by the generic conservation law where all concentrations are taken to be dimensionless ratios relative to a suitable reference, K i is clearance of the signaling species through chemical reaction, transmembrane ion pumps and advective-diffusive transport, S i is a dynamic source term for the signaling molecule and A i 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.

Stretch-Ca Dynamics
Since our focus is on the interactions between Ca 2+ and NO, we introduce a minimal representation of the Ca 2+ dynamics rather than a fully detailed model of Ca 2+ oscillations as may be found in the literature [10,16]. We retain the following features: i) at rest, Ca 2+ is at a low concentration in the cytoplasm of LMCs; ii) a contraction is initiated when Ca 2+ 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 Ca 2+ 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 Ca 2+ levels have returned to near equilibrium. As the Ca 2+ 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][20][21]. We mathematically express the release of Ca 2+ into the cytoplasm from intracellular stores and the extracellular fluid as the sum of a steady source S ca0 needed to maintain the baseline Ca 2+ concentration and a transient component that is sufficiently rapid to be modeled as an impulse function δ(t) where t is the time since the Ca 2+ concentration most recently passed the threshold necessary to trigger another contraction C CaThresh . This is expressed as, where S Ca1 /(1 + γC NO ) is the magnitude of a bolus of Ca 2+ with a possible reduction due to NO that is scaled by γ. We model the clearance of Ca 2+ from the cytoplasm with where K Ca is a rate constant and β scales the possible enhancement of Ca 2+ clearance attributed to NO [12,16]. At high concentrations of Ca 2+ 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 and dC Ca dt 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 S Ca1 and K Ca now include the adjustments due to NO introduced in Eqs 3 and 4.

Shear-NO Dynamics
In the previous section we developed the model so that it reproduces Ca 2+ 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 ¼ À 2pR dR dt which when integrated with respect to z along a vessel yields QðzÞ ¼ À 2pR dR dt z þ C 1 where the constant of integration depends on the end conditions for the lymphangion. When the segment is contracting dR dt < 0 and the upstream valve is closed [Q(0) = 0] we have QðzÞ ¼ À 2pR dR dt z. Alternatively, if the segment is expanding dR dt > 0 and the downstream valve is closed (Q(L) = 0), we obtain QðzÞ ¼ 2pR dR dt ðL À zÞ. We can express the mean flow along the length more compactly as " Q ¼ pRj dR dt jL. When additional flow Q 0 is imposed on the segment by an axial pressure gradient we have " 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 t ¼ 4m " Q pR 3 . The mean shear stress along the segment due to dynamic changes during contraction or expansion is therefore approximately t ¼ 4mL 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 where S NO has absorbed the remaining constants in the shear stress expression and S NO 0 represents NO released due to the through-flow term Q 0 or chronic sources of NO such as might arise during inflammation. The source term S NO 0 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. S NO 0 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 ¼ À 8mQ pR 4 . When contracting we have @p @z ¼ 16m dR dt z, which when integrated along the length gives pðzÞ ¼ 16m Averaging over the length of the segment yields " dR dt þ 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 dR dt , 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 F Ca and the elastic effects. Mathematically, for small αC NO , we can write this as By defining F NO F Ca C Ca0 α, we can write And the net force from Eq 8 becomes where the net contractile force is decomposed into a positive term set by the baseline Ca 2+ levels and a negative term that represents how the Ca 2+ levels are modulated by NO. Thus, the NO-dependent term is not a true outward force, but arises mathematically from a reduction in the Ca 2+ -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 Ca 2+ 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 S NO and F Ca C Ca0 α, 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][27][28][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 Ca 2+ [32]. Parameters such as the contractility, NO production and the mechanical stiffness appear to depend on anatomical location, species and age [13,28,[33][34][35][36], suggesting that the full range of possibilities realizable in vivo awaits further investigation.

Results
Here we first present analytical and numerical results for cases in which the stretch-Ca 2+ dynamics are active, but with constant NO levels. We will then present results of our model for when shear-NO dynamics are active, but Ca 2+ levels do not spike, but remain near baseline. Finally, we will present results of our model for fully coupled stretch-Ca 2+ and shear-NO processes.

Stretch-Ca 2+ Dynamics with Constant NO Levels
In the absence of dynamic activation, Eq 6 implies the Ca 2+ concentration during each contraction cycle will decay as C Ca ðtÞ ¼ S Ca1 e À K Ca t . Using this as an input to Eq 5, we find that the radius varies from its baseline value during each contraction cycle as where we see that the return to equilibrium depends on two characteristic times, one set by the rate of Ca 2+ clearance t Ca = 1/K Ca and the other by the mechanical lag t mech = D/E 1 which can include the lag between the concentration increase and force production. As a reference time scale we have selected t Ca = 1/K Ca = 1 s which is in the range observed by Shirasawa and Benoit [15]. They observed a similar lag between the rise in Ca 2+ concentration and the peak force generation. Here we use this lag as an estimate of t mech which we take to incorporate the viscoelastic 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 t c that determines the return to equilibrium. Experimental observations of the magnitude of radius change as a function of pressure show Offset pressure for wall stiffness p offset 100 Pa [15] Wall damping coefficient D 3x10 6 Pa-s/m [15] Baseline Ca concentration C Ca0 1 Baseline Ca release rate S Ca0 1 Ca release pulse size S Ca1 0.85 [15] Ca force production coefficient F Ca 100 Pa [15] NO source from shear stress coefficient that it decreases with increased internal pressure [21]. This phenomenon is reproduced by our model as the vessel wall stiffens (larger E 1 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 t Ca and the magnitude of the random activation term. In addition to the steady source of Ca 2+ that establishes the vascular tone, we include a random component A CaRand ð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 Ca 2+ 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 σ / p m ). This will lead to a higher frequency at larger radii as found experimentally [21]. Exponential decay of Ca 2+ near equilibrium leads to a latency period between contractions that varies as T = t c log (C max /σ) where C max is the magnitude of the Ca 2+ 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 E 1 ) leading to a smaller mechanical time constant (smaller t mech ) and reduced amplitude for change in the radius.
We note that our simple model of stretch-Ca 2+ 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 t c log(C max /σ). We also see that increased transmural pressure can reduce the interval between contractions by stiffening the wall (p m ") R 1 ") E 1 ") t mech #) t c #) T#) and also by increasing the sensitivity of the Ca 2+ channels by stretching the vessel wall (p m ") σ ") T #) yielding higher frequency contractions (compare Fig 2b and 2e).
We also find that our model of the stretch-Ca 2+ process readily synchronizes when we impose extrinsic rhythmic pace-making since only small variations in the Ca 2+ concentration relative to the threshold level are needed to initiate the next contraction cycle (Fig 3). Such small variations in Ca 2+ 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].

Shear-NO Dynamics with Ca 2+ Near Baseline Levels
The conditions for oscillations in radius to arise near baseline Ca 2+ levels in the absence of sharp spikes in Ca 2+ as considered in the previous section are available from linear stability analysis of the shear-NO process near a point R 1 which yields  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: 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 Ca 2+ concentration described in Eq 6. If the radius is large enough so that E 1 /D + K NO > S NO F NO /DR 2 , then the fixed point is inherently stable. If instead, the radius is small enough that E 1 /D + K NO < S NO F NO /DR 2 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 S NO F NO /DR 2 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 t c % t mech + t NO + t FNO where the three contributions arise from mechanical lag t mech as before, the clearance of NO t NO = 1/K NO , and the rate of force modulation by NO t FNO = S NO F NO /E 1 K NO R 2 . To a similar degree of approximation, the instability of the NO-shear dynamics requires t mech + t NO t FNO . 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 The arrows indicate the direction of increasing R 1 . The eigenvalues during contraction always have a negative real part indicating stability. When the baseline radius is large enough to move the eigenvalues beyond point B the system is inherently stable during dilation. At point B the system is marginally stable and will oscillate at frequency f = (E 1 K NO /D) 1/2 /2π. For smaller baseline radius, between points A and B, the response is unstable, but oscillatory. And when the baseline radius is smaller at point A, the dynamic component of the radius increases exponentially without oscillation until the radius is large enough to reach the range between A and B where oscillations can occur. stretch-Ca process. Here the period of the NO-induced oscillations varies as t c log(C max /σ) as before but t c and C max now refer to NO rather than Ca 2+ . As with Ca 2+ , 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, Thereafter, the trajectory begins a stable return (blue arrows _ R _ < 0) to equilibrium at point 0. NO reaches its peak concentration at point 1 before the radius reaches its maximum, but this point does not directly influence the stability of the system. 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.

Combined Stretch-Ca 2+ and Shear-NO Dynamics
Having explored the system dynamics when Ca 2+ 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 Ca 2+ as modeled by Eq 10, ii) modify the availability of Ca 2+ by 1/(1 + γC NO ) or iii) speed clearance of Ca 2+ by 1 + βC NO [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 Ca 2+ -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 = (E 1 K NO /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-Ca 2+ 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-Ca 2+ and shear-NO mechanisms, but will typically be dominated by the faster of the two processes. Long latency intervals between Ca 2+ -induced contractions can permit NO to produce an unstable dilation, whereas, short intervals due to Ca 2+ can suppress the autonomous oscillations possible through the NO mechanism. Interestingly, the published clearance rates for Ca 2+ 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 Ca 2+ 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.

Discussion
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 The form of the characteristic polynomial in Eq 13 implies that the shear-NO oscillator may be written as where the generic variable x fills the role of the radius in the shear-NO model, x 1 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Þ ¼ ðx 1 ; 0Þ, but with x > x 1 . Therefore, the second term in Eq 15 can be either i) always be positive (x 2 1 > 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 (x 2 1 > B). Furthermore, even when the radius is sufficiently small (x 2 1 < 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 x 2 1 < 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 t FNO = S NO F NO /E 1 K NO R 2 , 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

Effect of Noise on Period
Increasing period with smaller noise T % t c log(C Camax / σ) Increasing period with smaller noise T % t c log(C NOmax / σ) (higher E 1 ). 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-Ca 2+ 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-Ca 2+ process causes the vessel to contract rapidly and unstably toward a smaller radius and then to dilate stably. An essential feature of both the Ca 2+ 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-Ca 2+ 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-Ca 2+ 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.

Effect of Radius on Frequency
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-Ca 2+ 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 socalled 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 onesided 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 Ca 2+ and EDRFs but can be extended into other fields. The ability of the Ca 2+ 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, noisedriven oscillator can help to explain the diverse pumping behavior of lymphatic vessels.