A Reduced Drosophila Model Whose Characteristic Behavior Scales Up

Computational biology seeks to integrate experimental data with predictive mathematical models—testing hypotheses which result from the former through simulations of the latter. Such models should ideally be approachable and accessible to the widest possible community,motivating independent studies. One of themost commonlymodeled biological systems involves a gene family critical to segmentation in Drosophila embryogenesis—the segment polarity network (SPN). In this paper, we reduce a celebrated mathematical model of the SPN to improve its accessibility; unlike its predecessor our reduction can be tested swiftly on a widely used platform. By reducing the original model we identify components which are unnecessary; that is, we begin to detect the core of the SPN—those mechanisms that are essentially responsible for its characteristic behavior. Hence characteristic behavior can scale up; we find that any solution of our model (defined as a set of conditions for which characteristic behavior is seen) can be converted into a solution of the original model. The original model is thus made more accessible for independent study through a more approachable reduction which maintains the robustness of its predecessor.


Introduction
It is widely acknowledged that mathematical models can provide valuable biological predictions [1,2].As levels of biological data increase, the scale of representative mathematical models must increase accordingly.This increase in scale will add computational complexity, posing a significant challenge to in silico testing.Hence, there is an ever increasing need for approaches which can facilitate the efficient study of large-scale models.A suitable approach is to reduce the dimension of large-scale models [3].Predictions of a reduced model will ideally translate to the original model [4].In this paper, we demonstrate such an approach for a model system in developmental biology.
Over the course of Drosophila embryogenesis (when the formation of embryonic cells and organs occurs [5]), a cascade of maternal and zygotic segmentation genes establishes the bodily layout [6].Normal development of segmental structures in the Drosophila embryo relies on a selective expression of the segment polarity (SP) genes [7][8][9][10] (which are initiated by the previously expressed pair-rule genes [11] at the beginning of gastrulation [12]) along the anteriorposterior axis.Products of SP genes (e.g., mRNAs and proteins) form a complex network of interaction which gives cells their distinctive identities [13] and functions [14,15], refines pattern formation [16,17] (in the form of stripes [18]), and maintains boundaries [19] between parasegments (repeated units of genetic expression [20,21] which lay the foundation for later adult structures [22]).The extensive interaction of SP gene products regulates their synthesis, and the resulting regulatory network the segment polarity network (SPN) has seen much interest from modelers; for example, see [23][24][25][26][27][28][29][30].Amongst these, no model has been more influential than the work of von Dassow et al.
In [23], von Dassow et al. discuss their modeling of the SPN as a system of coupled nonlinear ordinary differential equations (ODEs).The model uses time-dependent variables to describe the concentrations of SP gene products, as well as parameters including affinities and rate constants (see Appendix A).Time-evolution of the SPN can be simulated by integrating the model's ODEs.

ISRN Computational Biology
The model can also be represented as a directed graph  = (, ), where  are the nodes representing gene products of the system and  are the edges describing directed interactions among components (see Figure 1(a)).Typically, the directed edges of  indicate the influence exerted by source nodes (nodes with outgoing edges [31]) over the concentrations of target nodes (nodes with incoming edges).This influence can be either positive (excitatory) or negative (inhibitory) and different types of edge are used accordingly.
Following a comprehensive study of their SPN model, von Dassow et al. reported a robustness in its characteristic behavior (a single-cell wide wingless () stripe confronts a single-cell wide stripe of engrailed () and hedgehog (ℎℎ) at parasegmental boundaries [32]) to variation of initial conditions (ICs) and parameter values.
von Dassow et al. developed and tested their SPN model (the ODE form of which is given here in Appendix A) with Ingeneue (a simulation program in Java code) [33], having encountered computational slowness when initially using Mathematica.Similarly, we found extensive simulations of their model to be time-consuming when using MATLAB.This computational slowness for popular platforms is a significant barrier to independent studies that could offer new insight into how the SPN achieves its characteristic behavior so robustly.We aimed to motivate such studies of this SPN model (and enable our own study) by facilitating its analysis on a widely-used platform.
von Dassow et al. found that characteristic behavior of the SPN arises more from the nature, rather than the strength, of interactions between its components.Hence a simpler model might prove to be equally representative [34] and would produce a reduction in computational time [35].Therefore the first challenge was to find an equally representative reduced model.We proposed that such a model should demonstrate the noted robustness of the original model (see Appendix A)-such robustness is vital to the proper functioning of biochemical networks [36] and the SPN is one of the most robust examples [37].
Such a model can facilitate a time-efficient study of the original model if solutions of the reduced model provide solutions of the original model (where we define solutions to be conditions (parameter values and ICs) under which the model can replicate the characteristic behavior of the SPN).Hence, the original model should be able to replicate the SPN's characteristic behavior under any conditions that produce characteristic behavior in the reduced model.In this piece, we present a reduction that satisfies the criteria above (see Figure 1(b) for a graphical representation and Appendix B for the corresponding ODE system).
The rest of this paper is organized as follows.Section 2 introduces our approach to model reduction, derives the reduced model, introduces a method which converts any solution of the reduced model into a solution of the original model, describes the simulations we performed in our testing, and addresses any shortcomings of our work.Section 3 presents and discusses our key results.Our conclusions are given in the final section.

Materials and Methods
There are many well-practiced approaches to model reduction.These are predominantly techniques based on timescales/sensitivities of the model or lumping methods [38].In recent years, reduction methodologies for biological regulatory networks have begun to appear more frequently (e.g., see [39][40][41][42] for techniques that can be applied to boolean models).The suitability of any reduction method depends on both the type of model and the objective of reducing; this is why so many different methods are used and why no method can be considered supreme [38].
The model under study can be reduced for our objective by eliminating selected time-dependent variables (and the ODEs which govern their time-evolution).A variable is only eliminated if that variable (a) does not affect the model's behavior or (b) does not affect the behavior of a revised version of the model.
Each of the model's time-dependent variables is represented by a node in Figure 1(a).Each variable's incoming interactions are described in the corresponding ODE and represented by the incoming edges of the corresponding node in .Hence, eliminating a variable corresponds to removing a node from  and eliminating a variable's ODE corresponds to removing a node's incoming edges from .The revised digraph will not therefore feature the node which represented the eliminated variable, nor the edge(s) which represented that variable's incoming interaction(s).
A variable with outgoing interactions does not satisfy our criteria for elimination (see (a) and (b) above).Hence, a node with outgoing edges is not removed (and the variable it represents is not eliminated) unless (c) the target nodes are also removed or (d) these edges can be rewired using the following approach which we call merging (a simple example of merging is given in Figure 2).Consider 2 nodes; say  and , which are connected by a single (, ) edge in some system where the values of  and  are, or can be treated as being, equal at steady state.Any outgoing edge of , say (, ), can be replaced by a new (, ) edge without affecting the system's steady state behavior.If all outgoing edges of  are replaced in this way, then  has no outgoing interactions and can be eliminated.

Model Reduction.
We use the following steady state properties of the original model to merge  and  to merge  and cellular  as part of our reduction.For the original model (given in Appendix A), at steady state (see Appendix C for derivation of these equations).Note that cellular  is denoted by  , = ∑ 6 =1  , , where  ∈ [1,4] and  ∈ [1,6] indicate the cell and cell-face locations, respectively (see Figure 3 The (, ) edge of (a) can be replaced by a new (, ) edge; that is, the value of  can depend on the value of  (in the same way the value of  did depend on the value of  in (a)) without affecting the system behavior.Values of , , and  are the same as in (a), but  now has no output and thus plays no role in the system's behavior.As such,  and its incoming interaction(s) can be removed (as shown in (c)) without affecting the system.We refer to this process as the merging of  and . the subscripts  and +3 denote "neighbouring cell" and +3 mod (7), respectively.
Steady state values of  and ℎℎ are functions of 's steady state value in the original model (see (A.8) and (A.11)); say () and (), respectively.This is reflected in the wiring of Figure 1(a), that is, by the outgoing edges of .
From (1) we can see that Therefore () and () could replace () and (), respectively, in the original model without affecting that system's steady state behavior.We make these replacements in our revised model.The time-evolution of  is governed by (4) in the original model and by (5) in our revised model: )) )) Similarly, the time-evolution of ℎℎ is governed by (6) in the original model and by (7) in our revised model: )) )) ] ℎℎ −ℎℎ  ) .
The digraph representation is rewired accordingly; the (, ) and (, ℎℎ) edges of Figure 1(a) are replaced by (, ) and (, ℎℎ) edges, respectively, in the revised model.Hence,  has no outgoing interactions in the revised model.Therefore,  cannot affect the behavior of our revised system [43].As such, we eliminate  as part of our reduction.
In the original model, cleaving of cubitus interruptus protein is a function of cellular  at steady state (see (A.9) and (A.10)); say ℎ( , ).This is reflected in the wiring of Figure 1(a), that is, by the outgoing edge of cellular  (an outgoing edge of the  node whose head connects to the edge from  to ).We shall refer to this edge by the notation (, (, )).From (2), we can see that  (9) in the original model and by (10) in our revised model: Similarly, the time-evolution of  is governed by (11) in the original model and by (12) in our revised model: The digraph representation is rewired accordingly; the (, (, )) edge of Figure 1(a) is replaced by an (, (, )) edge in the revised model.A subgraph of Figure 1(a) emerges from this rewiring.This subgraph is defined by the set of nodes  = {, , } and the set of directed edges  = {(, ), (, )} (see Figure 1(c)).No element of  (the set of nodes in our revised system which are not elements of ) is the target node of any interaction involving elements of  in our revised system.Therefore no element of  can influence the behavior of any element of  in our revised system.As such, we eliminate all elements of  as part of our reduction.
The only outgoing interaction of ℎℎ in the original system was represented by the (ℎℎ, ) edge of Figure 1(a).Having eliminated  as part of our reduction, ℎℎ has no outgoing interactions in our revised system.Therefore ℎℎ cannot influence the revised system's behavior.As such, we eliminate ℎℎ as part of our reduction.
Within (A.5), the ODEs which govern concentrations of  in the original model, there exists an intrinsic mathematical symmetry such that, in all cells (i.e., for all ), at steady state (see Appendix D for proof and Figure 3(a) for chosen numbering of cell faces here).Thus the six  nodes within each cell of the original model form three pairs where members of the same pair have equal values at steady state.As such, members of the same pair can be labeled identically in the ODE system (see Appendix E for details) for the purpose of a steady state analysis.This results in 12 pairs of ODEs associated with  where members of the same pair are identical.Thus one member of each pair is superfluous to a steady state analysis of the original model and is eliminated as part of our reduction.Hence the revised model has half as many  variables as the original model.The revised digraph has just three  nodes per cell (half as many as the original digraph).These are the only nodes which lie on the faces of cells in our revised model; therefore cells become triangular (see Figure 3(b)).
The reduced model which results from the steps described in this section can be represented graphically, as shown in Figure 1(b).By picturing  "below"  (as in Figure 1(b)) rather than "above"  (as in Figure 1(a)), we avoid an untidy cross-over of edges.Similarly, we adjust Figure 1(a)'s positioning of  for greater clarity in Figure 1(b).
The time-evolution of variables that were eliminated as part of our reduction was governed in the original model by (A.2), (A.7), (A.11), (A.12), and (A.13), respectively.These equations do not appear in Appendix B (the ODE system which governs the behavior of our reduced model) and 14 of the parameters they feature do not feature in any other ODE of the original model.As such, these 14 parameters do not feature in our reduced system.

Using Our Reduced Model to Study the Original Model.
The original SPN model (given in Appendix A) can be studied through our reduction (given in Appendix B) via a process we call lifting.Lifting converts any set of conditions (i.e., parameter values and ICs) that can be used in a simulation of our reduced model into a set of conditions that can be used in a simulation of the original model.Hence lifting assigns values to the parameters and time-dependent variables of the original model that do not feature in our reduction.Importantly, any solution of our reduction (which we define as a set of conditions for which wild-type ON-OFF patterns of both  and  are stably maintained in our model) can be converted into a solution of the original model (defined as a set of conditions for which wild-type ON-OFF patterns of , , and ℎℎ are stably maintained in the original model) through lifting.
Under the set of conditions which results from the conversion of a reduced model solution, the original model must therefore (i) be able to retain the stable wild-type ON-OFF patterns of  and  (that were seen before conversion in our reduced model) and (ii) be able to replicate stable wild-type ON-OFF patterns of ℎℎ (that were not seen before conversion in our reduced model; our reduced model does not feature ℎℎ).The following approach achieves this.
Let  and  equal the sets of parameters used in the reduced model and the original model, respectively.Then  ⊂ .Hence a parameter set of numerical values that can be used to simulate the reduced model, say R, can act as the subset of a parameter set of numerical values that can be used to simulate the original model, say P; this is our first step in converting any R.The second part of our conversion assigns numerical values to  \ , that is, the following 14 parameters which are used in the original model but not in our reduced model: and   .These parameters feature in the ODEs which govern the time-evolution of eliminated variables: ℎℎ, , ,  and .
An element of  \ , say , is assigned a numerical value which has been randomly chosen over its range unless there exists an "analogous" element; say  ∈ .In this case,  is assigned the numerical value assigned to  in P. We propose that  and  are "analogous" if they satisfy all of the following criteria: (1)  and  have the same range, (2)  and  have the same interpretation (e.g., "half-life" (inverse of degradation rate)), (3)  and  feature in the ODEs which govern the timeevolution of variables  and  respectively (where the value of  is driven by the value of  in the original model).
This approach was designed to facilitate the conversion of a reduced model solution into a solution of the original model.We shall discuss the shortcomings of our approach in a future section.Assigning values to the parameters which feature in ( 14), the ODE that governs time-evolution of ℎℎ in the original model, is evidently crucial in determining the ON-OFF patterns of ℎℎ.In [27], Ma et al. note that ℎℎ-expression depends entirely on  in the SPN.An equivalent observation for a steady state analysis of the original model is that "ℎℎexpression depends entirely on " since concentrations of  and  are equal at steady state (as (1) shows).The values we assign to the parameters which feature in (14) reflect this dependence.These parameters are assigned values taken by parameters which feature in (15), the ODE which governs the time-evolution of  (in both the original model and our reduced model).Consider the following: )) )) Parameters which feature in (14) are assigned the values taken by their analogous parameters in (15).Identifying analogous parameters is made intuitive by the shared structure of ( 14) and (15).Hence parameters on the left-hand side below are ISRN Computational Biology assigned their values for P using values assigned to the bold parameters (in both P and R) as follows: These are analogous pairs of parameters; both members of each pair (from top to bottom) represent "half-life" (inverse of degradation rate) of mRNA, cooperativity coefficient of upregulation of mRNA, cooperativity coefficient of downregulation of mRNA, potency coefficient of upregulation of mRNA, and potency coefficient of downregulation of mRNA (where the mRNA in all cases is ℎℎ for the left hand member and  for the right-hand member of each pair).Concentrations of , , and  are driven by those of , , and ℎℎ, respectively, in the original model.This is described in the right hand sides of the ODEs which govern , , and  (i.e.(A.2), (A.7), and (A.12), resp.) by the terms ( 17), (18), and (19), respectively.Consider the following: The values we assign to the three parameters which feature in (17), (18), and (19) (i.e.  ,   , and   ) reflect this driving influence.These parameters are assigned values taken by analogous parameters that feature in the ODEs which govern the time-evolution of ,  and ℎℎ, (A.1), (A.6), and (A.11), respectively.Hence parameters on the left hand side below are assigned their values for P using values assigned to the bold parameters (in both P and R) as follows: These are analogous pairs of parameters; both members of each pair represent the "half-lives" (inverses of degradation rates) of protein (left hand member) and mRNA (right hand member) forms of the same gene products (from top to bottom: engrailed, patched, and hedgehog).Note that the final pair in our list here says that the values of   and  ℎℎ are equal in P (as outlined previously, the value assigned to  ℎℎ in P is equal to the value of   in both P and R).
Six final elements of  \ , that is, are assigned values in P which result from random samplings over their individual ranges (see Table 1 for these ranges).

Simulations.
A focus on the system's most important behavior can improve our understanding of the system's design [44].The most important behavior of the celebrated SPN model under study here is its characteristic behavior.Hence, the primary objective of our work is to facilitate the discovery of solutions for the original model.Let  V and  V equal the sets of time-dependent variables used in our reduced model and the original model respectively.Then  V ⊂  V .Hence for any set of initial numerical values assigned to  V , say PV , we can define a set of initial numerical values that can be assigned to  V , say RV , such that RV ⊂ PV .In all testing of our reduced model we choose RV such that it is a subset of the PV termed "Close to target pattern" by von Dassow et al. in [23].These ICs yielded the highest frequency of solutions in the original model's testing (see [23]), making them ideally suited to our primary objective.Our choice of ICs is realistic; from the point at which SP genes are first expressed, they are stably expressed throughout a fruit fly's development [24].
Our approach to simulations largely mirrors the approach of von Dassow et al. for the original model.We used Matlab for our work.Each of our model's 34 parameters is assigned a value which results from a random sampling over its individual range (we use the ranges defined by von Dassow et al. which are given in Table 1) to define some set, R. If a parameter's range covers multiple orders of magnitude, values are sampled uniformly over a logarithmic scale (to avoid any bias toward the upper end of their range).Values of other parameters are sampled uniformly over a linear scale.The reduced model is then simulated for 1000 minutes under conditions defined by the sets R and RV defined above.If, after 1000 minutes, the ON-OFF patterns of  and  in our model replicate wild-type patterns, R is declared to be a solution of the reduced model.
Any such R is converted into some P for the original model using lifting.As described previously, 6 elements of P\ R result from random samples.These samplings follow the rules above (regarding parameters with(out) ranges covering multiple orders of magnitude).
Such P are then shown to be solutions of the original model; we use the approach of von Dassow et al. here.When simulating their SPN model under the PV termed "Close to target pattern" and some P, von Dassow et al. used two temporal checkpoints: the first of which was 200 minutes.If, after 200 minutes, the ON-OFF patterns of , , and ℎℎ in their model replicated wild-type patterns, the simulation of the model continued to a second checkpoint, 1000 minutes (whilst ending otherwise).If, after 1000 minutes, the ON-OFF patterns of , , and ℎℎ still replicated wild-type patterns, P was declared to be a solution of the original model.
In assessing the ON-OFF patterns of the gene products mentioned above, we distinguish "ON" from "OFF" as follows.Each of these gene products is considered to be ON if its concentration is above one tenth of its maximal equilibrium  [23]).
Although our reduced model is not intended to be a competing model, we anticipate comparisons being drawn with the original model and, given its extensive discussion in [23], the frequency of solutions is of interest.For the ICs mentioned above, von Dassow et al. found 464 solutions in the testing of 21526 randomly chosen P (a frequency ≈ 1 in 46) (although Ingeneue uses a random search more than any other strategy, a directed search strategy is also used (see supplement to [23])).The original model has 48 parametersthus a parameter which is assigned a randomly chosen value has chance  of being compatible with the ON-OFF patterns captured in solutions where Although the robustness of our reduced model was not our primary interest, we endeavor to enable a comparison between the original model and our reduction.In four consecutive (unfiltered) parameter surveys of our reduced model (each consisting of one hundred randomly chosen R), we found 1, 3, 2, and 2 solutions, respectively, (all of which could be converted into solutions of the original model).
Our reduced model has just 34 parameters, thus a parameter which is assigned a randomly chosen value has chance  of being compatible with the ON-OFF patterns captured in solutions where, Note that these trial runs indicate that solutions of the original model can be found through parameter surveys of our reduced model with an (aggregate) average frequency of 1 in 50 (under the choice of ICs mentioned above).In reducing the SPN model of von Dassow et al. the frequency of solutions therefore remains essentially unchanged.Moreover, the interpretation of this frequency with respect to robustness also results in a similar value (note the comparative values of  and  above).Finally, we compared the computational time for  simulations of both models (see Appendix F for details of our tests).We found our reduced model to be significantly quicker than the original model, even for smaller .For  = 5, the reduced model was ≈9.22 faster on average.For  = 10, the reduced model was ≈10.7 faster on average.For  = 15, the reduced model was ≈12.3 faster on average.These trial runs indicate that the more extensive the simulations, the greater the benefit of using our reduced model.

Shortcomings of Our Work.
The original model is rewired so that cleaving of cubitus interruptus protein is driven by the mRNA, rather than the protein, form of the patched gene in our reduced model, that is, by  rather than cellular .This step of our reduction is a simplification of the SPN as it nullifies the     [] 0 ∑ 6 =1  ,+3  , term in (2) (a term which could otherwise exhibit non-zero values).
In defending our approach, we first highlight the lifting process.It is clear that our lifting process does not seek to restrict any of the key values here.Therefore, although the value of     [] 0 ∑ 6 =1  ,+3  , in ( 2) is restricted to zero in reducing the original model, its value is not restricted in the solutions of the original model which result from lifting.We conclude that treating  and cellular  as interchangeable in this way does not restrict the solutions we find for the original model and it has significant benefits from a model-reducing perspective.The rewiring of cleaving results in a system where , ,  and ℎℎ no longer play any role.Hence, this step allows us to eliminate 76 time-dependent variables (over one half of all timedependent variables in the original model) and 13 parameters (over one quarter of all parameters in the original model).
Moreover, although the reduced model simplifies the SPN in treating  and cellular  as equal, we find that solutions of our model exhibit great variation in parameter values.These solutions can be converted into solutions of the original model via lifting.Therefore, treating   and  , even though they were equal at steady state does not deprive the reduced model of the original model's robustness to parameter variation, prevent the use of our reduction as a tool for finding solutions of the original model.Hence we do not regret this part of our approach.
However, we partly regret an aspect of our lifting process; 8 elements of each P \ R use the same values as 8 elements of R.This is a restriction; each parameter was assigned its own randomly selected value in the work of von Dassow et al.Our approach was motivated by a need to find solutions of the original model as swiftly as possible; we required data to investigate the original model's behavior.The higher the success rate of our lifting process; our means of obtaining solutions for the original model, the more swiftly we could obtain data for analysis.Our intention was to maximize this success rate through our approach, no false positives have been found; that is, we have found no single solution of our reduced model that could not be converted into a solution of the original model through lifting.
We regret that our approach places this restriction on results.However, we wish to emphasize that the solutions we have found for the original model through our approach exhibit great variation in values of all associated parameters.Furthermore, we stress that the 8 elements of each P \ R referred to in the previous paragraph are assigned values that were randomly chosen (albeit these values were randomly chosen for 8 elements of R originally).

Results and Discussion
The original model uses over three times as many timedependent variables as our reduced model (the original and reduced models use 132 and 40 time-dependent variables, resp.).In addition, the original model and reduced model use 48 parameters and 34 parameters, respectively, thus the number of parameters is reduced by over a quarter.These properties combine to make parameter surveys of our reduction significantly swifter.
von Dassow et al. define a solution of their SPN model to be a set of conditions (i.e., some P and PV ) under which their model can sustain the wild-type ON-OFF patterns of , , and ℎℎ simultaneously (which is the characteristic behavior of the SPN).We have found that any set of conditions (i.e., some R and a realistic choice of RV ) under which our reduced model can sustain the wild-type ON-OFF patterns of  and  simultaneously can be converted into a solution of the kind defined by von Dassow et al. via the process we call lifting.(In rare cases, the randomly chosen values which are assigned to six elements of each P \ R as part of each lifting process have to be rechosen in order to make this conversion successfully.) The original model's robustness to parameter variation was celebrated by von Dassow et al. in [23].Randomly sampled P were found to be solutions at a frequency which peaked at 1 in 46 for the PV termed "Close to target patterns".We achieved a similar result by testing our reduced model under a subset of these PV , an average of 1 in 50 randomly sampled R could be converted into solutions of the original model via the process we call lifting.At the time of writing we have found over 50 such solutions in this way.
Computational Biology uses in silico experiments, simulating predictive mathematical models.Useful biological insights can be gained from computational modeling [45].Thus there is a clear incentive to improve the accessibility of such models so that independent studies (which may yield fresh insights) might be carried out.Such studies rely on approachable models that can be simulated in a timeefficient manner for analytic purposes.In this respect, the accessibility of the SPN model proposed by von Dassow et al. is compromised, this model cannot be simulated extensively in a time-efficient manner using popular platforms such as Mathematica [23] or Matlab.
This SPN model can now be analyzed efficiently.Previously, extensive simulations of the model were required to find solutions and such simulations were not time-efficient.However, our reduced model can be extensively simulated in a time-efficient manner and such simulations provide solutions of the original model.Thus, our work facilitates independent studies of the original model (noting that the use of Java-based Ingeneue is not consistent with an independent study).
Reducing biological systems to tractable forms can result in poor representations of these systems.Therefore, our reduction is no small achievement.Our reduced model can exhibit the SPN's characteristic behavior robustly and, moreover, any reduced model solution can be converted into a solution of the original model (via lifting).Thus the findings of our reduced model can scale up which demonstrates a clear link between the original model and our reduction.Techniques which allow large-scale models to be reduced without loss of their main properties are much needed [46].Figure 4: Merging nodes can be viewed as an abridging of the translation process.Structures of type (a), where ACT, REP, and  are proteins and  is some mRNA ( and  are products of the same gene), are seen throughout Figure 1(a).When  and  are products of the patched or engrailed genes we propose that such structures can be replaced by a structure of type (b), where  here is the result of having merged the  and  nodes of (a).
The SPN model proposed by von Dassow et al. is an ideal system for testing model-reducing techniques since its properties suggest that a simpler version of the system could prove to be equally effective.Our techniques produce such a model, our (much simpler) reduction can replicate the SPN's characteristic behavior robustly.An exciting prospect is that this may prove to be commonplace, that is, many regulatory networks which show robustness (a property seen throughout biological systems [47]) could be representatively modeled more simply and, we anticipate, using techniques that we have demonstrated here for the SPN.
The main techniques we used to reduce the original system were elimination and merging (these techniques could be used to reduce other models, for example, Ingolia's model in [25]).The effect of merging two nodes is to abridge the connection between them.Here, we merge products of the patched and engrailed genes.Thus, from a biological point of view, our reduction indicates that the translation process can be abridged for both the patched and engrailed genes when the SPN is mathematically modeled (see Figure 4).In the resulting model,  and  perform the roles that were performed by  and cellular , respectively, in the original model.Hence, merging preserves the underlying architecture of the original network and has brought us closer to the core of the SPN.

Conclusion
This paper proposes techniques for reducing the dimension of large-scale mathematical models.Such models are commonly used to describe complex biological systems.However, highdimensional models are less approachable and less accessible than lower-dimensional equivalents.Therefore, modelreducing techniques can motivate new studies of complex biological systems.Such studies are important; the reproducibility of previous findings can be tested and fresh insights may emerge.
Here, we reduce the celebrated and influential SPN model of von Dassow et al.Our reduction techniques result in a model which is significantly more approachable than its predecessor (under a third of the variables, three quarters of the parameters, and more time-efficient simulations).However, our model is not merely a simpler version of a preexisting model.Unlike many reduced models, our reduction provides a screen for solutions of the original model, that is, solutions of our reduced model supply solutions for the original model (via the process we call lifting).Hence, our model is not an alternative to its predecessor but an accompanying tool; it allowed us to find many solutions of the original model and discover keys to its celebrated robustness.Often the complexity of large biological systems obstructs such discoveries [48], but reducing the model results in a removal of inessential factors-we have now identified the factors which are vital to the model's robustness (which add to those reported by Ingolia in [25]) and hope to report these in a future publication.

A. SPN Model of von Dassow et al.
Segments of fruit flies are four cells wide before cell proliferation and wild-type ON-OFF patterns of SP gene products are repeated every four cells [49].Here follows the system of ODEs proposed by von Dassow et al. in the supplementary information for [23] with corrections (with thanks to George von Dassow for verifying typing errors made in their original publication).Cells are indexed by  ∈ [1,4] in a four-cell repeat and  ∈ [1,6] refers to the cell face.Here  , = ∑ 6 =1  , ,  , =  ,−1 +  ,+1 ,  ,+3 = amount of  on the apposite cell face to  , , and  , = cumulative amount of  on cell faces apposite to the faces of cell .
)) )).The interpretation and ranges of all parameters used here are given in Table 1.These details were found through a combination of the supplementary material for [23], the Java-encoded program Ingeneue created and used by von Dassow et al. and through private communications with George von Dassow (with our thanks).
We use the term "half-lives" here (as used by von Dassow et al., e.g., see [33]) but note that these are not half-lives in the true sense as they do not represent the time required for the quantity of material to fall by a factor of two.Their ranges are measured in minutes.
Two additional parameters,  0 and ]  , take fixed values.von Dassow et al. wanted their model to be dimensionless, a change of variable  =  0  was used where  is time,  0 is a characteristic time constant (equal to 1 minute) and  is a dimensionless variable (see supplement to [23]).The parameter representing the cooperativity of 's induction of , ]  , is also fixed at one,  is ever-preset ( is basally expressed [23]).
The value of each time-dependent variable in ((A.1)-(A.13))describes a fraction of its maximal steady state concentration (see supplement to [23]).All variables, apart from secondary forms, take cellular values between zero and one.Each of the model's secondary forms, ,  and -are scaled according to their primary forms (,  and , resp.).As such, their cellular concentrations may exceed one hundred per-cent, that is, values of  , ,   and  , may each exceed one (with thanks to George von Dassow for confirming this in private communications).

B. Our Reduced SPN Model
The ODEs associated to  in our reduced model can be found in Appendix E.Here follows the remaining ODEs associated with our reduced SPN model, where  ∈ [1,4] refers to the cell number in a four-cell repeat and  ∈ [1,3] refers to the cell-face.All parameters here use the ranges defined in Table 1.Table 2 defines the modified interpretations of some parameters.The modified shape of our cells (this model uses triangular cells where the original model uses hexagonal cells) results in the following modified notations:  , = 2 ∑ 3 =1  , and  , = 2 × cumulative amount of  on cell faces apposite to the faces of cell : )) )) )) )).

C. Manipulations
Here follow the origins of our two identities, (1) and ( 2 Noting that the various subscript notations used here are defined in Appendix A. We can divide both sides above by  0 /  because it is always non-zero ( 0 equals 1 minute, where we have used the numbering shown in Figure 3(a) here.

E. ODEs for EWG in Our Reduced Model
The network of von Dassow et al. uses 24 different  nodes, that is, one for each face of the four hexagonal cells in their repeated segments.These nodes are represented in Appendix A by terms of the form  , where  ∈ [1,4] and  ∈ [1,6].At steady state these 24 terms can, given the symmetry proven in Appendix D, be relabeled as follows:

F. Comparison of Computational Times for Original and Reduced Models
Both models were tested five times, each test comprised five, ten or fifteen simulations of both models (see Table 3).The computational time taken for each of these five tests (in seconds, to 4 decimal places) is reported in Table 3.All tests were performed on the same Windows XP computer (using Matlab).Each simulation ran a model from  = 0 to  = 1000 under the initial conditions described in Section 2.3.For each simulation, a model's parameters were assigned values that had been randomly sampled over their individual ranges.

Figure 1 :
Figure 1: The segment polarity network (SPN) as modeled (a) by von Dassow et al.[23,24] and (b) in our reduction: A single cell is shown together with one of its cell faces and the face of a neighboring cell.All cells contain an identical network.Nodes are elliptical for mRNAs, rectangular for proteins, and hexagonal for protein-protein complexes (whilst a rhombus is used for the dummy node).Triangular-headed edges represent positive influence and round-headed edges represent negative influence between connected nodes.(c) Subgraph of partially reduced Figure1(a).

Figure 2 :
Figure 2: Merging for simple 3-node system in steady state.(a) Value of  depends on value of  which depends on value of .The values of  and  are equal.(b)The (, ) edge of (a) can be replaced by a new (, ) edge; that is, the value of  can depend on the value of  (in the same way the value of  did depend on the value of  in (a)) without affecting the system behavior.Values of , , and  are the same as in (a), but  now has no output and thus plays no role in the system's behavior.As such,  and its incoming interaction(s) can be removed (as shown in (c)) without affecting the system.We refer to this process as the merging of  and .

Figure 3 :
Figure 3: Numbering for Appendices A and B models.Grid of four-cell wide, one-cell high segments together with our chosen numbering for (a) the original system and (b) our reduced system.
if the value of     [] 0 ∑ 6 =1  ,+3  , (see Appendix A for ranges of values here) is zero.Therefore ℎ(  ) could replace ℎ( , ) in the original model without affecting that system's steady state behavior whenever     [] 0 ∑ 6 =1  ,+3  , is zero.Hence we treat this value as zero to facilitate reduction, replacing ℎ( , ) with ℎ(  ) in our revised model (the shortcomings of this step are discussed in a future section).The timeevolution of  is governed by