Cells change their sensitivity to an EGF morphogen gradient to control EGF-induced gene expression

How cells in developing organisms interpret the quantitative information contained in morphogen gradients is an open question. Here we address this question using a novel integrative approach that combines quantitative measurements of morphogen-induced gene expression at single-mRNA resolution with mathematical modelling of the induction process. We focus on the induction of Notch ligands by the LIN-3/EGF morphogen gradient during vulva induction in Caenorhabditis elegans. We show that LIN-3/EGF-induced Notch ligand expression is highly dynamic, exhibiting an abrupt transition from low to high expression. Similar transitions in Notch ligand expression are observed in two highly divergent wild C. elegans isolates. Mathematical modelling and experiments show that this transition is driven by a dynamic increase in the sensitivity of the induced cells to external LIN-3/EGF. Furthermore, this increase in sensitivity is independent of the presence of LIN-3/EGF. Our integrative approach might be useful to study induction by morphogen gradients in other systems.

(E) and (F) for wild-type animals and for lin-12(lf) and lin-3(++) mutants. The maximally induced mRNA level is defined as the mean of the mRNA in P6.p of all animals in the late induction stage, >9 hr after the start of L2. Error bars indicate the standard error of the mean. Each experiment represents data for animals. Even though in lin-3(++) animals internal Ras signaling in P6.p is expected to be hyperactivated, the measured values of the maximally induced mRNA levels in these mutants is similar to wild-type animals. This confirms our hypothesis that the measured mRNA levels in the late induction stage indeed correspond to the maximally induced mRNA levels.
(A) Mean lag-2 expression level (bars) and standard error of the mean (error bars) in all VPCs, P(3-8).p, in the lin-1(0) mutant ( animals). Shown are the expression levels at early induction (cyan bars, 1-4hr after the start of L2) and late induction (magenta bars, 9-12hr after the start of L2). (B) Same as (A) but for apx-1. (C) Overview of lag-2 and apx-1 expression in the lin-1(0) mutant. Columns correspond to individual VPCs. Rows represent different animals and are sorted according to increasing gonad length. Colors represent the mRNA level with the same scale as used in Fig. 2C and D in the main text. (E) Same as (A) but for the let-60(gf) mutant ( animals). (F) Same as (A) but for the lin-1(0);lin-3(e1417) mutant ( animals). (G) Same as (C) but for the let-60(gf) mutant. (H) Same as (C) but for the lin-1(0);lin-3(e1417) mutant.

Supplementary Figure 5. Overview of the best fits of Models A, B and C.
Best fits of the 8 models A1-A3 (red bars), B1-B3(green bars) and C1-C2 (blue bars) to the experimental data for lag-2 expression in P6.p in wild-type, lin-3(e1417), lin-3(++) and lin-1(0) animals as well as the model predictions for the lin-1(0);lin-3(e1417) mutants and a hypothetical mutant lacking the activator A. For the experiments, mean expression level (grey bar, animals) and standard error of the mean (black error bars) are shown. In contrast to Models A1 and B1, Model C1 is able to reproduce the experimental data for the lin-3(++) mutant, even though in this model the transition is due to a change in LIN-3 level. However, in contrast to Models A2-3, B2-3 and C2 this agreement only occurs for extremely precisely tuned model parameters and hence in the main text we already discard Model C1 based on lack of agreement with the lin-3(++) data.
(A) Single-cell correlation between apx-1 and lag-2 mRNA levels in P(5-7).p in wild-type animals (green, n=73) and P(4-8).p in lin-12(lf) animals (red, n=111). Line indicates the fit to Eq. (1) in the main text for wild-type animals with . (A) Overview of parameter combinations that yield good fits to the experimental data. Grey markers indicate parameter combinations that are good fits to the lin-12(lf) data, but not to the lin-3(++) and the lin-1(0) data. Blue markers indicate parameter combinations that give good fits to the lin-12(lf) and lin-3(++) data, but not to the lin-1(0) data. Green markers indicate parameter combinations that are good fits to the lin-12(lf) and the lin-1(0) data, but not to the lin-3(++) data. Finally, red markers indicate good fits to all three data sets. The overall best fit to the data, as characterized by the total RMSE, is indicated by the white marker. (A) Schematic description for Model A of the dependence of the lag-2 transcription rate r on the binding of the repressor LIN-1 to the binding site and the activator A to the binding site . Transcription only occurs, at the maximum rate r=1, when the activator A but not the repressor LIN-1 is bound to the promoter (solid box). (B) Schematic description for Models B and C of the dependence of the lag-2 transcription rate on the binding of the repressor LIN-1 to the binding site , the activator LIN-1-P to the binding site and the activator A to the binding site . When both activators LIN-1-P and A are bound transcription occurs at the maximum rate r=1 (solid box). When either LIN-1-P or A are bound, transcription occurs at a reduced rate (dashed box).

Supplementary Note 1. LIN-3 secretion and binding to LET-23
We assume that both the LIN-3 morphogen gradient dynamics and the LIN-3 binding kinetics are in steady state. With this assumption we arrive at the following diffusion equation for the free LIN-3 concentration , where is distance along the anteroposterior axis to the Anchor Cell (AC), the source of LIN-3: Here, is the LIN-3 diffusion constant, is LIN-3 clearance rate, is the total concentration of the EGF receptor LET-23 and and are the forward and backward rate of binding of LIN-3 to LET-23, respectively. In addition, secretion of LIN-3 by the AC at a rate is incorporated using the boundary condition: The steady state equation for the concentration of the activated LIN-3-LET-23 complex, , given by: Solving Eq. (3) yields: (4) Inserting this expression in Eq. (13) yields: Solving Eqs. (5) and (2) gives: where the decay length is given by: Using this expression, we find the following expression for : where (9) is a parameter that we label the LIN-3 input strength and that describes all extracellular dynamics of LIN-3 secretion and binding to the LET-23 receptor. Then, | | is the fraction of total free LIN-3 at position , with ∫ .

Supplementary Note 2. Induction of Ras signaling by activated LET-23
The Vulva Precursor Cell (VPC) integrates the concentration of activated LET-23 over the entire basolateral surface of the cell. Hence, in our model the Ras signaling strength for a VPC at distance to the AC is given by the total concentration of activated LIN-3-LET-23 complexes integrated over the body of the VPC: where is the length of the VPC body along the A-P axis. Using Eq. (8) we arrive at the following expression of the Ras signaling strength: We calculate analytically using the expression: where , and ( ).

Supplementary Note 3. Mathematical modeling of Notch ligand expression dynamics Model A
We assume that LIN-1 is phosphorylated in response to Ras signaling with rate , with given by Eq. (11), and spontaneously dephosphorylates at rate . The activator A is converted from its inactive form A* to the active form A with forward rate and converted back at rate . On the timescale of induction these reactions are in steady state, leading to the following equations for the activation of the transcription factors LIN-1 and A: (13)   (14) where [Lp] indicates the concentration of LIN-1-P, [L] the concentration of LIN-1, and . We assumed that the transcription rate of Notch ligand mRNA depends only on the probability of the transcription factors LIN-1 and A being bound to their respective binding sites and according to Fig. 8, constituting a so-c d "th r ody ic od " of g r ssio . I this ictur , th st dy st t bi di g rob bi iti s of th transcription factors LIN-1 and A are given by: where and are the concentration of LIN-1 and A at their binding sites and and are their respective binding constants. Transcription occurs at a rate when A but not LIN-1 is bound, reflecting their respective roles as activator and repressor. This gives rise to the following expression of the transcription rate: To calculate an explicit expression for the transcription rate, we use the following conservation equations: where and are the total concentrations of LIN-1 and A in the VPC. Using the above equations, we calculated the following expressions for the concentration of the different transcription factor species LIN-1, LIN-1-P and A: and the following expressions for the concentration of unbound promoter binding sites: yielding the following expressions for the concentration of bound promoter binding sites and : This results in the following expression for the transcription rate for model A: where:

Model B
In Model B, LIN-1 and A still act as a repressor and an activator, as in Model A, but now in addition LIN-1-P also activates Notch ligand expression (See Fig. 8). The concentrations of LIN-1, LIN-1-P and A are given by Eqs. (13)-(14). Similar to Model A, the steady-state concentration of transcription factors bound to their promoter binding sites are given by: We assumed that no transcription occurred when the repressor LIN-1 was bound (Fig.  S8, ), that transcription occurred at a low level when either LIN-1-P or A were bound and at the maximum level when both LIN-1-P and A were bound. This gives rise to the following expression for the transcription rate: where (32) and (33) For the concentration of unbound promoters we have the following expressions: leading to the following expression for the transcription rate for Model B: where:

Model C
Model C is very similar to Model B, with the sole change that now the conversion of the activator A from its inactive to its active form is dependen on Ras signaling, i.e. proceeds with a forward rate . This gives rise to the following modified expression for the activation of the transcription factors LIN-1 and A: Following the approach for Model B, we find the following expression for the transcription rate for Model C: where:

Supplementary Note 4. Summary of Models A, B and C for wild-type and mutants
In the table below, we summarize the expressions for the Notch ligand transcription rate in Models A, B and C. We also show the expressions for the transcription rate in the lin-3(0) mutant, corresponding to , as well as in the lin-1(0) mutant, corresponding to and the mutant lacking the activator A, corresponding to :

Supplementary Note 5. Overview of parameters in Models A, B and C
Parameter Definition Ras signaling strength Distance of VPC to AC LIN-3 gradient decay length LIN-3 input strength Sensitivity of Ras signaling pathway to external LIN-3 Binding affinity of repressive LIN-1 to lag-2 promoter Binding affinity of activating LIN-1-P to lag-2 promoter Transcription rate of either LIN-1-P or A bound relative to transcription rate when both LIN-1-P and A bound Binding affinity/amount of activator A

Supplementary Note 6. Fitting models A, B and C to lin-3(++) and lin-1(0) mutant data
To examine how the observed transition in Notch ligand expression is generated on the network level, we tested, for each of the Models A, B and C, three potential underlying mechanisms: first, the case where the transition was driven by a change in external LIN-3 level (Models A1, B1 and C1 in Fig. 5A in the main text). Second, the case where the transition was driven by a change in sensitivity of the Ras pathway to the external LIN-3 signal (Models A2, B2 and C2). Third, the case where the transition was driven by a change in the amount of activator A (Models A3 and B3). We implemented this in the context of the models A, B and C by assuming that all parameters except for a single parameter are constant in time. For Models A1, B1 and C1, then the only parameter changing in time is . For Models A2, B2 and C2 the only parameter changing in time is , rewriting Eqs. (11) and the expression for in Eqs. (27) or (38) so that , with . Finally, for Models A3 and B3 the only parameter changing in time is . Below, as an illustration we will outline the procedure for fitting Model B3, where varies in time, to the lin-3(++) and lin-1(0) data for lag-2 ( Fig.  5B-D in the main text). The procedure for fitting the other models is identical, apart from the fact that a different time-dependent parameter is used.
We perform our fitting procedure under the constraint that the model should reproduce the wild-type data for lag-2 and show no induction of lag-2 expression in the absence of LIN-3, i.e.
. The first constraint is given by: where and are the values for at the early (1-4hr) and late (9-12hr) induction and and are the mean lag-2 mRNA levels for early and late induction (Fig. 5B in the main text). In Model B3, we still have basal induction for , with the magnitude depending on the value of and . Here, we constrain so that in the absence of LIN-3 the lag-2 expression level is less than 1% of the fully induced level: . As an initial guess for the fitting procedure, we picked values for the parameters and , so that the constraints in Eqs. (43)and (44) were satisfied. We then fit the model to the lin-3(++) and lin-1(0) data using the Sequential Least Squares Programming (SLSQP) constrained optimization algorithm to minimize the following error function: where the index indicates the different mutants, indicates the different time points, is the mean lag-2 mRNA level of mutant at time point ( Fig. 5C and D in the main text) and is the lag-2 mRNA level at the late induction timepoint. The lag-2 transcription rate for the lin-3(++) mutant is given by -, corresponding to a tenfold increase in LIN-3 dosage. The lag-2 transcription rate for the lin-1(0) mutant is given by -, corresponding to the complete absence of LIN-1 and LIN-1-P. The constrained optimization algorithm ensured that the constraints in Eqs. (43) and (44) are maintained throughout the fitting procedure. In the case of Model B3 we found that we could simultaneously fit the wild-type, lin-3(++) and lin-1(0) mutant data (Fig. 5), with many parameter combinations giving equally good fits. We interpret this as Model B3 being consistent with the experimental data. However, for the other models we found that we were unable to fit the data for wildtype animals and one or both mutants simultaneously. We interpreted this as those models being incomplete or incorrect.

Supplementary Note 7. Regulation of lag-2 and apx-1 expression in model B3
The expression for the transcription rate for Model B3 can rewritten as: where: where is the only time-dependent variable. The mRNA level L is given by , where d is the mRNA degradation rate. For full induction, , the expression for the maximally induced mRNA level is given by , whereas in the absence of activator A, , the basal mRNA level is given by . We assumed that the only parameters that can control the difference in expression patterns observed between lag-2 and apx-1 are those governing the interaction between LIN-1, LIN-1-P and A and their binding sites in the lag-2 and apx-1 promoters, i.e. and . However, and only influence the basal, , and maximum expression level, , not how the expression level increases as a function of the amount of activator. Hence, to capture the difference in time dynamics between lag-2 and apx-1 expression we have to assume that the binding of A to the lag-2 and apx-1 promoter occurs with different affinity. For this we use the simplest assumption that for each Notch ligand i. Finally, to capture the difference in maximum expression level between lag-2 and apx-1, we use the simplest assumption that only one parameter, which we choose to be the LIN-1-P binding rate , differs between lag-2 and apx-1. Hence, we need to consider two values for each Notch ligand i, with the values of the remaining parameters and the same for both the lag-2 and apx-1 promoter.

Supplementary Note 8. Single-cell correlation between lag-2 and apx-1 expression
Using Eq. (46), we can write the expression level of Notch ligand as: Solving for gives: If the two ligands and are coregulated by same activator A, as we assume is the case for lag-2 and apx-1, we can insert the expression for into the expression for , yielding: Experimentally, we observe no expression of lag-2 and apx-1 in VPCs at the start of vulva induction, i.e. for i={lag-2,apx-1}. Hence, the parameters , and are tuned so that the basal expression levels vanish, , and we have: where we have , corresponding to Eq. 1 in the main text. In the context of Model B3, we interpret as follows: the conversion of active to inactive A, as given by the equilibrium constant , is the only process that depends on time. At the same time, the binding affinity of active A is constant in time, but different for each promoter. With these assumptions and using the expression for in Eq. (38), we get the following expression for : Thus, in this interpretation ̃ ̃ , which is the ratio between the thresholds for induction of Notch ligand i and j by the activator A. Here, would mean that Notch ligand i has a lower threshold to induction by Ras signaling than Notch ligand j.
Finally, we note that the maximally induced mRNA level in a specific VPC depends on the Ras signaling strength and hence on the distance of the VPC to the AC. This means that Eq. (51) is only expected to hold when we use the value of measured in that particular VPC. In Fig. 6 in the main text, we use for the maximally induced lag-2 and apx-1 mRNA levels observed in P6.p. Hence, a priori it is not expected that the single-cell correlation curve between lag-2 and apx-1 expression for cells other than P6.p should collapse on the same curve. Still, we find excellent overlap between curves for P6.p in wild-type animals and P4.p, P5.p and P6.p in lin-12(lf) animals. However, this could be responsible for the relative deviation from the theoretical curve observed for P7.p cells in the lin-12(lf) mutant.

Supplementary Note 9. Overview of parameters of the extended Model B3
Parameter Definition Ras signaling strength Distance of VPC to AC LIN-3 gradient decay length, time-dependent LIN-3 input strength Sensitivity of Ras signaling pathway to external LIN-3 Binding affinity of repressive LIN-1 to lag-2 and apx-1 promoter Binding affinity of activating LIN-1-P to promoter of Notch ligand i Transcription rate of either LIN-1-P or A bound relative to transcription rate when both LIN-1-P and A bound Binding affinity of activator A to the promoter of Notch ligand i, time-dependent Ratio of binding affinitys of activator A for the lag-2 and apx-1 promoter, lag-2 and apx-1 mRNA degradation rate Supplementary Note 10. Fitting Model B3 to experimental lag-2 and apx-1 expression data We found that the extended Model B3 has four independent global parameters: and . In addition, for each time point in Fig. 4A and B in the main text the model has two additional independent parameters, and . The remaining parameters and are fixed by the experimental data in the following way: no induction of Notch ligand expression should occur in the absence of LIN-3 or far away from the AC, i.e. when , even when A is fully activated, i.e. . This means that: where we use . Equation (53) uniquely fixes the value of . Because lag-2 is more highly expressed than apx-1 this value of also ensures no apx-1 expression occurs in the absence of LIN-3. However, not all combinations of and allow for a value of that satisfies Eq. (53). This only occurs when the following condition is met: The value of is fixed by the constraint that the model is able to reproduce the fully induced apx-1 expression level : However, also in this case not every combination of and allows for a value of that satisfies Eq. (55). This only occurs when (56) and (57) are both satisfied. We generated combinations of the parameters and , as outlined in the Methods section in the main text, that satisfy the above criteria. Then for each time point in Fig. 4A and B, we found values of and that minimize the SSE between the spatial expression profile predicted by Model B3 for lag-2 and apx-1 expression and the experimental data in Fig. 4A and B. In addition, we calculated for each parameter combination the SSE with respect to the lag-2 expression dynamics in the lin-3(++) and lin-1(0) mutants, as outlined in the section 'Fitting models A, B and C to lin-3(++) and lin-1(0) mutant data'.