Mating of budding yeast cells is a model system for studying cell-cell interactions. Haploid yeast cells secrete mating pheromones that are sensed by the partner which responds by growing a mating projection toward the source. The two projections meet and fuse to form the diploid. Successful mating relies on precise coordination of dynamic extracellular signals, signaling pathways, and cell shape changes in a noisy background. It remains elusive how cells mate accurately and efficiently in a natural multi-cell environment. Here we present the first stochastic model of multiple mating cells whose morphologies are driven by pheromone gradients and intracellular signals. Our novel computational framework encompassed a moving boundary method for modeling both a-cells and α-cells and their cell shape changes, the extracellular diffusion of mating pheromones dynamically coupled with cell polarization, and both external and internal noise. Quantification of mating efficiency was developed and tested for different model parameters. Computer simulations revealed important robustness strategies for mating in the presence of noise. These strategies included the polarized secretion of pheromone, the presence of the α-factor protease Bar1, and the regulation of sensing sensitivity; all were consistent with data in the literature. In addition, we investigated mating discrimination, the ability of an a-cell to distinguish between α-cells either making or not making α-factor, and mating competition, in which multiple a-cells compete to mate with one α-cell. Our simulations were consistent with previous experimental results. Moreover, we performed a combination of simulations and experiments to estimate the diffusion rate of the pheromone a-factor. In summary, we constructed a framework for simulating yeast mating with multiple cells in a noisy environment, and used this framework to reproduce mating behaviors and to identify strategies for robust cell-cell interactions.
Cell-to-cell signaling via diffusible molecules is an important mode of communication between cells in many mammalian systems such as neuron axon guidance [1], immune cell recognition [2], and angiogenesis [3]. These interactions involve sensing an attractant from the partner and responding by moving or growing in the appropriate direction (i. e. chemo-taxis/tropism), while secreting signaling molecules in a reciprocal fashion. This behavior is conserved in eukaryotes from fungi to humans [4,5]. The budding yeast Saccharomyces cerevisiae, undergoes a mating response that has served as a model system for studying cell-to-cell communication [6]. Yeast cells have two haploid mating types, a and α. By sensing the pheromone molecules (α-factor and a-factor), a- and α-cells detect the presence of a mating partner. These secreted peptides form a spatial gradient, bind to the pheromone-specific receptors, and elicit a response that includes cell-cycle arrest, gene expression, and formation of a mating projection (“shmoo”). Ultimately, the mating response results in the fusion of the two cells and nuclei to create an a/α diploid cell (reviewed in [7]). Mathematical modeling has provided a useful tool for studying cell-cell interactions. Previously, moving interface models have been used to investigate deforming the shape of eukaryotic cells. In [8], a 1D continuum model of cell motility in amoeboid cells based on a viscoelastic description of the cytoplasm was developed, and in [9], cells in a 2D domain were treated as a two-phase fluid. The immerse boundary and finite element based approach was developed to model the actin network and cell morphogenesis in [10], an evolving surface finite element method modeled cell motility and chemotaxis in [11], and the boundary tracking Lagrangian framework was used in [12,13]. Other models used agent-based frameworks such as the Potts model, which takes into account detailed chemical networks and moving cells [14]. Level set approaches have also been adopted [15,16] to simulate the cell membrane deformation coupled to chemistry reaction dynamics. Previous studies focused on the relationship between morphogenesis and its underlying biochemical or mechanical machinery. In this work, we extend this concept by including the molecular dynamics within the extracellular space to study multi-cell interactions. In laboratory yeast mating assays, wild-type cells mate with approximately 100% efficiency [17]. Genetic screens have identified mutants that mate at reduced efficiency [18]. One class of mutants prevents mating altogether. In addition, Hartwell and colleagues have modified the basic assay to investigate “three-way” mating between an a-cell that can mate with either an α-cell that makes α-factor or an α-cell that does not [19,20]. In this mating discrimination test, wild-type a-cells mate almost exclusively with α-factor producers. Mutations that affect the sensitivity of the system, such as the deletion of SST2 (a gene which downregulates signaling via the heterotrimeric G-protein) or the deletion of BAR1 (which encodes for an α-factor protease), dramatically reduce both mating efficiency and mating discrimination [20]. The communication between mating cells is mediated by the mating pheromones which bind their cognate G-protein-coupled receptors turning them on. Active receptor catalyzes the conversion of heterotrimeric G-protein into Gα-GTP and free Gβγ. The resulting Gβγ subunit can then recruit Cdc24 to the membrane where it activates Cdc42. Active Cdc42 is a master regulator of the cell polarity response orchestrating the cytoskeleton, exo/endocytosis, and signaling complexes [21,22]. All of these processes involve noise due to Brownian motion, stochasticity in gene expression or other intracellular fluctuations [23–26], which may affect cell assessment of signals and their responses [27]. In particular, the diffusion of ligand into the local neighborhood of the cell and the binding of ligand to receptor are thought to introduce significant stochasticity to gradient-sensing systems [24,28]. Therefore, it is necessary to consider the effects of noise when exploring cell behavior during mating. There has been extensive mathematical modeling of the yeast pheromone response system. The early models were non-spatial and emphasized signaling dynamics [29–31]. More recent modeling efforts have incorporated spatial dynamics, both deterministic [32–34] and stochastic [35–37]. Models have ranged from simple generic formulations to detailed mechanistic descriptions. Finally, we and others have modeled pheromone-induced morphological changes to cell shape [12,38]. In related research, Diener et al. employed a combination of image processing and computational modeling to describe the extracellular α-factor dynamics in a population of mating cells, and how those dynamics were altered by the protease Bar1 [39]. However, missing from the literature is modeling of the yeast mating process itself involving both a- and α-cells. In this paper, the goal was to construct the computational infrastructure for simulating the mating of two or more yeast cells, and then to investigate the factors responsible for robust mating behavior. We want to use our models to understand in greater detail the spatial dynamics that ensure efficient mating, and provide quantitative explanations and predictions on how perturbing these dynamics (e. g. mutants) disrupts the cell-cell interactions during mating. We succeeded in developing numerical methods for simulating yeast mating. Key elements include modeling the shape of the cell described by a moving boundary technique, and the extracellular diffusion dynamics of the pheromone ligand. Using this framework, we explored different model structures and parameters in a systematic fashion using generic models. We were able to simulate the high efficiency of mating among wild-type cells, and their ability to discriminate among partners that synthesized mating pheromone at different levels. Our simulations suggested that two critical factors ensuring robust mating under noisy conditions were the polarized secretion of mating pheromone, and the presence of the Bar1 protease. In addition, we demonstrated that supersensitive mutants disrupted both mating efficiency and discrimination, reproducing experimental data. More generally, this work makes progress toward the goal of a more detailed description of cell-cell interactions.
In this section, we describe the stochastic model for multi-cell mating systems in two-dimensional space. Cell shape is represented by a level set formulation to capture the deforming plasma membrane induced by pheromone signaling. As described in the Introduction, mating occurs when an a-cell and α-cell are in close proximity (Fig 1A). They sense the pheromone gradient generated by the partner and project toward the source. In Fig 1A, cells are labeled with a marker for the polarisome (a: Spa2-GFP or α: Spa2-mCherry), a cellular structure at the tip of the mating projection. From a simulation standpoint, this process can be broken down into a series of steps (Fig 1B) from the secretion and diffusion of pheromones to the resulting growth in the mating projection. Describing the mating process between two cells requires solving diffusion equations for the ligands in the extracellular space, which evolve according to the shifting positions of the pheromone sources. These sources in turn depend on the sensing of the ligand input and the morphological response. Thus, the cell boundary is evolved together with the molecular dynamics associated with the membrane for each cell. Unlike our previous model of a single polarizing cell which solves surface reaction-diffusion equations in Lagrangian coordinates to capture deformation of the cell membrane [12], here we apply the level set method [40], which can track the moving curve front implicitly by solving a Hamilton-Jacobi equation. In this way, it is easier to study the interactions of multiple cells, and it allows a straightforward extension to the case of multicellular interactions by introducing level set functions for each different cell, and inclusion of the surface diffusions for molecules on the cell membrane. With this methodology, we can distinguish between the intracellular and extracellular space, and couple extracellular pheromone diffusion with the intracellular reaction-diffusion dynamics. The numerical scheme is described in the Methods section. For simplicity, the cell is modeled as a two-dimensional (2D) circle with radius of 1 μm; the actual yeast cell is a three-dimensional (3D) sphere with radius 2 μm. The experimental mating assay involves placing the cells on a surface (i. e. paper filter) so that the mating reaction is effectively in two dimensions. The time unit is 100 seconds to approximate within an order of magnitude the growth velocity observed in experiments. In this model, the mating pheromone is denoted by f, which is the external cue of cell polarization. Two membrane-associated species, u1 and u2, initially are uniformly distributed and then undergo polarization upon sensing the pheromone signal. The system forms a two-stage cascade in which the output of the first stage (u1) is the input to the second stage whose output is u2. The species v1 and v2 provide negative feedback (integral feedback) to regulate u1 and u2 (see S1 Text). The cell grows in the direction determined by u2. This model is a generic model of the mating system and abstracts away the mechanistic details of yeast mating. As studied in the previous model for the two-stage yeast cell polarity system on a single cell [41], u1 represents the protein Gβγ, which is the output of the heterotrimeric G-protein system and the input to the Cdc42 system, and u2 represents active Cdc42, which is the master regulator of yeast cell polarization. Finally, the peak of the u2 distribution represents the polarisome which directs new secretion driving mating projection growth. To track morphological changes of multiple cells, we use a level set function, denoted by ϕ (x, t), to distinguish exterior and interior of one cell such that the domain D is decomposed into three regions: Γ = {x: ϕ (x, t) = 0} representing the plasma membrane, Ωin = {x: ϕ (x, t) < 0} corresponding to the intracellular space, and Ωex = {x: ϕ (x, t) > 0} is the extracellular space. The membrane is moving in a given velocity field V (x) which is described in the Methods section and represents growth of the projection. The stochastic dynamics of the diffusing pheromone ligand (fα represents α-factor and fa represents a-factor) are described in (1) and (2): ∂fα∂t=DαΔfα+Sα (x, t) −kαfα+κ1fα∂2W1 (x, t) ∂t∂x, x∈Ωex (1) ∂fa∂t=DaΔfa+Sa (x, t) −kafa+κ1fa∂2W2 (x, t) ∂t∂x, x∈Ωex, (2) where Sα (x, t) and Sa (x, t) denote the sources of pheromones, and they are either constant or localized Gaussian distributions with support on the membrane (see Methods). Each cell contains the membrane-associated species (uj, vj), j = 1,2, whose dynamics are described in Eqs (3–6), and the membrane velocity is described in (7): ∂u1∂t=DsΔsu1+k101+ (β1f˜) −q1+k111+ (γ1u1p1) −h1− (k12+k13v1) u1+κ2u1∂2W3 (x, t) ∂t∂x, x∈Γ (3) ∂v1∂t=k14 (u˜1−k1ss) v1+κ3v1∂2W4 (x, t) ∂t∂x, x∈Γ (4) ∂u2∂t=DsΔsu2+k201+ (β2u1) −q2+k211+ (γ2u2p2) −h2− (k22+k23v2) u2, x∈Γ (5) ∂v2∂t=k24 (u˜2−k2ss) v2, x∈Γ (6) V (x, t) =Vamp⋅u2⋅max (0, 〈n→, d→max〉), x∈Γ (7) f˜=fmaxx∈Sf+0. 1, u˜1=∫su1ds∫sds, u˜2=∫su2ds∫sds, p1=k101+ (β1f˜) −q1, p2=k201+ (β2u1) −q2. Note that Eqs (3–7) are restricted on the plasma membrane Γ, where Δs denotes the surface Laplace-Betrami operator for the lateral surface diffusion. In Eq (3), f˜ is associated with the pheromone factor from the opposite mating type; that is, if the cell is an a-cell, then f˜ is f˜α. In addition, instead of f, we use the normalized distribution f˜; in this definition, a constant is added to make the parameter consistent with [12]. This normalization represents the adjustable dynamic range mechanisms in the system designed to prevent the sensing from saturating (see Discussion). In Eqs (3–6), we ignore the advection terms ∇S⋅ (ujV) and ∇S⋅ (vjV), which describe the increased surface area where the membrane species reside, because this dilution effect is minimized by the integral control feedback in the model (see S1 Text). On the other hand, these terms necessitate a smaller time step because of the curvature appearing during the computation. Therefore, for numerical efficiency, we simulate Eqs (3–6) on the deforming membrane without the advection terms. In Eq (7), V (x, t) denotes the normal component of the growth velocity of the plasma membrane which describes the rate and direction of membrane movement (the tangential component is ignored by assuming it is small). Vamp is a constant specified with respect to the time scale, n→ is the unit outward normal vector, and d→max is the growth direction defined as the unit outward normal vector at the center of the polarisome. The normal velocity V (x, t) is assumed to be proportional to u2, i. e. , active Cdc42. The concentration and location of active Cdc42 determines the position of the polarisome, which directs the secretion [21]. We model the growth direction to be aligned with the normal direction at the polarisome. Multiplicative noise was adopted and each of the noise terms was weighted by a parameter κi, representing external or internal noise sources. The function Wi (x, t) is a random variable such that the white noise term ∂2Wi (x, t) ∂t∂x follows a normal distribution with variance the same as the time step according to the definition of a Wiener process in our simulations. For simplicity, we considered three noise effects in the simulations. One represents the diffusive noise of the extracellular ligands which is described in the pheromone equation. The second is associated with the dynamics of u1 (Gβγ) which represents noisy internal processes such as fluctuations in ligand-receptor binding and receptor activation of G-protein. The final noise effect represents noise in the regulatory feedback loop (v1). Noise introduced in the second stage of the model is ignored to focus on the sensing noise. In addition, we modified the definition of the velocity function in the stochastic model. Since the velocity depends on u2 which is fluctuating, it is necessary to apply filtering to smooth the dependence of the velocity function on u2 (see S1 Text). Finally, we can explore the deterministic dynamics simply by choosing zero for each κi. The default initial conditions for the simulations are two cells (one a-cell and one α-cell) whose centers are separated by 4 μm in which the membrane species ui is uniformly distributed on the cell surface, and thus no polarisome is formed in the beginning and initial pheromone secretion will be isotropic. Unless otherwise stated, no Bar1 (α-factor protease) is present, i. e. , cells are considered bar1Δ; there is no background α-factor source. At this point we note some of the limitations of the model which we expand upon in the Discussion. First the model is a generic representation of the system that lacks mechanistic detail. Second we employed a quasi-steady-state approximation of α-factor spatial dynamics to speed up the simulations. Third there was not rigorous fitting of the parameters to the experimental data but rather a sampling of different regions of parameter space that produced experimentally observed behaviors. In this section, we investigated the impact of noise in the context of exploring one specific parameter in the simulations, the a-factor diffusion constant. In the natural environment, yeast mating is efficient and robust to a variety of perturbations. In this section, we explored how features of the mating process could promote robust and efficient mating; we compared different mating scenarios by modifying the model parameters. A natural extension of two-cell mating simulations is three-cell mating simulations. In three-cell simulations, the set-up can be either two α-cells and one a-cell, or two a-cells and one α-cell. In the former case, if the two α-cells are equidistant from the a-cell, we found that in the absence of noise, the a-cell projected toward the middle in between the two α-cells. Interestingly, if the two α-cells are slightly offset (i. e. the a-cell is located 0. 1 microns below the middle line of two α-cells so that one is closer), then the a-cell still projected toward the middle (Fig 6A). If the a-cell is Bar1+, then the a-cell is able to gradually reorient to the closer mating partner. However adding noise to the simulations, the Bar1+ a-cell projected toward one or the other α-cell in a random fashion whether or not the cells were offset. Although the no-noise case is somewhat artificial, it indicates how Bar1 can improve the ability to detect the gradient direction in this idealized scenario. Alternatively, the simulation can be between one α-cell and two a-cells. If the two a-cells have different genotypes, then there is a competition between the two for the single α-cell. This corresponds to mating competition experiments, a second important type of mating assay [47], in which one mixes two a-cell genotypes with a limiting quantity of α-cells. We tested the importance of Bar1 using mating competition. In mating competition simulations between Bar1+ and bar1Δ cells we found that the Bar1+ cells mated with the single α-cell partner 20/20 times (Fig 6B). Snapshots of more simulations are provided in Fig K in S1 Text. Greater insight on why the Bar1+ cell has the advantage can be provided by the α-factor profiles for two sample simulations (Fig 6B, lower). The Bar1 helps to remove the excess α-factor so that the Bar1+ a-cell is able to sense the gradient from the α-cell. The bar1Δ cell is stuck in a region of high α-factor in which the gradient is shallower. The third mating arrangement is having a single a-cell choose between two α-cells of different genotypes. One specific scenario is having one α-cell make α-factor whereas the other α-cell makes less or no α-factor. Experimentally this simulation corresponds to a third important type of mating assay: mating discrimination in which the a-cell must discriminate between the α-cell mating partner secreting α-factor from α-cell decoys that do not [20,47]. This assay measures the ability of an a-cell to sense and respond accurately to a pheromone gradient. We imaged mating mixes using both Bar1+ and bar1Δ cells as well as a combination of the two. We found that mating was short-range when the a-cells were Bar1+, i. e. , both a- and α-cells made short projections (see S1 Text). With the bar1Δ a-cells, there was longer-range mating with only the a-cells forming longer projections. We hypothesized that degradation of α-factor by Bar1 resulting in short-range mating in the Bar1+ matings. The projection length in both simulations and experiments was defined by subtracting the initial cell radius from the distance between the center of the cell and the point that is farthest from the center on the cell membrane. The asymmetry in projections lengths in the bar1Δ matings was reminiscent of our simulations in which we varied the a-factor diffusion rate (Fig 2). In particular, as the a-factor diffusion rates became slower, the α-cell projection became shorter (and the a-cell projection became longer). We attributed this difference to the reduced spread of a-factor from its source when its diffusion constant is lower. To provide an estimate of the a-factor diffusion coefficient, we determined the relative length of the α-cell projection normalized by the total distance traveled by both projections, and plotted this α-cell length for both simulations and experiments in Fig 9. In the simulations we varied the a-factor diffusion rate from 0. 1 to 100. From this comparison we estimate that the a-factor diffusion rate is 1 μm2/s.
In this paper we performed computer simulations of the yeast mating process for the first time. The main advance was constructing a computational framework for yeast mating which we used to explore different model structures and parameters. We reproduced qualitatively the basic mating behaviors and calculated the simulated mating efficiency. In addition, we were able to model mating competition and mating discrimination which together with mating efficiency form the three basic assays of yeast mating [19,47]. From a computational perspective, we combined modeling the shape of the cell using a moving boundary technique with the extracellular diffusion of the pheromone ligands with a previously described minimal model of pheromone-induced cell polarity. The simulations were CPU intensive because of the multiple time-scales, the evolution of the level set function over the computational domain, and the calculation of the velocity field. Overall the simulation time depended on the number of cells, time step size, length of simulation, and α-factor diffusion rate. We examined for the first time the coupling among ligand secretion, ligand diffusion, and ligand-induced receptor activation which revealed new cell-cell interaction dynamics that could not be captured in single-cell simulations. We identified key factors that contributed to the efficiency and robustness of mating. First polarized secretion of mating pheromone resulted in higher mating efficiency than isotropic secretion. This finding is consistent with experimental data in which a-factor secretion through the Ste6 transporter is highly polarized [48]. It is likely that α-factor is secreted in a polarized fashion given the polarization of the secretory pathway during mating [22]. A second critical factor is the proper modulation of the sensitivity of the system. In experimental matings, strains that are “supersensitive” show considerably reduced mating efficiency and mating discrimination because they are unable to determine the pheromone gradient direction. By increasing the value of the parameter β1 in Eq (3), we were able to mimic the supersensitive phenotype, and the resulting mating simulations were defective. Finally, the presence of Bar1 helped cells to mate in the presence of background α-factor. Bar1 has been implicated to play an important role in modulating the pheromone dynamics [38,39]. Our results are consistent with the conclusions in [38] that Bar1 helps to shape the α-factor gradient for optimal mating. More specifically, both results show that Bar1 can create an α-factor sink that amplifies the α-factor gradient promoting gradient-sensing. The simulations in this work incorporated stochastic effects, a generic description of intracellular signaling that drives the cell membrane, and polarized secretion of both pheromone and Bar1. There are important limitations to this study. First, we did not attempt to present a detailed quantitative portrait of the mating process with mechanistic reactions. We employed a generic model of yeast cell polarity with a small number of variables for computational efficiency and to facilitate parameter exploration. Second, we employed a quasi-steady-state approximation of α-factor spatial dynamics, although we provide simulation data that this choice does not affect the basic results. Third, we employed mechanisms that only crudely approximate physical reality such as ligand normalization. Fourth, we have not attempted to fit the parameters to actual mating data; rather our approach was to test multiple parameters values to qualitatively explore different scenarios. We thus achieved our goal of constructing a computational framework that is capable of generating realistic-looking responses and reproducing basic behaviors. From a technical standpoint, one important future challenge is speeding up the simulations so that the boundary velocity can be reduced to a more realistic value. One possibility is to employ a quasi-steady-state approximation for the fast α-factor dynamics. For the model with multiple cells, each cell would be assigned with a level set function and a velocity field in our framework, and so there is the potential to improve the efficiency by performing parallel computation for different level set functions or representing all cells by one level set function with a mixed velocity field. Current simulations are all restricted to two-dimensional space. Theoretically it is feasible to extend this framework into three dimensions, although the computation could be very expensive because the computational cost increases exponentially with respect to dimensionality. In addition, experimentally the mating reaction occurs on a surface (i. e paper filter) which is effectively two-dimensional [17]. Importantly this research helps to identify the key processes to focus on for future work. The generic framework is easily extended, and we can incorporate more sophisticated and detailed mechanistic models. Because of the absence of mechanistic details, the models in this work can be thought of as “general mating models”, providing a generic description of gradient tracking informed by the yeast mating system. For example, we plan to replace the normalized f term with pheromone-induced Bar1. In the future, an important goal is to replace the generic terms with more mechanistic terms. With a more realistic mechanistic model of pheromone-induced cell polarity, we could attempt to simulate the mating defects of a variety of mutants. Numerous mutants have been isolated that affect mating efficiency and discrimination including fus1Δ, spa2Δ, etc. [49]. One goal would be to reproduce these mating phenotypes at a quantitative level; another goal would be to predict novel mutants that may affect mating.
An overview of the model including model equations is presented in the main text in Section 1. Here we present additional details. The evolution of the level set function Φ is governed by a Hamilton-Jacobi equation ϕt (x, t) +V|∇ϕ (x, t) |=0, x∈D, in which the velocity field V is defined in Eq (6). More information on how the boundary conditions are imposed on the computational grid as well as other technical details can be found in the S1 Text. The time step is set to be 4 × 10−4 for extracellular pheromone, and 0. 01 for membrane-associated dynamics. The simulations were performed with the authors’ original MATLAB codes, and they can be provided by the authors upon request.
One of the riddles of Nature is how cells interact with one another to create complex cellular networks such as the neural networks in the brain. Forming precise connections between irregularly shaped cells is a challenge for biology. We developed computational methods for simulating these complex cell-cell interactions. We applied these methods to investigate yeast mating in which two yeast cells grow projections that meet and fuse guided by pheromone attractants. The simulations described molecules both inside and outside of the cell, and represented the continually changing shapes of the cells. We found that positioning the secretion and sensing of pheromones at the same location on the cell surface was important. Other key factors for robust mating included secreting a protein that removed excess pheromone from outside of the cell so that the signal would not be too strong. An important advance was being able to simulate as many as five cells in complex mating arrangements. Taken together we used our novel computational methods to describe in greater detail the yeast mating process, and more generally, interactions among cells changing their shapes in response to their neighbors.
cell physiology
medicine and health sciences
enzymes
enzymology
cell polarity
physiological processes
model organisms
sex pheromones
cellular structures and organelles
research and analysis methods
proteins
cell membranes
biochemistry
biochemical simulations
cell biology
physiology
secretion
biology and life sciences
pheromones
yeast and fungal models
proteases
computational biology
Modelling of Yeast Mating Reveals Robustness Strategies for Cell-Cell Interactions