1. Introduction

DNA segregation is a complicated process that is critical for cell proliferation and survival [16,30]. Failures during segregation can result in aberrant DNA contents (aneuploidy), a phenomenon which is prevalent in human cancers [24,28]. The fidelity of chromosome segregation during cell division is monitored by control mechanisms called checkpoints, which ensure that particular criteria are met before moving on irreversibly to the next phase [30].

In mitosis, the Spindle Assembly Checkpoint (SAC; [33]) ensures that all chromosomes are properly attached to spindle microtubules via their kinetochores. Even a single unattached or misattached chromosome is sufficient to keep the checkpoint active and engaged [32,31]. (A human mitotic cell has 46 chromosomes and 92 kinetochores.) In budding yeast (Saccharomyces cerevisiae and Drosophila male), an additional control checkpoint exists to place the correct DNA into the right cell during asymmetric cell division. This regulation is known as the Spindle Position Checkpoint (SPOC), and delays mitotic progression until the spindle is correctly aligned along with the polarity axis [2].

Both SAC and SPOC have prominent similarities, although they constitute different mitotic checkpoints. They broadcast a 'wait' signal to the environment, and rely on turnover of the inhibitor and activator at an organelle (kinetochore for SAC and spindle pole body in SPOC). The SAC integrates signaling information about attachment of the individual kinetochores, which broadcast a 'wait'-signal unless a correct attachment is established. Many core SAC components like Mad2 are recruited to unattached kinetochores, and broadcast a nucleoplasm 'wait'-signal and inhibit Cdc20, the APC/C activator. Upon kinetochore-microtubule attachment, these components are rapidly removed from the kinetochores and APC/C:Cdc20 formation (SAC silencing) is turned on (see Fig. 1A). Similarly, central SPOC components Bf1:Bub2 are localized and regulated at the spindle pole bodies (SPBs), which are broadcasted through out the cytosol and inhibit the downstream pathway, Tem1. Signaling from the SPBs is shut down after correct spindle alignment with the polarity axis is achieved (see Fig. 1B).

The checkpoint mechanisms, SAC and SPOC, are hard to observe experimentally in living cells, due to a number of technical challenges. For example, even a low number of components can have various localizations and states upon which the interactions depend. Another issue is that the average diameter of many proteins at the kinetochore is about 40Å, making connections between them invisible to current microscopy techniques. Likewise, the SPOC protein Tem1 activity, whether it is GTP-or GDP-bound, and Bfa1 protein phosphorylation by various kinases (or lack thereof), are not observable in the wet lab. These limitations can be addressed by employing mathematical models and numerical simulations, which can improve our understanding of the mechanics of cell division. However, mathematical methods can be hindered by combinatorial explosion in the amount of intermediate components (complexes) and explicit representations. Also, the different components often interact nonlinearly in time and space; in the presence of various feedback loops, these interactions lead to phenomena that are difficult to predict [10,0,25,11,38]. A combination of experimental work and rigorous mathematical models was central to exclude some hypothesized SAC checkpoint architectures and to elucidate how the SPOC's elaborate system functions (e.g. [4]). All mathematical models extant in the literature studied either human SAC or yeast SPOC activation at a detailed molecular level or in abstract models to distinguish between different pathways [7,34,14,13,15,21,26,19,9,0,18,23,27,29]. However, none of these models addressed SAC or SPOC silencing. Moreover, no rigorous mathematical analysis of properties such as bifurcation, stability, or the existence of feedback loops has been performed for either SAC or SPOC to date. The research groups of Novak and Tyson work intensively on cell cycle-related modeling for yeast. Their recent small model consists of five reactions, five species, and two ODEs based on Michaelis-Menten kinetics with double-negative and two double-positive feedback loops[39]. Also, the smallest chemical reaction system with bistability in the literature contains four reactions, and two ODEs based on mass-action kinetics [41]. It has double positive and a single negative feedback loops. Unfortunately, no regeneration for the reactants is possible in this model.

The purpose of this paper is fourfold. First, a minimal bistable SAC model for activation and silencing was constructed. The model is based on mass-action kinetics and comprises four reactions, double-negative and two double-positive feedback loops, and two ODEs. It is structurally fully distinct to the known smallest biochemical model [41] and structurally comparable to the yeast mitotic model [39,17]. Second, the same model structure was applied to the SPOC, with both SPBs included. Subsequently, a one-parameter bifurcation was computed for these models, in order to demonstrate the realistic biochemical switches. Eventually, numerical simulations were carried out for the system as a system of ordinary differential equations (ODEs), and also as partial differential equations (PDEs; reaction-diffusion systems) with various parameters.

---

2. Methods

---

2.1. Reaction equations

The reaction rules governing the Spindle Assembly Checkpoint (SAC) system are (cf. Fig. 2A):

$\begin{align} &\text{Cdc20:C-Mad2}\xrightarrow{ k_{1}, \text{Unattached Kin }} \text{MCC} \label{eq:Mad2} \end{align}$

(1)

$\begin{align} &\text{MCC} + \text{APC/C} \mathop \rightleftharpoons \limits_{k_{-2}}^{k_2} \text{MCC:APC/C} \label{eq:MCCAPC} \end{align}$

(2)

$\begin{align} &\text{MCC:APC/C } \xrightarrow{k_{3}, \text{Attached Kin}, \text{APC/C}} \text{Cdc20:C-Mad2} + \text{APC/C}\label{eq:stripingMAD2} \end{align}$

(3)

The Spindle Position Checkpoint mechanism (SPOC) has very similar reaction rules (see below), but differs significantly in the initial concentrations, rate contacts, and signals of exactly two SPBs (cf. Fig. 4A):

$\begin{align} &\text{Bub2} \xrightarrow{ k_{1}, \text{misaligned SPB }} \text{Bfa1:Bub2} \label{eq:Bub2} \end{align}$

(4)

$\begin{align} &\text{Bfa1:bub2} + \text{Tem1} \mathop \rightleftharpoons \limits_{k_{-2}}^{ k_2 } \text{Bfa1:Bub2:Tem1} \label{eq:Bfa1Bub2} \end{align}$

(5)

$\begin{align} &\text{Bfa1:Bub2:Tem1 } \xrightarrow{k_{3}, \text{aligned SPB}, \text{Tem1}} \text{Bub2} + \text{Tem1} \label{eq:stripingBub2} \end{align}$

(6)

---

2.2. Ordinary differential equations model

By applying the law of mass-action kinetics, the reaction rules (Eqs.(1-3) and Eqs.4-6) can be translated into sets of time-dependent nonlinear ordinary differential equations (ODEs). The translation is done by computing $\textrm{d}S/\textrm{d}t = Nv(S)$ with state vector $S$, flux vector $v(S)$ and stoichiometric matrix $N$. The initial amounts for reaction species and rate constants are taken from the literature (cf. Table 1).

---

2.3. Partial differential equations model

Adding a diffusion term as a second spatial derivative transforms the system into one of coupled partial differential equations (PDEs), which is known as a reaction-diffusion system and has the following general form:

$\frac{\partial {[C_{i}]}}{\partial t} =\underbrace{\; {\text{D}_{i}}\nabla^2 {[C_{i}]\;}}_{\text{Diffusion}}+\underbrace{\;R_{j}\left(\{[C_{i}]\}; \text P\right)\;}_{\text{Reaction}}, \label{eq:RDS}$

(7)

where ${[C_{i}]}$ refers to the concentration of species $i = \{1,\ldots{},4\}$. The first term on the right hand side represents the diffusion and the second one represents the biochemical reactions $R_j = \{R_1,...,R_4\}$ where species $i$ is involved. The constant ${\text{D}_{i}}$ denotes per-species diffusion for species $i$. $t$ specifies time and $P$ symbolizes phenomenological parameters. The operator $\nabla$ refers to the spatial gradient in spherical coordinates $\left(\nabla\left(r, \theta,\varphi\right) = \frac{\partial}{\partial r}\vec{e_r} + \frac{1}{r}\frac{\partial}{\partial \theta}\vec{e_\theta} + \frac{1}{r\sin(\theta)}\frac{\partial}{\partial \varphi}\vec{e_\varphi}\right)$. Recent studies have shown that spherical symmetry is suitable and sufficient for use when modeling the SAC mechanism [21,12,13]; thus, we use the spherically symmetric form. Eq. (7) reduces to the following system of PDEs, which depend on $t$ and $r$:

$\frac{\partial {[C_{i}]}}{\partial t} = \underbrace{\; \frac{D_{i}}{r^2} \frac{\partial}{\partial r} (r^2 \frac{\partial {[C_{i}]}}{\partial r}) }_{\text{Diffusion}}+ \underbrace{\;R_{j}\left(\{[C_{i}]\};\text P\right)\;}_{\text{Reaction}}. \label{eq:RDSsymm}$

(8)

---

2.4. Numerical simulations of the ODEs model

The systems of ODEs were implemented in the freely-available software package XPPAUT[8], and integrated using the Rosenbrock method (stiff solver). The bifurcation analyzes and the related numerical integrations were conducted with AUTO [6] via an XPPAUT interface.

---

2.5. Numerical simulations of the reaction-diffusion PDEs model

For the spatial simulations, the mitotic cell is considered as a 3-sphere with radius $R$. The last unattached kinetochore is a small 3-sphere with radius $r$, which is located in the center of the cell. All boundary conditions are assumed to be reflective and the numbers of all interacting elements to be conserved. Additionally, all PDEs are assumed to be spherically symmetric (solely as a function of $t$ and $r$, cf. Eq. 8). All reactions are assumed to follow the mass-action-kinetics law. The kinetic rate constants are taken from the literature, as given in Table 1.

The reaction-diffusion system of PDEs were solved numerically using MATLAB (MathWorks), and integrated using its predefined function called pdepe-solver, which solves systems of parabolic and elliptic PDEs in one space variable $r$ and time $t$.

The pdepe-solver converts the PDEs to ODEs using a second-order accurate spatial discretization based on a fixed set of user-specified nodes [36]. This is done using piecewise non-linear Galerkin (regular case) and implicit Petrov-Galerkin (singular case, second-order accurate). The ordinary differential equations resulting from discretization in space are integrated via the multistep solver ode15s which is a variable order solver based on the numerical differentiation formulas (NDFs)using the numerical Gear method [35]. To check that our results were not influenced by the spatial discretization method used in pdepe, we repeated all simulations for 50,100 and 1000 grid cells.

---

3. Results

---

3.1. Mathematical framework of SAC

The wiring diagram of the SAC mechanism (Fig. 2A) was translated into a set of reaction equations (see Methods), which were then translated into a set of coupled ordinary differential equations (ODEs) under the assumption of mass action kinetics for all reactions. It is clear biochemically that the total concentration of Cdc20:Mad2 is constant in the system and can be expressed as [Cdc0:CMad2T] = [Cdc20:CMad2] + [MCC] + [MCC:APC/C]. The same is true for total APC/C, thus [APC/CT] = [APC/C] + [MCC:APC/C]. Also, for simplicity, the total amount of MCC was defined as: [MCCT]=[MCC]+[MCC: APC/C]. Under these assumptions, the reduced system can be easily written as the following nonlinear ODEs:

$\frac{\text{d[MCCT]}}{\text{d}t} = k_1. U(\text{[MadT]}-\text{[MCCT]})-k_{3}.A\text{[APC][MCC:APC/C]}) \label{eq:conmcc}$

(9)

$\frac{\text{d[APC/C]}}{\text{d}t} = -k_2\text{[MCC][APC/C]}+ (k_{-2} + A\text{[APC]}) \text{[MCC:APC/C] } \label{eq:conapc}$

(10)

The parameter $U$ refers to the number of unattached kinetochores, which vary from 92 to 0, while $A$ references the number of attached kinetochores and is defined by $A=92-U$. APC/C can be presumed in steady state to analyze the bifurcation for the kinetochore signal versus the total MCC. Also, plugging the constant concentration [APC/CT] into the system 9-10 produces

$\begin{align} &\frac{\text{d[MCCT]}}{\text{d}t}\nonumber\\ = &k_1. U(\text{[MadT]}-\text{[MCCT]})-k_{3}.A(\text{[APC/CT]-[MCC:APC/C]})\text{[MCC:APC/C]}) \label{eq:conmcc_b} \end{align}$

(11)

$a.\text{[MCC:APC/C]}^2 + b.\text{[MCC:APC/C]} + c = 0 \label{eq:conapc_b}$

(12)

where the parameters $a,b$, and $c$ are defined as follows: $a = -k_2 - A$, $b = k_2\text{[MCC]} + k_2\text{[APC/CT]} + k_{-2} + A.\text{[APC/CT]}$ and $c = -k2 \text{[MCC]}\text{[APC/CT]}$.

First, one-parameter bifurcation analysis was performed for the nonlinear system (Eq. 11 and 16). The aim is to demonstrate the bistable switch states influencing total MCC, while kinetochores are gradually attached. The simulations were conducted using AUTO software (see Methods). The results in Fig. 2B display a typical S-shape, representing the number of attached kinetochores versus the total concentration of the MCC inhibitor. The stable node points for steady states are shown in solid lines and the unstable saddle points are depicted with dashed lines. Stable and unstable steady states meet at saddle-node bifurcation points, which are indicated by solid circles. At the attached kinetochore number (91.98, i.e. nearly all), the SAC checkpoint switched off and APC/C activated rapidly. The total MCC lowers back to zero as the cell enters anaphase. The switch flipping from the SAC-active state to the SAC-inactive state as the number of attached kinetochores raises is indicated by a black dashed line. Slight movements of the bifurcation curve to the left or right are possible, depending on the values of the parameters $k_3$ and $k_2$. However, the earliest possible shift (to the left) can take place with high values of both $k_3$ (1.1) and $k_1$ (0.5).

Additionally, the dynamics as the change of concentrations over time were simulated and plotted (Fig. 2C) using XPPAUT (see Methods). The concentration of the APC/C component (Fig. 2C pink line) remains very low, as long as no kinetochores are attached. After approximately 25 minutes, the APC/C activity increases quickly to reach its maximum. This result is consistent with the experimental findings reported in the literature [1,27]. The inhibitor complexes MCC:APC/C that sequester APC/C display exactly the opposite behaviors compared to APC/C (Fig. 2C, brown line). MCC concentration behaves similarly to MCC:APC/C with respect to the difference in the initial amount (Fig. 2C, blue line).

Also interesting is the effect of the diffusion coefficient on the SAC system, particularly since the MCC diffusion constant is unknown. To investigate this, a second-derivative diffusion term was added to the original system (Eq.10-Eq.9). The resulting system of coupled PDEs is known as a reaction-diffusion system:

$\begin{align} &\frac{\text{d[MCCT]}}{\text{d}t}\nonumber\\ =& \frac{D}{r^2} \frac{\partial}{\partial r} (r^2 \frac{\partial {[MCCT]}}{\partial r}) +k_1. U(\text{[MadT]}-\text{[MCCT]})-k_{3}.A\text{[APC][MCC:APC/C]}) \label{eq:mccpde} \end{align}$

(13)

$\begin{align} &\frac{\text{d[APC/C]}}{\text{d}t}\nonumber\\ =& \frac{D}{r^2} \frac{\partial}{\partial r} (r^2 \frac{\partial {[APC/C]}}{\partial r}) -k_2\text{[MCC][APC/C]}+ (k_{-2} + A\text{[APC]}) \text{[MCC:APC/C] } \label{eq:apcpde} \end{align}$

(14)

This system was subjected to the initial conditions given in Table 1, with reflective (Neumann) boundary conditions and equivalent geometry details for the cell and kinetochore as specified in Methods. The reaction-diffusion system (Eq.13-Eq.14) was implemented in MATLAB and simulated with various diffusion constant values (Table 1). No qualitative changes were recorded using a wide range of diffusion coefficients for MCC from the literature and higher; all behave similarly as shown in the typical curves in Fig. 3A. To ensure that the assumptions used to reduce the model had no influence of this result, a full system including three PDEs was re-simulated using various diffusion values for MCC. Again, no effects were observed (Fig. 3B). We conclude that diffusion has no major influence on the SAC model, and that the use of ODEs is in principle sufficient. This is certainly not applicable to other model structures, and cannot be generalized, particularly to high-dimensional systems or using low diffusion constants.

---

3.2. Mathematical framework of SPOC

Following the same steps as for SAC analysis as in the previous section, the SPOC wiring diagram (Fig. 4A) can be translated into a set of coupled ODEs. Accordingly, the total concentrations of Bub2, Tem1 and Bfa1:Bub2 can be expressed as: [Bub2T] = [Bub2] + [Tem1] + [Bfa1:Bub2:Tem1], [Tem1T] = [Tem1] + [Bfa1:Bub2:Tem1], and [Bfa1:Bub2T] = [Bfa1:Bub2]+[Bfa1:Bub2:Tem1]. The SPOC system is governed by the following equations, under the aforementioned assumptions (see also Methods):

$\label{eq:conbfa1_b} \begin{split} \frac{\text{d[Bfa1:Bub2T]}}{\text{d}t} =& k_1. X(\text{[Bub2T]}-\text{[Bfa1:Bub2]})-\\ &k_{3}.Y(\text{[Tem1T]-[Bfa1:Bub2:Tem1]})\text{[Bfa1:Bub2:Tem1T]} \end{split}$

(15)

$a.\text{[Bfa1:Bub2:Tem1T]}^2 + b.\text{[Bfa1:Bub2:Tem1T]} + c = 0$

(16)

where the parameters $a,b$, and $c$ have the following meaning: $a = -k_2 - X$, $b = k_2\text{[Bfa1:Bub2]} + k_2\text{[Tem1T]} + k_{-2} + Y.\text{[Tem1T]}$ and $c = -k2 \text{[Bfa1:Bub2]}\text{[Tem1T]}$. $X$ refers to the number of misaligned SPBs (2 to 0), and $Y$ represents the number of aligned SPBs ($Y=2-X$).

Again, AUTO software was used to find the bifurcation curve (see Methods for details). The bifurcation curve (Fig. 4B) is shown for the number of misaligned SPBs versus the total concentration of the Bfa1:Bub2:Tem1 inhibitor. The system switches its bistable state at the value 1.99; subsequently, SPOC is turned off, and eventually Tem1 is rapidly activated.

The numerical simulations of the ODEs is depicted in Fig. 4C, using a stiff solver in XPPAUT. Tem1 (Fig. 2C pink line) is inactive until both SPBs are correctly aligned, which takes place after about 7 minutes. The inhibitor complexes Bfa1:Bub2 and also Bfa1:Bub2:Tem1 have asymmetric behavior compared to Tem1 (Fig. 4C, gray and brown lines).

---

4. Discussion

In eukaryotic cells, the mitotic control prevents DNA missegregation and aneuploidy. The evolutionarily conserved SAC mechanism guarantees that each chromosome has established its attachment to the spindle apparatus before commencing sister-chromatid separation, while the SPOC mechanism assures correct spindle alignment in some asymmetric cell divisions. The complexity of the mitotic control system arises from its fundamental spatial feature. In SAC, a single unattached or incorrectly attached kinetochore (out of 92) has to inhibit all APC/Cs of the cell and solely after last proper attachment the inhibitor has to be switched off rapidly. This behavior likely implies of a feedback loop contribution. The same is applied for SPOC mechanism with the distinguish in the signal that represents two SPBs.

Mathematical modeling can help to improve our molecular-level understanding of the interplay of SAC as well as SPOC components, allowing an understanding of the requirements that the system has to meet. To that end, a mathematical framework containing all signals in SAC (or SPOC) was studied. The network models were constructed based on biochemical reaction rules, and spatial characteristics, such as diffusion coefficients or cell size. Kinetochores and SPBs act as sensory-driven signals in SAC, and SPOC regulations, respectively. The simulation results (as ODEs or PDEs) were able to capture the desired behavior of both SAC and SPOC. Additionally to the simulations, a crucial feedback loop was built around the core components APC/C and MCC (or Tem1 and Bfa1 in SPOC) and their interplay. A one-parameter bifurcation diagram clearly depicts the bistability of the system and the realistic switch from active to inactive control states.

The presented mathematical models can be extended in future work to include various cell cycle checkpoints. Additionally, this approach will serve as a basis for designing experiments and evaluating novel hypotheses related to mitosis and cell cycle. The results provide systems-level details into the DNA segregation control mechanism and demonstrate that the combination of mathematical analysis with experimental data constitutes a powerful tool for investigation of complex biomedical systems.

---

Acknowledgments

The author gratefully acknowledges the visiting fund of the Institute for Numerical Simulation (INS) at Bonn University.

---