Introduction

Cyclin B–Cdk1 is the core regulator of meiotic and mitotic entry. In early embryos and oocytes, where cells can reach hundreds of microns to millimeters in size, Cdk1 activation does not occur uniformly throughout the cytoplasm. Instead, Cdk1 activity propagates as a large-scale cytoplasmic wave following nuclear envelope breakdown (NEB) Bischof et al. (2017); Chang and Ferrell Jr (2013); Puls et al. (2024) (Box 1,Fig. 1). Such waves provide a mechanism for coordinating cell cycle progression across distances that far exceed typical molecular diffusion lengths and therefore require mechanisms beyond simple diffusive equilibration. More broadly, the spatial coordination of biochemical oscillations via inter- and intracellular waves has emerged as a key focus of research, as wave-driven coordination is critical for maintaining proper cellular function and developmental processes Gelens et al. (2014); Beta and Kruse (2017); Deneke and Di Talia (2018).

Figure 1.

Cytoplasmic Cdk1 waves are directly coupled to the cell cortex. The cortex is a dynamic actomyosin network that controls membrane tension, surface deformation, and cell shape Bement et al. (2005); Nalbant et al. (2004); Machacek et al. (2009). Rho GTPases are central regulators of cortical patterning and their spatiotemporal organization underlies a wide range of cortical behaviors Bement et al. (2024); Beta et al. (2023). During meiosis and mitosis, spatial variation in Cdk1 activity is transmitted to the cell cortex through phosphorylation of the RhoGEF Ect2, which reduces its membrane affinity and inhibits Rho GTPase signaling Hara et al. (2006); Niiya et al. (2006); Su et al. (2011). The resulting cortical response manifests as surface contraction waves (SCWs), consisting of a relaxation phase (SCWa) followed by a contraction phase (SCWb), first described by Hara in amphibian eggs Hara (1971) and later characterized in multiple systems Rankin and Kirschner (1997); Bischof et al. (2017).

Box 1. Cyclin-Cdk1 dynamics and mitotic waves

Box 1— figure 1.

Despite extensive experimental characterization, the functional significance of SCWs and their directionality are still poorly understood. During mitosis, SCWs have been proposed to reduce sensitivity to cell shape during division plane positioning Minc and Piel (2012), a role for which directionality appears secondary. In contrast, during meiosis, SCW directionality correlates with polar body emission in some organisms. In starfish oocytes, SCW propagation and collapse have been associated with the site of polar body extrusion Satoh et al. (2013), although subsequent work has shown that cytoplasmic flows generated by cortical contractions are insufficient to reposition the spindle and that polar body emission can still occur when SCWb is inhibited Klughammer et al. (2018a). These findings suggest that SCWs may function primarily as a coordinated cortical response to an underlying cytoplasmic signal rather than as an autonomous cortical patterning mechanism Bement et al. (2024); Liu et al. (2025). Here, we build on and extend this interpretation by providing a mechanistic framework that directly links cytoplasmic wave geometry to cortical reactivation patterns, showing quantitatively how nuclear size, cell size, and inhibition kinetics together determine the directionality and mode of cortical response.

In this work, we test whether differences in SCW directionality can be explained by geometry-dependent propagation of cytoplasmic Cdk1 waves rather than by cortex-specific patterning mechanisms. We then examine how these cytoplasmic waveforms are transmitted to the cortex by coupling Cdk1 activity to an excitable Rho–actin network through Ect2 inhibition. In this framework, surface contraction waves emerge as downstream consequences of cytoplasmic signaling geometry and recovery dynamics.

Here we begin by analyzing spatiotemporal Cdk1 dynamics and asking how geometric parameters shape cytoplasmic Cdk1 waves (Fig. 1A,B). Using a spatially extended model of the Cdk1 oscillator Yang and Ferrell Jr (2013); Puls et al. (2024), we show that Cdk1 waves generically consist of two mechanistically distinct components. Activation fronts propagate as traveling waves with approximately constant velocity, consistent with trigger-wave behavior, whereas wave backs are controlled by Cdk1 gradients shaped by diffusion rather than by local bistable kinetics (Fig. 1B,C). As a consequence, wave fronts and wave backs can propagate at different speeds and, depending on nuclear Cyclin B–Cdk1 content and effective system size, even in opposite directions. These transitions between traveling-wave-dominated and gradient-dominated regimes provide a simple framework for interpreting differences in surface contraction wave directionality observed across biological systems (Fig. 1D).

Results

Geometry-controlled directionality of cytoplasmic Cdk1 waves

During both mitosis and meiosis, Cyclin B–Cdk1 accumulates in the nucleus prior to NEB. Upon NEB, active Cyclin B–Cdk1 is released into the cytoplasm, where it propagates spatially and temporally according to the spatiotemporal Cdk1 model described in Methods and Materials: The cytoplasmic Cyclin B–Cdk1 model. This cytoplasmic signal provides the upstream input for the cortical dynamics analyzed in the following section. Here, we first focus on the properties of the cytoplasmic Cdk1 waves themselves.

Wave fronts, wave backs, and directionality

Following NEB, the spatial Cdk1 profile exhibits a characteristic wave-like structure consisting of two distinct components: an activation front, marking the transition from low to high Cdk1 activity, and a wave back, corresponding to the decay of Cdk1 activity after the peak. Throughout this section, we quantify propagation by tracking the positions of the wave front and wave back in kymographs extracted along a one-dimensional radial coordinate extending from the nuclear region toward the cell periphery (see Fig. 1A, C, D).

From these trajectories, we define the front velocity v_F and back velocity v_B as the slopes of the corresponding features in space–time plots (see Supplementary Fig. S2). Positive velocities correspond to propagation away from the nucleus. While v_F is always positive, v_B may be either positive or negative, providing a quantitative measure of wave reversibility and defining the sign of back propagation used throughout this section.

To investigate how transitions between positive and negative v_B arise, we considered two control parameters. First, we designed initial conditions that allow controlled manipulation of the spatial Cdk1 distribution. Specifically, we initialized the system with a high total Cyclin B–Cdk1 concentration throughout the cell, that can be increased in the nucleus, while restricting high levels of active Cyclin B–Cdk1 to the nucleus. This was achieved by selecting two points along the limit cycle and introducing a nuclear scaling factor γ, which multiplicatively enhances the nuclear concentration of total and active Cyclin B–Cdk1 relative to the cytoplasmic baseline (see Supplementary Fig. S1). Second, we varied the diffusion coefficient, which is mathematically equivalent to changing the effective system size.

Baseline regime: traveling-wave–dominated propagation

We first establish a reference regime in which the cytoplasmic Cdk1 wave exhibits traveling-wave– like behavior. In systems with moderate diffusion coefficient and low nuclear scaling factor, the activation front propagates with approximately constant velocity, producing a linear trajectory in kymographs (constant-speed propagation) (see Fig. 2A (i)). In this regime, the wave back follows the same direction as the front, with v_B > 0, and the overall wave appears as a single traveling entity.

Figure 2.

Front–back asymmetry: traveling fronts and gradient-controlled backs

Across a broad range of parameters, we observe a qualitative asymmetry between the wave front and the wave back. While the activation front often propagates with nearly constant velocity, consistent with trigger-wave–like behavior, the wave back does not generically exhibit constant-speed propagation. Instead, back trajectories frequently display curvature or piecewise-linear behavior, indicating gradient-controlled dynamics rather than local bistable kinetics (see Fig. 2A (ii), (iii) and (iv) and Supplementary Fig. S4A). As a result, the front and back velocities respond differently to changes in nuclear Cdk1 content and system size. This asymmetry allows the wave back to slow down, stall, or reverse direction even when the activation front continues to propagate outward. This front–back asymmetry is the key mechanism underlying direction reversals in the regimes described below.

Although our results do not show substantial changes in the wave front, theory predicts that in the small-system limit the front behavior also changes, becoming purely diffusion-dominated. In this regime, the distance beyond which a traveling wave overtakes diffusive spreading can be estimated for an initially localized region of high concentration diffusing in space as

where [CyB-Cdk1]_N is the initial concentration in the nucleus, [CyB-Cdk1]_C is the intial concentration in the rest of the cell, and [CyB-Cdk1]_Th is the threshold concentration required to excite the system (see Supplemental Material). This diffusion-dominated limit appears as an initial transient even in large systems; however, when the system size is much larger than L_Crit, its effect on the measured front velocity is minimal. Consistent with this prediction, we observe an initial diffusion-dominated phase in Supplementary Fig. S2 before a constant-speed front is established. Accordingly, throughout the remainder of this paper we focus on regimes in which a well-defined traveling front exists.

Nuclear Cyclin B–Cdk1 content controls wave-back directionality

We tested the effect of nuclear accumulation of Cyclin B–Cdk1 by varying the control parameter γ. Keeping the diffusion coefficient fixed at 1.5μm²/s, we found that increasing the nuclear concentration caused the wave back to transition from positive velocity, through a regime of coexisting wave behaviors, to negative velocity (see Fig. 2D). In addition, increasing the amount of nuclear Cyclin B–Cdk1 either by raising its concentration or by enlarging the nuclear size produced equivalent effects (see Supplementary Fig. S3A). Therefore, these variables can be collapsed into a single control variable that is the total nuclear content of active Cyclin B–Cdk1 (see Fig. 2D and Supplementary Fig. S3B). During the transient regime of coexisting velocities, the positive wave-back velocity appears to diverge to infinity, while the negative velocity converged from infinity toward a value close to the wave-front velocity (see Fig. 2D). This apparent divergence arises because the spatial domain over which the corresponding velocity is defined shrinks to zero (see Supplementary Fig. S3C). Intuitively, when the total amount of Cyclin B–Cdk1 is low, the system supports a traveling wave. As this amount increases, two processes coexist: a propagating wave front and a slower diffusive spread. Because diffusion is slower, regions farther from the source —here the nucleus— begin to relax earlier, having not yet been reached by the diffusive component. This establishes a gradient that drives relaxation from the most distant regions back toward the nucleus, resulting in a backward-propagating wave back (see Fig. 2B).

System size and diffusion modulate gradient effectiveness

Because the proposed mechanism for wave-back reversal relies on diffusion-driven gradients, we hypothesized that it would be sensitive to diffusion timescale. To test this, we used the diffusion coefficient as our second control parameter, which is mathematically equivalent to changing the system size for fixed reaction rates. In small systems (large diffusion coefficients), large nuclear concentrations of active and total Cyclin B–Cdk1 diffuse rapidly throughout the domain, producing a nearly homogeneous decay of Cdk1 activity and eliminating directional wave-back propagation (see Fig. 2A (iv)). In larger systems (small diffusion coefficients), gradients persist over longer distances, allowing the wave back to exhibit hybrid behavior consisting of both traveling-wave– dominated and diffusion-driven components (see Fig. 2A (iii)). The range of total nuclear Cdk1 content over which such hybrid behavior occurs narrows as diffusion increases, consistent with the asymptotic behavior of the fits in Fig. 2D and Supplementary Fig. S4E.

Phase diagram of Cdk1 waves

Together, these results define a phase diagram of the cytoplasmic Cdk1 wave behavior as a function of the total nuclear active Cyclin B–Cdk1 concentration and the diffusion timescale (or system size). Distinct regimes can be identified based on the sign and spatial structure of the wave-back velocity: (i) a traveling-wave–dominated regime with v_B > 0, (ii) a counterpropagating-back regime with v_B < 0, (iii) a hybrid regime with coexistence of both behaviors, and (iv) a fast-diffusion regime in which no well-defined wave propagation occurs.

These regimes provide a simple geometric framework for interpreting organism- and stage-specific differences in observed wave directionality. Returning to the two prototypical cell types introduced in the top and bottom rows of Fig. 1B,C, we find that our results are consistent with the distinct characteristics of these systems. An embryo of Xenopus laevis during its first mitosis (approximately 1 mm in diameter) can be considered a very large system in comparison with the size of its nucleous. At the same time, a starfish oocyte in meiosis represents a smaller system (approximately 0.15 mm in diameter) with a much larger germinal vesicle. We now turn to the cortical network and examine how Cdk1 waveforms are transmitted to and converted into Rho–actin dynamics in the excitable cortex.

Cytoplasmic Cdk1 waves shape cortical dynamics

The cortex as a Cdk1-gated excitable system

Rather than treating SCWs as autonomous cortical phenomena, we view them as downstream dynamical consequences of cytoplasmic Cdk1 activity. The cortex functions as an excitable system whose dynamics is transiently gated by Cdk1-mediated inhibition. Cytoplasmic Cdk1 activity is transmitted to the cortex primarily through phosphorylation of the RhoGEF Ect2, which reduces its membrane affinity and suppresses Rho GTPase signaling Hara et al. (2006); Niiya et al. (2006); Su et al. (2011). As a consequence, changes in Cdk1 activity modulate the distance of the cortical network to a Hopf bifurcation separating quiescent and oscillatory regimes. SCWs thus emerge as phase waves (i.e., apparent wave propagation arising from spatially staggered threshold crossing, rather than local signal propagation) generated when different cortical regions cross this bifurcation at different times.

Box 2. Rho-actin dynamics and cortical waves

Box 2—figure 1.

Cortical excitability and control by Ect2

The cortical dynamics is described using an activator–depleted substrate–inhibitor (ADSI) model of the Rho GTPase network Bement et al. (2015); Michaud et al. (2022); Chomchai et al. (2024). The model incorporates Rho-GTP, Rho-GDP, and F-actin, alongside other regulators such as Ect2 and the GAP RGA-3/4—which are accounted for via parameters—and it captures transient oscillations as well as a wide range of spatial patterns. The effective Ect2 concentration α is a key control parameter that determines the dynamic regime of the model (see Box 2C). For low α, the cortex remains quiescent. Increasing α induces oscillations and spatial patterning (see Box 2). In the absence of inhibition by the Cyclin B-Cdk1 complex, the system expressing sufficiently high levels of Ect2 (0.59 < α < 1.26) does not spontaneously reset on the mitotic time scale (i.e., it remains in the active regime characterised by propagating waves).

Coupling cytoplasmic Cdk1 and the cortex dynamics

The interaction between cytoplasmic Cdk1 and the cortex is implemented by inhibiting Ect2 activity, using a reduction in the parameter α, which accounts for effective Ect2 concentration. Two types of inhibitory signals are considered. First, an idealized Heaviside-like inhibition, which approximates an abrupt on–off Cdk1 signal. Second, an inhibition derived directly from the cytoplasmic Cdk1 concentration profiles, which preserves the slow recovery and rapid drop of the Cdk1 wave-form following NEB. Both signals suppress cortical activity, but they differ in the way inhibition is released. The Heaviside signal crosses the Hopf bifurcation abruptly, whereas the Cdk1-derived signal exhibits a slow recovery followed by a fast drop. This difference can have a strong effect on how cortical activity reappears after inhibition, as we show later in this paper.

Besides the inhibitory function, we introduce two different scenarios depending on a choice of parameters. The functional term describing F-actin disassembly can be either homogeneous or heterogeneous, where heterogeneity accounts for additional factors not explicitly included in the model that also influence F-actin disassembly. As we show later, the way disassembly is implemented, together with the functional form of the inhibition, plays an important role in the resulting dynamics. Together, this defines a 2×2 comparison: (i) abrupt vs waveform-derived inhibition, each under (a) homogeneous vs (b) heterogeneous actin turnover.

Baseline cortical response to idealized inhibition

Under Heaviside inhibition, cortical activity is rapidly suppressed and the system approaches a homogeneous steady state. After inhibition is lifted, activity reappears as waves through two distinct mechanisms (Supplementary Videos 6,7; Fig. 3D,E). In heterogeneous systems, parameter inhomogeneities allow for an immediate recovery of complex dynamics, requiring only a single reactivation wave before spiral dynamics re-emerge. In contrast, in homogeneous systems, the recovery is delayed and proceeds through several wave cycles before spiral turbulence re-emerges. This behavior can be observed in Supplementary Videos 6 and 7 and is also reflected in the oscillation period of the system, since spiral dynamics in excitable systems is characterized by a reduction in the effective period Cebrián-Lacasa et al. (2024). This differences arise because, as the system approaches the stable state during inhibition, phase differences between the neighboring points — induced earlier by the pre-inhibition pattern — are reduced by diffusion. Upon reactivation, a homogeneous system therefore requires time to amplify residual concentration inhomogeneities, whereas in a heterogeneous system the imposed parameter variations immediately seed pattern formation, eliminating this delay.

Figure 3.

Cortical response to realistic Cdk1 waveforms in two dimensions

Under Cdk1 inhibition, cortical activity is also rapidly suppressed and the system approaches a homogeneous steady state. In a homogeneous system, the dynamics is equivalent to that observed under the Heaviside inhibition: several planar waves propagate through the system, leading to a slowing of the oscillations compared to the spiral regime (Supplementary Video 8; Fig. 3F). The same rationale as was used in the Heaviside case therefore applies here.

In contrast, heterogeneous systems exhibit a qualitatively different behavior (Supplementary Video 5). As shown in Fig. 3A, reactivation occurs through a propagating front of target patterns, or so-called bubble-like patterns (target-like oscillatory domains seeded by local pacemakers), rather than through planar waves. As a consequence, the oscillation period is set by the period of these target patterns, which is comparable to the period of wave cycles and to the intrinsic period of the system, rather than by the shorter period associated with spiral dynamics (see Fig. 3B,C). After a few oscillation periods, the bubble-like pattern evolves into a multi-spiral regime, recovering dynamics similar to that observed prior to inhibition and to the one of the corresponding Heaviside inhibition case. As we show below, this differences arise from the combined effects of slow recovery, parameter heterogeneities, and the manner in which the system crosses the Hopf bifurcation.

Cortical dynamics in a spherical geometry

We next extend our system to a three-dimensional cortical shell and use the different regimes described and shown in Fig. 2D as input cytoplasmic signals regulating Ect2 activity (see Fig. 4A; Supplementary Videos 9–12). The resulting cortical Rho GTPase signal exhibits the four behaviors previously described: a wave back propagating parallel to the wave front, a counter-propagating wave back, coexistence of these two regimes, and a flat recovery occurring simultaneously across the entire cell. The first and second behaviors correspond to the scenarios reported in experimental observations (Fig. 1D), whereas the latter may provide a good description of surface contraction waves once the embryo has divided multiple times and cell size is strongly reduced. Nevertheless, the mode of reactivation in these simulations does not fully match the experimental results reported by Bement et al. (2015), where cortical activity reappears after Cdk1 inhibition as an almost planar front (see Fig. 4D; Supplementary Video 18). This discrepancy suggests that not only the spatial ordering of inhibition, but also the temporal profile and strength of inhibition, control the reactivation mode — an issue we address next.

Figure 4.

The role of wave profile and inhibition time

We therefore returned to the results shown in Fig. 3 to understand why our simulated dynamics do not directly reproduce the experimental observations. In that figure, we showed that Heaviside inhibition is capable of producing planar waves upon reactivation. However, when F-actin disassembly is homogeneous and inhibition times are sufficiently long, several wave cycles can emerge even in the three-dimensional geometry, which again does not resemble the experimental behavior (see Supplemental Fig. S6A,B; Supplementary Video 13,14). In contrast, introducing heterogeneous F-actin disassembly is sufficient to reproduce a single reactivation wave followed by immediate spiral dynamics (see Supplemental Fig. S6C; Supplementary Video 15). Thus, long inhibition times in homogeneous systems favor multi-cycle recovery, whereas heterogeneity promotes rapid, single-wave reactivation.

Inhibition strength selects the reactivation mode

These observations lead us to conclude that the perfectly flat Heaviside inhibition profile can qualitatively bias reactivation compared to waveform-derived inhibition, which is characterized by a slow recovery followed by a fast decay. To test this, we varied the phosphorylation rate, or inhibitory strength (s), with which Cdk1 interacts with Ect2, and found that the strength of Ect2 inhibition determines how cortical activity reappears (Fig. 4B,C). Strong inhibition drives the system across the Hopf bifurcation during the fast recovery phase of Ect2, producing immediate planar phase waves followed by spiral turbulence, in agreement with experimental observations (Fig. 4C, right; Supplementary Video 17). In contrast, weaker inhibition crosses the bifurcation during the slow recovery phase, allowing spatial heterogeneities to trigger oscillations locally while other regions remain excitable. This results in bubble-like target patterns and transient oscillations before full pattern re-establishment (Fig. 4C, left; Supplementary Video 16). This type of dynamics is absent in idealized inhibition models.

Different reactivation modes have different activity properties

We next asked whether the pattern that characterizes reactivation can influence the mechanical feedback. To quantify this, we measured the activity proportion of the system, meaning the fraction of the area that exhibits high Rho GTP relative to the total area considered, as well as the maximum Rho GTP intensity (see Fig. 4C).

We performed this quantification in a section around the equator of the cell. Although the bubble-like pattern appears earlier during activation, both its intensity and the area it covers are much smaller than those observed for the spiral pattern. Only when the spiral pattern is reactivated does a larger fraction of the system display high Rho GTP levels, together with a higher maximum activity.

Alternative ways of controlling Rho GTPase activation

We conclude the Results section by showing that inhibitory profiles and cytoplasmic activity levels are not the only factors controlling cortical dynamics. As shown by Chang and Ferrell Jr (2013), since activity propagates as a concentric wave from the nucleus, varying the nuclear position alone is sufficient to alter not only the wave back but also the wave front. This is illustrated in Supplemental Fig. S7, where changing the position of the nucleus produces a range of behaviors, from a clearly defined phase wave when the nucleus is close to the cortex, to a completely flat activation profile in which Cdk1 reaches the cortex simultaneously when the nucleus is centrally located (corresponding to an effectively infinite apparent phase velocity).

Discussion

Surface contraction waves reflect how large cells coordinate cortical mechanics with cell cycle progression. A puzzling observation has been why these waves propagate in opposite directions in starfish oocytes Bischof et al. (2017) and Xenopus embryos Chang and Ferrell Jr (2013), two systems that share the same core biochemical machinery. Our results show that this difference does not require organism-specific cortical mechanisms. Instead, it follows directly from how nuclear geometry and effective cell size shape the cytoplasmic Cdk1 waveform through an asymmetry in the wave front and back.

The asymmetry arises because the activation front and the wave back are governed by fundamentally different physical mechanisms. Activation fronts propagate as trigger waves driven by local bistable dynamics Gelens et al. (2014); Beta and Kruse (2017); Deneke and Di Talia (2018), with approximately constant velocity set by the kinetics of the Cdk1 oscillator. Wave backs, by contrast, are not governed by local kinetics but by diffusion-driven Cdk1 gradients that develop after nuclear envelope breakdown. Because these two components respond differently to nuclear Cdk1 content and system size, they can propagate at different speeds or even in opposite directions. In large systems with a relatively small nucleus, such as a first-mitosis Xenopus embryo, the traveling-wave component dominates and front and back co-propagate outward. In smaller systems with a proportionally larger germinal vesicle, such as a starfish oocyte in meiosis II, the diffusion-driven gradient dominates the wave back and drives it inward, counter to the activation front. This single geometric framework, summarized in the phase diagram of Fig. 2C, accounts for both classes of behavior within a conserved network without invoking any cortex-specific asymmetries, and is consistent with earlier gradient-based interpretations of SCW directionality while extending them to explain both front and back propagation simultaneously.

The transmission of this cytoplasmic signal to the cortex adds a second layer of geometry-dependence. Cdk1-mediated phosphorylation of Ect2 provides the biochemical link between the cell-cycle oscillator and the Rho–actin cortex Hara et al. (2006); Niiya et al. (2006); Su et al. (2011), reducing Ect2 membrane affinity and transiently driving the cortical network below the Hopf bifurcation that separates quiescent and oscillatory regimes. As a consequence, cortical activity does not propagate by local activation along the cortex. Instead, different cortical regions independently cross the bifurcation threshold as inhibition is released, at times determined by the spatiotemporal profile of the cytoplasmic Cdk1 waveform. Surface contraction waves therefore correspond to phase waves, apparent propagation arising from spatially staggered threshold crossing, whose directionality is inherited from the underlying cytoplasmic signal.

The temporal structure of the inhibitory waveform, not only its spatial ordering, determines the mode of cortical reactivation. Rapid release of Cdk1-mediated inhibition synchronizes cortical recovery across the cell, producing a coherent planar reactivation wave consistent with experimental observations in starfish oocytes following Cdk1 inhibition Bement et al. (2015); Liu et al. (2025). Slower release, by contrast, exposes spatial heterogeneities in actin turnover. Regions with locally shorter oscillation periods cross the Hopf bifurcation earlier than their surroundings and act as pacemakers, generating transient bubble-like target patterns before spiral dynamics re-emerge. This class of dynamics is entirely absent from idealized on–off inhibition models, and arises specifically because the relaxation kinetics of the Cdk1 oscillator produce a slow recovery phase that amplifies pre-existing spatial variability in cortical excitability. The strength of Ect2 inhibition selects between these two reactivation modes: strong inhibition drives the cortex across the bifurcation during the fast recovery phase, synchronizing reactivation, while weaker inhibition allows heterogeneities to trigger local oscillations during the slow phase. Together, inhibition strength and recovery kinetics constitute a second geometric axis, one defined by the temporal rather than the spatial waveform, that shapes the pattern of cortical reactivation.

Our framework does not incorporate mechanical feedback. Cortical tension, cytoplasmic flows, and surface deformation have been proposed to influence SCW propagation Miller et al. (2018); Klughammer et al. (2018a); Yin et al. (2022), and more broadly mechano-chemical coupling has been identified as an important contributor to actin wave dynamics in many cell types Beta et al. (2023). Such effects may further modulate the dynamics we describe. However, our results indicate that much of the observed diversity in SCW behavior, including directionality, reactivation mode, and the transition between them, can be explained at the biochemical level, through the geometry of cytoplasmic signaling and the kinetics of cortical inhibition release. Mechanical coupling may refine these dynamics, for example by shifting the effective Ect2 control parameter or the degree of cortical heterogeneity, but it is not required to generate the primary phenomenology. Testing this will require experiments that decouple mechanical and biochemical contributions, for instance by comparing SCW directionality in cells where cortical tension is pharmacologically clamped against predictions from the phase diagram in Fig. 2C.

Taken together, our results establish that the contrasting wave behaviors observed across organisms reflect differences in nuclear-to-cell size ratio and nuclear Cdk1 content, not organism-specific cortical patterning mechanisms. Cell geometry, acting through a front–back asymmetry in the Cdk1 oscillator, is sufficient to account for this diversity within a single conserved biochemical network. The framework makes concrete, testable predictions: steeper or stronger Cdk1-dependent Ect2 inhibition should favor single coherent planar reactivation waves, while prolonged inhibition or increased spatial heterogeneity in actin turnover should favor pacemaker-driven recovery. More broadly, the front–back decoupling described here is not specific to the Cdk1 oscillator. Any excitable or oscillatory system in which a localized source drives waves across a large domain may exhibit similar geometry-dependent transitions Gelens et al. (2014); Beta and Kruse (2017); Deneke and Di Talia (2018), making this a potentially general principle for spatial coordination of the cell cycle and beyond.

Methods and Materials

The cytoplasmic Cyclin B–Cdk1 model

We model cytoplasmic Cyclin B–Cdk1 dynamics using the spatially extended oscillator Yang and Ferrell Jr (2013); Puls et al. (2024). The model tracks total Cyclin B–Cdk1 concentration, c(x, t), and active Cyclin B–Cdk1 concentration, a(x, t), and incorporates feedback through Cdc25, Wee1, and APC/C. The governing equations are

where D is the diffusion coefficient, assumed equal for active and total Cyclin B–Cdk1 for simplicity; allowing distinct diffusion coefficients does not qualitatively alter the results (data not shown).

The regulatory terms are given by ultrasensitive Hill functions,

All parameters are taken from Yang and Ferrell Jr (2013) and lie within experimentally measured ranges Tsai et al. (2014) (see Methods and Materials: Model parameters). The model exhibits relaxation oscillations characterized by rapid activation of Cyclin B–Cdk1 followed by slower decay driven by APC/C-mediated Cyclin B degradation.

To model cytoplasmic wave propagation following nuclear envelope breakdown (NEB), Eq. (2) was solved on a one-dimensional radial domain while retaining a three-dimensional diffusive Laplacian under the assumption of spherical symmetry. The origin of the coordinate system was placed at the center of the nucleus and no-flux (Neumann) boundary conditions were imposed at both the nuclear boundary and the cell cortex.

Initial conditions were chosen such that the nucleus contained elevated concentrations of both total and active Cyclin B–Cdk1, while the cytoplasm was initialized with low active and high total Cyclin B–Cdk1. Increased nuclear concentrations were modeled by multiplying the initial concentrations by a constant factor γ, applied only within the nuclear region.

Wave fronts were defined as the first threshold crossing of active Cyclin B–Cdk1 during activation, while wave backs were defined as the return crossing of the same threshold during decay (the threshold was set to the midpoint between the low and high activity steady states unless otherwise stated). Wave velocities were extracted from space–time kymographs using linear regression over the region where trajectories were approximately linear; in regimes with curved trajectories, piecewise fits were used as described in Supplementary Fig. S2.

System-size effects were explored by varying the diffusion coefficient D. Increasing D increases the characteristic diffusion length over a fixed timescale and is mathematically equivalent to decreasing system size for fixed reaction kinetics and domain geometry (i.e., after nondimensionalization, only the ratio D/T controls the effective length scale).

The cortical ADSI model and bifurcation structure

Cortical dynamics were described using an activator–depleted substrate–inhibitor (ADSI) model adapted from Michaud et al. Michaud et al. (2022). The model tracks active Rho GTPase (u), inactive Rho GTPase (v), and F-actin (w):

The reaction term is

The parameter α represents the effective Ect2 concentration and controls the distance to a subcritical Hopf bifurcation of the spatially homogeneous system. Once this bifurcation is crossed, oscillation amplitudes rapidly increase, leading to relaxation-like cortical dynamics that can support excitable wave propagation and spiral turbulence in space. All parameters are retained from Michaud et al. (2022) and lie within experimentally supported ranges (see Table 1). Spatial heterogeneities in actin turnover were introduced by adding Gaussian-correlated noise to the F-actin degradation rate k₁₀. The noise field was characterized by a fixed correlation length and amplitude (see Methods: Numerical details) and was kept constant in time (quenched disorder). These heterogeneities locally shift the Hopf bifurcation threshold, allowing some regions to enter oscillatory regimes earlier than others. This leads to the emergence of localized pacemaker regions and target-like patterns during cortical reactivation when the system is slowly driven across the bifurcation by Cdk1-dependent inhibition.

Table 1.

Coupling cytoplasmic Cdk1 to cortical inhibition

The interaction between cytoplasmic Cyclin B–Cdk1 and the cortex was modeled as an inhibitory effect on the effective Ect2 concentration α, consistent with experimental evidence that Cdk1 phosphorylation reduces Ect2 membrane association Hara et al. (2006); Niiya et al. (2006); Su et al. (2011). In this framework, cytoplasmic Cdk1 activity transiently shifts the cortical system across the Hopf bifurcation that separates quiescent and oscillatory regimes.

Inhibition derived directly from the cytoplasmic signal was implemented as

where s is the inhibition strength, a₀ is the basal active Cyclin B–Cdk1 concentration, and a(x, t) is the spatiotemporal cytoplasmic Cdk1 field obtained from Eq. (2). In practice, α was bounded below by zero to avoid unphysical negative effective Ect2 concentrations.

An idealized Heaviside inhibition was also considered,

where t₁ − t_NEB defines the inhibition duration. This idealization produces an abrupt crossing of the Hopf bifurcation, in contrast to the waveform-derived inhibition, which yields a slow recovery followed by a rapid decay in α(t).

In homogeneous cortical systems, longer inhibition drives the system closer to the stable fixed point, reducing residual spatial inhomogeneities. Upon release, this results in multiple planar traveling waves before spiral turbulence re-emerges. Shorter inhibition times or steeper release profiles synchronize reactivation and favor single planar phase waves.

In systems with heterogeneous actin degradation, localized pacemaker regions persist during inhibition, enabling immediate reactivation regardless of inhibition duration. When inhibition is released slowly, these regions cross the Hopf bifurcation earlier than their surroundings and seed target-like patterns. These behaviors are consistent with recovery dynamics in excitable media Cebrián-Lacasa et al. (2024).

Cortical dynamics were also simulated on a spherical surface representing an early embryonic cell. The cortical domain was implemented as a closed spherical shell with no-flux boundary conditions, and cytoplasmic inhibition was applied through the spatiotemporal Cdk1 field projected onto the cortex. Cytoplasmic inhibition was initiated from a localized nuclear pacemaker positioned off-center along the animal–vegetal axis.

Shifting the pacemaker position altered the relative timing of inhibition across the cortex. In the limit of a centrally located nucleus, all cortical regions are inhibited simultaneously, corresponding to a diverging apparent phase velocity (i.e., a spatially uniform phase reset). These effects reproduce experimentally observed transitions between U-shaped and V-shaped SCWs Chang and Ferrell Jr (2013).

Model parameters

In this study, we used two models and retained the parameter values from their original publications Yang and Ferrell Jr (2013); Michaud et al. (2022). These values are listed in Table 1 and were used throughout unless explicitly stated otherwise. In addition to the parameters listed in Table 1, the spatial parameters were specified as follows. For the cytoplasmic model, the system size was set to d = 125, μm, with a nuclear radius of 30, μm. For the cortical model, the system radius was R = 75, μm, and the nuclear position was set to 40, μm from the center. Heterogeneities in the cortical model were generated using spatially correlated Gaussian noise applied to the F-actin degradation rate, with a correlation length of 4, μm and an amplitude of 2 (variance set by the noise-generation procedure). The control parameter α was set to a constant value of 1 prior to nuclear envelope breakdown (NEB), which occurred at time t_NEB = 5, min.

Numerical methods

Both models were simulated using a second-order Runge–Kutta method for time integration, implemented manually. For spatial discretization, the cytoplasmic model and the two-dimensional cortical simulations used a finite-difference scheme, while the three-dimensional cortical simulations employed a finite-element method implemented using FEniCSx.

Unless otherwise stated, spatial and temporal discretization parameters (grid spacing, time step, and domain resolution) were kept fixed across simulations and are provided in the code referenced in the Code Availability section.

To account for dimensional effects in the cytoplasmic model, we used the Laplacian operator corresponding to three dimensions under the assumption of spherical symmetry:

where f (r) is a radially symmetric function. Under this assumption, the dynamics depend only on the radial coordinate, while preserving the effective three-dimensional diffusive geometry. The resulting solution can be projected onto any spherical surface representing the cell cortex. No-flux (Neumann) boundary conditions were imposed at all domain boundaries.

Supplemental material

The critical length at which wave propagation becomes more efficient than diffusion can be derived as follows. Consider a point source with a high initial concentration that evolves according to a Gaussian distribution,

Assuming that the concentration must reach a threshold value C_Th in order to excite the system, we can solve for the position x(t) at which C(x, t) = C_Th and obtain

The velocity of this diffusive front is obtained by taking the time derivative,

Since the wave velocity v_W is constant, while the velocity of the diffusive front decreases with time, we determine the crossover point by imposing v_W = v_D. This condition defines the critical time t_Crit,

Finally, using the distance traveled by the wave, L = v_Wt, and substituting t_Crit, we obtain the critical length,

which can be further corrected in the presence of a basal concentration level to recover the expression given in the main text.

Supplemental figures and videos

The initial conditions and the nuclear increase in total and active Cyclin B–Cdk1 are shown in Supplementary Figure S1. Total Cyclin B–Cdk1 displays a high concentration throughout the entire domain, both inside and outside the nucleus, whereas the active form is predominantly concentrated within the nucleus. The specific initial values were extracted from trajectories along the underlying limit cycle. The applied nuclear scaling factor rescales both total and active Cyclin B–Cdk1 proportionally, thereby preserving their relative ratio.

Supplementary Figure S2 shows the fits and residuals of the wave fronts and wave backs for different values nuclear scaling factor, corresponding to a subset of the velocities shown in Fig. 2C. The fits were performed using NumPy’s polyfit function with a first-order polynomial, applied either globally or piecewise (one or two segments) depending on the behavior.

Figure S1.

As expected, the wave fronts exhibit an approximately constant velocity, reflected in good linear fits. However, as the nuclear scaling factor increases, the residuals become bimodal. This reflects an initial phase dominated by diffusion rather than a traveling front. This also explains the slight variation in fitted velocity and suggests that the front dynamics can be decomposed into two phases: an initial diffusion-driven segment followed by a traveling-front-driven segment.

Figure S2.

In contrast, the wave backs are poorly approximated by a single linear fit. This is because the established Cdk1 gradient plays a dominant role in the dynamics. As nuclear scaling factor increases, the linear approximation becomes less accurate due to the increasingly diffusive character of the gradient, and two segments with opposite apparent velocities can be identified.

We next examined how this hybrid wave back, consisting of traveling-wave–dominated and diffusion-driven segments, varies (Supplementary Figure S3). First, we verified that the transition in wave-back behavior is controlled by the total amount of Cyclin B–Cdk1 by modifying the nuclear size instead of the nuclear overexpression. This produces equivalent behavior after renormalization, showing that both manipulations collapse onto a single control variable. We then quantified how the transition occurs by measuring the fraction of the domain governed by traveling-wave–dominated versus diffusion-driven dynamics.

To further complement Fig. 2, we introduce Supplementary Fig. S4A–C to illustrate the effect of system size—or, equivalently, the diffusion coefficient—on wave dynamics. A large diffusion coefficient (corresponding to a small effective system size) leads to a nearly homogeneous distribution of active Cyclin B–Cdk1, and as a result, no well-defined wave propagation is observed. In contrast, a small diffusion coefficient (i.e., a larger system size) produces a spatial gradient confined to part of the domain, giving rise to a wave back in which both counterpropagating and forward-propagating segments are present. As the system size increases further, the forward-propagating component becomes more dominant in shaping the overall system behavior.

Figure S3.

To characterize these transitions, we quantified the wave velocities for different nuclear scaling factors (see Supplementary Fig. S4D). Furthermore, we investigated how the region of coexistence of oppositely propagating wave-back segments narrows as the diffusion coefficient increases. This region is defined by the positions of the two asymptotes of the hyperbolic fit, which converge as the diffusion coefficient increases (see Supplementary Fig. S4E).

In Supplementary Fig. S5, we show how operating close to the Hopf bifurcation, in the presence of heterogeneous F-actin degradation, leads to target-pattern dynamics that we refer to as bubble-like patterns. This behavior arises from the excitable nature of the cortical system combined with spatial differences in local oscillation periods (Supplementary Fig. S5B). Regions with shorter intrinsic periods act as pacemakers that excite neighboring regions, generating outward-propagating waves from these localized sources (Supplementary Fig. S5A).

Supplementary Fig. S6 examines the effect of inhibition duration under idealized homogeneous actin degradation. In this case, spatial patterns emerge only from inhomogeneities in the dynamical variables, since the degradation term is uniform. During inhibition, the system relaxes toward its steady state. For short inhibition periods, residual spatial inhomogeneities remain large, and pattern reactivation occurs rapidly (Supplementary Fig. S6A). For longer inhibition periods, these inhomogeneities are further reduced by diffusion, and several planar waves appear before spiral or turbulent dynamics re-emerge (Supplementary Fig. S6B).

When F-actin degradation is heterogeneous, parameter inhomogeneities alone are sufficient to restart the pattern, leading to immediate reactivation regardless of inhibition duration (Supplementary Fig. S6C). In this case, the heterogeneities act as pacemakers that initiate and coordinate pattern reactivation, often before the inhibition phase has fully ended.

Figure S4.

Figure S5.

The position of the nucleus directly affects how the wave appears at the embryo surface, as shown in Supplementary Fig. S7. In panel C, we compare the front and back of the Ect2 wave and show that activation is more synchronous when the nucleus is closer to the center (z = 10 μm). The frames in panel D further illustrate that, in this configuration, cortical activity disappears and re-emerges more uniformly than in the off-center case shown in the main text.

Figure S6.

Figure S7.

Data availability

The numerical codes to reproduce the figures in this study are openly available in a GitLab Repository: "Geometry shapes cytoplasmic Cdk1 waves that drive cortical dynamics". GitLab; 2026.691 (https://gitlab.kuleuven.be/gelenslab/publications/coupledsystems_biomodels.692), and as an archived repository in RDR by KU Leuven [Upcoming]. The generated data is openly available on RDR [Upcoming].

Additional files

Video 1. Projection of Cyclin B–Cdk1 activity onto the surface of a sphere in the traveling-wave–dominated regime (i) (Fig. 2D).

Video 2. Projection of Cyclin B–Cdk1 activity onto the surface of a sphere in the front–back asymmetry regime (ii) (Fig. 2D).

Video 3. Transition from a leopard-like pattern to a spiral pattern for α = 0.8 (Fig. 3C).

Video 4. Bubble-like regime in an overall non-oscillatory system, with parameter heterogeneities acting as pacemakers for α = 0.55 (Supplementary Fig. S5).

Video 5. Rapid spiral dynamics induced by heterogeneities for α = 0.8.

Video 6. Spatiotemporal dynamics in a two-dimensional simulation with Heaviside inhibition and heterogeneous F-actin disassembly (Fig. 3D).

Video 7. Spatiotemporal dynamics in a two-dimensional simulation with Heaviside inhibition and homogeneous F-actin disassembly (Fig. 3E).

Video 8. Spatiotemporal dynamics in a two-dimensional simulation with Cdk1-derived inhibition and heterogeneous F-actin disassembly (Fig. 3F).

Video 9. Spatiotemporal dynamics on a spherical cortical shell with underlying regime (i) cytoplasmic dynamics (Fig. 4A).

Video 10. Spatiotemporal dynamics on a spherical cortical shell with underlying regime (ii) cytoplasmic dynamics (Fig. 4A).

Video 11. Spatiotemporal dynamics on a spherical cortical shell with underlying regime (iii) cytoplasmic dynamics (Fig. 4A).

Video 12. Spatiotemporal dynamics on a spherical cortical shell with underlying regime (iv) cytoplasmic dynamics (Fig. 4A).

Video 13. Effect of inhibition duration with homogeneous F-actin disassembly: short inhibition time (Supplementary Fig. S6A).

Video 14. Effect of inhibition duration with homogeneous F-actin disassembly: long inhibition time (Supplementary Fig. S6B).

Video 15. Effect of inhibition duration with heterogeneous F-actin disassembly: long inhibition time (Supplementary Fig. S6C).

Video 16. Spatiotemporal dynamics on a spherical cortical shell with underlying regime (i) dynamics and strong inhibition (s large) (Fig. 4B).

Video 17. Spatiotemporal dynamics on a spherical cortical shell with underlying regime (i) dynamics and weak inhibition (s small) (Fig. 4B).

Video 18. Experimental data from Bement et al. (2015); the bottom-left embryo is used for qualitative comparison with the model dynamics.

Video 19. Effect of pacemaker position on cortical reactivation dynamics (Supplementary Fig. S7).

Additional information

Funding

KU Leuven (Katholieke Universiteit Leuven) (GPUE/22/033)

• Daniel Cebrián-Lacasa