Whole-embryo modeling of early segmentation in Drosophila identifies robust and fragile expression domains

Abstract

Segmentation of the Drosophila melanogaster embryo results from the dynamic establishment of spatial mRNA and protein patterns. Here, we exploit recent temporal mRNA and protein expression measurements on the full surface of the blastoderm to calibrate a dynamical model of the gap gene network on the entire embryo cortex. We model the early mRNA and protein dynamics of the gap genes hunchback, Kruppel, giant, and knirps, taking as regulatory inputs the maternal Bicoid and Caudal gradients, plus the zygotic Tailless and Huckebein proteins. The model captures the expression patterns faithfully, and its predictions are assessed from gap gene mutants. The inferred network shows an architecture based on reciprocal repression between gap genes that can stably pattern the embryo on a realistic geometry but requires complex regulations such as those involving the Hunchback monomer and dimers. Sensitivity analysis identifies the posterior domain of giant as among the most fragile features of an otherwise robust network, and hints at redundant regulations by Bicoid and Hunchback, possibly reflecting recent evolutionary changes in the gap-gene network in insects.

Figure 1

Gap-gene-network model. (A) Experimentally measured maternal protein gradients (Bcd and Cad) and protein expression of the terminal system proteins (Tll and Hkb) are taken as time-varying inputs to the gap-gene network. (B) knirps mRNA at C14A, with the mesh defined by the positions of the nuclei. High and low expression are indicated by red and blue, respectively (the full data set used in this model is shown in Fig. S4). Data are taken from the BDTNP database. (C) Modeling geometry. A-P denotes the anterior-posterior axis and D-V the dorsoventral axis. The thick band along the A-P axis (from EL position 35% to 92%) shows the geometry considered in 1D models (23). Unless stated otherwise, all embryos will be presented in this orientation, and 1D plots will be along the A-P line. (D) The model. mRNA and protein expression levels for the gap genes hb, Kr, gt, and Kni are modeled on the embryo surface; all mRNA and protein species can diffuse (Methods and Fig. S10). Each gap-gene mRNA is transcribed according to a linear model of transcriptional influences (u is a linear combination of the protein expression levels, p_j). A nonlinear transfer function, g(u), models saturation of the polymerase. The proteins are translated from mRNA using a linear model, and all degradations are first-order processes.

Figure 2

Model for the gap gene network calibrated on the whole surface of the embryo. mRNA patterns. (A) Experimental data (red or dark gray) and best fit (green or light gray) are superimposed for three gap genes. Yellow indicates that data and simulations agree perfectly. Notice that the giant domains are much improved compared to Fig. S12; in fact, the agreement is very good up to 20 min and then deteriorates when the head patterns become very fine. (B) Distribution of best-fit parameters across 21 independent solutions (all solutions have RMS error $\sqrt{S\left(\theta \right)/N}<36$; the maximum expression is set globally to 255). Here, only the interaction matrix T^ij is given. Each column shows regulation of one trunk gap gene by the nine regulators. (C) Reconstructed gene network. Gap genes are shown as circles, where green (light gray) indicates self-activation. Input genes are shown above and below. Activating (respectively inhibitory) links are in green or light gray (respectively, red or dark gray). Darker shading indicates smaller errors on the parameters as computed from Hessian matrix (darker shading indicates smaller errors on the parameters) as computed from the Hessian matrix (Section 6 in the Supporting Material). The T matrix and other inferred parameters are given in Table S5 and Table S6.

Figure 3

3D model predicts patterns of gap-gene mutants. The line-graph data for Giant (A–C) and Knirps (D) patterns in mutants (Gt data) are quantified from Kraut and Levine (15) and Hülskamp et al. The times for the mutant stainings correspond to 30 min. The A-P projections of the modeled proteins at 30 min are also shown (Model), as is a simulation from a published model (Circuit 28008 in Jaeger et al. (18)) (Model 2). In the side views, simulated mRNA patterns for the wild-type (green or light gray) are shown together with the simulated mutant pattern (red or dark gray). Yellow or white indicates no change between the simulated wild-type and the mutant patterns. (A) hb null mutant. The absence of the anterior gt mRNA domain is correctly predicted up to t = 40 min. In a similar way, the posterior extension of the posterior domain is also consistent. (B) Kr null mutant. The gt pattern for this mutant is poorly predicted. (C) kni null mutant. The vanishing of the posterior gt domain is correctly predicted. (D) Kni in hb null mutant. This mutant is correctly predicted.

Figure 4

Gap gene network is robust and fragile. (A) Cluster of stiffest eigenvectors. The mean and standard deviation of the cluster are indicated. Arrows indicate the important parameters in that mode (1–6), from left to right (see also C): activation of Kr by Cad (1), repression of Kr by Tll (2), activation of gt by Cad (3), repression of gt by Tll (4), mRNA and protein production of Kr (5), and mRNA and protein production of gt (6). (B) Cluster of softest eigenvectors. The mean and standard deviation of the cluster are indicated. Arrows indicate, from left to right, activation of kni by Hkb (7), and repression of hb (8) and Kr (9) by Hkb. (C and D) Perturbing the network along stiff (C) and soft (D) directions. (C, left) Perturbation associated with the stiffest eigenvector (for the solution in Fig. 2) is shown on the network model. Numbers refer to important parameters defined in A. The perturbation is taken as ${\tilde{p}}_i=p_i+q_i$, where p_i is the optimum and q_i = εv_i is a relative perturbation along the eigenvector $\overrightarrow{v}$, with $\left|\overrightarrow{q}\right|=1$. Components of the eigenvectors determine the color intensity of the links. Green or light gray indicates that the perturbation increases the magnitude of the regulation (positive or negative) and red or dark gray for the opposite. This eigenvector involves mainly the control of Kr and gt by maternal inputs Tll and Cad. (C, right) Simulation of wild-type (green or light gray) and perturbed parameters, ${\tilde{p}}_i$ (red or dark gray), indicates the high stiffness of this mode, as seen by the loss of the posterior gt domain. Yellow or white indicates that the perturbation has no effect. Here, the mRNA at 20 min is shown. (D) Same as in C, but for the softest mode. The perturbation has no effect in this case, even though the perturbation, $\overrightarrow{q}$, has the same magnitude (norm) as in C. Numbers refer to parameters defined in B.