Stochastic reaction-diffusion-dynamics model of mitotic spindle assembly

Reaction-diffusion master equation formalism.

We employed the RDME approach to modeling biochemical reactions and molecular transport [30,31,33]. The volume (~850 μm³) of the ellipsoidal shaped cell is divided into a large number of subvolume elements (subcells) of dimension l_SV = 250 nm (Table A in S1 Text), and the biomolecules in the cell are distributed among the subcells (Fig 1). Biochemical reactions are allowed only between molecules within a subcell. The biomolecules can diffuse randomly between the nearest-neighbor subcells. The state of the cell X is specified by the number of biomolecules and biomolecular complexes (e.g. KT-MT attachments) x_j,v of each type j = 1, 2, …, J in each subcell ν = 1, 2, …, V. The time evolution of the probability for the cell to be in state X is given by the sum of contributions from biochemical reactions, described by the reaction operator R, and from molecular diffusion events, described by the diffusion operator D, i.e.

(1)

In the second line in Eq 1, which describes the rate of change of P(X,t) due to biochemical reactions, x_v is the column vector containing the number of molecules in the ν-th subcell, α_μ(x_v) is the reaction propensity for the μ-th reaction (μ = 1, 2,…, M) to occur in the ν-th subcell, and S_μ is the μ-th column of the J×M stoichiometry matrix S, which describes the changes in the number of molecules when the μ-th reaction occurs. In the absence of molecular transport, Eq 1 reduces to the spatially homogeneous chemical master equation (CME): , which is widely used to stochastically model chemical reactions in a well-stirred volume. Therefore, the RDME extends the CME to account for spatial degrees of freedom. In the third line in Eq 1, which describes the rate of change of P(X,t) due to molecular diffusion events, d_j is the diffusion propensity for a biomolecule of type j to move from subcell v to the next-neighbor subcell v+ξ, where ξ is a next-neighbor subcell in the ±x, ±y, and ±z direction (total of 6 next-neighbor subcells) denoted by the unit vectors i,j,k (see Fig A in S1 Text), x_j,ν is the number of biomolecules of type j in subcell v, and 1_j,ν represents a single molecule of type j in subcell v. In the absence of chemical kinetics, Eq 1 reduces to the diffusion equation, .

The RDME (Eq 1) was sampled numerically using the Gillespie approach, which is based on the propensities, rather than probabilities of chemical reactions used in more traditional Monte Carlo approaches. In the Gillespie approach, the probability that the μ-th reaction will occur within the next time interval between t+τ and t+τ+dt is given by P₀(t+τ)c_μh_μdt, where P₀(t+τ) is the probability that at time t+τ no reaction has occurred in the previous time interval (t, t+τ). The reaction propensity for the μ-th reaction is given by α_μ = c_μh_μ and total propensity for all M reactions is . For a unimolecular reaction μ = 1, 2,…, M (e.g. MT-Ndc80 complex dissociation, MT growth, MT shortening, MT catastrophe and MT rescue; see Tables 3 and 4) to occur in the ν-th subcell with the rate constant k = k_μ, the rate equation for a chemical species of type A is . We define c_μdt = cdt to be the probability that a particular combination of reactants will interact through the same reaction μ in the time interval dt. If h_uni is the total number of distinct molecular reactant combinations at time t, then for a single molecule of type A, c = k and h_uni = x_A, and the reaction propensity is α_uni = ch_uni = kx_A. For a bimolecular reaction (e.g., Ndc80 phosphorylation, dephosphorylation and MT-Ndc80 complex formation; see Table 4), the rate equation for species A and B is . If h_bi is the number of distinct molecular reactant combinations for a bimolecular reaction at time t, then for a single combination of A and (l_SV is the size of subcell) and h_bi = x_Ax_B, and the reaction propensity is . The diffusion propensity for a molecule of type j (j = 1, 2, …, J) to move from subcell v to the next-neighbor subcell v+ξ is given by , where D_j is the diffusion constant (see Table 3 and Table A in S1 Text).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Biochemical reactions at kinetochore-microtubule interface.

Enzymatic reactions (e.g., phosphorylation and dephosphorylation) and association-dissociation reactions, which involve the MT associated protein Ndc80 linking MTs with KTs, and the reaction rate constants and characteristic timescales. The subscript p = 0, 1,…,6 denotes the number of phosphate groups attached to Ndc80 and changes in the rate constants and reaction propensities for MT-Ndc80 complex dissociation.

https://doi.org/10.1371/journal.pcbi.1010165.t004

The RDME formalism asymptotically approximates the Smoluchowski diffusion-limited reaction method when the average time for the next reaction to occur in a subcell τ_R is much longer than the average time until the next diffusion event τ_D, i.e. τ_R≫τ_D [30–32,51]. In the SRDDM, the slowest diffusing particles are Aurora A and B, for which the diffusion timescale is τ_D = 4.3×10⁻⁴ s (Table 3), and the most rapid reaction is MT-Ndc80 complex formation, for which the characteristic time is τ_R = 2.6×10⁻³ s (Table 4). Therefore, the diffusion timescale is ~10-fold shorter than the characteristic reaction time. This large separation of timescales for kinetics and diffusion (τ_D≪τ_R) justifies our using the RDME formalism. We adapted the multi-particle diffusion (MPD) method [52] and the next-subvolume method (NSM) [51,53] for exactly sampling the RDME (Eq 1).

Force field for mechanically active components

Langevin dynamics formalism describing mechanically active components

We used the Langevin Dynamics (LD) approach to model the dynamic mechanical processes. The cell configuration r is specified by the positions of all mechanically active components r_i, i = 1, 2,…, N, where N is the total number of components. The cell dynamic evolution was followed by integrating the Langevin equations in the overdamped limit for each position r_i of each mechanically active component, (10)

In Eq 10, U(r) is the total potential (mechanical) energy function (see Eq 2), γ is the friction coefficient, and g_i(t) is the Gaussian zero-average random force with the variance σ². The Langevin equations were propagated forward in time with the timestep δt = 50 ps at room temperature (T = 300 K) using water viscosity (Table A in S1 Text). For the KT beads (R_CH = 362.5 nm, Table A in S1 Text), this corresponds to the friction coefficient γ = 6πηR_CH = 6.8×10⁶ pN ps/nm. For the MT beads (R_MT = 12 nm, Table A in S1 Text), this corresponds to the friction coefficient γ = 4.5×10⁵ pN ps/nm. The variance σ² is related to the diffusion constant (k_B is the Boltzmann’s constant), and so σ² = 2k_BTγ. In Eq 10, is the deterministic force and the second term contains the random force.

Dynamics of microtubules: growth, shortening, catastrophe, and rescue

The MT growth is described by the equation for the MT length l, where k_gr = 5.0 s^-1 is the rate of growth and Δl = 24 nm is the discrete increment of length. Similarly, the MT shortening is described using the equation , where k_sh = 18.6 s^-1 is the rate of shortening (Table 3 and Fig D in S1 Text). The reaction rate constants k_gr and k_sh are chosen to recover the experimental rates of MT growth and shortening, v_gr = 7.5 μm/min and v_sh = 27 μm/min (Table 3 and Table A in S1 Text). During the process of dynamic instability, MTs switch between phases of growth and shortening. The frequency of catastrophe ω_cat is set to the experimentally observed value [57], ω_cat = 2.5×10⁻³ s^-1; the frequency of rescue ω_res is set to the experimentally observed value [57] ω_res = 3.0×10⁻² s^-1 (Table 3 and Table A in S1 Text).

Dynamics of microtubules: interactions with chromosomes

In the SRDDM, MTs physically interact with CHs. These interactions are the following: 1) When an MT overlaps with a CH arm (or a centromere), the excluded volume interaction adds non-zero energy to the system (Eq 7) which generates a pushing force (the first term in the right-hand-side of Eq 10). The direction of this pushing force is along the line connecting the centers of two interacting spherical beads representing the MT and CH arm. This force is pushing the CH away from the MT, while also bending the MT away from the CH. As all MTs are either growing or shortening all the time, this excluded volume interaction lasts for a finite time interval. With the sizes of the beads and dynamic instability rates we use, the pushing process lasts for a short time rarely exceeding 2 s (see Fig E in S1 Text). After that, the MT effectively slides off the CH. The magnitude of the pushing force is determined by the parameter ε in the Lennard-Jones potential (Eq 7) chosen so that the average pushing force is in the 10-pN range (see Fig E in S1 Text) to conform with the experimentally established forces [46,58]. 2) When the plus-end bead of a growing MT is in the vicinity of a KT, which contains Ndc80 linker ends, the MT-Ndc80 bond formation takes place with the rate constant k = 3.8×10⁹ s^-1M^-1 (see Table 4). At the Ndc80 density used (~ 750 Ndc80 linkers per KT; Table 2 and Table A in S1 Text), this results in the formation of the Ndc80-MT attachment within, on average, 2.6×10⁻³ s (Table 4) as soon as the MT plus-ends bumps into the KT. 3) When an MT-Ndc80 attachment occurs, an immediate catastrophe takes place, and the MT starts to shorten. 4) During the KT-MT attachment time interval, the shortening MT stretches the Ndc80 spring thereby exerting the 10-pN average pulling force on the KT (see Fig E in S1 Text). The direction of the pulling force is along the line connecting the MT plus-end and the Ndc80 end bound to this plus-end.

Numerical implementation

The SRDDM was mapped into CellDynaMo package (CUDA language), fully implemented on a GPU. In LD, the particle-particle interactions (e.g. excluded volume interactions, stretching and bending within the filament structures, such as MTs and CHs, and KT-MT attachment) are the computational bottleneck. However, these interactions are described by the same empirical potential energy function (Eqs 2–9). Therefore, when running Langevin Dynamics on a GPU, it is then possible to execute the same operation, e.g. generation of random forces, calculation of the potential energy, evaluation of forces, integration of Langevin equations of motion, for many particles at the same time. When mapping the RDME, we implemented the next-subvolume method (NSM) extension [51] of the original Gillespie algorithm [59,60]. In CellDynaMo, we implemented the multi-particle diffusion (MPD) approach to the reaction-diffusion master equation [52]. Numerical routines for the generation of (pseudo)random numbers (Hybrid Taus) for RDME and LD are described in our previous publications [61,62]. To achieve top performance on a GPU, all the numerical algorithms implemented in RDME and LD have been recast into a data-parallel form so that the computational threads run the same instruction stream, but on different data sets (i.e., subcells and particles). We made the tasks compute-intensive so that, most of the time, the GPU performs computations rather than reading and writing data. These efforts enabled us to reach the biologically important timescales. For example, it takes ~72 hours of wall-clock time to generate a few ~30 min trajectories of cell dynamics on a contemporary graphics card GeForce GTX 1080.

Benchmark test simulations

To assess possible errors in the numerical implementation of the LD module for the mechanical components, and in the RDME description of chemical kinetics and molecular transport, we carried out benchmark test simulations. Here, we summarize the main results; technical details are given in S1 Text.

Langevin Dynamics: In CellDynaMo, the Langevin equations of motion in the overdamped limit for each position r_i of each mechanically active component are integrated numerically using the first-order integration scheme (Ermack-MacCammon algorithm) [63]: (11)

First, to access the accuracy of the numerical integration, we carried out short 1-μs simulation runs (20,000 integration steps) for a small system of three beads interconnected by the harmonic springs (Fig F-A in S1 Text). Benchmark simulations of the cell dynamics have been carried out at zero temperature with the time step δt = 50 ps (see Table A in S1 Text). For this three-body system, the equations for the forces and displacements for all beads are known exactly (Eqs 1–3 in S1 Text), and so the results obtained both numerically and analytically can be directly compared. The displacements Δx₁₃ and Δx₂₃ displayed in Figs F-A, F-B in S1 Text shows excellent agreement between the exact and numerical results. The numerical error does not exceed 1.5% for t>1 μs. For a 5-pN force on bead 3, we obtained Δx₁₃ = 6.62×10⁻¹ nm (simulations) vs. 6.66×10⁻¹ nm (exact), and Δx₂₃ = 6.21×10⁻¹ nm (simulations) vs. 6.27×10⁻¹ nm (exact). For a 50-pN force, we obtained Δx₁₃ = 6.32 nm (simulations) vs. 6.33 nm (exact), and Δx₂₃ = 5.91 nm (simulations) vs. 6.33 nm (exact). We also investigated the dependence of the relative error between the asymptotic values of the particle displacements obtained numerically (Δx_sim) and analytically (Δx_exact) for the 50-pN pulling force on the integration timestep δt = 5×10⁻²–5×10² ps (see S1 Text for more detail). Err(Δx) was found to be very low (<1.8%; see the inset to Fig F-B in S1 Text).

Next, we investigated the dependence of the relative error between the average displacement for N = 100 Brownian oscillators obtained numerically (〈Δx(t)〉_sim) and exactly (〈Δx(t)〉_exact) on the cytoplasmic viscosity η = 1, 5, and 10 cPs and on the integration timestep δt = 5×10⁻²–5×10² ps (see S1 Text for more detail). The numerical and analytical results practically collapse on the same curve (Fig F-C in S1 Text), and Err(〈Δx(t)〉) was found to be very low, i.e. <0.4% for variable solution viscosity (the inset to Fig F-C in S1 Text) and <0.3% for variable timestep (Fig F-D in S1 Text; see the inset to Fig F-D in S1 Text for average relative error). Hence, our choice of δt = 50 ps as the timestep for Langevin Dynamics is reasonable.

We also carried out benchmark test simulations for a single KT pair (centromere; CH) to estimate the coefficients of one-dimensional translational diffusion and one-dimensional rotational diffusion based on the numerical output from the Langevin Dynamics simulations. These were then be compared with the exact analytical results for the translational diffusion coefficients, , and rotational diffusion coefficient, , of a cylinder (see Fig G in S1 Text). In these simulations, we set the CH height to L_CH = 4R_CH = 1.450 μm, and the CH diameter to d_CH = 2R_CH = 0.725 μm. We performed these test runs with a short timestep δt = 5 ps to obtain more detailed trajectories with a large number of KT positions and KT orientations (see Fig G in S1 Text). By averaging over 3 runs (each of 2.5-min duration), we obtained D_x,sim = 0.154 μm²/s (from simulations) vs. D_x = 0.152 μm²/s (exact result), and D_θ,sim = 0.866 rad²/s (from simulations) vs. D_θ = 0.869 rad²/s (exact result).

Reaction-Diffusion Master Equation–Brownian diffusion: We assessed the accuracy of the numerical implementation of the diffusion part of RDME by performing simulations of molecular transport, for which the exact distribution of the particles’ displacements is described by Gaussian statistics (Eq 4 in S1 Text). We carried out several short (100 s) benchmark simulations, in which we placed 10⁴ molecules of Aurora B in the central subcell x₀ = 0 at time t = 0 (see Fig H in S1 Text). For Aurora B, we set the diffusion constant to D = 7.3×10⁷ nm²/s based on the Einstein-Stokes relationship for spherical particles (S1 Text) with the radius of R_A = 2.9 nm [64] (Table A in S1 Text). We ran the multi-particle diffusion algorithm (MPD; see Fig H in S1 Text) and observed spreading of molecules at later time points t = 1, 2, 5, 10, and 20 s, which correspond to 2.50×10⁴, 5.00×10⁴, 1.25×10⁵, 2.50×10⁵, and 5.00×10⁵ steps of iteration. The non-parametric density estimates of the distributions of particles’ displacements constructed based on the test runs are compared with the exact probability distribution curves in one dimension in Fig H in S1 Text, which shows excellent agreement between the numerical and exact results. We also compared the variance from the simulations with the exact values , for t = 1, 2, 5, 10, and 20 s. We found 1.4×10⁸ nm² (simulations) vs. 1.5×10⁸ nm² (exact) for the t = 1 s time point, 2.9×10⁸ nm² (simulations) vs. 2.9×10⁸ nm² (exact) for t = 2 s, 7.1×10⁸ nm² (simulations) vs. 7.3×10⁸ nm² (exact) for t = 5 s, 1.4×10⁹ nm² (simulations) vs. 1.5×10⁹ nm² (exact) for t = 10 s, and 2.8×10⁹ nm² (simulations) vs. 2.9×10⁹ nm² (exact) for t = 20 s. The numerical errors are all below 3%.

Reaction-Diffusion Master Equation–Biochemical kinetics: Next, we assessed the accuracy of the numerical implementation of the reaction part of the RDME. We carried out benchmark simulations for the consecutive two-step irreversible kinetics, , with species A, B and C and reaction rate constants for the first and second steps, k₁ and k₂, and for the single-step reversible kinetics, A⥨B, with the reaction rate constants for the forward (A→B, k₁) and backward (B→A, k₋₁) steps. For the two-step irreversible kinetics, the time-dependent populations p_A(t), p_B(t), and p_C(t) are given by Eqs 8, 9, and 10 in S1 Text. For the single-step reversible kinetics, the populations p_A(t) and p_B(t) are given by Eqs 13 and 14 in S1 Text. We compared the results of numerical calculations of p_A, p_B, and p_C with the initial conditions p_A(0) = 1, p_B(0) = p_C(0) = 0 and reaction rate constants k₁ = 1 s^-1 and k₂ = 2 s^-1 with the exact solutions given by Eqs 8, 9 and 10 in S1 Text. We compared the results of calculations of p_A and p_B with the initial conditions p_A(0) = 1 and p_B(0) = 0 and reaction rate constants k₁ = 1 s^-1 and k₋₁ = 3 s^-1 with the exact solutions given by Eqs 13, 14 in S1 Text. The results are presented in Fig I in S1 Text, which shows that the time-dependent populations obtained analytically and numerically practically collapse onto the same curves.

Optimal MT detachment rate maximizes correctness of KT-MT attachments

In the remainder of the paper, we apply our modeling toolkit to simulate the Search-and-Capture process for a single CH. The process of a single CH incorporation into the spindle is already so complex that in this paper we limit ourselves to this situation and leave the case of multiple CHs for the next study. Initially, there are two immobile CSs with dynamic MT asters and completely unconnected CH. We allow 30 min of biological time for MT and CH movements and attachment kinetics to play out, record the spindle angle (angle between the vector connecting sister KTs and the vector connecting the centrosomes/spindle poles), KT-KT distance, distance from the CH to the spindle equator, and number of MTs attached to KTs as functions of time. The results reveal complex and nuanced dynamics that results either in correct, amphitelic attachments, or, in a significant fraction of cases, in erroneous merotelic attachments or incomplete monotelic attachments. In a very few cases, the result is erroneous syntelic attachment. Analysis of the respective statistics gives unexpected insights about key feedbacks responsible for accuracy of the spindle assembly and pitfalls that the nascent spindle must circumvent.

We start with the ‘naked’ CH without chromosome arms. In this case, the CH is just a roughly spherical centromere with two sister KTs on the opposite sides of the centromere. We also start with the centromere spring to be very stiff, so that the KT-KT distance changes very little. This simplest case sheds light on one of the most basic questions–how the attachment accuracy depends on the rate of MT detachment from the KT. Previous modeling studies suggested that faster KT-MT detachment improves the accuracy, allowing destabilization of erroneous attachments. The main regulators of the KT-MT binding strength are members of the family of Aurora protein kinases and respective phosphatase enzymes. We model explicitly Aurora A (on the spindle poles) and Aurora B (on the centromeres), which reduce the KT-MT attachment strength by phosphorylation of the Ndc80 tails linking the MTs and KTs, and also the phosphatase enzymes, which make the KT-MT bond stronger via dephosphorylation of the Ndc80 tails. Therefore, it is only the number of Phosphatase and Aurora B (AB) molecules around the KT-MT interface that affect the turnover of individual KT-MT attachments. In the simulations, phosphatase, for simplicity, is uniformly distributed in the interior of the cell, while AB undergoes a confined diffusion and so is distributed with a spherically symmetric Gaussian distribution centered at the center of the centromere. The distribution width is such that AB spreads up to approximately 250 nm away from the centromere surface (see Fig 4A), where the 65-nm long Ndc80 linkers are located. The kinetic, force field and cell morphology parameters, including the Phosphatase-to-AB (P:AB) ratios, are given in Tables 1, 2, 3 and 4 and Table A in S1 Text.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Exploring the role of Phosphatase and Aurora B kinase on dynamics of kinetochore-microtubule attachments.

A) Initial position and orientation of the KT pair inside the cell. The KT pair is placed close to the equatorial plate, shifted by 1.5 μm below and along the axis perpendicular to the axis of the spindle. Orientation of the KT pair axis is along the axis perpendicular to the spindle axis. The AB particles are spatially distributed around the center between the KT pair as shown in the blowout. B) Time dependent evolution of KT-KT orientation angle (solid lines; left y-axis) and distance between KT-pair and equatorial plate (dash-dotted lines, right y-axis) for the Phosphatase to Aurora B (P:AB) ratio = 1:100 (blue), 1:10 (red), and 1:1 (black). Each curve shows results from a single simulation run. C) Probability to find each type of attachment for different P:AB ratio = 1:100 (blue bars), 1:10 (red bars), and 1:1 (black bars). Statistics was collected from n = 8 independent runs for each of these three cases. D) The P:AB ratio influences the average bond lifetime and the frequency of attachment-detachment switches. Shown is evolution of the maximum number of microtubules (MTs) attached to a single KT from a single centrosome for P:AB ratio = 1:100 (blue lines), 1:10 (red) and 1:1 (black) during 30 min of simulation. Changes in the KT-MT attachment status are shown in the inset. E) Snapshot of the final state of the system (red and cyan spheres are particles of Aurora B kinase and Phosphatase).

https://doi.org/10.1371/journal.pcbi.1010165.g004

To create an initial configuration, we placed the naked CH close to the equatorial plate, shifted it by 1.5 μm ‘downward’ along the axis perpendicular to the axis of the spindle, and oriented the KT-KT axis perpendicularly to the spindle axis (Fig 4A). This initial arrangement eliminates any positional and orientational bias, providing the same probability for MTs from both CSs to reach the same KT, which in turn should lead to a higher probability of forming merotelic attachments. In other words, these initial conditions, which we used in all simulations, are ‘making it hard’ for the CH to achieve the correct amphitelic state (see also S1 Movie).

Before we investigate how the Ndc80 phosphorylation state affects the KT-MT attachments, we first examine the effect of different P:AB ratios on the degree of Ndc80 phosphorylation N_phos. We ran 3 sets of simulations for a single naked CH without Langevin Dynamics (LD is switched off), only with kinetics and diffusion (RDME is switched on), for 3 different P:AB ratios for 15 min of biological time. The number of AB enzymes was fixed at about 100, while the number of phosphatase molecules was varied from 10 to 100 and to 1000. Thus, the smallest P:AB ratio (P:AB = 1:10) results in the largest degree of phosphorylation, N_phos =6.3±0.7, which is close to the highest possible value (7 phosphorylated sites) that would result in the least stable KT-MT attachments. For larger P:AB ratios, N_phos = 3.3±0.4 (P:AB = 1:1) and N_phos = 0.1±0.0 (P:AB = 10:1), which should result in more stable KT-MT attachments.

Interestingly, the characteristic number of KT-MT attachments to a CH that remains close to the spindle equator (two of the three cases analyzed in Fig 4B, 4C and 4D) grows roughly linearly with similar rates despite different KT-MT detachment rates (Fig 4D). The reason is that roughly the same set of MTs growing into a sector that allows MTs to reach the CH, interact with the KTs. Thus, a detached MT often shortens, gets rescued, grows in roughly the same direction, and gets attached again. Note that after 30 min of biological time, the characteristic number of KT-MT attachments approaches 15 (Fig 4D and 4E), in agreement with the experimental observations in eukaryotes [65,66].

Eight independent simulation runs carried out for each P:AB ratio, provided the statistics shown in Fig 4C, which illustrates that slow KT-MT detachment leads to very inaccurate spindles–only 25% of amphitelic attachments evolve, while the majority, 50%, of the attachments are merotelic. The reason is that when a KT is exposed to both poles, there is a large probability for MTs from both poles to attach to the same KT, thus, creating the merotelic attachment. If the detachment rate is too slow, there is very little chance to break the erroneous attachments. Note that a surprisingly noticeable fraction of attachments, 25%, are monotelic. A close look at the CH trajectory, orientation and KT-MT attachment numbers (black curves in Fig 4B and 4D) hints that the origin of such attachments is the ‘geometric vicious cycle’; that is, if one of the KTs ‘acquires’ significantly more attachments to one pole than its sister KT to another pole, then the mechanical pulling starts to bring the CH closer to the first pole. The closer the CH is to the first pole, the more attachments the proximal KT accumulates, the stronger the pull to the first pole. The distal pole is becoming farther away from the CH, and the sister KT does not make enough attachments to the distal pole to stop this process. Sometimes, we observe that the sister KT does not have time to make even a single attachment from the distal pole. More often though, initially monotelic attachments turn into either amphitelic ones, when a MT from the distal pole reaches the unconnected KT (as happens around 15 min in the case shown in Fig 4D, black curves) or merotelic ones when an MT from the proximal pole reaches the KT previously unconnected to this pole (as happens around 20 min in the case shown in Fig 4D, black curves). Syntelic attachments were observed very rarely because rare transient syntelic connections rapidly turn into merotelic ones: in the syntelic state, when both sister KTs are attached to one pole, the CH is oriented so that both KTs are also exposed to the opposite pole causing rapid merotelic connection.

Simulation with an intermediate detachment rate revealed a dramatic improvement in attachment accuracy: while the fraction of the monotelic connections, 25%, is unchanged, the percentage of the amphitelic connections, 62%, increased more than two-fold, while the percentage of the merotelic connections, 12%, decreased four-fold. The reason is that many merotelic connections transiently lose a minor number of ‘incorrect’ MTs connecting one of the KTs to the ‘wrong’ pole. When this happens, the ‘beneficial’ geometric feedback kicks in: transient amphitelic connections, by pulling sister KTs to the opposite poles, orient the CH so that the spindle angle approaches zero, which enables the centromere body to protect each of the KTs from the ‘wrong’ pole and accelerates formation of attachments from the ‘right’ pole. Thus, the correct attachments improve the spindle angle, while better angles promote accurate attachments. The geometric vicious cycle does not occur because amphitelic attachments do not allow the CH to deviate from the spindle equator too much (see also S1 Movie).

Interestingly, when we explored the fastest KT-MT detachment rate, the accuracy worsened: though the number of amphitelic attachments decreased just a little, from 62% to 56%, the number of merotelic attachments increased significantly, from 12% to 38%. The reason is that too rapid MT attachment dynamics does not stabilize the amphitelic state, allowing frequently detaching CHs to rotate slightly into unfavorable positions exposing the KTs to wrong poles and causing merotelic attachments. Note also the drastic decrease in the number of monotelic attachments in this case. The reason is that the rapid dynamics, while ruining the beneficial geometric feedback, also disables the geometric vicious cycle: when one of the KTs starts moving toward one of the poles, the number of MTs from that pole pulling that KT does not grow fast enough due to the rapid detachment, and the movement slows down giving time for MTs from the opposite pole to attach to the sister KT and to stop the divergence from the equator. In summary, we conclude from this case study that there is an optimal intermediate KT-MT attachment rate maximizing the spindle accuracy.

Aurora A kinase has a small positive effect on the attachment accuracy

Next, we investigated the influence of Aurora A kinase (AA) in the vicinity of the spindle poles (Fig 5A), which was switched off in the simulations for Fig 4 focused on the accuracy of KT-MT attachments. This case study was motivated by the possibility that AA could help in breaking the KT-MT attachments responsible for bringing the CH too close to one of the spindle poles. The idea is that near one of the poles, AA phosphorylating activity would increase the KT-MT detachment rate, which in turn would allow the CH to move further away from the pole and back to the center of the spindle, where the sister KT could then acquire proper KT-MT connections.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Exploring the effect of the Aurora A presence.

A) Snapshot of the final KT pair position and orientation (amphitelic attachment) for a simulation with AA. B) Probability to find each type of attachment for the system without AA (blue bars) and for the system with AA (red bars). Statistics were collected from n = 20 simulation runs for both cases. C) Addition of AA kinases to the system changes the Phosphatase to Aurora B ratio close to centrosomes (CSs). This influences the average bond lifetime and the frequencies of attachment/detachment switches. Blue line shows the evolution of total number of microtubules (MTs) attached to a single KT pair from both CSs for the case study without AA, red line shows the number of MTs for the case study with AA. D) An example of how adding AA kinases to the system corrects the trajectory of the KT pair. Lines show the 2D projection (xz-plane) of the KT pair trajectory during 30 min simulations. For the case study without AA (blue line), KT pair reaches one of the CSs and for the case study with AA (red line), KT pair stops right before the cloud of AA kinases.

https://doi.org/10.1371/journal.pcbi.1010165.g005

We carried out 20 independent simulations for the system without AA and then ran 20 simulations with AA present, whose Gaussian distributed concentration profile is centered on both poles (see Fig 5A and 5D), for 30 min of biological time. In these simulations, we used the highest P:AB ratio for AB, corresponding to the fastest KT-MT detachment kinetics; all the kinetics, mechanical and geometric parameters are given in Tables 1, 2, 3 and 4 and Table A in S1 Text. The addition of AA improved the accuracy (Fig 5B) but only slightly: the probability of formation of amphitelic attachments slightly increases, while the probability of formation of merotelic attachments decreases (Fig 5B). The AA effect on the total number of MTs attached to the naked CH can be seen in Fig 5C. In the absence of AA, the number of attached MTs starts to grow abruptly (a signature of the geometric vicious cycle) and changes little after the CH approaches one of the poles. In the presence of AA, the number of attached MTs grows more evenly, and, importantly, fluctuates significantly near the maximum, because AA frequently breaks the attachments to the proximal pole. As a result, occasionally, the effect of AA stopped the movement of the naked CH to one of the poles (Fig 5D). Equivalent behavior in the simulations without AA was not observed. In summary, AA does have a small but overall positive effect on correcting the KT-MT attachments and reversing incorrect motion of KT pairs toward the CS poles.

Chromosome arms and centromere softening improve the accuracy of the spindle assembly

For the purpose of simplicity, in our previous simulations (Figs 4 and 5), we used a stiff centromere without the CH arms. Here, we can use those simulations as benchmarks to examine how the addition of the CH arms and softening of the centromere affects the number and distribution of KT-MT attachment types and the final orientation of the CH. The mechanical effect of the soft centromere is hypothesized to be important for stabilizing the amphitelic attachments: when sister KTs are pulled in the opposite directions by amphitelic attachments, the sister KTs deform the centromere and shift away from the ‘Aurora B cloud’. This decreases the degree of phosphorylation of the Ndc80 linkers, decreases the KT-MT detachment rates and stabilizes the connections. The erroneous connections either do not generate the stretch (syntelic) or generate a smaller stretch (merotelic) on CH (Fig J-E in S1 Text), and the detachments are faster, which increases the probability of a correction. It is impossible to simply intuit what the effects of the CH arms on the accuracy are, and this is where the detailed 3D simulations become crucial.

We carried out 2 sets of simulations for 30 min of biological time, In the first set, we ran 20 independent simulation runs for the naked CH with stiff centromere, and in the second set we then repeated 20 simulation runs with an entire CH (centromere plus CH arms) and stiff centromere. All parameters are given in Tables 1, 2, 3 and 4 and Table A in S1 Text. Initially, we positioned the centromere as described above (see Fig 4A; see also S2 Movie) with straight CH arms perpendicular to the spindle axis and the path from the poles to the KTs. Indeed, the accuracy improved with the addition of CH arms and with softening the centromere (Fig 6A). The characteristic dynamics of one example of an optimal path to amphitelic attachment is shown in Fig 6C and 6D. There was never an erroneous attachment in this case study. The spindle angle rapidly turned to zero, the number of MTs attaching to the KTs increased steadily and equally for both KTs, thus increasing the stretching force and pulling the KTs apart (see Fig J-E in S1 Text). One of the reasons for this improved accuracy is the soft centromere spring, but we also expect that since CH arms are large structures, they might slow the CH movement. This would keep the CH near the spindle center (see Fig 6B), where there is a higher chance of formation of correct attachments from both poles. We will confirm this intuition in the case studies described below.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Describing chromosome arms and flexible corona surface in Stochastic Reaction-Diffusion-Dynamics Model.

A) Probability to find each type of attachment for a single chromosome (CH) with CH arms. Blue and red bars represent statistics for kinetochore pairs (KTs) without and with CH arms, respectively. The statistics are based on n = 20 simulation runs for each case. B) Snapshot after 30 min simulation of biological time shows a representative example of merotelic attachment. C) Changes in KT-KT orientation angle (solid lines; left y-axis) and distance between KT-pair and equatorial plate (dash-dotted lines, right y-axis) over time for the most representative simulation run for the case study with CH arms present. D) Number of MTs vs. time profiles for the same simulation run as in panel C showing proper amphitelic attachment. Abbreviations used: L1 and L2 denote the numbers of MTs from the left CS attached to the first KT and second KT, respectively; R1 and R2 are the numbers of MTs from the right CS attached to the first KT and second KTs, respectively. E) Probability to find each type of attachment for a rigid corona surface (blue) and for a flexible corona surface (red). The statistics are based on n = 20 simulation runs for each case study. F) KT-MT interface and corona flexing for a case of amphitelic attachment. G) KT-MT interface and corona flexing for a case of merotelic attachment.

https://doi.org/10.1371/journal.pcbi.1010165.g006

Small KT deformations do not significantly affect the accuracy of spindle assembly

The KT outer layer is likely not rigid but flexible, and it was hypothesized that this flexibility could play a positive role in correction of merotelic attachments [6]. If a ‘wrong’ MT connects to an otherwise amphitelically attached KT, this KT is already pulled by multiple amphitelic connections away from the AB cloud. However, the merotelic MT could then deform a part of the KT back into the AB cloud, which would accelerate the detachment of this MT. Therefore, we examined whether there was any effect of the KT flexibility on the number and types of KT-MT attachments. For this case study, we carried out simulations of a naked CH without CH arms. All parameters are given in Tables 1, 2, 3 and 4 and Table A in S1 Text. In this case study, we compared the results obtained for a rigid KT with the results for a soft KT. To model a soft KT, we decreased 40-fold the value of K_KT,r in Eq 4 for the potential energy going from a value of K_KT,r = 3.3×10³ pN/nm (rigid KT) to 82.5 pN/nm (soft KT). Also, we changed the cut-off radius R_Ndc−bond for the KT structure (with χ = 0.5 curvature) by decreasing it from 200 nm to 50 nm (see Fig 2A and 2B).

We carried out 20 independent simulation runs for 30 min of biological time both for rigid and deformable KTs. The results obtained for final KT-MT attachments observed at the end of simulations resulted in the statistics of KT-MT attachments shown in Fig 6E. We see that there is little difference between the results obtained for the rigid KT versus deformable KT. The reason is likely that the deformations brought about by forces exerted by the attached individual MTs are not large enough to exert the desired effect (see Fig 6F and 6G for the amphitelic and merotelic cases; see also S3 Movie for the amphitelic case).

Thermal noise worsens accuracy of KT-MT attachments for naked centromeres but improves accuracy of CH possessing CH arms

Thermal noise often plays important roles in cells, helping to achieve robust biological functions [67–75]. For example, higher levels of stochastic noise result in the increased robustness of cell polarization [76]. Previous research has not addressed the role of thermal noise in mitosis, and so here we explored this role computationally by switching on and off the random force exerted on each mechanically active component. In all case studies reported above, the random force modeling stochastic cell environment (σg_i(t) term in Eq 10) was switched on. For this case study, we performed independent simulations, which include: i) 8 trajectories for a naked CH with P:AB ratio = 1:100, 8 trajectories for a naked CH with P:AB ratio = 1:10, and 8 trajectories for a naked CH with P:AB ratio = 1:1; ii) 20 trajectories for a single naked CH with AA present; iii) 20 trajectories for a single CH with CH arms; and iv) 20 trajectories for a naked CH with flexible KT. The results of simulations with and without thermal noise are displayed in Fig 7.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Influence of stochastic noise on dynamics of KT-MT attachments.

Comparison of the types and numbers of KT-MT attachment from the simulations with the random force component (shaded bars) and without random force component (blank bars) for: A) A single KT pair with the Phosphatase to Aurora B (P:AB) ratio = 1:10, B) A single CH with CH arms (see also Fig J in S1 Text). The KT-KT orientation angle vs. time (solid lines) and distance to the equatorial plate vs. time profiles (dashed lines) for the case of a single KT pair with P:AB ratio = 1:10 (panel C) and for a single CH with CH arms (panel D). Black lines in panel C corresponds to simulation run without thermal fluctuations and demonstrate amphitelic attachment. Red and blue lines correspond to simulation runs with thermal fluctuations and demonstrate monotelic and merotelic attachments, respectively. In panel D, black and green lines represent simulation runs without and with thermal fluctuations, respectively. Panels E and F demonstrate the profiles for the number of MTs attached to KTs from both CS for a single KT pair with P:AB ratio = 1:10 and for a single CH with CH arms, respectively. In panels E and F, the assignment of curve color is same as in panels C and D. Abbreviation L1, L2, R1, and R2 are same as in Fig 6D.

https://doi.org/10.1371/journal.pcbi.1010165.g007

Two sets of simulations show significant effects of noise. Specifically, for the case of a naked CH pair with P:AB ratio = 1:10, neglecting the random force component dramatically increases (by 3-fold) the number of amphitelic attachments, while the numbers of merotelic and monotelic attachments decreases 4-fold and 2-fold, respectively. The hint to the reason for this effect came from characteristic dynamics of the monotelic and merotelic attachments (red and blue curves in Fig 7C and 7E, respectively). Let us start with the monotelic ones: for the first ~7 min, a small number of MTs from the ‘right’ pole attach to the ‘upper’ KT (the one ‘exposed’ to both poles). In the first 3 min, these MTs pull rapidly turning the CH, exposing the attached KT further to the right pole and exposing it to more MTs from this pole. Starting from 2 min, these MTs start to rapidly pull the CH away from the equator and closer to the right pole. This proximity and beneficial angle create an avalanche of KT-MT binding between 7 and 15 min (Fig 7E), and by 15 min the CH practically falls on the right pole. (Fig 7A). In the case of the merotelic attachment (blue curves in Fig 7C and 7E) the following events take place: almost instantly, one MT from the right pole and another from the left pole bind to the upper KT exposed to both poles. The MTs pull in the opposite directions, which keep the CH near the equator and not rotated. This keeps the sister KT hidden, and it never makes any KT-MT attachments. Although the number of KT-MT attachments for the opposite poles becomes different, this difference does not significantly change the angle and position of CH.

Next, we compared the movements of the CHs and MTs with and without thermal noise. We found the following: CHs are large, and thermal forces have virtually no effect on them. However, the MTs are thin elastic rods, and thermal forces are capable of bending the MTs, and so the latter start to undergo rapid small-amplitude undulations. What turned out to be important is rapid trembling of the MT plus-ends. It enables the plus-ends to rapidly find a KT when a growing MT passes by a KT. Thus, effectively, this increases the KT-MT attachment rate. At first glance, this is beneficial for rapid assembly. However, this is not beneficial for accuracy! Indeed, either numerous MTs rapidly attach to the upper KT from both poles, while the lower one is hidden, and then the CH never turns, thus keeping the lower KT hidden, or the first MT has time to rapidly turn the CH toward one of the poles, and then the geometric vicious cycle begins pulling the CH closer and closer to that pole and acquiring more and more MTs from that pole ensues. Of course, other scenarios exist (and were observed in the simulations), but these two cases illustrate the fundamental issue.

What improves the situation when the thermal noise is turned off is revealed by the amphitelic dynamics (shown in black in Fig 7C and 7E). The MTs do not tremble, there is a slower accumulation of attachments. Because of that, there is a sizable time interval between consecutive attachments, and the first attachment formed earlier has time to turn the CH before the second attachment from the opposite pole formed later in time preventing the CH turning. This CH turning makes the attached KT less visible to the opposite pole, while the sister KT is now more visible to the opposite pole, and so it is more likely to create an amphitelic attachment.

The surprise was that when we performed a comparison for the case of a CH with CH arms, the result was quite the opposite: thermal noise improved the accuracy. The probability of amphitelic attachment slightly decreased and the probability of merotelic attachment strongly increased (see Fig 7B). Dynamics shown in Fig 7D and 7F gives us a hint why this is the case. Green curves show the amphitelic attachment dynamics for the case with thermal noise. One can see that in this case the number of KT-MT attachments grows rapidly, to the correct poles. Indeed, the attachment rate is high because of the noise effect. The numbers of attachments are also asymmetric–the upper KT attaches to the left pole first, which rapidly turns the centromere, but does not rapidly move the CH. The reason is that due to CH size, it’s effective viscous drag is large. This prevents the geometric vicious cycle from starting. While the CH starts moving away from the equator very slowly, the rapid turn (turning of the small centromere is rapid, which can be seen from the movies; see S2 Movie) allows MTs from the opposite pole to attach to the sister KT, thus establishing the amphitelic connection. The characteristic merotelic connection in the case study without thermal noise is shown by black curves in Fig 7D and 7F. The MTs attach rarely (low attachment rate), CH does not shift much, and in the end the fatal merotelic attachment is formed.

The other two sets of simulations, i.e. with flexible KT and with AA, did not show significant differences between the simulations with and without thermal noise (Fig J in S1 Text).

3D model of search-and-capture process

We developed and tested the molecularly detailed, mechanochemical 3D Stochastic Reaction-Diffusion-Dynamics model (SRDDM) of the Search-and-Capture pathway of mitotic spindle assembly. This model extends beyond the previous modeling attempts [28] by taking into account realistic 3D geometry of CHs with deformable sister KTs of realistic size and shapes the flexible centromeric region, and elastic dynamic MTs. The model also simulates the diffusion of Aurora A and B kinase and incorporates the molecularly explicit kinetics of (de)phosphorylation of Ndc80 linkers, and resulting formation and rupture of Ndc80-mediated KT-MT connections. The SRDDM is based on solving the spatially inhomogeneous RDME for describing biochemical reactions and molecular transport, and LD to treat the force-dependent dynamic mechanical processes (Fig 3). The model combines the stochastic description of kinetic processes (enzyme catalysis, formation and rupture of protein-protein complexes) and molecular transport of chemical species in the cell interior (diffusion of Aurora A and B kinases and Phosphatase) with the stochastic description of force-generating processes, (CH pushing/pulling by growing/shortening MTs).

Several frameworks have been proposed for computational modeling of the cytoskeleton, and we briefly discuss the most relevant. One of the first computational platforms for numerical modeling of the cytoskeleton network was Cytosim [77], which has been very successful and widely used by the scientific community. This model describes the physics using Brownian dynamics and allows one to simulate interactions between flexible fibers and motors. Cytosim has been utilized to model a wide variety of cytoskeletal phenomena, including spindle dynamics. The main thrust of CellDynaMo is to integrate transport and biochemical reaction processes with cytoskeletal mechanics in 3D and to introduce more complex and molecularly detailed heterogeneity into modeling the intracellular environment, i.e. geometrically and mechanically realistic chromosomes and kinetochores. This builds on the Cytosim success and extends it further. An example of a more recent modeling environment is MEDYAN [78], a sophisticated model of the actomyosin network that includes a large number of molecular components and the cell membrane. As in CellDynaMo, in MEDYAN Langevin Dynamics is used to describe the force-dependent dynamic mechanical processes, and the Gillespie algorithm is implemented to treat the biochemical kinetics. MEDYAN is capable of simulating biological systems with intricate membrane-cytoskeleton interactions [79]. In Lattice Microbes [52], as in CellDynaMo, the next-subvolume extension [51] of the Gillespie algorithm [59,60] is used to sample the RDME. Odell and Foe developed a detailed 3D model of the spindle-related assemblies [23] based on a combination of mechanical interactions and transport equations. Another model is aLENS developed by Yan et al. [80]. This 3D model describes cytoskeletal assemblies and also includes MTs and crosslinkers. The aLENS employs the Brownian motion of MT filaments (with excluded volume interactions) and the kinetics of binding and dissociation of crosslinkers (Kinesin-5 and dynein motor proteins) described with the Monte Carlo algorithm. The aLENS model takes advantage of parallel simulations with MPI and OpenMP. Yet another cytoskeleton model is CyLaKS developed by Fiorenza et al. [81]. These authors also use the kinetic Monte Carlo algorithm to describe binding of motors and crosslinkers (kinesin-1 and PRC1) to MT filaments and Brownian Dynamics to update the filaments’ positions. In the CyLaKS model, the mechanochemistry (cross-communication between the kinetic and dynamic components) is realized as in CellDynaMo. A model of the cytoskeletal network developed by Belmonte et al. [82] uses the 2D Brownian Dynamics of MT filaments, and stochastic kinetics to account for binding and dissociation of molecular motors (myosins, kinesins, and dyneins) to and from the MT filaments.

CellDynaMo package

We mapped the SRDDM into the CellDynaMo package. The initial input for computer simulations with CellDynaMo is provided by a user in a set of configuration files, which specify the kinetic parameters, e.g. reaction rate constants, number of copies, the force field parameters, e.g. stretching and bending rigidities of MTs and CH arms, strength of KT-MT attachments, and the cell morphology parameters, e.g. number of chromosomes, number of kinetochore pairs, membrane shape, etc. These parameters specify the geometry of spatial arrangements of various cell components (Fig 1), biochemical kinetics (Tables 2 and 3 and 4), diffusion of biomolecular species (Table 3), force-generating properties, and mechanical interactions (Table 1) for cell components. For example, in the configuration file for cell morphology one can describe the cell geometry (shape and size), the number of chromosomes, the length of chromosome arms, the shape and size of kinetochore corona surface, the number of microtubules per centrosome, the initial configuration of all cell components, and the boundary conditions (cell membrane). Information about the viscosity and temperature of the cell cytoplasm is contained in the distribution of random forces. In the kinetics configuration file, the user can specify the reaction network, copy numbers of biomolecules, and kinetic rate constants for all biochemical reactions. In the mechanics configuration file, the user can specify the properties for all force-generating components, including stretching and bending rigidities, strength of KT-MT attachment, etc. A full list of parameters used in CellDynaMo to specify cell morphology, and to solve numerically the RDME and LD equations is provided in Table A in S1 Text.

Numerical algorithms for sampling RDME and LD are organized into the following workflow. For each molecular component in every subcell, the RDME module determines if the component diffuses out of that subcell, records what neighboring subcell it diffuses into, and performs this check for all components. Next, for each subcell the RDME module checks if any reaction occurs, performs this check for all reactions and updates the system state. Biological time is determined by the number of steps of the RDME algorithm. Next, the algorithm checks if there are changes in the chemical state (see Table 4) that might lead to changes in the mechanical state (Fig 3). The mechanochemical coupling is realized through the particle-particle interactions (e.g. KT-MT attachment/detachment). When the reaction and diffusion events lead to changes in the mechanical state (e.g. formation/dissociation of KT-MT bond, MT pulling or pushing on chromosomes), the RDME module switches off. At this time point, the LD switches on and brings the cell to a new state of mechanical equilibrium. For example, if a dissociation of an MT from the KT occurs, this information is passed to the LD module which switches off the excluded volume interactions. When a new state of mechanical equilibrium is reached, the LD module switches off, and the RDME module is back on (similar approach is implemented in MEDYAN [78]). CellDynaMo generates a simulation output, which includes coordinate and force files for mechanically active components, and file with subcell specific content (e.g., copy numbers of chemical species in each subcell).

The results of application of CellDynaMo show that contemporary GPUs can be utilized to generate multiple trajectories of subcellular dynamics on minute-to-hour timescales for a large system of ~10⁶ particles treated implicitly (RDME) and ~10⁴ particles described explicitly (LD) using the one-run-per-GPU approach [62]. These are the experimentally relevant timescales for many biological processes in eukaryotic cells, including mitotic spindle assembly. The 8–12 GB GPU on-board memory is adequate to describe cell dynamic processes involving ~10⁶ particles. For example, when running CellDynaMo on a single contemporary GPU (GeForce GTX 1080), it takes ~70 hours of wall-clock time to generate a single trajectory over 10⁶ steps of the RDME algorithm, which translates to ~30 min of biological time. It takes 3–5 days of computational time to probe the importance of Aurora A enzyme for proper biorientation, or the balance of Aurora B and Phosphatase enzymes on the strength of KT-MT linkages.

Correct versus incorrect KT-MT attachments

We employed the CellDynaMo framework to interrogate computationally several important aspects of the mitotic spindle assembly process and to revisit and assess the validity of the proposed Search-and-Capture paradigm [8–11,16,18–20]. The results are by no means complete or systematic: this study was mainly devoted to development and testing of the computational modeling framework, and to gaining a few insights about the workings of the Search-and-Capture process with realistic cell geometry and mechanics, as well as kinetics of the most relevant chemical reactions. Surprisingly, even such a limited exploration offers several unexpected insights. We found that, at best, the geometrically and mechanically realistic Search-and-Capture mechanism can result in 2/3 correct (amphitelic) KT-MT attachments (Figs 4–7). The amphitelic attachments usually have ~10 KT-MT linkages (Figs 4–7), which compares well with experimental observations. Cells in vivo perform much better [83,84] likely due to the processes we have not included yet into the model (see below); but in vitro a similar accuracy was observed [6]. The majority of incorrect attachments are predicted to be merotelic, which agrees with both previous predictions [15] and experimental observations [85]. The reason is that monotelically attached CHs tend to ‘acquire’ more MTs and convert into other KT-MT attachment types. Syntelic attachments, similarly, tend to acquire an MT from the ‘wrong’ pole and become merotelic, or else destabilize and become monotelic again. Merotelic attachments, on the other hand, are hard to correct, and additional attachments do not change their character.

One of the conclusions from this study is intuitively clear: an intermediate rate of KT-MT detachment is optimal for CH connection accuracy. When the MTs detach too rapidly, the amphitelic attachments do not have much chance to stabilize. However, when the MTs do not detach rapidly enough, then incorrect attachments are hard to correct. Our simulations also vividly illustrated the beneficial geometric feedback leading to the characteristic amphitelic attachment scenario: the first attachment rapidly turns the attached KT toward the respective pole, thus shielding this KT partially from the other pole while exposing the other KT more to the other pole. This leads to more attachments biased in the amphitelic way, which improves the spindle angle further. In addition, resulting pulling of the sister KTs apart brings the KT-MT linkers away from the Aurora B cloud, thus reinforcing the correct attachments. On the other hand, one of the greatest ‘dangers’ for a CH that our simulations have revealed is the geometric vicious cycle that occurs when an increasing number of MTs attach to a KT pull this KT toward a pole, thereby exposing it to even larger number of MTs and the other KT to smaller number of MTs. One of the unexpected insights from our study is that CH arms break this vicious cycle by making the CH ‘bulky’ and not allowing it to move rapidly away from the spindle equator, while also allowing both KTs to have a chance to acquire connections from the opposite poles.

Another nontrivial insight from this study is that thermal forces due to the stochastic environment in the cell cytoplasm are important despite bulky CHs: thermal trembling of the MT plus-ends effectively accelerates the KT capture, which poses no danger in the presence of flexible CH arms. It was shown that thermal fluctuations drive random pivoting of microtubules [86], so that the plus-ends of 1-2-μm long microtubules constrained at the minus-end are capable of spanning the 30-degree angle, and, as a result, accelerate the process of kinetochore capturing to just ~3–4 min [86]. In the SRDDM, the 8-μm long microtubules are relatively stiff, and changes in the bending angle do not exceed 2 degrees. However, we found that even this reduced microtubule flexibility is sufficient to enable kinetochore capturing within the first few minutes of computational experiment, in agreement with the effect posited in [86].

We found that the presence of Aurora A can improve the accuracy of formation of KT-MT attachments by breaking the incorrect attachments when a CH is too close to one of the poles, where Aurora A phosphorylates Ndc80 linkers, thus increasing the KT-MT detachment rates and reversing the CH convergence toward the pole. However, at this stage of development, the model is not able to completely reverse the CH movement back to the center of the spindle. This is because the growing MTs impinging on the CH arms bend and slide off the arms rapidly, not exerting any significant polymerization force. This is another valuable insight from our simulations: MT pushing is ineffective and chromokinesins, not yet incorporated into the model, are necessary to generate a pushing force. Lastly, we found that mechanical flexibility of the KTs has little effect on the statistics of KT-MT attachments. However, so far, we tested only relatively stiff KTs that deform only slightly. The logic of hypothesized effects of KT deformation is based on producing large KT deformations, which our simulations cannot currently handle.

A number of important questions remain for the future, including the effects of: 1) excluded volume obstruction of the cellular space by the CH arms, and 2) key molecular motors, such as chromokinesins, kinesin-5, and dynein. In particular, progress in this direction will allow us to explore the ‘anti-poleward wind’ concept and formation of the spindle interpolar bundles [87–89], not to mention the mechanisms proposed as alternatives to Search-and-Capture. A pioneering study [90] paved the way for understanding the roles of molecular motors in spindle maintenance, examining conditions under which antiparallel microtubules overlapping can stabilize and maintain a finite distance between the spindle poles [90]. As the spindle assembly process is complex enough, in this study we fixed the spindle pole positions. In the future, using SRDDM, we will be able to interrogate the role of molecular motors in full three-dimensional spindle evolution. Recent data suggests that KTs do not become captured by the end-on MT attachments right away; rather, long centrosomal MTs interact with CHs laterally by bringing CHs to the spindle surface, where the lateral attachments are converted into the end-on ones [91]. Utilizing the CellDynaMo package will enable us to directly simulate this phenomenon. When the formation of the spindle surface is disrupted, the assembly likely reverts to the Search-and-Capture mechanism, and so a thorough investigation of this mechanism in future reports is still valuable.