Reverse engineering of force integration during mitosis in the Drosophila embryo.

The mitotic spindle is a complex macromolecular machine that coordinates accurate chromosome segregation. The spindle accomplishes its function using forces generated by microtubules (MTs) and multiple molecular motors, but how these forces are integrated remains unclear, since the temporal activation profiles and the mechanical characteristics of the relevant motors are largely unknown. Here, we developed a computational search algorithm that uses experimental measurements to 'reverse engineer' molecular mechanical machines. Our algorithm uses measurements of length time series for wild-type and experimentally perturbed spindles to identify mechanistic models for coordination of the mitotic force generators in Drosophila embryo spindles. The search eliminated thousands of possible models and identified six distinct strategies for MT-motor integration that agree with available data. Many features of these six predicted strategies are conserved, including a persistent kinesin-5-driven sliding filament mechanism combined with the anaphase B-specific inhibition of a kinesin-13 MT depolymerase on spindle poles. Such conserved features allow predictions of force-velocity characteristics and activation-deactivation profiles of key mitotic motors. Identified differences among the six predicted strategies regarding the mechanisms of prometaphase and anaphase spindle elongation suggest future experiments.

Introduction

Mitosis, the process by which identical copies of the replicated genome are distributed to the products of each cell division, involves a highly dynamic sequence of coordinated motility events mediated by a bipolar mitotic spindle (Karsenti and Vernos, 2001; Pines and Rieder, 2001; Scholey ) (Figure 1A). The motility is driven by forces generated by multiple molecular motors—kinesins and dyneins—together with dynamic microtubules (MTs) (Figure 1B) (Kline-Smith and Walczak, 2004; Brust-Mascher and Scholey, 2007) whose activities are controlled by kinases, phosphatases and proteases (O'Farrell, 2001; Parry and O'Farrell, 2001; Peters, 2002). These force-generating and regulatory proteins form a vast network that coordinates spindle assembly, maintenance and elongation, as well as orchestrating chromosome segregation.

Figure 1

The spindle protein machinery. (A) A cartoon that shows all major components of the spindle. Four MT populations (astral (as), kinetochore (kt), chromosome arm (chr) and inter-polar (ip) MTs) extend from the poles creating the spindle. Molecular motors bind to MTs and either regulate their ends' kinetics, or slide them, or exert forces on the chromosomes and centrosomes. (B) Eight possible MT–motor combinations, with the respective velocities and forces acting on a single MT. asMT: cortical dynein pulling the MT generates an outward force F1 on the spindle pole. chrMT are anchored at the pole, while MT polymerization and chromokinesins generate a pushing force F2; a force F3 is associated with such an MT if depolymerases are activated at the MT's minus end in addition to motor activity at the MT's plus end. ktMT: an inward force, F4 is generated on an MT anchored at the pole while kt motors act on the MT plus end; modified force F5 acts on an MT depolymerized at the minus end in addition to the plus end motor activity; force F6 is exerted if ktMT is depolymerized at its minus end and anchored at its plus end. ipMT: an outward force, F7 results from the combination of kinesin-5 and kinesin-14 actions on the MT anchored at the pole, while a force F8 is exerted by these motors on an MT being depolymerized at its minus end. (C) The experimentally measured time series for spindle length (pole–pole distance) in wild type (WT) and inhibited spindles used in the optimization process (details in Supplementary Figure 1, referenced in Supplementary Table 2). Colors correspond to motor colors in the legend. Previous studies revealed that the double inhibition of kinesin-5 and kinesin-14 fully rescues metaphase spindle assembly (dashed blue line, not used here) (Sharp ). However, recent studies suggest that this effect results from the partial inhibition of kinesin-5, whereas a more complete inhibition leads to prometaphase spindle collapse even in kinesin-14 null mutants (green line, used here), suggesting that an additional unknown inward force also opposes kinesin-5 and contributes to the collapse (Brust-Mascher and Scholey, in preparation).

In the Drosophila syncytial embryo, hundreds of mitotic spindles progress synchronously through a well-defined and reproducible sequence of transitions, in which periods of rapid pole–pole separation are interspersed with quiescent pauses (Sharp ). Each mitosis begins with prophase when the centrosomes located on the nuclear envelope separate with roughly hyperbolic kinetics to reach a steady-state separation distance of about 6–8 μm. Then following nuclear envelope breakdown at the onset of prometaphase, the spindle elongates further to reach another steady-state length of 10 μm in metaphase. Sister-chromatid segregation occurs during anaphase A and then the spindle undergoes a final linear episode of elongation to reach a final length of 14 μm in anaphase B (Sharp ; Brust-Mascher and Scholey, 2002).

Mitotic spindle dynamics depends upon the combined effect of several distinct molecular processes including, for example, force-generating mechanisms, changes in the concentration of MT components, the presence or absence of centrosomes, the establishment of morphogen gradients, etc. (Sharp ; Karsenti and Vernos, 2001; Mitchison and Salmon, 2001; Mitchison ). The extent to which these distinct processes influence mitotic spindle behavior appears to differ in different systems. However, the important idea that balances of antagonistic forces contribute to mitosis is thought to apply to a broad range of mitotic spindles. This idea was originally proposed by Ostergren (1951) to explain chromosome positioning and provided a plausible explanation for spindle pole dynamics in experimentally perturbed diatom spindles (Leslie and Pickett-Heaps, 1983). On the basis of observations of interactions among mutant genes encoding members of the yeast kinesin-5 and kinesin-14 families (Lawrence ), it was proposed that the corresponding motor proteins could exert antagonistic outward and inward forces on spindle poles, respectively (Saunders and Hoyt, 1992; Hoyt and Geiser, 1996). The idea that these counterbalancing forces are generated by an ‘antagonistic sliding filament mechanism' was supported by biochemical studies showing that purified kinesin-5 is a slow, plus-end-directed bipolar homotetramer capable of crosslinking and sliding apart antiparallel MTs, whereas kinesin-14 is a minus-end-directed homodimeric MT bundling motor that could slide antiparallel MTs inwards (McDonald ; Walker ; Sawin ; Cole ; Kashina ; Kapitein ; Oladipo ; Furuta and Toyoshima, 2008). The hypothesis that pole–pole separation depends upon a balance of forces generated by such an antagonistic sliding filament mechanism is further supported by observations that purified kinesin-5 and kinesin-14 can antagonize and balance one another in motility assays (Tao ), although further work is needed to establish if and how this mechanism contributes to spindle pole dynamics.

In the Drosophila syncytial embryo, a motor-generated force balance produced by systems of complementary and antagonistic motors is proposed to play a dominant role in spindle assembly/elongation and chromosome segregation (Sharp ). For example, comparison of the temporal changes in Drosophila embryo spindle length in five mutant and biochemically inhibited spindles reveals characteristic defects in pole–pole separation compared to wild-type (WT) spindles. Each defect can be most naturally explained by a shifting force balance in the spindle resulting from the inhibition of specific molecular motors (Figure 1C). For example, when either dynein or kinesin-5, presumably pulling the astral (as) MTs outward or sliding the inter-polar (ip) MTs outward, respectively, are inhibited, spindle length decreases. On the other hand, both kinesin-13 and -14 are hypothesized to contribute to the inward force on spindle poles by shortening the MTs connecting the poles and chromosomes and by sliding inward ipMTs, respectively, and when either of these motors is inhibited, spindle length increases.

Recently, mathematical modeling (reviewed in Gardner and Odde, 2006; Mogilner ) was used to examine the role of forces and MT dynamics in spindle development (Nédélec, 2002; Cytrynbaum ; Brust-Mascher ; Gardner ; Civelekoglu-Scholey ; Burbank ; Cheerambathur ). For example, one of these models explained pre-nuclear envelop breakdown spindle elongation to the steady state during prophase as being the result of a balance between a constant cortical dynein-generated force pulling asMTs outward and an antagonistic spindle length-dependent Ncd-generated force pulling the ipMTs that link the poles inward (Cytrynbaum ). Another model (Brust-Mascher ) quantitatively explained how the persistent sliding apart of ipMTs by kinesin-5 motors, combined with changes in the activity of an antagonistic MT depolymerase on the spindle poles, produces poleward flux in pre-anaphase B spindles and drives anaphase B spindle elongation. These models provided a successful description of the experimental data because of the relative biochemical simplicity of prophase, during which most mitotic motors are sequestered in the nucleus and do not contribute to the pole–pole separation mechanism, and the structural simplicity of anaphase B, during which the ipMTs dominate spindle pole dynamics. However, even these relatively simple models were based on guessing (plausible values of) a formidable number of parameters. When model complexity grows beyond certain level, intuition alone is insufficient. Reverse engineering, an approach complementary to this ‘explicit' modeling of a simplified system, uses computational optimization to automatically identify the appropriate model parameter values by constraining the parameters with quantitative experimental data.

Now that the inventory of mitotic force-generating and regulatory molecules is close to completion (e.g. Bettencourt-Dias ; Goshima ), the task of elucidating the mechanism of action of the mitotic spindle at a systems biology level is becoming realistic. Part of this task is to ‘reverse engineer' the spindle—i.e. using experimental data to understand the temporal activation sequence and the mechanical characteristics of the force generators acting during mitosis—and reconstituting the spindle in silico. This presents a challenge that seems prohibitive: with more than 10 molecular motors being involved, each characterized by unknown mechanical, kinetic and regulatory parameters, and the structural complexity of the spindle, especially in metaphase, it is impossible to use intuition and traditional modeling to explain the dynamics associated with the sequence of transitions characteristic of mitotic progression. To address this challenge, here we develop and utilize a novel computational algorithm that automatically builds force balance models from a few MT–motor modules and uses quantitative experimental data to screen and optimize them. This algorithm ultimately identifies all plausible activation sequences and mechanical characteristics of the molecular motors that mediate spindle elongation in the Drosophila embryo.

Results

Force balance model

We divide the complex spindle machinery into its elementary structural components based on four distinct MT populations that are known to act within the spindle, namely astral (asMT) inter-polar (ipMT), chromosomal (chrMT) and kinetochore (ktMT) (Figure 1B). For each MT population, we first calculate the force acting on a single MT. For example, to obtain the force acting on a single chrMT, we consider two possible scenarios shown in Figure 1B. In the first case, the chrMT's minus end is anchored at the pole, while its plus end is connected to the chromosome through the chromokinesin motor on the chromosome arm. Then, the relative velocity between the motor and the chrMT is the difference between the pole velocity, Vpole, and the chromosome velocity, Vchr. We characterize each mitotic motor with an assumed linear force–velocity relation (Nédélec, 2002; Cytrynbaum ) characterized by two parameters—maximal stall force and free unloaded velocity (in this case, Fchr,mx and Vchr,mx, respectively): if the MT does not move relative to the motor, then the motor pushes the MT in the poleward direction with the maximal, stall, force. However, if the MT slides poleward with the rate (Vpole–Vchr) relative to the chromokinesin motor, then the force exerted by this plus-end-moving motor on the chrMT is lower,

(there is zero force if Vpole−Vchr=Vchr,mx). In principle, multiple force generators can contribute to this force on the chromosome arm, including MT polymerization. In the model, we combine them into one single ‘composite' motor and assume that it can be characterized by a linear force–velocity relation.

An alternative possibility is that the chrMT is not anchored at its minus end, but rather is being actively depolymerized by the pole-associated MT depolymerase (i.e. kinesin-13, or a combined effect of all the MT depolymerases; Figure 1B). In this case, we have to consider the force balance on this MT. Since the viscous drag for a single MT is negligible (Howard, 2001), we simply have to balance two motor forces (on the chromosome arm and at the pole: left- and right-hand sides of the equation, respectively), each characterized by its own force–velocity relation:

where Vchrspeckle is the velocity of the chrMT in the laboratory frame of reference (that could be observed as the velocity of a fluorescent tubulin ‘speckle'). Linear equation (2) is easily solved, giving us the chrMT velocity Vchrspeckle such that the motor forces acting on that MT are balanced. Then, the force pushing the spindle pole outward can be computed by substituting this Vchrspeckle into the negative depolymerase force:

To compute the total force on the chrMT population, we have to determine which motors are active and when. We do not explicitly model the mitotic regulatory network, but rather reduce it to the effective binary ‘on' and ‘off' ‘switches' of the molecular motors. We define the binary switch time-dependent parameters Pdep(t) and Pchr(t) (for depolymerase and chromokinesin) that are equal to 1 if the respective motor is engaged and generating force and either MT movement or growth/shortening occurs and to 0 otherwise. Using these switch parameters and equations (1) and (3), we obtain the total force on the entire population of chrMTs:

Equation (4) works as follows. Expression F=((1−Pdep(t))F2+Pdep(t)F3) represents the single MT force: if the depolymerase at the pole is active, Pdep=1, and F=F3, otherwise, Pdep=0, and F=F2. Nchr is the total number of chrMTs, so NchrF is the maximal possible total chrMT force. If the chromokinesin motor is active, Pchr=1, otherwise Pchr=0 and the whole force is zero. Finally, Achr, additional parameter, is the chromosome arm area, while the difference between S(t), the pole-to-pole distance, and D(t), the inter-sister chromatid distance, is the double pole-to-chromosome distance. The dimensionless geometric factor Achr/π(S−D)2 determines the fraction of the MTs impinging on the chromosome arm.

Similar arguments are used in the Supplementary information to calculate the total forces Fip, Faster and Fkt acting on the three other MT populations, and we summate them to obtain the total forces applied to the spindle pole, Fip+Fchr+Faster−Fkt, and to the chromosome, Fkt−Fchr−Fcohesion. (The results indicate that the kt forces are usually directed inward, while all others outward, hence the expression to the resulting net outward force on the pole. A chromosome is pulled outward by the kt force resisted by inward chromokinesin force and cohesion between the sister chromatids; the latter is approximated by a linear spring.) One additional twist of the model is that to calculate the ipMT forces, we calculate the dynamic overlap between the ipMTs at the spindle equator and assume that two active motors (kinesin-5 and -14) and passive crosslinkers generate additive sliding forces proportional to the overlap length. In the low Reynolds number environment of the cell, pole and chromosome separation velocities are determined by the balance of the total MT and cohesion forces and the effective viscous drag (Nédélec, 2002; Cytrynbaum ; Brust-Mascher ; Civelekoglu-Scholey ):

where μpole and μchr are the effective pole and chromosome drag coefficients, respectively. (The coefficients ½ are included to account for the fact that the pole and chromosome rates of movement away from the spindle equator are half that of the respective pole–pole and chromosome–chromosome separation velocities.) By solving equations (5) and (6), we recover the temporal dynamics of spindle poles and chromosomes.

Thus, our model has a modular character: we construct the ‘virtual spindle' (Figure 1A) from the ‘building blocks'—eight possible MT–motor configurations (Figure 1B)—using a combinatorial approach, in which the ‘building blocks' composition changes over time based on the 11 time-dependent binary switch parameters representing the activity timing of all major known force generators in the spindle (Figure 1). We allow each motor to switch its activity only once during mitosis: the ith motor (where i=1,…, 8) is active or not active all the time, or it is active (inactive) in the beginning and switches off (on) at a random time tswitch. In addition to these 8 switching times, each such virtual spindle is characterized by a random choice of 3 more switching times (regulating changes in MT numbers and chromosome cohesion), 16 stall forces and free unloaded velocity of 8 mitotic motors, and 8 other geometric, kinetic and mechanical parameters, 39 parameters in total. Each parameter's value can vary in certain wide range justified by available information (see the Supplementary information, Supplementary Table 1). This formulation allows us to define a 39-dimensional ‘model parameter space' and encode multiple models by points in this space corresponding to a large number of possible combinations of force–velocity relations, variations in the kinetics and timing of the regulatory switches.

The model is simple, based on drastic assumptions of perfectly symmetric, an effectively one-dimensional spindle, a homogenous distribution of motors, no dependence of motor affinity on the generated forces, only binary variations in the motor activities and MT numbers, additive multiple motor forces, etc. Even with this conceptual simplicity and a relatively small number of spindle structural elements, the total number of possible ways in which multiple mitotic motors of various characteristics can be integrated to build different mitotic spindles is astronomical. Thus, a straightforward scan of the entire model space is impossible, and we resort to the stochastic optimization process to identify the model parameters that obey experimental constraints.

Results of computational searches of the model space

To determine if the temporal regulation of mitotic motors' activity is essential, we searched for ‘virtual spindles' that could mimic the spindle elongation in WT embryos (Figure 1C) without any motor switching during mitosis, and we did not identify any good fit for the data under these conditions. We then searched for models that could explain the kinetics of spindle pole separation (Figure 1C) and chromosome motility (Supplementary Figure 4) in WT embryos, with motors allowed switching on or off only once. The number of models identified this way was very high (∼10 000). Figure 2A–D shows the time series for forces acting on four MT populations in the virtual spindle for just four out of the thousands of different models that agree with the WT data predicting almost exactly the time series for WT spindle elongation in Figure 1C. In the first two of these models (Figure 2A and B) that are very similar to each other, the spindle is governed by a large inward ktMT force being balanced by a large outward force generated by chrMTs. The ipMT and asMT forces in these two models are almost negligible. In the third, conceptually different, model (Figure 2C), all four MT populations participate in the force balance, but only up to the anaphase onset, after which all forces decrease by several orders of magnitude from hundreds to single picoNewton range. In the fourth, also very distinct, model (Figure 2D), the forces generated by all four MT populations are significant at certain time intervals (Supplementary information and Supplementary Figure 2).

Figure 2

Iterative elimination of models that do not fit experimental data. Initial analysis identified ∼10 000 models (belonging to ∼1500 groups) that agree with WT data. (A–D) Four examples of the forces on the four MT populations (asMT—blue; ipMT—green; chrMT—red and ktMT—cyan) in different models. (E) Crosses illustrate projections of the points in the parameter space corresponding to all identified models onto a two-dimensional manifold in the parameter space (explained in the Supplementary information). The models corresponding to (A–D) are marked with red circles. The iterative process of addition of experimental data to the search reduced the number of identified groups of models (F) and the area they occupy in the two-dimensional projection of model space (G) (same projection as that in (E)). The probability density distributions of the parameter values governing the timing of kinesin-14 activity and force for each iteration are shown in (H) with the same color coding as that in (F).

To analyze the multiple models statistically, we developed a quantitative distance measure that estimates how similar pairs of models are. This distance measure is based on differences in the magnitudes of forces and in the activity profiles for model pairs (Supplementary information). To illustrate the inter-model distances that are hard to visualize in high (39-) dimensional parameter space, we used this distance measure and projection onto a two-dimensional manifold in the model parameter space (Supplementary information; Figure 2E). In this figure, the metric two-dimensional distance between the model pairs reflects the multi-dimensional distance between full parameter sets characterizing the pairs of models. The positions of four sample models from Figure 2A–D in this projection, the quantitative similarity of models A and B, and significant differences among models B–D can be seen in Figure 2E.

The results of this first search lead us to the following conclusions: (i) there is a tremendous variety of plausible model parameters that can explain the WT behavior, whose combination can be complex and counterintuitive—there is no way to come up with a motor combination generating the force sequence shown in Figure 2D, for example. So, the WT data set, on its own, is not sufficient to discriminate between the multiple potential mechanisms of mitosis in this system. (ii) Many models are qualitatively different from each other, like those illustrated in Figure 2B–D, but some are very similar to each other, such as the pair shown in Figure 2A and B—the difference between those is in the small variation of a few parameters, and biologically, this pair describes basically the same molecular mechanism. This demonstrates the need for proper clustering of the models resulting from the computer search.

Clustering analysis (Supplementary information and Supplementary Figure 3) of the ∼10 000 models identified ∼1000 distinct model groups. Even this large number of possible distinct model types is an underestimate of the expected number of force integration scenarios. Convergence analysis (Supplementary information) showed that the search is still far from a complete exploration of total model space, and we estimate that there are ∼1500 possible model groups that can explain the WT data (Supplementary information).

To further constrain plausible models, we repeated the above search in an iterative manner, using at each iteration more of the experimental data on the dynamics of spindle length following inhibition of different motors (Figure 1C), in addition to the WT data. Specifically, we first used data for both WT spindles plus those with inhibited dynein, then WT, dynein- and kinesin-5-inhibited spindles, then we added data for kinesin-13-inhibited spindles, then for kinesin-14-null spindles and finally we incorporated the data for kinesin-5-inhibited/kinesin-14-null spindles (explained in Figure 1). The models were deemed successful if they predicted, with only a small error, the time series for spindle length dynamics both for WT, and for experimentally perturbed spindles. At the present, we do not have the data for the time-dependent inter-chromosomal distance for the perturbed spindles; as discussed in the Supplementary information, such data probably would not be of much use constraining plausible models. With each iteration, the search was conducted independently of previous iteration thereby not restricting the possible models to ones that were previously identified.

As expected, we saw that the addition of more experimental data decreased the number of identified model groups. Figure 2F shows the predicted number of groups from the convergence analysis as the experimental data accumulate. We saw that from ∼1500 model groups, when only WT data were used, there were only 6 different model groups ‘surviving' the scrutiny of the whole body of the experimental data after 5 iterations. Furthermore, with each additional iteration, the models that were in agreement with more experimental data occupied a lesser segment of the parameter space (shown in the two-dimensional projection (same as that in Figure 2E) in Figure 2G). This suggests that the additional experimental data can indeed constrain the number of viable spindle models. This trend can be seen further in the properties of individual components of the spindle machinery. For example, Figure 2H shows the probability density estimates of two parameters (switching time and force) that characterize the kinesin-14 motor, Ncd. The initially wide and almost unconstrained distributions of the possible parameter values obtained when only WT data are used are observed to narrow down significantly as the additional inhibition data are used. A similar trend was seen for parameters characterizing other force generators.

Comparison of conserved model features with experimental data

The final result of this iterative elimination process is a set of ∼1000 models that are clustered into six groups (Figure 3; the temporal changes in total force generated by each MT array and the ‘on–off' motor switching are shown for all these models in Figure 3A and for one selected model as an illustrative example in Figure 3B). Convergence analysis suggests that the search reached saturation and that all possible model types were identified (Supplementary information and Supplementary Figure 2). Interestingly, though the available data are insufficient to definitively narrow down the number of spindle model groups to one, thereby identifying the ultimate Drosophila spindle model, we discovered that several model properties were highly conserved among all ∼1000 models belonging to the six final groups, hinting that these properties are required for proper spindle design and mitotic progression in this system.

Figure 3

Cluster analysis of all identified models. (A) Results of the cluster analysis for models that fit all available experimental data. Dendrogram shows the hierarchical tree of all ∼1000 models. Each imaginary vertical line across the panel corresponds to a specific model fitting all available experimental data. Time series (represented on the y axis of each panel running from top (early) to bottom (late)) for the forces on the four MT populations and cohesion forces from prometaphase (t=0) till the end of anaphase B (t=278 s) follow immediately below the dendrogram. Six identified clusters within the tree are color-coded. The forces (in picoNewtons) are color coded according to the bar shown at the upper left corner (extreme red (blue) corresponds to 1500 pN (−1500 pN)). The forces are abbreviated as follows: Fas—total force on astral MTs; Fip—total force on inter-polar MTs; Fchr—total force on MTs reaching to the chromosome arms; Fkt—total force on the kinetochore fiber and Fcoh—forces induced by the cohesion complex that hold sister chromatids together. For the reference, the time series for the pole–pole distance are shown immediately below the force bar. Time series for 10 motor switches' activity follow immediately below the force time series. White and black correspond to active and inactive motors, respectively. The switches are: Pdep—pole depolymerizer; Pchr—chromokinesin; Pdyn—dynein; Pk5—kinesin-5 sliding motor; Pk14—kinesin-14; Pkt—combined kt motors; Pmt—MT plus end depolymerization activity at the kinetochore; Ppoly—MT plus end polymerization activity at the kinetochore; Pas—switch regulating the number of astral MTs; Povrlp—switch regulating the number of MTs at the overlap zone at the spindle equator. (B) The time series for the forces and switches predicted by a single selected model from group 1 as an illustrative example. The vertical bars (left) in the upper half of the panel represent the color-coded forces for each MT population, whereas the corresponding graphs (blue line, right) depict the change in magnitude of the corresponding force with time. The vertical bars in the lower half of the panel (left) represent the black–white switching ‘off' or ‘on' of the specific motors at the times indicated on the graphs (green line, right) for this particular model.

Thus, in all the models, the depolymerase(s) (Pdep in Figure 3), shortening the MTs at the pole thereby counteracting the outward thrust of other motors, must be active throughout mitosis to restrain spindle elongation, and must be switched off only at the onset of anaphase B to allow rapid pole–pole separation, in agreement with previous experimental data (Brust-Mascher and Scholey, 2002; Brust-Mascher ; Rogers ). The timing of kinesin-5 activity (Pk5 in Figure 3) is also highly conserved among all the models. This motor has to be switched on during prometaphase to exert outward forces on ipMTs until the end of anaphase B and it is the main force generator sliding the poles apart, again in agreement with experimental observations (Sharp ). Dynein (Pdyn in Figure 3) switches off uniformly before the end of metaphase and does not contribute to the outward force in the end of mitosis—this model prediction is surprising, since previous work suggested that cortical dynein contributes to spindle pole separation during late anaphase B (Sharp ). This (possibly transient) downregulation of dynein activity would be hard to predict using intuition without the system-level search, and it needs to be tested in the future.

Also, the forces exerted on asMTs and ipMTs are conserved: outward asMT forces of hundreds of picoNewtons pull the poles apart before metaphase and switch off afterward, and then ipMT sliding forces of hundreds of picoNewtons take over and push the poles apart during late prometaphase and metaphase and drastically decrease at the anaphase onset. All six model groups predict that the spindle is largely balanced by the outward ipMT forces, assisted at early stage by asMT forces and by inward ktMT forces. The magnitude of the forces and the timing of their action are model predictions that can be tested in the future. Other conserved mechanical features are discussed in the Supplementary information.

In addition to illuminating properties of the entire spindle, the modeling suggests that specific biophysical properties of the participating molecular motors are conserved. The plus-end-directed kinesin-5 and minus-end-directed kinesin-14 motors are proposed to act antagonistically on the antiparallel overlap zone of the ipMTs (Figure 1). Our search predicts that to reproduce the experimental results, kinesin-5 should be strong (great stall force) and slow (small unloaded velocity), while kinesin-14 should be weak (small stall force) and fast (great unloaded velocity). Figure 4 shows the predicted force–velocity curves for these two motors. Interestingly, this prediction was recently supported by in vitro biochemical studies: the experimentally measured unloaded velocities of the two motors are shown in the inset of Figure 4; indirect measurements reported in Tao are consistent with the notion that the stall force of kinesin-5 is greater than that of kinesin-14. Note, that the predicted motors' force–velocity relations (Figure 4) illustrate the robustness of the models to a few-fold parameter fluctuations, but not to order of magnitude changes. Qualitative feature of intersection of the force–velocity curves of the opposing motors is conserved as well. Predictions of other motors' force–velocity properties can be gleaned from the Supplementary information, Supplementary Table 3 and Supplementary Figure 6.

Figure 4

Constraints on the motors' mechanical properties. The predicted force–velocity relations of kinesin-5 (blue) and kinesin-14 (green). The main panel shows the force–velocity relation predicted in 50 resulting models randomly chosen from ∼1000 models that fit all available experimental data. The inset shows free load velocity measurements reported in Tao .

Open questions: differences between the six identified model groups

While many of the model features are conserved among all identified models, the six identified model groups have interesting biological differences among them summarized in Table I. Group 1 is unique in its difference from the other five groups in the following respect: outward forces resulting from chromosome arms, due either to chromokinesins or to MT polymerization, are downregulated during anaphase B. To maintain the balance of forces, the inward kt forces in this group are much smaller in magnitude than in the rest of the groups. This smallness is due to ktMT polymerization counteracting the effect of the inward-thrusting motors. Therefore, in this group of models, the balance of small, picoNewton-level forces is characteristic for the anaphase B spindle, unlike hundreds of picoNewton forces at the end of mitosis predicted by all other groups. The largeness of the parameter space corresponding to this group (Figure 3) could also mean that this is the most robust strategy of the spindle design, though at the present we cannot support or reject this statement. Below, we describe in detail the sequence of molecular events predicted by this group.

Table 1

Differences among six identified model groups

Group number	1	2	3	4	5	6
Strong cortical forces during prometaphase	−	+	−	−	−	−
Chromosome arm forces—active during prometaphase	+	+	+	−	+	+
Chromosome arm forces—active during anaphase B	−	+	+	+	+	+
ktMTs are in polymerization state during prometaphase	−	−	−	+	+	−
ktMTs are in depolymerization state during metaphase	−	−	+	−	−	−
ktMTs are in polymerization state during anaphase B	+	−	−	−	−	−

Group 2 is unique among the six groups since it is the only one that has a strong component of cortical forces that contribute to prometaphase elongation. Mechanistically, this is achieved by having a high number of asMT (∼200) with the cortical dynein active and pulling on the asMTs comparing to a low number of asMT (∼20) in the rest of the groups. Although in groups 2–6, forces on the chromosome arms are generated during anaphase B, in groups 2–3 and 5–6, these forces' contribution starts early in prometaphase, whereas group 4 is unique in that the chromosome arm forces are activated only during anaphase B. The rest of the differences among the remaining groups all result from differences in ktMT dynamics. Group 3 is unique since it is the only group where ktMT switch into depolymerization mode in metaphase. Groups 4 and 5 have ktMTs in polymerization state in early prometaphase, while in group 6, ktMTs are inactive until they start depolymerizing during anaphase B. Mechanistically, these differences in ktMT activity are achieved by differences in the switching times regulating activities of several participating motors and MT dynamics.

Predicted sequence of molecular events guiding spindle elongation

We chose one of the models from the largest model group 1 to describe the predicted sequence of molecular events guiding the spindle elongation (Figure 5). In this model, at the beginning of prometaphase, cortical dynein pulling forces supported by pushing forces from the chromosome arms are balanced by inward forces resulting from MT minus end depolymerization at the poles. Then, recruitment of kinesin-5 to the antiparallel ipMTs increases the net outward force, while dynein is switched off followed by the switching on of MT depolymerases at the kts. All these pre-metaphase events result in relatively small force balance changes, because strong antagonism between kinesin-5 and kinesin-13 on ipMTs largely determines the total force. Upregulation of kinesin-14 then reduces the net ipMT outward force to produce the metaphase steady state and the subsequent downregulation of the chromosome arm forces and degradation of cohesins substantially reduce the ktMT tension and mark the transition to anaphase A. Finally, anaphase B is the result of the switching off of MT depolymerization at the poles. Figure 5A shows the predicted time series for the spindle length, faithfully imitating the WT data. Figure 5B shows the forces acting on all four MT populations and generating the predicted spindle elongation, while Figure 5C illustrates the timing of switching of the mitotic motors responsible for generation of these forces. Prediction of this, likely most plausible, mechanistic scenario of the Drosophila embryo mitotic spindle elongation is arguably the most valuable modeling result, and would have been impossible to obtain without such a massive computer search.

Figure 5

Example of an identified model. The identified forces and main switches in one of the identified models (from group 1). (A) The pole–pole distance resulting from the net force acting on the spindle poles. (B) The MT populations' forces: asMT (blue), chrMT (red), ipMT (green) and ktMT (black). (C) The predicted sets of force generators acting in the spindle at sequential mitotic stages; roman numerals correspond to the stages in (A, B). (i) At the beginning of prometaphase, dynein pulling forces supported by pushing forces from the chromosome arms are balanced by inward forces resulting from MT minus end depolymerization at the poles. (ii) Recruitment of kinesin-5 to the antiparallel ipMT tips increases the net outward force. Then, downregulation of dynein (iii) followed by upregulation of kinesin-14 (iv) antagonizes kinesin-5, reduces the ipMT outward force (B) and decreases the total outward force resulting in the metaphase steady state. Upregulation of kt motors (v) changes the composition of the balance of forces (B) but has very small net effect on the total force (A). These changes followed by the degradation of cohesion (vi) are the only necessary transitions required for progression through anaphase A. Finally, anaphase B is the result of downregulation of chr-arm forces immediately followed by downregulation of depolymerization at the poles (vii).

Testable predictions generated by modeling

Variability of the multiple resulting mechanisms suggests that a few alternative activity profiles of mitotic force generators are possible and that further experimental work is needed to narrow down the number of groups and to identify the ultimate mechanism of spindle pole separation in Drosophila embryo mitosis. In other systems, it has been shown that chromokinesins are degraded in an APC-dependent manner during the metaphase-to-anaphase transition (Funabiki and Murray, 2000) and that this degradation is essential for entry into anaphase. It is yet to be determined whether this is the case in the Drosophila embryo system. The possibility that forces are generated by chrMT polymerization even after the chromokinesins are degraded further complicates the issue, so biophysical measurements of forces on the chromosome arms, as well as a chromokinesin-GFP-tagging experiment will be required to discriminate among models in group 1, group 4 and all other groups.

Time lapse movies of spindles with GFP-labeled tubulin (our unpublished data) together with immunofluorescence microscopy (Sharp ) suggest that the density of asMT is higher during later stages of mitosis than in early stages. This evidence does not support group 2 where the number of asMT in prometaphase is predicted to be much higher than in the later stages. The caveat of this observation is that it is impossible to determine the exact number of asMTs that are actively engaged with the actin cortex by using GFP-tubulin studies alone. Reconstruction of the MT array using serial section electron microscopy (EM) at different stages of mitosis would be required to resolve this issue. Finally, determining the state of ktMTs is very challenging since it is very difficult to distinguish between the different MT populations, let alone their dynamic state, by using fluorescence microscopy. The use of EM could be very helpful in this respect as well: it will be possible soon to determine the MT dynamic state in a fixed EM image based on the MT structure. Systematic characterization of kt fibers during consecutive mitotic stages would provide information necessary to distinguish among groups 1, 3, 2–6, and 4–5. In principle, if and when the data described in these two paragraphs become available, Table I demonstrates that such data will be sufficient to choose a single model group from the current six possibilities.

The computational modeling described above also allows computer experiments, suggesting additional future experimental studies. Previous studies revealed that the double inhibition of both kinesin-5 and kinesin-14 fully rescues metaphase spindle assembly (Figure 1C, dashed blue line) (Sharp ), but in recent studies we observed that the inhibition of kinesin-5 in null mutants totally lacking kinesin-14 function sometimes caused a collapse of the spindle (Figure 1C, green line). While this apparent discrepancy requires further experimental scrutiny, its analysis using the model suggested that the observed differences in the extent of inhibition of kinesin-5 activity could be a critical factor. Specifically, we ran simulations of the models from group 1 that were modified so that the kinesin-5 motor was partially inhibited—its maximal force was factored by a parameter less than unity. The simulations predicted that the partial inhibition of kinesin-5 to a concentration of ∼90% would result in a phenotype similar to that seen in embryos containing a WT level of kinesin-5. Further simulations predicted that the effect of more severe partial inhibitions of kinesin-5 would produce shorter metaphase spindle than those seen in WT spindles but should have little effect on anaphase B spindle elongation (Figure 6). This prediction is in contrast with previous models for spindle length carried out in Drosophila S2 cells (Goshima ), which suggested that metaphase spindle length is insensitive to the levels of kinesin-5. In addition, recent unpublished experiments suggest that the kinesin-5-dependent force balance that maintains the prometaphase spindle requires inward forces generated by kinesin-14 together with another unidentified factor. Our model further predicts that this additional inward force could be produced by a kinesin-13 depolymerase acting on ipMTs at the spindle poles and future experiments will be directed at testing this prediction. These examples illustrate the utility of our modeling strategy in identifying key, experimentally testable predictions.

Figure 6

Model prediction. Simulation result of dynamics of spindle length over time after complete (black line) and partial (color lines) inhibitions of kinesin-5. Partial inhibition is from 80% (blue line) to 0% (red lines).

Discussion

To summarize, we performed computer searches in the vast model space and identified six strategies for the temporal and structural organization of mitotic motors and MTs within the spindle, which quantitatively explain our observations of spindle elongation kinetics in WT, mutant and inhibited Drosophila embryos. Each of these strategies is characterized by specific activity timing and mechanical properties of each motor, MT number and a few other parameters. Note that the discovered models could not be obtained simply by combining earlier explicit models of specific mitotic stages, for example the inactivation of dynein (on asMTs) prior to the end of metaphase has not been considered nor suggested so far. Importantly, a number of features are conserved for all predicted models, including the timing of the activity of dynein and a few kinesins, as well as the forces and velocities of crucial mitotic motors. The search also revealed that large inward and outward forces in the range of hundreds of picoNewtons are closely balanced, which hints at a general design principle of mitotic spindle mechanics. In addition, the search also uncovered areas of uncertainties, mostly regarding the role of forces generated on the chromosome arms by MT polymerization and chromokinesins, the exact forces at the kts, and MT numbers and lengths. We make suggestions about future experiments that could help to resolve these uncertainties.

The use of an unbiased systematic computational search revealed a plethora of models producing basically the same overall phenotype. Even the most restricted search identified hundreds of model variants, all predicting almost identical spindle dynamics. Polymorphism is well known to exist among protein molecules such as those composing parts of the mitotic machinery. Our analysis suggests that a parallel polymorphism exists in terms of accurate model representation, so that there is not necessarily ‘one true model', but rather a set of models, all slightly different in their characteristics but capable of generating the same overall phenotype. This ‘model polymorphism' makes our computational search approach a useful ‘hypothesis generator' that can identify key experiments that are needed to further elucidate the mechanism of spindle elongation.

As is true for any model, our models depend on a set of assumptions, and this introduces a number of caveats into our modeling strategy. For example, (i) it is plausible that other, either known (e.g. kinesin-6), or yet to be identified, force generators act during mitotic progression; (ii) motor forces may not be additive; (iii) force–velocity relations may not be linear and (iv) motors may switch on and off more than once during mitosis. In addition, we assumed that the microinjected antibody inhibitors (Figure 1) serve only to reduce the effective concentrations of the target motor, but other effects (e.g. the antibody-induced inhibition of the mechanochemical cycle) could contribute to the observed inhibition (Ingold ). Such problems will be addressed in future applications of our modeling methodology.

Furthermore, the large-scale search comes with a price. It requires using an approximation similar to a mean-field approximation in theoretical physics and other simplifying assumptions to reduce the computation time for each possible set of parameters. In reality, the Drosophila embryo spindle is highly stochastic (Cheerambathur ), but this feature was ignored in the search. Incorporating realistic stochasticity, originating from MT dynamic instability and relatively small (in a thermodynamic sense) number of motors, into the entire search is computationally prohibitive, so instead we focused on investigating how the ‘winning models' behave following incorporation of simplified stochasticity, such as white Gaussian noise and spread of individual MT lengths as a result of dynamic instability. The simulation results (not shown) demonstrated that incorporating the MT dynamic instability and additional stochasticity in motor concentration would not change the predicted pole–pole separation profile, supporting the deterministic approximation used in the large-scale search.

This study builds on and extends several previous spindle models that were developed for the Drosophila embryo (Brust-Mascher ; Cheerambathur ) and other systems (Nédélec, 2002; Burbank ). We extended the idea of model screening and parameter search first proposed for spindle models by Nédélec (2002) and improved it from random sampling to more efficient search that uses repeated stochastic optimization. The improved efficiency is crucial since, unlike in the work of Nédélec (2002), in which the goal was to identify regions of parameter space that produce qualitatively interesting behavior, in this study the goal was to produce good fit to experimental data, similar to Gardner .

Some predictions derived from all six groups of models appear to contradict available data (Sharp ; Brust-Mascher and Scholey, 2002). For example, the models predict that kinesin-14 is only activated after prometaphase, that ktMT flux is virtually absent during metaphase (Supplementary Figure 5) and that the partial inhibition of kinesin-5 produces instability in the metaphase steady states. These discrepancies indicate areas of uncertainty that merit further attention. It is possible, based on preliminary estimates, that changing assumptions (i) and/or (iii) would lead to a correct prediction of the timing of kinesin-14 activation, while changing assumption (iv) would lead to a correct value for ktMT flux. (In the present form, the model has ktMTs permanently attached to the kts and, in fact, agrees with recent observations of DeLuca ). We see these inconsistencies as supporting the credibility and usefulness of the modeling, since they identify topics where further work is needed to reconcile theory and experiment, which is healthy for the discovery process.

The ultimate goal of systems biology is to construct comprehensive quantitative models for cellular function. One possible strategy is to ‘build with a scaffold' (Ideker and Lauffenburger, 2003) by initially constructing coarse grain models and later building on the results of these models to construct more elaborate detailed models. The construction of the initial models can be based on high-throughput data or, as in our case, reverse engineering of simplified models. A few pioneering studies have used similar reverse engineering approaches to different biological systems (reviewed in Ma'ayan ), mainly applying them to cell regulatory networks (Sachs ). Our study applied a reverse engineering approach that uses global indirect quantitative data to perform a comprehensive computational search to identify the mechanical design of the spindle that can explain such data. Our strategy allows us to examine numerous possible parameter values and alternative mechanisms using coarse grain models and later refine the ‘promising' models to include additional components with more detailed models. The suggested framework can be easily adapted to mitotic spindles in other organisms and in vitro (that may be designed differently) and, in fact, to many other biomechanical systems for which sufficient quantitative data exist.

Materials and methods

In brief, the analysis, performed automatically by the computer algorithm, proceeds as follows: (i) model parameters are chosen randomly from the allowed ranges; (ii) at each computational step, total MT forces for the four MT populations are found using equations similar to equation (4), the spindle geometry is updated by numerically solving equations (5) and (6), then the forces are updated for the updated geometry, etc., for the time t≈280 s. The solutions for each such parameter choice predict time series for the spindle forces and length. (iii) We collected several experimental data sets for time-dependent changes in spindle length in WT and experimentally perturbed spindles, and we averaged, smoothed and aligned the data (Supplementary information and Supplementary Figure 1; Figure 1C). The predicted time series for the spindle length are automatically compared with the corresponding data (Figure 1C), and a ‘score' is assigned to the model, so that if certain mathematical difference between the predicted and measured length time series is great, then the model's score is poor; while if the difference is small, the score is good (quantitative details in Supplementary information). The model is ‘screened' so that if the score is poor (good), the model is discarded (selected). (iv) Another set of model parameters is chosen randomly, until thousands of ‘good' models are amassed. (v) A genetic algorithm then automatically modifies or ‘mutates' the successful ‘selected' models, and repeated stochastic optimization results in the evolution of increasingly adequate multiple models, each producing excellent fits to all the available experimental data. (vi) Finally, a cluster analysis determines which of the adequate models are qualitatively different from each other, and which differ just by slightly distinct parameter values, and identifies minimal sets of successful strategies for force integration in mitosis. We conducted this optimization search process until a convergence analysis showed that the search had reached saturation and all possible model groups had been found. All optimizations and simulations were performed on an 11-node Linux cluster (each node had a 2 × Opteron-246 processor) using a combination of Matlab and C codes. Full description of the details of the numerical analysis and computational search is provided in the Supplementary information.