Constructing quantitative models from qualitative mutant phenotypes: preferences in selecting sensory organ precursors

(1)

Systems biology

Constructing quantitative models from qualitative mutant

phenotypes: preferences in selecting sensory organ precursors

Chao-Ping Hsu

1,2,

, Pei-Hsuan Lee

1

, Ching-Wei Chang

1

and Cheng-Tsung Lee

1

Institute of Chemistry, Academia Sinica, 128 Section 2 Academia Road, Nankang, Taipei 115 Taiwan and2Institute of Molecular Sciences, National Chiao-Tung University, 1001 Ta Hsueh Road, Hsinchu 300 Taiwan

Received on September 23, 2005; revised on March 1, 2006; accepted on March 3, 2006 Advance Access publication March 7, 2006

Associate Editor: Alfonso Valencia

ABSTRACT

Motivation: To study biology from the systems level, mathematical models that describe the time-evolution of the system offer useful insights. Quantitative information is required for constructing such models, but such information is rarely provided.

Results: We propose a scheme—based on random searches over a parameter space, according to criteria set by qualitative experimental observations—for inferring quantitative parameters from qualitative experimental results. We used five mutant constraints to construct genetic network models for sensory organ precursor formation in Drosophila development. Most of the models were capable of generating expression patterns for the geneEnhancer of split that were compatible with experimental observations for wild type and two Notch mutants. We further examined factors differentiating the neural fate among cells in a proneural cluster, and found two opposite driving forces that bias the choice between middle cells and the peripheral cells. Therefore, it is possible to build numerical models from mutant screening and to study mechanisms behind the complicated network.

Contact: [email protected]

Supplementary information: Supplementary data are available at Bioinformatics online.

INTRODUCTION

Living organisms can be viewed as complex biochemical systems that require enormous amounts of research regarding their essential components and interactions. Mathematical models and computer simulations are becoming increasingly useful tools in such studies (see, e.g. Amonlirdvimanet al., 2005; Barkai and Leibler, 1997; Bhallaet al., 2002; Chen et al., 2004; Hoffmann et al., 2002; Jaeger et al., 2004; von Dassow et al., 2000). In these studies, molecular production, degradation and interactions are modeled with chemical kinetics that are simulated using computers. This type of analysis allows discoveries for new intrinsic biological property and assists in experimental design.

Current computer technologies allow for large-scale dynamic biological network simulations, but to set up network models both qualitative and quantitative data are required. Qualitative information can be gathered from most molecular biological experiments, but quantitative data are much harder to come by. A project for taking large-scale metabolite measurements has

been established (Ishii et al., 2004), but the approach ignores differences among cells. When quantitative data are available, numerical algorithms can be used to find optimal parameters that minimize the distance between the model and the measured data (Zwolaket al., 2005). In a relatively small model, parameters can be estimated by trial-and-error. As an example, the dynamics and mutant phenotypes in circadian cycle dynamics is studied using this approach (Becker-Weimann et al., 2004). In constructing network models for dynamic simulations, estimating appropriate parameters is a challenge, making techniques for identifying para-meter values an important issue.

Network simulation without predetermined parameters

There are several approaches in network simulation that do not rely on predetermined quantitative parameters. One example is analyz-ing biochemical reaction fluxes in which known constraints are used in optimizing an output function (e.g. organism growth rate). This technique has been used for the genome-scale metabolism modeling of yeast (Famili et al., 2003) with genetic controls included as binary switches. For networks with complex genetic controls, find-ing all the thresholds of genetic switches that interface numerical fluxes and binary switches may not be easy. Another approach optimizes scores to obtain kinetic parameters, where a large penalty is added to the scores if the parameters cannot express mutant phenotypes (Amonlirdviman et al., 2005). This scheme requires a proper score function and a suitable optimization algorithm—both of which can directly affect computational complexity and numer-ical quality.

Computer-generated random searches are probably the easiest

approach to work without predetermined parameters. von

Dassowet al. (2000) used randomly generated parameter sets to study segment polarity networks. Parameter sets were accepted as ‘suitable’ whenever a network generated cell array patterns con-sidered compatible with those found in developing embryos. The likelihood of finding suitable parameters reflects system robustness, which Morohashi et al. (2002) believe serves as a measure of plausibility. Eldaret al. (2002) used randomly generated parameters to model bone morphogenic protein (BMP) gradient formation in Drosophila, where system properties were studied with a small number of parameter sets that were unusually robust in their gradient-forming capabilities.

The random screening approach is attractive because it allows us to take advantage of biological robustness and work with

To whom correspondence should be addressed.

at National Chiao Tung University Library on April 26, 2014

http://bioinformatics.oxfordjournals.org/

(2)

computer-generated parameters. However, there is always a chance of the correct output resulting from incorrect dynamics, thus making the establishment of search schemes that produce quality models very important. In this paper we address issues associated with screening strategies and subsequent analyses.

To uncover genetic control networks, scientists manipulate cer-tain genes and observe the resulting phenotypes. The appearance of mutant phenotypes not only indicates the presence of certain control pathways, but also implies that pathway strength exceeds a certain threshold in order for the effect to be observable. These clues can be used to build models. In this paper, we propose and test an approach in which the random generation scheme described above is used to screen for parameter sets that satisfy mutant conditions. Our goal is to demonstrate the possibility of obtaining such parameter sets and using them to study the system properties.

Lateral inhibition: Drosophila’s external sensory organ formation

Lateral inhibition is a process in which, through cell–cell inter-actions, a cell prevents adjacent cells from moving in the same developmental direction. In Drosophila, sensory organ precursor (SOP) formation is one well-understood example of lateral inhibi-tion via the Notch (N) protein signaling pathway. Mathematical modeling of the lateral inhibition process offers an ideal platform for determining how computational and theoretical studies can lead to increased understanding of system behavior and dynamics (Chang et al., 2003; Collier et al., 1996; Matsuno et al., 2003; Meiret al., 2002; Wearing et al., 2000).

The body of an adult fruit fly is covered with evenly spaced external sensory organs composed of four terminally differentiated cells derived from a single SOP. SOP formation involves expression of proneural proteins encoded byachaete (ac) and scute (sc) genes (Artavanis-Tsakonaset al., 1999; Campuzano and Modolell, 1992; Greenwald, 1998; Ghysen and Dambly-Claudiere, 1988). Initially, all cells in a proneural cluster are capable of forming a neural precursor, but only one cell actually forms an SOP. The N-Delta (Dl) signaling pathway allows for competition among neighboring cells. This has been shown by Hartenstein and Posakony (1990), Oellerset al. (1994) and Singson et al. (1994), and computationally demonstrated by Collieret al. (1996) and Meir et al. (2002).

Most pathways that affect SOP formation are able to amplify an inital bias among the proneural cells if one exists. However, it is unclear how a bias is formed at the beginning. Uniformly over-expressing the gene encoding signal-sending ligand, Dl, lead to a slightly affected phenotype (Seugnet et al., 1997; Doherty et al., 1997) and this result seems to contradict the N signaling model. Baker (2000) and Schweisguth (2004) have discussed several possibilities that may regulate N signaling activity and overcome the ectopic Dl expression by initiating the differentiation. While new bias-initiating pathways may be discovered in the future, it is important to understand the bias-initiating capability of the known N signaling pathway.

In this work, we first numerically built N signaling models using the mutant screening scheme. We demonstrated that mutant screen-ing has made a significant improvement on the quality of the models as shown by the high probability of generating correct expression patterns of Enhancer of split [E(spl)] proteins in wild type and twoN mutant models, and through several other analyses included in the

supplements. With the numerical models, we showed that a bias can arise from having different numbers and types of neighbors. We found that lateral inhibition prefers cells in contact with fewer signal-sending cells. On the other hand, cells in contact with epidermal cells may have a lower chance of becoming SOPs owing to an indirect mechanism involving the Dl’s autonom-ous regulation of N and the irreversible N–Dl interaction on cell surfaces. As a result, the chance of becoming an SOP is differen-tially affected by a cell’s position. While our deductions from numerical data are only hypotheses, we demonstrated that it is possible to construct numerical models from qualitative phenotypes, and such models offer a chance to study possible mechanisms hid-den in the complicated interactions and pathways.

METHODS

We modeled the mRNA and protein concentrations of each gene by integ-rating ordinary differential equations (ODEs) derived from a lateral inhibi-tion genetic network. Searching ranges were chosen so that corresponding rates, equilibrium constants or threshold concentrations were within typical reported ranges. In each trial, random values were generated to form a parameter set that was screened using wild type and mutant criteria. Para-meter sets that met all stated criteria were selected for further analysis; parameter sets that passed all but one mutant test were refined, with their values varied in a small range until the set in question met all mutant criteria (see Supplementary Tables S1 and S2 for details).

This parameter-inferring scheme was used to build aDrosophila SOP development network (Fig. 1A). At cell–cell contact regions, the interaction of N and Dl triggers a series of events, leading to the cleavage of N. The N intracellular domain, NICD, is then released as it enters the nucleus. Together with the Suppressor of Hairless [Su(H)] protein, NICDactivates the tran-scription of bHLH proteins encoded by theE(spl) gene complex. E(spl) proteins are repressors for the proneural genesac and sc. The proneural

Su(H) acorsc Ac orSc C H E(spl) E(spl) H NICD da Dl Da Dl N N C Dl N NICD Ac/Da or Sc/Da Neighbor H/Su(H) Su(H)/N C su(H) C A

A cell A neighbor cell

Cell membranes

1 2

4 3

B

Fig. 1. (A) Network model and (B) cell arrangement. In A, ovals represent mRNA, rectangles proteins and octagons protein complexes. ‘C’ denotes a constant supply of mRNA for ubiquitously expressed genes according to given parameters (independently varied for each gene). Arrows represent transcription activations (in this case, from a protein or protein complex to mRNA), translations (from mRNA to a protein), or complex formation (from several proteins to a single complex). Lines with circular dots denote repres-sion. In most cases a specific reaction is repressed, therefore dots were placed on corresponding arrows. The inhibitory effect of Su(H)/H onE(spl) tran-scription was modeled as a competitor with Su(H)/N. The seven cells in the middle of B represent a proneural cluster. The presumptive SOP cell is marked 1; the cell marked as 2 is one of the six adjacent cells. In a model with 6-fold symmetry, solving ODEs for 19 cells is the equivalent of solving ODEs for the cells marked 1–4.

(3)

genetic products Acheate-Daughterless (Ac/Da) and Scute-Daughterless (Sc/Da) complexes turn on the transcription ofDl, forming a mutual inhib-itory loop. Ac/Da and Sc/Da are able to turn on the transcriptions of pro-neural genesac and sc as well as their inhibitor E(spl), forming positive and negative feedback loops, respectively. The recently discovered ‘default repression’ effect of Hairless (H) onE(spl) transcriptions were included (Baroloet al., 2002) as a competition between NICD_{and H in binding}

with Su(H) and as a second competition to regulate theE(spl) gene tran-scription. Autonomous Dl regulation takes place during post-translational modifications of N (Schweisguth, 2004). In our model, the newly translated Dl inhibits N translation.

We note that our model is still far from complete, since there exist several other factors influencing the N signaling pathway. The influences of these factors, and the difficulties of including them in the current model, are discussed in the accompanying Supplement for interested readers.

NICDoverexpression leads to a loss of SOP phenotype—the result of N signaling (Rebayet al., 1993; Struhl et al., 1993). Temperature-sensitive N mutants (Nts) in non-permissive temperatures lead to a supernumerary SOP phenotype (Hartenstein and Posakony, 1990). TheH mutant also leads to a loss of SOP phenotype (Banget al., 1991). A simultaneous overexpression of genesSu(H) and H results in the production of supernumerary SOP cells (Morelet al., 2001), while an overexpression of the gene Su(H) only leads to a loss of SOP phenotype (Schweisguth and Posakony, 1994). Utilizing constraints derived from these mutants, it should be possible to obtain a balanced description of N, H and Su(H) activities.

Except when stated otherwise, each cell in our model had six adjacent cells and they were simulated in a group of 19 (Fig. 1B). In a 6-fold symmetrical configuration, only four cells need to be explicitly calculated. To account for the effect of a ‘positional cue’ from an unknown source (Koelzer and Klein, 2003), we used a mild prepattern for the initial amount of proneural gene products (mRNAsac and sc, and proteins Ac, Sc). The presumptive SOP cell in the center received twice as much proneural gene product as the six adjacent cells. All other cells received zero proneural genetic products, mimicking epidermal cells outside the cluster. Other details for model setting, parameter ranges and the ODEs used are presented in our Supplementary Material. Simulations were performed with a program written in Matlab run on either an Intel Pentium IV-based PC (Linux oper-ating system) or Power Mac G5 (OS 10.3).

RESULTS

Screening for mutants

Screening statistics for the parameter sets from the observed mutant conditions are presented in Table 1. Despite the low probability, we obtained two sets that met all five mutant criteria and 39 sets that met four out of five. After refining the latter parameter sets by varying their values within a narrow range, 29 of the 39 sets met all five criteria. We consequently obtained 31 parameter sets that satisfied all five mutant criteria.

E(spl) expression patterns

To determine whether the identified parameter sets actually medi-ated N signaling in SOP selection, we examined the expression patterns of the N signal’s direct target, the E(spl) gene, whose protein product represses the expression of theac and sc proneural genes. According to Jenningset al. (1995), E(spl) proteins can be detected in cells surrounding SOP cells, but not in SOPs. Further-more,E(spl) expression is greatly reduced in Notch mutants, while the ubiquitous expression of NICDleads to ectopic expression of E(spl)bHLH proteins. These results served as good tests for our models.

Among the 31 parameter sets obtained through mutant screening, 22 behaved similar to the wild-typeE(spl) gene expression—i.e. the amount of E(spl) protein in adjacent cells was at least double than that in the center cells. Furthermore, results from NICD

over-expression tests revealed that, in all of the 22 sets, the amount of E(spl) protein in every cell in the proneural cluster was uniformly higher than in any cell from the wild-type model. In N loss-of-function mutants, 20 of the 22 sets were capable of generating E(spl) expressions that did not exceed those in any cell in a wild-type cluster. These results, which are in agreement with observa-tions reported by Jenningset al. (1995), indicate that our mutant screening scheme is capable of identifying parameter sets that choose SOP cells by means of N signaling. In contrast, most of the parameter sets that only passed wild-type screening did not yield correctE(spl) expression patterns (Table 2).

Meiret al. (2002) used a second test with the same initial Ac and Sc protein concentrations across cells, without changing other initial concentrations. The test’s purpose was to avoid SOP generation via a simple threshold-cutting mechanism. In our model setting, the probability of passing this second test was about 30% for parameter

Table 1. Screening for parameter sets with mutant conditions

Genes Mutation Phenotypes Number of

parameter Sets (I) Screening

Number of parameter sets tested 2 852 053 Wild type Single SOP 15 259 Notch Loss of functiona _{All SOP}b ₂₂₃c

NICD _{Overexpression} _{No SOP}d _{4 398}

H Loss of function No SOPe _{2 194}

Su(H) Overexpression No SOPf ₆₀

Su(H) and H Overexpression All SOPg 443

All mutants and wild type 2

(II) Refinement

Number of parameter sets passing any four mutant tests

39 Number of parameter sets tested 120 866

All mutants and wild type 29

a_{Mutation tests performed for parameter sets that passed wild type single SOP test.} b_{Hartenstein and Posakony (1990).}

c_{Among those that passed wild-type screening.} d

Rebayet al. (1993); Struhl et al. (1993).

e_Bang_{et al. (1991).} f

Schweisguth and Posakony (1994).

g_Morel_{et al. (2001).}

Table 2. Number of parameter sets that produce expected E(spl) expressions in wild type and two different mutants

Tests performed

Wild type NICDOver N Loss Passed all mutant test 22/31a 22/22 20/22 Passed wild type test only 53/520 35/53 8/53

Passed a second testb 9/149 6/9 1/9

a_{Number of sets/total sets tested.} b_{As defined by Meit}_{et al. (2002).}

(4)

sets capable of generating the wild-type SOP phenotype—a much higher probability than for our mutant screening. However, for those parameter sets that passed this second test, we observed a very low probability of producing correct E(spl) expression patterns (Table 2).

Based on this requirement to generate a correct E(spl) pattern in wild type and two mutant settings, we ended up with 20 parameter sets which we studied further.

Factors influencing neural fates

Several pathways in the network are known to promote single SOP formation. The lateral inhibition pathway is the major component mediating intercellular signaling. The auto-activation of Ac/Da (which acts like a toggle switch) may play an auxiliary role (Becskei et al., 2001; Bhalla et al., 2002; Gardner et al., 2000). The regulatory effect of Dl on N production may provide extra stabilization for normal SOP formation. While it is difficult to elucidate the exact role of each process within the network by traditional means of genetic manipulation, in silico simulation can provide useful information.

Lateral inhibition prefers peripheral cells In a proneural cluster, the cells in the middle are in contact with more signal-sending cells than those on the periphery, implying that the middle cells may have less chance to become an SOP. By changing the number of con-tacting neighbors for the center cell, we tested for such a preference. In Figures 2A and B we depict new arrangements for a proneural cluster, where the center cell is in contact with four or five

neighbors. Time-course data for one representative parameter set are included in Figure 2.

In Figure 2, it is shown that reducing the number of contacting neighbors was accompanied by an increase in the proneural activity in the center cell (Ac/Da, Fig. 2C), and a decrease of inhibitory signals [NICDand E(spl), Fig. 2E and G]. For peripheral cells, the opposite changes were seen (Fig. 2D, F and H). The fewer signal-sending cells a cell is in contact with, the less NICDit receives, and the stronger Dl signal it is able to send. Therefore, with lateral inhibition, cells that contact with fewer signal-sending cells seem to be in a better position to become an SOP. Results for all 20 parameter sets are included in Figure S2 in the Supplementary material. All are qualitatively very similar to those in Figure 2.

To further explore if it is easier to build lateral inhibition models that select a peripheral cell rather than a center cell, we tested with a simplified model that contains only a ‘pure’ lateral inhibition path-way by removing all the branches and auxiliary pathpath-ways from the model (shown in Figure 3). We repeated the parameter selection scheme with this model to select an SOP in the middle (‘LI’) or on the periphery (‘LI-inv’). Since we have omitted the competitive repression that arises from H and Su(H), only two mutants are left for the screening tasks: N loss-of-function and NICD over-expression. In Table 3 we list the results.

Without auxiliary pathways, the chance of finding SOP forming parameter sets was very low (LI). In contrast, the chance increased by> 500-fold when we looked for the inverse final pattern—i.e. the peripheral cells became SOPs while the center cell was inhibited (LI-inv). Our data show that it is much easier to form a lateral inhibition module that selects a peripheral cell than one that selects a center cell.

This result is contrary to the fact that most SOPs arise from the central region of the clusters. Sometimes an SOP may arise from an eccentric position in a proneural cluster, but not at the edge of a cluster (Cubaset al., 1991). In the early development of a proneural cluster, the presumptive SOP does not seem to be very different from other cells in its shape, size or number of contacting neighbors (Cubaset al., 1991; Skeath and Carroll, 1991). Thus the preference towards cells in the periphery region from the mutually repressive lateral inhibition pathway may actually exist. This contradiction highlights the importance of other mechanisms that promote cells in the central region or suppress peripheral cells. In the

Fig. 2. (A) and (B) New arrangements for two proneural clusters tested. Cells marked 1, 2, 3 and 4 are equivalent to those similarly marked in Figure 1B, except that the hexagonal cell 1 is now a pentagon (A) or a square (B). (C–H) Time-dependent concentrations (in nM) of Ac/Da (C,D), NICD(E,F) and E(spl) (G,H) in the center cell (C,E and G) and the adjacent cells (D,F and H). Solid lines are from the original hexagonal model, while the dashed lines are from the cluster in A, and dotted lines are from the cluster in B.

Su(H) acorsc Ac orSc C E(spl) E(spl) NICD da Dl Da Dl N N C Dl N NICD Ac/Da or Sc/Da Neighbor Su(H)/N su(H) C

A cell A neighbor cell

Cell membranes C

Fig. 3. A simplified model containing only the N – E(spl) mediated lateral inhibition pathway. The gray dashed inhibition from Dl to N was removed in the LI and LI-inv tests, but included in LI+ test.

(5)

following, we discuss a mechanism for suppressing peripheral cells through existing pathways.

An auxiliary mechanism that suppresses peripheral cells One of the 20 parameter sets selected had an unusual capability of sup-pressing the neural fate of peripheral cells. When started with a reversed initial pattern (i.e. initial concentrations ofac (sc) and Ac (Sc) 10 times higher in the peripheral cells than in the center cell) the network still propagates to generate an SOP in the center. While the dynamics from this extreme case may not be realistic, a full analysis uncovers a possible mechanism that affects the neural fates. In Figure 4 we included the corresponding time course data.

Despite an initial 10-fold difference in the proneural gene pro-ducts, the difference inac mRNA and Ac/Da complex among the cells was quickly reduced (Fig. 4A and B). The higher concentration of proneural gene products in the peripheral cells (A and B, time regions marked in yellow shade), did not lead to a higher Dl protein concentration (C). Instead, the proteins NICD (D) and E(spl) (E) seemed to follow the proneural preferences set up by Dl, and they were followed byac and Ac/Da. The Dl translation of the two cells were similar (H), but the depletion due to N–Dl interaction showed a remarkable difference (I). At about 40–140 min the depletion in cell 2 was larger (more negative) than that in cell 1. (The red peak at 170 min is from the center cell sending massive signals to the adjacent cells—a green peak is also seen in NICD[D] at this time.) We further decomposed the depletion (J) and found that the deple-tion of Dl protein in cell 2 interacting with N protein on cells 3 (pink lines in J and K) and 4 (data not shown) is responsible for the higher depletion rate of Dl in the adjacent cell, suppressing its proneural activity. The ubiquitously expressed N is not regulated in the blank epidermal cells 3 and 4, since they do not express any proneural genes or theDl gene. However, in cells 1 and 2, N is regulated (F). Therefore, the abundant N protein on the surface of cells 3 and 4 offered a chance to decrease the amount of Dl protein in the adjacent cells (2), since N–Dl interaction leads to cleavage of N

and subsequent endocytosis. The mechanism is summarized in panel L of Figure 4.

In a separate test, we added the autonomous regulation of N to the simplified lateral inhibition model (‘LI+’ in Figure 3 with the gray inhibition indicated). Random screening results are included in Table 3 under the column ‘LI+’. The chance of finding normal SOP formation parameter sets with this model increased by a factor of 10, despite having two more parameters in the model. We further tested the screened parameter sets with a number of different initial concentrations. Among the 140 sets for the LI+ model, 45 could survive with a higher level of proneural gene products in the peri-phery cells, generating an SOP in the center at the end of the simulation. On the contrary, none of the 21 models obtained from the original LI model could survive when the adjacent cells received a larger amount of proneural gene products. (Full details and data are in the supplemental materials).

The Dl autonomous regulation of N can directly decrease the incoming NICD

signaling by regulating the amount of N on the cell membrane. Our results further indicated that this pathway is potentially able to bias the selection of an SOP via an indirect mechanism outlined in Figure 4L, such that cells in the periphery of a proneural cluster are less likely to become an SOP.

Effects of initial amounts of proneural genes

In lateral inhibition, a small initial difference leads to different cell fates. We tested to see whether the models possess this type of competitive nature, and whether the small initial bias used in the model setting is necessary.

We varied the initial concentrations of proneural genes on models derived from the mutant screening. In each test, the initial concen-tration of mRNAsac (sc) and proteins Ac (Sc) were varied inde-pendently using two scaling factors, one for the center cell and another for the adjacent cells. Two sets of results are shown in Figure 5. The diagonal line in each plot indicates cases in which the cells started with equal amount of proneural genes. The position of each dot in Figure 5 indicates the initial concentrations in the simulation, and the color specifies the phenotypes at the end of simulation.

In Figure 5A and B, results from two of the 20 mutant screened parameter sets are presented. The full results are included in Figure S4 A (Supplementary Material). A major characteristic of the results in Figure 5A and B, which is also seen in many datasets in Figure S4 A, is that the normal SOP foming region (filled by red dots) includes a boundary with a positive slope at its upper left side. Therefore, along this boundary, when the initial concentrations of the proneural genes are increased (or decreased) in all of the cells, the single SOP phenotype will still appear. This is a result of a competition mechanism, and it is in nice agreement with experi-mental observations (Chienet al., 1996), where overexpression of sc led to ectopic es organs that were solitary and evenly distributed. The details and shapes of the single SOP forming region vary for different parameter sets. In Figure 5B, the SOP region extends to the upper triangle, indicating that the model is able to suppress the neural fate of the adjacent cells even if the adjacent cells initially received higher amounts of proneural gene products. This result indicates that it is possible for the network to start from a uniform proneural gene distribution and evolve to a single SOP in a cluster. This parameter set has helped us to identify the auxiliary mechanism that suppresses neural fate in the peripheral cells discussed above.

Table 3. Random screening of a reduced model that contains only a lateral inhibition pathway (LI), with the same model but reversed SOP selection (LI-inv), and with a lateral inhibition plus an automonous regulation pathway (LI+)

Models LIa LI-invb LI+c

Number of parameter sets tested

3 585 062 2 036 429 2 467 460 Wild Type/NICD

Over/N Loss

21/21/21 2636/2635/2607 141/140/141

Passing all 21 2606 140

Number of parameters 36 36 38

Hit rates 5.86· 106 _1.28_{· 10}3 _5.67_{· 10}5

Inverse initial patternd 0/21 — 45/140

a

Model depicted in Figure 3, without the D1 to N autonomous regulation.

b_{With the same LI model but reversed wild type criteria: the parameter set was selected if}

it could generate SOPs at the peripheral site (#2 cell) and repress the sensory fate in the center cell. The initial conditions and the screening criteria remained unchanged.

c_{Model depicted in Figure 3, with the D1 to N autonomous regulation.}

d_{Starting from an inverse distribution of proneural gene products, the number of}

para-meter sets that are able to generate the normal phenotype is listed. The denominator is the total number of sets tested.

(6)

The proneural clusters are typically larger than 7 cells, as in the 19-cell model we simulated. Cells that are not in direct contact with the presumptive SOP were shown to be suppressed by the long range inhibition in which Dl is transported via filopodia (Joussineauet al., 2003; Renaud and Simpson, 2001). However, the indirect mechanism discussed above depends on the level of N on the cell surface, and therefore it is necessary to study a large cluster where the effect of proneural activity of cells that are not in direct contact with the center cell can be studied. In Figure 5 we included results of simulations extending to a 37-cell model, in which 19 cells in the middle form a proneural cluster and the 18 peripheral cells are blank cells forming an epidermal boundary. Extending to a larger cluster does not affect the major characteristic of the SOP-forming region in set I, but set II is sensitive to such a change. Results from most other sets (Fig.S4 E in the supplemental material) are largely the same as the original 19-cell model.

In set II, the tendency of becoming an SOP of the center cell is decreased when an extra layer of proneural cells are added. In the ‘no SOP’ region (pink dots, which were judged by proneural pro-ducts in cells 1 and 2), none of cells 3 and 4 could become an SOP. Even though parameter set II is not likely to be representative of a realistic situation, from this result we see that the auxiliary mech-anism suppresses peripheral cells which contact the epidermal cells. In experimental observations, some eccentric SOPs arose from cells that were about a cell’s distance from the boundary of the proneural clusters (Fig. 2 in Cubaset al., 1991), and such observations are consistent with our conclusion.

With the asymmetries in the numbers and types of neighbors, on which the cell–cell interaction strengths depend, it is possible to differentiate cells that are initially the same. In Table 4 we list statistics of searching for SOP-forming parameter sets starting with an equal amount of proneural genes. It is seen that, without an initial bias, it is quite unlikely to generate SOP-forming models from a pure lateral inhibition model (LI). With the auxiliary path-way added in the LI+ model, there exists a good chance to form models in which SOPs arise from the middle of a seven-cell cluster. Our results show that the pre-patterning used in simulation is not necessary, as long as there exist auxiliary pathways that help sup-press the neural fate of the peripheral cells.

A number of other numerical tests were performed to see the differences mutant screening made. The normal phenotype forming ranges of each parameter, robustness towards random parameter variation, and the complete results of varying initial amount of proneural gene products are included in the Supplementary Sections S2.1–S2.3.

DISCUSSION

Building models with mutant screenings

When constructing genetic networks, it is possible to impose a simple pattern formation condition and to obtain parameter sets using random search techniques. Requiring models to comply with the observed phenotypes of mutants significantly improves the chances of having them perform properly. For a 55-parameter network it is feasible to let a computer generate parameters accord-ing to constraints and ranges set by experimental results (Table 1). Accordingly, we observed that E(spl) expression patterns (a con-dition that was not included in the screening criteria) were mostly

Dl Dl N N C N NICD NICD Dl Dl N N C Dl N NICD Dl N C Ac/Da or Sc/Da Ac/Da or Sc/Da N C • Center cell: – N is regulated. – Initially, N on neighbors’ surface are regulated. – Dl is less depleted by N-Dl interaction. • Boundary cells:

– are epidermal cells. – N is not regulated.

There is abundant N on cell surface.

• Adjacent cells:

– N is regulated. – Some neighbors have

large amounts of N. – Dl depleted by N-Dl

interaction.

1 2 3

L

Fig. 4. Time-evolution in concentrations (A–F, in nM) or rates (G–J, in nM/min) calculated from the parameter set that can suppress the neural fate of the adjacent cells. Shown are the time-course data starting from a reversed pattern, in which the amount of proneural gene products given to the center cell was the same as in a standard test, but the adjacent cells received 10 times as much as the center cell. In panels A–I, red solid lines represent the con-centrations or production rates of the center cell, while green dashed lines are those of the adjacent cell. In F, concentration of N protein in the boundary cell #3 is included as a blue dotted line. The concentrations ofac (A), Ac/Da (B), Dl (C), NICD_{(D), E(spl) (E) and N (F) are shown. The yellow shaded regions}

in A–E indicate higher proneural activity in the adjacent cell (higher in ac, Ac/Da Dl and lower in NICD_{, E(spl)). Production rates for Dl are in G, which}

were further decomposed into translation (H) and depletion (I, J). The Dl depletion rates due to N-Dl complex formation are shown in I (for whole cells) and J (for each cell-cell contact). In J, the corresponding cell-cell contact regions are marked in the same color in the cluster arrangement in K. The solid red line in J, e.g. represents the depletion of Dl for cell 1 due to interaction with N on cell 2. L is a summary of the underlying mechanism. The dashed inhibitory curve from NICDto Ac/Da represents multiple steps of interaction in the model.

(7)

compatible with experimental results for the wild type and two mutants.

Mutant phenotypes carry pieces of quantitative information implicitly—in other words, the fluxes in pathways associated with perturbed genes must be above some unknown threshold for phenotypes to be observable. By requiring our network model to perform in a manner similar to observed mutant phenotypes, such implicit information is utilized in a simple and direct way. Since the literature contains a large amount of qualitative data, our scheme allows for the construction of numerical models to aid research in systems biology.

Since the mutants in our model share the same set of perturbed pathways and dynamics, the mutant screening events are not independent. Thus, the mutant screening task is not as time-consuming as one might expect. The probability of a parameter set that has passed the wild-type test meeting all mutant criteria is 2/15259’ 0.00013, which is three orders of magnitude larger than the probability of passing five independent mutant tests (4398· 223 · 2194 · 60 · 443/152595_{’ 6.91 · 10}8_{). The}

covari-ance matrix elements for the five mutant screening tests were all positive (data now shown), meaning that the odds of passing one test increases when another test is passed.

When it is necessary to add new elements, interactions or mutant constraints to a network, it is desirable to reduce the time and effort required for the model-building process. Our refinement strategy may be useful in this regard, since it includes a local screening of the region that surrounds the original parameter values. As shown in Table 1, whenever a parameter set met four out of five mutant criteria, a local search had a much higher probability of finding suitable parameters. We used Gaussian random searches (which are not strictly bound) for this task. Parameters that were considered more relevant to genes altered in the mutants were searched using larger ranges. This local search scheme allows for easy model extension and modification.

The random screening scheme allows for an even probe of para-meter space according to experimental observations. Despite the presence of uniform mutant phenotypes andE(spl) expression pat-terns, the 20 parameter sets we identified still differed in their detailed behaviors, reflecting the complexity involved in biological model construction (e.g. Figs S2–S4). The diversity of the models allows us to examine a large number of possibilities. We traced the source of differences among cells and identified two possible driv-ing forces that bias neural fates. In searchdriv-ing for possible outcomes from a complicated interaction network, we showed that the random screening scheme in conjunction with mutant constraints is very useful.

Neural fate determination

The N signaling pathway is seen in many tissues inDrosophila and other organisms. In different contexts, strength or timing require-ments may not be the same. Signaling networks need to survive these different cellular conditions and they may need to turn on different downstream targets. Utilizing functionally similar path-ways is a commonly seen strategy to ensure a proper outcome.

In selecting an SOP, all of the pathways studied are potentially able to amplify and stabilize a difference in the proneural activity, but little is known about their ability to create a difference. With our numerical analyses, it is shown that differences can be created within the N signaling network. With different numbers of signal-sending and non-signal-signal-sending neighbors, cell’s fates may be biased.

The role of the indirect mechanism cannot be confirmed until experimental evidence is presented. Since this pathway is heavily coupled with essential players in N signaling, experimental veri-fication may be rather indirect. One possibility is to use Notch mosaic clones, where regions lacking N are created. Careful stat-istics on the number of SOPs arising at different distances from the edge of mutant clones, when compared with similar statistics over normal tissue, can possibly provide supporting data.

CONCLUSION

In this paper we proposed and tested a random search scheme to obtain parameters from qualitative experimental results. We used the scheme to construct mathematical models for a developmental module ofDrosophila—the lateral inhibition network involved in SOP selection. We showed that most parameter sets identified in this manner can generate correctE(spl) expression patterns for wild type and twoN mutants. Further analyses revealed two opposite driving forces that bias the choice between middle cells and the peripheral cells in a proneural cluster. Our data show that lateral inhibition

Fig. 5. Phenotypes for varying initial concentrations of proneural genes in the center and adjacent cells. [red: single SOP (wild type); green: all SOP; pink: no SOP; yellow: inverse; blue: others.] Shown are the results of two parameter sets that passed all mutant tests. The initial concentrations used in mutant screening are at (1,0.5) in the scale used. A and B are results from a 19-cell model, C and D are from a 37-cell model, as indicated on the right. In C and D, red crosses are used to indicate cases where the central seven cells satisfy the single SOP condition, but with additional SOPs formed among cells 3 or 4.

Table 4. Random screening of LI and LI+ models, starting with equal amount of proneural gene products in all seven cells in the proneural cluster

Models LIa LI+b

Number of parameter sets tested 3 229 482 2 241 396 Wild Type/NICD_Over/_{N loss} _0/0/0 _23/23/23

Passing all 0 23

a_{Model depicted in Figure 3, without the D1 to N autonomous regulation.} b_{Model depicted in Figure 3, with the D1 to N autonomous regulation.}

(8)

favors cells that have contact with fewer inhibitory signal-sending neighbors. Through numerical analysis, we found an indirect mech-anism that arises from the autonomous regulation of N and the irreversible N-Dl interaction, which decreases the amount of Dl in the peripheral cells and thus suppresses neural fates in these cells. This new mechanism may contribute to a preferential selec-tion of center cells in proneural clusters. With these findings, we have demonstrated the potential of the computational approach in identifying working hypotheses and providing an analytical tool for complicated interactions.

ACKNOWLEDGEMENTS

The authors wish to acknowledge helpful discussions with

Drs Cheng-Ting Chien, (Institute of Molecular Biology,

Academia Sinica) and Haiwei Pi (Department of Life Sciences, Chang Geng University). The authors also acknowledge support from Academia Sinica and the National Science Council (Grant No. NSC 94-2627-M-001-004), Taiwan, ROC.

Conflict of Interest: none declared.

REFERENCES

Amonlirdviman,K.et al. (2005) Mathematical modeling of planar cell polarity to understand domineering nonautonomy.Science, 307, 423–426.

Artavanis-Tsakonas,S.et al. (1999) Notch signaling: cell fate control and signal integ-ration in development.Science, 284, 770–776.

Baker,N.E. (2000) Notch signaling in the nervous system. Pieces still missing from the puzzle.BioEssays, 22, 264–273.

Bang,A.G.et al. (1991) Hairless is required for the development of adult sensory organ precursor cells inDrosophila. Development, 111, 89–104.

Barkai,N. and Leibler,S. (1997) Robustness in simple biochemical networks.Nature, 387, 913–917.

Barolo,S.et al. (2002) Default represson and Notch signaling: Hairless acts as an adaptor to recruit the corepressors Groucho and dCtBP to Suppressor of Hairless. Genes Dev., 16, 1964–1976.

Becker-Weimann,S.et al. (2004) Modeling feedback loops of the mammalian circadian oscillator.Biophys. J., 87, 3023–3034.

Becskei,A.et al. (2001) Positive feedback in eukaryotic gene networks: cell differen-tiation by graded to binary response conversion.EMBO J., 20, 2528–2535. Bhalla,U.S.et al. (2002) MAP kinase phosphatase as a locus of flexibility in a

mitogen-activated protein kinase signaling network.Science, 297, 1018–1023. Campuzano,S. and Modolell,J. (1992) Patterning of theDrosophila nerves system: the

achaete-scute gene complex.Trends Genet., 8, 202–207.

Chang,C.-W.et al. (2003) Network modeling of Drosophila external sensory organ precursor formation: the role of recently studied genes.J. Gene. Mole. Biol., 14, 243–251.

Chen,K.C.et al. (2004) Integrative analysis of cell cycle control in budding yeast. Mol. Biol. Cell, 15, 3841–3862.

Chien,C.-T.et al. (1996) Neuronal type information encoded in the basic-helix-loop-helix domain of proneural genes.Proc. Natl Acad. Sci. USA, 93, 13239–13244. Collier,J.R.et al. (1996) Pattern formation by lateral inhibition with feedback: a

mathematical model of Delta-Notch intecellular signalling.J. Theor. Biol., 183, 429–446.

Cubas,P.et al. (1991) Proneural clusters of achaete-scute expressioin and the genera-tion of sensory organs in theDrosophila imaginal wing disc. Genes Dev., 5, 996–1008.

Doherty,D.et al. (1997) The Drosophila neurogenic gene big brain, which encodes a membrane-associated protein, acts cell autonomously and can act synergistically with Notch and Delta.Development, 124, 3881–3893.

Eldar,A.et al. (2002) Robustness of the BMP morphogen gradient in Drosophila embronic patterning.Nature, 419, 304–308.

Famili,I.et al. (2003) Saccharomyces cerevisiae phenotypes can be predicted by using constraint-based analysis of a genome-scale reconstructed metabolic network. Proc. Natl. Acad. Sci. USA, 100, 13134–131398.

Gardner,T.S.et al. (2000) Construction of a genetic toggle switch Escherichia coli. Nature, 403, 339–342.

Ghysen,A. and Dambly-Claudiere,C. (1988) From DNA to form: the achaete-scute complex.Genes Dev., 7, 723–733.

Greenwald,I. (1998) LIN-12/Notch signaling: lessons from warms and flies.Genes Dev., 12, 1751–1762.

Hartenstein,V. and Posakony,J.W. (1990) A dual function of the Notch gene in Drosophila sensillum development. Dev. Biol., 142, 13–30.

Hoffmann,A.et al. (2002) The IkB-NF-kB signaling module: temporal control and selective gene activation.Science, 298, 1189–1190.

Ishii,N.et al. (2004) Toward large-scale modeling of the microbial cell for computer simulation.J. Biotechnol., 113, 281–294.

Jaeger,J.et al. (2004) Dynamic control of positional information in the early drosophila embryo.Nature, 430, 368.

Jennings,B.et al. (1995) Role of Notch and achaete-scute complex in the expression of Enhancer of split bHLH proteins.Development, 121, 3745–3752.

Joussineau,C.D.et al. (2003) Delta-promoted filopodia mediate long-range lateral inhibition inDrosophila. Nature, 426, 555–559.

Koelzer,S. and Klein,T. (2003) A Notch-independent function of Suppressor of Hairless during the development of the bristle sensory organ precursor cell of Drosophila. Development, 130, 1973–1988.

Matsuno,H. et al. (2003) Boundary formation by notch signaling in Drosophila multicellular systems: experimental observations and gene network modeling by genomic object net.Pac. Symp. Biocomput, 152–163.

Meir,E.et al. (2002) Robustness, flexibility, and the role of lateral inhibition in the neurogenic network.Curr. Biol., 12, 778–786.

Morel,V. et al. (2001) Transcriptional repression by Supressor of Hairless involves the binding of a Hairless-dCtBP complex inDrosophila. Curr. Biol., 11, 789–792.

Morohashi,M.et al. (2002) Robustness as a measure of plausibility in models of biochemical networks.J. Theor. Biol., 216, 19–30.

Oellers,N.et al. (1994) bHLH proteins encoded by the Enhancer of split complex of Drosophila negatively interfere with transcriptional activation mediated by pro-neural genes.Mol. Gen. Genet., 244, 465–473.

Rebay,I.et al. (1993) Specific truncations of Drosophila Notch define dominant activ-ated and dominant negative forms of the receptor.Cell, 74, 319–329. Renaud,O. and Simpson,P. (2001) scabrous modifies epithelial cell adhesion and

extends the range of lateral signalling during development of the spaced bristle pattern inDrosophila. Dev. Biol., 240, 361–376.

Schweisguth,F. (2004) Regulation of Notch signaling activity. Curr. Biol., 14, R129–R138.

Schweisguth,F. and Posakony,J.W. (1994) Antagonistic activities of Suppressor of Hairless and Hairless control alternative cell fates in theDrosophila adult epi-dermis.Development, 120, 1433–1441.

Seugnet,L.et al. (1997) Transcriptional regulation of Notch and Delta: requirement for neuroblast segregation inDrosophila. Development, 124, 2015–2025. Singson,A. et al. (1994) Direct downstream targets of proneural activators in the

imaginal disc include genes involved in lateral inhibitory signaling. Genes Dev., 8, 2058–2071.

Skeath,J.B. and Carroll,S.B. (1991) Regulation of acheate-scute gene expression and sensory organ pattern formation ini theDrosophila wing. Genes Dev., 5, 984–995.

Struhl,G.et al. (1993) Intrinsic activity of the Lin-12 and Notch intracellular domains in vivo. Cell, 74, 331–345.

von Dassow,G.et al. (2000) The segment polarity network is a robust developmental module.Nature, 406, 188–192.

Wearing,H.J.et al. (2000) Mathematical modelling of juxtacrine patterning. Bull. Math. Biol., 62, 293–320.

Zwolak,J.W.et al. (2005) Parameter estimation for a mathematical model of the cell cycle in frog eggs.J. Comp. Biol., 12, 48–63.