CellDynaMo–stochastic reaction-diffusion-dynamics model: Application to search-and-capture process of mitotic spindle assembly

We introduce a Stochastic Reaction-Diffusion-Dynamics Model (SRDDM) for simulations of cellular mechanochemical processes with high spatial and temporal resolution. The SRDDM is mapped into the CellDynaMo package, which couples the spatially inhomogeneous reaction-diffusion master equation to account for biochemical reactions and molecular transport within the Langevin Dynamics (LD) framework to describe dynamic mechanical processes. This computational infrastructure allows the simulation of hours of molecular machine dynamics in reasonable wall-clock time. We apply SRDDM to test performance of the Search-and-Capture of mitotic spindle assembly by simulating, in three spatial dimensions, dynamic instability of elastic microtubules anchored in two centrosomes, movement and deformations of geometrically realistic centromeres with flexible kinetochores and chromosome arms. Furthermore, the SRDDM describes the mechanics and kinetics of Ndc80 linkers mediating transient attachments of microtubules to the chromosomal kinetochores. The rates of these attachments and detachments depend upon phosphorylation states of the Ndc80 linkers, which are regulated in the model by explicitly accounting for the reactions of Aurora A and B kinase enzymes undergoing restricted diffusion. We find that there is an optimal rate of microtubule-kinetochore detachments which maximizes the accuracy of the chromosome connections, that adding chromosome arms to kinetochores improve the accuracy by slowing down chromosome movements, that Aurora A and kinetochore deformations have a small positive effect on the attachment accuracy, and that thermal fluctuations of the microtubules increase the rates of kinetochore capture and also improve the accuracy of spindle assembly.


{
= 13 cos − 23 cos = 0 = − 13 sin − 23 sin = 0 (1) In Eqs 1 above, and are the − and − components, respectively, of the total force applied to bead 3 (tagged), 13 and 23 are the restoring forces with which the first and second springs pull on bead 3, is the angle formed by beads 2, 1, and 3, and is the angle formed by beads 1, 2, and 3 (see Fig F-A). The forces and displacements are related through Hooke's law: 13 = 13 Δ 13 and 23 = 23 Δ 23 (2) where 13 and 23 are the spring constants of the springs which connect beads 1 and 3, and beads 2 and 3. Δ 13 and Δ 23 are distances between beads 1 and 3, and beads 2 and 3, respectively. By substituting the expressions for 13 and 23 from Eqs 2 into Eqs 1, we obtain the exact expressions for the beads' displacements:

Dependence of numerical error on integration timestep:
To investigate the dependence of the numerical error on the integration timestep = 5×10 -2 -5×10 2 ps, we calculated and compared for the three-bead system ( Fig F) the asymptotic values of the particle displacements obtained for the applied 50-pN pulling force with the CellDynaMo implementation (∆ ) and using the exact analytical expression (∆ ) for Δ 13 and Δ 23 (Fig F). The relative error, (∆ ) = |∆ −∆ | ∆ was found to be very low (<1.8%) in the 5×10 -2 -5×10 2 ps range of the integration timestep (Fig F-B).

Dependence of numerical error on solution viscosity:
To investigate the dependence of numerical error on the cytoplasmic viscosity = 1, 5, and 10 cPs and on the integration timestep = 5×10 -2 -5×10 2 ps, we calculated and compared, for a large system of = 100 Brownian oscillators (with the spring constant = 1 pN/nm), the numerical solution of the average particle displacements obtained with CellDynaMo (〈∆ ( )〉 ) and the exact analytical expression (〈∆ ( )〉 ). The exact expression is given by the formula 〈∆ ( )〉 = (1 − − / ), where = 6 is the friction coefficient, is the viscosity, and is the particle size (see Fig F). First, we varied the viscosity but fixed the integration timestep = 50 ps. The numerical and analytical results practically collapse on the same curve (Fig F-C), and the relative error, , was found to be very low <0.4% (see the inset to Fig F-C) for all viscosity values. Next, we fixed the viscosity at = 1 cPs and varied the integration timestep over 4 orders of magnitude from 5×10 -2 to 5×10 2 ps. The relative error was found to be very low <0.3% (see Fig F-D), and the average relative error was <0.16% (see the inset to Fig F-D). Hence, our choice of = 50 ps as the timestep for Langevin Dynamics based description of the mechanical and forcedependent processes is reasonable.
Benchmark test for translational and rotational diffusion in Langevin Dynamics: To examine the system's thermal fluctuations, we compared the numerical implementations of the one-dimensional translational diffusion for a cylinder along the -axis (see the inset in Fig G) and one-dimensional rotational diffusion for a cylinder around the -axis (see the inset in Fig G)  , where is the Boltzmann's constant, = 300 K, is the solution viscosity set to be equal to the viscosity of water ( = 1 cPs), and = 4 = 1.450 µm is the length and = 2 = 0.725 µm is the diameter of the centromere. We carried out three 2.5-min independent runs for a single centromere. In these test simulations, the MT dynamics was turned off, and the integration timestep was decreased from 50 ps to 5 ps to make thermal fluctuations more noticeable. The time profiles of 〈∆ 2 〉/2 and 〈∆ 2 〉/2 are displayed in Fig G The obtained asymptotic values of 0.154 µm 2 /s and 0.866 rad 2 /s compare well with the exact analytical solutions for = 0.152 µm 2 /s and = 0.869 rad 2 /s, respectively (see Fig G).

Benchmark test for diffusion component of RDME algorithm:
To assess the accuracy of numerical implementation of the diffusion part of RDME, we placed 10 4 molecules of AB in the center of the cell (central subcell), 0 = 0, at the initial time = 0, as displayed in Fig A and observed spreading of AB particles at later time points = 1, 2, 5, 10, and 20 s. The non-parametric density estimates of the distributions of particles displacements in one dimension, constructed using the numerical output from CellDynaMo simulations, are compared with the theoretical curves of the exact probability distributions of particles' displacements, in Fig H, which shows excellent agreement between the numerical and exact analytical results. In Eq 4, = 7.3×10 7 nm 2 /s is the diffusion constant (Fig H) for AB molecules, which was estimated using the Einstein-Stokes formula, = 6 , where is the viscosity set to be equal to the viscosity of water ( = 1 cPs) and = 2.9 nm is the size of the AB molecule (see Table A). Therefore, the results obtained indicate excellent agreement between the exact description of the Brownian diffusion (Eq 4) and the numerical description of the Brownian diffusion implemented in the CellDynaMo package.
Benchmark tests for kinetics component of RDME algorithm: Consecutive two-step irreversible kinetics is described by the following scheme, 1 → 2 → , with species , and and the reaction rate constants for the first step ( ⟶ ) 1 and the second step ( ⟶ ) 2 . The populations , , and of species , and are described by the following coupled ordinary differential equations: The exact analytical expressions for the time-dependent populations ( ), ( ), and ( ) are the following: Single-step reversible kinetics is described by the scheme ⥨ , with the reaction rate constants for the forward step ( ⟶ ), 1 , and the backward step ( ⟶ ), −1 . The populations and , of states and are described by the following coupled ordinary differential equations: The exact analytical expressions for the time-dependent populations ( ) and ( ) are: Eqs 8, 9, and 10 for the consecutive two-step irreversible kinetics and Eqs 13 and 14 for the single-step reversible kinetics were used in benchmark simulations to compare the exact description of chemical kinetics with the numerical description of kinetics implemented in the CellDynaMo package.

Flexibility of chromosomes:
To model flexible chromosomes, in the CellDynaMo package the chromosome arms are represented as a collection of spherical beads (see Fig 6B in the main text and Fig  C) of 362.5-nm radius. Each pair of beads is connected by a harmonic spring characterized by the stretching rigidity , (Table A), and each triplet of connected beads is described by the bending rigidity , (Table A). We parameterized the chromosome flexibility by setting the stretching rigidity to , = 3.3×10 3 pN/nm (Table A) and by varying the bending rigidity ( , ) in the bending potential (see Eq 5 in the main text). In the test LD simulations, we set , to be equal to 2.5×10 2 kJ/mol·rad 2 (flexible chromosome arms), to 2.5×10 5 kJ/mol·rad 2 (semi-flexible chromosome arms), and to 2.5×10 8 kJ/mol·rad 2 (stiff chromosome arms), and we ran ~1,000 equilibrium LD trajectories for each value of , . Using the simulation output, we constructed the nonparametric density estimates (see Methods in the main text) of the distributions of bending angle Δ , (Δ ), for , = 2.5×10 2 , 2.5×10 5 , and 2.5×10 8 kJ/mol·rad 2 , which are compared in Fig B. Based on published experimental studies, which report the 1-5 degree range of bending-angle fluctuations [1], we selected , = 2.5×10 5 kJ/mol·rad 2 to characterize the flexural rigidity of the chromatids within the same chromosomes.
Bending constraints imposed by cohesin rings: To model cohesin rings constraining the movement of chromosome arms in the CellDynaMo package, we represented the cohesin rings by weak harmonic springs that link flexible chromatid arms together ( Fig 1B in the main text and Fig C-C). Because the chromatid arms are described using the bead-spring representation, these weak harmonic springs connect the corresponding beads in different chromatids ( Fig 1B in the main text, Fig C-C). We considered chromosomes with two different values of the contour lengths: 4 μm and 8 μm. We parametrized the fluctuations in the chromosome width by varying the stretching rigidity for cohesin rings ( ℎ ) in the stretching potential ℎ (Eq 5 in the main text). In the test LD simulations, we set ℎ, to be equal to 10 -3 kJ/(mol·nm 2 ) (soft cohesin rings) and to 10 3 kJ/(mol·nm 2 ) (stiff cohesin rings), and we ran ~1,000 equilibrium LD trajectories for each value of ℎ, . Using the simulation output, we constructed the nonparametric density estimates (see Methods in the main text) of the distribution of the separation distances , ( ), for ℎ, = 10 3 kJ/(mol·nm 2 ) and 10 -3 kJ/(mol·nm 2 ), which are compared in Fig C-A Figures   Fig A. Illustration of the numerical algorithms implemented in CellDynaMo for simulations of molecular transport: A) Brownian diffusion of biochemical species (AB and AA enzymes) is modeled by the exchange of particles between the next-neighbor subcell. For each particle, the probability of diffusion is compared with a (pseudo)random number r generated at each time step. If < the diffusion move is accepted. B) Each subcell is surrounded by 6 next-neighbor subcells. The value of determines the direction (forward, backward, up, down, left, and right) and the next-neighbor subcell for particle transfer.

Fig F. Benchmark tests Langevin Dynamics algorithm implemented in CellDynaMo: A)
The timedependent profiles of the spring extensions Δ 13 (red) and Δ 23 (blue) obtained numerically using CellDynaMo (dashed lines) and analytically using Eq 3 (solid black lines) for = 5 pN. The inset shows the three-bead system with the constrained beads 1 and 2 and tagged bead 3 in the -plane (blue balls), which are connected by the harmonic springs (shown in yellow) with the spring constants 13 and 23 , respectively. The tagged bead 3 is pulled with force . The pulling force and the components and of the restoring force are indicated by black arrows. Also shown are angle formed by beads 2, 1, and 3, and angle formed by beads 1, 2, and 3. B) The time-dependent profiles of the extensions Δ 13 (red) and Δ 23 (blue) obtained numerically using CellDynaMo (dashed lines) and analytically (solid black lines) using Eq 3 for = 50 pN. Relative error of calculation of asymptotic values of Δ 13 (red lines and data points) and Δ 23 (blue lines and data points),

Fig H. Benchmark test for diffusion component of RDME algorithm implemented in CellDynaMo:
A) Nonparametric density estimates (dashed lines), obtained using CellDynaMo and the exact theoretical profiles (thin curves), obtained with Eq 4, of the probability density function of the displacements of Brownian particles in one dimension, ( , | 0 ). The particle is positioned at the origin 0 = 0 at time = 0 are the particles' spreading is compared for the time points = 1, 2, 5, 10, and 20 s. The nonparametric density curves were constructed as described in Ref. [16]. Color denotation is presented on the graph. B) Snapshots of the spatial distributions of a total of 10 4 Brownian particles in three dimensions obtained for different time points, = 1, 2, 5, 10, and 20 s, which show spreading of the Brownian particles (AB molecules) due to free diffusion.

Fig I. Benchmark tests for kinetics component of RDME algorithm implemented in CellDynaMo:
A) Time-dependent profiles of the populations , , and of chemical species , , and obtained numerically with CellDynaMo and analytically using Eqs 5, 6, and 7 for the two-step consecutive irreversible kinetics. The initial conditions (0) = 1, (0) = 0, and (0) = 0, and the reaction rate constants 1 = 1 s -1 and 2 = 2 s -1 were used. B) Time-dependent profiles of the populations and of chemical species and obtained numerically with CellDynaMo and analytically using Eqs 11 and 12 for the single-step reversible kinetics. We used the initial conditions, (0) = 1 and (0) = 0, and the reaction rate constants 1 = 1 s -1 and −1 = 3 s -1 . Color denotation is presented on the graphs. C) Snapshots of the reaction volume for the two-step irreversible kinetics taken at different time points = 0, 0.7, and 5.0 s, which show the reaction progress (transformation of into and formation ofthe end product). D) Snapshots of the reaction volume for the reversible kinetics taken at different time points = 0 and 1.5 s, which shows the presence of both species and in the equilibrium mixture.