A network model for the specification of vulval precursor cells and cell fusion control in Caenorhabditis elegans

The vulva of Caenorhabditis elegans has been long used as an experimental model of cell differentiation and organogenesis. While it is known that the signaling cascades of Wnt, Ras/MAPK, and NOTCH interact to form a molecular network, there is no consensus regarding its precise topology and dynamical properties. We inferred the molecular network, and developed a multivalued synchronous discrete dynamic model to study its behavior. The model reproduces the patterns of activation reported for the following types of cell: vulval precursor, first fate, second fate, second fate with reversed polarity, third fate, and fusion fate. We simulated the fusion of cells, the determination of the first, second, and third fates, as well as the transition from the second to the first fate. We also used the model to simulate all possible single loss- and gain-of-function mutants, as well as some relevant double and triple mutants. Importantly, we associated most of these simulated mutants to multivulva, vulvaless, egg-laying defective, or defective polarity phenotypes. The model shows that it is necessary for RAL-1 to activate NOTCH signaling, since the repression of LIN-45 by RAL-1 would not suffice for a proper second fate determination in an environment lacking DSL ligands. We also found that the model requires the complex formed by LAG-1, LIN-12, and SEL-8 to inhibit the transcription of eff-1 in second fate cells. Our model is the largest reconstruction to date of the molecular network controlling the specification of vulval precursor cells and cell fusion control in C. elegans. According to our model, the process of fate determination in the vulval precursor cells is reversible, at least until either the cells fuse with the ventral hypoderm or divide, and therefore the cell fates must be maintained by the presence of extracellular signals.


INTRODUCTION
Caenorhabditis elegans is a nematode used extensively as a model organism for study in the areas of genomics, cell biology, neuroscience, aging, genetics, developmental biology, and cell differentiation (Hodgkin, 2005;Herman, 2006;Golden and Melov, 2007;Hobert, 2010). In particular, the vulva of C. elegans has been amply used in studies of organ formation, cellular fusion, and intracellular signaling (Sharma-Kishore et al., 1999;Sternberg, 2005;Félix, 2012). The vulva is a small organ with the main functions of copulation and egg laying. Anatomically, it is formed by a stack of seven different epithelial rings, namely (in ventral-todorsal order): vulA, vulB1, vulB2, vulC, vulD, vulE, and vulF, containing a total of 22 nuclei (Figure 1). Each of these rings is either a single tetranucleate syncytium, a binucleate syncytium (vulD) or two half-ring binucleate syncytia (vulB1 and vulB2). Despite its small size, this organ interacts with muscles, nerves, the gonad, and the ventral hypodermis (Lints and Hall, 2009).
Vulval development may be divided into three main stages; namely, formation of the VPCs (Figure 1, 0-12 h), cell specialization, and morphogenesis. During cell specialization (Figure 1, 25-36 h) the fate of the cells is determined by induction from the anchor cell (AC)-a gonadal cell located dorsally with respect to the cell P6.p-, by lateral signaling among the VPCs, and the concentration of Wnt ligands secreted by the AC and other cells near the tail of the worm. Then, during morphogenesis the vulval cells migrate towards the AC, and they fuse forming the seven rings that give the adult vulva its final shape. Furthermore, during this stage the AC breaks the membrane that separates the gonad from the epidermis, connecting both tissues and opening the vulval channel (Sharma-Kishore et al., 1999;Sternberg, 2005;Lints and Hall, 2009).
There are several models describing the process of cell specialization in the vulva of C. elegans (Félix and Barkoulas, 2012). The first developed models were diagrammatic (Sternberg andHorvitz, 1986, 1989), where a concentration gradient of the inductive signal determines the cell fate. Then, dynamical models were created to highlight the importance of the order in the sequence of signals (Fisher et al., 2005(Fisher et al., , 2007, while other models emphasized the importance of the inductive signal gradient (Giurumescu et al., 2006(Giurumescu et al., , 2009Hoyos et al., 2011). Furthermore, some models incorporated an evolutionary perspective (Giurumescu et al., 2009;Hoyos et al., 2011), while other models were developed to test new modeling techniques (Kam et al., 2003(Kam et al., , 2008Sun and Hong, 2007;Li et al., 2009;Fertig et al., high levels of active  In the second fate cells the transcriptional complex CSL-formed by LIN-12, SEL-8, and LAG-1-is active, and the lateral signaling targets LIP-1 and LIN-11 are present. Finally, the third fate cells have the same pattern of activation as a VPC. A comprehensive list of the observed patterns of activation is included as Figure A1, which includes the patterns of expression that have been reported for the VPCs at different stages.

Extracellular signaling
All vulval precursor cells are equivalent before induction. If a VPC is experimentally removed, another one -usually the nearest neighbor-acquires the fate that would correspond to the ablated cell in the wild type. Also, if all VPCs except P3.p are ablated, P3.p acquires the first fate (Sternberg and Horvitz, 1986). These results suggest that the extracellular concentration of three types of ligands, namely WNT, DSL, and EGF, determine the fate of VPCs, while the concentration of these ligands are determined by gradients determined by the relative positions of the VPCs and the AC (Sternberg and Horvitz, 1986;Wang and Sternberg, 1999).

Formation of the vulval precursor cells
During larval phases L1 and L2, before induction, canonical RTK/Ras/MAPK and Wnt signaling maintain the competence of VPCs (Myers and Greenwald, 2007). The presence of the Hox gene lin-39, together with the absence of the Hox genes mab-5 and ceh-13, are necessary for the formation and competence of the VPCs. Importantly, the expression of mab-5 and ceh-13 act as a boundary for the vulval equivalence group. As a result of the activation of Wnt and RTK/Ras/MAPK signaling cascades, the VPCs express LIN-39. This gene, together with its cofactors CEH-20 and UNC-62, activates the expression of ref-2, which inhibits the expression of the fusogen EFF-1. The posterior VPCs P7.p and P8.p express MAB-5, another Hox gene that activates the expression of ref-2. As a result, Wnt signaling mutants have a small penetrant effect on the posterior VPCs (Alper and Kenyon, 2002;Shemer and Podbilewicz, 2002;Shemer et al., 2004;Alper and Podbilewicz, 2008). During L1 and L2, the gradients of the Wnt ligands EGL-20 and CWN-1 are the most important factors controlling cell fusion. Moreover, these ligands may also be necessary for VPC competence, in a LIN- 39-independent mechanism (Pénigault and Félix, 2011b). Finally, their anterior to posterior gradients reach a critical concentration near the VPC P3.p, with the result that P3.p fuse with hyp7 in half of the organisms Pénigault and Félix, 2011a).
GTP-bound LET-60 binds to LIN-45 activating it, LIN-45 then binds to KSR-1 and KSR-2. The LIN-45/KSR-1/KSR-2 complex phosphorylates and activates MEK-2, which in turn phosphorylates MPK-1, which becomes active. MPK-1 then moves to the nucleus, where it phosphorylates target proteins, many of which are transcription factors like  Phosphorylated LIN-1 and LIN-31 are not able to form a complex that binds to the promoter -PJW-5-and inhibits the expression of lin-39. Instead, phosphorylated LIN-1 and LIN-31 activate the expression of LIN-39. Finally, phosphorylated LIN-39 activates its own expression Sundaram, 2006;Wagmaister et al., 2006a), and the transcription of lin-12 .

The canonical Wnt cascade
There are several Wnt ligands, CWN-1 and EGL-20 with penetrant phenotypes, and LIN-44, MOM-2, and CWN-2, with weak phenotypes. Also, there are several members of the Frizzled family of Wnt receptors, of which LIN-17, MIG-1, and MOM-5, are the most important during the vulva formation . A Wnt ligand binds to a Frizzled-family Wnt receptor, and this membrane complex binds MIG-5 and APR-1. APR-1 forms a complex with KIN-19, GSK-3, and PRY-1. This complex marks the β-catenins BAR-1, WRM-1, and SYS-1 for ubiquitination and degradation. Also, when APR-1 is bound to the Frizzled receptor the concentration of BAR-1 increases. BAR-1 forms a complex with POP-1 (TCF), and activates the transcription of lin-39 Wagmaister et al., 2006b).

Determination of the first fate
Twenty-five hours after birth, P6.p responds to the EGF signal LIN-3 activating the canonical RTK-Ras-MAPK cascade, which strongly activates the expression and activity of LIN-39 that in turn activates the transcription of egl-17 (Cui and Han, 2003). Both LIN-39 and egl-17 are markers of the first cell fate. The increased expression of LIN-39 in P6.p during this stage depends on SUR-2 and LIN-25, subunits of the Mediator complex Wagmaister et al., 2006b). The expression of the main components of the lateral signal emitted by the first cell fate, namely the DSL ligands LAG-2, DSL-1, and APX-1 is also SUR-2 dependent in the VPC P6.p. While the ligands LAG-2 and APX-1 are localized on the membrane of P6.p, DSL-1 is secreted . During the determination of the first fate, LIN-39 acts with its cofactor CEH-20 to activate the transcription of elt-5/egl-18 and elt-6. ELT-5 and ELT-6 inhibit eff-1 expression inhibiting the fusion of first fate cells with hyp7 Alper and Podbilewicz, 2008).

Determination of the second fate
NOTCH signaling is a key element for the determination of the second cell fate. Before the larval phase L3, LIN-14 inhibits LIN-12, preventing the second fate determination. Later, lin-4 RNA concentration increases, inhibiting the expression of lin-14 by binding to its mRNA and targeting it for degradation . Twenty-eight hours after eclosion, P5.p and P7.p cells respond to the lateral signal expressed by P6.p, due to the activation of the NOTCH signaling cascade by the DSL ligand . The signal activates , which undergoes a SUP-17 mediated cleavage at the extracellular site 2, and then is cleaved again at the transmembrane site 3 mediated by the γ -secretase protease complex conformed by SEL-12 (HOP-1, the catalytic subunit), APH-1, APH-2, and PEN-2. The resulting intracellular domain of LIN-12 is translocated to the nucleus where it binds to LAG-1 (CSL) and SEL-8, forming a complex that activates the transcription of the target genes ark -1, lip-1, dpy-23, lst-1, lst-2, lst-3, lst-4, and lin-11, among others. Importantly, LIP-1 and LIN-11 are markers of the second cell fate (Greenwald, 2005). Cells that acquire the second fate inhibit the expression of eff-1, very likely mediated by the LAG-1/SEL-8 complex represented by the CSL node since some regulatory regions of eff-1 contain candidate LAG-1/CSL binding sites and NOTCH signaling inhibits the fusogenic function of eff-1 during the formation of the digestive tract of C. elegans (Rasmussen et al., 2008). Finally, the second fate has at least two positive feedback circuits. First, LIN-12 activates the LAG-1/SEL-8 complex which in turn activates lin-12 transcription (Wilkinson et al., 1994). And second, LIN-12 activates mir-61 transcription, which causes VAV-1 down-regulation, and as a result promotes lin-12 activity (Yoo and Greenwald, 2005).
There is another mechanism involved in the determination of the second fate. It has been reported that a diminished LIN-3/EGF signal directly promotes the second fate (Sternberg andHorvitz, 1986, 1989;Katz et al., 1995Katz et al., , 1996. As a result, a model has been proposed  where a small concentration of LIN-3/EGF causes LET-60/Ras to activate RGL-1, which then activates RAL-1/RalGEF instead of Ras activating LIN-45/Raf, and thus inhibiting the first fate and indirectly promoting the second fate. The authors of the model mentioned above also proposed that RAL-1 may directly activate NOTCH signaling.

The first and second fates inhibit each other
The SUR-2 dependent LIN-12 endocytosis and/or degradation is activated in P6.p , thus promoting the first cell fate. In turn, the lateral signal targets inhibit RTK-Ras-MAPK signaling in the VPCs that adopt the second fate. Specifically, LIP-1 inactivates MPK-1 , while ARK-1 inhibits LET-23 in a SEM-5 dependent mechanism . Also, the lateral signal targets lst-1, lst-2, lst-3, lst-4, and dpy-23, inhibit first fate determination in P5.p and P7.p. The loss-of-function mutants of those same targets have phenotypes characterized by the ectopic expression of egl-17 .

The gradients of different WNT ligands determine the polarity of the VPCs
The ligand EGL-20 binds to the receptors CAM-1 and VANG-1 to establish the ground polarity of the VPCs. EGL-20 is secreted by the cells near the posterior end of the worm, so that its concentration is higher in the posterior part of all VPCs . During larval phase L3 the AC secretes the ligands LIN-44 and MOM-2, which bind to the receptors LIN-17 and LIN-18. After these two events, the concentration of polarizing WNT ligands is higher in P5.pp than P5.pa, and higher in P7.pa than in P7.pp.
In the cells with a higher concentration of polarizing WNT ligands, the β-catenins BAR-1, WRM-1, and SYS-1 are not degraded, WRM-1 forms a complex with LIT-1 and POP-1 that moves out of the nucleus , and thus the ratio of SYS-1 to POP-1 is sufficiently large to activate the transcription of some target genes .

THE REGULATORY NETWORK AS A DISCRETE DYNAMICAL SYSTEM
Boolean networks constitute the simplest approach for modeling the dynamics of regulatory networks. These networks consist of a set of nodes, usually representing genes or proteins, each of which may attain one of only two possible states; namely, 0 if the node is inactive and 1 if the node is active. The state of activation of the i th node is represented by x i , and is updated in discrete time steps according to a Boolean function F i such that is the state of the regulators of x i at time t. In the simplest case the functions F 1 , F 2 , . . . , F n are solved simultaneously, which is known as a the synchronous approach. In such case, the dynamics of the Boolean network is deterministic, and for any given initial state the network reaches either a fixed-point or a cyclic state. The set of all these asymptotic behaviors is known as the attractors of such network. Kauffman (1969Kauffman ( , 1993 proposed that the attractors of Boolean networks represent the experimentally observed gene expression profiles that characterize different cell types. We developed a discrete dynamical system of the network that controls fate determination, function, and polarity of the VPCs based on the molecular basis of the regulatory network. We used a generalization of a Kauffman network (Kauffman, 1969), where nodes may attain two or more states of activation (Mendoza, 2006;Sánchez et al., 2008;Schlatter et al., 2009;Franke et al., 2010). Specifically, in our model there are seven nodes with four possible levels of activation -LIN-3 * , LET-23, SEM-5, SOS-1, LET-60, LIN-45, and MEK-2-. The reason is that the VPCs P3.p, P8.p, and P9.p, which acquire the third fate, have no Ras activity (i.e., a level of 0); the VPCs P5.p and P7.p usually have a moderate level of Ras signaling which is sufficient to determinate the second fate (i.e., a level of 1); P6.p is characterized by a high level of Ras signaling (i.e., a level of 2), which is sufficient to determinate the first fate, but only in the absence of negative regulators like GAP-1, ARK-1 and LIP-1; and finally, in some experiments with worms that have more then one ACs, the level of Ras signaling is high enough to overcome the effects of the negative regulators (i.e., a level of 3). For the genes lin-3, lin-23 and lin-60, the phenotypic effects of some mutant alleles correspond to the different levels of activity described above. Moreover, five nodes need three levels of activation. MPK-1, LIN-39, and LIN-39a, which are at the end of the Ras signaling cascade or downstream from it, but they have no inhibitors to overcome, so only the levels 0, 1, and 2 are needed. Then, PJW-5 is a promoter of lin-39, which is completely inactivated by unphosphorylated LIN-1 (level 0), the loss of function of LIN-1 causes a moderate level of activity in PJW5 (level 1), and phosphorylated LIN-1 and LIN-31 completely activates PJW5 (level 2). Finally, EGL-20 * also requires three levels, the highest (level 2) is enough to activate MAB-5, a medium level (level 1) is enough to polarize the VPCs and activate POP-1, and its absence (level 0) is needed to allow VPCs to fuse with hyp7. For the rest of the nodes, the experimental evidence report either a full gain or total loss of function, and therefore only 2 levels of activation are necessary. The rules determining the state of activation of each node as a function of their regulatory inputs, as well as the references used to infer such rules are shown in the Table A1.
We obtained all the attractors of the discrete dynamical system with the use of GINsim (Gonzalez et al., 2006). The reader may find the GINsim implementation of our model as the Supplementary Materials vpcwt23h.ginml. Besides the attractors of the wild type model, we obtained the attractors for all possible single loss-and gain-of-function mutations, included as Supplementary Materials mutants.xls. Furthermore, we also obtained the attractors of the network after systematically removing one interaction at a time; such results are included as Supplementary Materials interactions.xls. Finally, we simulated the processes of fate determination, and the transitions from one cell type to another with the use of a script (Supplementary Materials vpc.py).

THE REGULATORY NETWORK
We were able to reconstruct the regulatory network that controls the VPC fate determination and cell fusion, made of 88 nodes and 126 regulatory interactions (Figure 2). Most interactions were inferred from the experimental data mentioned in the methodology. Four interactions, however, are predictions that are supported by our modeling effort; namely, (1) RAL-1 activates the protein complex that allows the lateral signal targets to be expressed, (2) the self activation of LIN-39 transcription requires LIN-39 to be phosphorylated by MPK-1, (3) the complex formed by LAG-1, LIN-12 and SEL-8 inhibits the transcription of eff-1 in the second fate cells, and (4) the Hox factors MAB-5 and CEH-13 inhibit vulval fate determination by rendering the cofactors UNC-62 and CEH-20 unavailable to  In our model, three or four levels of activation are necessary in some nodes (see Methods) to simulate the effect of additional ACs, with the result of a higher concentration of LIN-3, as well as to describe the reported effects of different mutant alleles.

STATIONARY PATTERNS OF ACTIVATION
The analysis of the dynamical behavior shows that the network has 11 fixed point attractors (Figure 3), all of which can be interpreted as the stable patterns of molecular activation possible for different vulval cells. According to their biological interpretation, these attractors can be grouped into three categories: first fate, second fate, and third fate/VPC.
The first group contains five attractors corresponding to the first cell fate. Attractors in this group are characterized by a high level of activation in the Ras/MAPK signaling pathway, a high level of activation of LIN-39, and the expression of apx- 1, dsl-1, egl-18, elt-6, and egl-17, which is used as a first fate molecular marker. Attractors 1-3 present the highest possible level of extracellular LIN-3. A high level of inductive signaling is enough to determine the first fate, even in the presence of NOTCH signaling, as observed in attractor 1. This particular pattern of activation recovers the expression observed in Sternberg and Horvitz (1989) and Wang and Sternberg (1999), where VPCs acquire the second fate, but they respond to LIN-3 by becoming first fate cells that express egl-17, and whose granddaughters divide transversely. Attractor 4 represents a case where a VPC acquires the first fate with an extracellular LIN-3 concentration as found in P6.p in the wild type, such pattern has already been reported Sternberg and Horvitz (1989). And finally, attractor 5 corresponds to the pattern of expression reported for P6.p after acquiring the first fate. The determination of the first fate in P6.p occurs in an extracellular environment with a high concentration on LIN-3 and no lateral signaling.
The second group, comprised of attractors 6-10, correspond to the second cell fate. These attractors are characterized by the activity of the LAG-1/SEL-8 complex, represented by the CSL node, as well as the expression of the lateral signal targets lst -1, lst-2, lst-3, lst-4, mir-61, dpy-23, lin-11, and lip-1. Importantly, lin-11 and lip-1 are routinely used as second fate markers. Attractor 6 has both a high level of extracellular LIN-3 and an active lateral signaling, as reported in Sternberg and Horvitz (1989). Attractor 7 presents a high level of the inductive signaling, together with active RGL-1 and RAL-1 but no lateral signaling. This pattern has not been reported in the literature and therefore constitutes a prediction of the model. Attractor 8 fits the pattern of activation and the extracellular conditions for P5.p and P7.p, which acquire the second fate in the wild type. These conditions, include a small extracellular concentration of LIN-3 and active lateral signaling. Attractor 9 has a low inductive signaling, active RGL-1 and RAL-1, and no lateral signaling. This kind of pattern has been proposed so as to explain the determination of the second fate in an extracellular environment lacking lateral signaling . Finally, attractor 10 has an active lateral signal, but no inductive signal, as reported in Wang and Sternberg (1999).
The third group contains only one attractor-number 11-which represents the third fate. This attractor is characterized by a low level of LIN-39, and well as the activation of lin -4, lin-12, and ref-2. Importantly, the very same pattern of activation is found in the VPCs after the L3 molt but before induction.

THE DIFFERENTIATION PROCESS
The dynamical modeling of the network is able to describe the capacity of the VPC cells to acquire the first, second, and fusion fates (Figure 4). We simulated the determination of the first fate starting from a VPC ( Figure A2A). While the pattern of a VPC is a stationary state, the sudden activation of LIN-3 to its highest level, which simulates the arrival of a high inductive signal from the AC, originates a cascade of activation that activates the Ras/MAPK signaling cascade, inducing the activation of LIN-39, EGL-17, ELT-6, and EGL-18. Another effect of the activation of Ras signaling pathway is the endocytosis of LIN-12 by the Mediator complex. While the VPC and first state attractors are stationary, the transition of the former to the latter is reversible ( Figure A2B). Indeed, the disappearance of the inductive signal, simulated by turning LIN-3 from its highest to its lowest value of activation, reverses the changes described above. Thus, our model shows that if a first fate cell is moved into an environment with WNT ligands but lacking LIN-3, such cell becomes a VPC; a similar effect has been observed experimentally (Euling and Ambros, 1996;Wang and Sternberg, 1999). The second fate may be acquired by a VPC in an environment containing DSL ligands ( Figure A3A), as described experimentally (Sternberg and Horvitz, 1989;Wang and Sternberg, 1999;Chen and Greenwald, 2004). The presence of the lateral signal activates the NOTCH signaling, which leads to the activation of the lateral signal targets ). In our model, a VPC also acquires the second fate in an environment with a moderate concentration of LIN-3, which is known to occur (Katz et al., 1995(Katz et al., , 1996;   . We propose that this effect may happen through the direct activation of the CSL transcriptional complex by RAL-1 ( Figure A3B). Also, in an environment with both DSL ligands and a moderate concentration of LIN-3 the cells P5.p and P6.p acquire the second fate (Sternberg and Horvitz, 1989), a behavior that is also captured by our model (Figure A3C). In this case the lateral signal targets are activated by the lateral signal before the activation of RAL-1. Finally, as with the first fate, the second fate also shows reversibility of differentiation ( Figure A3D). When a second fate cell is moved into an environment with WNT ligands but lacking LIN-3 and lateral signal, the cell would de-differentiate and become a VPC, reflecting an experimentally observed behavior (Euling and Ambros, 1996).

Hoyos
Our model predicts that a first fate cell in an environment with a lateral signal but no inductive signal would acquire the second fate ( Figure A4A). Conversely, when a second fate cell is moved into an environment with DSL ligands and a very high concentration of LIN-3 ( Figure A4B). The Ras/MAPK signaling activates the first fate markers, while the lateral signal targets remain active, reflecting an experimentally observed behavior (Sternberg and Horvitz, 1989;Wang and Sternberg, 1999).
Finally, it is known that in half of all worms P3.p fuses with hyp7 at about twelve hours after birth. This is determined by the concentration of EGL-20 and CWN-1, which are the only Wnt ligands reaching the cell at this stage of development (Myers and Greenwald, 2007;Pénigault and Félix, 2011a,b). In order to simulate the transition to the fusion fate it was necessary to change the parameters of the simulation, so as to reproduce an environment with no Wnt, LIN-3 or lateral signaling. The resulting dynamics is presented in Figure A5, showing that an environment lacking a sufficient concentration of WNT, DSL ligands, and LIN-3, the fusogen EFF-1 is activated and therefore the VPC acquires the fusion fate.

SIMULATION OF MUTANTS
Experimentally, mutations on genes that control the formation of the vulva in C. elegans may cause the following phenotypes: (1) a vulva-less hermaphrodite (Vul), (2) an hermaphrodite with multiple vulvae (Muv), (3) a worm with two incomplete vulvae (Biv), (4) a worm defective in egg laying (Egl), (5) a worm with at least one protrusive vulva (Pvl), and (6) a fertile individual whose eggs hatch inside the worm (Bag).
We simulated the effect of mutations that cause each molecule to stay at a fixed state of activation (Table A2), and it was possible to assign the development of Vul, Muv, and Egl phenotypes (Table 1). Accordingly, we compared the simulated mutants against the reported mutant phenotypes (Table A3). Importantly, the model was able to describe most of the reported mutants, namely, 15 of 19 type Vul, 11 of 17 Muv, 19 of 32 Egl, 8 of 8 Biv, and 24 out of 26 wild types. Only a few mutants where incorrectly classified, specifically 3 mutants as Vul, 1 as Muv, 3 as Egl, and 23 as wild type. Most of the discrepancies with the reported phenotypes are due to three reasons. First, considering RAL-1 activity as sufficient to determinate the second fate causes many NOTCH mutants to lose their effect. Second, some of the mutants have an effect at a later stage of vulval development not included in our model. And third, some mutants have an effect on less than half of the worms affected, and we could not include such effects due to the deterministic nature of our model. Additionally, the model also allows for the appearance of the Biv phenotype. Specifically, if the presence of second fate cells is possible, the Biv phenotypes arises if the basal polarity of the vulval precursor cells is lost (loss of function mutants for the genes egl-20 * , cam-1, vang-1), or the re-polarization by the ACs is affected in mutants with combinations of the loss of function of the genes lin-44, mom-2, lin-17, lin-18, wrm-1 or sys-1.

EFFECT OF REMOVING A SINGLE INTERACTION AT A TIME
Removing one of 35, out of a total of 115 interactions, cause the model to recover an incorrect dynamical behavior. This is shown in Table A4, which incorporates a partial robustness analysis. These 35 interactions comprise those that are part of the Ras signaling cascade, those that are downstream of the Wnt signaling cascade, those that are essential for ref-2 to inhibit cell fusion, and those which activate the second fate markers.

IMPLICATIONS OF THE MODEL ON SOME SELECTED MOLECULAR MECHANISMS
Our model is able to recover the second fate in an extracellular environment without DSL ligands and with a moderate concentration of LIN-3 (Figures 4, A3B). While this process does not fit the classical description outlined in the introduction, it has been observed experimentally nonetheless (Katz et al., 1995(Katz et al., , 1996. One possible explanation is that the low extracellular concentration of LIN-3 activates Ras signaling in such a way that LET-60 activates RGL-1 instead of LIN-45; RGL-1 then activates RAL-1 which may then cause the determination of the second fate
MAB-5 is usually expressed in P7.p and P8.p. Ectopic expression of MAB-5 in all VPCs reduces their sensitivity to the inductive signal, and in mab-5(lf) mutants, the 3 posterior VPCs are all very likely to acquire the first fate (Clandinin et al., 1997). In ceh-13(lf) hermaphrodites several anterior and posterior Pn.p cells remain unfused, and mab-5(gf);ceh-13(lf) mutants are almost wild type, suggesting that mab-5 may substitute for CEH-13 to control the fusion of Pn.p cells . Based on this information, in our model we propose that the presence of either mab-5 or ceh-13 decrease the amount of cofactors available for lin-39 activity. Not including either the MAB-5 to HOXCO or the CEH-13 to HOXCO interactions, would allow a VPC to acquire a vulval fate even in the presence of MAB-5 or CEH-13.
Another interesting property of our model is that the VPCs that acquire the third fate, do not fuse with hyp7 before dividing because the concentration of Wnt ligands required for the proper specification of the first fate, which are secreted by the AC is high enough to prevent cell fusion. After proper first fate specification and the first VPC division, the concentration of Wnt ligands must drop bellow the critical concentration required to prevent cell fusion.

LIMITATIONS OF THE MODEL AND SOME POSSIBLE IMPROVEMENTS
Our model has 88 nodes some of which have up to 4 possible levels of activation, and the state space of the model is very large 4 7 × 3 5 × 2 76 = 3.008194e + 29, exploring each possible initial state would take a very long time. The use of initial conditions that are biologically relevant, and the efficient algorithm used by GINsim allowed us to find all the fixed point attractors for all our different versions of the network. Simulating the loss of an interaction or a mutation produces a new network, thus, we decided to study only the most important paths that lead from one cellular fate to another one, using the known pattern of gene activity in the VPCs as a starting point. We decided to incorporate the regulatory interaction from RAL-1 to CSL because it improves the whole dynamical behavior of the wild type model. While the absence of such interaction would result in the worsening of the wild type behavior, it would improve the dynamics describing the following mutants: aph-1(0), aph-2(0), lag-1(0), lin-12(0), lin-14(1), lin-4(0), pen-2(0), sel-12(0), sel-8(0), and sup-17(0). A possible solution to reconcile this misbehavior would be to postulate that RAL-1 somehow activates the expression of DSL-1. We will explore this possibility in the near future.
There are two general changes that might improve the model. First, the implementation of the network as a stochastic model. And second, the inclusion of events that occur after the first longitudinal division of the VPCs. All these changes will be explored separately.

CONCLUSION
Our model is the largest reconstruction to date of the molecular network controlling the specification of vulval precursor cells and cell fusion control in Caenorhabditis elegans. According to our model, the process of fate determination in the vulval precursor cells is reversible, at least until either the cells fuse with the ventral hypoderm or divide, and therefore the cell fates must be maintained by the presence of extracellular signals, in agreement with a previous hypothesis (Euling and Ambros, 1996). Furthermore, our model predicts that the trans-differentiation from the first to second fate and vice versa are possible ( Figure A4).

REASONS FOR THE DISCREPANCIES BETWEEN OUR MODEL AND THE EXPERIMENTALLY REPORTED PHENOTYPES
In bar-1(0) mutants, P3.p and P4.p almost always adopt the F fate, while P5.p, P6.p, P7.p, and P8.p adopt the F fate less frequently. Our model is deterministic, and sys-1 may replace bar-1 functionally.
The double mutant elt-6(0);egl-18(0) is usually lethal, but when its effect is limited to the vulva, it causes a Vul phenotype on some of the worms, in our model ref-2 and the CSL complex control cell fusion in these mutants, leading to a wild type phenotype.
ksr-1(0);ksr-2(0) double mutants and lip-1(1) mutants have a Vul phenotype, which our model does not reproduce when a very high level of inductive signal is allowed.
Some ceh-20(0) mutant worms have VPCs to be induced (normally, they are not induced), and P5-7.p not to be induced, originating a Muv and Egl phenotype. In our model, ceh-20(0) causes a Vul phenotype since first fate cells do not form.
48% of dpy-23(0);gap-1(0) double mutants are Muv, but our model produces a wild type phenotype because the precise target for dpy-23 is unknown, and our model is deterministic.
lin-1(0) mutants produce a Muv phenotype, but the first fate does not express egl-17 as much as the wild type P6.p does. According to our model, in such mutants no cell acquires the first fate, but the VPCs may acquire the second or third fate (Case 7 in Table A2); this may cause a Vul or a Muv phenotype. To determine which phenotype may arise, our model would need to cover the events after the first VPC division. Our model indicates, however, that the vulva will not be functional.
lin-31(constitutive phosphorylation) causes a Muv phenotype. In our model this mutation is equivalent to allowing pjw5 to have only two levels of activity, namely 1 or 2, which could cause a Vul, or even a wild type phenotype. In order to simulate the Muv phenotype, our model would need to cover the events after the first VPC division.
lin-12 (1);vav-1(0) increases the penetrance of the lin-12(1) Muv phenotype. In our model, both lin-12_2(1) and lin-12_3 (1), cause a Muv phenotype, and due to the deterministic nature of the model, the effect is not enhanced by vav-1(0). pry-1(0) cause a Muv phenotype in a small percent of the mutants due to defective Wnt signaling regulation, but our model predicts a wild type. In order to simulate this effect, our model would need to be stochastic, and include the events after the first division of the VPCs.
Only a small percentage of mig-5(0) worms present a Muv phenotype or have an ectopic fusion of P5.p or P7.p with hyp7, most likely because dsh-1 and dsh-2 may functionally replace mig-5. In our model we did not include dsh-1 and dsh-2 because their function during the development of the vulva is not known. According to our model losing mig-5 function will cause a Vul phenotype.

FIGURE A3 | The determination of the second fate is reversible.
Only the nodes whose activity changes dynamically are included, the time at which they change is highlighted in yellow and the pattern, which corresponds to the second fate is colored in salmon.