Abstract
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 the most commonly modeled 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.
1. 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–10] (which are initiated by the previously expressed pair-rule genes [11] at the beginning of gastrulation [12]) along the anterior-posterior 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–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.
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.
(a)
(b)
(c)
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.
2. 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–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).
(a)
(b)
(c)
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.
2.1. 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 , where and indicate the cell and cell-face locations, respectively (see Figure 3(a)). Finally, denotes on apposite cell face to ; that is, the subscripts and denote “neighbouring cell” and mod (7), respectively.
(a)
(b)
Steady state values of and are functions of ’s steady state value in the original model (see and ); 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 and ); 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
if the value of (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 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 time-evolution of is governed by (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 , 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 , , , , and , 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.
2.2. 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 , can act as the subset of a parameter set of numerical values that can be used to simulate the original model, say ; this is our first step in converting any . 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 . 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 time-evolution 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 assigned their values for using values assigned to the bold parameters (in both and ) 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. , , and , 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 , , , and , respectively. Hence parameters on the left hand side below are assigned their values for using values assigned to the bold parameters (in both and ) 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 (as outlined previously, the value assigned to in is equal to the value of in both and ).
Six final elements of , that is,
are assigned values in which result from random samplings over their individual ranges (see Table 1 for these ranges).
2.3. 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 and equal the sets of time-dependent variables used in our reduced model and the original model respectively. Then . Hence for any set of initial numerical values assigned to , say , we can define a set of initial numerical values that can be assigned to , say , such that . In all testing of our reduced model we choose such that it is a subset of the 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, . 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 and defined above. If, after 1000 minutes, the ON-OFF patterns of and in our model replicate wild-type patterns, is declared to be a solution of the reduced model.
Any such is converted into some for the original model using lifting. As described previously, 6 elements of result from random samples. These samplings follow the rules above (regarding parameters with(out) ranges covering multiple orders of magnitude).
Such 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 termed “Close to target pattern” and some , 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, 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 value, that is, if the value of its time-dependent variable is greater than 0.1. This ON/OFF threshold mirrors the approach crudely taken by von Dassow et al. (see supplement to [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 (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 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
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 ), 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 , the reduced model was ≈9.22 faster on average. For , the reduced model was ≈10.7 faster on average. For , 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.
2.4. 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 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 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 time-dependent 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 use the same values as 8 elements of . 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 referred to in the previous paragraph are assigned values that were randomly chosen (albeit these values were randomly chosen for 8 elements of originally).
3. Results and Discussion
The original model uses over three times as many time-dependent 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 and ) 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 and a realistic choice of ) 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 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 were found to be solutions at a frequency which peaked at 1 in 46 for the termed “Close to target patterns”. We achieved a similar result by testing our reduced model under a subset of these , an average of 1 in 50 randomly sampled 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 time-efficient 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].
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.
(a)
(b)
4. 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, high-dimensional models are less approachable and less accessible than lower-dimensional equivalents. Therefore, model-reducing 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.
Appendices
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 in a four-cell repeat and refers to the cell face. Here , , = amount of on the apposite cell face to , and cumulative amount of on cell faces apposite to the faces of cell .
where and . 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, and , take fixed values. von Dassow et al. wanted their model to be dimensionless, a change of variable was used where is time, 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 refers to the cell number in a four-cell repeat and 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: and × cumulative amount of on cell faces apposite to the faces of cell : where and .
C. Manipulations
Here follow the origins of our two identities, (1) and (2), which are derived from (A.2) and (A.7), respectively. At steady state (A.2) becomes
We can divide both sides of (C.1) by , because it is always nonzero ( equals 1 minute, takes values between 5 and 100 minutes), giving . Therefore , that is, we have obtained identity (1). To obtain identity (2), we consider (A.7) at steady state; that is,
This identity is satisfied for every cell face since each cell face in the model of von Dassow et al. has its own distinct . Since each cell face satisfies (C.2) at steady state, the following must be satisfied for each cell at steady state (as each cell has 6 such faces):
Noting that the various subscript notations used here are defined in Appendix A. We can divide both sides above by because it is always non-zero ( equals 1 minute, takes values between 5 and 100 minutes). Thus we obtain identity (2); that is,
D. Proof of Intrinsic EWG Symmetries in Original Model
The 24 ODEs which govern in Appendix A take the following form at steady state:
where both sides of (D.1) can be divided by ( is always non-zero). Hence a simple rearrangement gives
and, by dividing both sides of (D.2) by (which is always non-zero), can be made the subject:
Thus we have a matrix equation of the form
where and are both column vectors and is a matrix of coefficients. Left-multiplication of both sides of (D.4) by the inverse of (which exists) gives
Multiplying out the right hand side of (D.5) gives 24 linear equations of the form , where , , , and are functions of system parameters. Therefore any two which share the same , the same , the same and the same must be equal at steady state. Hence we used Maple to show that , , and 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 and . At steady state these 24 terms can, given the symmetry proven in Appendix D, be relabeled as follows: , , , , , , , , , , , and . Using this relabeling in the 24 ODEs which govern in Appendix A (see (A.5)) gives 12 distinct pairs of ODEs. Both members of each pair are identical. We give here a single member from each pair,
where , , , and .
Hence ((E.1)–(E.12)) describe the time-evolution of all the twenty four in the model of von Dassow et al. for a steady state analysis. Thus half of the 24 ODEs associated with in Appendix A are superfluous to a steady state analysis. As such, we use just the 12 ODEs above to monitor the time-evolution of in our reduced model, employing a second relabeling such that , , , , , , , , , , , and .
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 to 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.
Acknowledgments
The author would like to thank Professor Paul Glendinning for his invaluable contributions to this work as well as for his encouragement and kindness. The author is also grateful to the referees for their help in improving this paper.