# A ReducedDrosophila Model Whose Characteristic Behavior Scales Up

Authors: Andrew David Irving
Journal: ISRN Computational Biology
 (2013)

---

## 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 inDrosophila 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.


---

## Body

## 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 ofDrosophila 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 graphG=(V,E), where V are the nodes representing gene products of the system and E are the edges describing directed interactions among components (see Figure 1(a)). Typically, the directed edges of G 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.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 Figure 1(a).
(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 widewingless (wg) stripe confronts a single-cell wide stripe of engrailed (en) and hedgehog (hh) 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 AppendixA) 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 Figure1(b) for a graphical representation and Appendix B for the corresponding ODE system).The rest of this paper is organized as follows. Section2 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 Figure1(a). Each variable’s incoming interactions are described in the corresponding ODE and represented by the incoming edges of the corresponding node in G. Hence, eliminating a variable corresponds to removing a node from G and eliminating a variable’s ODE corresponds to removing a node’s incoming edges from G. 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 Figure2).Merging for simple 3-node system in steady state. (a) Value ofc depends on value of b which depends on value of a. The values of a and b are equal. (b) The (b,c) edge of (a) can be replaced by a new (a,c) edge; that is, the value of c can depend on the value of a (in the same way the value of c did depend on the value of b in (a)) without affecting the system behavior. Values of a, b, and c are the same as in (a), but b now has no output and thus plays no role in the system’s behavior. As such, b 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 a and b.
(a)
(b)
(c)Consider 2 nodes; sayx and y, which are connected by a single (x,y) edge in some system where the values of x and y are, or can be treated as being, equal at steady state. Any outgoing edge of y, say (y,z), can be replaced by a new (x,z) edge without affecting the system’s steady state behavior. If all outgoing edges of y are replaced in this way, then y has no outgoing interactions and can be eliminated.


### 2.1. Model Reduction

We use the following steady state properties of the original model to mergeen and EN to merge ptc and cellular PTC as part of our reduction. For the original model (given in Appendix A),
(1)eni=ENi,(2)PTCi,T=ptci-HPTCkPTCHH[HH]0∑j=16‍HHn,j+3PTCi,jat steady state (see AppendixC for derivation of these equations). Note that cellular PTC is denoted by PTCi,T=∑j=16‍PTCi,j, where i∈[1,4] and j∈[1,6] indicate the cell and cell-face locations, respectively (see Figure 3(a)). Finally, HHn,j+3 denotes HH on apposite cell face to PTCi,j; that is, the subscripts n and j+3 denote “neighbouring cell” and j+3 mod (7), respectively.Numbering for AppendicesA 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.
(a)
(b)Steady state values ofci and hh are functions of EN’s steady state value in the original model (see (A.8) and (A.11)); say f(EN) and g(EN), respectively. This is reflected in the wiring of Figure 1(a), that is, by the outgoing edges of EN. From (1) we can see that
(3)f(EN)=f(en),g(EN)=g(en).Thereforef(en) and g(en) could replace f(EN) and g(EN), 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 ci is governed by (4) in the original model and by (5) in our revised model:
(4)dciidτ=T0Hci×(Bi(κENciνENci/(κENciνENci+ENiνENci))νBciκBciνBci+Bi(κENciνENci/(κENciνENci+ENiνENci))νBci-cii),(5)dciidτ=T0Hci×(Bi(κENciνENci/(κENciνENci+eniνENci))νBciκBciνBci+Bi(κENciνENci/(κENciνENci+eniνENci))νBci-cii).Similarly, the time-evolution ofhh is governed by (6) in the original model and by (7) in our revised model:
(6)dhhidτ=T0Hhh(ENi(κCNhhνCNhh/(κCNhhνCNhh+CNiνCNhh))νENhhκENhhνENhh+ENi(κCNhhνCNhh/(κCNhhνCNhh+CNiνCNhh))νENhh-hhi),(7)dhhidτ=T0Hhh(eni(κCNhhνCNhh/(κCNhhνCNhh+CNiνCNhh))νENhhκENhhνENhh+eni(κCNhhνCNhh/(κCNhhνCNhh+CNiνCNhh))νENhh-hhi).The digraph representation is rewired accordingly; the(EN,ci) and (EN,hh) edges of Figure 1(a) are replaced by (en,ci) and (en,hh) edges, respectively, in the revised model. Hence, EN has no outgoing interactions in the revised model. Therefore, EN cannot affect the behavior of our revised system [43]. As such, we eliminate EN as part of our reduction.In the original model, cleaving ofcubitus interruptus protein is a function of cellular PTC at steady state (see (A.9) and (A.10)); say h(PTCi,T). This is reflected in the wiring of Figure 1(a), that is, by the outgoing edge of cellular PTC (an outgoing edge of the PTC node whose head connects to the edge from CI to CN). We shall refer to this edge by the notation (PTC,(CI,CN)). From (2), we can see that
(8)h(PTCi,T)=h(ptci)if the value ofHPTCkPTCHH[HH]0∑j=16‍HHn,j+3PTCi,j (see Appendix A for ranges of values here) is zero. Therefore h(ptci) could replace h(PTCi,T) in the original model without affecting that system’s steady state behavior whenever HPTCkPTCHH[HH]0∑j=16‍HHn,j+3PTCi,j is zero. Hence we treat this value as zero to facilitate reduction, replacing h(PTCi,T) with h(ptci) in our revised model (the shortcomings of this step are discussed in a future section). The time-evolution of CI is governed by (9) in the original model and by (10) in our revised model:
(9)dCIidτ=T0HCI(cii-CIi)-T0CCICIi(PTCi,TνPTCCIκPTCCIνPTCCI+PTCi,TνPTCCI),(10)dCIidτ=T0HCI(cii-CIi)-T0CCICIi(ptciνPTCCIκPTCCIνPTCCI+ptciνPTCCI).Similarly, the time-evolution ofCN is governed by (11) in the original model and by (12) in our revised model:
(11)dCNidτ=T0CCICIi(PTCi,TνPTCCIκPTCCIνPTCCI+PTCi,TνPTCCI)-T0CNiHCN,(12)dCNidτ=T0CCICIi(ptciνPTCCIκPTCCIνPTCCI+ptciνPTCCI)-T0CNiHCN.The digraph representation is rewired accordingly; the(PTC,(CI,CN)) edge of Figure 1(a) is replaced by an (ptc,(CI,CN)) edge in the revised model.A subgraph of Figure1(a) emerges from this rewiring. This subgraph is defined by the set of nodes S={PTC,PH,HH} and the set of directed edges C={(PTC,PH),(HH,PH)} (see Figure 1(c)). No element of S¯ (the set of nodes in our revised system which are not elements of S) is the target node of any interaction involving elements of S in our revised system. Therefore no element of S can influence the behavior of any element of S¯ in our revised system. As such, we eliminate all elements of S as part of our reduction.The only outgoing interaction ofhh in the original system was represented by the (hh,HH) edge of Figure 1(a). Having eliminated HH as part of our reduction, hh has no outgoing interactions in our revised system. Therefore hh cannot influence the revised system’s behavior. As such, we eliminate hh as part of our reduction.Within(A.5), the ODEs which govern concentrations of EWG in the original model, there exists an intrinsic mathematical symmetry such that, in all cells (i.e., for all i),
(13)EWGi,1=EWGi,4,EWGi,2=EWGi,3,EWGi,5=EWGi,6at steady state (see AppendixD for proof and Figure 3(a) for chosen numbering of cell faces here). Thus the six EWG 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 EWG 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 EWG variables as the original model. The revised digraph has just three EWG 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 Figure1(b). By picturing en “below” CN (as in Figure 1(b)) rather than “above” CN (as in Figure 1(a)), we avoid an untidy cross-over of edges. Similarly, we adjust Figure 1(a)’s positioning of ptc 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.


### 2.2. Using Our Reduced Model to Study the Original Model

The original SPN model (given in AppendixA) 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 en and wg 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 en, wg, and hh 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 ofen and wg (that were seen before conversion in our reduced model) and (ii) be able to replicate stable wild-type ON-OFF patterns of hh (that were not seen before conversion in our reduced model; our reduced model does not feature hh). The following approach achieves this.LetR and P equal the sets of parameters used in the reduced model and the original model, respectively. Then R⊂P. 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 P∖R, that is, the following 14 parameters which are used in the original model but not in our reduced model: HEN, Hhh, κENhh, κCNhh, νENhh, νCNhh, HHH, HPTC, rLMxferPTC, rLMxferHH, [HH]0, [PTC]0, kPTCHH, and HPH. These parameters feature in the ODEs which govern the time-evolution of eliminated variables: hh, HH, EN, PTC and PH.An element ofP∖R, say m, is assigned a numerical value which has been randomly chosen over its range unless there exists an “analogous” element; say q∈P. In this case, m is assigned the numerical value assigned to q in P~. We propose that m and q are “analogous” if they satisfy all of the following criteria:(1)
mand q have the same range,(2)
m and q have the same interpretation (e.g., “half-life” (inverse of degradation rate)),(3)
m and q feature in the ODEs which govern the time-evolution of variables M and Q respectively (where the value of M is driven by the value of Q 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 hh in the original model, is evidently crucial in determining the ON-OFF patterns of hh. In [27], Ma et al. note that hh-expression depends entirely on EN in the SPN. An equivalent observation for a steady state analysis of the original model is that “hh-expression depends entirely on en” since concentrations of EN and en 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 en (in both the original model and our reduced model). Consider the following:
(14)dhhidτ=T0Hhhdhhidτ=×(ENi(κCNhhνCNhh/(κCNhhνCNhh+CNiνCNhh))νENhhκENhhνENhh+ENi(κCNhhνCNhh/(κCNhhνCNhh+CNiνCNhh))νENhh-hhi),(15)denidτ=T0Hendenidτ=×(EWGn,T(κCNenνCNen/(κCNenνCNen+CNiνCNen))νWGenκWGenνWGen+EWGn,T(κCNenνCNen/(κCNenνCNen+CNiνCNen))νWGen-eni).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 P~ using values assigned to the bold parameters (in both P~ and R~) as follows:
(16)Hhh=Hen,νENhh=νWGen,νCNhh=νCNen,κENhh=κWGen,κCNhh=κCNen.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 ishh for the left hand member and en for the right-hand member of each pair).Concentrations ofEN, PTC, and HH are driven by those of en, ptc, and hh, respectively, in the original model. This is described in the right hand sides of the ODEs which govern EN, PTC, and HH (i.e. (A.2), (A.7), and (A.12), resp.) by the terms (17), (18), and (19), respectively. Consider the following:
(17)T0HEN(eni-ENi),(18)T0HPTC(ptci6-PTCi,j),(19)T0HHH(hhi6-HHi,j).The values we assign to the three parameters which feature in (17), (18), and (19) (i.e. HEN, HPTC, and HHH) reflect this driving influence. These parameters are assigned values taken by analogous parameters that feature in the ODEs which govern the time-evolution of en, ptc and hh, (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:
(20)HEN=Hen,HPTC=Hptc,HHH=Hen,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 HHH and Hhh are equal in P~ (as outlined previously, the value assigned to Hhh in P~ is equal to the value of Hen in both P~ and R~).Six final elements ofP∖R, that is,
(21)rLMxferPTC,rLMxferHH,[HH]0,[PTC]0,kPTCHH,HPHare assigned values inP~ which result from random samplings over their individual ranges (see Table 1 for these ranges).Table 1
Interpretation of parameters in SPN model of von Dassow et al.
Parameter
Range
Meaning
κ
W
G
e
n
1
0
-
3
–
1
Potency ofEWG’s induction of en
κ
C
N
e
n
1
0
-
3
–
1
Potency ofCN’s inhibition of en
κ
W
G
w
g
1
0
-
3
–
1
Potency ofIWG induction of wg
κ
C
I
w
g
1
0
-
3
–
1
Potency ofCI’s induction of wg
κ
C
N
w
g
1
0
-
3
–
1
Potency ofCN’s inhibition of wg
κ
C
I
p
t
c
1
0
-
3
–
1
Potency ofCI induction of ptc
κ
C
N
p
t
c
1
0
-
3
–
1
Potency ofCN’s inhibition of ptc
κ
B
c
i
1
0
-
3
–
1
Potency ofB’s induction of ci
κ
E
N
c
i
1
0
-
3
–
1
Potency ofEN’s inhibition of ci
κ
P
T
C
C
I
1
0
-
3
–
1
Potency ofPTC’s cleaving of CI
κ
E
N
h
h
1
0
-
3
–
1
Potency ofEN’s induction of hh
κ
C
N
h
h
1
0
-
3
–
1
Potency ofCN’s inhibition of hh
ν
W
G
e
n
1–10
Cooperativity ofEWG’s induction of en
ν
C
N
e
n
1–10
Cooperativity ofCN’s inhibition of en
ν
W
G
w
g
1–10
Cooperativity ofIWG’s induction of wg
ν
C
I
w
g
1–10
Cooperativity ofCI’s induction of wg
ν
C
N
w
g
1–10
Cooperativity ofCN’s inhibition of wg
ν
C
I
p
t
c
1–10
Cooperativity ofCI’s induction of ptc
ν
C
N
p
t
c
1–10
Cooperativity ofCN’s inhibition of ptc
ν
E
N
c
i
1–10
Cooperativity ofEN’s induction of ci
ν
P
T
C
C
I
1–10
Cooperativity ofPTC’s cleaving of CI
ν
E
N
h
h
1–10
Cooperativity ofEN’s induction of hh
ν
C
N
h
h
1–10
Cooperativity ofCN’s inhibition of hh
H
e
n
5–100
Half life ofen
H
E
N
5–100
Half life ofEN
H
w
g
5–100
Half life ofwg
H
I
W
G
5–100
Half life ofIWG
H
E
W
G
5–100
Half life ofEWG
H
p
t
c
5–100
Half life ofptc
H
P
T
C
5–100
Half life ofPTC
H
c
i
5–100
Half life ofci
H
C
I
5–100
Half life ofCI
H
C
N
5–100
Half life ofCN
H
h
h
5–100
Half life ofhh
H
H
H
5–100
Half life ofHH
H
P
H
5–100
Half life ofPH
r
E
n
d
o
W
G
1
0
-
3
–
1
Rate ofEWG endocytosis
r
E
x
o
W
G
1
0
-
3
–
1
Rate ofIWG exocytosis
r
M
x
f
e
r
W
G
1
0
-
3
–
1
Rate ofEWG cell-to-cell exchange
r
L
M
x
f
e
r
W
G
1
0
-
3
–
1
Rate ofEWG intramembrane diffusion
r
L
M
x
f
e
r
P
T
C
1
0
-
3
–
1
Rate ofPTC intramembrane diffusion
r
L
M
x
f
e
r
H
H
1
0
-
3
–
1
Rate ofHH intramembrane diffusion
α
C
I
w
g
1–10
Saturability ofCI’s wg-induction
α
W
G
w
g
1–10
Saturability ofIWG’s wg-induction
[
P
T
C
]
0
1
0
3
–
1
0
6
Max. cellular concentration ofPTC
[
H
H
]
0
1
0
3
–
1
0
6
Max. cellular concentration ofHH
k
P
T
C
H
H
1
0
-
5
–
1
0
-
3
Second-order binding rate
C
C
I
1
0
-
3
–
1
Rate at whichPTC cleaves CI


### 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.LetRv and Pv equal the sets of time-dependent variables used in our reduced model and the original model respectively. Then Rv⊂Pv. Hence for any set of initial numerical values assigned to Pv, say P~v, we can define a set of initial numerical values that can be assigned to Rv, say R~v, such that R~v⊂P~v. In all testing of our reduced model we choose R~v such that it is a subset of the P~v 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 Table1) 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 R~v defined above. If, after 1000 minutes, the ON-OFF patterns of en and wg in our model replicate wild-type patterns, R~ is declared to be a solution of the reduced model.Any suchR~ 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).SuchP~ 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 P~v 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 en, wg, and hh 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 en, wg, and hh 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 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 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 parameters—thus a parameter which is assigned a randomly chosen value has chance X of being compatible with the ON-OFF patterns captured in solutions where
(22)X48=46421526⟹X=(46421526)1/48⟹XX48=0.9232(to4  decimal  places).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 chosenR~), 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 Y of being compatible with the ON-OFF patterns captured in solutions where,
(23)Y34=8400⟹Y=(150)1/34⟹YY34=0.8913(to  4  decimal  places).Note that these trial runs indicate that solutions of theoriginal 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 X and Y above).Finally, we compared the computational time forn 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 n. For n=5, the reduced model was ≈9.22 faster on average. For n=10, the reduced model was ≈10.7 faster on average. For n=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.


### 2.4. Shortcomings of Our Work

The original model is rewired so that cleaving ofcubitus interruptus protein is driven by the mRNA, rather than the protein, form of the patched gene in our reduced model, that is, by ptc rather than cellular PTC. This step of our reduction is a simplification of the SPN as it nullifies the HPTCkPTCHH[HH]0∑j=16‍HHn,j+3PTCi,j 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 ofHPTCkPTCHH[HH]0∑j=16‍HHn,j+3PTCi,j 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 ptc and cellular PTC 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 PTC, PH, HH and hh 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 treatingptc and cellular PTC 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 ptci and PTCi,T 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 eachP~∖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 eachP~∖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).


## 2.1. Model Reduction

We use the following steady state properties of the original model to mergeen and EN to merge ptc and cellular PTC as part of our reduction. For the original model (given in Appendix A),
(1)eni=ENi,(2)PTCi,T=ptci-HPTCkPTCHH[HH]0∑j=16‍HHn,j+3PTCi,jat steady state (see AppendixC for derivation of these equations). Note that cellular PTC is denoted by PTCi,T=∑j=16‍PTCi,j, where i∈[1,4] and j∈[1,6] indicate the cell and cell-face locations, respectively (see Figure 3(a)). Finally, HHn,j+3 denotes HH on apposite cell face to PTCi,j; that is, the subscripts n and j+3 denote “neighbouring cell” and j+3 mod (7), respectively.Numbering for AppendicesA 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.
(a)
(b)Steady state values ofci and hh are functions of EN’s steady state value in the original model (see (A.8) and (A.11)); say f(EN) and g(EN), respectively. This is reflected in the wiring of Figure 1(a), that is, by the outgoing edges of EN. From (1) we can see that
(3)f(EN)=f(en),g(EN)=g(en).Thereforef(en) and g(en) could replace f(EN) and g(EN), 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 ci is governed by (4) in the original model and by (5) in our revised model:
(4)dciidτ=T0Hci×(Bi(κENciνENci/(κENciνENci+ENiνENci))νBciκBciνBci+Bi(κENciνENci/(κENciνENci+ENiνENci))νBci-cii),(5)dciidτ=T0Hci×(Bi(κENciνENci/(κENciνENci+eniνENci))νBciκBciνBci+Bi(κENciνENci/(κENciνENci+eniνENci))νBci-cii).Similarly, the time-evolution ofhh is governed by (6) in the original model and by (7) in our revised model:
(6)dhhidτ=T0Hhh(ENi(κCNhhνCNhh/(κCNhhνCNhh+CNiνCNhh))νENhhκENhhνENhh+ENi(κCNhhνCNhh/(κCNhhνCNhh+CNiνCNhh))νENhh-hhi),(7)dhhidτ=T0Hhh(eni(κCNhhνCNhh/(κCNhhνCNhh+CNiνCNhh))νENhhκENhhνENhh+eni(κCNhhνCNhh/(κCNhhνCNhh+CNiνCNhh))νENhh-hhi).The digraph representation is rewired accordingly; the(EN,ci) and (EN,hh) edges of Figure 1(a) are replaced by (en,ci) and (en,hh) edges, respectively, in the revised model. Hence, EN has no outgoing interactions in the revised model. Therefore, EN cannot affect the behavior of our revised system [43]. As such, we eliminate EN as part of our reduction.In the original model, cleaving ofcubitus interruptus protein is a function of cellular PTC at steady state (see (A.9) and (A.10)); say h(PTCi,T). This is reflected in the wiring of Figure 1(a), that is, by the outgoing edge of cellular PTC (an outgoing edge of the PTC node whose head connects to the edge from CI to CN). We shall refer to this edge by the notation (PTC,(CI,CN)). From (2), we can see that
(8)h(PTCi,T)=h(ptci)if the value ofHPTCkPTCHH[HH]0∑j=16‍HHn,j+3PTCi,j (see Appendix A for ranges of values here) is zero. Therefore h(ptci) could replace h(PTCi,T) in the original model without affecting that system’s steady state behavior whenever HPTCkPTCHH[HH]0∑j=16‍HHn,j+3PTCi,j is zero. Hence we treat this value as zero to facilitate reduction, replacing h(PTCi,T) with h(ptci) in our revised model (the shortcomings of this step are discussed in a future section). The time-evolution of CI is governed by (9) in the original model and by (10) in our revised model:
(9)dCIidτ=T0HCI(cii-CIi)-T0CCICIi(PTCi,TνPTCCIκPTCCIνPTCCI+PTCi,TνPTCCI),(10)dCIidτ=T0HCI(cii-CIi)-T0CCICIi(ptciνPTCCIκPTCCIνPTCCI+ptciνPTCCI).Similarly, the time-evolution ofCN is governed by (11) in the original model and by (12) in our revised model:
(11)dCNidτ=T0CCICIi(PTCi,TνPTCCIκPTCCIνPTCCI+PTCi,TνPTCCI)-T0CNiHCN,(12)dCNidτ=T0CCICIi(ptciνPTCCIκPTCCIνPTCCI+ptciνPTCCI)-T0CNiHCN.The digraph representation is rewired accordingly; the(PTC,(CI,CN)) edge of Figure 1(a) is replaced by an (ptc,(CI,CN)) edge in the revised model.A subgraph of Figure1(a) emerges from this rewiring. This subgraph is defined by the set of nodes S={PTC,PH,HH} and the set of directed edges C={(PTC,PH),(HH,PH)} (see Figure 1(c)). No element of S¯ (the set of nodes in our revised system which are not elements of S) is the target node of any interaction involving elements of S in our revised system. Therefore no element of S can influence the behavior of any element of S¯ in our revised system. As such, we eliminate all elements of S as part of our reduction.The only outgoing interaction ofhh in the original system was represented by the (hh,HH) edge of Figure 1(a). Having eliminated HH as part of our reduction, hh has no outgoing interactions in our revised system. Therefore hh cannot influence the revised system’s behavior. As such, we eliminate hh as part of our reduction.Within(A.5), the ODEs which govern concentrations of EWG in the original model, there exists an intrinsic mathematical symmetry such that, in all cells (i.e., for all i),
(13)EWGi,1=EWGi,4,EWGi,2=EWGi,3,EWGi,5=EWGi,6at steady state (see AppendixD for proof and Figure 3(a) for chosen numbering of cell faces here). Thus the six EWG 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 EWG 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 EWG variables as the original model. The revised digraph has just three EWG 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 Figure1(b). By picturing en “below” CN (as in Figure 1(b)) rather than “above” CN (as in Figure 1(a)), we avoid an untidy cross-over of edges. Similarly, we adjust Figure 1(a)’s positioning of ptc 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.


## 2.2. Using Our Reduced Model to Study the Original Model

The original SPN model (given in AppendixA) 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 en and wg 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 en, wg, and hh 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 ofen and wg (that were seen before conversion in our reduced model) and (ii) be able to replicate stable wild-type ON-OFF patterns of hh (that were not seen before conversion in our reduced model; our reduced model does not feature hh). The following approach achieves this.LetR and P equal the sets of parameters used in the reduced model and the original model, respectively. Then R⊂P. 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 P∖R, that is, the following 14 parameters which are used in the original model but not in our reduced model: HEN, Hhh, κENhh, κCNhh, νENhh, νCNhh, HHH, HPTC, rLMxferPTC, rLMxferHH, [HH]0, [PTC]0, kPTCHH, and HPH. These parameters feature in the ODEs which govern the time-evolution of eliminated variables: hh, HH, EN, PTC and PH.An element ofP∖R, say m, is assigned a numerical value which has been randomly chosen over its range unless there exists an “analogous” element; say q∈P. In this case, m is assigned the numerical value assigned to q in P~. We propose that m and q are “analogous” if they satisfy all of the following criteria:(1)
mand q have the same range,(2)
m and q have the same interpretation (e.g., “half-life” (inverse of degradation rate)),(3)
m and q feature in the ODEs which govern the time-evolution of variables M and Q respectively (where the value of M is driven by the value of Q 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 hh in the original model, is evidently crucial in determining the ON-OFF patterns of hh. In [27], Ma et al. note that hh-expression depends entirely on EN in the SPN. An equivalent observation for a steady state analysis of the original model is that “hh-expression depends entirely on en” since concentrations of EN and en 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 en (in both the original model and our reduced model). Consider the following:
(14)dhhidτ=T0Hhhdhhidτ=×(ENi(κCNhhνCNhh/(κCNhhνCNhh+CNiνCNhh))νENhhκENhhνENhh+ENi(κCNhhνCNhh/(κCNhhνCNhh+CNiνCNhh))νENhh-hhi),(15)denidτ=T0Hendenidτ=×(EWGn,T(κCNenνCNen/(κCNenνCNen+CNiνCNen))νWGenκWGenνWGen+EWGn,T(κCNenνCNen/(κCNenνCNen+CNiνCNen))νWGen-eni).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 P~ using values assigned to the bold parameters (in both P~ and R~) as follows:
(16)Hhh=Hen,νENhh=νWGen,νCNhh=νCNen,κENhh=κWGen,κCNhh=κCNen.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 ishh for the left hand member and en for the right-hand member of each pair).Concentrations ofEN, PTC, and HH are driven by those of en, ptc, and hh, respectively, in the original model. This is described in the right hand sides of the ODEs which govern EN, PTC, and HH (i.e. (A.2), (A.7), and (A.12), resp.) by the terms (17), (18), and (19), respectively. Consider the following:
(17)T0HEN(eni-ENi),(18)T0HPTC(ptci6-PTCi,j),(19)T0HHH(hhi6-HHi,j).The values we assign to the three parameters which feature in (17), (18), and (19) (i.e. HEN, HPTC, and HHH) reflect this driving influence. These parameters are assigned values taken by analogous parameters that feature in the ODEs which govern the time-evolution of en, ptc and hh, (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:
(20)HEN=Hen,HPTC=Hptc,HHH=Hen,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 HHH and Hhh are equal in P~ (as outlined previously, the value assigned to Hhh in P~ is equal to the value of Hen in both P~ and R~).Six final elements ofP∖R, that is,
(21)rLMxferPTC,rLMxferHH,[HH]0,[PTC]0,kPTCHH,HPHare assigned values inP~ which result from random samplings over their individual ranges (see Table 1 for these ranges).Table 1
Interpretation of parameters in SPN model of von Dassow et al.
Parameter
Range
Meaning
κ
W
G
e
n
1
0
-
3
–
1
Potency ofEWG’s induction of en
κ
C
N
e
n
1
0
-
3
–
1
Potency ofCN’s inhibition of en
κ
W
G
w
g
1
0
-
3
–
1
Potency ofIWG induction of wg
κ
C
I
w
g
1
0
-
3
–
1
Potency ofCI’s induction of wg
κ
C
N
w
g
1
0
-
3
–
1
Potency ofCN’s inhibition of wg
κ
C
I
p
t
c
1
0
-
3
–
1
Potency ofCI induction of ptc
κ
C
N
p
t
c
1
0
-
3
–
1
Potency ofCN’s inhibition of ptc
κ
B
c
i
1
0
-
3
–
1
Potency ofB’s induction of ci
κ
E
N
c
i
1
0
-
3
–
1
Potency ofEN’s inhibition of ci
κ
P
T
C
C
I
1
0
-
3
–
1
Potency ofPTC’s cleaving of CI
κ
E
N
h
h
1
0
-
3
–
1
Potency ofEN’s induction of hh
κ
C
N
h
h
1
0
-
3
–
1
Potency ofCN’s inhibition of hh
ν
W
G
e
n
1–10
Cooperativity ofEWG’s induction of en
ν
C
N
e
n
1–10
Cooperativity ofCN’s inhibition of en
ν
W
G
w
g
1–10
Cooperativity ofIWG’s induction of wg
ν
C
I
w
g
1–10
Cooperativity ofCI’s induction of wg
ν
C
N
w
g
1–10
Cooperativity ofCN’s inhibition of wg
ν
C
I
p
t
c
1–10
Cooperativity ofCI’s induction of ptc
ν
C
N
p
t
c
1–10
Cooperativity ofCN’s inhibition of ptc
ν
E
N
c
i
1–10
Cooperativity ofEN’s induction of ci
ν
P
T
C
C
I
1–10
Cooperativity ofPTC’s cleaving of CI
ν
E
N
h
h
1–10
Cooperativity ofEN’s induction of hh
ν
C
N
h
h
1–10
Cooperativity ofCN’s inhibition of hh
H
e
n
5–100
Half life ofen
H
E
N
5–100
Half life ofEN
H
w
g
5–100
Half life ofwg
H
I
W
G
5–100
Half life ofIWG
H
E
W
G
5–100
Half life ofEWG
H
p
t
c
5–100
Half life ofptc
H
P
T
C
5–100
Half life ofPTC
H
c
i
5–100
Half life ofci
H
C
I
5–100
Half life ofCI
H
C
N
5–100
Half life ofCN
H
h
h
5–100
Half life ofhh
H
H
H
5–100
Half life ofHH
H
P
H
5–100
Half life ofPH
r
E
n
d
o
W
G
1
0
-
3
–
1
Rate ofEWG endocytosis
r
E
x
o
W
G
1
0
-
3
–
1
Rate ofIWG exocytosis
r
M
x
f
e
r
W
G
1
0
-
3
–
1
Rate ofEWG cell-to-cell exchange
r
L
M
x
f
e
r
W
G
1
0
-
3
–
1
Rate ofEWG intramembrane diffusion
r
L
M
x
f
e
r
P
T
C
1
0
-
3
–
1
Rate ofPTC intramembrane diffusion
r
L
M
x
f
e
r
H
H
1
0
-
3
–
1
Rate ofHH intramembrane diffusion
α
C
I
w
g
1–10
Saturability ofCI’s wg-induction
α
W
G
w
g
1–10
Saturability ofIWG’s wg-induction
[
P
T
C
]
0
1
0
3
–
1
0
6
Max. cellular concentration ofPTC
[
H
H
]
0
1
0
3
–
1
0
6
Max. cellular concentration ofHH
k
P
T
C
H
H
1
0
-
5
–
1
0
-
3
Second-order binding rate
C
C
I
1
0
-
3
–
1
Rate at whichPTC cleaves CI


## 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.LetRv and Pv equal the sets of time-dependent variables used in our reduced model and the original model respectively. Then Rv⊂Pv. Hence for any set of initial numerical values assigned to Pv, say P~v, we can define a set of initial numerical values that can be assigned to Rv, say R~v, such that R~v⊂P~v. In all testing of our reduced model we choose R~v such that it is a subset of the P~v 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 Table1) 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 R~v defined above. If, after 1000 minutes, the ON-OFF patterns of en and wg in our model replicate wild-type patterns, R~ is declared to be a solution of the reduced model.Any suchR~ 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).SuchP~ 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 P~v 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 en, wg, and hh 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 en, wg, and hh 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 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 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 parameters—thus a parameter which is assigned a randomly chosen value has chance X of being compatible with the ON-OFF patterns captured in solutions where
(22)X48=46421526⟹X=(46421526)1/48⟹XX48=0.9232(to4  decimal  places).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 chosenR~), 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 Y of being compatible with the ON-OFF patterns captured in solutions where,
(23)Y34=8400⟹Y=(150)1/34⟹YY34=0.8913(to  4  decimal  places).Note that these trial runs indicate that solutions of theoriginal 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 X and Y above).Finally, we compared the computational time forn 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 n. For n=5, the reduced model was ≈9.22 faster on average. For n=10, the reduced model was ≈10.7 faster on average. For n=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.


## 2.4. Shortcomings of Our Work

The original model is rewired so that cleaving ofcubitus interruptus protein is driven by the mRNA, rather than the protein, form of the patched gene in our reduced model, that is, by ptc rather than cellular PTC. This step of our reduction is a simplification of the SPN as it nullifies the HPTCkPTCHH[HH]0∑j=16‍HHn,j+3PTCi,j 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 ofHPTCkPTCHH[HH]0∑j=16‍HHn,j+3PTCi,j 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 ptc and cellular PTC 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 PTC, PH, HH and hh 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 treatingptc and cellular PTC 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 ptci and PTCi,T 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 eachP~∖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 eachP~∖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).


## 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., someP~ and P~v) under which their model can sustain the wild-type ON-OFF patterns of en, wg, and hh 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 R~v) under which our reduced model can sustain the wild-type ON-OFF patterns of en and wg 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 P~v termed “Close to target patterns”. We achieved a similar result by testing our reduced model under a subset of these P~v, 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 usesin 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-basedIngeneue 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, en and ptc perform the roles that were performed by EN and cellular PTC, 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.Merging nodes can be viewed as an abridging of the translation process. Structures of type (a), where ACT, REP, andG are proteins and g is some mRNA (g and G are products of the same gene), are seen throughout Figure 1(a). When g and G are products of the patched or engrailed genes we propose that such structures can be replaced by a structure of type (b), where g here is the result of having merged the g and G nodes of (a).
(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.



---

*Source: 756829-2013-11-14.xml*