Understanding the mechanisms that regulate the dimensions of subcellular compartments to match cell volume and metabolic needs remains a fundamental challenge in cell biology. Current evidence suggests that organelle dimensions are governed by a combination of physical limitations and autonomous self-organization [1]. Researchers have broadly categorized these regulatory strategies into several main types, notably quantal synthesis, molecular rulers, and dynamic balance [2].

The principle of dynamic balance dictates that the continuous turnover of structural components can yield a stable size, provided that the rates of building and dismantling the structure eventually equalize. Crucially, if either of these kinetic rates is intrinsically tied to the structure’s current dimensions, the system naturally contracts to a steady-state size. This is exemplified by microtubules in eukaryotic flagella, which experience ongoing turnover at their distal ends. A uniform breakdown rate is counteracted by an elongation rate that scales with length—a consequence of a finite pool of motor proteins delivering building blocks from the soma. This equilibrium ultimately dictates the final length of the flagellum [3–5].

Similar equilibrium models have been proposed for actin-rich protrusions like inner ear stereocilia [6]. In these systems, actin filaments undergo continuous treadmilling: they are degraded at their proximal base while growing at their distal tips. Because the addition of new material relies on the passive diffusion of actin monomers to the far end, the growth rate inherently decreases as the structure lengthens. A comparable diffusion-limited growth model has also been hypothesized to govern hook length in bacterial flagella [7].

Conversely, yeast cells utilize a distinct balancing strategy to manage microtubule length. In this case, specific kinesins travel processively to the ends of microtubules to actively promote depolymerization. Because longer microtubules offer a larger surface area to capture these depolymerizing kinesins from the surrounding cytosol, the rate of structural breakdown increases directly with length. Pairing this length-dependent degradation with a constant rate of polymerization ensures the microtubule reaches a definitive steady-state size [8–13]. Other related models of microtubule regulation focus on adjusting their dynamic balance, specifically through the targeted control of catastrophe frequencies [14, 15].

Here, we develop a multiscale model that focuses on the dynamic balance of microtubules. Microtubules are directionally-polarized filaments with biophysically distinguishable (+) and (–) ends [16]. The principal functions of microtubules are twofold: (1) maintain structural integrity of cell shape, and (2) aid in transport of vesicles and organelles throughout the interior of a cell. They also play a key role in cell motility during mitosis [16, 17].

There are two broad categories of microtubules: centrosomal and non-centrosomal. Centrosomal microtubules emerge from the pericentriolar material surrounding the centrioles of the centrosome. Their (–) ends are anchored to the centrosome and are thus not subject to dynamic growth and shrinking. The (+) ends extend out radially in the cytoplasm. This type of microtubule is prominent in most dividing animal cells. Non-centrosomal microtubules emanate from other microtubule organizing centers, such as the Golgi apparatus, nuclear envelope, or even from other microtubules via branching nucleation or severing. Both (+) and (–) ends are free and subject to growth and shrinking. They are prominent in polarized cells, such as neurons and other specialized cells, and support local organization rather than global radial structure [18, 19]. In this work, we focus on non-centrosomal microtubules, as they exhibit a variety of dynamical behaviors, such as growth, shrinking, and treadmilling.

Transport and delivery of vesicles and organelles is carried out by molecular motors, such as kinesin and dynein. Although the precise mechanism by which these motors collectively transport vesicles is debated, it is known that the polarity at a given end of the microtubule dictates what kind of molecular motor will travel along the microtubule in a given direction. For example, kinesin motors generally walk in the (+) direction along microtubules, whereas dynein motors tend to walk in the (–) direction. The transport and delivery process enabled by molecular motors is crucial for cellular and organism health. Breakdown in these processes in neurons has been implicated in neurodegenerative diseases such as Alzheimer’s and Huntington’s [20–22].

The question of how microtubule length is controlled therefore reflects an interesting balance between distinct effects. On the one hand, microtubules are continuously switching between states of growth and shrinking to facilitate cytoskeleton rebuilding and robustness to injury; on the other hand, they are required to maintain a semblance of constant length to facilitate robust vesicle transport throughout a cell. Understanding the mechanisms that generate specific microtubule length distributions is therefore of interest to both theorists and experimentalists.

Situations like this are precisely where mathematical models can play a key role in advancing biological understanding—namely to theorize outcomes for hypothesized mechanisms underlying a particular biological phenomenon. Indeed, a plethora of mathematical models exist that investigate some aspect of microtubule length control. However, these models tend to focus on the dynamics at one end of centrosomal microtubules (the (+) end) and are typically framed in the context of in vitro experiments. Several biological and mathematical models have been specifically proposed to address an important growth control mechanism for microtubules, where filament assembly is limited by length- or age-driven microtubule disassembly [23, 24]. While an exact molecular basis for this length-dependent catastrophe mechanism has not been fully established, experiments have supported several hypotheses. One is the ‘antenna’ mechanism, where certain types of depolymerizing kinesin motors accumulate at microtubule (+) ends and promote catastrophe [9, 25]. Another is a mechanism that proposes that, as microtubules grow, their tips become more ragged due to the variable protofilament lengths and are thus destabilized [26, 27].

Here, we develop a model of microtubule length dynamics at both ends, with parameters informed by in vivo data, coupled with a model of molecular motor kinetics, parametrized using in vitro measurements. We assume a length regulation paradigm where shrinkage is dependent on the density of depolymerizing molecular motors attached to microtubules. We explicitly model the dynamics of two types of molecular motors in the kinesin family and link their kinetics to overall microtubule length. We use a 1D macroscopic parabolic partial differential equation to model kinesin motor density dynamics. Such an equation can be derived in a principled way from a more detailed biophysical model formulated in terms of a stochastic hybrid system where a particle switches according to a Markov process between states where it moves ballistically in the (+) direction, ballistically in the (–) direction, or remains stationary [28, 29]. The domain of the advection-diffusion equations describes the length of the microtubule and thus changes through time (figure 1).

Figure 1. Schematic showing the bidirectional motor-based depolymerization of microtubules. These dynamics provide the foundation of our theoretical model.

Download figure:

Our approach is reminiscent of the earlier model [25], which proposed the first mathematical formulation for the antenna model of MT depolymerization due to motor accumulation. This model tracked the density of depolymerizing Kip3p (kinesin-8) motors along a microtubule and showed that different motor binding models could not be distinguished. As they were interested in different phases of MT depolymerization, this study did not model individual microtubule end dynamics and instead focused on reproducing the decrease in microtubule length due to the depolymerizing motors. In addition, our approach is most similar to the model in [15], which also describes the density of depolymerizing motors along a single microtubule. As in [25], they consider a single class of motors that depolymerizes microtubules from the (+) end, and do not distinguish between bound and unbound motor states. This study then couples the microtubule length to the density of motors at the last site on the filament. A more complex version of the model they propose includes dynamic instability at the (+) end and incorporates a catastrophe rate that varies with length. By modulating the catastrophe frequency as a function of the density or the flux of depolymerizing motors, [15] makes analytical approximations of the mean MT length in these two theoretical scenarios. However, this study does not intend to model the dynamics of non-centrosomal microtubules, which undergo growth and shrinking at both filament ends, and does not distinguish between the contribution of depolymerizing motors and other mechanisms to the length-dependent catastrophe of microtubules.

Our model therefore goes beyond these prior studies by incorporating dynamics at (+) and (–) ends of microtubules, informing parameters from in vivo and in vitro data whenever possible, and assessing the contribution of motor-driven and non-motor-driven length-dependent disassembly. Exploring the parameter space of this model allows us to investigate the qualitative differences of the system’s behavior in various regimes. In the adiabatic regime, where motor dynamics are assumed to occur faster than microtubule length dynamics, we derive a reduced model that is amenable to analytic techniques. Using homotopy continuation, we analyze equilibria of the reduced model and describe regions of parameter space where different types of microtubule behavior (e.g. growth, disassembly, treadmilling) can be observed. Specifically, our reduced model provides a simple relationship between model parameters to distinguish between growth and disassembly in the adiabatic regime.

In the unrestricted parameter regime, we perform numerical simulations of the model on a growing domain and are also able to identify qualitatively distinct microtubule behaviors in various parameter regions. We are thus able to assess the impact of motor kinetics and different length-dependent disassembly mechanisms on microtubule dynamics. The main contribution of our work is thus the development of a flexible model that can shed light on the regimes of biophysical parameters that characterize microtubule length regulation. Such quantitative tools are especially useful given that ascertaining characteristics of in vivo behavior of microtubules remains challenging.

As mentioned above, one of the ways in which microtubules achieve length regulation is through a length-dependent catastrophe mechanism, which has been observed both in vitro and in vivo [24]. Microtubules have been found to experience increased catastrophe frequency as they grow, and a leading mechanism for this phenomenon is the so-called ‘antenna’ mechanism, where kinesin motor proteins walk to microtubule ends and stimulate catastrophe. Studies on the kinesin-8 motor (Kip3P in the budding yeast) found that more motors accumulate at the (+) ends of longer microtubules, which in turn leads to further destabilization of these ends [9, 25]. Kinesin-8 was also found to be a slow, highly processive motor in [9], leading to the analogy of the microtubule filament to an antenna that directs Kip3p molecules to its (+) end and thus leads to a concentration gradient of the motor along the microtubule. These findings hold in vitro and in vivo and have positioned kinesin-8 as a length-dependent depolymerase that controls microtubule length from the (+) ends [9].

On the other hand, kinesin-13 motor proteins (MCAK and Klp10A) were found to depolymerize microtubules from both ends, though depolymerization was faster from the (–) ends [9, 30]. Unlike kinesin-8, kinesin-13 is known to start depolymerization by bending the terminal subunits of microtubules, leading to their loss of lateral interactions [25, 30]. The MCAK kinesin-13 protein is known to diffuse in a 1D random walk along the microtubule filament, as opposed to using directed motion like kinesin-8 [25, 31]. These studies have thus established kinesin-13 as a depolymerase that uses rapid diffusion to get to both ends of the microtubule, and which is able to cover short distances rapidly [31]. The protein Patronin was found to help (–) ends resist to complete depolymerization by kinesin-13, and the regulation of dynamic instability by these motors was thus found to be responsible for controlling mitotic spindle length in vivo [30].

Motivated by the biological evidence outlined in section 2, we develop a model that describes the motion of both (+) and (–)-end depolymerizing motors along microtubules. Neurons are long structures and, since we are interested in the dynamics of microtubules along such structures, we model space as one-dimensional, implicitly assuming all significant dynamics occur along the length of the neuron. Our goal is to model the impact of (+) and (–)-end motors on microtubule length. This setting yields the following coupled system of partial differential equations:

where is the density of MT-bound (moving) (+)-end motors, is the density of unbound (+)-end motors, is the density of MT-bound (moving) (–)-end motors, and is the density of unbound (–)-end motors. The motors move with speeds and diffuse with diffusion coefficients along microtubules. Motors diffuse with diffusion coefficient D when unbound from microtubules. The binding rates of the motors to microtubules are and their unbinding rates are given by . We denote the MT length by L, which extends from to .

The advection-diffusion equations in equations (1) and (3) are effective mean-field reduction equations of a more detailed stochastic biophysical model describing the transport of a population of motors of a given type, which randomly switch between a motile state moving ballistically in the (+) direction, a motile state moving ballistically in the (–) direction, and a stationary state. Assuming the transitions between states are fast relative to the transport dynamics, one can employ a quasi-steady state approximation upon the associated Chapman–Kolmogoroff equations to obtain an effective partial differential equation for the dynamics. The diffusivity describes the modification to motor-cargo dynamics arising from the intermittent stationary state [28, 32].

We need to specify boundary conditions at the (+) and (–) ends for both the motor densities and for the MT ends. At the (–) end, we model a constant source of moving and diffusing (+)-end motors, as well as of diffusing (–)-end motors:

while at the (+) end, we model a constant source of moving and diffusing (–)-end motors and of diffusing (+) end motors as follows:

The boundary conditions for the moving (+)-end motors at the right edge of the domain and for the moving (–)-end motors at the left end of the domain assert that the Fickian flux of motors departing the domain is zero. The justification for imputing zero Fickian flux upon the motors is again derived from the underlying stochastic model of the advection-diffusion equations. In the stochastic model, motion across the boundary of the domain occurs in the motile, ballistic states [28, 33]. Thus, we allow for convective flux but set Fickian flux to naught.

We assume that the dynamics of the (+) and (–) ends of the MTs is given by:

where the first term corresponds to intrinsic growth of the MTs, the second term denotes intrinsic length-dependent MT disassembly, and the third term models MT shrinking driven by the (+) and (–) end motor populations, respectively. These equations do not explicitly model the switch from growth to shrinking of MTs, but rather the overall dynamics at the (+) and (–) ends. We therefore incorporate the length-dependent shrinking mechanism in the second term of these equations.

Equations (1) and (7) describe a model that can be derived from first principles, and solving them numerically provides insight into how solutions behave. However, they are analytically difficult to parse. We thus invoke some approximations to obtain analytic insights from the model. Specifically, we reduce the model in equations (1) and (7) into a coupled set of ODEs to explore the parameter space of the model and examine the biophysical insights that emerge.

We assume that switching occurs on a faster timescale than the advective and diffusive processes so that and are larger than and , where is the equilibrium microtubule length attained in the absence of motor influence. This has been observed experimentally in the peripheral sensory neurons of Drosophila [34]. We make this explicit with the mapping , where . The effective dynamics then manifest as a perturbation away from the stationary measure of the process governing the switching dynamics. We show the derivation explicitly for (+) moving motors and the derivation for (–) moving motors is nearly identical.

We rewrite the dynamics of the (+) moving motors as

where

and

The co-kernel of A is spanned by the vector and the kernel of A is spanned by the vector

so that . Indeed u is the eigenvector corresponding to the zero eigenvalue of A and represents the stationary density of a Markov process formed by the motor switching rates.

Let and , so that q is proportional to the component of p in the co-kernel of A and w is in the orthogonal complement. Applying vT to both sides of equation (8) gives

Substituting into equation (8) and using equation (9) and the fact that , we obtain

where is the identity matrix. We introduce the expansion

and substitute it into equation (10). Collecting terms gives , since we require . In this case, higher order terms are not needed to extract a reasonable effective equation, so we do not compute higher order terms.

Substituting the leading order term for w () into equation (9) gives the effective advection-diffusion equation for , the total density of (+) moving motors at position x at time t:

where

We refer to as the effective velocity of the (+) motors, and to as the effective diffusivity of the (+) motors. After performing a similar approximation on the model of (–) moving motors, we obtain the following system of parabolic equations describing motor dynamics of (+) moving motors and (–) moving motors :

where

These are coupled with the boundary conditions:

We note that the Robin conditions arise from the fast switching between Neumann and Dirichlet conditions [35].

We make the further assumption that effective motor velocities are greater than the growth rates of the microtubules, to ensure that motors can reach the end of the microtubules and induce motor-dependent shrinking before the microtubule reaches its basal equilibrium length. Implementing this approximation will allow us to investigate and make predictions regarding those regimes where such assumptions hold. Acquiring such data about an application system is challenging, however this is exactly where mathematical modeling is especially useful. We can leverage it to make predictions about what could happen in the particular parameter regimes that correspond to the approximation assumptions here.

We assume the motor dynamics reach equilibrium before studying microtubule length dynamics. The equilibrium solutions are readily obtainable from equations (12) and (13), yielding

Substituting these into equations (7) gives

Defining yields the 1D flow for the length of the microtubule:

Equation (15) yields a unique equilibrium microtubule length provided This is straightforward to see by rewriting equation (15) as

where

Since and , , if (i.e. if ), then the Intermediate Value Theorem guarantees that f and g intersect for some L > 0. The intersection point is the unique microtubule length, (see figure 2).

Figure 2. Dynamics of equation (15) in different parameter regimes. (A) If , the model predicts the existence of a unique equilibrium microtubule length, ; (B) If , then the microtubule shrinks.

Download figure:

An interesting aspect of the reduced model is the prediction that the basal microtubule shrinking rates, , do not affect the existence of an equilibrium length but only its value (see figure 3). Only the basal growth rates of the microtubules and the motor-dependent shrinking rates are responsible for establishing whether a microtubule reaches an equilibrium length or shrinks to zero.

Figure 3. Dependence of upon various parameters in the reduced model in equation (15).

Download figure:

Treadmilling. A generic equilibrium resulting from equation (15) results in a treadmilling equilibrium, meaning but . This is another key prediction of the model: equilibrium microtubule lengths are effectively only maintained dynamically with growing and shrinking tips rather than fixed tips.

Equilibrium tip positions can be obtained in the reduced model, (15), but they occur in very constrained regimes of parameter space. To have an equilibrium point of equations (14), we must have and . Consider the nullclines of the system obtained by setting each equation individually equal to zero. The nullclines and