Robust polarity establishment occurs via an endocytosis-based cortical corralling mechanism

In neuronal synapses, zones of endo- and exocytosis are segregated laterally on the cortex such that a central exocytic zone is bracketed by compensatory endocytosis (Haucke et al., 2011). In tip-directed growing pollen tubes, zones of endo- and exocytosis are also spatially segregated; however, endocytosis occurs at the apex, whereas exocytosis occurs in an adjacent zone (Zonia and Munnik, 2009). In budding yeast, endo- and exocytic vesicles localize predominantly to growth sites in the bud (Kilmartin and Adams, 1984; Pruyne et al., 1998). Given the antagonistic relationship of these processes, this arrangement could be incompatible with the membrane efflux required for polarized growth. To address this issue, we investigated the dynamics of these trafficking compartments using high resolution in vivo imaging and in silico modeling.

The Rho GTPase Cdc42 is a key regulator of polarity in eukaryotes. In budding yeast, Cdc42 activates formins at growth sites, promoting the nucleation of actin cables along which exocytic vesicles are transported (Evangelista et al., 1997; Schott et al., 1999). A Cdc42-GTP–scaffold–guanine nucleotide exchange factor (GEF) complex triggers the autocatalytic activation of Cdc42, promoting the recruitment of more exocytic vesicles and polarity establishment (Butty et al., 2002). Endo- and exocytosis have been suggested to maintain polarity by counteracting diffusion via the recycling of polarity factors such as Cdc42 (Marco et al., 2007; Slaughter et al., 2009). Recent studies also indicate a negative role for endocytosis in polarity, possibly by removing polarity factors from the pole (Layton et al., 2011; Savage et al., 2012).

However, polarization of vesicle trafficking and Cdc42 are likely to be coordinated; here, we identify the salient properties of this intricate system in silico using a new stochastic mathematical model, and we study how the spatial organization of endo- and exocytosis arise dynamically in vivo during polarity establishment.

We modeled stochastic biochemical reactions between polarity proteins, their diffusion in the cytosol and on the plasma membrane of a spherical in silico cell, together with vesicle trafficking events using time-dependent simulations (see Materials and methods [M&M]). The model contains three modules, interconnected in a polarity factor delivery and recycling circuit: a previously introduced Cdc42 autoamplification module and endo- and exocytosis modules (Fig. 1). Polarization of the Cdc42 module was previously attributed to a Turing-type mechanism. However, when a Cdc42 diffusion rate of 0.036 µm²s⁻¹ is incorporated, as has been measured in vivo, Turing-based models predict significant loss of polarity (Goryachev and Pokhilko, 2008; Howell et al., 2009; Savage et al., 2012). Without excluding the possibility that a Turing-based mechanism may contribute to polarity establishment, we present here an alternative regimen of the same autoamplification module where polarity emerges from a multistable reaction network (M&M; Semplice et al., 2012). The autoamplification module stimulates exocytosis in two ways. First, active Cdc42-GTP promotes actin polymerization via formin activation, enabling vesicle delivery. Second, the Cdc42 module recruits the exocyst complex, facilitating tethering to the plasma membrane and ensuing SNARE-mediated fusion (Zhang et al., 2001; Zajac et al., 2005; France et al., 2006). Exocytic vesicles are assumed to contain Cdc42, establishing a trafficking-based positive feedback loop, and proteins such as clathrin, seeding the early steps of endocytic vesicle formation (Wedlich-Soldner et al., 2003; Boyd et al., 2004; Newpher et al., 2005). Together, endo- and exocytosis contribute to Cdc42 recycling (Valdez-Taubas and Pelham, 2003; Marco et al., 2007). In an initial simplified model, details of the endocytic pathway were neglected; each endocytic event was instantly triggered by endocytic proteins diffusing on the membrane. In 90% of wild-type (WT) in silico cells, Cdc42-GTP polarized to a region on the cortex during the simulation, whereas endo- and exocytic events clustered in and around the Cdc42-GTP pole. The simulated fluorescence of Cdc42-GTP, exocytic, and endocytic events was plotted as a kymograph by drawing a line around the cell cortex (y axis) projected over time (x axis; Fig. S1 A, i).

Disrupting either Cdc42-GTP autoamplification (Fig. S1 A, ii) or both actin-dependent and -independent Cdc42 recycling mechanisms had drastic effects on polarity establishment (Fig. S1, A [iii] and C). Moreover, polarity was dependent on the localization of endocytosis relative to the exocytic pole where early endocytic proteins were delivered. Colocalization of endo- and exocytosis within the Cdc42 cluster resulted in severe polarity defects (Fig. S1 B, i and iii). In our model, restoring the spatial segregation of endo- and exocytic events rescued polarity to almost WT levels (Fig. S1, B [ii and iv] and C). However, when incorporated into the model proposed by Layton et al. (2011), the spatial segregation of endo- and exocytosis was not sufficient to generate a polarized state (unpublished data).

To explore the role of endocytosis in polarity establishment, we modeled a more detailed endocytic pathway, encompassing distinct phases of coat recruitment to early endocytic patches, actin-based vesicle formation, and subsequent scission (Kaksonen et al., 2003, 2005; Weinberg and Drubin, 2012). In addition to coat proteins, the cargos of this endocytic pathway include Cdc42 and v-SNARES (Valdez-Taubas and Pelham, 2003; Marco et al., 2007). The complete model shown schematically in Fig. 1 (parameters listed in Table 1) is used in the remainder of this paper.

As in the simplified model (Fig. S1), Cdc42-GTP polarized within 20 min of starting the simulation (Fig. 2 A, top; and Video 1). Meanwhile, exocytosis focused to a vertex, around which endocytic events clustered (Fig. 2 A, bottom; and Video 2). These results were robust to individual parameter changes (Fig. S3).

To compare in silico results with in vivo data, endo- and exocytic trafficking compartments were monitored simultaneously in live cells by near-total internal reflection fluorescence microscopy (TIRFM) that facilitated long-term imaging every second for 20 min or more. Many vesicular events were illuminated in the evanescent field by near-TIRFM because post-Golgi secretory vesicles traverse actin cables under the plasma membrane, whereas endocytic sites develop at the plasma membrane (Yu et al., 2011). Endo- and exocytic vesicles were monitored with actin binding protein Abp1-RFP and the Rab GTPase marker GFP-Sec4, respectively.

Trafficking compartments were spread over the cortex in most unbudded WT cells (Fig. 2 B, left nonpolarized cell). However, in a minority of cells, they were confined to a polarized cortical region (Fig. 2 B, right polarized cell), as seen previously (Kilmartin and Adams, 1984; Layton et al., 2011). Kymographs of the endo- and exocytic domains were generated as for in silico data, providing a view of membrane trafficking over a narrow region of the entire cortex (Fig. 2 B). During polarity establishment, and in agreement with our in silico model, an endocytic ring bracketed the constricting exocytic pole (Fig. 2 B and Video 3; McCusker et al., 2012). The diameter of the endo- and exocytic regions diminished abruptly, both in vivo and in silico (Fig. 2 C).

The spatial reorganization of endocytosis correlated with a temporal change in the endocytic frequency during polarization (Fig. 3 A). Endocytic dynamics were analyzed along the time axis of the endocytosis kymograph near the pole, where consecutive endocytic events were observed (M&M). The endocytic signal was differentiated and smoothed to discriminate individual endocytic events from noise (Fig. S2 A and M&M). In a polarized WT cell, endocytic events within the ring became more regular, displaying a constant frequency and amplitude, in contrast to nonpolarized cells (Fig. 3, B and C). The constant amplitude of the differentiated Abp1-RFP signal intensity in polarized WT cells reflected constant quanta of Abp1 molecules being internalized from the endocytic zone (Fig. 3 B). The time interval between consecutive endocytic events at the same position on the cortex was longer in nonpolarized cells (79 s; SD = 50 s) than in polarized cells (53 s; SD = 28 s). The lower SD of time intervals in polarized cells reflected endocytic events becoming regular, or periodic. The variance of the distribution of time intervals was significantly different between nonpolarized and polarized cells (Fig. 3 C). These results indicate that during tightening of the exocytic pole, the endocytic system reaches an equilibrium at which a stable, periodic endocytic flux operates. This behavior, also borne out in silico (Fig. 3 C), appears to be a characteristic signature of the endocytic system in polarized cells.

The observation that regular endocytic events increased in frequency, encircling the exocytic zone while both domains abruptly constrict, led us to hypothesize that endocytosis may confine or corral the exocytic zone. In this scenario, endocytic corralling could limit the spreading of polarity proteins outside the exocytic zone, stabilizing the polarity axis. Endocytic mutants in which corralling is perturbed might therefore display polarity defects. To test this, in silico simulations were run in which coating of early endocytic vesicles and subsequent loading with Cdc42/v-SNARE were impaired. Under these conditions, endocytic dynamics may resemble a sla2Δ mutant, the Hip1R (Huntington Interacting Protein 1) homologue in budding yeast (Holtzman et al., 1993). In these simulations, endocytic events were depolarized over the cortex and exocytic clusters were unstable, consistent with endocytic dynamics contributing to a robust polarity axis (Fig. 4 A, top; and Video 4). In sla2Δ mutants in vivo, the endocytic marker Abp1-RFP was depolarized over the cortex and severe polarity defects were observed (Fig. 4 A, bottom; and Video 5; Hervás-Aguilar and Peñalva, 2010). As predicted in silico, kymographs of trafficking compartments in polarized sla2Δ cells in vivo showed multiple exocytic foci that often disintegrated over time, consistent with the model that endocytic corralling stabilizes the exocytic pole (Fig. 4 A, bottom kymograph; and Video 5). The internalization of Abp1-RFP was less efficient in sla2Δ cells, resulting in longer Abp1-RFP residency times at the cortex (22 ± 12 s) compared with WT cells (9.9 ± 3.6 s; Fig. 4 B). Consistently, the patch residency time of RVS167-GFP, an additional endocytic marker, was also significantly increased in sla2Δ cells (39.3 ± 36 s) compared with WT cells (5.5 ± 3.5 s; Fig. 4 B). The characteristic endocytic signature of polarized WT cells was abolished in sla2Δ cells (Fig. 4, C and D).

We next tested whether other endocytic mutants display polarity defects. Endocytic vesicle formation occurs sequentially via the recruitment of coat, actomyosin, and scission modules to the incipient endocytic site (Kaksonen et al., 2003, 2005; Weinberg and Drubin, 2012). We therefore altered in silico endocytic vesicle maturation rates to recapitulate the effect of sla2Δ, sla1Δ bbc1Δ, and rvs167Δ rvs161Δ mutations on endocytic dynamics (M&M). In silico sla2Δ, sla1Δ bbc1Δ, and, to a lesser extent, rvs167Δ rvs161Δ mutants disrupted the endocytic signature, resulting in irregular endocytic events (Fig. 4 D). These effects were also evident from in vivo analyses: the time interval between endocytic events was longer and irregular in sla2Δ mutant cells (76 s; SD = 66 s), in contrast to polarized WT cells (53 s; SD = 28 s). Similar defects were observed in sla1Δ bbc1Δ (111 s; SD = 83 s) and rvs167Δ rvs161Δ mutants (74 s; SD = 48 s; Fig. 4 D).

To study the effect of modifying endocytic dynamics on exocytic polarization, we analyzed the diameter of the exocytic pole in WT and endocytic mutants. As predicted in silico, sla2Δ, sla1Δ bbc1Δ, but also rvs167Δ rvs161Δ mutants displayed wider exocytic poles in vivo (1.88 ± 1.1, 1.92 ± 1.4, and 1.39 ± 0.82 µm, respectively), in contrast to the focused poles in WT cells (0.75 ± 0.17 µm; Fig. 5 A). Wider exocytic poles were also observed in additional endocytic mutants including ede1Δ (1.21 ± 0.6 µm) and clc1Δ (2.09 ± 1.4 µm). The wider diameter of the exocytic pole in endocytic mutants was also evident from their respective kymographs (Fig. 5 B). These results indicate that a focused polarity axis requires endocytosis-based cortical corralling.

What is the role of endocytosis during polarity establishment? Endocytosis has been proposed to enhance polarization by initiating the recycling of polarity factors (Valdez-Taubas and Pelham, 2003; Marco et al., 2007; Slaughter et al., 2009; Yamamoto et al., 2010; Orlando et al., 2011). In contrast, recent modeling studies have predicted polarized endocytosis to be detrimental to polarity (Layton et al., 2011). To understand how polarized endocytosis contributes to robust polarity, we designed a stochastic mathematical model where trafficking and polarity pathways self-organize. We found that whereas Cdc42-GTP autoamplification drives the clustering of exocytic activity to discrete sites, endocytic corralling ensures the selection of a unique, focused cluster for robust polarity establishment. In agreement with a positive role for endocytosis in polarity, we observed an increased stabilized endocytic frequency within the corralling region as polarization proceeded in vivo, as predicted by in silico simulations. Together with the in silico and in vivo studies in WT and mutant cells, we show the utility of concentrating the cell’s endocytic activity in the corral for maintaining directed exocytosis during robust polarity establishment.

The yeast strains used in this study are listed in Table 2. DDY904 was a gift from D. Drubin and Y. Sun (University of California, Berkeley, Berkeley, CA). The CEN LEU2 GFP-Sec4 plasmid was obtained from R. Collins (Cornell University, Ithaca, NY). PCR-based homologous recombination at the endogenous gene locus was used to tag ORFs at the C terminus. Cells were grown in minimal media supplemented with dextrose at 25°C.

The output from a laser bench consisting of two solid-state lasers emitting at 491 and 561 nm (Cobolt) was fiber coupled to an inverted microscope (Axiovert 200; Carl Zeiss). The microscope was equipped with a 100× TIRFM objective (oil, NA 1.46, plan Apo), a mercury lamp (HBO100; Carl Zeiss), a motorized TIRF slider, a Dualview beam splitter (Photometrics), and a sensitive electron multiplying charge coupled device camera (Evolve; Photometrics). For all measurements presented in the current paper, exocytic vesicles were marked by a GFP fusion to the N terminus of the Rab GTPase Sec4, whereas endocytosis-associated actin patches were marked by Abp1-RFP, unless specified otherwise. The filter set 76 HE (Carl Zeiss) allowed simultaneous excitation of both wavelengths (excitation 406/493/561 nm, dichroic 427/503/578 nm, emission 460/525/608 nm), enabling visualization of exo- and endocytic dynamics. Additional filters were used in the Dualview beam splitter: a dichroic 565dcxr and two single-band pass filters—510/30m and 630/75m for GFP and mRFP channels, respectively. The illumination settings (low excitation power and well separated filter settings) prevented channel bleedthrough. An additional 1.6× Optovar lens provided 100-nm/pixel sampling, compatible with Nyquist-Shannon criterion. The optical system was calibrated daily using beads of 100 nm (Tetraspeck) to align the Dualview channels. A workstation (T7500; Dell) installed with MetaMorph 7.7 software (Molecular Devices) was used for image acquisition, image analysis, and creating movies.

Because of the presence of the cell wall in yeast, it was not possible to image the cells in complete TIRF. Therefore, the angle of incidence was optimized to obtain the best signal to noise ratio to visualize the trafficking events at the plasma membrane. In this case, the angle of incidence was set slightly oblique (5–8°) to the sample plane in near-TIRFM mode, as used previously (Layton et al., 2011; McCusker et al., 2012). The motorized TIRF arm of the microscope was controlled via the MetaMorph software for changing the excitation angle of the laser lines.

For in vivo imaging, cells were grown to log phase in selective minimal yeast media and imaged at room temperature (25°C), whereas fixed cells were imaged in PBS. The low bleaching and high sensitivity of the system enabled long-term acquisition up to 1 h with an exposure time of 100 ms captured every second.

Kymographs were used to analyze the dynamics of endo- and exocytic compartments. A median filter with a filter size of 3 × 3 was applied to the images obtained from the endo- and exocytic domains for signal averaging, before the construction and analysis of kymographs. The kymographs were constructed by selecting a 10-pixel (1-µm)-wide region of interest around the cell cortex for the TIRFM images (Figs. 2 B and 4 A; MetaMorph 7.7).

Intensity line scans along the time axis (x axis) of the endo- and exocytic kymographs were extracted at 50-s time intervals, providing intensity profiles along the cell cortex (y axis) at different time points (MetaMorph 7.7). For each time point, the endo- and exocytic intensity profiles were best fit with a Gaussian function (Origin Pro7; Origin Laboratory). The full width at half maximum (FWHM) was extracted from the fits (as 2.3 times the SD of the Gaussian) and plotted over time, providing the diameter of the exo- and endocytic zones over time during the polarization process (Fig. 2 C).

A 5-pixel (0.5-µm)-wide line scan was created along the time axis (x axis) of the kymograph of the endocytic channel near the pole, where consecutive endocytic events were observed. The line scan on the kymograph was precisely chosen at the position where the intensity of the endocytic marker began to diminish, marking an individual endocytic event (Kaksonen et al., 2003, 2005). Fluctuations of the endocytic intensity profiles were analyzed by differentiating and smoothing the data using the Savitzky-Golay algorithm (Origin Pro7; Savitzky and Golay, 1964). The algorithm used a first order polynomial to calculate the derivatives at each time point and smoothed the data over ten time points. This procedure enabled individual intensity peaks to be distinguished from background, while preserving the features of the original distribution. Endocytic events were sufficiently distant in time compared with the speed of acquisition that differentiation provided changes in intensity corresponding to individual endocytic events. It also enabled analysis of the frequency of endocytic events, as discussed in Calculation of endocytic time intervals.

In the differentiated data (dI/dt), the intensity distribution of each endocytic event is represented by a sine wave encompassing a crest and a trough (Fig. S2 A). The local maxima of dI/dt corresponds to the time point when the slope of the endocytic intensity rises, whereas dI/dt tends to zero as the endocytic intensity peaks. Different controls were performed to study the effect of differentiation on noise and to assess the risk of generating false signals. The background noise was differentiated by selecting a region of interest outside the cell and compared with differentiated data from the cortex. The average intensity collected from the entire cell was also differentiated and compared likewise. The temporal change in intensity was negligible for the data collected from both of these regions in comparison to the data obtained from the cell cortex, confirming the characteristic endocytic signature to be specific for the cortex (Fig. 3 A, black lines compared with red lines in graphs). Also, the time intervals between the peaks of the differentiated data from the cortex matched the intervals between individual endocytic peaks of the original data, confirming that there were no false signals (Fig. S2 A). To confirm that the Abp1 signals corresponded to individual endocytic events, an additional control was performed where the dynamics of Abp1 and the scission protein Rvs167 were compared (Fig. S2 B). Both proteins colocalized before their disappearance, consistent with their cointernalization and the assertion that Abp1-RFP signal intensity changes are a reliable marker of endocytic events.

The differentiated data highlight intensity peaks marking sequences of endocytic events. In contrast to nonpolarized cells, polarized WT cells displayed a constant amplitude of dI/dt, marking internalization of a constant fraction of Abp1 molecules during each endocytic event. To analyze the frequency of these events with constant amplitude, a threshold of 50% of the maximum intensity of each endocytic sequence was set. Only the peaks with amplitude higher than the threshold were selected and the time interval between two consecutive events was calculated and presented as a scatter dot plot. An F-test was used to verify that the variance of distributions between the strains was significantly different (Figs. 3 C and 4 D; GraphPad Prism 5; GraphPad Software). Similar analysis was performed on polarized and nonpolarized WT cells as well as for endocytic mutants. This analysis revealed the emergence of a specific endocytic signature for polarized WT cells with relatively constant frequency and amplitude (Fig. 3, B and C), which was abolished in endocytic mutant cells (Fig. 4, C and D).

The stringent threshold criterion set for the selection of endocytic events limited the number of selected events, especially in polarized WT cells because the time window between the unbudded stably polarized state and bud emergence was short. In endocytic mutants such as sla2Δ and sla1Δ bbc1Δ cells, individual endocytic events were more difficult to discriminate because of the clustering of Abp1 molecules at the cortex (Video 5). The analysis took into account all of the Abp1-containing patches leaving the selected region on the membrane and did not discriminate whether they were internalized completely or retained within the cytosol beyond the evanescent field. These events were still considered endocytic events in the analysis.

To calculate the residency times of Abp1-RFP and/or Rvs167-GFP in the immobile actin patches, images were segmented using the “à-trous” wavelet-based segmentation algorithm to identify individual actin patches (Holschneider et al., 1989; Shensa, 1992; Racine et al., 2007; Fig. 4 B). Fluorescent images were first decomposed into the two first wavelet maps using a B-spline function. Then, actin patches were segmented by applying a threshold of six to eight times the noise SD of the image on the second wavelet map. The residency time was defined as the period between the appearance of the protein and its disappearance, as previously defined (Kaksonen et al., 2003). The residency times were presented as a scatter plot, and an unpaired two-tailed test with Welch’s correction was used to calculate whether the mean residency times of the proteins differed significantly between WT and mutants (Fig. 4 B; GraphPad Prism 5). The residency times of Abp1-RFP and Rvs167-GFP in WT cells were consistent with previous studies, indicating that the imaging conditions and quantitative analyses used here recapitulate the dynamics of these proteins (Kaksonen et al., 2005). In the sla2Δ mutant, clusters containing multiple actin patches were observed to have a cortical residency time exceeding a few minutes. These clusters were not taken into account for the calculation of the patch residency times. Similarly, in the sla2Δ mutant, GFP-Rvs167 containing patches formed stable clusters that were also not included in the calculation.

WT cells and endocytic mutants expressing GFP-Sec4 were fixed and actin was stained with Alexa 543 phalloidin. Z stacks of 3 µm were acquired using a motorized Z Piezo scanner (PI France SAS) with a step size of 0.15 nm. For the analysis of pole diameter, the sagittal plane of the exocytic channel for polarized cells was chosen. A 5-pixel-wide line was selected around the cell cortex for this plane and the exocytic intensity was measured using a line scan. The obtained intensity profile was best fit with a Gaussian function and the FWHM was calculated as the exocytic pole diameter (Origin Pro7). The exocytic pole diameters were scatter plotted and an unpaired two-tailed test with Welch’s correction was used to calculate whether the mean exocytic diameter was significantly different between strains (Fig. 5 A; GraphPad Prism 5).

We modeled a stable Bem1–Cdc24 GEF complex that cycles between the cytosol, where it diffuses very rapidly (cytosolic diffusion constant D_c ≈ 5–10 µm²s⁻¹), and the membrane, where it diffuses with D_m = 0.036 µm²s⁻¹ (Brown and Kholodenko, 1999; Marco et al., 2007; Goryachev and Pokhilko, 2008). The complex binds and dissociates from the membrane with rates k₁ and k₋₁, respectively, and activates Cdc42 with rate k₂. This catalytic activity is increased to k₃ > k₂ upon formation of a Bem1–Cdc24–Cdc42-GTP trimolecular complex (formation rates k₄ and k₇ from membrane-bound and cytosolic Bem1–Cdc24, respectively; dissociation rate k₋₄). Cdc42 is inactivated with rate k₋₂, and inactive Cdc42 binds and dissociates from the membrane with rates k_b and k_u. This approach is equivalent to the Rdi1-mediated circuit used in previous modeling, provided that k_u/k_b = (η + k₋₅/k₅)k₆/k₋₆ and overall binding/dissociation kinetics are similar (Goryachev and Pokhilko, 2008; Howell et al., 2009; Savage et al., 2012). This biochemical network is assumed to follow mass action kinetics, and therefore molecular species concentrations obey the following reaction-diffusion system:

where $Cdc42m*,$ Cdc42_m, and Cdc42_c denote the local concentrations of membrane-bound Cdc42-GTP, Cdc42-GDP, and cytosolic (presumably Rdi1-bound) Cdc42-GDP, respectively, and $GEFm*,$ GEF_m, and GEF_c represent concentrations of membrane-bound trimolecular, bimolecular, and cytosolic GEF complex, respectively. η is the ratio between the volume V_m of the neighborhood of the plasma membrane (η = 0.01; Goryachev and Pokhilko, 2008; Howell et al., 2009; Savage et al., 2012) and the cell volume V_c. Fast diffusion results in a well-mixed cytosol, with uniform molecular concentrations.

Conservation of the total amount of Cdc42 and GEF during the evolution of the system imposes

where N_Cdc42 and N_GEF are the total populations of Cdc42 and the GEF complex, respectively.

Mathematical analysis by Goryachev and Pokhilko (2008) revealed that this autoamplification module exemplifies a Turing-type mechanism of cell polarity establishment. The mechanism relies on a significant difference in the diffusion of membrane-bound versus cytosolic Cdc42 species. However, unlike integral membrane proteins, peripheral membrane-associated Cdc42 diffuses rapidly in WT cells in vivo (D_m = 0.036 µm²s⁻¹; Marco et al., 2007). Moreover, mutants that fail to extract Cdc42 from the membrane ought to display severe polarity defects, as the model predicts that such mutants would disrupt the balance between fast and slow diffusing Cdc42 populations. In contrast, deletion of the sole Rho dissociation inhibitor in yeast (rdi1Δ strains) has subtle effects on polarity (Slaughter et al., 2009). For these reasons, we explored the possibility that other mechanisms may underlie polarity establishment.

Although we do not exclude that Turing-type mechanisms are involved in polarity establishment, we were unable to identify parameters that allow such mechanisms to break symmetry under fast Cdc42 diffusion (D_m = 0.036 µm²s⁻¹) and low Cdc42 copy number (N_Cdc42 = 3,000). However, the analysis of the biochemical reaction network described in the previous section revealed the potential for an alternative mechanism for spontaneous symmetry breaking. This mechanism is based on multistability, the existence of multiple stable uniform steady-states (Semplice et al., 2012). In the (uniform) steady-state, the system of equations for Cdc42 autoamplification reduces to the following:

where GEF_m and $GEFm*$ are functions of $Cdc42m*$ only:

and

We solved this system numerically using Matlab (Mathworks; see Matlab script in online supplemental material). Each equation defines a curve in the plane (⁠$Cdc42m*$ and Cdc42_m), and uniform steady-states are defined by the intersections of these two curves. Under adequate parameter regimes (in silico WT cells; Table 1), the system shows multiple solutions. Successful polarization of the autoamplification module with fast Cdc42-GTP diffusion (Fig. S1, A [i] and B [iv]) supports multistability as a potential mechanism for symmetry breaking and polarity establishment. Remarkably, the autoamplification module approaches multiple uniform steady-states using previously published parameters. However, reaching the multistable regimen using 3,000 Cdc42 molecules (based on the expression levels of other Rho GTPases in the Saccharomyces Genome Database) required subtle adjustments of parameter values (Table 1).

In the presence of Cdc42 fluxes to and from the plasma membrane caused by exo- and endocytosis, respectively, steady-state equations are only quantitatively modified, and the autoamplification module remains multistable within comparable ranges of parameter values.

Cdc42-GTP promotes the nucleation of actin cables, with rate k_nuc. Actin cables grow within d_nuc seconds (Yang and Pon, 2002). To account for the small number of cables observed in individual cells, we assume that the nucleation rate decreases exponentially with the total number of cables n_c (factor e^−βtot*nc). Cables detach with the rate k_det, locally reduced by a factor e^−β*nCdc42 in the presence of n_Cdc42 Cdc42-GTP molecules on a patch of membrane (Marco et al., 2007). Parameters were chosen to yield cable dynamics resembling dynamics observed previously (Yang and Pon, 2002). Actin-related parameters could be varied over wide ranges with little disruption of polarity (Fig. S3, A and D). In the absence of cables, all membrane sites receive exocytic vesicles with the same probability. In the presence of cables, the local exocytic rate is proportional to the local density of cables. The overall secretion rate from an unspecified internal buffer is set to K_ex. Exocytic vesicle cargos include Cdc42, v-SNAREs, exocyst complex subunits, and endocytic coat proteins, such as clathrin (Gall et al., 2002; Gurunathan et al., 2002). An upper limit is chosen for the amount of exocytic cargo on a single individual vesicle (see Table 1). If the internal buffer becomes depleted of a cargo protein, this protein will be absent in the exocytic vesicles until endocytic events refill the internal buffer.

After actin cable-directed secretion, vesicles are tethered to the membrane through the assembly of the exocyst complex. The membrane-bound subunit Sec3 is recruited by Cdc42-GTP (Zajac et al., 2005; France et al., 2006; Zhang et al., 2008). We assumed that the local exocyst assembly rate is the product of the local Sec3 density, proportional to the local density of Cdc42-GTP relative to the total membrane-bound Cdc42-GTP, and a constant k_exocyst. Subsequently, vesicles are docked to the membrane via assembly of the SNARE complex with rate n_vk_snare, where n_v is the number of v-SNARES carried by the vesicle. Exocytosis concludes with the fusion of vesicles with the plasma membrane at the rate k_fus (Liu et al., 2005; Domanska et al., 2009; Karatekin et al., 2010). After fusion, the exocytic cargo is released on the plasma membrane, where v-SNARES diffuse slowly (D_mslow = 0.0025 µm²s⁻¹ [Valdez-Taubas and Pelham, 2003]). Diffusion of peripheral membrane -associated Cdc42 molecules, GEF complexes, and early endocytic proteins is significantly faster D_m = 0.036 µm²s⁻¹ (Marco et al., 2007). In agreement with a multistability-based mechanism, polarity weakly depends on diffusion (Fig. S3 C), although fast diffusion yields larger polar clusters. The probabilities, per second, of these processes are:

where AC(r) is the number of actin cables binding to a site r on the plasma membrane. AV(r), EXO(r), and SNARE(r) are binary variables that equal 1 if there is an exocytic vesicle approaching the site, tethered to the site via the exocyst complex, or fusing to the site via the SNARE complex, respectively. Otherwise, these variables equal 0. N_AC, N_sec3, and n_v are the total number of actin cables at a given time, the total number of Sec3 exocyst subunits (Table 1), and the number of v-SNARES attached to the vesicle being delivered at the location r, respectively. The binary variable χ is 1 if there is already an endocytic vesicle maturing at the same site or an exocytic vesicle being delivered, and 0 otherwise, indicating that vesicle fusion/fission events are mutually exclusive.

In the implementation of the model with the Gillespie algorithm, fluxes of individual Cdc42-GTP, Cdc42-GDP, v-SNARES, and endocytic proteins delivered by exocytic vesicles and internalized in endocytic vesicles are taken into account together with the biochemical reactions describing the autoamplification module.

To analyze the model behavior and inherent mechanisms (Fig. S1), endocytic vesicles are assumed to be instantly internalized with a probability proportional to the local density of endocytic proteins relative to their total, cell-wide density. This preliminary endocytic pathway is described by the following set of equations:

where the first equation represents the diffusion of early endocytic proteins of membrane density EP, and the second equation gives the probability, per second, that an endocytic event is triggered at the location r on the plasma membrane. <EP(r)>_cell is the mean endocytic protein density, and the overall internalization rate K_en = 1.67 s⁻¹ was chosen in agreement with previous studies (Kaksonen et al., 2003, 2005; Layton et al., 2011). Each vesicle is assumed to carry the same amount of Cdc42-GTP, Cdc42-GDP, v-SNARES, and endocytic proteins as exocytic vesicles (see Table 1). The binary variable χ is 1 if there is already an endocytic vesicle maturing at the same site, or if an exocytic vesicle is being delivered, and 0 otherwise.

As modeling the clustering of early endocytic patches is beyond the scope of this article, they were modeled as individual entities that diffuse on the plasma membrane until the coat is recruited. Coat recruitment is triggered at the rate k_coat and completes within the duration d_coat. Then, actomyosin module recruitment is triggered at the rate k_growth. Concomitantly, Cdc42-GDP, Cdc42-GTP, and v-SNARES are loaded into the forming endocytic vesicle at the rates k_load, $kload*,$ and k_loadv, respectively. This loading and the endocytic tubule extension phase completes within the duration d_growth (Liu et al., 2009; Carroll et al., 2012). Finally, the vesicle is internalized at the rate k_sc, and its cargo is delivered to the internal buffer for subsequent exocytosis. The scission rate can be varied over orders of magnitude without significantly affecting polarity (Figs. S3 B and 4 D, rvs167Δ rvs161Δ mutants). These processes are described by the following equations:

where the first equation reflects the diffusion of early endocytic patches and the remaining equations give the probabilities of transition through the steps of endocytic vesicle maturation; e.g., patch immobilization (endocytic patch EP becomes a recruiting patch EP_rec), cargo recruitment and vesicle maturation (recruiting patch become a mature patch EP_mat), and scission (mature patch is internalized). Similar to the various steps of exocytosis, EP_rec and EP_mat are binary variables that describe the current state of vesicle maturation at a given time (if the vesicle is recruiting cargo, EP_rec = 1 and EP_mat = 0; if the vesicle is mature, EP_mat = 1 and EP_rec = 0; otherwise both equal 0). During the recruitment phase, cargo is loaded with cargo protein local density-dependent rates $kload*Cdc42m*,$ k_loadCdc42_m, and $kloadv*v-SNAREm,$ where the v-SNARE density obeys the standard diffusion equation

with a lower diffusion constant (see Table 1).

Wherever possible, parameter values were chosen in agreement with available published data (see Table 1 and references therein). Other parameters (e.g., k_coat, k_growth, and K_ex) were chosen to yield overall endo- and exocytic rates and dynamics in accordance with previous studies (Kaksonen et al., 2003, 2005; Layton et al., 2011). The maximal load of v-SNARE exocytic cargo was set to eight molecules per vesicle (Domanska et al., 2009; Karatekin et al., 2010). In the absence of similar data for Cdc42 cargo, we used the same value for Cdc42-GTP and Cdc42-GDP. Loading rates have little influence on polarity except in the regimes where Cdc42 has a very high affinity for endocytic vesicles (high values of k_load and $kload*;$ Fig. S3 B) or a very low affinity for exocytic vesicles (Fig. S3 D). We assumed that exocytic vesicles carry sufficient endocytic coat proteins to form a single coat. Duration of coat and actomyosin module recruitment phases were set in agreement with measured Abp1-RFP residency times (Fig. 4 B) and recent observations (Liu et al., 2009; Carroll et al., 2012).

Protein populations were estimated from the Saccharomyces Genome Database. The number of Bem1–Cdc24-GEF complexes was limited to the Cdc24 copy number (see Table 1). Owing to a limited number of early coat proteins (e.g., clathrin), the maximal number of early endocytic patches that may simultaneously develop was estimated by dividing the total number of clathrin molecules in a cell by the typical number of clathrin molecules in a coat (Agrawal et al., 2010). The number of Cdc42 molecules in a cell was set to 3,000, in agreement with amount of other Rho GTPases in yeast.

Sla2 is thought to enable binding of the actin cytoskeleton to clathrin coats (Weinberg and Drubin, 2012). sla2Δ mutants were therefore modeled using decreased k_coat and k_growth. Sla1 is an early endocytic protein that contributes to the formation of endocytic patches. Deletion of SLA1 leads to a delay in patch formation and maturation (Kaksonen et al., 2005). We have therefore modeled the sla1Δ mutant using a decreased k_coat. Bbc1 is an inhibitor of Las17, a protein that promotes actin branching via activation of the Arp2/3 complex. Deletion of BBC1 leads to fast growth of actin structures that drive endocytic tubule invagination, uncoupling the invagination and scission machineries. This loss of efficiency in vesicle scission in bbc1Δ strains is further amplified in sla1Δ bbc1Δ and rvs167Δ rvs161Δ double mutants (Kaksonen et al., 2005). We have modeled the bbc1Δ mutant using a decreased k_sc. Rvs167 and Rvs161 proteins are involved in the scission of endocytic vesicles. rvs167Δ rvs161Δ mutants were modeled using a decreased k_sc. Modified rates for mutants are listed in Table 1.

The stochastic evolution of the Cdc42 autoamplification module was modeled using the standard Gillespie algorithm. The algorithm translates the deterministic reaction-diffusion equations of the GTPase module into a probabilistic language. Individual reaction rates were calculated by combining the rates defined in the previous section with reaction stoichiometry, following the standard procedure (Gillespie, 1976). Diffusion of reacting molecular species on the membrane was accounted for in this framework using a method introduced previously (Bernstein, 2005). The plasma membrane was triangulated into elementary patches, and each molecule was labeled according to the patch it belonged to. Then, diffusion (i.e., molecular displacement from one patch to a neighbor) was modeled as a reaction modifying the patch label. Diffusion in the cytosol was supposed to be instantaneous. Endo- and exocytic vesicles at different stages of their evolution and actin cables were described as different reacting species and progression through vesicle events as reactions between these species, with reaction rates given by the probabilities. We implemented long-term reactions, such as the formation of an actin cable, endocytic coat recruitment, or actomyosin-based vesicle growth using a delayed version of the Gillespie algorithm (Bratsun et al., 2005; Barrio et al., 2006).

We implemented this novel algorithm within the programming environment of the free software Surface Evolver and ran numerical computations on a 6-µm-diameter sphere sampled into 320 (simplified model, Fig. S1) or 1,280 (complete model, Figs. 2–5) elementary triangles (Brakke, 1992). The triangle size was chosen to be comparable with the typical size of the membrane region involved in each individual event. 100 Cdc42-GTP molecules, all early endocytic protein patches, and half of the v-SNARES were initially randomly seeded on the membrane. The remaining molecules were shared between the cytosol and the internal buffer. In Figs. 2 A and 4 A, the model was evolved from this initial state for 2 min in the absence of autoamplification (e.g., k₃ = 0) to simulate the nonpolarized state transitioning into the polarized state. In other figures, autoamplification was included from the beginning of the simulations.

We used the Surface Evolver binary image export tool to produce the images of the simulated cell (Figs. 2 A and 4 A) and to produce Videos 1, 2, and 4 at one frame per second. Movies were converted to the .avi format using screen capture with the software VirtualDub, and then to the .mov format using h.264 compression in QuickTime. The kymographs were obtained by creating a contour of the cell along the equator crossing the pole with the highest Cdc42-GTP concentration, similar to kymograph construction for in vivo data. The Cdc42-GTP concentration and endocytic/exocytic events along the cell contour (y axis) were plotted as a function of time (x axis) using ImageJ.

The width of endo- and exocytic regions was obtained by taking a line scan along the y axis (cell contour) of the kymographs and calculating the FWHM of endo- and exocytic signals along this scan. However, because of the stochastic and binary nature of vesicular events in the simulations, line scans performed at a single time point were very irregular and variable, making it difficult to define the FWHM. To smooth the exocytic signal, we calculate the mean over 20 s (Fig. 5 A). This time was chosen because it was long enough to smooth out the stochasticity of the exocytic signal, while being sufficiently short to monitor the slower dynamics of polarizing exocytic vesicles. Because of their wider dispersion over the membrane, endocytic events were less frequent than exocytic events in the polar region, necessitating a longer averaging time (63 s) to obtain a smooth endocytic signal (Fig. 2 C). The endocytic pole size did not show a significant dependence on the averaging time (Fig. S2 C). The exocytic pole size as a function of time was determined directly from the kymograph (Fig. 2 C). Time intervals between consecutive endocytic events in the corral were obtained by taking a line scan of kymographs along the x axis (time) in a window originating at the pole and covering 7.5% of the selected plasma membrane contour (half size of the corral according to Fig. 2 C). Endocytic signal along this line showed peaks corresponding to individual endocytic events, separated by variable time intervals that were plotted as distributions in Figs. 3 C and 4 D.

Fig. S1 shows that stable polarity requires both autoamplification and the segregation of endo- and exocytic events. Fig. S2 shows differentiating Abp1-RFP intensity fluctuations to monitor individual endocytic events and a graph showing that endocytic pole size is independent of the signal averaging time. Fig. S3 shows the robustness of the multistability-based polarization mechanism to changes in model parameters. Video 1 shows that during polarity establishment in in silico WT cells, active Cdc42 is focused to a small region on the plasma membrane. Video 2 shows the reorganization of trafficking domains during polarity establishment in in silico WT cells. Video 3 is of a polarizing WT cell in vivo displaying endocytic corralling of the exocytic zone. Video 4 is of a sla2Δ cell in silico showing perturbed endocytic corralling. Video 5 is of a polarizing sla2Δ cell in vivo showing perturbed endocytic corralling and multiple exocytic poles. The Matlab script calculating the concentration of Cdc42-GTP solves the steady-state equations of the autoamplification module. Additional data are available in the JCB DataViewer at http://dx.doi.org/10.1083/jcb.201206081.dv.

We thank D. Drubin, Y. Sun, and R. Collins for reagents and advice. We also thank the McCusker lab and A. Royou for discussion.

Work in the McCusker lab is funded by Seventh Framework Programme (FP7) Marie Curie grant IRG249298, Agence Nationale de la Recherche grant 2010 JCJC 1210 01, Fondation pour la Recherche Medicale (FRM; grant INE20100518678), Agence pour la Recherche sur le Cancer (ARC; grant PDF20120605172), Centre National de la Recherche Scientifique (CNRS), University of Bordeaux Segalen (UBS), and Conseil Régional d’Aquitaine Volet Recherche 20091301015. M. Jose, S. Tollis, and D. Nair are funded by FRM, ARC, UBS, Region Aquitaine, CNRS, and FP7 Marie Curie IEF Fellowships.