Global asymptotic stability of the active disassembly model of flagellar length control

As described in the Introduction, the process of IFT involves the transport of various biomolecules between the flagellar base and tip, and is believed to be closely linked with the dynamics of flagellar length. In the development of the mathematical model described next, we conceptually simplify this process by capturing only the essential biomolecular constituents for length control. In particular, we restrict attention to the tubulin that makes up the axoneme structure, the molecular motors that power transport along the axoneme, and the depolymerase conjectured to be necessary for length control. We assume that both assembly and disassembly are in continual competition so that the overall growth rate is the net result of these processes. Next, we describe the assembly and disassembly rates that give rise to differential equations for the flagellar lengths.

Let $L_{1}(t)$ and $L_{2}(t)$ denote the lengths of two flagella on a single Chlamydomonas organism at a particular time $t$. The flagellum is composed of tubulin filaments, and we assume there is a total amount $T$ of tubulin available. The units of tubulin are given in microns, so that the pool size and flagellar lengths have the same units and the amount $T_{f}$ of free tubulin at the flagellar base is given by $T_{f}=T-L_{1}-L_{2}$.

We assume there are a total number $M$ of anterograde molecular motors, of which a number $M_{f}$ are free in the basal pool while the others are undergoing IFT. In line with the pattern of motion found experimentally for kinesin [18], the motors move ballistically in the anterograde direction with velocity $v$ and diffusively in the retrograde direction with diffusion coefficient $D$.

The injection flux $J$ of IFT particles into the flagella satisfies mass action kinetics, so that

There is a separation of timescales between the motion of IFT particles, which traverse the flagellum in several seconds, and the flagellar length dynamics, which varies over a timescale of several minutes [19]. As shown in Appendix A, upon combining (1) with conservation of protein number, this separation of timescales leads to a quasi-steady state approximation for the flux:

The flux $J$ represents the number of IFT particles injected per second into the flagella, and these particles are assumed to carry an amount of tubulin proportional to the free tubulin in the basal pool. Therefore, the overall speed of assembly is equal to $\gamma J\left(T-L_{1}-L_{2}\right)$, where $\gamma$ is a proportionality constant that represents the fraction of tubulin arriving at the flagellar tip that is incorporated into the flagellum.

Disassembly is taken to be a linear function of the depolymerase concentration at the flagellar tip, i.e. each depolymerizing enzyme located near the flagellar tip removes a constant number of tubulin subunits per second. This implies that for the $i$^th flagellum for $i=1,\,2$, the speed of disassembly is equal to $d_{0}+d_{1}\frac{JL_{i}}{D}$, where $d_{0}$ and $d_{1}$ are the coefficients of the first-order Taylor series approximation. As explained in detail in Appendix A, this follows from the assumption that there is a depolymerase which has the same ballistic-to-diffusive pattern of motion as kinesin: the disassembly rate equals $d_{0}+d_{1}c_{d}(L_{i})$, where $c_{d}(L_{i})$ is the concentration of depolymerase at the flagellar tip.

Setting the growth rate to be the difference between assembly and disassembly speeds leads to a system of coupled nonlinear differential equations:

where $J\equiv J(L_{1},L_{2})$ satisfies

Next, we turn to the mathematical analysis of the dynamical equations. To simplify notation, we drop tildes and use the dimensionless formulation for the remainder of the article. We will first consider the case of $N=2$ flagella, and later on show that the key results generalize to arbitrary flagellar number. The active disassembly model eqs. 6 and 7 is not a gradient system, i.e. it cannot be written as the gradient of an energy function $E(L_{1},L_{2})$. This is proven by contradiction. Suppose there is an energy function such that $L_{1}^{\prime}=\partial E/\partial L_{1}$ and $L_{2}^{\prime}=\partial E/\partial L_{2}$. Then, by equality of mixed partials

i.e.

However, the left-hand side of the equality is equal to

whereas the right-hand side is equal to

Because $J(L_{1},L_{2})=J(L_{2},L_{1})$, we have $\partial J/\partial L_{1}=\partial J/\partial L_{2}$. Assuming $\pi_{1}\neq 0$, it follows that there is equality if and only if $L_{1}=L_{2}$. We have thereby reached the contradiction since $L_{1}\neq L_{2}$ in general. It follows that no energy function exists.

However, using a technique applicable to many nonlinear dynamical systems, we may construct a Lyapunov function to show that the steady state is globally asymptotically stable. First, we specify the subspace $\Omega$ with $(L_{1},L_{2})\in\Omega\subset\mathbb{R}^{2}$ over which the model is defined. Because we must have $L_{1},\,L_{2}\geq 0$ for the variables to have physical meaning as lengths, $\Omega$ is contained within the first quadrant. Moreover, the sum of lengths cannot exceed the total pool size so that $L_{1}+L_{2}\leq 1$. To summarize, $\Omega$ is defined to be the triangular region in the first quadrant bound by the lines $L_{1}=0$, $L_{2}=0$, and $L_{1}+L_{2}\leq 1$.

We define the deviations from steady-state $\Delta L_{1}:=L_{1}-L_{\text{ss}}$ and $\Delta L_{2}:=L_{2}-L_{\text{ss}}$ on $\Omega$ and express the dynamics in terms of $\Delta L_{1}$ and $\Delta L_{2}$:

or equivalently

By symmetry, the steady-state lengths $L_{1,\text{ss}}$ and $L_{2,\text{ss}}$ are equal. We denote the common steady-state length by $L_{\text{ss}}$. It is determined by setting the derivative equal to zero in (21) or (22), which yields the identity

where

(Eq. (23) yields a quadratic equation for $L_{\text{ss}}$, which has a unique positive root given by (9).)

We next use the steady-state relationship (23) to rewrite eqs. 21 and 22 in a different form. Before we proceed, we first note the identities given by

where the first line follows from adding and subtracting $J_{\text{ss}}\left(1-2L_{\text{ss}}\right)$ and the second line follows from adding and subtracting $\pi_{1}J_{\text{ss}}L_{\text{ss}}$. Plugging the above expressions into eqs. 21 and 22 and using the steady-state relation (23) yields

Adding and subtracting the above equations yields

It follows that

where we have used $\mathrm{d}L_{i}/\mathrm{d}t=\mathrm{d}(L_{\text{ss}}+\Delta L_{i})/\mathrm{d}t=\mathrm{d}\Delta L_{i}/\Delta t$ for $i=1,\,2$ and

We use eqs. 31 and 32 to form a Lyapunov function for the dynamics. Let

Next, we consider whether $\Omega$ is a trapping region closed to the dynamics. That is, we must consider whether any trajectories originating in $\Omega$ may eventually leave $\Omega$. First, we show that no trajectories cross the $L_{1}+L_{2}=1$ boundary. To show that all trajectories move inward from this boundary, we calculate

This proves that trajectories flow inward cannot leave $\Omega$ across the $L_{1}+L_{2}=1$ boundary. See Figure 4 for an illustration of $\Omega$ as well as the slope field and contours of $\phi(\Delta L_{1},\Delta L_{2})$.

In addition, for $\Omega$ to be a trapping region, no trajectories should escape the first quadrant. In terms of $\Omega$, this means that no trajectories cross the $L_{1}$ axis on the interval $0<L_{1}\leq 1$, and no trajectories cross the $L_{2}$ axis on the interval $0<L_{2}\leq 1$. However, the model eqs. 6 and 7 does not satisfy this requirement. To the contrary,

and upon solving the quadratic equation resulting from $J(T-L_{1})-\pi_{0}=0$ one finds that there is a unique positive root

Therefore, to prevent trajectories that cross the $L_{1}$-axis resulting in negative lengths, we construct a hybrid system by setting $L_{2}^{\prime}=0$ when $L_{1}>L_{\text{crit}}$ on $L_{2}=0$. This results in trajectories that flow down to the $L_{1}$-axis, move horizontally according to this modified dynamics with $L_{2}^{\prime}=0$ until they reach $(L_{1},L_{2})=(L_{\text{crit}},0)$, at which point they rejoin the full dynamics and move toward steady state. The dynamics are modified symmetrically to prevent trajectories from crossing the $L_{2}$-axis as well.

To characterize the trapping region of the full dynamics, we therefore define the trajectory $\Gamma_{1}(t)$ going backward in time from the point $(L_{\text{crit}},0)$, as defined by

up until it hits the diagonal $L_{1}+L_{2}=1$. We define another trajectory $\Gamma_{2}(t)$ similarly as the one going backward in time from the point $(0,L_{\text{crit}})$. The points from which the trajectories $\Gamma_{1}(t)$ and $\Gamma_{2}(t)$ originate are illustrated by red diamonds and the trajectories are shown as red curves in Fig. 4.

As a consequence of the modified dynamics, trajectories originating in the corner regions of phase space bounded by $\Gamma_{1}(t)$ and $\Gamma_{2}(t)$ consist of three distinct phases: first, there is flow down to the $L_{1}$ or $L_{2}$-axis, followed next by horizontal or vertical motion to the point $(L_{\text{crit}},0)$ or $(0,L_{\text{crit}})$, respectively, and finally leading to recovery by the full dynamics back to steady state. One interesting consequence of this behavior is its interpretation in the context of severing experiments. We find that the typical severing protocol, in which the two flagella first reach steady-state before severing takes place, does not lead to states within the corner region (purple trajectory in Fig. 5).

Fig. 5 shows a representative simulation of a severing experiment. The code used to generate this simulation and the corresponding figures is included in our open source code repository located at https://github.com/thomasgfai/StabilityFlagellarLengths.

If the system is initially out of steady-state, however, severing may result in a state outside the trapping region. This yields two phases of flow after severing: first, a horizontal flow to the accumulation point $L_{\text{crit}}$, and second a flow by the full dynamics to reach the global steady-state (green trajectory in Fig. 5). As shown in the inset of Fig. 5(c), this first phase of horizontal flow leads to an apparent delay in the response of the cut flagellum.

As mentioned above, the corner regions that lead to negative lengths in the original active disassembly model are not encountered in typical biological scenarios. One may nevertheless ask how a model derived from physical principles could lead to negative lengths. In our case, this behavior arises because of the basal disassembly rate $d_{0}$ that appears in the model, which implies a nonzero rate of disassembly even when the flagellum is very small. In the corner regions of the phase space, this basal term is greater than the rate of assembly because the tubulin pool is nearly depleted. This leads to a negative growth rate at zero length. Note that this issue does not arise in agent-based models in which each tubulin monomer is treated explicitly, since in that case disassembly only occurs when the flagellum contains at least one monomer. Within the ODE model, the modified dynamics ensure that disassembly only occurs when there are tubulin monomers remaining on the flagellum. Alternatively, one could rectify this issue by defining $d_{0}=d_{0}(L)$ to be a step function of length that turns off once a flagellum is sufficiently small.