The precise results of non-linear interactions with positive and negative feedback, such as the genetic circuit described above, are difficult to predict by intuition alone. Therefore we took advantage of mathematical modeling to determine whether the interactions that we described could achieve reliable conversion of a graded signal into a step function. Details of the postulated network evolved during attempts to simulate our experimental results based on the modeling, including the mutant situations. Deviations between the observations and the simulations allowed us to refine the model.
Equations such as those presented in are easy to understand conceptually even for non-specialists. Each equation describes the concentration change of a component in a short time interval as function of the concentrations of the other players. These changes depend upon each protein’s production or activation rate and its decay rate, and in the case of UPD, upon gain or loss by diffusion to neighboring cells. The rates depend on regulatory influences of the other components. In the computer simulation these changes are added to an initial concentration, providing the distributions at a slightly later time point. Subsequently, the next concentration changes are calculated, added to the updated concentrations, and so on. In this way, the total time course of the system is calculated. The initial condition is that the polar cells are specified and produce a polar cell activator (P), which is set arbitrarily to one. It remains unchanged. All other concentrations are initially set to zero. First consider the simplest case, that for the morphogen UPD (U) in Equation 2. UPD is produced in proportion to P, and has a first order decay, i.e., the number of molecules that disappear per time unit is proportional to the number of molecules present. The last term in this equation reflects the diffusion of the molecule where x represents distance. In Equation 1, JAK/STAT activity (J) is increased in proportion to the UPD signal (U). Activity saturates at high levels, giving us 1/0.1 J2, which allows achievement of a stable high level. APT (A) attenuates JAK/STAT activity and thereby appears in the denominator. The downregulation of APT on JAK/STAT is antagonized by SLBO (S). To reflect the empirical observation that there is no JAK/STAT activity detectable in the polar cells, we assume that the polar cell activator (P) attenuates signaling in the polar cells. The decay of JAK/STAT activity is described by the second term in the equation. Equation 3 represents the change in APT concentration. APT (A) is produced in response to JAK/STAT activity (J), but this effect saturates. APT is also downregulated by SLBO (S). The last term in the equation accounts for the JAK/STAT-independent production of APT in response to EYA. In Equation 4, we see that JAK/STAT activates SLBO expression non-linearly, but APT downregulates it. By including a second-order dependence of SLBO on JAK/STAT (J2), the distribution of SLBO expression becomes steeper than that of APT (which is merely a function of J). SLBO also has a linear decay rate.
Figure 4 Differential equations describe the establishment of biological patterns. (A) Differential equations describing the changes in concentration over time of each pathway component as a function of the concentrations of other components. See text for details. (more ...)
Since we do not have precise kinetic data, numerical values in the equations are estimates and can vary without changing the outcome, as long as the APT:SLBO balance is preserved. In this model, we achieve a steady state at approximately 4,500 iterations. Since the cell-specification process occurs over about 6 hours, each iteration accounts for about 5 seconds. In , the final stable concentrations are plotted as a function of the distance from the polar cells. An advantage of this model system is the small number of cells involved and their arrangement in a single layer. For simplicity they can be considered as a linear array (although, in reality, the signaling occurs radially out from the polar cells in three dimensions). The model clearly recapitulates the signaling pattern across the field of cells, with the polar cells in the center, and the border cells on either side.
The computer simulation accurately mimics the wild-type case as well as several mutant situations that were tested; in each case the outcome matched our empirical observations. In wild-type, SLBO is activated at high levels of STAT activity and antagonizes APT. Near the polar cells, APT activity is reduced, leading to an increase in STAT activity, which leads to an even higher level of SLBO, and thus lower APT activity so that STAT and SLBO predominate. Further from the polar cells, APT outcompetes STAT and SLBO due to its distribution, which is a result both of the saturation of APT production at relatively low concentrations of STAT and of the pattern of EYA expression. Thus, the ON-OFF behavior (to become a border cell or not) that we see in wild-type does not appear to depend on JAK/STAT auto-regulation but instead requires the switch-like effect of two mutually repressing genes, SLBO and APT. The model also describes correctly how the step-like JAK/STAT distribution evolves: first a graded distribution emerges. Subsequently, the cells distant to the polar cells downregulate STAT activity. In other words, cells that are not exposed to a threshold morphogen concentration switch off the migratory cell fate. This is very different from other threshold systems in which the ON-switch only occurs if the signal concentration is above a certain threshold and is maintained even in the absence of the morphogen (see ). This property may be important as the anterior follicle cells also require JAK/STAT activity for their fate specification. These cells experience not only a lower level of signaling compared to the border cells but also a transient signal. Moreover border cell fate requires continuous JAK/STAT signaling.66
If STAT activity is lost after they initiate migration, the cells stop migrating and turn on expression of at least one marker of nurse cell associated follicle cell fate. This system is considered therefore to be monostable. The only stable state is off. In order to keep the signal on, continuous input is required.
The model also simulates mutant phenotypes correctly. For example, if APT is absent (), the JAK/STAT distribution remains graded similar to its activator, UPD. Thus more cells achieve high levels of STAT signaling and acquire border cell-like properties. In contrast, if APT is over-produced, JAK/STAT activity eventually becomes repressed everywhere and no cells migrate. We tested the effect of co-overexpression of both APT and SLBO in silico, and the computer simulation suggested that the migration defect caused by overexpression of either gene alone would be rescued, as we observed in vivo.
A short PC-program including well-commented source code is available in the supplementary material. The program allows the simulation of wild type and mutant situations. Since the program can be recompiled with an open-source compiler, alternative molecular interactions or other ranges of parameters can be easily tested.
Several interesting concepts emerge from the mathematical modeling. First of all, it demonstrates that our minimal four-component genetic interaction network is at least theoretically sufficient to describe the stable outcomes we see in vivo with respect to cell fate decisions, and thereby, cell behaviors. Secondly, it accurately predicts the outcomes of several experiments in which we varied the levels of the components. But perhaps most interestingly, it hints at slightly different biological mechanisms than we expected. For example, we found that STAT does not require auto-regulatory activity for border cell specification, even though STAT activity is known to enhance its own expression in the embryo.41,42
Early versions of our model did incorporate a positive feedback factor for STAT activity (), and inclusion of such a term can simulate wild-type patterning and produce a switch-like behavior (). However, this model does not accurately simulate the apt
mutant phenotype, as it shows expanded signaling but still retains the ON/OFF property (). There are several possible explanations for this apparent conundrum. STAT auto-regulation may be tissue-specific and not occur in border cells, or may occur at low levels. Alternatively, this auto-regulation may occur more slowly than the critical time frame between when the border cells are specified and when they move away. Finally, in the ovarian epithelium, the localized and continuous production of UPD from the polar cells defines a domain of maximal activation, which reduces a need for auto-activation for the generation of pattern. Thus, the STAT positive feedback loop may be important in signaling during embryogenesis, but may be insufficiently rapid or simply unnecessary for specifying the border cell population in the presence of APT action.
Another surprise was that cells below a threshold become deactivated, in contrast to many other reactions in which activation is irreversible. In the case of border cells, and also in Drosophila testis stem cells,67,68
it is important that those cells receiving a low level of UPD signal quickly extinguish STAT activity. It may be that the detailed properties of the JAK/STAT signaling pathway make it particularly suitable to this type of regulation, which is also known as monostability. Therefore, the mechanisms we are studying may reveal fundamental properties of JAK/STAT signaling and monostable circuits in general.
Some features of the genetic circuit still require an explanation. For example, the apt loss of function phenotype is incompletely penetrant. That is, 80% of apt mutant egg chambers have additional migratory cells and expanded STAT signaling. But the other 20% appear normal. A possible reason could be that the mechanical separation of the border cell cluster from the epithelium, which physically removes the source of activating UPD signal, may sometimes be sufficient to result in loss of low STAT activity. Alternatively, or in addition, STAT autoregulation in the absence of APT could play a role. Experiments are underway to test these alternatives by genetically blocking border cell movement and measuring the effects on the STAT activity gradient. Note that in the mathematical models, one goal is to find conditions that lead to a steady state. In the case of the border cells, the activating signal is physically moving, but signaling properties altered by the movement are not factored into the simulations. In light of this, models incorporating the putative non-linear STAT self-enhancement may be correct, but the graded signals observed in vivo reflect the early iterations, prior to the achievement of a steady-state, (i.e., the first panel in ).
In conclusion, our molecular genetic investigation of border cell determination together with the modeling revealed a new type of all-or-nothing determination under the influence of graded morphogenetic signal. The modeling may provide help to understand other systems that still appear to be puzzling. Further investigation will reveal whether related mechanisms are also employed in other biological situations.
Program code and parameters:
CASE 246 ‘GE-BISTABLE STATE, FOR POLAR-BORDER CELLS
‘SLBO does not only hinder the production but the effect of APT on J/S
‘a = STAT, activated by UPD (term ba * b), inhibited by APT; thin inhibition is undermined by SLBO
‘J/S is excluded by pole cells (term cc * c)
‘b = UPD, produced by pole cells (=c)
‘c = pole cells, defined by initial condition, unchanged
‘d = APT is quenching the signal (a) that is produced J/S (=> a); is inhibited by SLBO (term cd * e * e)
‘and by itself (term sd * d * d), both terms are non-linear; a base-line APT production exists (se)
‘e = SLBO, produced by j/S (a), inhibits and is inhibited by APT (d)
‘e.g., The term ‘olddecaydiffA’ is the new concentration of a as given by the change due to decay (and diffusion)
‘e.g., olddecaydiffA = a *(1 - ra) where a is the removal per iteration
‘DA, DB… diffusion; RA, RB… removal rates, BA, BB… production rates, others used for coupling
FOR i = ja TO js; GOSUB olddecay
axt(1, i) = olddecaydiffA + ba * B/(1 + cc * c * c + sa * A * A + cb * d/(1 + cd * e * e))
axt(2, i) = olddecaydiffB + bb * C
axt(4, i) = olddecaydiffD + bd * (a/(1 + sd * d * d + sc * e) + se)
axt(5, i) = olddecaydiffE + be * a * a/(1 + ce * d* d)
DA = 0.0000; RA = 0.0100; BA = 0.0400; SA = 0.1000; CA = 0.0000; AA = 0.0000; GA = 0.0000;
DB = 0.1000; RB = 0.0060; BB = 0.0100; SB = 0.0000; CB = 33.0000; AB = 0.0000; GB = 0.0000;
DC = 0.0000; RC = 0.0000; BC = 0.0000; SC = 1.0000; CC = 10.0000; AC = 0.0000; GC = 0.0000;
DD = 0.0000; RD = 0.0001; BD = 0.0001; SD = 4.0000; CD = 15.0000; AD = 0.2000; GD = 0.0000;
DE = 0.0000; RE = 0.0005; BE = 0.0020; SE = 0.1000; CE = 0.5000; AE = 0.0000; GE = 0.0000;
(PC software for such simulations can be found on a CD included in the book: Hans Meinhardt ‘The Algorithmic Beauty of Sea Shells’
edition, 2003); Springer, Heidelberg, New York. See also http://www.eb.tuebingen.mpg.de/meinhardt