Abstract
During cell division, the mitotic spindle moves dynamically through the cell to position the chromosomes and determine the ultimate spatial position of the two daughter cells. These movements have been attributed to the action of cortical force generators which pull on the astral microtubules to position the spindle, as well as pushing events by these same microtubules against the cell cortex and plasma membrane. Attachment and detachment of cortical force generators working antagonistically against centring forces of microtubules have been modelled previously (Grill et al. in Phys Rev Lett 94:108104, 2005) via stochastic simulations and mean-field Fokker–Planck equations (describing random motion of force generators) to predict oscillations of a spindle pole in one spatial dimension. Using systematic asymptotic methods, we reduce the Fokker–Planck system to a set of ordinary differential equations (ODEs), consistent with a set proposed by Grill et al., which can provide accurate predictions of the conditions for the Fokker–Planck system to exhibit oscillations. In the limit of small restoring forces, we derive an algebraic prediction of the amplitude of spindle-pole oscillations and demonstrate the relaxation structure of nonlinear oscillations. We also show how noise-induced oscillations can arise in stochastic simulations for conditions in which the mean-field Fokker–Planck system predicts stability, but for which the period can be estimated directly by the ODE model and the amplitude by a related stochastic differential equation that incorporates random binding kinetics.
Similar content being viewed by others
Avoid common mistakes on your manuscript.
1 Introduction
Embryos develop, on the most basic level, as a result of one cell dividing into two cells. In a tissue, the orientation of cell division is an important factor in determining either the outcome for the daughter cells (e.g. cell fate due to distribution of intracellular components or the daughter cell local environment) or the tissue as a whole (e.g. building tissue and organ architecture by tissue stratification or spreading) (Bergstralh and St Johnston 2014; Morin and Bellaïche 2011). Cell division orientation is determined by the mitotic spindle, the large microtubule-based structure which forms in the cell and segregates genetic material into two discrete daughter cells (Karsenti and Vernos 2001; Mitchison and Salmon 2001). Prior to anaphase, where the chromosomes are pulled apart to opposite ends of the cell, the mitotic spindle is positioned translationally and rotationally (Fig. 1).
Key to spindle positioning is the pushing and pulling of astral microtubules which reach between the spindle pole and the cell cortex (Dogterom et al. 2005; Burakov et al. 2003; Zhu et al. 2010; Pecreaux et al. 2006; Howard 2006; Okumura et al. 2018; Bosveld et al. 2016; Pecreaux et al. 2016). Pulling is mediated at the cell cortex through interactions between astral microtubules and the motor protein dynein, which is anchored at the plasma membrane through its association with a tripartite complex consisting of NuMA, LGN and G\(\alpha \)i (Okumura et al. 2018; Bosveld et al. 2016; Pecreaux et al. 2016). Mathematical models have been used to investigate how astral microtubule pushing and pulling can drive spindle positioning. For example, minimising the calculated torque, created by pulling forces along the microtubule length (Minc et al. 2011) or by concentrated populations of dynein-associated proteins at the cell periphery (Théry et al. 2007; Bosveld et al. 2016), has been shown to predict the cell division orientation in sea urchin zygotes (Minc et al. 2011), micropattern-adhered HeLa cells (Théry et al. 2007) and the Drosophila pupal notum epithelium (Bosveld et al. 2016). These models highlight the importance of cell geometry (Minc et al. 2011) and the localisation of dynein (Théry et al. 2007; Bosveld et al. 2016) in spindle orientation and consequently cell division orientation. In a different approach, Li and Jiang (2017) used a stochastic model to describe the interactions of microtubules with chromosomes, motor proteins and boundaries to create self-assembled spindles within cells. This model was adapted to investigate spindle orientation: microtubules and chromosomes self-assemble into a mitotic spindle and orient within the simulated cells as a result of a combination of microtubule pushing forces and dynein-mediated pulling forces at both the cortex and within the cytoplasmic domain (Li et al. 2019), resulting in spindles which align with sites of localised dynein similarly to what has been shown in simpler models (Théry et al. 2007; Bosveld et al. 2016). Interestingly, the simulated spindles were shown to form already in line with their final division axis, with no notable movements of the spindle once assembled.
However, the mitotic spindle has been observed to approach its final destination less directly after its assembly. During the first division of the C. elegans fertilised egg, the posterior spindle pole undergoes a defined oscillation as the mitotic spindle is asymmetrically positioned in the cell to produce two daughter cells of unequal sizes (Pecreaux et al. 2006, 2016). Similarly, the rotational movements of the spindle in Xenopus epithelial tissue have been shown to be dynamic, culminating in oscillations of the spindle angle immediately prior to anaphase (Larson and Bement 2017). In the developing airway epithelium of mice, a subset of cells have been identified which continuously change their mitotic spindle angle throughout metaphase (Tang et al. 2018). Different strains of nematodes related to C. elegans show robust spindle oscillations, albeit with inter- and intra-specific variations in dynamical features(Valfort et al. 2018).
In the context of highly coordinated spindle movements, such as the oscillation of the C. elegans zygote posterior pole, contributions from both microtubule pushing and dynein-mediated pulling have been shown mathematically to produce the observed oscillatory dynamics (Grill et al. 2005; Pecreaux et al. 2006; reviewed in Beta and Kruse 2017). This mathematical model describes changes in the position of the mitotic spindle pole in 1D as a result of pulling by cortical force generators, combined with microtubule-based restoring forces (Grill et al. 2005). The relative simplicity of this model compared with that used by Li and Jiang (2017) lends itself more readily to an investigation of the primary parameters giving rise to dynamic movements and, crucially, replicates the oscillations seen in C. elegans. Wu et al. (2024) also describe the C. elegans zygote posterior pole oscillation, using a 2D model omitting pushing from microtubules. This omission is congruent with their analysis of subcellular fluid flows, which suggests that cortical pulling forces dominate to drive movement. The region of the cortex subtending the microtubule array as the spindle pole approaches and recedes from the cortex is highlighted as a key factor for creating oscillations in the pole position. As a result, pushing forces are unnecessary to create a reversal of the spindle pole velocity, as the direction of pulling by force generators is re-distributed across angles away from the cortex upon approach of the spindle pole (Wu et al. 2024). Strikingly, the resulting oscillations are nonlinear, in contrast with those demonstrated by Grill et al. (2005), though the significance of these nonlinear oscillations has yet to be explored. The correct balance of microtubule pushing and dynein-mediated pulling forces are likely the drivers for producing spindle movements and subsequently for determining the division orientation.
Other cells demonstrate spindle movements which are more complex. By mathematical amplification of pulling forces on the mitotic spindle from discrete cortical locations, rotational spindle dynamics have been simulated to match those observed in HeLa cells (Corrigan et al. 2015). Stochastic switching between active and non-active cortical cues simulates noisy rotation toward the long axis of the cell as defined by anisotropy in the placement of the cortical locations (Corrigan et al. 2015), highlighting both the importance of cortical cue elements in spindle orientation and the possibility of stochasticity in creating dynamic movements of the spindle. Stochastic processes can result in behaviours which are not captured by deterministic models due to processes such as stochastic resonance (Erban and Chapman 2020). Indeed, the addition of noise inherent to biological systems (Tsimring 2014) should not be discarded in considerations of dynamic behaviour.
In Xenopus embryo epithelial tissue, mitotic spindles have been shown to undergo both a net rotation towards the final division axis and a stereotypical oscillation prior to anaphase onset (Larson and Bement 2017) (see Online Resource 1). Figure 1b–d illustrates such oscillations in Xenopus animal cap epithelial tissue, obtained by tracking the movement of the metaphase plate, with which mitotic spindle movements are highly correlated. Oscillations are noisy with a nonlinear structure suggestive of relaxation oscillations (with rapid reversals of direction, Fig. 1d), a feature not reportedly observed in the C. elegans spindle oscillation (Pecreaux et al. 2006). The factors which affect the structure of oscillations in the mitotic spindle, specifically the nonlinear structure identified in Fig. 1d, have not yet been fully described, although spindle movements driven by cortical force generators in the absence of microtubule pushing forces appear to create nonlinear oscillations (Wu et al. 2024), suggesting that relaxation oscillations may arise in pull-dominated systems. Furthermore, spindles which do not oscillate are also present within the same tissue (Fig. 1e, f), in contrast to the defined and characteristic spindle behaviour of the C. elegans zygote (Pecreaux et al. 2006, 2016). It is unclear how more complex tissue environments, such as is found in the Xenopus epithelium, may affect the ability of mitotic spindles to oscillate, or the non-linearity of the oscillatory spindle movements, motivating the present study of spindle dynamics.
In this paper, we revisit the mathematical model presented by Grill et al. (2005), investigating factors which promote relaxation and noise-driven oscillations. The model is outlined in Sect. 2 and stochastic simulations are presented in Sect. 3. Representations of solutions using mean-field Fokker–Planck equations which account for noise in the random walking of force generators, but which incorporate a deterministic representation of binding kinetics, are given in Sect. 4; these are then reduced to a set of ordinary differential equations (ODEs) using systematic asymptotic analysis in Sect. 5 assuming slow binding kinetics. The ODEs turn out to represent a special case of the (less formally derived) ODEs proposed by Grill et al. (2005). A stability analysis gives predictions for the period of oscillation at the onset of neutral oscillations, as well as the position of the neutral oscillation boundary in parameter space, in agreement with Fokker–Planck predictions. Further asymptotic reduction of the ODE model in the limit of small pushing forces (Sect. 5.2) yields a single algebraic equation which describes the structure of nonlinear relaxation oscillations. We provide numerical evidence that smaller-amplitude irregular oscillations, characteristic of observations (Fig. 1), can be induced by noise associated with random binding kinetics through stochastic resonance. This is supported by analysis in Sect. 6 of a stochastic differential equation, derived from the ODE model, that seeks to estimate the amplitude and spectrum of the noise-induced oscillations.
2 The Model of Spindle Pole Dynamics
In the 1D model of spindle pole dynamics proposed by Grill et al. (2005), a spindle pole at position \({\bar{z}}\left( {\bar{t}}\right) \) at a time \({\bar{t}}\) moves along an axis \({\hat{\textbf{z}}}\) spanning a cell according to
The parameter \({\bar{\xi }}\) models viscous drag on the pole from the cytoplasm. The stiffness parameter \(k_{\textrm{MT}}\) represents a restoring force towards the cell midplane, arising from dynamic instability and bending of astral microtubules that emanate from the spindle pole and extend to the cell cortex (Grill et al. 2005; Pecreaux et al. 2006, 2016; Howard 2006; Rubinstein et al. 2009). The spindle is pulled towards either side of the cell under fluctuating forces \({\bar{F}}^\pm ({\bar{t}})\). Pulling arises from individual force generators which lie at the cell cortex and bind to astral microtubules. For simplicity, the model considers two opposing populations of force generators which sit in an ‘upper’ and ‘lower’ cortex, labelled ± hereafter. The force generators comprise a motor protein head connected to the cortex via an elastic linker of stiffness \(k_{\text {g}}\) (Fig. 2). The motor protein head can be considered to be dynein, which binds to and walks along the microtubules towards the spindle pole. The two populations of N force generators are assumed to exert pulling forces toward their respective cortex with a magnitude
where \({\bar{y}}^{\left( n\right) \pm }_{\text {b}}({\bar{t}})\) is the extension of the elastic linker of bound (subscript b) force generator n. In (2) it is assumed that the pulling force due to an individual linker is proportional to its length. Unbound (subscript u) force generators of length \({\bar{y}}^{\left( n\right) \pm }_{\text {u}}\) (Fig. 2c) are not connected to the spindle pole and are unable to provide any forcing. The superscript (n) in (2) is used to label the \(n_{\textrm{b}}^\pm \) linkers that, in any short interval \(({\bar{t}}, {\bar{t}}+\delta {\bar{t}})\), are bound to a microtubule, where \(0\le n_{\textrm{b}}^\pm \le N\).
The forcing in (1) fluctuates because the linkers bind and unbind randomly. The movements of the spindle pole are tightly coupled to the individual extension lengths of the linkers via (2), and by the fact that spindle motion influences the length of bound linkers. The motor protein heads have walking velocities given by
Here \(f_0\) is the stall force of a force generator, i.e. the force required to bring the motor protein head to rest relative to the spindle pole; it is assumed that the unloaded walking velocity \(v_0\) is reduced in proportion to \(k_{\text {g}}{\bar{y}}^{(n)\pm }_{\text {b}}/f_0\), the tensile force acting upon the motor protein head by the elastic linker scaled relative to \(f_0\). Equivalently, \(y_0\equiv f_0/k_g\) is the extension of a linker at which it stalls. In (3a), the spindle pole velocity term \(\mp {\text {d}{\bar{z}}}/{\text {d}{\bar{t}}}\) arises due to the force generator being connected to the moving spindle pole via the microtubules. Thus, as the spindle pole moves towards a bound force generator it will compress the elastic linker by pushing on the bound motor, reducing its relative walking velocity (Fig. 2). After a linker detaches from a microtubule, becoming unbound, it contracts with velocity
where \({\bar{\xi }}_g\) is a drag coefficient of unbound dynein. Using the largest proteins in the force generator complex (dynein, length approximately 50 nm (Trokter et al. 2012), and NuMA, length approximately 210 nm (Compton and Cleveland 1993)), a force generator has a Stokes radius of order \(10^{-1}\) smaller than the spindle pole, so that \({\bar{\xi }}_{\text {g}}\approx {\bar{\xi }}\times 10^{-1}\). The superscript (n) in (3b) is used to label the \(n_{\textrm{u}}^\pm \) linkers that, in the short interval \(({\bar{t}}, {\bar{t}}+\delta {\bar{t}})\), are unbound, where \(0\le n_{\textrm{u}}^\pm \le N\) and \(n_{\textrm{u}}^\pm +n_{\textrm{b}}^\pm =N\)
The model is closed by relating linker lengths to linker velocities, incorporating noise in the linker dynamics through effective diffusion coefficients \({\bar{D}}_{\textrm{u}}\) and \({\bar{D}}_{\textrm{b}}\), and by modelling the transitions between bound and unbound states as random events taking place at rates \({\bar{\omega }}_{\textrm{on}}\) and \({\bar{\omega }}_0 \exp [{\bar{\gamma }}{\bar{y}}_{\textrm{b}}^{(n)\pm }]\) respectively. Here \({\bar{\gamma }}\) parametrizes the slip-like manner in which dynein detaches from microtubules under loading (Ezber et al. 2020). Estimated values of dimensional parameters are summarized in Table 1.
Scaling lengths on the stall length \(y_0\) and time on the stall time \(y_0/v_0\) (so that \({\bar{z}}=y_0 {z}\), \({\bar{t}}=(y_0/v_0){t}\), etc.), the force balance on the spindle pole (1, 2) becomes, in dimensionless form,
\({\xi }={\bar{\xi }}v_0/(k_{\text {g}}y_0)\) and \(K=k_{\text {MT}}/k_{\text {g}}\) are dimensionless drag and stiffness parameters respectively. The velocities of the bound and unbound generators (3a) become
where \(\Gamma =f_0/({\bar{\xi }}_g v_0)\). This parameter measures \(k_g/{\bar{\xi }}_g\), the retraction rate of unbound linkers, relative to \(v_0/y_0\), the stall rate. Dimensionless counterparts of the stochastic parameters are diffusion coefficients \({D_{\textrm{b}}}\) and \({D_{\textrm{u}}}\) describing the mobility of bound and unbound linkers respectively, and transition rates \({\omega }_{\textrm{on}}\) and \({\omega }_0 e^{\gamma {y}_{\textrm{b}}^{(n)\pm }}\) respectively. Dimensionless parameters are summarised in Table 2.
Over any short interval, the populations of bound and unbound linkers have a distribution of lengths. Average extensions are defined by
3 Stochastic Simulations
To capture the discrete interactions between a small number of force generators and the spindle pole, we discretize \({y}^{\left( n\right) \pm }_{\text {b}}\) in increments of \(\Delta {y}\) and use a Gillespie algorithm to model the stochastic extensions and retractions of bound and unbound linkers and the stochastic transitions of the binding state of the force generators, as they bind and unbind from microtubules. As explained in Appendix B, the extension and retraction of the force generators are treated as 2N biased random walks with drift \({v}_{\text {b(u)}}^\pm \), diffusion \({D_{\mathrm {b(u)}}}\) and state change (between bound and unbound states), all coupled to displacement of the spindle.
Figure 3 presents a simulation displaying the emergence of spontaneous oscillations of the spindle pole, using the parameters shown in Table 1, with \(N=15\) linkers at either cortex. The spindle location (Fig. 3a), numbers of bound and unbound force generators \(n^\pm _{\text {b(u)}}\) (Fig. 3b) and average extensions \(\langle {y}^\pm _{\text {b(u)}}\rangle \) (6) (Fig. 3b, c) show noisy but oscillatory dynamics. The average extensions of the bound force generators \(\langle {y}^+_{\text {b}}\rangle \) and \(\langle {y}^-_{\text {b}}\rangle \) (6) oscillate in anti-phase to one another (Fig. 3b). The average extension of unbound force generators \(\langle {y}^\pm _{\text {u}}\rangle \) remains close to 0 following initial transients (Fig. 3b). This can be explained by considering the movement of the spindle pole through one cycle of oscillation (Fig. 3a) and \(\langle {y}^\pm _{\text {b}}\rangle \left( {z}\right) \) (Fig. 3c), discussed further below. Apparent gaps in the \(\langle {y}^\pm _{\text {b}}\rangle \) plots (Fig. 3c) occur where there are no bound generators from which to extract an average (where \(n_{\text {b}}\) = 0 in Fig. 3b).
Consider the following phases of movement identified by coloured symbols in Fig. 3a, c.
-
1.
Spindle moving away from the upper cortex (green to cyan) At the peak of the spindle pole oscillation, movement of the spindle is dominated by the microtubule restoring force. The bound generators are extended equally in the upper and lower cortices (\(\langle {y}^+_{\text {b}}\rangle \sim \langle {y}^-_{\text {b}}\rangle \) at the green timepoint (Fig. 3c)) though there are a greater number bound in the upper cortex rather than the lower (\(n_{\text {b}}^+ > n_{\text {b}}^-\), comparison in Fig. 3(bi) vs (bii)). The restoring force (\(-K{z}\)) is greater than the net upward pulling force provided by this unbalanced population ratio. As the spindle pole moves towards \({z}=0\), this restoring force decreases while the increasing spindle pole velocity results in a net compression of the elastic linkers on the lower cortex, due to a switch in the sign of \({v}^-_{\text {b}}\left( \langle {y}^-_{\text {b}}\rangle \right) \). Additionally, the spindle pole velocity increases the relative velocity of the force generators in the upper cortex, resulting in an extension of the elastic linkers at the upper cortex (Fig. 3c), and shortening of the linkers at the lower cortex. Due to the tension-sensitive unbinding rate \({\omega }_0e^{\gamma {y}^+_{\text {b}}}\), this results in a gradual decrease in the number of upper bound force generators as \(\langle {\omega }_0e^{\gamma {y}^+_{\text {b}}}\rangle \) increases in value, while the number of bound force generators in the lower cortex increases due to a constant binding rate and a decreased unbinding rate (Fig. 3b).
-
2.
Spindle moving through the centre of its oscillating range, toward the lower cortex (cyan to yellow) As the spindle moves through \({z}=0\) the restoring force steadily increases from 0 to \(-K{z}\). This slows the movement of the spindle such that the velocity of the force generators in the lower cortex may become positive \({v}^-_{\text {b}}\left( \langle {y}^-_{\text {b}}\rangle \right) >0\) which allows these elastic linkers to extend (Fig. 3c), decreasing the relative velocity of the remaining upper force generators, the average extension of which is also reduced due to the unbinding of those with larger extensions and binding of force generators with reduced extensions (Fig. 3c). The number of bound generators in the lower cortex also begins to decline as they extend due to the increased unbinding rate (Fig. 3(bii)).
-
3.
Spindle moving away from the lower cortex (yellow to magenta) This phase replicates the first phase, but with the behaviours of upper and lower cortex reversed. The motion away from the cortex due to the restoring force results in a compression of the upper elastic linkers and an extension of the lower elastic linkers (Fig. 3c), and a corresponding decrease in the absolute number of bound force generators in the lower cortex as opposed to the increased binding observed in the upper cortex (Fig. 3b).
The closed loops in \((\langle {y}^\pm _{\text {b}}\rangle , {z})\) space are traced anti-clockwise in the lower cortex and clockwise in the upper cortex (Fig. 3c). At the stall force (when \(\langle {y}^\pm _{\text {b}}\rangle \approx 1\)), the direction of the solution loop is determined by the direction of acceleration of the spindle pole with respect to the cortex. That is, a force generator in the lower cortex whose elastic linker is at \({y}^{(n)\pm }_{\text {b}}=1\) will be decreasing its extension as the spindle pole accelerates toward it (negative acceleration, green point in Fig. 3a, c) and increasing as the spindle pole accelerates away (positive acceleration, yellow point in Fig. 3a, c).
Removing the tension sensitivity of unbinding by setting \(\gamma =0\) results in less well-defined oscillations (Fig. 3d) of reduced amplitude relative the baseline case shown in Fig. 3a. Thus the tension-sensitive unbinding rate appears to promote coherent oscillations of the spindle pole, although fluctuations persist due to the stochastic binding and unbinding. Similarly, increasing the restoring force by increasing the parameter K reduces the deviation in the position of the spindle pole from the centre (Fig. 3e), but also leads to a marked reduction in the coherence of the spindle motion (compare Fig. 3a and e).
4 A Fokker–Planck Description
Simulating the system stochastically reveals the role of noise in individual realisations of spindle dynamics. To explore properties of the model over multiple realisations in a more computationally efficient manner, we turn to a system of partial differential equations (PDEs) for probability density functions (pdfs) \(P_{\textrm{b}(\textrm{u})}({y},{t})\) for the extensions of bound and unbound linkers at either cortex, where the elastic linker extension y is now considered as a continuous variable. The model may be written as
where
Equation (7) is a nondimensional version of the mean-field Fokker–Planck equations proposed by Grill et al. (2005). The continuous velocities \({v}^\pm _{\text {b(u)}}({y})\) in (8) evolve as in (5). The pulling force toward each cortex (2) is calculated as
modifying the force balance on the spindle (4), which becomes
The boundary conditions
ensure conservation of total probability
Given some initial conditions \({P}_{\text {b}}^\pm \left( {y},0\right) ={P}_{\text {b}0}^\pm \left( {y}\right) \), \({P}_{\text {u}}^\pm \left( {y},0\right) ={P}_{\text {u}0}^\pm \left( {y}\right) \), and \({z}\left( 0\right) ={z}_{0}\), the system (7–11) may be solved in time to return the dynamics of the spindle pole and the populations of cortical force generators, represented as probability densities over multiple realisations of the system. We computed numerical solutions using the method of lines.
The underlying stochastic system (Appendix B) combines two sets of random processes: binding and unbinding of force generators; and the random walk of force generators along microtubules. The description of the system provided by (7) can be described as mean-field in the sense that it combines a deterministic model of binding/unbinding kinetics (via the reaction terms in (7)) with a stochastic description of force-generator motion (via the diffusive terms in (7)). We shall revisit the role of randomness in binding kinetics in Sect. 6 below.
The solutions of (7)–(12) presented in Fig. 4 show an oscillating spindle displacement z(t) corresponding to fluctuations in \({P}^{\pm }_{\text {b}}\left( y,{t}\right) \) and \({P}^{\pm }_{\text {u}}\left( y,{t}\right) \). For large \(\Gamma \), (7b) is dominated by the advective term which sweeps any unbound force generators with a non-zero extension down toward \({y}=0\). As there is no flux through this boundary by (11), \({P}_{\text {u}}^\pm \) has a defined peak at \({y}=0\) which decays with y over the diffusive lengthscale \(D_{\textrm{u}}^{1/2}\). For the bound pdfs \({P}_{\text {b}}^\pm \), the location \({y}_{\text {c}}^\pm \) and amplitude \({P}_{\text {b}}^{\pm ,\text {max}}\) of the maximum of the pdf oscillate concurrently with z (Fig. 4c, g), mirroring the behaviour of the average extension \(\langle {y}^\pm _{\text {b}}\rangle \) and number of bound force generators \(n^\pm _{\text {b}}\) in the stochastic simulation (Fig. 3b). Variations of the initial conditions \({P}^{\pm }_{\text {b0}}\) and \({P}^{\pm }_{\text {u0}}\) had no effect on the final solutions following initial transients (data not shown).
Decreasing \({D_{\textrm{b}}}\) and \({D_{\textrm{u}}}\) by a factor of 10 results in taller and narrower pdfs (Fig. 4d, h), confining \({P}_{\text {b}}^\pm \) to a region of y which is spatially separated from \({P}_{\text {u}}^\pm \) at all times. In this limit we can partition the y-domain into three distinct regions: I, of width \({\mathcal {O}}(D_{\textrm{u}}^{1/2})\), encompassing the peak of \({P}_{\text {u}}^+\); III, of width \({\mathcal {O}}(D_{\textrm{b}}^{1/2})\), encompassing the peak of \({P}_{\text {b}}^+\); and II between them, which remains distinct throughout an entire oscillation (Fig. 4h, see Online Resource 2c). As well as modulating the shape of the pdfs, \({D_{\textrm{b}}}\) and \({D_{\textrm{u}}}\) also affect the resulting dynamics of the spindle pole. Decreasing \({D_{\textrm{b}}}\) and \({D_{\textrm{u}}}\) results in an increased period, T, of oscillation (\(T\approx 890\) increases to \(T\approx 1000\) upon a decrease in \({D_{\textrm{b}}}\) and \({D_{\textrm{u}}}\) by a factor of 10), a decrease in the amplitude of the oscillation (Fig. 4a, e) and longer transients (Fig. 4a, e).
The solution in Fig. 4a–d was run with parameters matching those in the stochastic simulation in Fig. 3, except that \(N=25\) in the former and \(N=15\) in the latter. Nevertheless, PDE predictions show a comparable period, without capturing the detailed fluctuations in an individual realisation. To obtain a broader view of parameter dependence, solutions of (7–11) for a range of values of N and \({\omega }_{\text {on}}\) are reported in Fig. 5a. For the baseline value \(\omega _{\textrm{on}}=0.003\), sustained oscillations arise in the PDE model with \(N=25\) and low diffusivities (as in Fig. 4e–h) but not \(N=15\) (see Fig. 5c). Increasing diffusivities with \(N=15\) leads to the sustained oscillations seen in the PDE model in Fig. 5d. Noise, therefore, is likely to play a role in promoting oscillatory dynamics.
The PDE stability boundary for low diffusivities is mapped out in (N, \({\omega }_{\text {on}}\))-space (Fig. 5a), distinguishing oscillatory from non-oscillatory solutions. The period of oscillation for neutrally-stable disturbances increases with \(\omega _{\text {on}}\) (Fig. 5b). Reduction of the number of force generators, leading to a decrease in pulling forces, results in a cessation of oscillations (Fig. 5a, c). For large N, two thresholds exist for values of \({\omega }_{\text {on}}\) at which oscillations arise, with the oscillatory section of parameter space forming a wedge shape (Fig. 5a). We explore the origins of these thresholds in more detail below. This wedge-shaped parameter space was described previously by Grill et al. (2005) through analysis of a reduced model where it was assumed that unbound force generators instantaneously relax down to zero extension. The presence of a threshold between oscillatory and non-oscillatory solutions has been experimentally validated in C. elegans embryos (Pecreaux et al. 2006).
Significantly, the oscillatory behaviour reported in stochastic simulations for \(N=15\) (Fig. 3a), lies in a regime in which the Fokker–Planck model predicts steady distributions of \(P_{b(u)}^{\pm }\) with \(z\rightarrow 0\) at large times (Fig. 5c). We provide evidence below that the sustained oscillations in Fig. 3a are noise-driven (or a form of stochastic resonance (Erban and Chapman 2020)), with noise arising from stochastic binding kinetics, a feature missing from the mean-field Fokker–Planck model.
Additional PDE solutions with low diffusivities are reported in Fig. 6a, c, illustrating respectively the impact of increasing N and decreasing the strength of the restoring force K. The latter leads to larger-amplitude oscillations with a relaxation structure, characterised by periods of approximately uniform spindle velocity, interspersed with rapid changes in direction (Fig. 6c). Correspondingly, the oscillations combine slow phases in the time-evolution of \({P}^\pm _{\text {b}}\) and \({y}^\pm _{\text {c}}\) (Fig. 6c) in which z is approximately linear in t, with short intervals in which the rapid change in the direction of motion of the spindle pole coincides with fast changes in the extension of the force generator bound probability centre \({y}^\pm _{\text {c}}\) and amplitude \({P}^\pm _{\text {b}}\). We explore the origins of this strongly nonlinear behaviour below, through comparison to a simplified model reported in the remaining panels of Fig. 6.
The Fokker–Planck description (7–12) reveals many of the characteristics promoting spindle-pole oscillations, but still requires extensive computation. We now reduce this model to a system of ODEs by asymptotic analysis, allowing us to more fully explore the relationships between the most important factors promoting oscillations. Rather than follow the heuristic approach proposed by Grill et al. (2005), we seek a systematic reduction valid at an appropriate distinguished limit in parameter space.
5 Asymptotic Reduction to ODEs
When diffusivities are small, the PDEs reveal distinct regions of y space where the pdfs \({P}^\pm _{\text {u}}\) and \({P}^\pm _{\text {b}}\) have most of their mass (Fig. 4h). While varying the diffusive terms has an impact on the oscillations, with larger diffusive terms promoting oscillations (Fig. 5d), the amplitudes and periods of the more and less diffusive solutions are still of a similar order. We now pursue the behaviour of the model with lower diffusivity to create a system of ODEs.
We develop an approximation to the oscillating spindle system in a distinguished limit for which \({\omega }_{\text {on}}\sim {\omega }_0\sim {D_{\textrm{b}}}^{1/2}\sim {D_{\textrm{u}}}^{1/2}\ll 1\) (where \(\sim \) means ‘scales like’) and rescale the time and spindle position parameters in the PDE problem (7), (8) and (10) by \({t}={\tilde{t}}/{\omega }_{\text {on}}\) and \({z}={\tilde{z}}/{\omega }_{\text {on}}\). In other words, we assume that binding kinetics happens slowly relative to movement of force generators, and that fluctuations of such movement are weak. The range of extension values y are split into the three regions identified in Fig. 4h: I, over which \({P}_{\text {u}}^\pm \) is peaked around \({y}=0\) with a width \({D_{\textrm{u}}}^{1/2}\); III, over which \({P}_{\text {b}}^\pm \) is peaked with a width of \({D_{\textrm{b}}}^{1/2}\) but whose centre moves as \({y}_{\text {c}}=1\mp {\tilde{z}}_{{\tilde{t}}}\); and II, where advective terms dominate.
In Appendix C, we express the governing equations in rescaled coordinates, expand in powers of \(\omega _{\textrm{on}}\), solve for \({P}_{\text {u}}^\pm \) in region I and \({P}_{\text {b}}^\pm \) in region III, and then match asymptotic limits across region II. This procedure yields the three coupled ODEs
Here \({\hat{B}}^\pm \left( {\tilde{t}}\right) =\sqrt{2\pi {D_{\textrm{b}}}}B^\pm \), where \(B^\pm \) approximates the amplitude of the peak of \({P}^\pm _{\text {b}}\). The parameters in (13) are
Formally, the parameters (14) are assumed to remain \({\mathcal {O}}(1)\) in the limit \({\omega }_{\text {on}}\sim {\omega }_0\sim {D_{\textrm{b}}}^{1/2}\sim {D_{\textrm{u}}}^{1/2}\ll 1\). \(\rho \) is the binding affinity (under zero load) of force generators for microtubules; \({\hat{\xi }}\) measures the spindle drag (assuming the spindle moves at the walking speed of a linker) relative to the stall force generated by the full population of linkers; \({\hat{K}}\) measures the restoring force, driving the spindle to the centre of the cell (assuming a displacement of the spindle comparable to the distance walked by a linker), relative to the stall forces generated by the full population of linkers. We solved (13) numerically, imposing initial conditions \({\tilde{z}}_0\) and \({\hat{B}}^\pm _0\).
The ODEs defined by Grill et al. (2005) may be re-written in the notation used above as
where \({\tilde{y}}^\pm \), the typical length a linker extends before it detaches, is determined from
Assuming \(\omega _0\ll 1\) in (15c), then \({\tilde{y}}\simeq 1\mp {\tilde{z}}_{{\tilde{t}}}\) and we recover (13). Our asymptotic reduction therefore recovers a special case of the heuristically-determined ODEs presented by Grill et al. (2005).
In Fig. 6a–d, two solutions of the Fokker–Planck system are compared with solutions of the ODEs (13) for equivalent parameters, taking \({y}_{\text {c}}^\pm =1\mp {\tilde{z}}_{{\tilde{t}}}\). The spindle pole dynamics and associated force generator behaviours in both cortices are fully captured in the ODE model in both cases.
Relaxation oscillations arise when pulling forces dominate over restoring forces, through a reduction of \({\hat{K}}\). Figure 6c, d shows a strong match between the relaxation oscillations returned by the PDE and ODE models for an equivalent reduction in K, with similar period and amplitude as well as shape. Likewise applying an increase in N (Fig. 6e), the ODE model predicts a strong relaxation structure (with very sharp changes in \({y}_{\text {c}}^\pm \) at the peaks of oscillation). The amplitude of this relaxation structure is, however, larger than would be viable within a biological cell. This model description omits explicit boundaries, and we would expect that within a cell the presence of a boundary would introduce a nonlinear contribution to the restoring force close to the cell edge. Despite this, we conclude that the balance of pulling to restoring forces controls the general structure of the oscillation of the spindle pole.
Interestingly, the oscillations in the lower-\({\hat{K}}\) and high-N relaxation oscillation are slightly different. For low \({\hat{K}}\), the peak of the bound pdf hits its maximum at the same time as the spindle pole experiences its maximum velocity (when \({y}_{\text {c}}^\pm \) is at its minimum value, Fig. 6c, d). Alternatively, when N is increased (leading to changes in \({\hat{K}}\) and \({\hat{\xi }}\)), the maximum of the peak of the bound pdf lags behind the spindle pole velocity (Fig. 6e). This lag represents a delay between the binding of the force generators and the movements of the spindle pole, likely due to there being a greater number of force generators in the system to bind to the microtubules before saturation of the force generators onto the microtubule. Indeed when \({\omega }_{\text {on}}\) is small, the oscillations of the spindle pole are more non-linear (data not shown), as the number of bound force generators takes longer to saturate.
5.1 Stability Analysis
The simplicity of the ODE model (13) lends itself to stability analysis. Linearising about the steady state
assuming that small disturbances are proportional to \(e^{s{\tilde{t}}}\), yields the relationship
Setting \(s=\mu +i\Omega \) and collecting real and imaginary parts defines the growth rate
and frequency of oscillation
Setting \(\mu =0\) at the onset of neutral oscillations identifies the frequency
and the stability threshold
Both (19b) and (19a) provide good predictions of the stability boundary identified by PDE solutions (Fig. 5a) and the period of oscillations at the stability boundary (Fig. 5b).
The frequency of oscillation determined by Grill et al. (2005) via (15) may be rewritten using the notation above as
at the stability threshold
Equations (20a) and (20b) are equivalent to (19a) and (19b) respectively when \(\omega _0\sim \omega _{\text {on}}\ll 1\). (19b) and (20b) are compared in Fig. 5a, showing a near-perfect match despite the additional terms present in (20b); both predictions bound almost perfectly the oscillatory region observed by individual solutions of the PDEs.
We highlight two limits of (19b). First, taking \({\hat{K}}\rightarrow 0\), leaving \({\hat{\xi }}={\xi }/N\) as the only parameter which depends on N, (19b) reduces to
provided the denominator \(2[\gamma (\lambda -1)-\lambda ]\) is positive, i.e.
In this limit, the period of oscillation \(T=2\pi /\Omega \) is
\(\omega _{\textrm{on}}\rightarrow \omega _{\textrm{on}}^\dag \) also appears in the large-N limit of (19b), taking \({\hat{K}}\sim {\hat{\xi }}\ll 1\). In addition, taking \(N\gg 1\), \(\omega _{\textrm{on}}\ll 1\) with \(N \omega _{\textrm{on}}=O(1)\), we recover the additional approximation (provided \(\gamma >1\))
for which
Thus, for large N, the upper branch of the stability boundary defined by (19b) in Fig. 5a approaches \(\omega _{\textrm{on}}=\omega _{\textrm{on}}^\dag \) in (22), confirming that a necessary condition for oscillations is that the tension-sensitivity parameter satisfies \(\gamma >1\), i.e. that linkers exhibit slip-bond behaviour. Indeed, removal of the tension-sensitivity of the unbinding rate in the stochastic simulations leads to a reduction of the coherence of the oscillatory behaviour of the spindle pole (Fig. 3d). The upper-branch asymptote \(\omega _{\textrm{on}}=\omega _{\textrm{on}}^\dag \) appears to be shared also by PDE solutions (which suggests an upper stability threshold between \(0.006<{\omega }_{\text {on}}^\dag <0.007\) for \(N\le 80\), within \(80\%\) of \({\omega }_{\text {on}}^\dag = 0.0074\)). Also in the large-N limit, the lower branch of (19b) is captured by (24), consistent with PDE solutions in this limit. This limit shows explicitly how increasing the restoring force K has a stabilising effect.
We also recall that, in the Fokker–Planck model, decreasing the restoring force parameter K promotes oscillations at smaller N (Fig. 6c, where \(N=15\)). This behaviour is conserved in the ODE system, where the low-\({\hat{K}}\) approximation (21) shown in Fig. 5a, predicts oscillations in a greater region of the (\(N, {\omega }_{\text {on}}\))-plane. Evaluating \(\textrm{d}\omega _{\textrm{on}}/\textrm{d}N\) using (21) gives
Thus for the neutral curve to lie in \(N>0\) requires
providing a lower bound on the number of linkers needed for oscillations in terms of the walking speed and stall force of a linker, and the drag on the spindle.
The period of oscillations along the neutral stability curve predicted using (19a) increases as K decreases (Fig. 5b); thus a reduction of restoring forces corresponds to longer periods of oscillation. The rapid increase of the period as \({\omega }_{\text {on}}\rightarrow {\omega }_{\text {on}}^\dag \) coincides with \(N\rightarrow \infty \). (19a) is well matched with (20b) determined by Grill et al. (2005), as well as with the periods along the approximate stability curve identified by numerical solutions of the PDEs.
5.2 The Structure of Relaxation Oscillations
A further simplification to the model can be implemented by exploiting \({\hat{K}}\) as a small parameter. The approximately linear sections of \({\tilde{z}}\) (e.g. Fig. 4c) scale like \({{\hat{K}}}^{-1}\) in both time and amplitude, and are interrupted by rapid changes in spindle direction. Re-scaling \({\tilde{t}}=\tilde{{\tilde{t}}}/{\hat{K}}\) and \({\tilde{z}}=\tilde{{\tilde{z}}}/{\hat{K}}\) such that \({\tilde{z}}_{{\tilde{t}}}=\tilde{{\tilde{z}}}_{\tilde{{\tilde{t}}}}\), the ODEs (13) describing the slower phases of the dynamics become
Posing expansions \({\hat{B}}^\pm = {\hat{B}}_0^\pm + {\hat{K}}{\hat{B}}_1^\pm +\cdots \) and \(\tilde{{\tilde{z}}} = \tilde{{\tilde{z}}}_0+{\hat{K}}\tilde{{\tilde{z}}}_1 +\cdots \nonumber \), to leading order (28) becomes
We may rewrite (29b) as
with \({\hat{B}}_0^{\pm }\) defined by (29a).
The displacement-velocity reationship (30) approximates the slow phases of the limit cycle in (\(\tilde{{\tilde{z}}}_0\), \(\tilde{{\tilde{z}}}_{0,\tilde{{\tilde{t}}}}\))-space as \({\hat{K}}\rightarrow 0\). Recalling that \({y}^\pm _{\text {c}}=1\mp {z}_{{t}}\), then following a parameter rescaling, (30) can also be used to describe the limit cycle in (\({z}_0\), \({y}^\pm _{\text {c}}\)) (black curve in Fig. 6c–e). The limit cycles obtained by solving the ODE and PDE systems with equivalent parameters are shown to closely match with this expected limit cycle (Fig. 6c–e). These cycles show the fast phases of the relaxation oscillation as the spindle pole changes its direction of motion (the approximately vertical sections at the maximum and minimum values of \(\tilde{{\tilde{z}}}\)). The maximum amplitude of oscillation can be estimated by the roots of G, which can be determined by solving
for roots \(G_{\text {max}}\) and \(G_{\text {min}}\). Then the amplitude of oscillation during relaxation oscillations can be estimated by
Thus the amplitude of oscillation can be estimated from the ratio of pulling to pushing (\({\hat{K}}\)), the effective drag (\({\hat{\xi }}\)), the ratio of the unbinding to binding rates (\(\rho \)) and the tension sensitivity of unbinding (\(\gamma \)).
This approximation also illustrates how the tension-sensitivity of the cortical force generators, mediated by \(\gamma \), is key for oscillations. Setting \(\gamma =0\) in (29a) uncouples the values of \({\hat{B}}^\pm _0\) from the spindle pole dynamics, thus \({\hat{B}}^+_0={\hat{B}}^-_0\) and (29b) becomes \( \tilde{{\tilde{z}}}_0 = -\left( {\hat{\xi }}+2\left( 1+\rho \right) ^{-1}\right) \tilde{{\tilde{z}}}_{0,\tilde{{\tilde{t}}}}\), giving a linear relationship between \(\tilde{{\tilde{z}}}_0\) and \(\tilde{{\tilde{z}}}_{0,\tilde{{\tilde{t}}}}\) and eradicating the limit cycle. Thus the coupling of the populations of bound force generators through the tension-sensitive unbinding rate is required for oscillations, as was shown by stability analysis of the ODE system (22).
5.3 Testing the Accuracy of the ODE System
For \(\omega _0=0.001\), comfortably satisfying the condition \(\omega _0\ll 1\) that allows the reduction of the Fokker–Planck system (7–12) to the ODEs (13), the latter make accurate predictions for the onset of sustained oscillations (Fig. 5a). In this limit, the stability threshold (19b) is almost indistinguishable from that of the heuristic model (20b) proposed by Grill et al. (2005). To test the robustness of each approximation, Fig. 7 shows predicted stability thresholds for \(\omega _0=0.1\), against solutions of (7–12). A marked difference in the minimum value of N required to cross the neutral curve is observed, with the threshold (19b) underestimating the lowest value of N required to elicit oscillations in the PDE solutions (by 45% for \(\omega _{\textrm{on}}=0.3\)). Alternatively, (20b) overestimates the threshold value of N (by 77% for \(\omega _{\textrm{on}}=0.3\)), demonstrating a modest benefit of the rigorously derived ODE system in this parameter regime.
6 Characterising Noise-Induced Oscillations
We now use the ODE system (13) to provide further evidence that the oscillations in Fig. 3a are noise-induced. Despite lying outside the neutral curve (Fig. 5a), the period of the stochastic oscillations is well approximated by (18b) (Fig. 3a, red bar), indicating that noise due to the binding and unbinding of a relatively small number of linkers may be sufficient to overcome the damping evident in the Fokker–Planck description (Fig. 5c) and in the ODE model. As explained in Appendix B, the Fokker–Planck system (7, 10) proposed by Grill et al. (2005) is a simplified form of the high-dimensional chemical Fokker–Planck equation associated with the full stochastic model; we attribute the failure of (7, 10) to predict the oscillations in Fig. 3a to this simplification, and show below how reintroducing stochastic effects associated with binding kinetics can explain some features of observations.
To do so, we adopt the framework outlined by Boland et al. (2008) to estimate the amplitude of small-amplitude noise-induced oscillations. At small amplitudes, it is appropriate to linearise the exponential term in (13) using \(\textrm{exp}(\pm \gamma {\tilde{z}}_{{\tilde{t}}})\approx 1\pm \gamma {\tilde{z}}_{{\tilde{t}}}\). We then treat the ODE model (13) as a chemical kinetic system with eight reactions, written as
where \(a\equiv \gamma ({\hat{B}}^+-{\hat{B}}^--{\hat{K}}{\tilde{z}})/({\hat{\xi }}+{\hat{B}}^++{\hat{B}}^-)\). The columns \(\varvec{\nu }_i\) of the \(3\times 8\) stoichiometric matrix can be assembled with the reaction rates \(a_i\) (\(i=1,\dots ,8\)) to form the correlation matrix \({\textsf{D}}=\tfrac{1}{2}\sum _i \varvec{\nu }_i\varvec{\nu }_i^\top a_i\), where
and \(\lambda \equiv 1+\rho e^\gamma \). Evaluated at the equilibrium point (16), \({\textsf{D}}\) simplifies to
Linearising (33) about the equilibrium yields the Jacobian matrix \({\textsf{J}}^*\) satisfying
The eigenvalues of \({\textsf{J}}^*\) satisfy (17). For a particular set of parameters, including \({\hat{K}}_c\), \({\hat{\xi }}_c\) satisfying (19b), at which \({\textsf{J}}^*={\textsf{J}}^*_c\) (say), the Jacobian has one real negative eigenvalue and a complex conjugate pair with zero real part and frequency \(\Omega _c\) satisfying (19a). Moving away from neutral stability by changing N (i.e. moving horizontally in Fig. 5a) can be represented by setting \({\hat{K}}={\hat{K}}_c(1-\epsilon )\), \({\hat{\xi }}={\hat{\xi }}_c(1-\epsilon )\) for some \( \epsilon \equiv \delta N/N\ll 1\) (using (14)), with \(\lambda \) and \(\gamma \) remaining fixed. Perturbing (17) in this way leads to complex eigenvalues
confirming instability (\(\textrm{Re}(s)>0\)) for an increase in N (\(\delta N>0\)). For \(\epsilon <0\) (associated with a small reduction in N from the neutrally stable case), we propose that the small negative growth rate (37) balances the noise from stochastic forcing to determine the amplitude of noisy oscillations.
Writing \({\textbf{x}}({\tilde{t}})=({\hat{B}}^+,{\hat{B}}^-,{\tilde{z}})^\top \), the stochastic differential equation that generalises (33) to describe small-amplitude noise-driven oscillations can then be written
Here \({\textbf{f}}=\sum _i \varvec{\nu }_i \sqrt{a_i}\) and \(W({\tilde{t}})\) is a Wiener process. Following Gardiner (1985), when \(\delta N<0\) (so that the eigenvalues of \({\textsf{J}}^*\) have negative real part), the stationary covariance \(\varvec{\sigma }\equiv \langle {\textbf{x}}(t),{\textbf{x}}^\top (t)\rangle \) satisfies the Lyapunov equation \({\textsf{J}}^*\varvec{\sigma }+\varvec{\sigma }{\textsf{J}}^{*\top }=-2{\textsf{D}}^*\). (Equivalently, the stationary distribution of the Fokker–Planck equation associated with (38) is proportional to \(\textrm{exp}[-\tfrac{1}{2} {\textbf{x}}^\top \varvec{\sigma }^{-1}{\textbf{x}}]\).) Writing
the coefficients satisfy
where \(\alpha \equiv (\lambda -1)\gamma -\lambda (2+\lambda {\hat{\xi }})\). We can use \(\sqrt{f}\) to estimate the amplitude of noise-driven oscillations in Fig. 3a lying outside the neutral curve. Solving (40) for f gives
Perturbing about the neutral curve, f simplifies to
in terms of the eigenvalues (37). Clearly this estimate of f is unbounded as \(\vert \Re (s)\vert \rightarrow 0\), violating the small-amplitude assumption and suggesting that (42) is best considered as an approximate upper-bound of the true amplitude.
Using the parameters for the stochastic simulation in Fig. 3a (a further example is shown in Fig. 8c), \(\sqrt{f}/\omega _{\text {on}}\approx 510\). Approaching the neutral curve by increasing N from 15 to 18 (Fig. 8d) increases \(\sqrt{f}/\omega _{\text {on}}\) by 35% to approximately 690. Figure 8e, f shows that the interquartile range of the density of z values of stochastic oscillations increase by 38% (from 59.5 when \(N=15\) to 82.1 when \(N=18\)). \(\sqrt{f}\) therefore captures the trend in amplitude but overestimates its magnitude by approximately a factor of 8. Discrepanices may arise from a number of sources, including linearisation of the exponential term in (13), linearisation leading to (42), interaction with random motion of force generators and insufficient sampling of stochastic time series.
The spectrum matrix \({\textsf{S}}\) (the Fourier transform of the time correlation matrix in the stationary state) satisfies (Gardiner 1985)
The \({\tilde{z}}{\tilde{z}}\) component of \({\textsf{S}}\) is plotted in Fig. 8a using the parameters for the stochastic simulation in Fig. 3a (\(N=15\)), and for \(N=18\) (approaching the neutral curve). Equation (43) predicts a sharpening of the spectrum with an approximate doubling of the root-mean-square amplitude near the resonant frequency; the maximum of \(\sqrt{{\textsf{S}}_{3,3}}/\omega _{\text {on}}\) corresponds to \(z\approx 360\) and \(z\approx 700\) for \(N=15\) and \(N=18\) respectively, broadly consistent with values of \(\sqrt{f}/\omega _{\text {on}}\). The relation \(\varvec{\sigma }=\int _{-\infty }^\infty {\textsf{S}}(\omega )\,\textrm{d}\omega \) emphasises the contributions of a narrower range of frequencies to the noise-induced oscillation as N increases. Accordingly, Fig. 8c, d shows a more coherent oscillation for larger N. Figure 8b illustrates similar increases in the predicted amplitude of stochastic oscillation as the neutral curve in Fig. 5a is approached along different values of \(\omega _{\textrm{on}}\); again this is best interpreted as a likely upper bound.
7 Discussion
We have investigated the factors promoting relaxation and noise-driven oscillations of the mitotic spindle identified experimentally (Fig. 1), by revisiting the mathematical model proposed by Grill et al. (2005). To this end, we used stochastic simulations to demonstrate the effect of noise on 1D pole movement (Sect. 3, Fig. 3); this involved 2N random walkers (linkers) switching between bound and unbound states, with their motion coupled via an ODE to that of the spindle. The corresponding mean-field Fokker–Planck equations (Sect. 4, Fig. 4), involving four PDEs coupled to an ODE, were reduced systematically (Sect. 5) to three nonlinear ODEs (13). When binding kinetics is slow (\(\omega _0\ll 1\)), these turn out to be a special (and simpler) case of the ODE system presented by Grill et al. (2005), and both systems show close agreement with the Fokker–Planck solutions in predicting conditions necessary for the onset of self-excited oscillations (Fig. 5a). (The ODE systems deviate as \(\omega _0\) increases, with (13) being marginally more accurate than the Grill et al. model for \(\omega _0=0.1\) (Fig. 7).) Further asymptotic reduction of (13) revealed the single algebraic equation describing the slow dynamics of the nonlinear relaxation oscillations and the associated amplitude of oscillation (Sect. 5.2, Fig. 6c–e).
While there is consistency between the descriptions in many respects, a striking feature is the appearance in stochastic simulations of noise-induced oscillations in a regime predicted to be linearly stable by the mean-field Fokker–Planck description (and the associated ODE system). The oscillations arise close to a stability boundary, and their period is well predicted by analysing the three linearised ODEs (the green circle in Fig. 5a highlights the position in parameter space occupied by the stochastic solution in Fig. 3a; the period prediction is given by (18b)). By restoring a representation of stochastic binding/unbinding kinetics (Sect. 6), we provide evidence that the amplitude of the oscillations is likely regulated by the noise associated with binding kinetics; the approximate SDE (38) captures amplitudes to within an order of magnitude (Fig. 8). The relationship between cell shape and division orientation was first described in the 1880s (Hertwig 1884) but it is inherently noisy (e.g. Nestor-Bergmann et al. (2019); Lam et al. (2020); Bosveld et al. (2016)). This makes trying to understand more complex regulation of spindle orientation, such as regulation by external forces, very challenging and can require experimental analysis of many spindles to obtain meaningful results. Our study reveals two underlying mechanisms of oscillation: small-amplitude spindle oscillations driven by noisy binding kinetics; and the previously described larger-amplitude oscillations driven by noise in force generator motion. Together, these may explain some of the noise seen in spindle and division orientation. For example, the predicted period of experimentally-observed oscillation can be used to infer (or at least constrain estimates of) some of the parameters relevant to the Xenopus system illustated in Fig. 1. For parameters as in Table 1, for a period of \(T\approx 100\) s as seen experimentally (Fig. 1d) then \(N=175\) force generators would be required. However, if restoring forces were reduced to \(K=0.005\) then \(N=21\) force generators would be required to achieve the same period. These would result in oscillations with an amplitude of \(\approx 15~\upmu \)m which is comparable with the typical size of a cell (\(\approx \)20 \(\upmu \)m diameter) though approximately three to four times larger than what was recorded experimentally (Fig. 1d, \(\approx 2\)–\(5\,\upmu \)m). Further work is needed to refine assessments of parameters to allow more direct comparison between theory and experiment.
Figure 5c, d shows how increasing demographic stochasticity by increasing the value of diffusive terms \(D_{\text {b}}\) and \(D_{\text {u}}\) can promote oscillations. Thus, noise associated with movement of force generators increases the ease with which oscillations are sustained. Expansion of regions of parameter space giving rise to oscillatory solutions under the addition of noise has been seen elsewhere, in studies of oscillations in protein expression. For example, Phillips et al. (2016) show that stochasticity increases the robustness of the oscillatory phenotype of gene expression resulting in the correct timing of cell differentiation. Indeed, if oscillations of the spindle pole play an important role in correctly orientating the mitotic spindle, then the inherent stochasticity of biological systems due to fluctuations in protein levels or ATP availability (driving dynein movement) would aid the robustness of correct spindle orientation.
Relaxation oscillations arise when restoring forces decrease relative to pulling forces (Fig. 6c–e). Biologically, this may be achieved by reducing the restoring force, for example by hinging of microtubules at the spindle pole (Howard 2006; Rubinstein et al. 2009), or by increased numbers of dynein linkers at the cell cortex (denoted in this work by an increase in N). The resulting linear sections of the spindle pole oscillation align temporally with slow phases in the time-evolution of \(P_{\text {b}}^\pm \) and \(y_{\text {c}}^\pm \) (Fig. 6c–e), until the spindle pole is sufficiently displaced from the centre for the restoring force to create a rapid reversal in the spindle pole velocity \(z_{t}\). This change in the spindle pole velocity results in a rapid increase in the value of \(y_{\text {c}}\) on the opposite cortex, which in turn creates a rapid decrease in the value of \(P_{\text {b}}^\pm \) due to the tension-sensitivity of the unbinding rate. \(P_{\text {b}}^\pm \) and \(y_{\text {c}}^\pm \) have amplitudes that are self-limited by the tension-sensitive unbinding, and amplified by their connection to the motion of the spindle pole. We have shown that \(\gamma >1\) is required for oscillations to occur at all (e.g in (22), (24)); dynein’s slip–bond (Ezber et al. 2020) is crucial for oscillatory dynamics, but the sensitivity of the slip–bond to tension may affect the nonlinearity of the oscillation.
In the limit of small restoring forces, the model may be simplified to a single algebraic equation describing the slow phases of the limit cycle in (\(\tilde{{\tilde{z}}}_0, \tilde{{\tilde{z}}}_{0,\tilde{{\tilde{t}}}}\)) and subsequently in (\(z_0,y_{\text {c}}^\pm \)), where the maximum amplitude of oscillation can be estimated using (32). This relation has similarities with the bistable force-velocity relation derived by Schwietert and Kierfeld (2020) in their model of kinetochore-chromosome dynamics, which underlies relaxation oscillations in that system. Re-dimensionalisation of the amplitude seen during relaxation oscillations (Fig. 6c–e) results in an oscillation with an amplitude of order \(\approx 0.1\) mm (for chosen baseline parameters), which is an order of 10 larger than the typical size of a cell in the Xenopus epithelium (\(\approx \)20 \(\upmu \)m diameter). This large amplitude is an artefact of the linear force-displacement law (1) and the 1D description; the imposed geometry necessary for a 2D description would allow the redistribution of pulling forces away from the direction of motion upon close proximity of the spindle pole as in Wu et al. (2024). Indeed, without pushing forces a redistribution of pulling forces is sufficient to cause reversal of spindle motion and relaxation oscillations (Wu et al. 2024). This effect could be described by a nonlinear restoring force at the boundary in 1D.
While we have studied a 1D model, imaged mitotic spindles in epithelial cells show 2D dynamics (Fig. 1; Online Resource 1), with forces acting on both spindle poles originating from the entire cell periphery. To properly consider the spindle dynamics and infer Xenopus system specific parameters, the relative motion of the two spindle poles must become part of the equation. The simplification of the present model to ODEs or low-dimensional SDEs is a key step to fully modelling 2D movements of a full mitotic spindle. From there, we may begin to piece together the full processes by which the mitotic spindle is positioned and orientated in tissue-based cells. In doing so, we must remain mindful that inherent stochasticity may increase the mobility of the spindle.
8 Supplementary Information
Online Resource 1: Time–lapse of a dividing cell in the epithelium of a Xenopus laevis embryo at stage 10. The mitotic spindle is seen in green (GFP-\(\alpha \)-tubulin) and metaphase plate in magenta (mCherry-Histone 2B). Images taken every \(t=6.0~\)s.
Online Resource 2: a) Movie of the evolution of fluxes \(J^+_{\text {b}}\), \(J^+_{\text {u}}\) and their sum \(J^+_{\text {b}}+J^+_{\text {u}}\). b) as in a) with a truncated y-axis to better demonstrate the dynamics of \(J^+_{\text {b}}\) and \(J^+_{\text {b}}+J^+_{\text {u}}\). c) The evolution of the pdfs \(P^+_{\text {u}}\) and \(P^+_{\text {b}}\). Parameter values as in Table 1 but with \(D_{\text {b}}=0.008\), \(D_{\text {u}}=0.004\) and \(N=25\).
Availability of Data and Materials
Biological data available upon request.
Code Availability
References
Belyy V, Hendel NL, Chien A et al (2014) Cytoplasmic dynein transports cargos via load-sharing between the heads. Nat Commun 5(1):1–9
Bergstralh DT, St Johnston D (2014) Spindle orientation: what if it goes wrong? In: Seminars in cell developmental biology. Elsevier, pp 140–145
Beta C, Kruse K (2017) Intracellular oscillations and waves. Annu Rev Condens Matter Phys 8:239–264
Boland RP, Galla T, McKane AJ (2008) How limit cycles and quasi-cycles are related in systems with intrinsic noise. J Stat Mech Theory Exp 09:P09001
Bosveld F, Markova O, Guirao B et al (2016) Epithelial tricellular junctions act as interphase cell shape sensors to orient mitosis. Nature 530(7591):495–498
Burakov A, Nadezhdina E, Slepchenko B et al (2003) Centrosome positioning in interphase cells. J Cell Biol 162(6):963–969
Compton DA, Cleveland DW (1993) NuMA is required for the proper completion of mitosis. J Cell Biol 120(4):947–957
Corrigan AM, Shrestha R, Draviam VM et al (2015) Modeling of noisy spindle dynamics reveals separable contributions to achieving correct orientation. Biophys J 109(7):1398–1409
Dogterom M, Kerssemakers JW, Romet-Lemonne G et al (2005) Force generation by dynamic microtubules. Curr Opin Cell Biol 17:67–74
Erban R, Chapman SJ (2020) Stochastic modelling of reaction–diffusion processes, vol 60. Cambridge University Press, Cambridge
Ezber Y, Belyy V, Can S et al (2020) Dynein harnesses active fluctuations of microtubules for faster movement. Nat Phys 16(3):312–316
Gardiner CW (1985) Handbook of stochastic methods. Springer, Berlin
Gillespie DT (1977) Exact stochastic simulation of coupled chemical reactions. J Phys Chem 81(25):2340–2361
Goddard GK, Tarannum N, Woolner S (2020) Applying tensile and compressive force to xenopus animal cap tissue. Cold Spring Harbor Protoc 3:pdb-prot105551
Grill S, Kruse K, Jülicher F (2005) Theory of mitotic spindle oscillations. Phys Rev Lett 94(10):108104
Harborth J, Weber K, Osborn M (1995) Epitope mapping and direct visualization of the parallel, in-register arrangement of the double-stranded coiled-coil in the NuMA protein. EMBO J 14(11):2447–2460
Hargreaves D (2023) Mitotic spindle dynamics in stretched epithelial tissue in vivo and in silico. PhD thesis, University of Manchester, UK
Hertwig O (1884) Das Problem der Befruchtung und der Isotropie des Eies: eine Theorie der Vererbung, vol 18. Verlag von Gustav Fischer
Howard J (2006) Elastic and damping forces generated by confined arrays of dynamic microtubules. Phys Biol 3(1):54–66
Joshi SD, Davidson LA (2010) Live-cell imaging and quantitative analysis of embryonic epithelial cells in Xenopus laevis. JoVE 39:e1949
Karsenti E, Vernos I (2001) The mitotic spindle: a self-made machine. Science 294(5542):543–547
Lam MS, Lisica A, Ramkumar N et al (2020) Isotropic myosin-generated tissue tension is required for the dynamic orientation of the mitotic spindle. Mol Biol Cell 31(13):1370–1379
Larson ME, Bement WM (2017) Automated mitotic spindle tracking suggests a link between spindle dynamics, spindle orientation, and anaphase onset in epithelial cells. Mol Biol Cell 28(6):746–759
Li J, Jiang H (2017) Geometric asymmetry induces upper limit of mitotic spindle size. Biophys J 112(7):1503–1516
Li J, Cheng L, Jiang H (2019) Cell shape and intercellular adhesion regulate mitotic spindle orientation. Mol Biol Cell 30(19):2458–2468
Milo R, Phillips R (2015) Cell biology by the numbers. Garland Science, New York
Minc N, Burgess D, Chang F (2011) Influence of cell geometry on division-plane positioning. Cell 144(3):414–426
Mitchison T, Salmon E (2001) Mitosis: a history of division. Nat Cell Biol 3(1):E17–E21
Morin X, Bellaïche Y (2011) Mitotic spindle orientation in asymmetric and symmetric cell divisions during animal development. Dev Cell 21(1):102–119
Nestor-Bergmann A, Goddard G, Woolner S et al (2018) Relating cell shape and mechanical stress in a spatially disordered epithelium using a vertex-based model. Math Med Biol 35(Supplement-1):i1–i27
Nestor-Bergmann A, Stooke-Vaughan GA, Goddard GK et al (2019) Decoupling the roles of cell shape and mechanical stress in orienting and cueing epithelial mitosis. Cell Rep 26(8):2088–2100
Okumura M, Natsume T, Kanemaki MT et al (2018) Dynein-dynactin-NuMA clusters generate cortical spindle-pulling forces as a multi-arm ensemble. Elife 7:e36559
Pecreaux J, Röper JC, Kruse K et al (2006) Spindle oscillations during asymmetric cell division require a threshold number of active cortical force generators. Curr Biol 16(21):2111–2122
Pecreaux J, Redemann S, Alayan Z et al (2016) The mitotic spindle in the one-cell C. elegans embryo is positioned with high precision and stability. Biophys J. 111(8):1773–1784
Phillips NE, Manning CS, Pettini T et al (2016) Stochasticity in the mir-9/hes1 oscillatory network can account for clonal heterogeneity in the timing of differentiation. Elife 5:e16118
Rubinstein B, Larripa K, Sommi P et al (2009) The elasticity of motor-microtubule bundles and shape of the mitotic spindle. Phys Biol 6(1):016005
Schneider CA, Rasband WS, Eliceiri KW (2012) Nih image to imagej: 25 years of image analysis. Nat Methods 9(7):671–675
Schwietert F, Kierfeld J (2020) Bistability and oscillations in cooperative microtubule and kinetochore dynamics in the mitotic spindle. New J Phys 22(5):053008
Tang Z, Hu Y, Wang Z et al (2018) Mechanical forces program the orientation of cell division during airway tube morphogenesis. Dev Cell 44(3):313–325
Théry M, Jiménez-Dalmaroni A, Racine V et al (2007) Experimental and theoretical study of mitotic spindle orientation. Nature 447(7143):493–496
Trokter M, Mücke N, Surrey T (2012) Reconstitution of the human cytoplasmic dynein complex. Proc Nat Acad Sci 109(51):20895–20900
Tsimring LS (2014) Noise in biology. Rep Prog Phys 77(2):026601
Valfort AC, Launay C, Sémon M et al (2018) Evolution of mitotic spindle behavior during the first asymmetric embryonic division of nematodes. PLoS Biol 16(1):e2005099
Woolner S, Miller AL, Bement WM (2010) Imaging the cytoskeleton in live Xenopus laevis embryos. Cytoskelet Methods Protoc 586:23–39
Wu HY, Kabacaoğlu G, Nazockdast E et al (2024) Laser ablation and fluid flows reveal the mechanism behind spindle and centrosome positioning. Nat Phys 20(1):157–168
Zhu J, Burakov A, Rodionov V et al (2010) Finding the cell center by a balance of dynein and myosin pulling and microtubule pushing: a computational study. Mol Biol Cell 21(24):4418–4427
Acknowledgements
The authors have applied a Creative Commons Attribution (CCBY) licence to any Author Accepted Manuscript version arising. This work was supported by the Wellcome Trust (098390/Z/12/Z and 225408/Z/22/Z). DH was supported by a Wellcome Trust PhD studentship (220054/Z/19/Z). OEJ and SW acknowledge support from the Leverhulme Trust (RPG-2021-394). The authors are grateful to Simon Cotter for technical advice and to anonymous reviewers for helpful suggestions.
Author information
Authors and Affiliations
Contributions
All authors contributed to the study conception and design. Mitotic spindle movie acquisition was performed by SW. Imaging of the metaphase plate and data analysis was performed by DH. The first draft of the manuscript was written by DH and all authors commented on subsequent versions of the manuscript. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethical Approval
All work with Xenopus laevis was performed using protocols approved by the UK Government Home Office under the Home Office Project Licence PFDA14F2D (Holder: Professor Enrique Amaya) and Home Office Personal Licences held by SW and DH.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
Below is the link to the electronic supplementary material.
Supplementary file 1 (avi 753 KB)
Supplementary file 2 (mp4 3733 KB)
Appendices
Appendix A: Acquisition of Biological Data
Xenopus laevis Xenopus laevis male and female frogs were housed within tanks maintained by the in-house animal facility at the University of Manchester. Female frogs were pre-primed 4–7 days in advance of embryo collection by injection with 50 U of pregnant mare serum gonadotrophin (Intervet UK) into the dorsal lymph sac. One day prior to embryo collection, male and pre-primed female frogs were primed by injection with 100 U (male) and 200 U (female) of human chorionic gonadotrophin (hCG; Chorulon, MSD) into the dorsal lymph sac. 2–5 h ahead of embryo collection, primed male and female frogs were transferred into the same tank for amplexus. Embryos were collected over 1 h time periods. Embryos were dejellied with 2% L-cysteine solution (Sigma Aldrich, #168149-100 G) in 0.1X Marks Modified Ringers (MMR) [1X MMR: 100 mM NaCl, 2 mM KCl, 1 mM MgCl2, 1 mM CaCl2, 0.5mM EDTA and 5 mM HEPES, pH7.8], rinsed with \(0.1\times \) MMR and incubated at room temperature (RT) to reach two-cell stage.
All Xenopus work was performed using protocols approved by the UK Government Home Office and covered by Home Office Project License PFDA14F2D (License Holder: Professor Enrique Amaya) and Home Office Personal Licenses held by SW and DH.
1.1 A.1 Whole Embryo Movies, Spindle and Metaphase Plate
For live imaging of mitotic spindles in Xenopus embryos (Fig. 1a, b), both cells of two-cell embryos were microinjected with 5 nl of mRNA for GFP-\(\alpha \)-tubulin (needle concentration of 0.5 g/l) and mCherry-Histone 2B (0.1 g/l), to highlight spindle microtubules and chromosomes, respectively. Embryos were incubated for 20 h post fertilization at \(\hbox {16}\,^\circ \)C and then mounted for live imaging in \(0.1\times \) Marks Modified Ringers (MMR) [\(10\times \) solution: 1 M NaCl, 20 mM KCl, 10 mM \(\hbox {MgCl}_2\).6\(\hbox {H}_2\)O, 20 mM \(\hbox {CaCl}_2\).2\(\hbox {H}_2\)O, 1 mM EDTA disodium salt, 50 mM Hepes, up to 5 L with distilled water], using a ring of vacuum grease to contain the embryos and support a glass coverslip as in Woolner et al. (2010). Imaging took place at developmental stages 10–11. Single focal plane live-cell images of spindles were collected at RT (\(\hbox {21}\,^\circ \)C) every 6 s using a confocal microscope (FluoView FV1000; Olympus) with FluoView acquisition software (Olympus) and a \(60\times \), 1.35 NA U Plan S Apochromat objective. Time-lapse videos were constructed from the single focal plane images using ImageJ.
1.2 A.2 Animal Cap Movies, Metaphase Plate Only
For imaging in the animal cap to follow dynamics of spindle movements (Fig. 1c–f), both cells of two-cell embryos were microinjeceted with mRNA for mCherry-Histone 2B (0.1 g/l; to highlight chromosomes) and BFP-CAAX (0.1 g/l; to highlight cell edges) using a Picospritzer III Intracel injector (Parker instrumentation). Injected embryos were washed in 0.1% MMR and incubated in fresh 0.1% MMR overnight at \(16\,^\circ \)C. Animal cap explants were prepared from the injected embryos at the early gastrula stage (stage 10). The embryos were transferred to Danilchik’s for Amy explant culture media (DFA) [53 mM NaCl, 5 mM \(\hbox {Na}_2\) \(\hbox {CO}_3\), 4.5 mM Potassium gluconate, 32 mM Sodium gluconate, 1 mM \(\hbox {CaCl}_2\), 1 mM \(\hbox {MgSO}_4\), up to 1 L with MilliQ water, pH 8.3 with Bicine] in 0.1% BSA (Sigma, A7906). The vitelline membranes were removed from the embryos using forceps, and the explant removed by incisions with the forceps around the animal pole resulting in separation of the animal cap tissue from the embryo (Joshi and Davidson 2010). The animal caps were then transferred onto a fibronectin-coated PDMS membrane with the basal side in contact with the membrane to prevent balling up and held in place with a coverslip. The caps were then left to recover for 2 h at \(18\,^\circ \)C before imaging (Goddard et al. 2020).
The PDMS membranes were prepared as described previously (Goddard et al. 2020). The PDMS membrane was mounted onto a stretch apparatus (custom made by Deben UK 722 Limited) and subjected to a uniaxial stretch of 0.5 mm displacement to ensure that the membrane remained taut under gravity and the weight of the animal cap (Nestor-Bergmann et al. 2018). Images were acquired every 5 s on a Leica TCS SP8 AOBS upright confocal using a \(20\times \) dipping objective at 2X confocal zoom. The confocal settings were as follows: pinhole 1.9 airy unit, 600 Hz bidirectional scanning, format \(1024\times 1024\). Images were collected using hybrid detectors with the detection mirror settings for red and blue at 600–690 nm, and 415–516 nm respectively, using the white light laser with excitation at 586 nm, and 405 nm laser lines. The images were collected non-sequentially with a z-spacing of 10 \({\upmu }\)m between sections. Images were taken continuously without resetting for drifting in order to ensure no missing data points for dividing cells. The maximum imaging duration per animal cap was 2 h.
1.3 A.3 Spindle Movement Analysis
All images were processed and analysed on ImageJ (Schneider et al. 2012). All measurements were taken from maximum intensity projections of the z-stack images. Each end of the metaphase plate was tracked in each frame using the ImageJ multi-point tool, returning an x and y coordinate for each point. The centre of the cell at the beginning of metaphase, \({\textbf{R}}_1\), and the end of metaphase, \({\textbf{R}}_2\), were used to create a linear correction to the metaphase plate position across metaphase time. Cell edges and tricellular vertices were manually traced at the beginning and end of metaphase using the ImageJ ‘Paintbrush tool’ (brush width = 1 pixel). The manual traces were processed using in-house Python scripts to return the cell centre position (Nestor-Bergmann et al. 2018, 2019). The measurements were made based on the polygonised cell according to the positions of the tricellular vertices. Oscillations were detected from signals using a periodogram.
Appendix B: Implementation of the Gillespie Algorithm
The extension of an elastic linker is discretised into states \(y^{\left( n\right) \pm ,i}_{\text {b(u)}}\) with \(i=0,1\ldots M\), separated by a fixed distance \(\Delta {y}\) such that \(y^{\left( n\right) \pm ,i+1}_{\text {b(u)}} = y^{\left( n\right) \pm ,i}_{\text {b(u)}} + \Delta {y}\). Each force generator n has identifiers which denote the associated cortex (±), the current extension state (i), and the binding state (\(\text {u}\) for unbound, \(\text {b}\) for bound). The binding state will be identified in the subscript and written as \(\text {b(u)}\), referring to a subscript b or u.
At any time, a generator may
-
retract: \({y}^{\left( n\right) \pm ,i}_{\text {b(u)}}\rightarrow {y}^{\left( n\right) \pm ,i-1}_{\text {b(u)}}\) with probability \(r^{\left( n\right) \pm ,i}_{\text {b(u)}}\);
-
extend: \(y^{\left( n\right) \pm ,i}_{\text {b(u)}}\rightarrow {y}^{\left( n\right) \pm ,i+1}_{\text {b(u)}}\) with probability \(f^{\left( n\right) \pm ,i}_{\text {b(u)}}\); or
-
switch between bound and unbound states: \({y}^{\left( n\right) \pm ,i}_{\text {b}}\leftrightarrow {y}^{\left( n\right) \pm ,i}_{\text {u}}\) with probability \(s^{\left( n\right) \pm ,i}_{\text {b(u)}}\).
These state-changing events are illustrated graphically in Fig. 9a.
Probabilities \(r^{\left( n\right) \pm ,i}_{\text {b(u)}}\), \(f^{\left( n\right) \pm ,i}_{\text {b(u)}}\), and \(s^{\left( n\right) \pm ,i}_{\text {b(u)}}\) are related to model parameters as follows. The switching probabilities were chosen such that an unbound generator may switch to become a bound generator within a short time \(\tau \) with a probability \(s^{\left( n\right) \pm ,i}_{\text {u}}=\tau \omega _{\text {on}}\) for a constant binding rate \(\omega _{\text {on}}\). In a short time \(\tau \), a bound generator may unbind with probability \(s^{\left( n\right) \pm ,i}_{\text {b}}=\tau \omega _0e^{\gamma y^{\left( n\right) \pm ,i}_{\text {b}}}\).
In order to obtain expressions for \(r^{\left( n\right) \pm ,i}_{\text {b(u)}}\) and \(f^{\left( n\right) \pm ,i}_{\text {b(u)}}\), consider
as an effective drift speed and diffusion coefficient for force generators respectively. These arise from considering extension or contraction of each linker as a biased random walk.
We summarise these transition probabilities as
No flux conditions were enforced by setting \(r^{(n)\pm ,i=0}_{\text {b(u)}}=0\) and \(f^{(n)\pm ,i=M}_{\text {b(u)}}=0\).
The Gillespie algorithm (Gillespie 1977) stipulates that the probability of a state-changing event (extension, retraction, or switch) happening within a short time \(\tau \) is exponentially distributed with rates \(r_{\text {b(u)}}^{\left( n\right) \pm ,i}/\tau \), \(f_{\text {b(u)}}^{\left( n\right) \pm ,i}/\tau \), and \(s_{\text {b(u)}}^{\left( n\right) \pm ,i}/\tau \), which sum together to give a total rate
Here 2N is the total number of force generators within the system (N per cortex), each of which is associated with either the upper (+) or lower (−) cortex, has an extension state i, and is either bound (b) or unbound (u). We assume that only one event for one force generator may occur, removing the possibility of simultaneous events. As \(r_{\text {b(u)}}^{\left( n\right) \pm ,i}\), \(f_{\text {b(u)}}^{\left( n\right) \pm ,i}\) and \(s_{\text {b(u)}}^{\left( n\right) \pm ,i}\) are proportional to the short time \(\tau \) (B2), the rates \(r_{\text {b(u)}}^{\left( n\right) \pm ,i}/\tau \), \(f_{\text {b(u)}}^{\left( n\right) \pm ,i}/\tau \), and \(s_{\text {b(u)}}^{\left( n\right) \pm ,i}/\tau \) (and thus R) are independent of \(\tau \). A random variable \(\zeta _1\) is chosen from a uniformly random distribution between 0 and 1 (\(\zeta _1\sim {\mathcal {U}}\left[ 0,1\right] \)) and the time to the next event is calculated using \(\tau ={R}^{-1}\log {\left( 1/\zeta _1\right) }\). The rescaled rates \(a^{(n)\pm ,r}_{\text {b(u)}}(i) = R^{-1}r_{\text {b(u)}}^{(n)\pm ,i}/\tau \), \(a^{(n)\pm ,f}_{\text {b(u)}}(i) = R^{-1}f_{\text {b(u)}}^{(n)\pm ,i}/\tau \) and \(a^{(n)\pm ,s}_{\text {b(u)}}(i) = R^{-1}s_{\text {b(u)}}^{(n)\pm ,i}/\tau \) are concatenated in triplets for each force generator n, giving a list of potential states \(a^j\) with \(j\in [1,6N]\) which sum together to give \({{\sum }}_{j=1}^{6N}a^{j}=1\) (Fig. 9b). Choosing an independent random variable from a uniformly random distribution, \(\zeta _2\sim {\mathcal {U}}\left[ 0,1\right] \), the next state-changing event is determined as the first j such that \({{\sum }}_{j'=1}^ja^{j'}>\zeta _2\). Force generators in the upper (\(n^+\)) and lower (\(n^-\)) cortex have corresponding events \(a^j\) where \(j\in [1,3N]\) and \(j\in [3N+1,6N]\) respectively.
In order to calculate the spindle pole position, we implement a forward Euler approximation of (4) which may be used to calculate the pole position at a time \(t+\tau \),
Here \(n'\) and \(N'\) are the equivalent of n and N, introduced to separate the upper and lower cortex in this expression.
In summary, the state of the system at any instant can be described by a vector \({\textbf{X}}(t)\) of size \({\mathcal {N}}=1+4N(M+1)\), comprising the spindle location plus the occupancies of 4N linkers (N at each cortex, in bound or unbound states) in states of different lengths (over a scale discretized into \(M+1\) elements). In principle such a system can be represented (Erban and Chapman 2020) by a chemical master equation for \(p({\textbf{x}},t)\), the probability that \({\textbf{X}}(t)={\textbf{x}}\), coupled (Langevin) stochastic differential equations for the elements of \({\textbf{X}}\) and a chemical Fokker–Planck equation for \(p({\textbf{x}},t)\). The latter is a PDE of \({\mathcal {N}}+1\) dimensions. This is distinct from the heavily reduced Fokker–Planck system (7,10) proposed by Grill et al. (2005) that motivates the stochastic system illustrated in Fig. 9.
Appendix C: Reducing the Fokker–Planck Equations to ODEs
To reduce the Fokker–Planck model to a system of ODEs, we rescale using the motor-protein-to-microtubule binding rate, writing \({t}={\tilde{t}}/{\omega }_{\text {on}}\) and \({z}={\tilde{z}}/{\omega }_{\text {on}}\). Then (10) and (7) become
We develop an approximation to the oscillating spindle system for which \({\omega }_{\text {on}}\sim {\omega }_{0}\sim {D_{\textrm{b}}}^{1/2}\sim {D_{\textrm{u}}}^{1/2}\ll 1\). To minimise the introduction of further notation, we expand our solutions in terms of the small order parameter \({\omega }_{\text {on}}\) and remain mindful that these parameters are taken to be of similar order. The range of extension values y is split into three regions (Fig. 4h): region I over which \({P}_{\text {u}}^\pm \) is peaked around \({y}=0\) with a width \({D_{\textrm{u}}}^{1/2}\); region III over which \({P}_{\text {b}}^\pm \) is peaked with a width of \({D_{\textrm{b}}}^{1/2}\) but whose centre moves as \({y}_{\text {c}}=1\mp {\tilde{z}}_{{\tilde{t}}}\); and region II where advective terms dominate and the asymptotic limits of I and III are matched. Solutions for \({P}_{\text {u}}^\pm \) and \({P}_{\text {b}}^\pm \) will be determined in regions I and III respectively, followed by matching their asymptotic limits in region II to reveal the ODE system which governs the time evolution of the parameters.
1.1 C.1 Region I
In Region I, we seek solutions \({P}_{\text {u}}^\pm \sim {P}_{\text {u}0}^\pm +{\omega }_{\text {on}}{P}_{\text {u}1}^\pm +\cdots \) where \({P}^\pm _{\text {u}0}\) is a quasi-static solution whose shape is static but whose amplitude varies slowly in time. We assume further that \({P}^\pm _{\text {b}}\sim {\omega }_{\text {on}}\) in this region (Fig. 4h). Here, it is observed that \({P}_{\text {u}}^\pm \) is sharply peaked about \({y}=0\) over a diffusive length-scale \({D_{\textrm{u}}}^{1/2}\) (Fig. 4h). Thus, setting \({y}={D_{\textrm{u}}}^{1/2}Y\) in (C5c) gives
with the boundary condition \({J}_{\text {u}}^\pm \left( {\tilde{t}},0\right) =0\). This boundary condition therefore becomes
where
which are both individually zero at \(Y=0\) due to (C7). To leading order in \({\omega }_{\text {on}}\), (C6) becomes
which may be integrated to give \(\Gamma \left[ Y{P}^\pm _{\text {u0}}+{P}^\pm _{\text {u0},Y}\right] ^Y_0=0\). Thus, due to boundary condition (C7),
which gives
as a solution, with \(A^\pm \left( {\tilde{t}}\right) \) an amplitude which varies slowly in time.
At \({\mathcal {O}}\left( {\omega }_{\text {on}}\right) \), (C6) becomes
We highlight here the absence of the \({\omega }_0e^{\gamma {y}}{P}_{\text {b}}^\pm \) term as we have assumed that \({P}_{\text {b}}^\pm \sim {\omega }_{\text {on}}\) in this region. Then (C12) may be integrated to
using the boundary conditions on flux (C7) at \(Y=0\). As the \(Y{P}_{\text {u1}}^\pm \) term will dominate the left-hand side as \(Y\rightarrow \infty \), we may write
Then,
Re-substitution of \(Y={D_{\textrm{u}}}^{-1/2}{y}\) gives
when \({D_{\textrm{u}}}^{1/2} \ll {y}\).
Now that we have an expression for how the shape of the unbound force generator pdf varies in time in region I, we seek similar solutions for the bound generator pdf in region III.
1.2 C.2 Region III
In region III we seek solutions of the form \({P}_{\text {b}}^\pm \sim {P}_{\text {b}0}^\pm +{\omega }_{\text {on}}{P}_{\text {b}1}^\pm +\cdots \). Here, the pdf \({P}_{\text {b}}^\pm \) is sharply peaked about \({y}_{\text {c}}=1\mp {\tilde{z}}_{{\tilde{t}}}\) over a diffusive length-scale \({D_{\textrm{b}}}^\frac{1}{2}\) (Fig. 4h). Both \({D_{\textrm{b}}}^{1/2}\) and \({\omega }_{\text {on}}\) are assumed to be small parameters of similar order. Thus, in this region about the peak of \({P}_{\text {b}}^\pm \), we set \({y}=1\mp {\tilde{z}}_{{\tilde{t}}}+{D_{\textrm{b}}}^\frac{1}{2}{\hat{Y}}\). Noting that
(C5b) becomes, to leading order,
Recalling (8), (C18) may be further simplified to
Substituting the expansion \({P}_{\text {b}}^\pm \approx {P}_{\text {b0}}^\pm +{\omega }_{\text {on}}{P}_{\text {b1}}^\pm +\cdots \), to first order (C19) becomes
taking \({{\omega }_{\text {on}}}/{{D_{\textrm{b}}}^{1/2}}\) to be \({\mathcal {O}}\left( 1\right) \). By integration,
for \(C\left( {\tilde{t}}\right) \) some constant of integration. The boundary condition \({J}^\pm _{\text {b}}\left( {\tilde{t}},{y}={y}_{\text {max}}\right) =0\) becomes \({J}^\pm _{\text {b}}\left( {\tilde{t}},{\hat{Y}}\rightarrow \infty \right) \rightarrow 0\), while \({J}^\pm _{\text {b}}\left( {\tilde{t}},{y}=0\right) =0\) becomes \({J}^\pm _{\text {b}}\left( {\tilde{t}},{\hat{Y}}\rightarrow -\infty \right) \rightarrow 0\). Then
where
which must both separately also tend to zero as \({\hat{Y}}\rightarrow -\infty \). Using this, (C21) may be rewritten as
Making the assumption that \({P}^\pm _{\text {b}0}\rightarrow 0\) as \({\hat{Y}}\rightarrow -\infty \), which enforces (C22) at leading order, then \(C\left( {\tilde{t}}\right) =0\) and
Integrating (C25) gives the solution
for some \({\tilde{B}}^\pm \left( {\tilde{t}}\right) \) and \(B^{\pm }\left( {\tilde{t}}\right) \), where \(B^\pm \left( {\tilde{t}}\right) \) describes the amplitude of the peak of the pdf (subject to smaller corrections) which varies in time. Returning to (C19), with \({P}^\pm _{\text {u}}\ll {P}^\pm _{\text {b}}\), then to \({\mathcal {O}}\left( {\omega }_{\text {on}}\right) \)
which may be rewritten as
Thus using (C26) in (C28) and integrating gives
By assuming that \({P}^\pm _{\text {b}1}\rightarrow 0\) as \({\hat{Y}}\rightarrow \infty \), which enforces boundary condition \({J}^\pm _{\text {b1}}\rightarrow 0\) as \({\hat{Y}}\rightarrow \infty \), then
The left-hand side (LHS) may be rewritten
and so in the limit \({\hat{Y}}\rightarrow -\infty \), (C29) is dominated by the \({\hat{Y}}{P}^\pm _{\text {b}1}\) term. Rearranging the right-hand side gives, in this limit,
The second integral vanishes, while the first integral can be evaluated and thus, as \({\hat{Y}}\rightarrow -\infty \),
The asymptotic limits (C16) and (C30) as \(Y\rightarrow \infty \) and \({\hat{Y}}\rightarrow -\infty \) respectively will now be matched inside region II.
1.3 C.3 Region II
In region II, advection terms dominate. These ‘sweep’ the bound force generators toward the peak of \({P}^\pm _{\text {b}}\) such that bound force generators will tend to have elastic linkers with an extension \({y}_{\text {c}}=1\mp {\tilde{z}}_{{\tilde{t}}}\), and the unbound force generators toward the peak of \({P}^\pm _{\text {u}}\) such that unbound force generators will tend to have an elastic linker with zero extension. Given that the pdfs are peaked in regions I and II, \({P}^\pm _{\text {b},{y}}\) and \({P}^\pm _{\text {u},{y}}\) are both relatively small in region II, being given by the small correction terms \({\omega }_{\text {on}}{P}^\pm _{\text {u1}}\) and \({\omega }_{\text {on}}{P}^\pm _{\text {b1}}\), expressions for which we have determined in the limits \(Y\rightarrow \infty \) and \({\hat{Y}}\rightarrow -\infty \) respectively. Then together, using (7a), \( {J}_{\text {b}}^\pm ={v}^\pm _{\text {b}}{P}_{\text {b}}^\pm -{D_{\textrm{b}}}{P}^\pm _{\text {b},{y}} \approx {v}^\pm _{\text {b}}{P}_{\text {b}}^\pm \), and so substitution of \({P}^\pm _{\text {b}}\approx {\omega }_{\text {on}}{P}^\pm _{\text {b1}}\) when \({\hat{Y}}\rightarrow -\infty \) and (8) returns
Continuing, using (7b), \( {J}_{\text {u}}^\pm = -\Gamma \left( {y}{P}_{\text {u}}^\pm +{D_{\textrm{u}}}{P}^\pm _{\text {u},{y}}\right) \approx -\Gamma {y}{P}^\pm _{\text {u}}\), and so substitution of \({P}^\pm _{\text {u}}\approx {\omega }_{\text {on}}{P}^\pm _{\text {u1}}\) when \(Y\rightarrow \infty \) returns
By the form of (C31) and (C32), \({J}^\pm _{\text {b},{y}}=0\) and \({J}^\pm _{\text {u},{y}}=0\) to leading order, and therefore (C31) and (C32) are valid across the whole of region II. It follows that
is a constant across region II. As demonstrated in Online Resource 2(a,b), detailed balance (\(J_{\textrm{b}}^\pm +J_{\textrm{u}}^\pm =0\)) does not hold in this region.
1.4 C.4 Matching Solutions: Regions I-II
In region I it was assumed that \(P_{\text {b}}^\pm \) was sufficiently small (\({\mathcal {O}}\left( \omega _{\text {on}}\right) \)) that its dynamics could be neglected to leading order. By (C31) it can be estimated that \(P^\pm _{\text {b}}\) is of magnitude \(\omega _{\text {on}}\sqrt{D_{\text {u}}}\), where it has been assumed that \(D_{\text {b}}\sim D_{\text {u}}\). Then letting \(P^\pm _{\text {b}}=\omega _{\text {on}}\sqrt{D_{\text {u}}}{\hat{P}}_{\text {b}}^\pm \) in (C5c) when \(y=\sqrt{D_{\text {u}}}Y\) results in expressions for \(P_{\text {u0}}^\pm \) and \(P^\pm _{\text {u1}}\) which are unchanged from (C11) and (C16) respectively. However, (C5b) becomes
To leading order, (C34) becomes
Then, integrating over Y from \(Y=0\) to \(Y\rightarrow \infty \),
Assuming that \(v_{\text {b}}^\pm \approx 1\mp {\tilde{z}}_{{\tilde{t}}}\) to leading order, then
Since \(J_{\text {b}}^\pm \) is independent of y in region II (Online Resource 2a,b), then we match (C31) to (C37), resulting in
We now perform the same analysis on the boundary of regions II and III.
1.5 C.5 Matching Solutions: Regions II-III
In region III it was assumed that \(P_{\text {u}}^\pm \) was sufficiently small that its dynamics could be neglected to leading order. By (C32) it can be estimated that \(P^\pm _{\text {u}}\) is of magnitude \(\omega _{\text {on}}\sqrt{D_{\text {u}}}\). Then letting \(P^\pm _{\text {u}}=\omega _{\text {on}}\sqrt{D_{\text {u}}}{\hat{P}}_{\text {u}}^\pm \) in (C5b) when \(y=1\mp {\tilde{z}}_{{\tilde{t}}}+\sqrt{D_{\text {b}}}{\hat{Y}}\) results in expressions for \(P_{\text {b0}}^\pm \) and \(P^\pm _{\text {b1}}\) which are unchanged from (C26) and (C30) respectively. However, (C5c) becomes
To leading order,
where we have used that \(e^{\gamma \left( 1\mp {\tilde{z}}_{{\tilde{t}}}+\sqrt{D_{\text {b}}}{\hat{Y}}\right) }=e^{\gamma \left( 1\mp {\tilde{z}}_{{\tilde{t}}}\right) }e^{\gamma \sqrt{D_{\text {b}}}{\hat{Y}}}\approx e^{\gamma \left( 1\mp {\tilde{z}}_{{\tilde{t}}}\right) }\) as \(\sqrt{D_{\text {b}}}\) is a small parameter. Integrating from \({\hat{Y}}\rightarrow -\infty \) to \({\hat{Y}}\rightarrow \infty \) gives
Taking \(y\approx 1\mp {\tilde{z}}_{{\tilde{t}}}\) to leading order, then
Since \(J_{\text {u}}^\pm \) is independent of y in region II (Online Resource 2a,b), then we match (C32) to (C42) resulting in
1.6 C.6 Combining the Whole System
We may now use expressions (C11), (C16) for \({P}_{\text {b}}^\pm \) and (C26), (C30) for \({P}_{\text {u}}^\pm \), and their coupling in region II (C38, C43) to close the system. Recalling (12), then to leading order
The first term of this integral is easily evaluated, while the second term is more complex. Consider only the leading-order terms of the exponent, due to \({\omega }_{\text {on}}\) being a small order parameter. Then
which we know is a peak contained within region III. That is, we integrate over the Gaussian, which does not intersect \({y}=0\) with any value of significance at leading order. Using this logic, the integral (C45) may be evaluated and thus
This can be used to eliminate \(A^\pm \) from (C38) to give
(C46) may similarly be used in (C43) to return (C47).
Equation (C47) predicts that \(\sqrt{2\pi {D_{\textrm{b}}}}B^\pm \) relaxes to \({{\omega }_{\text {on}}} / \{{\omega }_{\text {on}}+{\omega }_{0}e^{\gamma \left( 1\mp {\tilde{z}}_{{\tilde{t}}}\right) }\}\) and (C46) predicts that \(\sqrt{{\pi {D_{\textrm{u}}}}/{2}}A^\pm \) relaxes to \({{\omega }_{0}e^{\gamma \left( 1\mp {\tilde{z}}_{{\tilde{t}}}\right) }}/\{{\omega }_{\text {on}}+{\omega }_{0}e^{\gamma \left( 1\mp {\tilde{z}}_{{\tilde{t}}}\right) }\}\) provided \({\tilde{z}}\) does not change too rapidly.
Further to this, \({P}_{\text {b}}^\pm ={P}_{\text {b0}}^\pm +{\omega }_{\text {on}}{P}_{\text {b1}}^\pm +\cdots \) can be put into (C5a) to obtain a leading-order equation for the motion of the spindle pole. This requires the evaluation of
where we assume \({y}_{\text {max}}\) is sufficiently large that it exceeds the bounds of region III and can thus be taken as \({y}_{\text {max}}\rightarrow \infty \). Again we let \({y}={y}_{\text {c}}+{D_{\textrm{b}}}^{1/2}{\hat{Y}}\), which is where \({P}^\pm _{\text {b0}}\) has a significant value. Then
which to leading order becomes
Recalling that \({y}_{\text {c}}=1\mp {\tilde{z}}_{{\tilde{t}}}\), (C5a) becomes
and thus
This can be rewritten as
where \({\hat{B}}^\pm =\sqrt{2\pi {D_{\textrm{b}}}}B^\pm \), \({\hat{K}}={K}/\{{\omega }_{\text {on}}N\}\), and \({\hat{\xi }}={{\xi }}/{N}\). Recalling (C47) which may be alternatively written as
where \(\rho ={\omega }_0/{\omega }_{\text {on}}\), the coupled system (7a), (7b), and (10) is reduced to solving (C53, C54) along with initial conditions on \({\tilde{z}}_0\) and \({\hat{B}}^\pm _0\).
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Hargreaves, D., Woolner, S. & Jensen, O.E. Relaxation and Noise-Driven Oscillations in a Model of Mitotic Spindle Dynamics. Bull Math Biol 86, 113 (2024). https://doi.org/10.1007/s11538-024-01341-w
Received:
Accepted:
Published:
DOI: https://doi.org/10.1007/s11538-024-01341-w