A three-way decision landscape

Dynamical Systems Theory allows us to frame the Waddington landscape idea in mathematical terms. The state of a cell at a particular time t will be represented by the position on the landscape, x(t) = (x(t), y(t)). The evolution of x(t) in time, representing the differentiation of a cell, will be related to the gradient or steepness of the landscape represented by a potential function V(x), and will move downhill in this landscape until it reaches a local minimum (an attractor). Each attractor of V(x) corresponds to a stable fate that a cell can adopt. Given an attractor x*, its basin of attraction is comprised of all the points with a trajectory that ends in that given attractor. Finally, the lowest point of the barriers that separate these basins are saddle points of the potential function (Fig 1B).

Now, in order for the trajectory to move from one attractor to another one, the barrier between them needs to disappear. The induction of a new fate is biologically controlled by cues, known as morphogens or signals, that the cell receives and produce a change in its gene expression, pushing the cell towards a new differentiated state. Translating this idea to mathematical terms, the shape of the landscape must be controlled by these signals and its shape will change as these signals change in time. Therefore V(x) also depends on some parameters ω, which will relate to the biological signals, becoming the family of functions V(x, ω) = V_ω(x).

The most important changes in the landscape are the ones in which an attractor is created or destroyed as the parameters change, and these are called bifurcations. These normally happen when an attractor collides with a nearby saddle. An example of this is shown in Fig 1B. From left to right, by changing the parameter ω of the landscape, the basin of attraction corresponding to the blue attractor, A, gradually becomes flatter as the attractor approaches the saddle point. These two points eventually merge together and disappear. This bifurcation would shift the trajectory of a cell starting in the blue attractor towards either the red or green one. The values of the parameters at which bifurcations occur form the bifurcation set of the potential function V. The bifurcation set of a potential defines regions in the parameter set that correspond to different stability configurations in the landscape, or fate map (Fig 1B).

Choice of landscape.

The problem now is to find the simplest possible parameterized family of landscapes that can quantitatively reproduce the substantial set of experimental results summarized in Table 1. A good starting point is to note that some of the landscapes in this family must contain three attractors as, during the process of C. elegans vulval development, a VPC must be able to adopt one of the three fates introduced above: primary, secondary and tertiary. If we demand simplicity by not permitting any repelling points in the landscape, there are essentially only two landscapes with three attractors (Fig 1C1 and 1C2). An example of a more complex landscape with a repeller is shown in Fig 1C3, other, even more complex examples, also exist.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Table of experimental data obtained from the literature [5].

https://doi.org/10.1371/journal.pcbi.1009034.t001

As discussed above, cell state transitions occur because changes in the signals a cell receives destabilize the attractor where the cell state lies allowing it to transition to a new one. In considering these landscapes it is important to keep in mind that since the attractors correspond to biological cell states, there must be a clear correspondence between each attractor and a particular state, and this relationship must be tracked as the attractors move due to signal changes. We label them here A, B and C. The first landscape (Fig 1C1), which is related to Thom’s butterfly catastrophe [23], differs from the second (Fig 1C2) in that, with regard to decision-making, it is more restrained. It is characterized by the fact that, in the parameter region where there are three attractors, the identity of the central state B cannot change. The only transitions allowed are A → B (by bifurcating A with the top saddle), B → A (by bifurcating B with the top saddle), B → C (by bifurcating B with the bottom saddle) and C → B (by bifurcating B with the bottom saddle). This means that cells need to pass through state B in order to go from A to C. This landscape, in fact, corresponds to the morphogen model, in which a VPC would stay in tertiary fate (A) under a low EGF signal, would transition from tertiary fate to secondary fate (A → B) in response to a middle EGF signal, and to primary fate (A → B → C) in response to a high EGF signal. This landscape, however, is not consistent with the AC ablation experiments in [22], where early AC ablation experiments show an equal escape of cells to secondary and primary fates (Table 1) at early AC ablation times, since this means that cells take the same time to go from A to B than from A to C.

On the other hand, for the alternative landscape (Fig 1C2), the fact that the unstable manifold of the top saddle A can flip from B to C via a heteroclinic connection allows extra transitions and for either of B and C to be the central attractor (Fig 1D). Therefore, in this case, the transition A → C is also allowed and cells can go to the state C without passing through B, which is consistent with the AC ablation experiments mentioned earlier. We call the landscape in Fig 1C2 the binary flip with cusp landscape.

The more complex landscape containing a repeller shown in Fig 1C3 is essentially that used in [3, 5], where it was shown that it can reproduce the data we are concerned with. In this model, all attractors are equivalent to each other and all bifurcations are allowed. Biologically, this means that all cell states are equivalent and cells can transition to any fate at all times.

Our aim is to show that the same data set can also be explained by the considerably simpler binary flip with cusp landscape. This results in a different biological interpretation of the decision-making process. In this picture, not all the fates are equivalent, but the tertiary fate is a special fate. It is the default state for a cell if it does not receive any signal; and cells do not return to it once partially induced, even if the signals are switched off (as AC ablation experiments show in Table 1 [22]). Transitions from fate 3° to fate 1°, and fate 3° to fate 2°, are allowed [24]. Also, the transition from fate 2° to fate 1° is achieved by increasing the EGF signal and turning Notch signal off [17]. These reasons lead us to hypothesize that the topology of decisions is represented by that of the binary flip with cusp landscape. First, cells will decide whether to differentiate into a vulval fate or stay non-vulval (determined by a bifurcation between the blue attractor A and the top saddle in Fig 1C2) and, if they do so and surpass the top saddle, then they will decide whether to become 1° or 2° fated cells (determined by bifurcations between the bottom critical points in Fig 1C2).

Building the landscape.

In order to find a representative model for this landscape, we could attempt to use CT to search for a parameterized set of functions such that the parameterized landscape is given by the gradient flows of these functions. However, this approach is not practicable for very complex landscapes. Therefore we employ an alternative method that can be used in this more general context. We use catastrophe and bifurcation theory to analyze the local structure and then use ideas from differential topology to glue the local models together. In our case this involves composing two simple catastrophes: a fold (commonly also known as a saddle-node bifurcation) that will explain the non-vulval vs vulval bifurcation (bifurcation between attractor A and top saddle in Fig 1C2), and a cusp catastrophe that will explain the primary vs secondary fate bifurcation (bifurcation between attractors B and C and bottom saddle in Fig 1C2). For gradient and gradient-like systems these are the only local bifurcations that occur generically and are universal in one and two-parameter systems [25].

The fold or saddle-node is the only bifurcation that appears generically in one-parameter gradient-like systems of any dimension n, and it has the normal form (i = 2, …, n), where (x₁, x₂, …, x_n) is a choice of coordinate system, and are parameters (Fig 2A). Any system displaying a saddle-node bifurcation can be transformed into this [25]. Here, when the parameter c < 0, the system has two equilibria, an attractor and a saddle, and when c > 0 it has none (see S1 Appendix for more details). Thus, it can represent the disappearance of the tertiary fate upon receiving signals. The distance between these two equilibria grows like . We will use it to define the attractors corresponding to 3° fate and the point of transition between this non-vulval state and the other two vulval states (bifurcation between 3° attractor and top saddle in Fig 2C).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Binary flip with cusp model of vulval development in C. elegans.

(A) The fold catastrophe, controlled by the parameter c, defines the stability of the attractor representing 3° fate. The fate map is represented at the bottom, showing that the fold has a bifurcation at c = 0 (in orange), separating two stability regions: one containing 3° fate (for c < 0) and one where the fate disappears (c > 0). At the top, from left to right, different examples of landscapes for different values of c. For c < 0, the potential has two critical points in the state space: a minimum (blue circle) and a maximum (white triangle), such as the one represented on the top panels. Close to the bifurcation point, c = 0, the two critical points get near each other to merge together. For values of c > 0, the potential does not contain any critical point. (B) The cusp catastrophe, controlled by the parameters a, b, defines the stability of the attractors representing 1° and 2° fates. The fate map, at the center, shows that the cusp has a bifurcation set determined by Δ = 8a³ + 27b² = 0 (purple cusp), separating three stability regions: one containing 1° and 2° fate (for Δ < 0) and two where only one of the two fates is present (Δ > 0). Examples of landscapes determined by different values of a, b are shown; in each case, 1° and 2° are represented by red and green circles, respectively, and the saddle between them is represented by a white triangle. On the cusp mid-line (gray dotted line in the fate map) the two fates are equally probable. Getting closer to the cusp favors one fate over the other, by the saddle approaching one of the attractors. Outside the cusp, the landscape only has one attractor, either 1° or 2°. (C) The two catastrophes are combined to generate a 2-dimensional flow with two saddles (white triangles) and three attractors representing each fate (blue, green and red circles represent 3°, 2° and 1° fates, respectively). By mapping the signals θ_E, θ_N onto the control parameters a, b, c, a fate map is obtained, defining the available fates in the landscape for different signaling profiles. The presence of 3° fate is controlled by the fold line, as shown in (A), and the stability of the 1° and 2° fates is controlled by the position of the signaling profile with respect to the cusp line, as shown in (B).

https://doi.org/10.1371/journal.pcbi.1009034.g002

For 2-parameter gradient-like systems the only other generic bifurcation is the cusp. Similarly, the relevant normal form is (i = 2, …, n), where (x₁, x₂, …, x_n) is a choice of coordinate system, and are parameters. In this case, when the discriminant Δ = 8a³ + 27b² < 0, the normal form system contains two minima and a maximum; when Δ = 0, it contains one minimum and a degenerate point; and if Δ > 0, it contains just one attractor (see Fig 2B and S1 Appendix for more details). We will use it to define the attractors corresponding to fates 1° and 2° and the saddle between them (Fig 2B and 2C).

We merge the fold and the cusp catastrophe models into a two-dimensional dynamical system (Fig 2C), the mathematical details of which are in Box 1 and S1 Appendix.

Box 1: Mathematical details of the model

We merge the fold and the cusp catastrophes into the following model, which will be used to model the evolution of the state of a single VPC in time, represented by x(t) = (x(t), y(t)): (1) where H is a Heaviside function that is equal to 0, if y is less than 0, and equal to 1 otherwise. Strictly speaking to get a smooth dynamical system we should smooth the Heaviside function but, as will become clear and is shown in S1 Appendix, this does not affect the analysis shown here.

The attractors of this system are obtained by solving f₁ = f₂ = 0. Since f₂ does not depend on x, one can quickly solve it to obtain the y-coordinates of the critical points, which are given by , and . If f₁ = 0 is now solved with each of the possible values of y obtained, the x-coordinate of the critical points with or is equal to 0; while the x-coordinates of the critical points with are given by the zeros of the cusp, i.e. f_cusp(x*, a, b) = 0. The parameter M controls the relative position of top and bottom saddles in Fig 2C.

The stability of these points can be obtained by looking at the Jacobian matrix of the system. It can be shown that the point is stable (blue attractor in Fig 2C, assigned to represent 3° fate), the point is unstable (white top triangle in Fig 2C), while the points laying on the x-axis can be either two stable points (representing 1° and 2° fates) separated by a saddle (as in Fig 2C), a stable point and a degenerate point or just a stable point. We refer the reader to S1 Appendix for more details and for more examples of landscape configurations.

Taken together, the parameter c controls the existence and positions on the y-axis of critical points with positive y coordinate (which can be two, one or none), while the parameters a, b control the the positions on the x-axis and existence of critical points with zero y-coordinate (which can be three, two or one). With this model, 3° attractor can be bifurcated away and be removed from the landscape, however, either 1° or 2° fates need to be present and both cannot be bifurcated away to give a landscape where only 3° fate is present, which is intrinsically different from the model in [3, 5].

Since the 3° attractor is controlled by the fold we know that if c < 0, the attractor corresponding to this fate will be present, no matter the values of a and b. On the other hand, the presence of the attractors corresponding to 1° and 2° fates depends only on the values of a and b, specifically on the value of the discriminant Δ = 8a³ + 27b² as explained above, independently from the value of the parameter c, and therefore their presence will be controlled by the value of Δ.

One can then compute the bifurcation set of the landscape in Eq 1, and obtain: (2) where and are represented in orange and purple colors, respectively, in Fig 2. Bifurcations will happen at and only at the values of the parameters in . This bifurcation set divides the control space into regions with common stability, some of which are represented in S1 Appendix. More details can be found in S1 Appendix.

The steady states of this system will lie either on the x-axis, in which case their x-coordinates will be given by the steady states of the cusp (i.e. there will be three, two or one steady states on the x-axis); or they will have y-coordinates which are the zeros of a fold, (there will be two, one or none depending on c) in which case their x-coordinates will be 0 (see Box 1, Fig 2C and S1 Appendix for more details).

We choose the attractor with zero x-coordinate to represent the 3° fate, the attractor with zero y-coordinate and negative x-coordinate to represent the 1° fate and the attractor with zero y-coordinate and positive x-coordinate to represent 2°, as shown in Fig 2C. Taken together, parameter c will control whether 3° fate is present in the landscape, and a, b will control whether the attractors representing 1° and 2° fates will be present in the landscape.

The advantage of building the landscape from these two catastrophes is that, firstly, mathematically one has the universality property mentioned above and further described in S1 Appendix and, secondly, that the position and stability of the various restpoints and their dependence on the parameters is transparent and, as shown in Box 1, one has complete understanding of the regions of the parameter space with common stability, being able to characterize the shape of the landscape for any parameter value.

Parameter estimation: ABC SMC implementation

Now that a parameterized model has been developed, the next step is to find whether there are parameters that allow the model to reproduce the experimental data in Table 1 and to determine the extent of such parameters. Table 1 contains the probabilities of each VPC (P4.p, P5.p, P6.p) acquiring each fate (1°, 2°, 3°) under twenty one experimental conditions (i.e. 171 data points or probabilities in total). From this data set, we will use a subset of nine experimental conditions as training data set for the model to fit, and the remaining ten as validation results (TD and VD in Table 1, respectively).

Let us denote by p_e,f,c the experimental probability of cell c becoming fate f under experimental condition e, the values of which are given in Table 1. Similarly, represents the simulated probability of cell c becoming fate f under experimental condition e given parameters θ. Taking a Bayesian approach, we treat the parameters of the model as random variables with a probability distribution. Our goal is to estimate the posterior distribution of the parameters that accurately reproduce the training data set X₀ = (p_e,f,c)_{e ∈ TD}, and check that they also reproduce the validation data set.

If we denote a parameter vector by θ, the posterior distribution satisfies: (8) where π(θ) is the prior distribution of the parameters θ, is the likelihood of the data X₀ given the parameters θ and π(θ∣X₀) is the posterior distribution of the parameters θ given the data X₀.

Since it is not possible to find an analytical expression for the likelihood of the model we propose, we approximate it using approximate Bayesian computation (ABC). ABC methods, also known as likelihood-free methods, have been developed for inferring the posterior distribution of the parameters when the likelihood function is either too complex or too expensive to compute but observations of the model can be simulated fairly easily. They have been successfully applied to a wide range of intractable likelihood problems over the past twenty years [26] including population genetics [27], pathogen transmission [28], reaction networks models [29, 30] and epidemic modeling [31], among others. ABC methods use a comparison between the experimental and simulated data to measure the goodness of fit, instead of using the likelihood function, allowing for more flexibility.

In particular, here we take advantage of Sequential Monte Carlo ABC (ABC SMC) to explore the parameter space and determine the posterior distribution of the parameters. Various ABC SMC algorithms have been proposed in the literature [27, 32–34], but here we focus on the fairly general ABC SMC algorithm in [34]. The ABC SMC sampler methodology approximates the posterior distribution by sequentially sampling from a sequence of intermediate distributions, or approximate posterior distributions, {π_t}_0≤t≤T, that increasingly resemble the posterior distribution given in Eq 8: (9) where d is a metric that compares the experimental and simulated data, and {ε_t}_0≤t≤T is a decreasing sequence of distance values. This ABC method is particularly useful in cases where the number of parameters and the data are high dimensional, as it can be easily parallelised. It can also be tuned to efficiently explore the parameter space maintaining areas of high likelihood via the choice of its perturbation kernel. A summary of how this is done is shown in Fig 3.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Application of ABC SMC to the fitting of the binary flip with cusp model.

(A) First step of the algorithm. Having defined priors for the fitted parameters and constraints between them, as well as an initial threshold ε₁, particles are sampled from the priors and accepted if the distance between the simulated data and the experimental data is less than the initial threshold. This is repeated until N particles are accepted. (B) An initial approximation of the posterior distribution is obtained, generated by the N particles accepted in (A). (C) The algorithm then proceeds in a sequential manner. In each step t of the algorithm, a new threshold ε_t is defined, in our case, the 0.3 quantile of the previous distribution of distances ε_t−1. Particles are sampled from the last approximate posterior distribution and perturbed using a Markov Kernel obtained from the Optimal Local Covariance Matrix (OLCM) [35] (black ellipses in the figure represent confidence intervals), (as described in S1 Appendix). This is repeated until N new particles that satisfy the distance threshold are accepted. Also, each new particle is assigned a weight proportional to its prior probability and inversely proportional to the Markov Kernels evaluated at this particle, which control for efficient exploration. (D) After each step, a new improved approximate posterior distribution is obtained, which restricts the values of the parameters to a restricted region of the parameter set.

https://doi.org/10.1371/journal.pcbi.1009034.g003

Our model depends on a total of 22 parameters: 1 time constant that controls the velocity of the trajectories in the landscape, 10 landscape parameters that map the signaling values into the shape of the landscape and 11 signaling parameters that define the signaling profile that each cell is exposed to in time (see S1 Appendix). Both the model’s definition and the experimental data constrain the range of the landscape parameters and inform the choice of their priors, for which we choose flat priors on their support. For the rest of parameters we choose fairly non-informative priors that still reflect our knowledge about their ranges. Moreover, taking into account the dimensionality of the parameter space, we choose to approximate the posteriors by sampling N = 2 × 10⁴ particles at each step of the algorithm. More details can be found in S1 Appendix.

In order to compare the experimental and simulated data we define a distance d that measures the level of similarity between two data sets. If is the subset of data corresponding to experiments in the training data set (see Table 1), we define the distance between and the corresponding simulation X^TD(θ) as (10) where E is the number of experiments in the training data set, p_e,f,c and are the experimental and simulated probabilities of cell c becoming fate f in experimental condition e, respectively, and is the proportion of times that our model could not assign a fate to cell c when simulating experiment e (see S1 Appendix for more information). The distance function penalizes parameter values for which our model cannot assign a fate, since we would like to avoid these parameter values. More details about our implementation of the fitting algorithm such as the choice of thresholds sequence and perturbation kernel function can be found in S1 Appendix.

Our choice of experiments for training was based on trying to find a maximally informative set while constraining the number of experiments. For example, we included all experimental conditions (a total of 6) with fully penetrant phenotypes resulting from different signaling combinations, which explored different regions of the signal space. With regards to partially penetrant phenotypes, we included one experimental condition involving EGF overexpression. We expected that including these experiments would inform the fitting of the values of the scaling parameters s and l and other parameter values related to the strength of the signals and the crosstalk between them. To constrain the temporal features of the model, we also added two AC ablation conditions which we considered the most informative ones, expecting that the model would be able to reproduce the rest. The particular details of the fitting can be found in S1 Appendix.

A description of the methodology, software implementation in MATLAB, and instructions for use are publicly available in the following GitHub repository: https://github.com/ecamacho90/VulvalDevelopment.

Fitting the model to the training data: A two-step decision

The results after 14 generations of ABC SMC are shown Figs 4 and 5, S1 and S2 Figs, and S1–S13 Videos, which show excellent agreement overall between the simulations and the data. We are able to reproduce both fully and partially penetrant phenotypes.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Fitting of the training data.

(A) The patterns of the 6 fully penetrant phenotype perturbations considered in the training data set are correctly fitted by the binary flip with cusp model. (B) The pattern of the partially penetrant phenotype given by the EGF overexpression perturbation considered in the training data set is also correctly fitted by the model. On the left, the experimental patterns. On the right, the mean simulated patterns and trajectories on the landscape. In both the experimental and simulated patterns, primary, secondary and tertiary fates are represented by the red, green and blue colors, respectively, and proportions higher than 90% have been rounded. On the landscape, the mean simulated trajectories for the particle with best overall approximation for P4.p, P5.p and P6.p are colored in blue, green and red, respectively.

https://doi.org/10.1371/journal.pcbi.1009034.g004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. AC ablation fitting and validation.

(A) Patterns observed when the AC is ablated at different developmental stages. (B) Mean simulated patterns for the different AC ablation conditions. Only proportions for P5.p and P6.p are available, which are presented in blue, green or red representing tertiary, secondary or primary fate, respectively.

https://doi.org/10.1371/journal.pcbi.1009034.g005

S1 Video shows how the simulated VPCs pattern in the WT case. The two-step decision logic of non-vulval vs vulval followed by primary vs secondary fates can be clearly seen. As P6.p starts receiving EGF signal, the cell leaves the attractor corresponding to 3° fate, i.e. the first decision is made, and moves into the regions of the fate map and landscape containing vulval fates. The increase in EGF signal and lack of Notch signal from its neighbors, positions its signaling profile in a region of the fate map where only primary fate is stable, and therefore, its trajectory on the landscape moves towards this attractor. A similar effect happens for P5.p but, due to the higher distance to the AC resulting in lower EGF signal, but higher Notch signal from the neighboring P6.p differentiating to primary fate, the cell differentiates towards 2° fate. Finally, the signals received by P4.p are not enough for it to escape the attractor corresponding to 3° fate, and therefore the cell differentiates into the non-vulval fate. S2–S7 Videos show the differentiation of VPCs in mutants 2–7 in Table 1.

It is worth mentioning that, in our model, due to the lack of information about the exact timings of the stages at which AC ablations were performed, relative to each other, and also about the EGF and Notch signal dynamics, the increase in EGF and Notch signals are modeled by monotonically increasing functions, which are initially hypothesized to be sigmoidal functions of the simulation time. We find that the fitting of the data suggests a slight modification of the EGF dependence on time (see S1 Appendix) which then allows for the model to fit the AC ablation data (Fig 5 and S8–S13 Videos). This is probably due to the fast dynamics around the bifurcation point between the attractor corresponding to tertiary fate and the saddle. More details can be found in S1 Appendix.

Analysis of the evolution of the approximate posterior distributions of the parameters at each step of the ABC algorithm shows that most of the parameters are very well constrained by the data (S1 and S2 Figs), with the exception of the scaling parameters s and l that correlated with the parameter q₃ of the linear transformation that determines the position of the fold line. Also, the parameter l_d, which controls a direct coupling between EGF and Notch pathways, was not crucial to reproduce the data, and this is especially important because, as we will show in the following section, the model is able to reproduce epistatic effects without this coupling. Finally, the exponential decays λ_E and λ_N were only constrained such there was sufficient time to allow the system to return to zero-signal condition (S1 Appendix). Interestingly, the noise magnitude represented by the parameter σ_dif was strongly constrained, suggesting that the data strongly constrained noise-driven fluctuations. Moreover, the fitting did not show strong parameter correlations with the exception of parameters H_E and M_E, defining the sigmoidal increase of EGF signal in time. Also, it is important to note that the parameters γ and n₁ that control the balance between the morphogen and sequential models were strongly constrained. In particular, γ = 0.22±0.06 < 1 and n₁ pointed towards first fate, so we can infer from the data that both morphogen and sequential features were necessary.

Finally, we also observe that the results of fitting suggest that the geometry of the fate map is strongly constrained by the data, the details of which are shown in S1 Appendix. Taken together, our results suggest that aspects from both the morphogen and the sequential model are important for the correct patterning, which combine into a two-step decision logic determined by the fitted fate map and controlled by the dynamics of the EGF and Notch signals received by the cells.

Validating the model: The model reproduces epistasis between EGF and Notch

There is an extensive amount of experimental data available with different combinations of signaling perturbations. In particular, [21] provides quantitative signaling data for many perturbation lines. We compiled a set of mutants from the literature and tested whether our model was able to reproduce this data (Table 1). An important feature of the model developed in [3, 5] is that it can reproduce epistatic effects between the signals, i.e. a mutation in EGF signal can alter the effect of Notch and vice versa. Here we show that our model can also explain these effects without including any direct coupling between the pathways.

As shown in [21], an EGF overexpression perturbation of 2.75-fold with respect to the WT (based on measured lin-3 mRNA levels), named JU1107, increased P5.p induction towards primary fate, and P4.p towards secondary fate. This is an epistatic event, as an increase in EGF signal increases the secondary fate in P4.p, promoted by Notch signaling. Our model correctly reproduces these features as shown in Fig 6A and Table 2. This is achieved because a 2.75-fold increase in the EGF signal of P5.p locates its signaling profile close to the cusp mid-line in the signal space, where primary and secondary fates are equally probable. In turn, primary-fated P5.p cells signal through Notch to their P4.p neighbors which are receiving more EGF than the WT P4.p cells. This positions P4.p closer to the fold line, in a region where secondary and tertiary fates are equally probable (Fig 6A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Comparison between experimental and simulated data for the validated experimental conditions.

https://doi.org/10.1371/journal.pcbi.1009034.t002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Model’s further validations.

(A) Left: Representative mean trajectories on the landscape and signal space for simulated WT vulva and different signal perturbations included in the validation data set. Blue, green and red trajectories represent the evolution of the state of P4.p, P5.p and P6.p, respectively. Right: Mean simulated and experimental patterns for the different perturbations shown on the left. Blue, green and red represent tertiary, secondary and primary fates, respectively. (B) Representative mean trajectories on the landscape and signal space for simulated WT vulva and EGF hypomorph perturbations. Blue, green and red trajectories represent the evolution of the state of P4.p, P5.p and P6.p, respectively. (C) Mean simulated and experimental patterns for the EGF hypomorph perturbations. Blue, green and red represent tertiary, secondary and primary fates, respectively.

https://doi.org/10.1371/journal.pcbi.1009034.g006

Reduced Notch perturbations do not have a strong phenotype unless crossed with an EGF overexpression mutant (see Table 2), in which case the probability of P5.p differentiating into secondary fate is slightly reduced. We checked whether our model was capable of reproducing these features. Since the level of reduction of Notch is not quantified in the experimental settings, we first fitted a multiplicative magnitude of Notch reduction that reproduced the data. With a Notch reduction of 0.4 × WT, our model is able to reproduce both mutant JU2039, where EGF is WT and there is no strong phenotype, and mutant JU2113, where EGF is slightly increased by 1.25-fold and secondary fate in P5.p is destabilized, increasing the probability of primary and tertiary fates. This fitted fold increase is, in fact, similar to the one considered in [5]. The model predicts that, under JU2039 signaling regime, P5.p would stay in a region of the signal space where secondary fate is still predominant, and therefore the pattern would not change. Although our model is not capable of reproducing a small number of tertiary fated P5.p cells under the JU2113 perturbation, since the signaling profile moves them slightly outside the tristability region, it does reproduce the increase in primary fated P5.p, by approaching the cusp mid-line in the signal space while not leaving the bistability region inside the cusp (Table 2 and Fig 6A).

We also tested whether the model is able to reproduce mutants with ectopic Notch expression. Similar to the previous subset of data, these perturbations are silent unless crossed with EGF overexpression perturbations. As in [5], we modeled these perturbations as an additive constant for the Notch signaling parameter as, in this case, the perturbation is independent of the state of the cell and its neighbors (see S1 Appendix). We fitted this constant to the experimental outcome of the JU2092 mutant, where ectopic Notch is combined with a 1.25-fold EGF overexpression. This fitting resulted in a constant equal to 0.12. With this, we were able to fit the three crosses shown in Table 2. Interestingly, in the strongest EGF overexpression perturbation (mutant JU2092), P4.p was positioned close to the right of the intersection between the cusp mid-line and the fold line, giving a mix of tertiary and secondary fates, as observed in the data (Fig 6A).

Finally, we simulated EGF hypomorph mutants. Interestingly, these perturbations show that a strong decrease of EGF signals gives a mix of primary and tertiary fated P6.p cells while not affecting the WT pattern of the rest of VPCs. Moreover, a cross with a mutant with mild ectopic Notch activity, such as the one considered above, promotes primary fated P6.p, showing again epistasis between the signals. To simulate these perturbations, we fitted a multiplicative downregulation of EGF signal to the EGF hypomorph data, CB1417. With a reduction of 0.36 × WT EGF, our model can reproduce the mix of tertiary and primary fated P6.p observed in the EGF hypomorph mutant (CB1417) (Fig 6C). This is achieved because a reduction in EGF positions the signaling profile of P6.p close to the fold line while still being to the left of the cusp-midline (Fig 6B). Simulating the cross of this perturbation with ectopic Notch (mutant JU2095) moves P6.p signaling profile further from the fold line, promoting vulval fates. However, it also positions it closer to the cusp mid-line and therefore, also promotes secondary fates in P6.p in our simulations (Fig 6B). The model is also not able to reproduce the high probability of P5.p staying in tertiary fate. We believe that the differences between simulations and experiments observed in P6.p could be fixed by adding the perturbations to the training data set. Reproducing the differences observed in P5.p would likely require a non-linear mapping of the fold into the signal space.

Model predictions

Having reproduced a large set of data, we explored whether the model could predict interesting outcomes.

We tested whether crossing two silent perturbations in the training data could show a new epistatic event. Halving EGF ligand or Notch receptors does not affect the wild type phenotype (Table 1 (3–4)). In our model, this is because making either of these two perturbations leaves cells within the same stability regions of the fate map. However, our model suggests that, if both perturbations are crossed, we observe epistasis, where induction of P5.p to secondary fates is strongly reduced by almost half (Fig 7A). This is consistent with predictions shown in [3]. Biologically, this prediction suggests that EGF promotes secondary fates which, at first, could seem counterintuitive.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Model’s predictions.

(A) Predicted epistasis under the cross of reduced EGF ligand and Notch receptors perturbations. Top: Simulated proportions under the corresponding signaling profiles. Bottom: Schematics of the simulated positions of P5.p in the signal space under the different signaling profiles. (B) Predicted pattern proportions for different AC ablation times under a reduced EGF ligand perturbation. t* is a predicted time between 3° divided and 2-cell stages. The stages for which simulations are not shown are colored in gray. (C) Predicted pattern proportions for different AC ablation times under a reduced Notch receptors perturbation. The names of the stages for which simulations are not shown are colored in gray.

https://doi.org/10.1371/journal.pcbi.1009034.g007

AC ablation experiments are very insightful, as they are the only experiments that provide information about the trajectory of the cells in the landscape. However, AC ablation data is only available under a WT signaling profile. As we mentioned earlier, a wild type pattern is achieved even under either half dose lin-3 or lin-12. However, it is not known whether the pattern is dynamically formed in the same way as under WT signaling. We used our model to test this and observed an interesting effect.

Simulating AC ablation under a half dose of lin-3, the EGF ligand, showed that the process followed the same order as under WT signal, but proceeded more slowly (Fig 7B). Under this condition, our model predicts (1/3, 1/3, 1/3) probabilities for P6.p becoming one of the three fates if AC is ablated at a time between 3° divided stage and 2-cell stage. However, under a half dose of lin-12, a Notch receptor, our model predicted that there is very low probability of P6.p becoming a secondary fate under any ablation time, as observed in the WT and lin-3 mutant (Fig 7C). Instead, there is a direct transition from tertiary to primary fate. This is also consistent with the effect observed in [3].

Finally, we explored whether there was a prediction that could highlight the differences between the model presented here and that proposed in [3, 5]. An important difference between the two models is the fact that the binary flip with cusp model is not symmetric, in that once cells have sufficiently left the basin of attraction of the tertiary fate they cannot go back for any value of the external signals. For example, if isolated cells (with no autocrine or paracrine Notch signaling) are exposed to a short EGF pulse of the same WT strength (from time t₁ = 0 to t₂ = 0.4), our model predicts that most of them would leave the basin of attraction of tertiary fate and commit to either the primary or secondary fate (S3(A) Fig). Moreover, due to the directions of the flow lines, there is no signaling profile that would be able to return cells back to tertiary fate. For example, if this short EGF pulse is followed by a very short Notch signaling pulse (t₂ = 0.4 to t₃ = 0.62), cells remain outside the basin of attraction of the tertiary fate and would differentiate into the secondary fate (S3(A) Fig). However, this is not the case in the model proposed by [3, 5]. As in our model, isolated cells exposed to the same short EGF pulse leave the tertiary fate and commit to primary fate, resulting in 90% of primary-fated cells (S3(B) Fig). Strikingly, if this is followed by the short Notch signaling pulse, the cells return to the basin of attraction corresponding to the tertiary fate, resulting in 71% now differentiating to the tertiary fate and only 19% differentiating to the primary fate (S3(B) Fig). Performing this or a similar experiment would be able to differentiate between the two models.