Modeling protein dynamics in Caenorhabditis elegans embryos reveals that the PLK-1 gradient relies on weakly coupled reaction-diffusion mechanisms.

SignificanceIntracellular gradients have essential roles in cell and developmental biology, but their formation is not fully understood. We have developed a computational approach facilitating interpretation of protein dynamics and gradient formation. We have combined this computational approach with experiments to understand how Polo-Like Kinase 1 (PLK-1) forms a cytoplasmic gradient in Caenorhabditis elegans embryos. Although the PLK-1 gradient depends on the Muscle EXcess-5/6 (MEX-5/6) proteins, we reveal differences in PLK-1 and MEX-5 gradient formation that can be explained by a model with two components, PLK-1 bound to MEX-5 and unbound PLK-1. Our combined approach suggests that a weak coupling between PLK-1 and MEX-5 reaction-diffusion mechanisms dictates the dynamic exchange of PLK-1 with the cytoplasm, explaining PLK-1 high diffusivity and smooth gradient.

Protein concentration gradients transmit spatial information and are crucial for cell and developmental biology. At the organismal scale, morphogen gradients can form over long spatial ranges, affecting several cells, and are critical for patterning cell identity across tissues (1–5). Morphogen gradients can form through source–sink mechanisms: a localized source generates diffusing molecules that move away and are eventually degraded (2, 6). At the cellular level, gradient formation is important for different aspects of cell division, including asymmetric cell division and cell-fate determination (7–11). At these smaller scales, the source–sink mechanism cannot generate a gradient, as diffusion would rapidly dissipate asymmetries (2, 12). Time- and space-dependent protein modifications that regulate a protein’s diffusion state sustain intracellular gradient formation (2, 12–14).

The one-cell Caenorhabditis elegans embryo is a ∼50-μm-long cell that provides a beautiful example of cytoplasmic gradient formation. Cell polarity is established shortly after fertilization, with the conserved PARtitioning defective proteins (PARs) localizing to the opposite embryonic poles: PAR-3, PAR-6, and PKC-3 accumulate at the anterior and PAR-1 and PAR-2 at the posterior cortices (15–23). Concomitantly, the RNA-binding protein MEX-5 (Muscle EXcess-5) and the redundant MEX-6 (hereafter referred to as MEX-5/6) segregate to the anterior cytoplasm, whereas the germline proteins PIE-1 and POS-1 accumulate in the posterior cytoplasm (24–29). Formation of the MEX-5 gradient depends on the dynamic switching between fast and slow diffusive states. At the posterior, phosphorylation by PAR-1 leads to an increase in MEX-5 mobility (fast MEX-5, MEX-5f); a uniformly distributed PP2A phosphatase counteracts PAR-1 phosphorylation by dephosphorylating MEX-5 (30). As a consequence, the majority of dephosphorylated MEX-5 will be at the anterior, increasing the probability of binding RNA complexes and resulting in a decreased MEX-5 mobility (slow MEX-5, MEX-5s) (29–31). The reaction–diffusion model developed by Lipkow and Odde (32) can be applied to explain the formation of MEX-5 cytoplasmic gradient in the one-cell embryo. The solution of differential equations describing the dynamics of the two MEX-5 species and the unevenly distributed reactions (phosphorylation and dephosphorylation) that determine the switch between the species demonstrate that a nonzero concentration gradient is established when the two diffusion coefficients (Dc) differ from each other (32). At the steady state, the relative concentrations of fast and slow MEX-5 along the embryo axis result in a gradient pattern mirroring the one of PAR-1 activity.

Formation of PIE-1 and POS-1 posterior gradients also depends on spatial regulation of their diffusivity, as they derive from the segregation of slow-diffusing complexes to the posterior cytoplasm (33, 34). The retention of POS-1 complexes is regulated by the Polo-Like Kinase 1 (PLK-1) (35). PLK-1 is an important mitotic regulator essential for many aspects of cell division (36–41). The cytoplasmic pool of PLK-1 accumulates in the anterior, and this enrichment depends on MEX-5/6 (42). After priming of their polo-docking site by the DYRK kinase MBK-2, MEX-5/6 bind by the polo-box domain of PLK-1, driving the accumulation of PLK-1 in the anterior. Depletion of MEX-5 reduces PLK-1 gradient, and the same result is obtained, to a lesser extent, with the depletion of the MEX-5 redundant homolog, MEX-6 (42). After depletion of both MEX-5/6, PLK-1 anterior enrichment is abolished (42). Despite the dependency of PLK-1 anterior enrichment on MEX-5/6, the biophysical mechanisms behind PLK-1 gradient establishment and how they compare to the dynamics of MEX-5 gradient formation are not known.

In this work, we combine experimental and modeling approaches to characterize how the MEX-5 and PLK-1 gradients form and evolve over the first asymmetric division in the embryo. We find that the PLK-1 gradient is shallower and forms later than the gradient of MEX-5. The two proteins also differ in their diffusivity, consistent with a substantial cytoplasmic pool of PLK-1 not bound to MEX-5. To characterize how the coupling between MEX-5 and PLK-1 reaction–diffusion mechanisms leads to patterning at the cellular scale, we developed an in silico approach in which the localization and interaction of MEX-5 and PLK-1 are modeled at the single molecule level. We present a three-dimensional computational model to study protein interactions and gradient formation of MEX-5 and PLK-1 over the first asymmetric cell division. The framework relies on Monte Carlo simulations, which take as inputs the measurements of protein diffusivity and reaction rates for both MEX-5 and PLK-1. Our simulations suggest that PLK-1 gradient formation is driven by weakly coupled reaction–diffusion mechanisms between PLK-1 and the slow MEX-5 component and recapitulate critical aspects of PLK-1 dynamics in the cytoplasm. This study provides a unique theoretical and quantitative model for understanding how dynamic interactions in the cytoplasm can give rise to complex patterning and can be extended to the reproduction of other polarity-related mechanisms in which cytoplasmic and cortical factors are involved.

Results

The Gradients of MEX-5 and PLK-1 Differ in Their Morphology and Dynamics.

PLK-1 becomes enriched in the anterior cytoplasm of one-cell embryos and, at the two-cell stage, is enriched in the anterior cell (called AB) relative to the posterior cell (called P1) (Movie S1) (42–47). PLK-1 anterior accumulation depends on MEX-5 (42). However, whether PLK-1 and MEX-5 gradient formation follows the same dynamics is not known. We measured and compared the morphology and evolution of the MEX-5 and PLK-1 cytoplasmic gradients. We used strains in which PLK-1 and MEX-5 were endogenously fused to mCherry (mCh::plk-1, this work, and mCh::mex-5, ref. 48). These fusion proteins recapitulate the anterior cytoplasmic enrichment of both PLK-1 and MEX-5 (Movie S1 ), and PLK-1 is properly recruited to different mitotic structures (i.e., kinetochores, centrosomes, nuclear envelope, and midbody), as previously shown (42, 49).

The average strength of MEX-5 and PLK-1 gradients became steeper during the first asymmetric cell division (Fig. 1) as a result of the anterior enrichment of both proteins over time. The average slope of MEX-5 gradient was at least twofold higher than that of PLK-1. This was confirmed in a plk-1::sgfp strain ().

Fig. 1.

PLK-1 and MEX-5 evolve with different dynamics. (A) Average absolute values of the gradient in the mCh::mex-5 and mCh::plk-1 strains, extrapolated at different stages of the first asymmetric cell division. The stages used for quantifications in this study are shown on the left of the y-axis. Statistical significance among the strains, for each separate stage, was obtained with two-tailed unpaired t test. (B) Time after which the maximum gradient is achieved for early stage embryos. Statistical significance was obtained with two-tailed unpaired t test. (C) Average experimental curves for MEX-5 and PLK-1 gradient as a function of time, for the recorded early embryos. The number N of analyzed embryos is reported next to each curve. The curves were linearly fitted for the time gap (8 to 70) frames for mCh::mex-5 and (35 to 90) frames for mCh::plk-1, and the absolute value of the linear coefficient m from the linear is reported in the graph. (D) Relative increase of protein intensity in AB with respect to P1, calculated using . Statistical significance was obtained with unpaired two-tailed Mann–Whitney U test. In the bar graphs of all main figures, bars represent mean values, dots represent single embryo measurements, and error bars represent the SDs. The number of analyzed embryos is specified for each condition. ns: nonstatistically significant; *P < 5 × 10−2; **P < 10−2; ***P < 10−3; and ****P < 10−4. The statistical analysis for all figures is summarized in .

PLK-1 took a longer time than MEX-5 to form the maximum gradient (refer to Fig. 1 and for the plk-1::sgfp strain). The gradient of PLK-1 was indeed negligible (comparable to zero) until the ruffling stage, suggesting that PLK-1 anterior enrichment was delayed compared to MEX-5 (Fig. 1 and ). To confirm this, we plotted the average experimental values of MEX-5 and PLK-1 gradients obtained from recordings of early embryos. While MEX-5 gradient increased steeply (in absolute values) from the beginning of the process, PLK-1 showed a flat plateau around the value of −0.045 until timeframe 40 (refer to Fig. 1 and for the plk-1::sgfp strain). We quantified the steepness of the curves in their decreasing part (Fig. 1), and the values, extrapolated as the linear coefficient m from the fit (black lines and m values in Fig. 1), showed that PLK-1 anterior enrichment was not only delayed but also characterized by slower dynamics compared to MEX-5.

A less-steep gradient should result in a more homogeneous protein distribution in the daughter cells. We next asked whether the difference in the strength of the MEX-5 and PLK-1 protein gradients was reflected at the two-cell stage. PLK-1 concentration in AB was about 30% higher than in P1 (Fig. 1 and ), whereas MEX-5 concentration in AB was almost 75% higher than in P1 (Fig. 1), consistent with the steepness of the gradients measured in one-cell embryos.

Our data indicate that although PLK-1 binds to MEX-5 (42), its gradient initiates later and is significantly weaker than the gradient of MEX-5, suggesting that a significant fraction of cytoplasmic PLK-1 does not follow MEX-5 dynamics.

PLK-1 Diffuses Faster than MEX-5.

The behavior of MEX-5 and PLK-1 gradients suggests that the biophysical properties of these proteins are different. To estimate their overall mobility in the anterior and posterior cytoplasmic compartments, we used Fluorescence Recovery After Photobleaching (FRAP) microscopy ().

We fitted FRAP recovery curves of MEX-5 to a one-component model as such for a single diffusion species to extrapolate an apparent Dc. While a one-component model does not reflect the presence of fast- and slow-diffusing species, it provides a description of the overall mobility of MEX-5. At prepolarization, MEX-5 Dc was not significantly different between anterior and posterior (Fig. 2 and ). At nuclear envelope breakdown (NEBD), when polarity establishment is completed in control embryos, MEX-5 Dc in the posterior had nearly doubled, whereas the value of the MEX-5 Dc in the anterior remained similar to prepolarization (Fig. 2 and ), consistent with published data using transgenes (24, 30). This indicates that MEX-5 gradient formation coincides with an increase in MEX-5 mobility in the posterior.

Fig. 2.

Measurements of MEX-5 and PLK-1 apparent diffusion coefficient (Dc) in control and RNAi-treated embryos. MEX-5 (A) and PLK-1 (B) Dc at prepolarization stages (light blue) and at polarization (NEBD, orange). Measurements were performed at the anterior and posterior of the embryo. (C) PLK-1 Dc after MEX-5/6, MBK-2, and PAR-1 depletion compared to control (empty vectors, L4440 on the left of the dotted bar and T444T on the right; ). The measurements were performed at the anterior and posterior of the embryos at NEBD. The statistical significance in all panels was obtained with two-way ANOVA test, with Tukey’s multiple comparison. In C, only the significance among the PLK-1 diffusivity was tested. ns: nonstatistically significant; *P < 5 × 10−2; **P < 10−2; ***P < 10−3; and ****P < 10−4. The significance levels discussed in the text are reported in the figure. The full results can be found in .

Similar to MEX-5, FRAP measurements of PLK-1 showed that its Dc was not statistically different between anterior and posterior at prepolarization. The absolute values of PLK-1 Dc () were threefold higher than those of MEX-5 (Fig. 2 and ). At NEBD, the anterior PLK-1 Dc had decreased significantly, whereas the posterior Dc remained similar to prepolarization (Fig. 2 and ). Therefore, PLK-1 gradient formation coincides with a decrease in PLK-1 mobility in the anterior. The Dc of PLK-1 was higher than the one of MEX-5 also at NEBD. PLK-1 Dc at NEBD measured in the plk-1::sgfp strain did not show a statistically significant difference from the mCh::plk-1 strain ().

These results show that at NEBD, the average diffusivity of both MEX-5 and PLK-1 is higher in the anterior than in the posterior. There are two major differences in the behavior of these proteins: 1) PLK-1 diffusivity is higher than MEX-5, and 2) polarization corresponds to an increase in MEX-5 mobility in the posterior and to a decrease in PLK-1 mobility in the anterior.

Binding to MEX-5/6 Slows Down PLK-1 Diffusion.

We then asked how PLK-1 dynamics are regulated by MEX-5. We performed FRAP measurements of PLK-1 Dc at the NEBD, after different RNA interference (RNAi) conditions. Consistent with the results in Fig. 2 , the diffusivity of PLK-1 was higher than the one of MEX-5 (Fig. 2 and ). The difference in the absolute values compared to Fig. 2 was due to the different growth condition ().

The values of PLK-1 Dc in mex-5/6(RNAi) embryos, in which PLK-1 gradient is abolished (42), were similar in the anterior and posterior at NEBD (Fig. 2 and ) and almost twice the values of the anterior Dc in control.

We then measured PLK-1 diffusivity after depleting the DYRK kinase MBK-2 (42). MBK-2 phosphorylates MEX-5/6, thereby priming the interaction with PLK-1 (42). In mbk-2(RNAi) embryos, the anterior and posterior values of the Dc were not significantly different between each other (Fig. 2 and ) and not significantly different from the Dc measured in mex-5/6(RNAi) embryos (). This indicates that the interaction of PLK-1 with MEX-5/6 limits the Dc of PLK-1.

We then characterized PLK-1 dynamics in par-1(RNAi) embryos. Since PAR-1 increases MEX-5 mobility (30), MEX-5 is symmetrically distributed with uniformly low mobility in par-1(RNAi) and par-1 mutant embryos (24, 25, 30, 50). To achieve strong PAR-1 depletion, we performed par-1(RNAi) in a mCh::plk-1; par-1(zu310ts) strain at restrictive temperature (). In the control, the Dc of MEX-5 and PLK-1 was higher in the posterior than in the anterior (). PAR-1 depletion resulted in a uniform PLK-1 Dc along the anterior–posterior axis (Fig. 2 and ), consistent with published data showing a homogeneous PLK-1 localization (42, 43, 45). The Dc values of PLK-1 in PAR-1-depleted embryos are lower than in either mex-5/6(RNAi) or mbk-2(RNAi) embryos. In the absence of PAR-1 in which the majority of MEX-5 was slow (30), if the majority of PLK-1 was bound to MEX-5, PLK-1 Dc should have been more similar to the anterior Dc in control. On the contrary, the Dc of PLK-1 was more similar to the Dc of posterior PLK-1 in control. This suggests that in par-1(RNAi) embryos, the unbound PLK-1 pool strongly contributes to the diffusivity or that the binding of PLK-1 to MEX-5 is reduced ( and ).

These results are consistent with a model in which the binding between PLK-1 and MEX-5s is primed by MBK-2 and reduces PLK-1 mobility in the anterior.

Monte Carlo Simulation of MEX-5 Dynamics and Gradient Formation.

To understand how the MEX-5 and PLK-1 reaction–diffusion mechanisms are intertwined, we developed an in silico model of their dynamics during gradient formation in the one-cell embryo. We started by modeling the formation of the MEX-5 gradient, since this has extensively been described biologically and mathematically (24, 30). Compared to previous models, we introduced a time-resolved analysis of how the concentration profiles evolve during the whole first asymmetric cell division.

We used an approach that relies on Monte Carlo simulations, which allow modeling of the evolution of stochastic phenomena in both space and time at the single particle level. Detailed description of the code implementation is in . To reproduce in silico both the biological and physical mechanisms behind the establishment of gradients, protein dynamics was followed in a three-dimensional geometry, updating the protein state and localization at each simulated timepoint dt. We implemented in the simulations a software replica of the chosen detection technique (time-lapse fluorescence microscopy) and of the image analysis pipeline ().

The approach combined the physical properties of the protein (diffusivity) and biochemical reactions that are known to change MEX-5 state (phosphorylation/dephosphorylation).

We used the experimental apparent MEX-5 Dc obtained from our FRAP studies to calculate the diffusivity of the MEX-5f and MEX-5s species (Fig. 2 and ). The obtained values ( μm2/s and μm2/s) were in agreement with data from fluorescence correlation spectroscopy (30). We used these values to tune MEX-5 particle velocity in the embryo (): the velocity for MEX-5s was set to = 0.45 μm/s and the one for MEX-5f to = 4.9 μm/s.

To check for the validity of this calculation, we extracted the overall mobility of MEX-5 as a function of time and position along the embryo axis (Fig. 3). We obtained from the simulation the average velocity of MEX-5, calculated as the mean of the MEX-5f and MEX-5s particles in a specific voxel. The simulated average velocity of MEX-5 linearly increased at the posterior pole from prepolarization (red) to steady state (blue), leading to a twofold ratio between posterior and anterior (Fig. 3 and ), consistent with the ratio obtained from the measurement of the Dc in Fig. 2 and .

Fig. 3.

Comparison between simulation and experimental results of MEX-5 gradient and mobility. (A) Average simulated MEX-5 velocity as a function of the position along the embryo axis in the central volume slice Δz. Different timepoints are shown, as indicated by the color legend, from prepolarization (red) to maximum gradient (light blue). (B) Sensitivity analysis on the reaction rates and . Comparison of MEX-5 gradient steepness as a function of time when reaction values are varied from the published ones (30). (C) Comparison of MEX-5 experimental results (curve in Fig. 1, translated to zero by using the average of the values from timepoints 0 to 8, green) with the best-fitting simulation results, obtained using = 0.005 s−1 and = 0.001 s−1 to 0.0055 s−1 (violet). The points represent the experimental and simulation data. The fits obtained using the analytical three-parameter function reported in Eq. are shown (solid lines). In A–C, results are average values from 10 simulation runs. In B and C, the shaded area represents the SD. (D) Relative difference between experimental (green) and simulated (violet) gradients, calculated using Eq. 2, as a function of time. (E) Concentration of slow (blue), fast (red), and total (green) MEX-5 as a function of the position along the embryo axis at the steady state. The concentration was calculated as the number of particles of each MEX-5 species in each voxel of the central volume slice Δz (Movie S2) over the total number of simulated particles. Results are average values from 20 simulation runs, and the shaded area represents the SD.

We introduced MEX-5 phosphorylation by PAR-1 and dephosphorylation by a PP2A phosphatase (24, 30, 32) by means of the kinase rate and the phosphatase rate (). Concerning , we took into account a time-delayed gradient formation for PAR-1 activity throughout the cytoplasm ().

To test whether the values of the and rates influenced the characteristics of the MEX-5 gradient, we tuned these parameters with respect to published values by Griffin et al. (30) and measured how the variations influenced the model (parameter sensitivity analysis). The final gradient steepness did not vary significantly if the parameters (including the lower and upper limits for ) were scaled by several factors of 10 (1/10 pink, 1/20 violet, and 1/30 green in Fig. 3). When all the reaction rates were decreased by a factor of 50 (light blue), the dynamics started to slow down considerably. When the reduction was of two orders of magnitude (dark blue), the change in gradient strength monotonously decreased over the investigated time range, without reaching a steady-state value (Fig. 3). In conclusion, the absolute value of MEX-5 gradient at the steady state depends on the ratio between and , and changing their values accordingly only modifies the dynamics of the process.

We then compared the experimental values (green in Fig. 3) to the simulation curve obtained using 1/20th of the reference and rates (30) (i.e., linearly increasing from 0.001 s−1 to 0.0055 s−1 and = 0.005 s−1) (violet in Fig. 3 ). The simulation converged to the experimental values, with a maximum gradient of −0.45 obtained around 70 timeframes, corresponding to 10 min. These results agree with the observations in Fig. 1 .

To further characterize how well our simulation of MEX-5 gradient formation corresponds to the observed MEX-5 gradient formation, the plots in Fig. 3 were fitted with the following analytical functions (Fig. 3, continuous lines):where t is time and A, B, and C the free parameters extrapolated from the fit. The values extrapolated from the fit along with a 95% CI and the goodness of the fit are reported in . The analytical function in Eq. reproduced the trend of the experimental data. From the fits, we extracted the parameter A (negative asymptote) as indicator of the absolute value of the maximum gradient at steady state. This was 0.44 0.01 for the experimental curve and 0.43 0.01 for the simulations. We then quantified the relative difference between the fit curves for the simulation (violet) and the experimental (green) datapoints over time using the formula:

We considered the experimental values from the fit curve as the reference, and therefore, the relative difference was equal to zero (Fig. 3). For early timepoints, there was 25% difference between experiments and simulations. This was explained by the variability scored among imaging experiments and because the simulations needed time to reach an equilibrium. The difference steadily decreased over time to reach values of the order of 5% (starting from timeframe 40 to 50), showing agreement between experimental data and the in silico behavior. We concluded that the reaction rates = 0.001 s−1 to 0.0055 s−1 and = 0.005 s−1 can reproduce the MEX-5 gradient and they were therefore kept fixed for the following simulation runs.

To test whether we reproduced the characterized properties at steady state (24, 30) under these initial conditions, we analyzed the contributions of MEX-5s and MEX-5f to total MEX-5 gradient. MEX-5s displayed a 3-fold gradient from anterior to posterior (concentration of voxel−1 at the anterior versus voxel−1 at the posterior; Fig. 3). MEX-5f was instead preferentially produced in the posterior, but it quickly diffused to create a homogenous cytoplasmic background (Fig. 3). At the anterior, the concentration of MEX-5s was three times higher than the concentration of MEX-5f, while at the posterior the ratio slow-to-fast MEX-5 dropped to 1:1. For this reason, the computations showed that accumulation of MEX-5s in the anterior underlay the MEX-5 asymmetry, with an almost twofold gradient ( voxel−1 at anterior versus voxel−1 at posterior; Fig. 3). Simulated time-lapse movies of the distribution of total MEX-5, MEX-5s, and MEX-5f components are shown in Movie S2 .

In conclusion, we have established a time-resolved Monte Carlo simulation that can successfully reproduce the process of MEX-5 dynamics and gradient formation in silico.

Simulation of the PLK-1 Reaction–Diffusion Mechanism.

To address how PLK-1 and MEX-5 reaction–diffusion mechanisms are coupled, we extended our Monte Carlo approach to include PLK-1. PLK-1 particles were simulated with initial velocity based on PLK-1 Dc at prepolarization ( = 5.3 μm/s), under the assumption that at that stage, PLK-1 is mainly unbound from MEX-5. It is not known whether PLK-1 preferentially binds the fast or slow component of MEX-5 nor how stable the interaction is. We investigated whether different binding and unbinding kinetics impacted PLK-1 gradient formation. We allowed PLK-1 to either bind to MEX-5s, MEX-5f, or both (MEX-5f+s) (). The interaction was stochastic according to a binding rate and such that once bound to MEX-5, PLK-1 inherited the velocity of the MEX-5 partner (). We allowed the proteins to dynamically unbind through one of two mechanisms (). The first mechanism was regulated by a constant off-rate , which we will refer to as “free unbinding.” It takes into account a bidirectional, dynamic exchange between MEX-5 and the PLK-1 cytoplasmic pools (, blue panel). In the second mechanism, PLK-1 dissociated from MEX-5 each time that MEX-5 switched between the fast- and slow-diffusing states as regulated by PAR-1 and PP2 phosphorylation/dephosphorylation (hereafter referred to as “imposed unbinding”; , violet panel). This mechanism considers a scenario in which there is a conformational change in MEX-5 that causes unbinding of PLK-1.

We established a model in which the interaction between PLK-1 and MEX-5 is dictated by the rates and :

We tested different binding modalities: the binding rate was constant from anterior to posterior at prepolarization, due to the uniform concentration of both MEX-5 components before polarity establishment (0.045 s−1). Then, as a function of time, increased at the anterior (0.08 s−1) while depleting the posterior (0.03 s−1), coherently with MEX-5 gradient formation (). When PLK-1 was only allowed to bind MEX-5f, it was also homogeneously distributed (green curve in Fig. 4). When PLK-1 was only allowed to bind MEX-5s, PLK-1 was enriched in the anterior and formed a gradient with a steepness similar to the observed gradient (Fig. 4, red curve). The last scenario, in which PLK-1 could bind indiscriminately to MEX-5f+s, did not result in a significant difference compared to binding to MEX-5s only (Fig. 4, black curve). Since there are no data indicating preferential binding between PLK-1 and MEX-5s, we considered for our subsequent simulations the scenario in which PLK-1 bound to MEX-5f+s.

Fig. 4.

Comparison of PLK-1 simulations and experimental curve suggests the existence of a large fraction of unbound PLK-1 and a delay in binding to MEX-5. (A) Results for PLK-1 gradient formation as a function of time, simulated for the three different binding scenarios presented in . Both unbinding modalities were allowed, with an off-rate = 0.1s−1. (B) Dependency of PLK-1 gradient on the free unbinding rate . Simulations were performed by using settings in which binding to MEX-5f+s and both unbinding modalities (free and imposed) were allowed. The unbinding rate was tuned within a range from to 0.00 s−1 (purple) to 0.50 s−1 (blue). The average PLK-1 experimental observation (translated to zero using the average of the values from timepoints 0 to 35 in all panels, orange) is reported as reference. The points represent the experimental and simulation data, and the shaded area represents the SD from the experiments and from 10 simulation runs. The fits obtained using the analytical three-parameter function reported in Eq. are shown (solid lines). (C) Quantification of steady-state concentrations of unbound PLK-1 (light blue) and PLK-1 bound to MEX-5s (blue) and to MEX-5f (red) when binding to MEX-5f+s species was allowed. Total PLK-1 concentration (green) along the embryo axis is also shown. The concentrations are calculated as number of particles of each PLK-1 species in each voxel of the central slice Δz (Movie S3) over the total number of simulated particles. (D) Average simulated PLK-1 velocity as a function of the position along the embryo axis. Different timepoints are shown, as indicated by the color legend, from prepolarization (red) to steady state (light blue). (E) Results for PLK-1 gradient formation when the percentage of MEX-5 available for binding to PLK-1 was increased linearly from 0 to 100% in the time gap T1 to T2 (reported in the plot). The average experimental curve is reported for comparison. (F) Results for PLK-1 gradient formation when combining an initial time lap with no binding (from 0 to T1) between the two proteins and a linear increase in MEX-5 availability to bind (from T1 to T2, as reported in the plot). The average experimental curve is reported for comparison. In all panels, simulation results are average values from 10 simulation runs, and the shaded area represents the SD.

We performed a sensitivity analysis on : since the strength of the coupling between the MEX-5 and PLK-1 reaction–diffusion mechanism is unknown, we tuned the value of the free unbinding rate (Fig. 4). When no free unbinding is allowed ( = 0 s−1), PLK-1 unbinding from MEX-5 only occurred when MEX-5 switched diffusive states, as regulated by the PAR-1 and PP2 reaction rates (imposed unbinding). In this scenario of almost “perfect coupling,” the unbound PLK-1 pool was depleted, leading to PLK-1 dynamics similar to MEX-5 and thus to a similar gradient (purple curve, Fig. 4). The simulated PLK-1 gradient was significantly stronger than the experimental one, suggesting that in embryos, the PLK-1 and MEX-5 gradients were not perfectly coupled. With increasing , PLK-1 gradient became smoother, as the pool of unbound PLK-1 increased. The curves in Fig. 4 were fitted using Eq. and compared to the experimental curve for PLK-1 (orange, Fig. 4). The values of the free parameters A, B, and C extrapolated from the fit, along with the 95% CI and the goodness of the fit, are reported in . The comparison of the fit parameters A (negative asymptote, ) for the different conditions led to the conclusion that the reproduction of PLK-1 gradient at steady state required a significant amount of PLK-1 unbound from MEX-5. Given the chosen value, this pool was sufficiently replenished only with a uniform = 0.1 s−1. This parameter was fixed as the reference condition for our subsequent simulations.

Under these conditions, we then asked which pool of PLK-1 created the homogeneous background that flattens its concentration gradient. The simulations predicted the unbound pool of PLK-1 to be uniformly distributed in the cytoplasm and its abundance did not vary during the first division (∼6 10−5 voxel-1, , Upper). The concentration of PLK-1 bound to MEX-5s increased during gradient formation in silico at the anterior from prepolarization to polarization (from 1.8 10−5 to 4 10−5 voxel−1), while it dropped to 1 10−5 voxel−1 in the posterior (, Middle). The concentration of PLK-1 bound to MEX-5f was, instead, uniformly distributed and increased from prepolarization to polarization of 67% (from 6 10−6 to 1 10−5 voxel−1; , Lower). From this analysis, we concluded that the concentration of unbound PLK-1 (light blue, Fig. 4) was higher with respect to the one of the bound pools: it was almost twice the concentration of PLK-1 bound to MEX-5s (blue curve, Fig. 4) at the anterior and more than four times the amount bound to MEX-5s and MEX-5f (red, Fig. 4) at the posterior. The evolution of the concentration profiles for total PLK-1 and its different pools over time is shown in Movie S3 . These results show that the coupling between PLK-1 and MEX-5 is weak.

Before polarization, the relative concentrations of unbound PLK-1 (with = 5.3 μm/s) and PLK-1 bound to the uniformly distributed MEX-5s (with = 0.45 μm/s) resulted in an average PLK-1 velocity of 4.2 to 4.3 μm/s throughout the embryo (Fig. 4 and ). This predicted velocity of PLK-1 was higher than that of MEX-5 at the same stage, in agreement with the FRAP measurements (Fig. 2 ). During polarization, PLK-1 velocity decreased in the anterior due to the preferential binding to MEX-5s and slightly increased at the posterior ().

To conclude, the slow component of MEX-5 is essential for PLK-1 gradient formation, and the PLK-1 gradient is properly established by assuming a two-component model (a pool of PLK-1 binding to MEX-5 and one unbound from it). The two compartments are not stagnant, as a continuous dynamic exchange between PLK-1 species, the MEX-5-unbound and the MEX-5-bound, is required to recapitulate the observed PLK-1 gradient dynamics.

Modeling the Delay in PLK-1 Gradient Initiation.

The simulations presented so far reproduced PLK-1 gradient steepness and PLK-1 diffusivity at steady state but failed to reproduce the initial delay in PLK-1 gradient establishment (Fig. 1). We hypothesized that PLK-1 binding to MEX-5 is controlled over time (depending, for example, on MBK-2 phosphorylation of the polo-docking site of MEX-5) (42). In this scenario, the binding could either increase as a function of time or even happen discontinuously. To computationally address the source of this plateau, we modeled two scenarios:

We linearly increased the percentage of MEX-5 available for the binding, from 0 to 100%, in precise time gaps (Fig. 4). These conditions increasingly slowed down the dynamics of PLK-1 gradient formation but did not reproduce the almost two-step function observed experimentally (orange in Fig. 4).

We imposed a fixed time delay of around 35 frames (350 s) during which binding was not occurring before applying the algorithm above in point 1. By doing this, the simulations reproduced the observed initial plateau in gradient formation. After this time gap, we let the MEX-5 availability for binding increase to 100% for different amounts of time (Fig. 4), as in the former scenario, to vary the velocity with which the final PLK-1 concentration gradient was established.

All the data points in Fig. 4 were fitted with the analytical function in Eq. , and the free parameters A, B, and C extracted from the fit are reported in . To compare the overall temporal evolution of gradient formation among these fit curves, we compared the flex points, defined as the time frames at which the curves reached 10% and 90% of the maximum gradient (plateau value 0.15 for simulation datasets, and 0.17 for the experimental one) (). We considered these frames to be the points at which the bending to the decreasing trend and the steady-state plateau were respectively reached. The best fit of the experimental results was obtained when T1 = 35 and T2 = 80 (350 s and 800 s, respectively; Fig. 4 and ), for a total time gap of 450 s. We finally compared this curve with the one fitting the average experimental results over time by calculating the relative difference (Eq. ). The simulation results displayed good agreement after timeframe 50 (), while the initial large variability was explained by how the relative difference was calculated (division for values close to zero) and was not representative of the real discrepancy between the fit curves, which were instead in very good agreement in the initial plateau phase (Fig. 4).

These results show that the proposed Monte Carlo code allows full reproduction of the PLK-1 gradient, not only at steady state but also in its full temporal dynamics. The biological findings have been interpreted by tuning two main parameters: 1) the free unbinding rate (best fitting value 0.1 s−1) and 2) the time gap for the delay in MEX-5 priming (best fit for T1= frame 35) and the subsequent gap for its linear increase (till T2= frame 80).

Prediction of Gradient Behavior after Varying Protein Concentrations.

Our modeling approach described the behavior of both MEX-5 and PLK-1 during gradient formation, assuming that the number of MEX-5 and PLK-1 molecules was the same. We next tested whether the simulation predicted MEX-5 and PLK-1 gradient characteristics when the concentration of MEX-5 relative to PLK-1 was changed. To achieve this, we integrated the relative concentration of the two proteins in each volume voxel to determine the likelihood of binding. The decision of whether PLK-1 bound still happened stochastically () but depended on the availability of at least one MEX-5 unbound protein in the interrogated voxel. If all MEX-5 particles in the voxel were bound to PLK-1, no further binding was allowed ().

We first simulated the MEX-5 gradient when varying the number of MEX-5 particles. The results are shown in Fig. 5. The gradient of MEX-5 remained unaltered when its concentration was either increased or decreased till a factor of 10, indicating that MEX-5 gradient formation is relatively insensitive to MEX-5 concentration. For decreasing concentrations, the gradient steepness progressively weakened, as a result of the low MEX-5 particle concentration inside the embryo volume. This resulted in a homogeneous but sparse distribution of MEX-5 for a very low number of particles (103, gray).

Fig. 5.

Modeling can predict gradient and diffusivity behaviors after change in protein concentration. (A) Simulation of MEX-5 gradient for different initial MEX-5 particle numbers. (B) Simulation of PLK-1 gradient for different relative concentrations of MEX-5 with respect to PLK-1, shown in the plot as “MEX-5:PLK-1 ratio.” The run-down condition 1:100 was achieved by increasing of an order of magnitude PLK-1 particle number, not to lose statistics on MEX-5 gradient. (C) Relative difference of MEX-5 intensity between AB and P1 in the mCh::mex-5 strain, normalized to AB, after mex-5(RNAi). Statistical significance was obtained with unpaired two-tailed Mann–Whitney U test. (D) Average experimental curves for PLK-1 gradient as a function of time, for the recorded early embryos after mex-5(RNAi), in the mCh::plk-1 strain. The number N of analyzed embryos is reported next to each curve. (E) Simulation of MEX-5 and PLK-1 gradients following PAR-1 depletion. A 30:70 MEX-5f:MEX-5s concentration ratio was used. In A, B, and E, the points represent the simulation data, and the shaded area represents the SD from 10 simulation runs. (F) PLK-1 simulated average velocity following PAR-1 depletion at the steady state, as a function of the position in the embryo, and comparison with the control condition (data elaborated from Fig. 4). Points represent the average velocity values at the anterior and at the posterior from single runs. The bars represent the mean and the errors represent the SD from the 10 simulation runs. ns: nonstatistically significant; *P < 5 × 10−2; **P < 10−2; ***P < 10−3; and ****P < 10−4.

We then probed how different MEX-5 concentrations influenced the PLK-1 gradient in silico. Increasing MEX-5:PLK-1 ratio till a factor of 10 did not disrupt PLK-1 gradient formation (Fig. 5). Reducing the MEX-5:PLK-1 ratio to 1:2, 1:5, 1:10, and 1:100 caused progressive weakening of the PLK-1 gradient (Fig. 5) even though the MEX-5 gradient was not weakened at ratios of 1:2, 1:5, and 1:10. These results argue that there is a threshold level of MEX-5 required to establish the PLK-1 gradient, which we predicted to be CMEX-5 > 1/2 CPLK-1.

To experimentally test these predictions, we used RNAi to partially deplete MEX-5 in the mCh::mex-5 and mCh::plk-1 strains. Partial depletion of MEX-5 was highly variable, as scored by quantifying the intensity of mCh::mex-5. For our analysis of the MEX-5 and PLK-1 gradients, we grouped embryos from the condition in which MEX-5 was depleted between 45 and 75% compared to control (). We then quantified the asymmetry of MEX-5 at the two-cell stage. Consistent with the modeling results, MEX-5 asymmetry was decreased by partial MEX-5 depletion (Fig. 5).

Quantifications of the PLK-1 gradient showed that beginning at pseudocleavage, the PLK-1 gradient was increasingly weaker in mex-5(RNAi) embryos, with a final steepness of almost 1/2 with respect to control embryos (Fig. 5 and ). Consistent with a weakened gradient of PLK-1 prior to cytokinesis, PLK-1 asymmetry at the two-cell stage was reduced up to a 48% in the mCh::plk-1 strain ().

Our modeling approach can therefore predict the behavior of the PLK-1 gradient when the ratio between MEX-5 and PLK-1 is decreased.

We then used simulations to guide the interpretation of the results on the diffusivity of PLK-1 after par-1(RNAi) (Fig. 2). Consistent with the absence of gradient formation, the diffusivity of PLK-1 was constant in PAR-1-depleted embryos but was more similar to the posterior value of the control than the anterior (Fig. 2), despite the fact that the majority of MEX-5 is in the dephosphorylated state (MEX-5s) in this condition (30).

We simulated par-1(RNAi) condition by setting the value of kPAR-1 = 0 s−1. Since Griffin et al. (30) have estimated MEX-5s being about 70% of total MEX-5 also in the par-1 mutant, we assumed this concentration ratio within the embryo. The results concerning gradient formation for both MEX-5 and PLK-1 are shown in Fig. 5. The gradients of MEX-5 and PLK-1 were not established, with average gradient values around 0. This was expected, since PAR-1 asymmetric contribution to the phosphorylation of MEX-5 is absent, in agreement with published data on par-1(RNAi) and mutants (e.g., refs. 24, 25, 30, 51) and our data on the constant Dc (Fig. 2). At steady state, in par-1(RNAi), the average velocity of PLK-1 extracted from the simulations was uniform in the embryo and assumed the value of 4.2 to 4.3 μm/s (Fig. 5 and ). This value was higher than the anterior value of control embryos and was due to the contribution of MEX-5f to the overall velocity of PLK-1 particles. This result was in agreement with the experimental observations in Fig. 2 ().

To conclude, our modeling can predict the behavior of the gradient when the concentration of one component is increased or decreased or when polarity is affected.

Discussion

In this work, we have combined microscopy measurements and Monte Carlo simulations to quantify PLK-1 and MEX-5 gradient establishment during the first asymmetric cell division of the C. elegans embryo. Our data reveal that PLK-1 gradient formation does not follow the same dynamics and characteristics of the MEX-5 gradient, despite the fact that it depends on MEX-5. Specifically, we find that the PLK-1 gradient is less steep and forms later than the gradient of MEX-5. Our experimental and computation results support a weak coupling between the MEX-5 and PLK-1 gradients due to a substantial pool of PLK-1 that is unbound from MEX-5.

We find that the apparent diffusion coefficient of PLK-1 follows the same trend as MEX-5, with a Dc at prepolarization not significantly different between anterior and posterior. At steady state, once polarity has been established, the Dc of PLK-1 is lower at the anterior, consistent with the trend of MEX-5, which accumulates at the anterior in its dephosphorylated slower form. However, the Dc of PLK-1 is always higher than the one measured for MEX-5, and PLK-1 gradient coincides with a decrease of the diffusivity in the anterior.

To explain the difference in gradient steepness and in the diffusivity of the two proteins, we propose a 2-component model for PLK-1: the first component is represented by a pool binding to MEX-5, the second one is not bound to MEX-5 and has its own Dc throughout the cytoplasm.

We integrated the experimental observations with Monte Carlo simulations within a modeling framework that is able to reproduce protein interactions and gradient formation during the whole process of cell division. It furthermore recreates the experimental setup used for the readout, along with the strategy adopted for the image processing and analysis. We used this framework to test different biological scenarios concerning the binding with MEX-5 (MEX-5s, MEX-5f, MEX-5f+s) and the unbinding modalities (free and imposed unbinding).

Our simulations suggest that the critical parameters for PLK-1 gradient formation are the binding with the slow component of MEX-5 (Fig. 4) and the dynamic exchange between the unbound and MEX-5 bound pools to PLK-1 (free unbinding rate of 0.1 s−1) (Fig. 4 ). With this, we show that the coupling between PLK-1 and MEX-5 reaction–diffusion mechanisms has to be rather weak to guarantee a smooth PLK-1 gradient (): PLK-1 therefore continuously and dynamically exchanges between the bound and unbound pools, as dictated by the rates and .

This hypothesis is consistent with the fact that PLK-1 is recruited to mitotic structures at different phases during mitosis and with the existence of a pool that is highly dynamic and available to perform other functions than the one related to polarity establishment.

Simulation results also predict that the bound component is less abundant than the unbound PLK-1, with unbound PLK-1 being 60% of total PLK-1 concentration at the anterior and 78% at the posterior at steady state (Fig. 4). Indeed our simulations predict that when PLK-1 is perfectly coupled to MEX-5 (i.e., no unbound PLK-1), PLK-1 gradient matches MEX-5's and it is significantly stronger than the one experimentally measured. This demonstrates that binding to MEX-5 is the main driving mechanism to PLK-1 gradient formation, and that any other mechanisms regulating PLK-1 cytoplasmic localization should result in either a gradient similar to the one of MEX-5 or in a uniform distribution.

We moreover show in this work that the temporal evolution of the PLK-1 gradient is different from the one of MEX-5. The formation of PLK-1 gradient initiates later, and the curve displays an initial plateau with gradient values close to zero (Fig. 1). We hypothesize that either binding of PLK-1 to MEX-5 increases linearly or that PLK-1 is not able to bind to MEX-5 for a specific amount of time. Simulations show that when PLK-1 is not able to bind MEX-5 for the first 350 s, the gradient is comparable to the experimental data. This suggests that in vivo, the interaction of PLK-1 and MEX-5 cannot occur immediately and is consistent with the fact that MEX-5 is primed by the kinase MBK-2 to interact with PLK-1. MBK-2 activation could be delayed in time after symmetry-breaking, likely taking time to establish a pool of primed MEX-5. In the future, it will be interesting to further test this finding by constructing a MEX-5 mutant that can constitutively bind PLK-1 and ask whether this makes the timing of PLK-1 gradient formation similar to MEX-5.

The proposed simulation framework can be easily extended for the reproduction of other cytoplasmic polarity-related factors by accordingly changing boundary conditions for the protein velocity and localization. While the distribution of different factors is known to persist after establishment of the steady state, many of the physical rationales behind their accumulation remains unknown. This tool offers the possibility to guide the formulation of biological hypotheses, possibly underlying complex cellular mechanisms, that are hard to test within an in vivo system. The foreseen scenarios can eventually drive the design and interpretation of experiments to target specific biological questions in a very iterative process. For example, we have used the modeling to change the relative concentration between MEX-5 and PLK-1. This has revealed that reducing MEX-5 results in a less-steep gradient of PLK-1, something that was confirmed by the experimental data. We moreover simulated the condition of depleted PAR-1 to interpret our FRAP results following this treatment. We found that the Dc of PLK-1 in the absence of PAR-1 is not different between anterior and posterior, consistent with the fact that the gradient does not form. However, we found that the Dc is high when compared to the value of the anterior Dc in control embryos. This was an unexpected result, since in this condition, MEX-5 is slow. Griffin et al. (30) have shown that even in the absence of PAR-1, the ratio between fast and slow component of MEX-5 remains 30:70. However, in par-1 mutants, the 70% slow component will not be enriched at the anterior but distributed throughout the cytoplasm. Indeed, by simulating this situation with a 30:70 ratio of MEX-5 fast and slow component throughout the embryo, the obtained velocity of PLK-1 was higher than the anterior one in control. One possibility to explain this result in vivo is that the amount of PLK-1 bound to the 70% MEX-5 slow, which is distributed throughout the embryo, may not be sufficient to overcome the Dc of the unbound PLK-1. Another not-mutually exclusive possibility is that in the absence of PAR-1, PLK-1 does not bind efficiently to MEX-5, resulting in an increase of the unbound pool and, therefore, in a Dc similar to the posterior rather than the anterior one of control embryos. Being able to measure PLK-1 binding to MEX-5 in one-cell embryos in control and polarity mutants using single-cell biochemistry (52) may help to answer this question.

Reproduction of events at small spatial and temporal scales via Monte Carlo allows to shed light on system properties that are otherwise undetectable, due to, for example, fast dynamics and intricated interaction networks. Our approach can be easily extended to reproduce results deriving from different laboratory procedures or experimental techniques (e.g., confocal microscopy), thus offering a powerful mean to assist the interpretation of specific biophysical mechanisms.

Materials and Methods

C. elegans strain-culturing, RNAi, and microscopy experiments are reported in the . Furthermore, detail on the analysis algorithms and pipelines and on the statistical analysis and the Monte Carlo modeling adopted in this work are described in the .