Computational model of dynein-dependent self-organization of microtubule asters

Dynamic polymers known as microtubules (MTs) provide lines of transport, communication and control in the cell and organize cell movements, in part by serving as tracks for molecular motors (Bray, 2001). The most prominent motor families are the cytoplasmic dyneins, which glide toward MT minus ends (Holzbaur and Vallee, 1994), and kinesin motors, which are mostly plus-end-directed (Goldstein and Philip, 1999). MTs are commonly organized into polar asters with minus ends gathered at the center and plus ends extending outward (Kellogg et al., 1994). The free plus ends are dynamic, switching between periods of growth and shortening. This behavior, called dynamic instability (Mitchison and Kirschner, 1984), plays an important role in the exploration of intracellular space (Holy and Leibler, 1994). During interphase, most cells use MT astral structures to establish the spatial organization of organelles as well as to organize transport (Lane and Allan, 1998). In mitosis, focal points of two MT asters serve as poles for the bi-polar mitotic spindle, which plays a crucial role in the separation of the chromosomes (Sharp et al., 2000). Asymmetric MT asters interact with the actin cytoskeleton in a complex and poorly understood way to guide cell migration (Waterman-Storer and Salmon, 1999).

MT aster formation is normally attributed to the capacity of centrosomes to nucleate and stabilize the MT minus ends (Schiebel, 2000). Molecular motors also play an important yet poorly understood role in the organization of MT arrays (Sharp et al., 2000). Remarkably, polar MT arrays can self-organize in the absence of centrosomes (McNiven and Porter, 1988; Maniotis and Schliwa, 1991; Verde et al., 1991). For example, in mitotic extracts, aggregation of MT minus ends is accomplished by large complexes consisting of multiple cytoplasmic dynein motors, the dynein-activator dynactin and the large protein NuMa (Verde et al., 1991). The model of aster formation suggested previously (Verde et al., 1991) is based on the ability of the multivalent minus-end-directed motor complexes to associate with a few MTs and to stay attached to MT minus ends. The model asserts that MT minus-end focusing is achieved by the simultaneous motor driven transport of each MT to the minus ends of the other MTs attached to the same motor complex. Indeed, in mitotic cell extracts, short MT seeds were observed to be transported by cytoplasmic dynein toward spindle poles (Heald et al., 1996). Multivalent motor complexes also form radial MT arrays in purified in vitro systems (Nedelec et al., 1997). These studies emphasized the role of the physical movement of MTs as opposed to the view of MTs as immobile tracks.

Other studies of self-organization in MT/motor systems are based on in vivo observation of MT aster formation in cytoplasmic fragments of melanophores (Rodionov and Borisy, 1997a; Rodionov and Borisy, 1997b; Vorobjev et al., 2001). In melanophore cells, thousands of pigment granules are associated with minus-end-directed cytoplasmic dynein motors and plus-end-directed kinesin-related motors (Tuma and Gelfand, 1999). The granules thus have a motor-mediated ability to aggregate quickly to the centrosome-organized cell center, or to disperse uniformly throughout the cell. In centrosome-free fragments of black tetra melanophores, dynein motors are activated upon exposure to adrenalin and, as a consequence, random MT networks with uniformly distributed pigment granules rapidly rearrange into a radial array with pigment granules and MT minus ends focused at the center and plus ends at the fragment boundary (Rodionov and Borisy, 1997a). The experiments with melanophore fragments suggest a different model of MT aster formation than the multivalent motor transport model because recent studies showed that MTs are not transported through the cytoplasm (Vorobjev et al., 2001). Therefore, melanophore fragments provide a simple and specialized experimental system for studying a form of MT/motor self-organization that relies on the traditional view of MTs as immobile tracks.

Although the centrosome normally dominates MT organization, the self-organization phenomenon is important because it elucidates the intrinsic properties of MT/motor systems, and because there exist polarized MT structures that are not associated with centrosomes (notably, in meiosis). Also, the MT/motor self-organization phenomenon is likely to contribute to the centrosome-governed organization. In this paper, we develop and analyze a computational model of the dynein-mediated phenomenon of pigment aggregation and MT aster self-organization in melanophore fragments, based on earlier observations and measurements (Rodionov and Borisy, 1997a; Rodionov and Borisy, 1997b; Vorobjev et al., 2001). These studies generated data that call for computational modeling. By simulating the self-organization process on a computer and comparing the results with the experiments, the modeling provides valuable insight. Recent modeling of MT/motor systems (Cytrynbaum et al., 2003; Joglekar and Hunt, 2002; Maly and Borisy, 2002; Nedelec, 2002; Tran et al., 2001) proved to be an indispensable tool complementary to experimental studies.

In fragments excised from melanophore cells, MTs lose their organization in space and orient randomly (Fig. 1A), and pigment granules disperse uniformly throughout the fragment. Upon treatment with adrenalin (Appendix 1), dynein motors associated with the granules are stimulated leading to pigment aggregation and MT re-organization into a radial array (Fig. 1B) such that MT minus ends are embedded in the granule aggregate and plus ends extend toward the fragment boundary. Formation of the MT aster and pigment aggregation always occur simultaneously. Both MT dynamics and motor activity are crucial for this process: taxol treatment stabilizing MTs or depletion of granules from the fragment prevent self-organization (Rodionov and Borisy, 1997a).

Cell fragments used in the experiments (Rodionov and Borisy, 1997a; Rodionov and Borisy, 1997b; Vorobjev et al., 2001) (this paper) have a characteristic size of ∼30-50 μm (Fig. 2). Within a few minutes of adrenalin treatment, multiple local pigment aggregates emerge throughout the fragment. During the next few minutes, the local aggregates merge into a single focus. At the same time, pre-existent MTs turn over, and newly nucleated and assembled MTs emerge organized into a polar aster. The diameter of the pigment aggregate is in the order of 10 μm.

After re-organization is complete, the total quantity of MT polymer increases roughly two-fold (Vorobjev et al., 2001). Most MTs (∼80%) have their minus ends embedded in the aggregate, and plus ends at the periphery of the fragment (Rodionov and Borisy, 1997b) (Fig. 3A). About 10% of the MTs have their minus ends embedded in the aggregate with their plus ends growing toward the fragment boundary at a rate of a few microns per minute. The plus ends of another 10% of MTs are stalled at the fragment boundary while their minus ends disassemble at a rate similar to the rate of assembly of the growing plus ends. Rare MTs (∼1%) treadmill. There is no physical transport of the MTs (Vorobjev et al., 2001): lateral and longitudinal displacements of all MTs (in those treadmilling, or with plus, minus, or both ends stabilized) are orders of magnitude less than the displacements of the plus and minus ends. Most likely, the absence of MT transport is due to MTs being cross-linked into the actin cytoskeleton. The phenomenon of dynamic instability is not observed in this assay (Rodionov and Borisy, 1997b).

Consistent with the 80-10-10-1% partitioning of the MT population states as described earlier is the observation that after nucleation, a growing MT plus end reaches the fragment boundary in ∼5 minutes (a fragment radius of ∼20 μm divided by a growth rate of ∼4 μm/minute) (Rodionov and Borisy, 1997b). Minus ends stay in the aggregate at the center for an average of 40 minutes (hence, the ratio 8/1 = 40/5 minutes). Released minus ends reach the boundary in ∼5 minutes, because shortening occurs at a rate close to the growth rate of plus ends. In rare cases, a minus end is released from the aggregate before the corresponding plus end reaches the boundary explaining the ∼1% treadmilling MTs.

Pigment granules glide toward the MT minus ends at a rate of a few microns per second. The rate of movement is consistent with measurements of movements generated by dynein motors (Ashkin et al., 1990; Yamada et al., 1998). The granules detach from the MTs or pause frequently (with a rate of ∼1/second), and then associate with a MT fiber and glide again (Gross et al., 2002) (V. Rodionov, unpublished). Observations show that Brownian motion of granules is very slow, probably because their diameter is comparable with the mesh size of the cytoskeleton. When the granules are depleted partially, individual granules dislocate by a few microns over a few minutes (Rodionov et al., 1998). This indicates that the effective granule diffusion coefficient is ∼ 0.01 μm²/seconds. On the time-scale of seconds, a granule can move over ∼0.1 μm by thermal diffusion, which is important (see below). On the time-scale of minutes, however, this movement (microns) is negligible in comparison with MT/motor-mediated movements (tens of microns).

Previous studies (Rodionov and Borisy, 1997a; Rodionov and Borisy, 1997b; Vorobjev et al., 2001) suggest that some factors associated with the pigment granules have the ability both to nucleate MTs and to stabilize (decrease the net disassembly rate) their minus ends. The following scenario can explain the observed sequence of events after dynein stimulation (Fig. 2). First, the granules glide toward pre-existing MT minus ends, or local clusters of minus ends, thereby creating multiple local pigment aggregates. The local aggregates enhance nucleation of MTs, some of which grow through other local aggregates and establish immobile tracks for transport of granules from the latter aggregates to the former ones. This causes merging of the aggregates reinforced by more MT nucleation at the larger aggregates until eventually all pigment granules merge into a single aggregate. Simultaneously, all pre-existing MTs turn over and are replaced by MTs nucleated at the aggregate. Thus, the positive feedback loop based on dynein-mediated minus-end-directed transport and the ability of granules to nucleate and stabilize MTs seems to be sufficient for self-organization. Note that motor transport, MT kinetics and fragment geometry, but not MT/motor force generation, play a role in self-organization.

The described hypothetical MT/motor dynamics and interactions are simple, but they lead to complex behavior. In this study, we derive, analyze and simulate the mathematical model of the self-organization phenomena. Using this model, we address several issues, quantitative in nature, in an effort to validate the qualitative hypothesis. The model allows us to answer, within a precise quantitative framework, the following questions:

1. What is the minimal system that has the observed behavior? More specifically, is it possible that the nucleating activity of granule associated factors alone is sufficient for self-organization? Or is it possible that minus-end stabilization alone without additional nucleation can be responsible for pattern formation? If both nucleation and minus-end stabilization are necessary, are they equally important (in some quantitative sense)?

2. What is the nature of the breaking of symmetry? How does aggregation start from the homogeneous MT/granule distribution?

3. What determines the time-scale (minutes) of self-organization? Is there a dependence of pattern formation time on fragment size?

4. How does fragment shape affect pattern formation?

5. How do MT nucleation and stabilization rates depend on granule concentration? What determines the number and length distribution of stable and treadmilling MTs?

In this section we describe a computational model based directly on the qualitative hypothesis described above, the precise mathematical details of which are given in the Appendices. The following model assumptions stem from observations, some of which are described above.

1. The MTs are straight and do not move. (As seen in Fig. 1, the MTs do bend but not significantly.)

2. The MT plus ends grow at a constant rate v_p.

3. Any plus end reaching the fragment boundary stops there. When the corresponding minus end reaches the boundary, the MT disappears.

4. The minus ends shorten at a rate v(g) which is a decreasing function of the local granule density, g. In the absence of the granules, the MTs treadmill: v(g=0)=v_p.

5. The MT nucleation rate, n(g), is an increasing function of the local granule density. The MTs are nucleated isotropically, at random angles.

6. Each granule is either not associated with the MTs and immobile, or is associated with a single MT and glides toward the corresponding minus ends with constant speed, v_g. Gliding granules detach frequently at a constant rate, k_off. Immobile granules attach to a MT with a rate proportional to the local concentration of MTs of a given orientation. The proportionality coefficient is k_on.

Note that the MTs are very likely to be crosslinked into the actin meshwork. Therefore, the force generated by motors moves granules, while the reactive force applies to the whole actin/MT ensemble, the movement of which can be neglected.

Fig. 4 illustrates possible granule density dependencies of the net minus-end disassembly rate v(g) and of the nucleation rate n(g). For modeling disassembly, we assume that minus ends are in one of two states: capped and disassembling. If the capping dynamics are sufficiently rapid, the nature of disassembly can be summarized by function v(g), which monotonically decreases to zero as the granule density increases and which depends on three parameters. A derivation of this functional form is given in the Apppendix but here, for simplicity, we present the range of possibilities in qualitative terms. We investigate two possibilities: (1) a linear or non-cooperative dependence; and (2) a threshold or cooperative behavior where minus-end shortening is not density dependent at low densities. The nucleation rate includes both a spontaneous (granule independent) rate and a density dependent rate. We consider three possible types of density dependence for the nucleation rate: (i) quasi-linear increase with growing density; (ii) quasi-linear increase at small density and saturation at greater density; and (iii) threshold behavior at small density and saturation at greater density. Simulations of the model help to choose which of the described possibilities are more likely.

We consider two mathematical implementations of the model. First, we consider a 1D model. Biologically, this is equivalent to carrying out the experiments in a narrow elongated fragment in which the average MT length is much greater than the fragment width (see Experimental verification of the model predictions). This case is much simpler to treat mathematically than the realistic 2D situation and it gives insight useful in understanding more complex cases. To address the issue of the role of the shape and area of the fragment in pattern formation, we explore a realistic 2D model. The models are explained below and described in detail in Appendices 2 and 5, respectively. The model parameters and variables are listed in Tables 1 and 2; some of them are known from experimental measurements, others can be derived indirectly from the observations, and the remaining few are determined from comparisons between the theoretical and experimental results (Appendix 3).

Fig. 5 illustrates the 1D implementation of the model. In this case, the boundary of the fragment consists of two ends of the fragment at x= ±L. The MTs can be separated into two dynamic sub-populations characterized by opposite orientations. Each population is described by the dynamic densities of the plus ends [p_r,l(x)] and minus ends [m_r,l(x)]. The index r (l) stands for the right- (left)-oriented MTs, which have their minus ends to the right (left) of their plus ends. This notation is chosen because the pigment granules slide to the right (left) on right- (left)-oriented MTs. Another important characteristic of the MTs is the polymer density N_r,l(x) defined as the number of MTs (expressed as a density) passing through the cross-section of the fragment at coordinate x (Fig. 5). The granules can be described by three densities – those gliding to the right (g_r) and left (g_l) with speed v_g on the right- and left-oriented MTs, respectively, and the density of static granules (g_s) dissociated from the MTs (Fig. 5). The static granules associate with the right- (left-) oriented MTs with rates proportional to the local polymer densities of the respective fibers, k_onN_r,l(x). This model is deterministic and does not consider stochastic effects. The model equations are introduced in Appendix 2.

The computational structure of the 2D model is a combination of a discrete stochastic simulation of individual MT fibers and numerical solution of continuous deterministic equations for granule density.

A deterministic continuous description of the MTs on a 2D domain is not reasonable. First, mathematically, such a description presents a multi-dimensional problem that is difficult to formulate and time consuming to simulate. Second, the total number of MTs in a fragment (100-200) is large, but the number of MTs passing through any local region in the fragment is not so large. Thus, a discrete, stochastic representation of MTs is both more convenient and more appropriate. Therefore, we implement a discrete stochastic model: individual MTs are nucleated with random orientation and position, and their dynamics are tracked in a manner consistent with the 1D model. They appear, in time, according to a Poisson distribution and, in space, according to the local nucleation rate, n(g), determined by the pigment distribution. Plus ends grow at a constant rate v_p and are stabilized when they encounter the cell membrane. Minus ends shorten at a rate dependent on the local pigment concentration, v(g).

In Appendix 4, we demonstrate that the rapid movements and attachment/detachment kinetics of the pigment granules can be successfully approximated with a combination of diffusion and advection of the granule density. The diffusion arises from the random walk undertaken by pigment granules as they attach to and detach from MTs pointing in various directions and has nothing to do with Brownian motion, which is several orders of magnitude smaller and hence irrelevant. The local rate of advection depends on the presence of anisotropy in the local MT distribution – an isotropic MT distribution generates no advection whereas any bias causes movement in the direction of the bias.

We use the following scheme to calculate the velocity field associated with a MT. Each MT generates an effective velocity field in a rectangular domain of influence (Fig. 6). The corresponding velocities are minus-end-directed and decrease away from the MT according to the formula in Appendix 5. To compute the influence of MTs on the local granule advection, we derive an expression best described as scaled linear superposition (see Appendix 5 for details). This scaled linear superposition describes a sampling of the local MT network by pigment granules rapidly attaching to and detaching from the network allowing for each granule to be influenced by all MTs present.

The described extrapolation procedure is consistent with the 1D model in the sense that a few parallel and anti-parallel MTs placed close together generate the same velocity field as a similar arrangement would in the 1D model. Furthermore, if the local MT distribution is isotropic then the corresponding average advection is zero, and the granules disperse by an unbiased random walk, which also conforms to the microscopic observations.

In Appendix 4, we show that the granule density equilibrates rapidly to the current MT distribution. That is, at any given moment, the MT array is changing sufficiently slowly that pigment granules achieve a pseudo steady state distribution based on the current state of the MTs. In the simulations, the MT distribution generates the velocity field at each time step according to the described procedure and we adjust granule density instantly according to the described combination of diffusion and local advection at each time step (see Appendix 5 for 2D equations and numerical procedure). We update the MT distribution according to the described discrete stochastic algorithm using the current granule density. These calculations are repeated at each computational step.

We simulated the 1D and 2D models numerically. The dynamic behavior of the models can best be appreciated by viewing Movies 1-7. Figs 7, 8, 9 and 12 show frames from these movies. The model simulates a broad range of features of the self-organizing MT/motor system. In addition, theoretical analysis of the 1D model is given in Appendix 4.

Our explanation of the nature of the symmetry break at the onset of pattern formation is similar to the detailed analysis of the same 2D system in a disk-shaped fragment (Maly and Borisy, 2002). As shown below, in the fragment, where MTs are nucleated uniformly and isotropically, most microtubules are directed with their minus ends toward the interior of the fragment. The effective granule velocity field generated by the simulated MT ensemble is almost radially symmetric, directed toward the center and increasing away from the center, in complete agreement with the earlier analysis (Maly and Borisy, 2002).

The 1D model provides additional insight into the mechanism of symmetry breaking. The left-oriented MTs nucleate everywhere within the fragment with equal likelihood, grow in length at a constant rate and treadmill to the right. Thus, treadmilling shifts longer left-oriented MTs to the right. The same argument applies to the right-oriented MTs. Thus, even though MTs are nucleated uniformly throughout the fragment, the array that forms is not uniform but rather generates a bias in pigment movement. Indeed, near the right (left) edge of the fragment, almost all MTs are leading attached granules to the left (right), while at the center of the fragment, the numbers of oppositely oriented MTs are equal, and there is no bias (Fig. 5).

Fig. 7A shows analytical solutions of the 1D model equations (see Appendix 4) illustrating the MT and granule distributions at the onset of the self-organization. As one goes from the left (right) to the right (left) edge of the fragment, one finds a linear increase in the number of left- (right)-oriented MTs. Rigorous calculation (Appendix 4) shows that the effective velocity of the granules is an odd function of the distance from the center: increasing from the center toward the edges, zero at the center, and directed toward the center. Thus, within a few tens of seconds after stimulation of dynein, rapid transport of pigment granules leads to the formation of a loose aggregate shown in Fig. 7A.

Using computer simulations in both 1D and 2D, we experimented with a lower nucleation rate. It was found that aggregation and aster formation still took place, at roughly the same time-scale, but to a lesser extent: the aggregate was looser, and the minus ends were distributed widely in the fragment. A similar effect was observed when the ability of pigment granules to stabilize minus ends was attenuated. In the latter case, the time-scale of pattern formation also increased. These results conform with the observations reported previously (Rodionov and Borisy, 1997a).

Movies 1 and 4 show the simulated time course of self-organization in 1D and 2D fragments, respectively. In 1D, a loose aggregate at the center nucleates more MTs, the minus ends of which stay in the aggregate while the plus ends grow outward. This accelerates granule aggregation, which causes the nucleation of more MTs that remain anchored in the aggregate. This positive feedback loop leads to a tight granule aggregate at the center with most MT minus ends focused in the pigment aggregate (Fig. 7B, Appendix 4).

In 2D, the granules rapidly aggregate into a few local foci (Fig. 8B) and the local aggregates coalesce into the single group (Fig. 8C). At the same time, the re-organization of the MTs starts: MTs with minus ends outside the aggregate treadmill toward the boundary and disappear. Finally, the pigment aggregate tightens and the MT aster emerges (Fig. 8D). This simulated pathway reproduces faithfully the experimentally observed sequence of events (Figs 1, 2; Movie 6). Moreover, the simulations illustrate transient distributions of MTs that are hard to observe in the experiments. They demonstrate that indeed the initial local aggregates emerge near local MT minus-end clusters, and that mini-asters hypothesized in the qualitative model do not have time to form.

Numerical simulations also give an estimate for the time-scale of pattern formation that is close to that observed (∼10 minutes in both experimental and theoretical systems) in fragments a few tens of microns across. In Appendix 3, we demonstrate that the relevant time-scale is determined by the characteristic time of treadmilling across the fragment, ∼L/v_p ∼20 μm/5 μm/minute = 4 minutes. Simulations show that pattern formation takes place in roughly three of these time units. The observation that self-organization takes roughly three time units is indicative of the fact that the MT array must reassemble roughly three times before aggregation and aster formation are complete.

Movie 1 and Fig. 7 demonstrate that, in 1D, aggregation proceeds directly to a single granule aggregate at the center without the appearance of transient local aggregates. However, Movie 2 starts with an inhomogeneous initial MT array. This is intended to represent the random variation in MT densities in the fragment without explicitly introducing stochasticity to the model. Initially, the pigment forms three aggregates near the respective initial minus-end clusters. Subsequently, the three aggregates merge into two and finally into a single aggregate.

To explore the role of stochasticity in a more realistic setting, we simulated the aggregation pathways in the 2D domain at various nucleation rates. Fig. 9 illustrates that at a moderate nucleation rate (reflected in the low number of MTs), a few transient aggregates emerge near local clusters of minus ends. At nucleation rates an order of magnitude higher, much greater numbers of MTs are associated with a very regular distribution of minus ends. This leads to a smooth and solitary aggregate (Fig. 9B). These results suggest that local inhomogeneities in initial minus-end distribution in the fragment are responsible for local transient aggregates. These inhomogeneities are due to stochasticity of MT nucleation when the nucleation rate is such that 100-200 MTs are present in a fragment a few tens of microns in size. Varying the nucleation rate through a full order of magnitude had almost no influence on the time course of aggregation. This last fact is a theoretically predictable feature of the model.

Using computer experiments, we established that there should be no threshold in the density dependence of the minus-end disassembly rate (dashed curve in Fig. 4). The reason is that if the minus-end disassembly does not slow down when the granule density increases, then the positive feedback loop leading to self-organization is not triggered. The self-organization pathway is not sensitive to the presence or absence of either a threshold or saturation in the density dependence of the nucleation rate.

Our computational model is not appropriate to specify the quantitative dependencies of the rates of MT nucleation and minus-end disassembly on granule density. The reason is that the model neglects the effect of the steric repulsion of granules in the aggregate: in reality, there is a maximal granule density in the aggregate due to crowding, while in the model, the density is not limited from above.

This result means that 10% of MTs in the fragment are nucleated spontaneously, and another 90% are formed by granule-dependent nucleation. In the aggregated state, the percentages are similar; however, owing to minus-end stabilization, virtually all MTs in the aster are granule nucleated. The density dependence of the nucleation rate is sub-linear and mildly saturating. Before stimulation with adrenalin, the net rate of minus-end disassembly is roughly three times less than the rate of plus-end growth. After aggregation, the net rate of disassembly of the minus ends embedded in the aggregate is approximately 30 times less than the rate of plus-end growth. These quantitative predictions of the model can be tested.

We computed the total length of MT polymer and the partitioning of the MTs in the presence of a tight aggregate. The results are shown in Fig. 3A. The increase in total length of MTs is due to both greater number and greater length of individual MTs. The total length of MT polymer increases by a factor of three, compared with roughly a two-fold increase observed experimentally (Vorobjev et al., 2002). However, the experimental number probably underestimates the total because portions of MTs embedded in the aggregate are not visible. The computations predict the 80-10-10-1% partitioning of the stable – minus-end disassembling, plus-end assembling – treadmilling MTs, in close agreement with the observations reported previously (Rodionov and Borisy, 1997b).

As noted above, the time-scale of self-organization is determined by the ratio L/v_p. Therefore, the model predicts that the pattern formation time-scale is linearly proportional to the size of the fragment.

To elucidate the role of fragment shape in self-organization, we performed computer simulations in 2D domains with a characteristic bi-lobed shape (Fig. 10; Movie 7). The simulations suggest the following sequence of events. Initially, very few MTs pass through the corridor, so there is little communication between the lobes. Thus, self-organization proceeds in the lobes almost independently, according to the scenario described above for regularly shaped fragments. However, after two polar asters are organized in the adjacent lobes, there is an increased number of MTs transiently anchored by their minus ends in one of the granule aggregates extending through the corridor and passing through the other aggregate. These MTs establish tracks for granule transport, so that granule density in the corridor increases. This augments the nucleation of MTs with their minus ends in the corridor extending outward and thus passing through both granule aggregates. This accelerates directional granule traffic into the corridor further enhancing the formation of a polar MT aster in the corridor and depleting granules from the initial aggregates. This positive feedback loop leads to the final centered aggregation.

Thus there are three possible aggregate/aster formations, one stable and the other two only transiently stable. The two aggregates that form in the lobes are stable only as long as there is no communication (i.e. MTs) between the lobes. The time-scale of their `stability' depends on the expected time before a MT grows through the neck, predicting an inverse relationship between the thickness of the neck and the time required for the two aggregates to fuse in the middle of the neck. Furthermore, a neck which is sufficiently curved to prevent growth of MTs from one lobe into the other should indefinitely sustain two stable aggregates in the lobes. Note that our results complement previous findings (Nedelec, 2002), who demonstrated theoretically that MTs and minus-end-directed motors cannot support stable multi-aster structures.

It is difficult to obtain two fragments of different size and similar shape, and so we tested Prediction 1 semi-quantitatively. We observed pigment aggregation in 20 fragments of different size excised from nine cells (for a sample cell with two fragments see Fig. 2; Movie 5). In the sample cell, the large and small fragments were roughly 70 and 40 μm wide and had areas 3800 and 1300 μm², respectively. At each time step, we measured the state of aggregation by taking phase contrast images of the fragments and calculating the second moment of the pigment distribution about the centroid in both the x- and y-direction and taking the square root of their product. (We imported the images into Matlab™ and used this software to measure and quantify the distribution of pixels.) This measure of the granule dispersal has the dimension of length. Before adrenalin treatment, when the granules were dispersed homogeneously, the measure of granule dispersal defined in this way could be used to quantify the size of the fragment. A circle with radius twice this value should contain roughly 95% of all pigment granules so that four times this value is roughly the diameter of the fragment. To find the associated time constant of aggregation for each fragment, a decreasing exponential function was fit to the temporal sequence of second moments. After a slower initial stage of aggregation, decreasing exponential functions provided excellent fits for the linear granule dispersals as functions of time. Fig. 11 shows the observed relationship between space and time constants for the 20 fragments. A straight-line fit to the data (least squares sense) is superimposed with dashed lines representing one standard deviation above and below. Note that we do not assume that the line goes through the origin, a feature that is required to satisfy the scaling prediction. Rather, this feature, which is evident in the figure, simply emerges in the least squares fit. The slope and intercept of the fit were found to be 26.1±7.1 second/μm and 10.5±54.8 seconds, respectively.

Note that the scaling law (T=L/v_p) offers a prediction for v_p, in particular, v_p=0.038 μm/second=2.3 μm/minute, which is reasonably close to the directly observed experimental value (4-8 μm/minute) (Rodionov and Borisy, 1997b; Vorobjev et al., 2002). These observations and the image analysis confirm that the pattern formation time-scale is proportional to fragment size.

We also made two long narrow fragments, observed the pigment aggregation process in them and compared the observations to the 1D model analysis. Owing to technical problems, it is hard to make such fragments, therefore we were not able to produce fragments of various length and test the model quantitatively. The results are shown in Fig. 12. During the aggregation process, local aggregates form, as seen in this figure. Finally, a single tight aggregate forms as the system approaches steady state. The comparison with the 1D model shows qualitative agreement of the experiment and theory.

To explore the relationship between pattern formation and domain shape, we prepared a number of fragments with a characteristic bi-lobed shape and observed self-organization in such fragments (Fig. 10; Movie 6). In complete agreement with the theory, we observed that pigment aggregation proceeded separately in each of the lobes during the first few minutes after adrenalin treatment, much like the aggregation process seen in regularly shaped fragments. A large aggregate developed in each lobe, somewhat closer to the corridor between the lobes than to the centers of the lobes. In the next few minutes, the aggregates converged and finally merged into the single aggregate in the corridor between the lobes. Fluorescent images of the MTs (not shown) also confirmed the theoretical scenario.

We have developed a computational model of the MT/dynein system that demonstrates in silico aggregation of granules coated with dynein motors into a single focus and formation of a polar MT aster, thereby reproducing the phenomena experimentally observed in melanophore fragments. The model accounts for the self-organization phenomenon by assuming that (1) there is no physical MT transport, (2) dynein motors transport granules toward MT minus ends, (3) the MT nucleation rate is an increasing function of granule density, and (4) the rate of MT minus-end disassembly is a decreasing function of granule density. The model explains self-organization as a positive feedback loop based on the mutually enhancing processes of motor transport to MT minus ends leading to motor concentration and motor-mediated nucleation and stabilization of the MT minus ends.

Faithful quantitative reproduction of the observed aggregation and aster formation phenomena requires the existence of granule-associated factors able to both nucleate MTs, and slow down their minus-end disassembly. The model supports the qualitative assumption of Vorobjev et al. (Vorobjev et al., 2001) that granule-mediated nucleation and disassembly inhibition are equally important. Although there are a number of possible corresponding molecular mechanisms, we favor the hypothesis that dynein motors themselves are able to enhance the MT nucleation rate, perhaps by binding tubulin dimers and assembling a template for filament growth, and to transiently cap the MT minus end.

Some of the model parameters are available, either directly or indirectly, from published data. We estimated the rest of the parameters by fitting the theoretical results to the quantitative experimental data. In order to test the model, we made theoretical predictions for the dependence of the time of aggregation on fragment size and for pattern formation in the bi-lobed fragments and compared them with experimental observations. The semi-quantitative agreement between experiment and theory lend additional support to the model. Detailed quantitative measurements in fragments with altered MT dynamics and pigment density are necessary to test a few other model predictions: (1) decreased granule density leads to a more diffuse aggregate; (2) increased granule density leads to tighter aggregation into a single focus without transient aggregates; (3) specific constraints on the functional dependencies of the MT nucleation and minus-end disassembly rates on granule density.

The model has a number of limitations stemming from the constraints of mathematical solvability. These are: (1) inadequate treatment of the effect of steric repulsion of the pigment granules and of the observed process of granules being squeezed upward into the third, vertical dimension in tightly packed aggregates; (2) neglecting the possibility of granule transport toward the MT minus ends continuously, without numerous dissociations; (3) underestimating the probability of a granule associating with the same MT it just dissociated from. More importantly, the model does not capture the self-centering phenomenon that we describe briefly below.

In a different experimental preparation described elsewhere (Rodionov and Borisy, 1997b), in which granules are dispersed but MTs are pre-organized by the centrosome of the mother cell, the granules rapidly aggregate to the edge of the fragment where the minus ends are concentrated. The aggregate subsequently shifts to the center of the fragment (Rodionov and Borisy, 1997b). Movie 3 demonstrates the model-predicted pattern formation corresponding to this assay in which granules are initially distributed uniformly with all MT minus ends at the left and their plus ends at the right. The initial stage of aggregation is as observed in experiments but the aggregate stops far from the center, near the left edge of the fragment. Similarly, the 2D simulations (see Movie 4; Figs 8, 9) show that the aggregate is positioned close to but not exactly at the center of the fragment. Aggregation exactly at the center in some simulations (Fig. 7) is the consequence of the perfectly symmetric initial conditions and absence of stochastic effects. Therefore, there are molecular mechanisms responsible for the self-centering phenomenon that are not accounted for in our model. We hypothesize that force generation and length dependent differential buckling of MTs are essential for self-centering (Tran et al., 2001). We will address this problem in detail elsewhere. Despite these limitations, the model elucidates a pathway for aster formation in MT/motor systems.

Tissue cultures of black tetra melanophores were prepared as described previously (Rodionov et al., 1994). To prepare fragments, melanophore processes were dissected with microneedles with a 0.1 μm tip diameter. Aggregation of melanosomes in the fragments was triggered with 0.5 μM adrenalin. Phase contrast images of fragments were obtained using a Nikon TE300 microscope equipped with a Watek high-resolution video camera. Fragments with labeled MTs were obtained by injection of parental cells with Cy3-labeled tubulin (Rodionov et al., 2001) and subsequent dissection of fragments. Images of MTs were captured with a Photometrics CH350 back-illuminated cooled CCD camera.

In the absence of granules, the net minus-end depolymerization rate is equal to the plus-end growth rate, v_p. The function v(g) decreases monotonically with increasing granule density. A simple model explains the form given above. Suppose a MT minus end can be in one of two states, capped (C) or disassembling (D). The uncapping rate is the product of the uncapping rate constant and the density of capped MTs, k_uC. However, the reverse reaction depends on the local pigment concentration and the structural nature of capping. If s pigment granules are required to cap a minus end, then the transition probability is k_cg^s. At equilibrium, if D is the fraction of disassembling minus ends, then C=(k_c/k_u)g^sD is the fraction of capped minus ends. Assuming the disassembly rate is constant (v_p), the model predicts that the net disassembly rate is v(g)=v_pk_u/(k_u+k_cg^s). We have renamed and rescaled the parameters in the expression above for reasons of interpretation and nondimensionalization.

The parameter s is refered to as the Hill coefficient, which is a measure of the cooperativity required of pigment granules in stabilizing the minus ends. Large values of s correspond to the presence of a threshold density below which granules have little or no effect on minus-end depolymerization. The parameter g_v is the value of the threshold density at large s, or, for s∼1, the value of the half maximum density at which the net depolymerization rate is half what it is in the absence of granules. Fig. 4B illustrates two qualitatively different cases. The solid curve corresponds to s=1 (no threshold/cooperativity), while the dashed curve corresponds to s=3 (threshold/cooperativity behavior). Comparison of the theoretical and experimental data shows that the first case is the best fit: s=1, g_v=0.5. (We normalize g_v by the value of the granule density in the untreated fragment.) At values of s∼1, the model is insensitive to varying parameter g_v. However, increasing the Hill coefficient s introduces sensitivity to changes in g_v. This occurs because large values of s introduce threshold behavior to the stabilization phenomenon and, when the threshold is too high, the positive feedback mechanism is not triggered.

We choose half the size of the fragment, L, as the unit of length measurement, and the constant net polymerization rate v_p as the unit velocity so that the unit of time is L/v_p. We choose the (uniform) density of pigment in the untreated fragment, ḡ, as the scale for granule density. We define the characteristic nucleation rate in the fragment n̄=n(ḡ), where n(g) is given by Eqn 5. The characteristic scales of MT plus- and minus-end densities, p̄ and m̄, are determined by the product of the characteristic nucleation rate and the characteristic time for a MT end to travel across the fragment: p̄=μ=n̄L/v_p.

The model parameters are listed in Table 1. Some of them (v_p,v_g, L, k_off) are known from experiments. For others, we rely on numerical simulations and dimensional analysis combined with indirect experimental evidence. The behavior of the model depends crucially on three non-dimensional combinations of parameters ϵ=v_p/v_g, k₁=k_offL/v_g, and k₂=k_onL³n̄/(v_gv_p). The small parameter ϵ ∼0.067 is the ratio between the fast time-scale associated with granule transport to the slow time-scale associated with MT treadmilling. The fact that this parameter is so small means that granule density rapidly equilibrates to a slow changing MT distribution. This fact is conceptually an important feature of the model and simplifies both analytic and numerical calculations. The parameter k₁=5 is the ratio of the characteristic time of granule transport, L/v_g ∼5 seconds, to the average time of gliding before detachment, 1/k_off ∼1 second. Its value, relatively large compared with unity, corresponds to (and legitimizes) the claim that granules attach and detach frequently, rather than glide to the minus end of MTs directly upon attachment.

The value of parameter k₂ can be estimated from the following argument. Note that k₂ can be represented as k₂=k_on × [(n̄L/v_p) × L] × (L/v_g). Here n̄L/v_p is the scale of the MT plus- and minus-end concentrations and the expression in the square brackets is the characteristic MT polymer density. Therefore, k_on × [(n̄L/v_p) × L] is the characteristic frequency of granule attachment, which means k₂ is the ratio of the characteristic granule transport time, L/v_g, to the characteristic time in detached state, v_p/(k_onn̄L²). It was observed in the experiments with nascent fragments (Rodionov and Borisy, 1997b) that granules were rapidly transported across the fragment, almost at the same rate as the rate of gliding. This indicates that the characteristic time in the detached state is much less than the characteristic granule transport time and, therefore, parameter k₂ must be much greater than unity. Simulations show that, provided k₂ is large compared with unity, the behavior of the model is insensitive to variations in that parameter. The derivation of the 2D model equations, which do not depend on the value of k₂, relies on this fact. In simulations of the 1D model we used k₂=100.

The time-scale of nucleation, 1/(n̄L) (the average time required for the nucleation of a MT anywhere in the fragment), can be estimated by noting that the MT array in a fragment is capable of turning over completely in a period of about 10 minutes. Assuming the MT array consists of 100-300 MTs, we infer a time-scale for the nucleation of a single MT to be few seconds. Using these indirect estimates for the non-dimensional combinations of parameters, we calculated the dimensional values of the model parameters and listed them in Table 1. In particular, we predict that n̄∼0.001 μm² second) and that k_on ∼1/second.

Finally, let us note that although we use the value k₁=5 in 1D, this value turns out to generate very diffuse aggregates in 2D. Good results are achieved at a value of k₁=500, two orders of magnitude greater than our original prediction. This modification is discussed and justified in Appendix 6, below.

Asymptotic analysis of the 1D model provides insight into the behavior of the system, allows us to find analytical solutions to the model equations and, more importantly, suggests a generalization of the 1D model to the more realistic 2D case.

The effective diffusion is due to the random motion of granules on the MT array. As for the advection term, as granules come on and off quickly, they essentially sample the local MT population. Thus, their velocity is given by the difference between the probability of attaching to a left mover, N_l/(N_r+N_l), and that of attaching to a right mover, N_r/(N_r+N_l).

To implement the 2D model, at each step we (1) nucleate new MTs and eliminate any MTs whose minus ends reach the boundary, as described above; (2) update the positions of the plus and minus ends of the MTs using equations (Eqn 14); (3) find the velocity field using Eqn 15; and finally, (4) solve the diffusion-advection equation (Eqn 13) using the updated velocity field and no-flux boundary conditions. The latter equation is solved on a rectangular grid using standard finite difference numerical methods (Garcia, 1994).

In the bi-lobed geometry, we must account for the fact that essentially no MTs enter the neck until late in the process. The first obvious correction to try is to reduce the diffusion coefficient (1/k₁) to a fraction of its original value in the neck (setting it to zero causes severe computational difficulties). Unfortunately, this approach fails because of the pseudo steady state assumption – in the model, pigment instantly equilibrates to the stable steady state associated with the current MT array. Through the neck, this instantaneous equilibration requires a unrealistically high rate of flux. Therefore, we solve the pigment equation independently on each half of the domain, applying a no flux boundary condition along the midline of the neck. This guarantees that there is no transport of pigment across the neck so that each half has a conserved quantity of pigment. This simplification can be justified by the fact that in the narrow neck the steric repulsion of the granules prevents large fluxes. We allow MTs to nucleate and grow everywhere in the fragment.

To obtain quantitative estimates of the nucleation and minus-end stabilization behavior, we assume that the packed pigment aggregate is fixed at the center of the fragment. The aggregate is circular with radius r=0.15 and has a constant granule concentration. Both the area and initial constant granule density are normalized by unity, so g_a is simply the total number of granules, 1, divided by the area of the aggregate: g_a=1/(πr²). Under these conditions, only three values of the granule density are present: 1, before aggregation; g_a, in the aggregate; and 0, outside the aggregate. Therefore, we have to define three values of the functions n(g) and v(g): n(0),n(1),n(g_a),v(0),v(1),v(g_a) that provide the best fit for the available data.

The following argument provides a rough estimate for these values. The minus end of a MT nucleated in the aggregate takes around r/v(g_a) time units to escape the aggregate. Meanwhile, the plus end takes about 1 time unit to reach the fragment edge. Assuming 1<r/v(g_a), the typical MT spends 1 time unit growing with its minus end in the aggregate followed by about [r/v(g_a)–1] time units with its plus ends at the fragment edge and its minus end in the aggregate. Once free of the aggregate, the minus end reaches the fragment edge about 1 time unit after escaping the aggregate. The proportion of time growing is given by 1/[r/v(g_a) +1]. Both ends are stable for a proportion of time given by [r/v(g_a) –1] / [r/v(g_a) +1]. Finally, shortening occurs for a proportion of time given by 1/[r/v(g_a) +1]. The experimentally measured proportions are approximately 0.1 (growing), 0.8 (both ends stable) and 0.1 (shortening). This means that r∼9v(g_a). Estimates of the aggregate size from the experimental data give the value of r∼0.15, and therefore v(g_a)∼0.0167. Under the assumption that MTs treadmill in the absence of granules, we set v(0)=1. The choice of a value for v(1) is determined qualitatively from the behavior immediately following adrenalin treatment. Before treatment, a sufficient MT array must develop in order to get a loose aggregate but not so much so that stochasticity is lost. A value of v(1)∼0.3 fits these criteria.

The value of n(0) represents the level of spontaneous nucleation in the absence of granules and must be small compared with 1. Otherwise, the number of treadmilling MTs is too high (recall that treadmillers represent ∼1% of MTs). The value of n(1)=1 is determined by the scaling. The value of n(g_a) has to be found numerically by matching the observed increase in total length of MT polymer after aggregation.

We thank G. Borisy for fruitful discussions and help at early stages of this work, and J. P. Keener for useful suggestions on modeling issues. A.M. and E.C. are supported by the NSF and HIH grants DMS-1097746 and NIGMS GM068952-01, respectively, and by the Burroughs Welcome Fund through the Program in Mathematics and Molecular Biology. V.R. is supported by the NSF and NIH grants NCRR RR13186 and NIGMS GM-62290, respectively.