A dual-clock-driven model of lymphatic muscle cell pacemaking to emulate knock-out of Ano1 or IP₃R

The lymphatic system vasculature is a network of vessels, nodes, and valves that returns interstitial fluid from between cells back to the blood circulation, operating in parallel with the venous system (Moore and Bertram, 2018). It also plays a vital role in the immune system as the returning fluid is filtered via lymph nodes for a range of particles, bacteria, and viruses. Lymphatic system defects are involved in numerous diseases, such as obesity, cardiovascular disease, inflammation, and neurological disorders such as Alzheimer’s disease (Mortimer and Rockson, 2014; Oliver et al., 2020). Lymphoedema is a long-term dysfunction of the lymphatic system wherein the fluid is not properly returned to the circulatory system and instead accumulates in the interstitium (Rockson, 2001); there is currently no treatment for the underlying causes.

The return of interstitial fluid via the lymphatic vascular system is driven by both intrinsic and extrinsic pumping, where intrinsic pumping results from periodic contractions of the lymphatic muscle in the vessel walls and accounts for about two-thirds of lymph flow in the extremities at rest (Olszewski and Engeset, 1980), with the extrinsic remainder a result of passive squeezing of lymph vessels (Breslin et al, 2019). Only the intrinsic contractile activity can be regulated, and so it is essential for the matching of lymph transport to demand as represented by the filling state of the initial lymphatics (Scallan et al., 2012), itself a result of the extent of blood capillary extravasation of fluid. The extent of the required contractile activity is also a function of posture, which may set up adverse hydrostatic pressure gradients (Zawieja, 2009; Li et al., 2022). Higher pressure in the contractile collecting lymphatics leads to more frequent contractions; this is pressure-dependent chronotropy (McHale and Roddie, 1976; Benoit et al., 1989; Zawieja et al., 1993). The frequency of contractions is determined by what is at present believed to be the autonomous pacemaking ability of lymphatic muscle cells (LMCs) themselves. Consequently, the study of the LMC pacemaking ability investigates a function that is vital to fluid homeostasis in general.

An incomplete understanding of the mechanisms involved in the oscillatory excitation for intrinsic pumping restricts the potential development of pharmacological or other treatments (Lee et al., 2022) to attempt to deal with system dysfunction, whether this be genetic or acquired. The dysfunction may be mechanical or innate to the pacemaker itself. A prominent example of the former is the loss of contractile ability in lymphatic beds that have been irradiated during cancer treatment. The latter is exemplified by Cantú syndrome, which Davis et al. (2023) have shown in a mouse model stems from a gain-of-function mutation in one particular ion channel in the LMC pacemaker. The resulting loss of pumping capability leads to lymphoedema in >50% of Cantú syndrome patients. In another example (Kim et al., 2023), long-QT type II patients may have mutations causing ERG channel loss-of-function or impaired trafficking of the channel to the cell membrane, again leading to lymphatic contractile dysfunction. Due to the complexity of the pacemaking system, modelling is critical to uncovering the function and failure of the underlying oscillatory mechanisms. While LMCs are almost uniquely complex¹ in mounting both tone and short-lived contractions (Muthuchamy et al, 2003), for the inbuilt pacemaker behind the latter they can draw extensively from models of oscillations in other cell types. Such models in a range of cell types have historically been based on either a so-called “membrane clock” (M-clock) as originated by Hodgkin and Huxley (1952) or a “calcium clock” (C-clock); see, e.g., Goldbeter et al. (1990). The M-clock is composed of periodic flows of ion currents in the cell membrane, resulting in an oscillating membrane voltage output consisting of rhythmic action potentials (APs) separated by a period of relative quiescence. Propagation of the AP along the fibre as a signal is the purpose of the nerve, as studied by Hodgkin and Huxley, and in the conduction system of the heart. The C-clock represents the oscillating levels of calcium in the cell, and the models are composed of the calcium concentrations in the main part of the cell (cytosol) and in the calcium store (endoplasmic reticulum, ER), and mathematical descriptions of the calcium flux across the ER and cell membranes. In muscle, where the store is termed the sarcoplasmic reticulum, the periodic peaking² of [Ca²⁺], giving rise to mechanical contraction, is the purpose. Models of pacemaking in cardiac cells were based around the M-clock (McAllister et al., 1975), while a range of models of other cell types have focussed on oscillations in C-clocks (Dupont et al., 2016). More recently there has been a shift, backed by experimental data, in the models of cardiac pacemaking to dual-coupled clocks, where neither clock is dominant in triggering APs (Yaniv et al., 2015).

Lymphatic muscle has so far seen limited modelling of pacemaker mechanisms. Hald et al. (2018) created a two-variable model of a voltage oscillator (i.e., M-clock only) based on modifications to a barnacle muscle model (Morris and Lecar, 1981). This model was created in order to study communications between neighbouring LMCs (and between neighbouring endothelial cells), and had limited detail for internal mechanisms. Imtiaz et al. (2007) added M-clock components to an earlier calcium-based model by Dupont and Goldbeter (1993). In this model, the C-clock drove the oscillations and triggered the AP in the M-clock. Like Hald et al. (2018), Imtiaz et al. (2007) created a model of multicellular AP conduction.

In a previous work (Hancock et al., 2022), we extended the single-cell model of Imtiaz et al. (2007) by modifying mechanisms in the M-clock: specifically, by altering the model components for the L-type calcium channel. These modifications enabled the M-clock to oscillate independently and drive oscillations in the C-clock, which was necessary for the model to match experimental data showing survival of APs after knock-out of anoctamin-1 (Ano1) calcium-activated chloride–ion channels. These channels form a key M-clock mechanism of coupling from the C-clock.

In this paper, we direct attention to the inositol 1,4,5-triphosphate (IP₃) receptor, which is known (Hille, 2001; Foskett and Mak, 2010; Emrich et al., 2021) to be the main source of calcium ion release from the store to the cytoplasm in smooth muscle types, including lymphatic muscle. As such it is essential for the LMC’s C-clock, and its numerical knock-out automatically stops pacemaking in any model such as that of Imtiaz et al. (2007), which is C-clock driven. However, this is not the case experimentally. Fig. 1 shows a comparison of APs recorded from mouse LMCs under matched control conditions and after deletion of the dominant IP₃ receptor 1 (IP₃R1-KO)³. Clearly, APs survive, albeit with reduced frequency. The systolic plateau phase is lengthened after IP₃R-KO and the total minimum-to-maximum excursion of the membrane potential is less.

The extent of the frequency reduction after IP₃R-KO is variable; Fig. 2 shows another example. Whereas in Fig. 1 the frequency reduces to 46% of control, in Fig. 2 it reduces to <37% of control⁴. Again, the systolic plateau phase is pronouncedly lengthened after IP₃R-KO, and the total excursion is less. In Fig. 2, the total excursion is greater under control conditions than in Fig. 1 because of both a lower minimum immediately after repolarisation and a higher maximum early-systolic spike, whereas the IP₃R-KO waveform in Fig. 2 displays a similar excursion to that in Fig. 1, but the whole waveform is raised relative to its Fig. 1 counterpart so that the membrane potential is considerably higher than control during the diastolic slow depolarisation phase.

In this paper, by modifying the C-clock mechanisms, we extend our existing model to the case where the M-clock and the C-clock are both driving oscillations. We first analyse our existing model to show how that model fails to emulate the observations of knock-out of the IP₃ receptor. We then analyse a special case of the proposed model to illustrate two different means by which the C-clock can control oscillation frequency in the M-clock: first, by varying the frequency of C-clock oscillations, and second, by varying the [Ca²⁺] level of a non-oscillating C-clock. Finally, we analyse the proposed model for the general case and show that the model largely matches features of the IP₃R-KO experimental data. We use time traces and phase-plane analysis to show that the C-clock and the M-clock are both driving their synchronised oscillations.

The rest of the paper is set out as follows. The Background section provides background information on the existing model that we extend in this paper. The Materials and methods section describes the new model and methodology used in this paper. The Results section presents the results from the proposed model, including the case where both clocks are driving oscillations and phase-plane analysis explanations. The Discussion section provides a discussion of the results, and the Conclusion draws brief conclusions.

In this section, we provide a background on our existing model (see Hancock et al. [2022] for details). That model matched the experimentally observed membrane potential oscillations under both control and Ano1-KO conditions, with Ano1-KO leading to larger APs and a much longer period, as seen in Fig. 3. However, Fig. 3 also shows the membrane potential computed with the IP₃ receptor disabled (IP₃R-KO) by setting the maximum IP₃R flux J_x3 to zero (see table of notation and equations of Hancock et al. [2022]). The added traces show that, under IP₃R-KO, once start-up transients decay, the existing model produces APs which are similar in all respects to those produced under control conditions. This is not consistent with the observed IP₃R-KO behaviour as shown in Figs. 1 and 2. The data for this latter numerical case were not part of the original study and are a focus of the work in the current paper.

The existing model is composed of an M-clock and a C-clock using four time variables. The C-clock is represented by two calcium variables: the cytosolic calcium concentration Z and the ER-stored calcium concentration Y. The M-clock is represented by two membrane variables: the membrane potential V and an inactivation gating variable h for the L-type calcium channel; see Eqs. 1, 2, 3, and 4, and Fig. 2, of Hancock et al. (2022).

Hancock et al. (2022) used phase-plane analysis to analyse the oscillation of the clocks. Phase-plane analysis shows the relationship between two different time variables on a single graph of their combined trajectory rather than plotting them individually against time; see Strogatz (2015) for further details. To simplify the analysis, the phase-plane analysis for the two clocks was separated, as shown in Fig. 4. The trajectories of the two clocks were plotted on their respective phase planes. The nullclines of the M-clock were also plotted, whilst assuming that the (C-clock-dependent) calcium concentration Z was constant—see Fig. 4. Similarly, the nullclines of the C-clock were plotted assuming that (M-clock-dependent) L-type calcium current I_L was held constant. A nullcline of a variable shows all points where that variable does not change with respect to time. It is useful for interpreting the phase plane since it divides the plane into regions where the time derivative of the variable has different signs.

In the existing model, the M-clock drove oscillations of both [Ca²⁺] and membrane potential. This was required to match the Ano1-KO data. In the model of Imtiaz et al. (2007), the M-clock was driven by oscillations arising in the C-clock. Under these circumstances, oscillations of membrane potential (APs) were impossible with Ano1-KO because the latter condition removed the mechanism for the C-clock to drive the M-clock. In Fig. 4 c, elevated Z depresses the peak of the V-nullcline. This can happen dynamically, as is exemplified in Fig. 4 by moving from panel d to panel c. Consequently, the M-clock orbital trajectory in Fig. 4 c is relatively small, indicating increased M-clock frequency and decreased peak membrane potential V. Conversely, Ano1-KO, which for the M-clock amounts to Z = 0, leads to the maximum-sized orbit, indicating decreased AP frequency accompanied by increased peak V and lowered minimum V.

In the existing model, the C-clock had the means to oscillate independently when uncoupled from the M-clock, but this did not occur when the M-clock was connected. In Fig. 4, a and b, the Z-nullcline is shown for two fixed values of I_L, the first (I_L = 0) pertaining to when the C-clock was disconnected from the M-clock, and the second (I_L = −10 pA) exemplifying values occurring when the clocks were coupled. Fig. 4 a shows trajectories from an arbitrary starting location for the C-clock independent of the M-clock. It can be observed that the C-clock oscillated between relatively low values of Z when I_L = 0 but approached a stable equilibrium at a considerably higher value of Z when I_L = −10 pA. Fig. 4 b shows the orbital trajectories when the C-clock was coupled to the M-clock. In Fig. 4 b, Z oscillates between extremes of 2.9 and 4.0 μM under control conditions (solid black line), and between 1.5 and 4.0 μM with Ano1-KO (dashed black line). Under both conditions, the orbit surrounds a nullcline intersection point that is stable (Fig. 4 a), showing that the C-clock is not oscillating in its own right. For this case, the oscillation is driven by changes in I_L and thus by the M-clock.

In this paper, we focus on modifying and analysing the C-clock model, as modelling IP₃R mechanisms is central to fitting the IP₃R-KO data. To achieve this, we pose and answer two questions. First, can the C-clock model be improved by using more modern mechanistic approaches (see Dupont et al. [2016] for background on calcium modelling)? Second, can the C-clock and M-clock each drive oscillations in the other, rather than just the M-clock driving the C-clock?

A moderately detailed description of the biological methods for the experiments in general and for securing genetic ablation or pharmacological inhibition of calcium-activated chloride channel Anoctamin 1 channels, i.e., Ano1 knock-out, was given by Hancock et al. (2022). All membrane potential recordings come from murine inguinal-axillary lymphatic vessels, cannulated and pressurised ex vivo. Here, genetic means were again used to obtain IP₃R knock-out using Myh11CreER^T2;IP3R1^fl/fl (IP3R1ismKO) mice and equivalent non-activated Cre controls. The companion paper to this one (Zawieja et al, 2023) gives full details.

In this section, we introduce a model for LMC oscillations by improving on our previous model of LMC pacemaking (Hancock et al., 2022). We modify the model of the C-clock to alter the coupling between the M- and C-clocks and the effect of IP₃. This modification incorporates a previous IP₃R model by De Young and Keizer (1992), as simplified by Li and Rinzel (1994).

In the newly adopted IP₃R model, calcium release from the store is a function of store [Ca²⁺] and cytosolic [IP₃] and [Ca²⁺], the latter property representing calcium-induced calcium release (CICR). While the former IP₃R model included less exact representations of all these, the new model, via the time-dependent properties imbued by the gating variable n, also emulates the sequential activation and inactivation of the IP₃R (Dupont et al., 2016).

The above equations were coded in Matlab (RRID:SCR_001622; The MathWorks, Inc.) and solved by finite differences on a single processor of a standard PC. Nullclines were found, in most cases iteratively, by solving the appropriate equation from Eqs. 1, 2, 3, 4, and 5 with the time derivative set to zero. All results shown are from times much later than required for start-up transients to decay, apart from those for the existing model in the left panel of Fig. 3. Matlab code for solving the equations is available from an open-access repository; see Data availability at the end of the paper.

The above analysis applies to the two-dimensional M-clock of this paper, i.e., it can be used to determine whether an intersection between the nullclines dV/dt = 0 and dh/dt = 0 is stable or not. However, the C-clock is three-dimensional, so in general, the analysis cannot be applied. We portray the C-clock herein by way of two-dimensional phase-plane slices of the (Z,Y,n)-space. Provided a single such slice, e.g., a phase plane in (Z,n)-space with Y = const., is considered, the analysis can be applied to the nullcline intersection(s) found therein.

Videos 1 and 2 animate the control cycle and the IP3R-KO cycle, respectively, for the conditions which will be shown in Fig. 12 below. Fig. S1 and S2 are related to Videos 1 and 2, respectively.

In this section, we present time traces and phase-plane analysis for the proposed models of oscillations in membrane potential and calcium concentration. We show different methods whereby the C-clock can control the oscillations of the M-clock. We also show that the proposed model qualitatively matches the experimental control, Ano1-KO, and IP₃R-KO data, and analyse the coupling of the M-clock with the C-clock.

We first analyse the two different ways that the C-clock can bring about AP initiation in the M-clock. To analyse this triggering, we mathematically removed the effect of the M-clock on the C-clock while leaving the C-clock effect on the M-clock. We set the dimensionless binary constant c = 0 in Eq. 1 above, which removes the effect of the L-type calcium channel on the cellular calcium concentration [Ca²⁺] while keeping its effect on voltage. This modification is equivalent to the L-type calcium channel passing not calcium but instead another bivalent cation, which changes the membrane potential V but not the cytosolic [Ca²⁺] (algebraically, Z). While the L-type Ca²⁺ channel acts on both clocks, it is an integral part of self-driven M-clock oscillation. In the C-clock, it is only a mechanism whereby the M-clock acts on the C-clock, and it is this coupling that we now remove. The modification thus allows us to study the temporal dynamics by which the M-clock is coupled to the C-clock.

We investigate two cases where the C-clock can modify the oscillations of the M-clock. In the first (c = 0), which will be the subject of Figs. 6 and 7, the isolated C-clock can be at a quiescent steady state. In the second, which will be the subject of Figs. 8, 9, and 10, the isolated C-clock can spontaneously oscillate and thereby manipulate the M-clock. The first case is idealised as it could not occur for the real case of c = 1 if there are APs. However, this case is informative as an approximation and a conceptual illustration.

Fig. 6 shows superimposed time traces from the simulated system with c = 0 in the three situations of (1) control conditions (all parameters take their default values), (2) Ano1 channel knock-out (zero Ano1 conductance), and (3) IP₃ receptor knock-out (zero Ca²⁺ efflux from the store via the IP₃ receptor). The case of IP₃R-KO displays constant [Ca²⁺], i.e., Z, after start-up transients decay, whereas the control case exhibits time-varying Z. Ano1-KO also displays time-varying Z, but will be regarded as a reference case in the sense that it is equivalent to Z = 0 as input to the Ano1 channel, thus disabling the mechanism by which the C-clock triggers APs in the M-clock.

Both IP₃R-KO and control conditions modify the calcium oscillation from the Ano1-KO reference case (upper panel of Fig. 6). For IP₃R-KO, the constant but non-zero [Ca²⁺] results in an increase in the frequency and decrease in amplitude of the membrane potential oscillation. In the control condition, the oscillating [Ca²⁺] also results in an increase in the frequency and decrease in amplitude of the voltage oscillation. In this case, the frequency of the voltage oscillation is identical to that of the calcium oscillation (cf., the two panels of Fig. 6). In contrast, the M-clock and C-clock oscillatory frequencies are not the same in the Ano1-KO case. In the control condition, oscillatory Z applied to the Ano1 channel has caused the M-clock to synchronise with the oscillation in the C-clock. Thus, the AP frequency is determined by the calcium flash frequency in the control condition, while it is determined simply by the steady [Ca²⁺] level in the IP₃R-KO case.

Phase-plane analysis was used to investigate why this system behaviour occurs. Fig. 7 shows the orbital trajectories following the left and right arms of the V-nullcline, as discussed in the Background. The AP is triggered when these trajectories reach the peak of the V-nullcline’s left arm and jump to the right arm (the default value of K_h = 4 µM for the half-saturation [Ca²⁺] of the h_∞Z-dependence ensures that the h-nullcline is stationary). The position of this V-nullcline peak (and whether there is one) varies depending on the instantaneous [Ca²⁺], Z(t). But there are two cases in which this occurs. First, in both Ano1-KO and IP₃R-KO, the AP occurs when the trajectory surmounts a fixed V-nullcline peak due to fixed Z (zero for Ano1-KO, non-zero for IP₃R-KO). Constant values of Z greater than zero result in a lower peak and consequently lead to higher AP frequencies and smaller AP amplitudes. But in the control condition, the AP occurs when the trajectory surmounts a V-nullcline peak that is changing with time due to changing Z. In this case, the cyclic Z(t) increases until the V-nullcline peak is depressed enough for an AP to occur. So, the AP always occurs on the upstroke of the Z waveform and is thus synchronised with the C-clock oscillation.

We next return to the full model where the M-clock and the C-clock are bi-directionally coupled together. To achieve this, we set c = 1 in Eq. 1 above so that the L-type calcium channel directly alters both the calcium concentration and the membrane potential.

Fig. 8 shows superimposed time traces from model runs under IP₃R-KO, Ano1-KO, and control conditions. All three oscillations have different forms. IP₃R-KO leads to oscillation of membrane potential with larger amplitude and period than control. Ano1-KO oscillation has yet larger amplitude and period. The relative order of these oscillation periods qualitatively matches the experimental data shown in Figs. 1 and 2, resolving the issue with the existing model (Hancock et al., 2022), which was discussed above.

In Fig. 8, we can study the AP trigger by comparing the oscillation in the M-clock with that in the C-clock. Under Ano1-KO, the [Ca²⁺] oscillation, Z(t), has a different frequency⁷ to that of the membrane potential, V(t). The C-clock is not triggering oscillation in the M-clock because the Ano1 coupling of the C-clock to the M-clock is missing. The peak of the [Ca²⁺] flash (peak Z) varies from one flash to the next. It is particularly high when there is an AP, this being caused by the inrush of calcium ions via the L-type calcium channel. In contrast, under control conditions, the oscillations of V(t) and Z(t) are synchronised; the calcium oscillation affects the triggering of the AP by the outrush of chloride anions via the Ano1 channel.

Under IP₃R-KO, the oscillations of V(t) and Z(t) are again synchronised, as in the control case. However, as can be observed in the M-clock phase plane in Fig. 9 a, the AP is triggered near the V-nullcline peak, which occurs when Z takes its minimum value of Z = 0.368 µM. To that extent, the orbital trajectory for IP₃R-KO resembles that for Ano1-KO and differs from that for control conditions. In the control situation, the AP is launched when h(t) is around 0.21, whereas the left-arm peak of the V-nullcline reaches h = 0.77 at another time during the cycle. Thus, in control conditions, rising Z(t) causes the V-nullcline peak to descend to meet the slowly rising operating point during the diastolic depolarisation phase of the cycle. The peak descent is in this case greater than the operating-point rise, whereas under IP₃R-KO, for this particular set of parameter values, the peak descent before the AP starts is negligible.

Fig. 9 a also shows that the IP₃R-KO cycle involves a slightly increased “resting” membrane potential, relative to both the control and Ano1-KO conditions. This is seen both in the V-nullcline, which at minimum Z intersects the axis h = 0 at a slightly higher V, and in the orbital trajectory, which never falls below −64 mV (see also Fig. 8). The slightly increased resting or diastolic membrane potential in IP₃R-KO is a consequence of the fact that Ano1 current continues throughout diastole; under control conditions, it decays to zero during diastole and in Ano1-KO it is of course always zero. The continuing Ano1 current for IP₃R-KO is a result of cytosolic [Ca²⁺], i.e., Z, remaining substantially higher through diastole, at around 0.4 µM, than is the case under control conditions (again, see Fig. 8). That is, in turn, a consequence of high leakage of Ca²⁺ ions from the store as a result of the extreme elevation of store [Ca²⁺], i.e., Y, when there is no efflux via IP₃ receptors and only leakage available to limit the concentration Y.

Fig. 9 b shows the (Z,n)-space phase plane⁸, including the static n-nullcline and the two extreme positions of the dynamic Z-nullcline. The position of the nullcline intersection varies with the Z-nullcline over this range. To test the stability of the intersection point, a reduced model was built in which Y was reduced to a linear function of Z, fitting the close approximation to such linear dependence shown by the trajectory in Fig. 9 c. This reduces the C-clock to two dimensions, thus enabling the use of the test outlined in Nullcline intersection stability section. For the extreme positions of the Z-nullcline in (Z,n)-space in the reduced model, the nullcline intersection was oscillatory but stable for the smaller magnitude I_L and oscillatory and unstable for the larger-magnitude I_L, showing that the reduced C-clock was capable of oscillating autonomously. Furthermore, when the full three-dimensional C-clock was isolated by setting a fixed value of I_L corresponding to the mean of the values traversed cyclically when I_L was set by the M-clock (−2.49 pA), trajectories diverged from the Z- and n-nullcline intersection point. The C-clock is thus operating in a region where it can drive its own oscillation. This contrasts with the trajectories from the existing model in Fig. 4 b, which circle around a nullcline intersection far to the right which, by the Nullcline intersection stability section test, can be shown to be a stable equilibrium point. The C-clock in that model was therefore not self-sustaining oscillation, and so the observed C-clock oscillation was being driven by the M-clock. This observation matches with the simulations in Fig. 6 showing C-clock oscillation in Z for Ano1-KO that is independent of, and thus not driven by, the M-clock oscillation in V. So, in this new model the C-clock and M-clock oscillators are coupled to and driving each other, unlike the existing-model case shown in Fig. 4, in which the M-clock was the sole driver.

Fig. 9 c shows that the orbital trajectory in the (Z,Y)-space phase plane traces almost a straight line, with Y decreasing as Z increases, and vice versa. The slope dY/dZ of this line is −5.0, i.e., the magnitude is almost the same as the ratio v_c/v_s = 5.5 (Table 1). This shows that, at least for the parameters used in this simulation, the total amount of calcium in the cell is relatively constant and is an approximately conserved quantity. The dynamics of intracellular [Ca²⁺] are thus dominated by the calcium passing between the cytosol and the ER. Mathematically, this relationship also suggests that for the parameters used in this simulation, with Y essentially determined by Z, the (Z,n)-space phase plane is the appropriate one to analyse to understand the C-clock behaviour. The Z,Y plane also shows that the Y- and Z-nullclines are very similar, as is to be expected given the almost linear relationship between Y and Z. The trajectories in this plane are driven by changes in the nullclines due to changes in the IP₃R gating variable n.

Largely keeping the same parameter values as used for Figs. 8 and 9, we can also enquire of the model what happens when Ano1 channels are blocked (by setting G_xC = 0) and simultaneously IP₃ receptors are blocked (by setting k_f = 0). The outcome is illustrated in Fig. 10, which compares the combined knock-out situation with all the other conditions, for waveforms of Z, V, I_L, and I_Ano1. The simulated combined knock-out produces APs that are identical to those shown for Ano1-KO alone in Fig. 8. However, the intracellular [Ca²⁺] varies cyclically in a way that is different from either the Ano1-KO trace or the IP₃R-KO trace in Fig. 8: the level, amplitude, and waveform shape are essentially those of the IP₃R-KO trace, but the frequency of calcium flashes is reduced to that of Ano1-KO so that each AP is accompanied by a calcium flash. Combined KO of Ano1 and IP₃R leaves I_L dynamics intact, but the store release of Ca²⁺ is now solely by leakage. The comparison of these four situations shows an intriguing difference in relative timing: under all three knock-out conditions (Ano1, IP₃R, and both), the cytosolic [Ca²⁺] peaks just before the descending phase of the AP steepens, whereas in control conditions peak [Ca²⁺] occurs earlier, near the middle of the sloping plateau phase of the AP.

One of the goals of creating models of LMC pacemaking is to understand the mechanisms whereby the frequency of spontaneous contractions increases with vascular distending pressure (e.g., see Fig. 2 of Bertram et al. [2017] for data from an isolated segment of rat mesenteric lymphatic vessel). In lymphatic endothelial cells, a mechanosensitive ion channel, Piezo1, mediates the sensing of oscillatory fluid shear stress (Choi et al., 2019), but Piezo1 deletion from vascular SMCs does not affect pressure sensing (see Figs. S3 and S4 of Retailleau et al. [2015]; we confirm the same in LMCs in unpublished work). Thus, the ion channels responsible for this pressure/frequency dependence have yet to be identified. However, variations in IP₃ concentration in the LMC lead to variations in the frequency of pacemaking in models, including this one. The [IP₃] variations are primarily driven by the quantity of external agonist (a neurotransmitter or hormone) that binds to receptors in the cell membrane. In a cascade of subsequent reactions, G-protein and phospholipase C are activated, ultimately causing IP₃ to be released into the cytosol (see Dupont et al. [2016] for details). However, it is known that at least some cell membrane G-protein-coupled receptors (see Fig. 5) are also stretch-dependent (Mederos y Schnitzler et al., 2008; Yasuda et al., 2008; Pires et al., 2017; Erdogmus et al., 2019), so this pathway provides one possible mechanism of pressure dependence.

In our model, cytosolic [IP₃] has so far been taken to be constant. If it increases, more calcium ions are released from the SR, all other things being equal, and the C-clock is thus impelled to speed up. Assuming the stretch dependence of the G-protein-coupled receptors, we can therefore use cytosolic [IP₃] as a surrogate for distending pressure. Fig. 11 shows how the frequency of APs and the associated calcium flashes varies with cytosolic [IP₃] in our model. The figure also shows equivalent data for the model published by Hancock et al. (2022) and for the resynthesised version of the model by Imtiaz et al. (2007) that was implemented for that paper.

It is seen (Fig. 11, left panel) that the frequency range of IP₃ concentrations and frequencies over which the reinterpreted Imtiaz model provided pacemaker oscillation was very limited⁹. Pacemaking over a much wider range of [IP₃] was provided by the model of Hancock et al. (2022), but the range of pacemaker frequencies was still limited. In contrast, under control conditions as for Figs. 8 and 9, the present model ranges between a minimum frequency of 4.1/min at P = 0 µM and a maximum of 16.8/min at P = 2 µM. For comparison, the ex vivo data of Bertram et al. (2017) display a yet wider frequency range, from <2/min at zero transmural pressure, to some 28/min at 50 cmH₂O. But the frequency–[IP₃] pattern shown here, saturating at high [IP₃], is suggestively similar to the frequency–pressure pattern found experimentally, saturating at high distending pressure. This suggests that the mechanism of GPCR mechanosensitivity and resulting [IP₃] change alone may be sufficient to explain most of the observed extent of frequency change in response to pressure. The right panel shows that, for all IP₃ concentrations from 0.2 to 2 µM, in the new model, Ano1-KO reduces the AP frequency to between 24 and 30% of the corresponding frequency under control conditions. This finding exactly parallels what we found experimentally in mouse lymphatic collecting vessels (Zawieja et al., 2019): deleting Ano1 using three different smooth-muscle-specific Cre lines, a similar result was found in each, i.e., the basal pacemaking rate was reduced to <25% of control and pressure-induced chronotropy was blunted or abolished. The right panel also shows how AP frequency is reduced in IP₃R-KO in the new model: for all IP₃ concentrations from 0.2 to 0.8 µM (beyond which IP₃R-KO oscillation ceases unless other parameters are adjusted), AP frequency is reduced to between 34 and 41% of control. At [IP ₃] = 0, the control and IP₃R-KO traces become one. This matches fairly closely the experimental results (Zawieja et al, 2023), which show that contraction frequency in IP₃R1-KO vessels was between 17% and 27% of IP₃R1 fl/fl control counterparts across the pressure range that was studied. Specifically, at the distending pressure of 3 cmH₂O, IP₃R1-KO frequency was 26.9% of the control (4.45 versus 16.53/min).

In all the new-model simulations shown thus far, the Z-dependence of the gating variable h(V,Z) has been suppressed for the sake of simplicity by setting a value for K_h not reached during the Z oscillation. Finally, this dependence is now reintroduced by setting K_h = 0.6 µM. Fig. 12 shows how this change affects the waveforms of Z(t) and V(t) and the M-clock phase plane (the C-clock phase planes are essentially unchanged from Fig. 9).

The effects are especially marked during Ano1-KO. Insofar as Z(t) is concerned, the lack of cyclic regularity, earlier footnoted, is now highly evident. With Ano1-KO removing clock coupling, so that the M-clock frequency is untied from that of the C-clock, the two clocks oscillate at incommensurate rates. However, the size of the [Ca²⁺] flash is affected by APs, but to an extent that depends on the timing of the AP relative to the flash. At the beginning of the time window shown (long after numerical start-up transients have decayed), the AP occurs at the same time as the Z flash, causing the Z peak to be higher than usual. Over the subsequent few Z flashes, the peak height decays back to its unperturbed value. But the next AP occurs slightly later than the corresponding Z flash, the peak of which is accordingly barely higher than the preceding one. But the AP causes a late rebound in [Ca²⁺] so that the overall effect is to lengthen the flash.

With the reactivation of the Z-dependence of the L-type gate function h, the effects of the C-clock oscillation are also evident in the membrane potential trace for Ano1-KO, i.e., there is C- to M-clock coupling. The frequency of Ano1-KO APs is now even further reduced relative to both control conditions and the corresponding frequency seen in Fig. 8. As with the time course of intracellular [Ca²⁺], the interaction of the two clocks via oscillations of incommensurate frequency means that there is not a strictly repetitive V(t) waveform, but a calculation based on the timing of the two Ano1-KO AP peaks shown yields a pseudo-frequency of 1.33/min. The corresponding frequency for the (properly repetitive) Ano1-KO APs in Fig. 8 is 3.74/min.

Of even greater interest is the fact that, during the diastolic slow depolarisation phase, the membrane potential trace for Ano1-KO in Fig. 12 b now shows perturbations arising from the [Ca²⁺] oscillations. In other words, the operation of the C-clock can be observed via the M-clock’s output. This is a direct result of the effect of changing Z on the h-nullcline (Fig. 12 c), the left arm of which now moves during the cycle between h = 1 and h = 0.015. This opens the possibility of observing the same thing experimentally. Fig. 13 shows experimental data of membrane potential which include a qualitatively similar phenomenon occurring during the slow depolarisation phase of Ano1-KO, induced by administration of the pharmacological agent Ani9 to inhibit Ano1 instead of by genetic knock-out. As in the numerical trace of Fig. 12 b, diastolic oscillatory perturbations of increasing amplitude occur. While it is not possible to corroborate this speculation in the absence of a simultaneous recording of [Ca²⁺] flash frequency (not technically possible in our laboratory at this time), it is possible that the frequency of these otherwise unexplained experimental oscillations may represent the operation of the real C-clock in this case. If so, then as in the numerical simulation, the C-clock frequency in Ano1-KO greatly exceeds that of the M-clock and can be approximated from the trace as 11.5/min. The frequency of the lesser Ano1-KO Z peaks (between the two that are affected by APs) in Fig. 12 a is 12.0/min.

In this work, our overarching goal is to understand the key mechanisms and system effects that drive pacemaking in LMCs. To achieve this goal, we have focused on determining the underlying oscillators in a single LMC¹⁰ and how they are coupled. The focus of this interest has been on the M-clock and C-clock, and determining which of the clocks drives oscillation and how the clocks are coupled. An important part of this determination is studying how individual mechanisms such as the IP₃ receptor affect the individual clocks and the system as a whole.

We have made a number of changes to our existing LMC model (Hancock et al., 2022), which focus on modifying the mechanisms in the C-clock model. These changes had two motivations. First, the experimental data used in this paper involve the knock-out of a key C-clock mechanism (IP₃R), so developing an accurate model of the C-clock was important for these data. Second, our existing model was based on a model by Imtiaz et al. (2007) which added dependent M-clock components to an existing C-clock model by Dupont and Goldbeter (1993). The Dupont–Goldbeter model relied on store depletion to reduce IP₃R current, contrasting with more recent models that incorporate late IP₃R inhibition. The latter is now viewed as a crucial feature of calcium oscillation and so the Dupont–Goldbeter model is no longer widely used (Dupont et al., 2016). As such, the IP₃R mechanism in the model (DeYoung and Keizer 1992) was brought in line with more recent C-clock models in other cell types. The result of the IP₃R channel modification is a model structure that has similarities with the existing model, but with an extra gating variable for the IP₃ receptor. This addition takes the number of variables required for the C-clock model to three¹¹. Apart from this change, the ion channels represented are equivalent to the existing model, with only minor modifications to kinetics. However, the added complexity of the IP₃R model can dramatically change the behaviour of the C-clock. In particular, analysis of the self-driven oscillation in the C-clock is now centred on the intracellular calcium concentration Z and IP₃R gate state n, rather than Z and the intrastore calcium concentration Y as in our previously published model. This can be observed in the phase planes, where C-clock oscillation is best viewed in the (Z,n)-space phase plane rather than in (Z,Y)-space (see Fig. 9).

The proposed model is a system with dual coupled clocks, where both are active and both are important for the initiation of the AP. To obtain dual-clock-driven oscillation, we also modified the values of various parameters relative to those used by Hancock et al. (2022). The two clocks synchronise in frequency. Dual-clock coupled behaviour has been found essential to explain pacemaking in cardiac muscle cells (Maltsev and Lakatta, 2009, 2013), where experimental measurements are much more advanced, but has not previously been associated with LMCs. The modelling for cardiac cells involves highly detailed mechanistic models with large numbers of equations; see, e.g., Noble et al. (1998). In contrast, the model presented here is relatively low order, while obtaining qualitatively similar results. This lower order enables us to gain insight from analysis that would not be possible with a more complex model. In particular, it enabled the use of phase-plane analysis to understand the coupling between the two clocks. To our knowledge, this is the first time that a dual-active-clocks pacemaker has been elucidated using minimal models and the clear language of phase-plane analysis.

Despite the change to the C-clock model, the calcium concentration Z can still be used to represent the C- to M-clock coupling, which was an important property of our existing model (Hancock et al., 2022). Similarly, the L-type current I_L can still best represent the M–to–C-clock coupling. These properties enable each coupling direction to be represented by a single variable, simplifying the models and enabling insights to be drawn from the analysis. The alternative is to use h and V together to represent the coupling from the M- to the C-clock, which is unnecessarily complex and obscures the real mechanism.

We have demonstrated that coupling can occur between the clocks by multiple different means. The C-clock can affect the oscillations in the M-clock either via a constant level of calcium concentration Z setting the height of a static V-nullcline threshold for the operating point to surmount or via oscillations in the calcium concentration moving the AP threshold dynamically. Both methods act via the Ano1 mechanism, but their system effects are different. In the first case, the higher the steady-state Z, the lower the level of the gating variable h at which the M-clock self-triggers an AP. In this case, C-clock oscillation is not required for APs. In the second case, the rising phase of the cyclic Z oscillation progressively reduces the M-clock’s h-level trigger point for the AP. If this point falls below the current state of h, then an AP is initiated. Thus, the C-clock effectively triggers the M-clock’s AP rather than changing the self-trigger. The resulting upswing in oscillating calcium concentration then necessarily coincides with the AP in the M-clock, synchronising the two clocks.

The second of these cases represents a situation that goes beyond the traditional physiological view of the AP happening at such time as slow diastolic depolarisation causes the membrane potential to surpass a static voltage threshold. Instead, we show that, whenever calcium flashes and APs are synchronised, the barrier or threshold itself is dynamically falling away at the instant when the AP is triggered. This threshold-lowering effect, rather than what it does to the membrane potential, is the major influence of CICR-induced Ano1 current. Furthermore, we see that, in both cases, the M-clock threshold is far more effectively viewed as a level of the L-type gating variable h rather than as a level of the membrane potential V. Obviously, for the experimentalist, the only measurable quantity remains the membrane potential¹², and V(t) and h(t) are indeed linked monotonically albeit nonlinearly during the cycle phase leading up to the next AP. But, expressed in terms of membrane potential, these thresholds and threshold motions are a small fraction of the total excursion, whereas in terms of gating variable h they are order-1 quantities. Thus, the closest physiologically observable quantity, slow diastolic depolarisation, while important (and reproduced in the model), is determined by many complex coupled factors. There is also coupling from the C-clock to the M-clock via inactivation of the L-type calcium channel, but this mechanism has no effect on the triggering in this model.

The M-clock coupling to the C-clock occurs via a large inrush of calcium ions during the AP via the L-type calcium channel. If the C-clock is not oscillating independently then this periodic inrushing current causes the C-clock to oscillate passively in response. Alternatively, if the C-clock oscillations trigger an AP in the M-clock then the inrush of current amplifies the increase in cytosolic calcium concentration. This amplification causes a [Ca²⁺] flash of larger amplitude and longer duration in the C-clock. Overall, the C-clock can initiate an earlier AP in the M-clock, while the M-clock causes amplification of the calcium uptake. Combined, this bidirectional coupling can cause the two to synchronise in frequency.

Comparing the performance of the model with the experimental findings in the accompanying paper (Zawieja et al, 2023), the model performs particularly well with respect to the AP frequencies observed in control, Ano1-KO, and IP3R-KO conditions. However, in other respects, the agreement is less good. Experimentally, the AP in IP₃R-KO has a longer systolic plateau than the control one, and the total peak-to-trough waveform excursion is less. In contrast, the simulated AP in IP₃R-KO has the same form and brevity as that both observed and simulated for Ano1-KO, albeit with a somewhat lower peak. Furthermore, the model predicts that the waveform excursion of the calcium flash is reduced by a factor of almost six in IP₃R-KO, whereas unchanged or slightly increased excursion is observed. In fact, with reference to the very few calibrated measurements of calcium flash size in lymphatics (Zawieja et al., 1999; Stolarz et al., 2019), it is likely that the model predicts an excessively large calcium flash under control conditions, and possibly under Ano1-KO (there are no data). However, on the plus side, the model predicts a shape for the calcium flash which is in noticeably good agreement with that observed. Furthermore, it seems possible on the basis of our (M.J. Davis, S.D. Zawieja) recent and unpublished experiments that the observed maintenance of calcium flash excursion may relate to the recruitment during IP₃R1-KO of normally dormant ryanodine receptors (RyRs)¹³ to sustain store Ca²⁺ release. The model presented here does not contain any representation of RyRs; thus, the IP₃R-KO that we model leaves only store leakage as a mechanism for Ca²⁺ release.

The presence or absence of an AP plateau is intimately linked to the shape of the V-nullcline in (V,h)-space and how this shape changes during the cycle. As shown in Fig. 9 a and Fig. 12 c, under control conditions, the V-nullcline moves between a curve having two arms of positive slope joined by a negative-slope region and a simple curve of continuously positive slope, according to the instantaneous value of Z. The tripartite shape corresponds to low values of Z and the monotonically positive-slope shape to high values. If the V-nullcline still has this latter shape when the trajectory of the decaying AP approaches the nullcline intersection at or near h = 0, the trajectory cannot proceed past this point until Z falls sufficiently that the V-nullcline reverts to the tripartite form. When this occurs, the operating point moves to attain the left-hand positive-slope arm of the V-nullcline at a rate that depends on its V-distance from the nullcline; this is the steeply descending final portion of the AP. Under Ano1-KO, which amounts to tying Z to zero insofar as the V-nullcline is concerned, this nullcline is always tripartite, and so a systolic plateau is impossible; this explains the extreme brevity of the Ano1-KO AP as both observed and simulated.

In the case of the control AP here, the gradual increase in descent rate from the sloping plateau occurs while the V-nullcline is still of monotonically positive slope, but as Z reduces, the low-h part of the V-nullcline moves leftwards, with the operating point following. The V-nullcline becomes tripartite again when V reaches about −50 mV. The operating point is then attracted to the left-hand positive-slope arm of the V-nullcline, but since the V-distance between the operating point and the V-nullcline remains relatively small, the rest of the descent occurs at a rate scarcely greater than before. See Video 1, presented as supplemental information.

Again referring to Fig. 9 a or Fig. 12 c, we see that, in IP₃R-KO as simulated here, the V-nullcline barely achieves the high-Z positive-slope-only form at all, and consequently reverts to the tripartite shape too quickly to delay the end of the AP. What would be needed for the model to emulate the relatively prolonged systolic plateau of the observed IP₃R-KO AP is for the V-nullcline to spend relatively long in the positive-slope-only form, trapping the operating point and preventing the end of the AP. This in turn implies excursion to higher values of Z, as during the simulated control calcium flash (Fig. 10), but the model as at present constituted does not predict such values during IP₃R-KO. See Video 2 as supplemental information.

More generally, the model we present is the first to deal with IP₃R-KO in the context of lymphatic pacemaking. Although greatly simplified relative to the in vivo situation, where a number of additional channel types have been shown to exist¹⁴, some of which have importance only when normally dominant ones are inhibited, it nevertheless has a rather large number of settable parameters. In future work, we plan to assess whether alternative tunings of these parameters address any of the issues identified above or whether more fundamental modifications to the model are needed. These could include the addition of further mechanisms of M- and C-clock coupling that have been omitted here; an example is the dependence of transmembrane calcium entry on cytosolic [Ca²⁺]. Obviously, not all discrepancies between experiment and model can be addressed in a simplified model. However, many of the data that would be needed for a mechanistically detailed model have simply not yet been measured for the LMC. Whole topics have yet to be contemplated experimentally, such as the possible dependence of LMC pacemaker parameter values on matters affecting the mechanical load that the LMC experiences, including gravity (Olszewski and Engeset, 1980; Gashev et al., 2006; Li et al., 2022), location relative to the high hydraulic resistance presented by a lymph node (Browse et al., 1984), etc. Our approach in this paper and previously has been to add mechanistic detail as required based on the given experimental data. The minimal model presented here is a useful starting point for more detailed mechanistic models, as well as providing important insight in its own right.

Modelling is an important tool to help understand pacemaking in LMCs, and consequently, it is important for understanding lymphatic system dysfunction. In this paper, we have extended the state-of-the-art to pacemaker modelling for LMCs. We have introduced a model of LMCs where the M-clock and C-clock are synchronised and both drive the APs and associated calcium concentration surges. This change resulted from changes to the mechanism of the IP₃R channel in the C-clock and enables the model to match experimental data for IP₃R knock-out. We have also provided phase-plane analysis explaining the result. The model, results and analysis are important for advancing comprehension of pacemaking and its dysfunction in LMCs.

Animation Video 1 and Video 2 are available for download. Matlab code for the solution of Eqs. 1–19 has been deposited in the open access Sydney eScholarship Repository at https://hdl.handle.net/2123/31688 (doi: 10.25910/vjxw-j379).

Joseph A. Mindell served as editor.

This research, E.J. Hancock, and in part C.D. Bertram, were supported by National Institutes of Health (NIH) grant R01-HL-122578 to M.J. Davis. S.D. Zawieja acknowledges NIH grant R00-HL-143198.

Author contributions: The project was jointly conceived by C.D. Bertram and M.J. Davis. The numerical models were developed by E.J. Hancock and C.D. Bertram, with help from C. Macaskill. The experimental data were recorded by S.D. Zawieja, with help from M.J. Davis. E.J. Hancock and C.D. Bertram contributed equally to the writing of the paper. All authors reviewed and contributed to the manuscript and all approved the final version.