Development and Analysis of a Quantitative Mathematical Model of Bistability in the Cross Repression System Between APT and SLBO Within the JAK/STAT Signaling Pathway

Cell migration is a key component in development, homeostasis, immune function, and pathology. It is important to understand the molecular activity that allows some cells to migrate. Drosophila melanogaster is a useful model system because its genes are largely conserved with humans and it is straightforward to study biologically. The well-conserved transcriptional regulator Signal Transducer and Activator of Transcription (STAT) promotes cell migration, but its signaling is modulated by downstream targets Apontic (APT) and Slow Border Cells (SLBO). Inhibition of STAT activity by APT and cross-repression of APT and SLBO determines whether an epithelial cell in the Drosophila egg chamber becomes motile or remains stationary. Through mathematical modeling and analysis, we examine how the interaction of STAT, APT, and SLBO creates bistability in the Janus Kinase (JAK)/STAT signaling pathway. In this paper, we update and analyze earlier models to represent mechanistically the processes of the JAK/STAT pathway. We utilize parameter, bifurcation, and phase portrait analyses, and make reductions to the system to produce a minimal three-variable quantitative model. We analyze the manifold between migratory and stationary steady states in this minimal model and show that when the initial conditions of our model are near this manifold, cell migration can be delayed.


INTRODUCTION
The acquisition of cellular migration plays a critical role in both normal and pathological development. A better understanding of the processes cells undergo as they transition from a stationary state to a migratory state is thus of broad interest. Studying the mechanics of how cohorts of cells move together introduces additional complexities, and existing models of collectively migrating cells differ greatly (Peercy and Starz-Gaiano, 2020;Stuelten et al., 2018;Aman and Piotrowski, 2010;Saadin and Starz-Gaiano, 2016;Olson and Nechiporuk, 2018;Leonard and Taneyhill, 2020;Friedl and Mayor, 2017). To study cell migration, some scientists turn to an experimental model system amenable to both genetic analysis and live imaging. In Drosophila melanogaster during oogenesis, a set of cells called border cells develop within a layer of follicle cells and later become migratory, leaving the nearby epithelial cells behind as they move to the oocyte (Montell et al., 2012;Saadin and Starz-Gaiano, 2016). Experimentalists have discovered much of the molecular regulation that governs this process and what causes some border cells to become motile while others remain stationary, including the primary biochemical and molecular pathways. We are interested in advancing a mathematical model of these pathways, which could have implications on acquisition of cell motility in animals in general.
The Janus Kinase/Signal Transducer and Activator of Transcription (JAK/STAT) signaling pathway has been shown by previous studies to be crucial in the motility of border cells, as well as in stem cells and immune response (Montell et al., 2012;Arbouzova and Zeidler, 2006;Trivedi and Starz-Gaiano, 2018;Amoyel and Bach, 2012;Amoyel et al., 2014). The JAK/STAT pathway is well-conserved from fruit flies to humans. Anterior polar cells in the Drosophila egg chamber (see Figures 1, 2A) secrete the cytokine Unpaired (UPD), which acts as the ligand for the transmembrane Domeless receptor in neighboring follicle cells. UPD is predicted to form a gradient across the adjacent cells (Van De Bor et al., 2011;Xi et al., 2003;Starz-Gaiano et al., 2008). The binding of UPD to Domeless activates JAK, leading to the phosphorylation of the activated JAK/receptor complex. The activated complex then recruits and phosphorylates STAT. The phosphorylated STAT dimerizes, moves to the nucleus, and acts as a transcription factor for specific target genes.
A key target of STAT is the C/EBP transcription factor Slow Border Cells (Slbo) (Montell et al., 1992;Starz-Gaiano et al., 2008). Cells specified by high levels of SLBO to become the border cells respond to chemoattractants that activate two receptor tyrosine kinases (RTKs), which turn on signaling cascades that promote directional movement (Montell et al., 2012;Stuelten et al., 2018;Saadin and Starz-Gaiano, 2016;Friedl and Mayor, 2017). Downstream of the RTKs, guided motility is governed in border cells largely through regulation of the actinomyosin cytoskeleton via the Rho GTPase, Rac, and myosin phosphorylation states, with coordinated changes in cell-cell adhesion mediated by E-cadherin (Stuelten et al., 2018;Saadin and Starz-Gaiano, 2016;Montell et al., 2012;Chen et al., 2020). The follicle cells receiving sufficiently high levels of UPDnot necessarily just those in closer proximity to the polar cells depending on the extracellular geometry (Manning et al., 2015)turn on a higher level of STAT activity and become motile border cells while nearby cells with lower levels do not. Interestingly, high STAT activity in ovarian follicle cell is sufficient to induce motility in usually stationary lateral follicle cells (Silver and Montell, 2001).
We focus on the protein products from two genes activated by STAT: Apontic (APT) and SLBO. SLBO promotes migratory behavior and an insufficient amount of SLBO prevents motility (Montell et al., 1992). APT protein is a transcription factor that downregulates the function of JAK/STAT and SLBO, and thus inhibits migration (Starz-Gaiano et al., 2008). APT acts as a feedback inhibitor on the JAK/STAT pathway, and this process is mediated by APT's activation of a microRNA that reduces STAT protein and activity (Yoon et al., 2011). APT also activates expression of Socs36E, which downregulates STAT signaling via a degradation pathway (Monahan and Starz-Gaiano, 2015). APT and SLBO exhibit cross-repressional behavior (Starz-Gaiano et al., 2008, 2009. APT directly represses slbo transcription while SLBO only decreases the level of expression of APT protein. The dominating protein in a given cell determines the cell fate: stationary or motile. This creates what appears to be a bistable system, as cells that receive intermediate amounts of STAT have the potential to reach either fate (Starz-Gaiano et al., 2009;Rorth, 1994). While it is reasonable based on the non-linearities in the system to model this cell fate regulation as bistability, in vivo experiments to support this are challenging and we are not aware of any particular experiments that have been done to confirm bistability. By identifying the conditions under which bistability occurs, modeling can help to design experimental protocols to confirm bistability. We base our model on the mechanistic model built by Ge et al. (2012). Focusing on the cross-repression system of APT and SLBO, they built a mathematical model using elementary interactions to identify which components of the system are sufficient for bistability. Depending on the strength of the UPD signal and thus the level of STAT activity, each border cell can become motile (SLBO dominates) or remain stationary (APT dominates). Cells with an intermediate level of STAT activity may fall above or below the threshold necessary for mobility.
Ge and Stonko created a 15-variable model including many mechanisms of STAT regulation (see Figure 2A) as well as the cross-repression system of APT and SLBO (see Figure 2B) with sufficient elements, specifically cooperativity in SLBO repressing apt mRNA translation, to cause bistability. We do not know the mechanism for this, but we suspect SLBO activates expression of a microRNA that mediates the effect. Ge et al. (2012) identified several miRNAs in the Drosophila genome that have upstream binding sites for SLBO activation and seed sequences that would target the apt mRNA 3 ′ untranslated regions.
To describe briefly the full model, in Equations (1) and (2) JAK (J) is altered to an enzymatically active state, J * , in the presence of UPD (U). Further in Equation (1) J * enzymatically activates two STAT monomers (S) to create the activated STAT dimer (S * 2 ) shown for STAT variables and complexes in Equations (3)-(6) with c 1 denoting the complex between activated JAK (J*) and STAT (S). Equation (5) represents APT sequestering S * 2 in c 2 . In Equations (7) and (8) APT (A) and SLBO (B) are produced from their mRNA and degraded, while A binding with S * 2 is accounted for. The mRNA production of stat (m σ ), apt (m α ), and slbo (m β ), is shown in Equations (7)-(9) based on a constant basal production and transcribing or probabilistic fraction based on non-transcribing stat (σ ), apt (α), and slbo (β) respectively, along with degradation. In Equation (9) B cooperatively enhances m α degradation, while in Equations (10) and (11) A enhances m σ and m β degradation. The gene state dynamics are shown in Equations (12)-(15), with S * 2 inducing transcription of σ , α, and β, while binding with A puts slbo into a repressed state (β R ). Model variables for Equations (1)-(15) from Ge et al. (2012) are in Table 1.   The 15-variable model is: In this paper we analyze and adapt the Ge and Stonko model so that a minimal reduction retains the dynamics of the detailed model. In section 2 we describe the methods we use to establish parameters and for bifurcation analysis. In section 3 we develop the minimal reduced model. In section 4 we analyze bifurcation results, identify the critical stable manifold separating the migratory and stationary steady states, validate experimental results, and compare our model to previous models. We apply this model in the interesting case of controlling microRNAmediated degradation of stat mRNA via APT and show that delays in STAT activation, even to the point of activation failure within a biophysical time span, are possible due to the proximity of the critical UPD level to a limit point bifurcation. We conclude with possible experimentation that could test and improve our model.

MATERIALS AND METHODS
In this section, we describe the methods used to analyze our minimal reduced model developed below. We use bifurcation and parameter analysis and identify the stable manifold between the stationary and motile steady states.

Time Course and Bifurcation Analysis
We use XPPAUT (Ermentrout, 2002) to create time course simulations and bifurcation diagrams and to analyze the bistability of the reduced model. We use the stiff integration method to solve our system of ODEs. A bifurcation occurs when a small change in parameter values results in a qualitative change in a system. Bifurcation analysis allows us to identify limit points treating UPD as a parameter, and helps to identify the separatrix between our steady states.

Establishing Parameters
In order to establish a more biophysically realistic model for this study, we researched existing literature to find data to establish parameter values. For some parameters we were able to find data specific to the JAK/STAT pathway or Drosophila. For other parameters, related pathways were used to obtain data relevant to this model. We were able to identify published values for general translation and transcription rates and applied the established rates to the lengths of JAK (encoded by hopscotch), STAT (encoded by Stat92E), APT, and SLBO genes and proteins (Hargrove et al., 1991;Lewin, 2004). The lengths are found in the fly genome database, Flybase (Thurmond et al., 2019). Protein and mRNA degradation rates have a wide range of average values so JAK, STAT, APT, and SLBO are assumed to conform to this range (Guido et al., 2006;Harris et al., 2011;Nicholson and Nicola, 2013). The protein to DNA binding and dissociation rates of JAK and STAT have been identified in general but not specifically for APT and SLBO (Halford and Marko, 2004;Parsaeian et al., 2013;Yang et al., 2002;Nicholson and Nicola, 2013).
Many signaling pathways, including the JAK/STAT pathway, are regulated by microRNAs (miRNAs) (Lui et al., 2015;Yoon et al., 2011). Ge et al. (2012) identified several microRNAs annotated in the genome (Thurmond et al., 2019;Betel et al., 2007) that have seed sequences corresponding to apt and slbo, suggesting they are also regulated by this mechanism. This regulation can occur through mRNA degradation, translational inhibition, or other means, making the parameters corresponding to miRNA kinetics difficult to assign. Ge and Stonko addressed this by condensing the various processes by which miRNA can affect mRNA into one degradation rate. We identified a parameter value for this combined effect through information from the model established in Yoon et al. (2011). The rate that slbo transitions in and out of its repressed state was also hard to identify due to lack of data, so we utilized the rate given in Ge and Stonko which was adapted from the rate for a different repressor in Harris et al. (2011) and information from Rorth (1994). The base levels of STAT and the total amount of JAK present in border cells were estimated from Yoon et al. (2011), Starz-Gaiano et al. (2008), and McGregor et al. (2002. Lastly, for the rates of STAT independent mRNA production, we again used the original parameter values from Starz-Gaiano et al. (2009). For slbo, this rate is most likely negligible. However, apt can be activated by means other than STAT. The protein Eyes Absent (EYA) can also activate transcription of apt (Starz-Gaiano et al., 2009).

Manifolds Separating Cell Fate Basins of Attraction
One goal in developing the three-variable model was to be able to fully understand the manifold that separates the steady states in the model. For a level of UPD in the bistable region, cells can either become motile or remain stationary depending on the initial conditions of the system.
We visually identified the manifolds by labeling initial conditions according to the steady state to which they converge. This allows us to see the basins of attraction for each steady state. These two stable steady states, one where SLBO dominates and the cell becomes motile and one where APT dominates and the cell remains stationary, are listed in the appendix with values for each variable. In our three-dimensional system we are able to approximate the manifold by fitting a surface to the boundary between the basins of attraction.
A graphical representation of the boundary manifold was achieved by creating a 3-d grid of initial conditions and determining to which steady state each converged within 500 min, a time deemed reasonable from experimental data. Any initial conditions that did not converge by this time were identified as lying near the manifold. We then used the MATLAB curve fitting toolbox to fit a surface to these points, creating an approximation of the manifold between the steady states.

DEVELOPING THE REDUCED MODEL
We began by reducing the fixed STAT cross-repression system of APT and SLBO (Equations, 7-10, 12-14) to a two-variable model. In the process of researching biophysically realistic parameters we discovered that the binding rate of STAT to target genes appears to be several orders of magnitude faster than any other process in the system, as seen in Table 2 (Halford and Marko, 2004;Parsaeian et al., 2013;Yang et al., 2002;Nicholson and Nicola, 2013;Karsten et al., 2006;Ekas et al., 2010). For example, k f α , k f β , and k f β R are at least three orders of magnitude faster than the translation and transcription kinetics. Additionally, α, β, and β R reach equilibrium significantly faster than the other variables. We used time-scale analysis to reduce the system. We made a quasi-steady state approximation for α, β, and β R and set those derivatives equal to zero. This allowed us to solve Equations (12), (13), and (14) for α * = 1 − α and β * = 1 − β − β R in terms of APT protein and STAT dimer: Our parameter values indicate that the mRNA processes occur at least twice as fast as the protein processes. This makes a quasi-steady state approximation for m α and m β plausible. Thus, we have a two-variable model where only APT and SLBO are dynamic.
The two-variable model is: After establishing the cross-repressional system of APT and SLBO, we reintroduced STAT dynamics to the model. Now we reduce the STAT activation system (Equations 1-6, 11, 15) through a number of assumptions. First, we ignored the theoretical APT-STAT complex (c 2 ) as its effects of APT sequestering STAT do not affect the steady state structure of the model, which we proved analytically. Then  Manning et al., 2015;Ghiglione et al., 2002;Wright et al., 2011;Ward et al., 1995;Hilton and Nicola, 1992 Dissociation rate of UPD to JAK k b UJ 0.1 nM· min −1 Ghiglione et al., 2002;Wright et al., 2011;Ward et al., 1995;Hilton and Nicola, 1992 Binding rate of J * complex to 2 STAT monomers  Harris et al., 2011;Rorth, 1994 we used the Michaelis-Menten approximation for the activated JAK (J * ) conversion of two STAT molecules to an activated STAT dimer (S * 2 ). This eliminates the JAK-STAT complex (c 1 ) and condenses the conversion. We assumed conservation of JAK to eliminate unactivated JAK (J) by defining a constant total JAK as J T = J + J * . We also assumed that UPD (U) activation of JAK and STAT activation are fast so J * and S * 2 can be solved for by quasi-steady state approximations. Lastly, Frontiers in Physiology | www.frontiersin.org similar to the reduction made for apt and slbo, the inactive stat gene state (σ ) and stat mRNA (m σ ) were solved by quasi-steady state approximations. These assumptions gave us the following additional approximations: Thus, producing our minimal three-variable model in APT, SLBO, and STAT: 4. RESULTS

Bifurcation Analysis
Time courses of STAT, APT, and SLBO show the difference between the motile and stationary steady states in Figures 3A,B. Figures 3C,D shows the comparable results for the original 15variable model. The two models converge to slightly different steady state values. This is due to the methods used to reduce the  STAT dynamics in the 3-variable model. Exact steady state values can be found in the Supplementary Materials. A bifurcation diagram of APT against UPD revealed a nontrivial state when U = 0 with a high level of APT, as seen in Figure 4A. This can be interpreted as the system being predisposed to the stationary cell fate until UPD and thus STAT is high enough, at UPD > 11.24 nM. As UPD approaches 0 on the upper branch, APT stops at a value of 44.73 nM. Since APT is also positively affected by STAT via UPD, very low UPD will decrease the steady state levels of APT but not to zero, leaving the stationary state the only steady state at very low UPD. Here the system encounters a limit point bifurcation and with a higher UPD level would transition to low APT and a motile cell fate.
Additionally, with the same parameters as Figure 4A, the STAT bifurcation ( Figure 4B) shows that if UPD begins at a high level in a cell, the cell will remain in the motile state even as UPD decreases to a very low level. This can be seen in experiments that dissociate polar cells and border cells. The border cells continue to migrate until their UPD level presumably drops below the threshold level, i.e., as they get too far away from the polar cells (Starz-Gaiano et al., 2008;Cai et al., 2014). This threshold also exists in the opposite direction-as the level of UPD increases in an epithelial cell, it becomes a migratory border cell once the threshold is reached (Manning et al., 2015).
The bistability in the model depends on two mechanisms. We show that the predominate non-linearity is in SLBO cooperative repression of apt mRNA translation. We also show the repressed state for the slbo gene induced by APT contributes to bistability (Starz-Gaiano et al., 2008). Bifurcation diagrams for STAT, APT, and SLBO with each combination of these criteria are presented as Figures S1-S3. As quadratic non-linearity in SLBO is made linear and the ability to reach the slbo gene repressed state is eliminated, the stationary basin of attraction becomes smaller until bistability is lost.

Parameter Values
Throughout the research on parameter values, the goal was to develop a range of realistic parameters to test in our model. There are two reasons why a range of values is desirable. First, from the biological point of view, many biological processes do not occur at a constant rate. Second, heterogeneity in cellular parameters will likely lead to some parameter variety. A range of average values is thus both more appropriate and more consistent with experimental data. We were able to test the robustness of the dynamics over the ranges of parameters to see if the model behavior matches experimentally observed outcomes.
The parameter values established through research from a variety of experimental systems and testing for bistability are shown in Table 2. The range of values identified for each parameter demonstrates a robust region of bistability. This adds confidence in the decision to use some parameter values that were established from a range of possible values. Figure 5 shows the percent change on a log 10 scale in each parameter that maintains bistability. The exact range for each parameter is listed in Table 3. The most sensitive parameter is δ Aσ . This makes sense biologically, as δ Aσ controls APT induced miRNA degradation of stat. Changing this rate directly affects whether a presumptive border cell produces enough STAT activity to pass the threshold needed for motility given a particular amount of UPD. APT has multiple levels of control on STAT, but this may be the most sensitive because there are multiple other regulators at the protein level. We explore the effects of this parameter on the outcome of the system below. The range of possible values for each parameter supports the robustness of the JAK/STAT pathway as a bistable system.

Manifolds in Three-Variable Model
We found that our three-variable minimal model retains the bistability and dynamics of the 15-variable model. Therefore, we use the three-dimensional manifolds in this model to better understand the behavior of the full model.
The stable manifold appears to be near-planar in the STAT-APT-SLBO phase space. The UPD value determines the position of the plane, with the manifold shifting from the STAT-SLBO plane to the STAT-APT plane as UPD increases, thus increasing STAT and SLBO production. Figure 6 shows how the manifold Any value in this range will maintain bistability in the model with other parameters held at their baseline values from Table 2. shifts as UPD increases from 0.0133 to 4, values that cover much of the bistable range of UPD.
To visualize the dynamics for each of the three values of UPD in the bistable region used in Figure 6, we plot trajectories to demonstrate how the manifold affects the outcomes of different initial conditions. Initial conditions that lie just below the manifold are attracted to the stable manifold near the unstable steady state and then repelled toward the stationary steady state. Initial conditions that lie just above the manifold behave similarly but converge to the motile steady state. Figure 7 depicts this behavior, with initial conditions that progress to the stationary steady state and the corresponding trajectories in blue and initial conditions that progress to the motile steady state and the corresponding trajectories in red. Three initial conditions and trajectories for each steady state are plotted. The motile steady state is in red, the stationary steady state is in blue, and the unstable steady state is in black. The unstable steady state lies on the 2D stable manifold.
When UPD is low, such as in Figure 7A, initial conditions in most of the phase space will result in the stationary steady state. When UPD is high, such as in Figure 7C, initial conditions in most of the phase space will result in the motile steady state. Figure 7D combines Figures 7A-C to show how the shifting manifold affects the system. As the value of UPD increases and the manifold moves below additional initial conditions, the motile basin of attraction expands and the motile and stationary steady states shift slightly. Movies showing a 3D rotation of each space and time series of each trajectory are available in the Figures S4-S8.

Delay From miRNA
The rate at which APT-induced miRNAs lead to degradation of stat mRNA is controlled in the model by the parameter δ Aσ . Figure 8A shows time courses of S in the three-variable model with UPD = 4 for different values of δ Aσ . We established a typical value of δ Aσ = 0.05 that allows S to equilibrate around 200 min, which is consistent with the time it takes for border cells to respond to UPD (Starz-Gaiano et al., 2008;Manning et al., 2015). The STAT level when δ Aσ = 0.05 is well above the threshold needed for the cell to become motile. An increase to δ Aσ = 0.17382 delays STAT convergence by a significant amount of time, about 600 min. If we increase δ Aσ to just 0.18, S never elevates within 1,000 min. Thus, the cell remains stationary. Each time course in Figure 8 represents a different outcome of the JAK/STAT pathway: GO, where SLBO dominates (δ Aσ = 0.05); STOP, where APT dominates (δ Aσ = 0.18); and SLOW, where the transition to motility is delayed (δ Aσ = 0.17382). The mechanisms of how APT-activated miRNAs affect STAT are still being analyzed, but initial findings support the idea that changes in miRNA activity can cause delays similar to those seen in our model (Yoon et al., 2011;Luo and Sehgal., 2012;Monahan and Starz-Gaiano, 2016;Sun et al., 2015). Delays in STAT activation and failure of activation are possible within a realistic time frame.
For basal levels of δ Aσ the system has normal STAT activation. However, as δ Aσ is increased the manifold separating the GO and STOP basins of attraction gets close to the initial condition, slowing the rise time to activation. If δ Aσ is raised enough the manifold crosses the initial condition and the trajectory is attracted to the STOP fate. Figure 8B further illustrates this delay. The 3D plot shows that at t = 300 the SLOW cell fate still has very low STAT and SLBO production. At the same time the GO and STOP cell fates have almost converged to their respective steady states. This change in outcome for one initial condition occurs because the increase in δ Aσ causes the manifold to shift just above the initial condition. This can be seen in Figure 9, where the initial condition lies between the SLOW manifold and the STOP manifold.

Experimental Tests
We tested our models to confirm if they reproduced certain behaviors identified in various experiments. It has been shown that if STAT activity is blocked by stage 9 of cell migration, the level of APT protein drops by about half (Starz-Gaiano et al., 2008, 2009Yoon et al., 2011;Monahan and Starz-Gaiano, 2015). It is also known that increasing the initial condition of APT should decrease the level of STAT protein and activity (Starz-Gaiano et al., 2008). These behaviors should be achievable by our models.
The three-variable model (Equations 29-31) reproduces the experimental behavior of APT when STAT is knocked down, as shown in Figure 10A. STAT initially converges to a value of about 4. Setting k S = 0 causes STAT to be constantly zero. This causes APT to converge to an equilibrium roughly half that of the stationary steady state. APT acts as an inhibitor to STAT activity, so for higher APT initial conditions we should see lower levels of STAT. This also occurs in the three-variable model ( Figure 10B).

Comparison to Previous Models
Previous mathematical models of the JAK/STAT pathway in Drosophila have focused on matching the qualitative behavior of the system. We compare our mechanistic model to the models developed in Starz-Gaiano et al. (2008) and Yoon et al. (2011) in Table 4. We rewrite our equations and use composite parameters to allow for easier comparison. The names of variables for STAT and SLBO differ across the models.
Our model and the earlier models share the same basic structure. The equations for STAT, APT, and SLBO all have a production term and a degradation term. The production terms include the cross-reactions of the proteins. However, our model describes these interactions in more detail, accounting for more of the known molecular interactions and biochemical production/degradation rates. For example, our STAT equation includes the feedback inhibition from APT. Our model also includes STAT independent production rates for both APT and SLBO, while the earlier models only include this for APT. By including the details of these molecular interactions, we gain the ability to mathematically examine the bistability of the JAK/STAT pathway.
Dynamic UPD is a feature of earlier models. Earlier models also include that levels of active STAT are positively regulated by production of UPD. Our model can be combined with the UPD dynamics developed in Manning et al. (2015). They used a partial differential equation to model the change in UPD concentration over time in three-dimensional extracellular space (Manning et al., 2015). This more detailed version of UPD dynamics includes the spatial element of how UPD diffuses from the polar cells into the border cells. It also helps to explain how border cells farther from the polar cells, which receive less UPD, can be delayed in becoming motile. Additionally, modeling The motile steady state is in red, the stationary steady state is in blue, and the unstable steady state is in black. t = 300 on each trajectory is marked in purple.
Frontiers in Physiology | www.frontiersin.org   in Peercy and Starz-Gaiano (2020) gives details of how the geometry of the egg chamber influences cell fate and migration dynamics. They also examine models of velocity and migratory cohort size as well as how clusters of border cells function to guide directional movement. Together our mechanistic model, the spatial dynamics of UPD from Manning et al. (2015), and models of collective migration from Peercy and Starz-Gaiano (2020) help us gain a better understanding of the threshold for motility in border cells.

DISCUSSION
We reduced the mechanistic Ge and Stonko model in stages to arrive at our minimal three-variable model. This final model allows us to better understand bistability in the JAK/STAT pathway. Our minimal model retains the parameters from the more complex model, allowing easier analysis but retaining critical properties. We discovered that non-linearity in SLBO repression of apt mRNA translation (Equation 9), and to a lesser extent APT repression of the slbo gene is required for bistability. The model displays the bistability of the stationary and motile cell states expected from experimental data for medium saturation of STAT activity. It has been shown that modulating UPD expression affects the number of migratory cells (Manning et al., 2015;Van De Bor et al., 2011;Silver and Montell, 2001;Xi et al., 2003;Grammont and Irvine, 2002). If UPD is low there are no migratory cells and if UPD is high there are many. This can induce normally stationary cells far away from the anterior to become motile (Manning et al., 2015). We established parameters that showed bistability was obtainable under realistic conditions. Every parameter in the 15variable model was found to have a wide range of values that guaranteed bistability. The robustness with respect to parameter values also suggests the biophysical utility of the model. Cell migration is an essential biological process, so the JAK/STAT pathway necessitates robustness under a range of parameter values, which the model supports. It makes sense that the JAK/STAT pathway would be able to operate successfully under a range of parameter values.
Ge and Stonko assumed that a key aspect necessary for bistability in APT and SLBO is cooperativity in SLBO repressing apt mRNA. In our dynamic STAT model, we found that the feedback inhibition of APT on STAT is necessary for bistability in STAT activity coming from the bistability between APT and SLBO. We also found that the initial amount of APT present in a cell is a major factor in whether or not the cell will become motile through affecting the level of STAT activity. Our minimal reduced model also reproduces behavior seen in experiments, such as declines in APT expression when STAT is knocked down and how initial levels of APT affect final levels of STAT (Starz-Gaiano et al., 2008).
We also showed that delays in STAT activation and failure of activation are possible within a realistic time frame. By controlling the degree of feedback inhibition of APT on STAT we can induce a delay in the transition to cell motility or cause the cell to remain stationary. This result is due to the proximity of the UPD level to a limit point bifurcation.
Assumptions made in the development of this model lead to some limitations. Specific parameter values for the exact reactions occurring in the egg chamber are not known, so we performed sensitivity analysis to show the robustness in parameter values. The STAT dynamics were reduced under a number of assumptions. One was the decision to ignore the c 2 equation, which allows APT to act as a buffer on STAT. Since the action of APT on STAT is known to function by preventing activation through miRNAs (Yoon et al., 2011) as well as limiting active protein present Starz-Gaiano, 2015, 2016), this assumption might oversimplify the larger system. Additionally, the quasi-steady state assumptions made in reducing STAT dynamics may also oversimplify the model.
Further study into the various methods by which APT inhibits STAT would enable us to improve how the model captures this interaction. Experiments that change the level of APT or STAT may give more information about how the delay in specification affects the migration of the cell cluster. Experimenting with apt mutants could tell us more about how the bistability we show in cell specification ultimately affects cluster migration and the further development of the egg chamber. When STAT is abnormally high more intermediate border cells trail behind the polar cells (Silver and Montell, 2001;Silver et al., 2005;Monahan and Starz-Gaiano, 2015;Starz-Gaiano et al., 2008). Socs36E limits STAT signaling and higher Socs36E reduces motile cell number. Both of these cases may suggest how a delay in cell specification affects cell migration.
It may also be important to consider how APT and SLBO may interact through miRNAs. Currently there is little data beyond the existence of these miRNAs (Ge et al., 2012). miRNA interactions directly affect bistability in the model, so greater detail could improve the model. Experiments that control the level of UPD secreted could identify how quickly the signal enters border cells and the levels of UPD which correspond with the activated STAT threshold to induce motility. To further gauge the STAT levels it would also be valuable to know the level of APT present in border cells through activation by EYA, prior to STAT activation and upregulation by STAT.
The JAK/STAT signaling pathway is known to be wellconserved. Specifically, it seems to be comparable in Drosophila and in humans (Trivedi and Starz-Gaiano, 2018;Arbouzova and Zeidler, 2006;Amoyel and Bach, 2012;Amoyel et al., 2014). Thus, as the model we have developed helps to explain cell motility in Drosophila, it may prove useful to our understanding of the process in humans. The "decision" for tumor cells to become motile is often a turning point for cancer progression. Additionally, STAT signaling is also well-known for controlling stem cell division decisions. Thus, the decisiveness of STATbased signaling may have a variety of roles in different cell types, but the model and its bistability shown in this paper may help to explain how this STAT signaling operates in different situations.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
BP and MS-G contributed conception and methodology, and reviewed and edited the manuscript. AB was the main writer the manuscript, and developed the reduced model, performed analysis, and ran simulations. All authors contributed to manuscript revision, read, and approved the submitted version.

FUNDING
This work was supported in part by National Science Foundation grant IOS-1656550 to MS-G. AB was funded in part by a summer research grant from the Dean of the College of Natural and Mathematical Sciences at UMBC.

ACKNOWLEDGMENTS
We would like to acknowledge the previous work done within the NSF-funded UBM program at UMBC that forms the basis of this research.