An agent-based model contrasts opposite effects of dynamic and stable microtubules on cleavage furrow positioning

It has long been known that (but not how) all animal cells use their mitotic apparatus to assemble their cytokinetic ring in the right place to bisect the mitotic apparatus. We propose and test a mechanism for how the mitotic apparatus determines furrow position. Recent work has identified the small GTPase, Rho, as essential for communication from the mitotic apparatus to the cell surface during cytokinesis (for reviews see Piekny et al., 2005; Wadsworth, 2005; Bement et al., 2006). Active Rho promotes both F-actin assembly at, and recruitment of activated myosin II to, the cortex (Kimura et al., 1996; Watanabe et al., 1997; Kawano et al., 1999; Afshar et al., 2000; Kosako et al., 2000; Royou et al., 2002; Dean and Spudich, 2006; Foe and von Dassow, on p. 457 of this issue). During cytokinesis, three essential proteins regulate Rho activity: a guanine exchange factor (GEF) that promotes Rho activation (called ECT2 in mammals, Pebble in flies, and LET-21 in worms; Prokopenko et al., 1999; Tatsumoto et al., 1999; Kimura et al., 2000; Kim et al., 2005; Yuce et al., 2005), a GTPase activating protein (GAP; called MgcRacGAP in mammals, Tumbleweed/RacGAP50C in flies, and CYK-4 in worms; Afshar et al., 2000; Hirose et al., 2001; Mishima et al., 2002, 2004; Somers and Saint, 2003; Zavortink et al., 2005), and a kinesin-6 motor (KRP110 in urchins, MKLP1 in mammals, Pavarotti in flies, and ZEN-4 in worms; Adams et al., 1998; Raich et al., 1998; Matuliene and Kuriyama, 2002). Henceforth, we use MKLP1 (the name for the founding member of this group) to subsume the inconveniently different organism-specific names for this kinesin-6 motor.

During anaphase, cell cycle biochemistry induces a potentially global cortical stimulation, acting through the Rho pathway (Niiya et al., 2005; Kamijo et al. 2006) that, in the absence of microtubules (MTs), actually occurs globally, but which the mitotic apparatus normally restricts to the equator (Canman et al. 2000; Foe and von Dassow, 2008). Activation at anaphase onset causes MKLP1 and the RhoGAP to form a heterotetramer known as “centralspindlin,” which binds and appears to activate the RhoGEF and, when concentrated near the cortex, promotes Rho activation (Jantsch-Plunger et al., 2000; Mishima et al., 2002; Somers and Saint, 2003; Yuce et al., 2005; Zhao and Fang, 2005; Chalamalasetty et al., 2006; D'Avino et al., 2006; Kamijo et al., 2006; Nishimura and Yonemura, 2006). Because the mystery we wish to solve is what role the mitotic apparatus MTs play in focusing Rho activation to a cortical zone where the contractile ring forms, we focus on the interaction between MKLP1 and MTs because the MKLP1 is a plus-end directed MT motor, potentially able to move this Rho-activating complex along MTs to the furrow site.

Recent studies reveal that two furrow specifying signals emanate from the mitotic apparatus (Bringmann and Hyman, 2005). The first causes polar relaxation, seen as loss of myosin and F-actin from the cortex near mitotic asters, which are composed predominantly of dynamic MTs (Werner et al., 2007; Foe and von Dassow, 2008; Murthy and Wadsworth, 2008). The second elevates localized myosin activation and contractility where the plus ends of stable MTs (astral or interpolar) underlie the cortex (Canman et al., 2003; Shannon et al., 2005; Foe and von Dassow, 2008; Murthy and Wadsworth, 2008). In echinoderm zygotes, a subpopulation of nocodazole-resistant astral MTs forms during anaphase, with those that aim toward the future furrow growing longest. The elevation of active myosin at the cell equator to a high level correlates with when and where the tips of these nocodazole-resistant (i.e., stable) astral MTs reach the cortex (Foe and von Dassow, 2008; see also Inoue et al., 2004). In contrast, beginning earlier in anaphase and continuing during furrowing, myriad dynamic astral MTs contact the cortex everywhere without causing myosin activation. When drugs shrink the mitotic apparatus so that the only MTs reaching the cortex are dynamic astral MTs contacting the poles, the poles remain devoid of activated myosin, whereas the rest of the cortex undergoes a global low-level recruitment of active myosin (Foe and von Dassow, 2008; see also Werner et al., 2007). These observations support a hypothesis first advanced by Canman et al. (2003): MTs have, simultaneously, opposite effects, with stable MTs boosting, and dynamically unstable MTs suppressing cortical activation during anaphase/telophase.

Treating echinoderm zygotes or cultured mammalian cells with concentrated nocodazole as they enter anaphase led to myosin and F-actin enrichment over the entire cortex (Foe and von Dassow, 2008; Murthy and Wadsworth, 2008). Or, when nocodazole treatment slightly later selectively abolished dynamic MTs while leaving stable microtubules in place, myosin was recruited to a wider-than-normal equatorial zone and to well above the high level typical of a normal furrow (Foe and von Dassow, 2008). This linkage suggests that dynamic and stable MTs likely exert their opposite effects via the same pool of Rho/myosin regulator (i.e., centralspindlin).

To explain how furrow specification works in light of the above experimental observations, we propose the following mechanistic hypotheses. H1: centralspindlin component MKLP1 binds to and walks along all MTs, but rarely reaches the cortex on dynamically unstable MTs because they continually depolymerize before MKLP1s arrive. Dynamically unstable MTs, aiming everywhere, thus sequester centralspindlin away from the cortex, suppressing Rho activation globally (with greatest sequestration nearest the highest MT density). H2: some astral MTs, aimed chiefly toward the spindle midplane equator, stabilize during anaphase and thereby become effective rails along which MKLP1 delivers centralspindlin, which activates Rho specifically at the furrow zone.

We use agent-based simulations to test whether, and how robustly, H1 and H2 can account for the myosin activation zones observed in cleaving echinoderm zygotes. We model 33,000 MKLP1 motors as surrogates for the Rho activating cargos they can deliver. We do not represent explicitly the cortical accumulations of activated myosin imaged, e.g., in Foe and von Dassow (2008). We assess only how two asters containing a subset of stable MTs can produce zones of elevated MKLP1 concentrations at the cell surface, which we interpret as source zones of Rho activation. From zones where MKLP1s concentrate, diffusion will carry the molecules of the signal transduction cascade away as they interact, bind and release the cortex, and eventually activate myosin. Such diffusive broadening must occur, but how much depends on rate and binding constants we do not know. Cortical contraction of the furrow may offset this diffusive broadening by an amount we also do not know. Modeling only MKLP1 concentration zones introduces slippage between our model's predictions of furrow width and the empirical results in Foe and von Dassow (2008). We did not try to disguise that slippage by adding to our model a many-step signaling cascade involving GAPs and GEFs, and Rho and myosin II phosphorylation, then adjusting rate and binding constants to make model predictions match empirical myosin activation data.

Our simulations show that the proposed mechanism is physically plausible, but only if the MKLP1 motors have a particular (experimentally verifiable) qualitative trait, and we predict that this mechanism is robust to perturbations of motor numbers, motor speeds, motor processivity, and to cell size variation (as long as cell diameter remains above ∼3.5 μm).

Our model represents the run-up to first cleavage in a rigid-cortex 3D model of a spherical 80-μm Strongylocentrotus purpuratus zygote with two large centrosomes (whose diameters we measured from micrographs of late anaphase cells). Each centrosome has a fixed statistically isotropic set of randomly aimed MT “sockets.” It nucleates new MTs (determined by a mean nucleation rate parameter) at random times from randomly chosen sockets until its total number of MTs reaches a number we prescribe. When an MT depolymerizes completely, the centrosome eventually nucleates another MT to take its place from a randomly selected socket. There are so many sockets that each is usually empty. MTs, with minus ends stabilized by attachment to the centrosome, aim away from the centrosome's center. Each in silico MT comprises a sequence of short straight segments, linked by springs, while other “coil watch springs” exert torques at hinge points, acting to straighten any MT that bends. The plus end segment's tip can elongate at a speed determined by a polymerization on-rate and the nearby concentration of GTP-charged tubulin dimers. The parameter values and tubulin dimer concentrations we used let MTs elongate at a maximum speed of ∼0.45 μm/s when no barrier interferes and GTP-tubulin dimers are plentiful. Simulated cells contained enough tubulin dimer (if it all polymerized) to make 1,500 MTs, each long enough to reach from a centrosome to the cell surface).

Our individual agent-based model uses a few hundred thousand differential equations to keep track of how, via Newton's and Fick's laws, myriad individual cytoskeletal agents (e.g., MT “segments” and individual MKLP1 motors) move, continuously with time, through a spatial continuum of viscous cytoplasm in response to forces these parts exchange after collisions between them. The growth of MTs by polymerization, as well as random thermal forces agitating every part, cause these collisions. As described in Materials and methods, our model also keeps track of diffusion of, and biochemical reactions among, soluble factors that form the individual agents, e.g., the MTs and the MKLP1s motoring along them. We will describe this model, implemented in Java computer code, in more detail elsewhere. Here, we describe the rules governing the biochemistry and molecular interactions this model implements and then summarize how the outcomes change when we vary some MKLP1 parameters, which this study identified as having a strong influence on how well the proposed mechanism works.

Extrapolating from in vitro findings (Ligon et al., 2006), our computer simulation assumes Clip170 binds to soluble GTP-tubulin dimers, copolymerizes with them to elongate MT tips, then falls off the tubulin lattice after tubulin-GTP hydrolysis (which occurs here immediately as the MT tip elongates). We assume the catastrophe rate increases as the amount of tip-bound Clip170 falls. This is consistent with Komarova et al. (2002) but differs from the interpretation made therein that tip-bound CLIP-170 boosts rescue rates. Clip170 in our model is a surrogate for any MT rescue or anti-catastrophe factor that recruits to plus ends, then departs by falling off the MT lattice. When an in silico MT tip hits a physical object, we assume (in agreement with Janson et al., 2003) that the resulting compressive force on the MT tip slows its growth (the growth rate, determined by on-rate and local free GTP-tubulin concentration, is attenuated by a factor that shrinks exponentially as compressive force rises). Even with copolymerized Clip170 continuously falling off, the tip concentration of Clip170 stays high as an MT elongates rapidly, but diminishes whenever a compressive force (or anything else, e.g., a GTP-tubulin shortage) slows tip growth rate. The indirect consequence of an MT tip hitting an obstacle is therefore to boost the catastrophe rate by reducing the tip concentration of Clip170. That leads promptly, stochastically, to a catastrophe. Usually, after an MT hits the cortex but before a catastrophe shrinks it, the compressive collision force buckles the MT slightly. Then, the MT shrinks until and unless rescue gets it growing back toward the cortex.

MT tip transitions between polymerizing/depolymerizing, an MKLP1 motor binding/falling off an MT, or a centrosome “socket” nucleating a new MT, etc., all occur as stochastic events.

We do not here address why some MTs stabilize long before they reach the cortex or why the longest aim toward the future furrow. To explore how this empirically observed pattern of MT stability (Foe and von Dassow, 2008) influences MKLP1 patterning, we stabilize a random subset of those MTs that aim, accidentally, toward the future furrow site when centrosomes nucleate them. We stabilized ∼10% of all MTs in this way. In our model, stabilized MTs do stop growing after a catastrophe but never shrink, whereas unstable MTs do shrink. In every other respect, including catastrophe and rescue kinetics, stable and unstable MTs are identical in our model, so MKLP1s interact identically with the two kinds. In Figs. 1, 3, and 4, and in the online supplement videos , green MT segments have high Clip170 levels (and so polymerized recently); orange ones have lost CLIP170 (and polymerized at least 25 s ago). Thus, 25 s after reaching the cortex, stable MTs turn, and remain, orange.

The assumptions and parameter values governing MT dynamics we used constitute one way to generate asters of MTs that rapidly elongate to the cortex and then stop growing shortly after collision with the cortex, resulting in asters that resemble what we see in micrographs of fixed anaphase cells: asters with a subset of stabilized equatorially aimed MTs. In the model, MTs serve only as rails, most of which are transient, which MKLP1s can bind and move along. Alternate rules that produced the same MT kinematics would likely make no difference to MKLP1 traffic. One MT kinematics issue that does matter is how far back toward the centrosomes depolymerizing MTs shrink before being rescued; different retreat distances lead to quantitative (not qualitative) differences in the spatial MKLP1 concentration patterns. Our parameter values let dynamic MTs shrink, on average, ∼10 μm back toward the centrosome, with a standard deviation of ∼8 μm. The mean time MT tips spent in contact with the cortex was 4.8 s (±3.6 s). The mean time between rescue and catastrophe was 5 s (± 5.4 s), and the mean time between catastrophe and rescue was 4.8 s (±4.3 s). In all the ∼100 simulations on which we based our conclusions, the parameters controlling MT dynamics were the same.

Each simulation run starts with 33,000 randomly located MKLP1 motors (tiny white dots in the figures). Model MKLP1s do not interact with each other, so changing their number does not change what the typical MKLP1 does. Before MT nucleation starts, the MKLP1 motors diffuse for 10 s into a uniform distribution. An MT segment can bind nearby soluble-form MKLP1s; stochastic events convert MKLP1s from soluble form to individual agents bound to the MT, which then walk along the MT at a fixed speed: different speeds between 0.1 μm/s and 0.6 μm/s in different simulations, and 0.2 μm/s in most cases (only the 0.6 μm/s speed outruns MT polymerization). MKLP1s fall off MTs stochastically with a half-life that we varied from 2.5 to 40 s. MT segment depolymerization casts attached MKLP1s into the local soluble MKLP1 pool.

Each computer run recorded the location and state of every MKLP1 agent and all other parts at each of 1,650 snapshots 0.3 s apart, which our software renders into videos that can exhibit any subset of the parts, colored as we like, from any vantage point. Figs. 1, 3, and 4 show simulation results at selected time points as static pictures, which are frames from our online supplemental videos, which convey our model's dynamic behavior.

We expected H1 and H2 to work. To our surprise, when we assumed that MKLP1 motors just walked off MT tips, we could not find any reasonable combination of parameter values for which H1 and H2 caused a significant buildup of MKLP1s at the furrow zone. MKLP1 motors did arrive faster and in greater numbers at stable MT tips near the equator than at the tips of dynamically unstable MTs, but diffusion scattered them as soon as they walked off the MT tips (Fig. 2 H). The modest concentration of MKLP1 motors in an equatorial band in Fig. 2 H is not, in our judgment, sufficiently higher than cortical MKLP1 concentrations elsewhere to make a biologically robust mechanism that could withstand diffusive attenuation during the signal transduction cascade that leads to Rho activation. These initial failures highlight a physical challenge any cell, especially a small cell, must meet to use MT motors to set up a sharp concentration pattern: how can MT motors generate any cortical pattern at all when diffusion can undo that pattern in a fraction of the time it takes motors to travel from centrosome to cortex? Adding another hypothesis to our proposed mechanism meets this challenge as follows. H3: upon reaching an MT plus end, an MKLP1 motor doesn't walk off, but stalls at the tip.

Our simulations show that H1–H3 lead robustly to sharp cortical equatorial concentrations of MKLP1s, with equator-to-pole concentration ratios exceeding 10 for broad ranges of parameter values. Images of dividing HeLa cells show that MKLP1 and its RhoGAP binding partner do indeed accumulate at MT tips in blobs resembling the heads of strike-anywhere kitchen matches (Nishimura and Yonemura, 2006).

Nedelec and Surrey (2001) anticipated that large consequences could result from the apparently minor difference of motors remaining attached to, versus motors walking off, MT tips. Our assumption that MKLP1s stall at MT tips means that in silico, they stop at tips until they fall off stochastically, according to their processivity half-life, and identically on stable and unstable MTs. If some other mechanism made tip-stalled MKLP1s stay attached for longer, MKLP1 concentrations at the furrow site would rise dramatically.

Each single multithreaded simulation of 8 real minutes took ∼30 h on a four-processor AMD Athlon Linux computer (2.33 Ghz clock). The parameters governing this simulation are so numerous that, given the run time to test each set of parameter values, we could not make any systematic exploration of parameter space. Instead, we sought parameter values, consistent with measured values in the literature, that yielded MT kinematics consistent with our observations of MTs in echinoderm embryos. With those fixed, we then tested values for various parameters controlling the MKLP1 motors. We found that, if and only if we assume that MKLP1s stall at the MT tips, then varying parameter values around those we used in Fig. 1 for MKLP1 properties and MT kinematics robustly yields sequestration of MKLP1 on dynamic MTs. This greatly attenuates cortical MKLP1 concentration everywhere, except where MKLP1s move to the furrow zone along stable MTs.

We also held the MKLP1 parameters fixed at Fig. 1 values and varied parameters controlling MT numbers and kinematics by factors of ±2, and we found that the furrow zone buildup of MKLP1s persists robustly. Too many MTs or too many MKLP1s so densely filled the video frames that the observer could no longer see clearly what was going on, so, in light of the model's robustness to MT kinematics parameter values, we made each centrosome nucleate only about half as many MTs as real urchin centrosomes do (compare Fig. 1, C and D, with the single confocal section in Foe and von Dassow, 2008, Fig. 2, D–G), and we held the MKLP1 count to 33,000. Fewer MKLP1s and MTs saved computer time. With only ∼100 stable MTs (10% of all MTs) occupying the furrow zone, there are very few at each latitude. Reduced MT and MKLP1 numbers cause jagged radial density plots. In test runs (unpublished data) with many more or fewer MTs or MKLP1s, the spatial patterning of MKLP1 was very similar to simulation results shown here.

To facilitate quantitative comparison of the outcomes with different parameter values, Figs. 1–4,²³⁴ show radial density plots of MKLP1 concentrations. These graph the number of MKLP1 molecules per volume of cytoplasm as a function of latitude. In all panels of all figures, red curves show the concentration of only freely diffusing MKLP1s (not those attached to MTs) averaged throughout a wedge-shaped conical shell, intersecting the cortex at each latitude in the graph (there are 200 such shells and each is the volume in the cell from latitudes Ø to Ø + ∂Ø). Shells near the cell poles have small volumes and contain few MKLP1s, making their computed concentrations noisy, so we omitted them.

The gray and white curves show densities of all MKLP1s, both bound to MTs and diffusing, in the outer 1.5-μm rim of the above conical shells (i.e., in the shells, but within 1.5 μm of the cortex, where, presumably, the centralspindlin these motors are part of would be close enough to the cortex to cause activation of Rho, hence of myosin, there). The gray and white curves do not change much qualitatively when we vary the rim thickness between 0.5 μm and 2.5 μm. The jagged gray profiles show instantaneous, high-spatial resolution MKLP1 concentrations. The smoother white profiles are low-pass spatially filtered, 25-s fading memory versions of the gray profiles, to mimic, respectively, short-range diffusion of Rho regulators and the likelihood that after Rho activation, deactivation occurs slowly.

The scales for the red versus the gray and white graphs differ by a factor of ten: ten radial grid lines for a red curve correspond to one grid line for the gray and white curves. This scaling allows representation of the outcome that convection of MKLP1 along MTs, especially stable ones, can build up much higher densities near the cortex than diffusing MKLP1s ever rise to anywhere.

Fig. 1 F gathers together and straightens out the five radial density plots that wrap around the spherical cell in the five video frames of Fig. 1. The collage summarizes the time course of MKLP1 patterning during the 8 min our simulations span (in real purple urchin embryos at 12°C, furrowing commences ∼6 min after anaphase onset). Fig. 2 shows eight such collages, each equivalent to Fig. 1 F, for different MKLP1 binding half-lives, motor speeds, etc.

Fig. 1 shows five video frames from a simulation (Video 1) in which MKLP1 motors move at 0.2 μm/s, have an MT binding half-life of 20 s, and, therefore, have a 4-μm mean run. Fig. 1 C, at t = 340 s, shows a strong cortical MKLP1 accumulation near the equator, with accumulations of MKLP1s at the tips of stable MTs. We know that the 10% of stable and equatorially aimed MTs are responsible for the dramatic accumulation of MKLP1s near the cell equator because this is much weaker when MKLP1s walk off MT tips (Fig. 2 H) and does not occur at all absent stable MTs (Fig. 2 G). Occasional match head clusters of MKLP1s do show up on tips of unstable MTs, but only transiently. The video also makes clear that MT perdurance alone acts to concentrate MKLP1s.

Fig. 1 (D and E) shows the consequence of simulated nocodazole eliminating all but stable MTs and casting most MKLP1s adrift to diffuse. In silico nocodazole treatment clearly amplifies the equatorial concentration of MKLP1s (see Fig. 1 E at t = 480 s). In real cells, nocodazole increased both the concentration and the width of the band of active myosin at the cell equator (see Foe and von Dassow, 2008, Fig. 4, A–C). Our simulations show increased concentration but little widening of the furrow because, as explained above, we have not implemented diffusion of the various players in the Rho activation and myosin phosphorylation pathways (which would unavoidably result in furrow widening).

In Fig. 2, we juxtapose panel Fig. 1 F (relabeled as Fig. 2 D) with seven other similar collages of radial density plots, made at the same time points and including nocodazole treatment from seven other simulation runs in which we varied parameters controlling MKLP1 behavior. In all these simulations, astral MTs are sufficiently numerous that most MKLP1s bind to MTs. Also, because dynamically unstable MTs outnumber stable MTs, the unstable MTs sequester most MKLP1s away from the cortex. We now describe salient points from each simulation summarized in Fig. 2 (see Video 2).

Using an MKLP1 speed of 0.2 μm/s and MT binding half-lives of 2.5, 5, 10, and 20 s, we had mean run lengths of single motors along MTs of 0.5, 1, 2, and 4 μm, which are shown, respectively, in panels A, B, C, and D of Fig. 2. In Fig. 2 F, the motor speed is 0.1 μm/s and the MT binding half-life is 40 s, yielding the same 4-μm mean run length as in panel Fig. 2 D. Fig. 2 (D and F) exhibits quantitatively similar equatorial buildups of MKLP1s. This illustrates our general finding that the amplitude of equatorial buildups of MKLP1 at equilibrium is proportional to mean run length but roughly independent of either motor speed or MT binding half-life separately.

In Fig. 2 E, with MKLP1s moving at 0.6 μm/s and the same 20-s MT binding half-life as in Fig. 2 D, the mean run length is 12 μm, and this produces a threefold greater equatorial buildup than in Fig. 2 D. Here, the motor speed exceeds the 0.45-μm/s MT polymerization rate. This lets some MKLP1s reach the cortex via unstable MTs, but only transiently. MKLP1s outrunning MT polymerization yield a lower “signal-to-noise” ratio because of low peaks of cortical MKLP1 concentrations all over the cortex.

Unfortunately, motor speeds for KRP110 (the purple urchin homologue of MKLP1) have not been measured. Using gliding assays, Mishima et al. (2004) found speeds of 0.16 μm/s for the mammalian homologue of MKLP1 monomers and 0.4 μm/s for MKLP1 dimers, speeds that lead to strong, high signal/noise ratios in our simulations. That is, we have explored the furrow-specification mechanism we propose using a low estimate of the speed of MKLP1 motors; higher motor speeds make the mechanism work better.

Fig. 2 G shows what happens, namely nothing but global sponging up of MKLP1s, when there are no stable MTs. We note that Foe and von Dassow (2008) show that, independent of MTs, the cell cycle gates global myosin activity levels (see also Foe et al., 2000; Royou et al., 2002; Werner et al., 2007); during late anaphase/early telophase, Foe and von Dassow (2008) observed intermediate level myosin accumulation where neither stable nor dynamic MTs impinge on the cortex, possibly because of an affinity of RhoGEF itself for the plasma membrane (Chalamalasetty et al., 2006).

In Fig. 2 H, MKLP1s walk off MT tips rather than stall, as in all other panels. Only a feeble equatorial concentration builds up before nocodazole treatment. After attenuation and spreading by diffusion of factors in the signal transduction cascade, we think that feeble equatorial elevation in MKLP1 concentration represents failure, underscoring how important the accumulation of MKLP1 on MT tips discovered by Nishimura and Yonemura (2006) is.

In all simulations summarized in Fig. 2, we used a diffusivity of 0.5 μm²/s (equal to 5 × 10⁻⁸ cm²/s) for soluble MKLP1 motors (in centralspindlin complex form). Greater or lesser diffusivities caused lower or higher equatorial peaks of MKLP1 buildup.

The simulation results we describe, and some 100 others not depicted, cannot begin to explore this model's parameter space fully to find the limits on all the rate/binding constants within which the model succeeds. Our results do establish that, in silico, H1–H3 can succeed robustly. MKLP1s, which have realistic speeds and processivities, and stall at MT plus ends, rapidly accumulate at MT tips to much higher levels on stable ones than on dynamically unstable ones. If stable MTs aim chiefly at an equatorial ring, whereas the dynamically unstable ones point everywhere, this causes sequestration of most MKLP1s away from the cortex on the great mass of dynamically unstable MTs, whereas the stable MTs focus a ring-shaped elevation of cortical MKLP1.

We call this intracellular, pattern-formation, reaction–diffusion–convection mechanism a “push–pull mechanism.” The “push” component involves some MKLP1s binding to stable MTs and motoring to their tips; because the stable MTs that reach the cortex do so chiefly at the spindle midplane, this push causes what Rappaport (1986) called “equatorial stimulation.” The “pull” component, which we interpret as causing classical polar (actually global) relaxation (Wolpert, 1960), involves most MKLP1s binding to the great mass of dynamically unstable MTs but never getting to the cortex along them because the MTs melt first. Diffusion acts against both push and pull components to attenuate the pattern by scattering MKLP1s whenever they don't bind MTs. Absent MTs, diffusion wins, causing globally uniform cortical myosin activation (and no cleavage). Even though the push and pull components act best synergistically, each alone has some capacity to specify furrow location.

We do not know yet whether echinoderm blastomeres at the 16, 32, … 1,024 cell stage have the same kind of stable MT array as the one-cell zygote has. Regardless, the mechanism we propose for large echinoderm zygotes is of general interest only if it scales robustly to operate across the large range of sizes spanned by the cells of many organisms, e.g., as the zygote partitions its cytoplasm into thousands of progressively smaller cells. After 10 divisions (there is no exact 1,024-cell stage due to unequal divisions and asynchrony), purple urchin blastomeres have a mean diameter of ∼8 μm. Will H1–H3 work robustly as cells become this small or smaller?

The “Peclet number” (Pe) is a dimensionless ratio measuring whether convective transport dominates diffusive transport when both mechanisms move a substance over a distance (L; e.g., the radius of a cell). If a molecule, e.g., centralspindlin, has a molecular diffusivity D, then the time it takes to move L is ∼L²/D by diffusion, and ∼L/v by convection (i.e., motoring) at speed v. If the ratio of these times, Pe = Lv/D, exceeds 1, convection works faster than, and should dominate, diffusion, and we can expect an equatorial buildup of MKLP1 concentration. This Peclet number falls below 1 as L shrinks below D/v, which is 2.5 μm for our estimated 0.5 μm²/s molecular diffusivity of centralspindlin in cytoplasm, and a v of ∼0.2 μm/s MKLP1 speed.

The classical Peclet number derivation assumes that convection and diffusion occur simultaneously, all the time, as for a solute diffusing in a flowing stream. Here, the MKLP1s move either by diffusion or convection, but never both at once. With f equal to the fraction of time moving by convection, we must modify the Peclet number to Pe = fLv/[(1 − f)D]. From our simulations, ∼75% of MKLP1s are bound to MTs at any moment after t = 200 s. With f = 0.75, this modified Peclet number falls below 1 as L shrinks below (1 − f)D/(fv) = 1.65 μm, so our mechanism should work for cells bigger than 3.3 μm in diameter. If we use this Peclet number calculation to ask how large the diffusivity of the molecule being concentrated could become before the mechanism fails in a 5-μm diameter cell, the answer is 1.5 μm²/s. The caveat to these Peclet number considerations is that nothing in them reveals, as our simulations did reveal, the apparent importance of MKLP1s stalling at, rather than walking off, MT tips.

Changing cell size had multiple unforeseen consequences that the calculation of a Peclet number does not account for, so this verification by a simulation involves subtleties. Fig. 3 shows H1–H3 succeeding in an 8-μm diameter cell, but with the following adjustments to compensate for unforeseen consequences. To make the small cell do anything interesting, we had to boost the concentrations of tubulin and MKLP1s in 8-μm, relative to 80-μm, diameter model cells. The large cells have 1,000 times greater volume than the small ones. If we keep concentrations the same, the small cell will have only 33 MKLP1s, not enough to form a recognizable pattern. To fix this, we put 3,000 MKLP1s in the small cell. If the small cell polymerized all of its 1,000-fold less tubulin, each centrosome could nucleate only 15 MTs long enough to reach the cortex, and by the time each centrosome made four such MTs, the tubulin concentration would have fallen too low to support polymerization. We boosted the tubulin concentration in the small cells 100-fold relative to large cells (so the small cells have one tenth the tubulin that the large cells have). Note that the 1,000-cell embryo has 2,000 copies of each gene supplying the same total volume of cytoplasm the diploid zygote managed with two copies, so synthesis from zygotic genes should be able to supply the ∼100-fold boosts we used for the small cell's tubulin and MKLP1 concentrations. We adjusted the centrosomes' nucleation capacity so that each made ∼80 MTs in the small cells. With the higher density of MTs per volume, an MKLP1 that falls off one MT soon collides with another nearby MT to bind to. This resulted in almost all MKLP1s being bound to MTs almost all the time in the small cell. This made f so close to 1 that the mechanism worked too well, for the trivial reason that the MKLP1s spent almost no time diffusing. Although this may be what actually occurs in tiny cells, we lowered the binding affinity for diffusing MKLP1s to attach to MTs until the fraction of time the typical MKLP1 is bound to MTs fell to ∼0.75 (making H1–H3 less likely to work). To resolve spatial heterogeneities of soluble factors in small cells, we reduced the diameter of cytoplasmic domains in them to 1.2 μm (compared with 5.4 μm in the large cells). We also reduced the diameters of centrosomes in small cells because centrosomes of the size we see in zygotes hardly fit into the small cells.

A single example here of how concentrations might change with cell size demonstrates the furrow specification mechanism we propose can operate robustly across a wide range of cell sizes. All variations of that example we tried came out similarly. Our need to scale up the tubulin and MKLP1 concentrations as the in silico cell becomes smaller hints that whether the furrow-specification mechanism we propose will or will not work robustly in real cells across a size range from hundreds of micrometers down to 5 μm hinges partly on how cells use the same number of gene copies to populate hugely different cytoplasmic volumes with key cytoskeletal proteins, and on how the capacity of centrosomes to nucleate MTs varies with cell size.

If future empirical studies fortify our verdict that the interplay between stable and dynamic MTs positions and shapes the cytokinetic furrow, which this in silico proof shows is effective and robust to cell size changes and other quantitative perturbations, we will have simply replaced an old puzzle with new ones: what causes some MTs to stabilize and why do the stable ones grow longer and predominantly toward the spindle midplane equator? Much additional barking will still be necessary, but we will at least now be barking up the right tree.

Our model keeps track of spatio-temporal concentration fields of soluble factors, and how soluble proteins interact with each other and with formed cytoskeletal elements such as MTs. For example, GDP tubulin dimers, liberated when MTs depolymerize, increase the GDP tubulin concentration locally near depolymerizing MT tips. The resulting GDP tubulin heterogeneities drive diffusive flux of GDP tubulin. Meanwhile, a biochemical reaction converts GDP tubulin dimers to GTP-charged dimers, generating spatial gradients, and thus diffusive flux, of GTP dimers. An MT plus end tip elongates only by incorporating (thus depleting) nearby GTP-tubulin dimers. We use continuum mass action kinetic equations controlled by rate constants, and a continuum Fickian representation of diffusive flux, to model temporal changes in 3D concentration fields of GDP-tubulin, GTP-tubulin, MKLP1, and CLIP170. This uses less computer time than representing myriad individual molecules as they diffuse and react. We implemented our mean field representations of soluble factors using a kind of Lagrangian coordinate system, which we now describe.

In most applications of our model, organelles (such as centrosomes, nuclei, yolk particles, etc.) occupied some fraction, and possibly a large fraction, of a cell's volume. We filled all the rest of the cell's interior with small spherical “cytoplasmic domains,” mathematical control volumes resembling Lagrangian elements of classical continuum mechanics. Cytoplasmic domains are small material parcels, which move continuously, carrying their soluble contents with them, and they number among the autonomous individual agents in our model. They exclude each other and exclude organelles by exchanging forces to repel one another whenever they collide. All cytoplasmic domains have the same diameter, a parameter of the model we adjust empirically to be big enough to minimize the number of domains to minimize computer time but small enough that the simulation outcomes do not change appreciably when we make the cytoplasmic domains even smaller. In the simulations shown in Figs. 1, 2, and 4, cytoplasmic domains had a radius of 2.7 μm, and in Fig. 3, 0.6 μm. Each domain is a “reaction vessel” agent, which interacts with other agents in the simulation. For example, MT agents pierce them and exchange viscous drag forces with them proportional to their relative velocity. Each cytoplasmic domain agent can implement any number of mass-action equations to account for biochemical reactions among the arbitrary number of soluble factors in it. Moving cytoplasmic domains carry the soluble factors they contain with them; this models convection. Random thermal agitation and collisions with centrosomes, nuclei, the cortex, or MTs move them. Regardless of how cytoplasmic domains move, each soluble factor can diffuse between any pair of neighboring cytoplasmic domains with a flux density proportional to the concentration difference times a molecular diffusivity coefficient and inversely proportional to the distance between those two neighboring domains. To facilitate efficient calculation of these fluxes, each of the cytoplasmic domain agents keeps a list of pointers to neighboring domains. Domains that move past each other, past organelles, or past MTs exchange fluid viscous drag forces. This exchange of viscous drag forces, together with exchange of forces that prevent domains from encroaching into each others' space, make the ensemble of cytoplasmic domains and organelles act like an incompressible viscous fluid; the implementation of mass action kinetic equations in each domain, and diffusive flux equations between adjacent domains, is an agent-based method for solving reaction–diffusion–convection conservation equations in that viscous fluid.

There are many variants of particle system models (which Reeves [1983] introduced) in use in computer graphics special effects algorithms to animate motion of water, fire, etc. Animators concoct equations of motion for each of myriad particles that make them move in apparently realistic ways and then render each particle with colors and texture that depend on scalar properties (like temperature) they carry. Many variants of particle systems are also used in scientific fluid dynamics simulations, where equations of particle motions account for real physics. Our model, which, with its cytoplasmic domains, is yet another particle system variant wherein equations of particle motion use forces to represent viscous drag when domains slide past each other (or past other organelles) and to represent collision forces that prevent domains intruding into each other, gets the physics approximately right. Steele et al. (2004) provides a good introduction, with references, to particle systems using equations of motion similar to ours. Chemical reactions among, and diffusive exchange of, soluble factors the particles contain do not feed back to influence particle motions, so we can simply add on equations to account for them.

Our individual agent-based simulation model is more complicated than an approach by Nedelec, Surrey and Maggs (2001), who matched experiments they performed with a continuum model of motors interacting with astral MT arrays in vitro (not involving cytokinesis). Both model and experiments were confined between coverslip planes, sufficiently closer together than the MT lengths, so that asters were essentially 2D. They used partial differential equations to represent the temporal development of spatial patterns of the concentrations of motors (both bound to MTs and diffusing). Using ∼1,000 50-μm long MTs, motor diffusivity coefficients, speeds on MTs, binding affinities, etc., similar to values we use, their models' predictions agreed well with the patterns of motor density they measured. They showed that motors would concentrate strongly farthest from the aster's center. The Peclet number consideration we discuss here applies directly to the situation they analyzed and affirms their result. Our model's complexity allows us to add several details important to account for dynamic asters with a patterned mix of stable MTs and unstable MTs (whose ends alternate between catastrophe and rescue), and to resolve whether MKLP1s will accumulate not just far from aster centers but specifically at the spindle midplane very close to the cortex. As is apparent in Fig. 4, predictions of their much simpler model, of our Peclet number, and of our agent-based model all agree to the extent possible given structural differences between the MT arrays in the two models.

Using myriad cytoplasmic domains may seem to be a more complicated way of solving reaction–diffusion–convection equations numerically than using a fixed-in-space Eulerian mesh to discretize a standard partial differential equation model, but the added computational cost is small. Regardless of how we deal with soluble factors, our general agent-based simulation must have a collision detection scheme. Any part can move anywhere and collide with any other part. We must therefore detect the many collisions that actually do occur between cytoskeletal parts and organelles as polymerization elongates MTs, and thermal agitation forces move all the parts around. We devised a collision detection scheme, which squanders computer memory to make computation time scale approximately linearly with the number of parts involved. Because the number of MT segments and MKLP1 motors (in the tens of thousands) greatly exceeds the number of cytoplasmic domains we require (about 1,000), the extra time required to detect cytoplasmic domain collisions is an insignificant increment. Also, the number of differential equations (three per domain) needed to keep track of the position of each domain is a small fraction of all the ordinary differential equations in the simulation. This means that the computational infrastructure already in our model can compute the trajectories of all the cytoplasmic domains for a small additional computational cost.

By keeping track of the locations of myriad close-packed cytoplasmic domains, and the numbers of soluble factor molecules in each one, we have a way to determine the concentration of any soluble factor near any 3D point, x, in the cell; we take it to be the concentration in the cytoplasmic domain nearest to x. Thus, for example, when x is the current location of a polymerizing MT tip, we can compute the concentration of GTP-tubulin near x, and thus compute how fast that tip should elongate (until it hits some obstacle). Similarly, for each MT segment, we can compute the concentration of soluble MKLP1 nearby, and compute from that the probability (during each time step) that an MKLP1 will bind to that MT segment (becoming a newly created individual agent that will walk along the MT and depleting, abruptly, the MKLP1 concentration in the cytoplasmic domain it came from). When, later on, that MKLP1 falls off its MT, it will disappear as an individual agent and boost abruptly the concentration of soluble MKLP1 in the nearest cytoplasmic domain (Fig. 4). Binding constants (for soluble MKLP1s binding MTs) and MKLP1-MT dissociation rates, etc., quantify all such stochastic transitions between soluble, and individual-agent forms of all cytoskeletal parts. Other forward and reverse rate constants control, via deterministic (i.e., nonstochastic) mass-action kinetic ordinary differential equations, the biochemical reactions among soluble factors in each cytoplasmic domain. Each of these mass-action equations has terms representing rates of loss or gain by diffusive flux between adjacent domains.

These Lagrangian cytoplasmic domains open modeling possibilities, which we do not exploit here, but which future work requires. For simplicity, the spherical cortex of the model cell in this paper is rigid. In other future applications, we plan to use cell cortices capable of moving and contracting in response to concentrations of regulators carried by cytoplasmic domains and by motors on MTs. This moving cortex will comprise a mesh of many small, contractile, 2D tiles (agents) that are linked together. Our general-purpose collision detector will detect collisions between these tiles and other agents the cell contains, in particular collisions with cytoplasmic domains. When a tile of a contracting cortex tries to push a cytoplasmic domain inward, the domain will push back with an equal and opposite reaction force. The result will be that our agent-based model will let us solve, approximately, the notoriously difficult fluid dynamics “free-boundary” problem. When, in classical fluid dynamics, a deformable boundary encloses a fluid, the partial differential equations of motion for the boundary and the partial differential (Navier-Stokes) equations of motion for the fluid contained by the boundary are coupled intimately. The boundary conditions on the fluid (what the stresses are and where they are applied) depend upon the unknown location and deformation rate of the boundary. How the boundary moves is, likewise, dependent on the movement of the fluid it contains. The best-known successful scheme for implementing this coupling in a numerical algorithm for solving the equations of motion is the one pioneered by C. Peskin (see Mittal and Iaccarino, 2005) and is used in many biologically relevant fluid-dynamical simulations. The agent-based simulation scheme we are developing attempts to solve this moving deformable boundary problem in a different, computationally tractable way.

The simulations were performed on AMD Athlon computers with 2.33 Ghz clocks and a variety of slower computers with fewer processors, all using Linux and Java.

Video 1, from which the panels in Fig. 1 were taken, shows our 3D computer simulation of MKLP1s interacting with stable and unstable MTs to elevate MKLP1s in the future furrow zone of an 80-μm zygote. Video 2, a video version of Fig. 2, shows the continuous time evolution of MKLP1 radial density plots for different assumptions and different motor speeds, and MT binding half-lives. Video 3, from which the panels in Fig. 3 are taken, shows the proposed mechanism succeeding robustly in a small 8-μm-diameter cell. Video 4, from which the panels of Fig. 4 are taken, shows the concentrations of soluble MKLP1s changing through time in cytoplasmic domains. Videos 1 and 4 are from the same simulation. Videos 5–10 show the individual simulations whose radial density plots are collected into Fig. 2 and Video 2. Video 5 corresponds to Fig. 2 A. Video 6 corresponds to Fig. 2 B. Video 1 corresponds to Fig. 2 D. Video 7 corresponds to Fig. 2 E. Video 8 corresponds to Fig. 2 F. Video 9 corresponds to Fig. 2 G. Video 10 corresponds to Fig. 2 H. The videos show only one half of the entire cell to make the compressed video files fit Journal of Cell Biology size limitations.

We thank Yvonne Beckham for help in manuscript preparation and George von Dassow, Bill Bement, and anonymous reviewers for suggesting improvements.

We thank The Seaver Institute, whose support enabled us to begin this work (grant “Using 3-D Microscopic Imaging To Study Mechanisms Of Cell Division, Cell Shape Change, And Cell Movement”), and the National Institute of General Medical Sciences, whose support of the Center for Cell Dynamics (1 P50 GM066050) allowed us to bring it to fruition.