 Research
 Open Access
 Published:
Modeling and analysis of the DeltaNotch dependent boundary formation in the Drosophila large intestine
BMC Systems Biologyvolume 11, Article number: 80 (2017)
Abstract
Background
The boundary formation in the Drosophila large intestine is widely studied as an important biological problem. It has been shown that the DeltaNotch signaling pathway plays an essential role in the formation of boundary cells.
Results
In this paper, we propose a mathematical model for the DeltaNotch dependent boundary formation in the Drosophila large intestine in order to better interpret related experimental findings of this biological phenomenon. To achieve this, we not only perform stability analysis on the model from a theoretical point of view, but also perform numerical simulations to analyze the model with and without noises, the phenotype change with the change of Delta or Notch expression, and the perturbation influences of binding and inhibition parameters on the boundary formation.
Conclusions
By doing all these work, we can assure that our model can better interpret the biological findings related to the boundary formation in the Drosophila large intestine.
Background
The large intestine of Drosophila embryo is a middle and large region of the Drosophila hindgut and is subdivided into dorsal and ventral domains, between which a onecellwide strand of boundary cells forms bilaterally for wildtype embryos [1–4]. The large intestine is a multicellular system that involves a number of cells composed of three cell types, dorsal, ventral and boundary cells, organized in a singlelayered epithelial tube [5]. For such developmental patterning problems, different kinds of computational strategies [6] have been proposed, e.g., signaling gradients [7] and activatorinhibitor systems [8], and many computational techniques are adopted, e.g., ordinary/partial differential equation [9] and colored Petri nets [10]. See [11] for a review. However, for the boundary formation in the Drosophila large intestine, the mechanism has been widely explored in vivo, e.g., [1–4, 12], but rarely from the computational point of view, e.g., [13].
It has been shown that the DeltaNotch signaling pathway plays an essential role in the boundary formation of the Drosophila large intestine [1, 12]. In fact, the DeltaNotch pathway is considered as one of the six major signaling pathways in cells, which is active in developing embryos at different phases [14]. Both Notch and Delta proteins are transmembrane proteins, where Notch proteins act as receptors and Delta proteins as ligands. When Delta ligands in a cell bind to Notch receptors in neighboring cells, all the cells in a system may evolve and finally form different types of patterns [15–18].
In order to understand the mechanism of the DeltaNotch pathway, there have been some mathematical models proposed. For example, Collier et al. proposed a simple ordinary differential equation (ODE) model for the DeltaNotch signaling pathway, and discussed numerical simulation of multiplecell systems [19]. Boareto et al. devised a theoretical framework that includes a couple of ODEs to explore the effects of Jagged in cellfate determination [20]. Sprinzak et al. gave a model of a set of ODEs for describing mutual inactivation of Notch and Delta and used this model to illustrate how cisinteractions between Notch and Delta generate mutually exclusive signaling states [21]. Specifically, Matsuno et al. analyzed the mechanism of the Notchdependent boundary formation in the Drosophila large intestine and built a hybrid Petri net model to numerically explore how the DeltaNotch pathway affects the boundary formation in twodimensional space [13].
In this paper, we propose a mathematical model for the DeltaNotch dependent boundary formation in the Drosophila large intestine based on the work of [13], aiming at better interpreting related biological findings and further making predictions. Compared with the existing work in this area, our work has the following main contributions.
(1) We give a mathematical model for the DeltaNotch dependent boundary formation in the Drosophila large intestine, which can better interpret relevant biological findings, so far obtained in the lab. Moreover, this model can be numerically simulated efficiently. To make the model both mathematically and biologically sound, we do the following analysis work.
(2) We perform local stability analysis on the twocell model to make the model mathematically sound. The analysis confirms that the model would reach an equilibrium, which corresponds to that the system would result in a stable pattern in either normal or mutant conditions. The model is also stable when some conditions are satisfied, which means even if there are small disturbances of parameters (e.g., small environmental noises), the system would still converge to a stable state after a period of time. We use the local stability analysis result to determine appropriate parameter values that make the system stable.
(3) We perform numerical simulations with and without noises of Notch expression, which shows that both deterministic and random simulation results are consistent with experimental observations. We further analyze how phenotypes change with the change of Delta or Notch expression.
(4) We analyze the perturbation influences of binding and inhibition parameters on the boundary formation. As the boundary formation is affected by many environmental factors or noises, to make the model more realistic, we need to consider noises into the model by finding those parameters that have significant influences on the dynamics of the model when they vary. According to the analysis result, we add appropriate random noise items to key parameters of the model.
This paper is structured as follows. In the section of methods, we introduce relevant biological background of the DeltaNotch dependent boundary formation in the Drosophila large intestine, and present a mathematical model for this biological phenomenon. In the section of results and discussions, we give stability analysis, simulation analysis, and parameter perturbation analysis of the model. Finally, the conclusion is given.
Methods
In this section, we introduce relevant biological background and give a description of the mathematical model we develop, as well as the simulation and analysis methods which we use.
Boundary formation in the Drosophila large intestine
The large intestine of Drosophila embryo is a major middle region of the Drosophila hindgut and is subdivided into dorsal and ventral domains (see Fig. 1), each of which is characterized by different cell types, dorsal or ventral. Between these two domains, a onecellwide strand of boundary cells forms in wildtype embryos.
The DeltaNotch signaling pathway plays an essential role in the formation of boundary cells [1, 12]. The main processes [1] are shown in Fig. 2. When the Delta ligands on the ventral cell surface bind to the Notch receptors of a neighboring Deltanegative dorsal cell, a DeltaNotch signaling cascade is activated. First the site 2 cleavage of the Notch proteins occurs to generate a transmembrane form (N^{ΔE}), followed by the Presenilindependent site 3 cleavage, producing an active Notch intracellular fragment (N^{intra}). N^{intra} then activates the target genes, inducing boundary cell differentiation. Moreover, Delta autonomously blocks the Presenilinmediated site 3 cleavage and thus inhibits Notch signal transduction within Deltapositive ventral cells.
In this paper, we aim to interpret the publicated findings about the boundary cell patterning in the large intestine [1, 4, 13] (see Fig. 3) and hypothesize the following scenarios for our model.

Scenario 1. In a wildtype embryo, a onecellwide strand of boundary cells forms at the interface of dorsal and ventral domains. See Fig. 3a for an illustration.

Scenario 2. In an embryo where overexpression of Notch or Delta proteins happens, the number of boundary cells would change. If we fix overexpression of Notch proteins, the number of boundary cells would increase with the decrease of Delta background. See Fig. 3 b1 to b2 for an illustration.

Scenario 3. However, if we fix the background of Delta proteins, the number of boundary cells would increase with the increase of overexpression of Notch proteins. See Fig. 3 c1 to c2 for an illustration.

Scenario 4. During the boundary formation in the Drosophila large intestine, environmental factors such as temperature may heavily affect the boundary formation, thus resulting in different (random) phenotypes. Therefore, the model to be built should incorporate such random noises into some parameters of interest.
In what follows, we will construct and analyze our model step by step.
Mathematical model
Celltocell interaction mediated by the DeltaNotch signaling pathway plays an essential role in the development of a multicellular organism such as the Drosophila large intestine. In this work we propose a mathematical model for the core part of the DeltaNotch dependent boundary formation in the Drosophila large intestine (see Fig. 4 for a onecell model) based on the work of [13]. That is, we only consider two key species, Delta (D) and Notch, and Notch can be in the inactive state (N) or active state (A). Inactive Notch can be converted into the active state when coupled with Delta from neighboring cells. Active Notch can inhibit the production of Delta of the same cell through target genes. Besides, Delta also inhibits the production of Notch of the same cell.
We arrange all cells in a regular M×N lattice (see Fig. 5 for an illustration), each site (cell) being a hexagon and having at most six neighboring cells. That is, we have N C=M×N/2 cells for an M×N lattice.
The mathematical model is then given as follows:
Here, 1≤i≤N C. D _{ i }, N _{ i } and A _{ i } denote the concentration of Delta proteins, inactive and active Notch proteins, respectively, in the i ^{th} cell. In \(\frac {\lambda }{1+\Delta \cdot {{A}_{i}}}\), λ represents the production rate of Delta proteins, and Δ the inhibition coefficient of active Notch proteins A _{ i }. d _{1} represents the degradation rate of Delta proteins. \(\sum \limits _{NG(i)}{{{D}_{i}}}\) represents the concentration of the Delta proteins (in the i ^{th} cell) that bind to the Notch proteins of the contacting cells of the i ^{th} cell, and N G(i) denotes all the neighbors of the i ^{th} cell. f _{1} represents the binding rate of Delta ligands to Notch receptors. λ _{ N } and d _{2} represent the production rate and the degradation rate of inactive Notch proteins, respectively. f _{2} represents the binding rate of Delta ligands to Notch receptors. In \(\frac {a{{N}_{i}}}{b{{D}_{i}}+{{N}_{i}}}\), a represents the conversion rate of Notch proteins from the inactive state to the active state, while b describes the inhibition effect of Delta on Notch in the same cell i. d _{3} is the degradation rate of active Notch proteins.
Simulation and analysis methods
We encode our mathematical model and perform all simulations with Matlab 2016. Besides, we employ the following analysis methods for our model.
Local stability analysis
We conduct local stability analysis to explore stability conditions of our model, with which we further determine values of parameters.
A set of autonomous ordinary differential equations can be written in the following vector form:
where x=(x _{1},x _{2},…,x _{ n }) is the state vector and f=(f _{1},f _{2},…,f _{ n }). The Jacobian matrix is:
Assume x ^{∗} is an equilibrium point, i.e., f(x ^{∗})=0. J ^{∗} is the Jacobian matrix evaluated at x ^{∗}. The equilibrium point x ^{∗} is stable if all the eigenvalues of the characteristic equation of J ^{∗} have negative real parts [22].
Sensitivity analysis
We conduct sensitivity analysis to investigate the significance of model parameters. Sensitivity analysis studies the uncertainty of the output of a model caused by the uncertainty of the inputs of the model [23, 24]. Assume the output (Y) of a model can be represented as the following equation:
where x _{ i } (i=1,2,…n) represent the input variables or factors. Each time, we analyze the sensitivity of one factor (or parameter) by assigning a perturbation term ε β to the factor, e.g., x _{ i }, where ε denotes the disturbance intensity and β is a random number that satisfies the uniform distribution between [1, 1]. We randomly sample values from the interval [x _{ i }−ε,x _{ i }+ε] and then compute the mean and variance of the output as follows:
where Y _{ k } (k=1,2,…,K) are values of the output Y.
Results and discussion
In what follows, we present analysis and simulation results for the model.
Local stability analysis
System (1) gives the dynamic equations of NC cells, which is impossible to be analyzed for a big NC in theory. However, due to the locality of the DeltaNotch interaction, we can still obtain valuable insight into the whole system by only considering two cells to explore their dynamic behavior such as equilibria and stability [25–27] from the theoretical point of view. The mathematical model of two cells is given as follows:
In the following, we will take λ (the production rate of Delta proteins) as an example to study the equilibria and stability of System (2).
(1) When λ=0, which corresponds to misexpression of Delta, we obtain the equilibrium
where \(D_{1}^{0}=0\), \(N_{1}^{0}=\frac {{{\lambda }_{N}}a}{{{d}_{2}}}\), \(A_{1}^{0}=\frac {a}{{{d}_{3}}}\), \(D_{2}^{0}=0\), \(N_{2}^{0}=\frac {{{\lambda }_{N}}a}{{{d}_{2}}}\) and \(A_{2}^{0}=\frac {a}{{{d}_{3}}}\).
(2) When λ>0, which corresponds to normal or overexpression of Delta, we obtain the equilibrium
where \(D_{1}^{1}=\frac {\lambda }{\left (1+\Delta A_{1}^{1}\right)\left ({{d}_{1}}+{{f}_{1}}\right)}\), \(N_{1}^{1}=\frac {b{{d}_{3}}\lambda A_{1}^{1}}{\left (1+\Delta A_{1}^{1}\right)\left ({{d}_{1}}+{{f}_{1}}\right)\left (a{{d}_{3}}A_{1}^{1}\right)}\), \(D_{2}^{1}=\frac {\lambda }{\left (1+\Delta A_{1}^{1}\right)\left ({{d}_{1}}+{{f}_{1}}\right)}\), \(N_{2}^{1}=\frac {b{{d}_{3}}\lambda A_{1}^{1}}{\left (1+\Delta A_{1}^{1}\right)\left ({{d}_{1}}+{{f}_{1}}\right)\left (a{{d}_{3}}A_{1}^{1}\right)}\) and \(A_{2}^{1}=A_{1}^{1}\). \(A_{1}^{1}\) is the solution of
Now, let
Then Eq. (3) becomes
Because \({{m}_{1}}=d_{3}^{2}({{d}_{1}}+{{f}_{1}})\Delta \ne 0\), we further obtain
Thus, we have
Please note that the equilibria of System (2) we obtain above mean that the system will converge to a stable state after a period of time.
Now, we further study the stability of System (2) at equilibrium E _{0}. The Jacobi matrix [27] at E _{0} is computed as follows:
Assume S is the eigenvalue of the characteristic equation \(\phantom {\dot {i}\!}{{J}_{{{E}_{0}}}}\). Then the characteristic equation is calculated as follows [27]:
After some intermediate computation steps, we have
It is obvious that all of the six eigenvalues of Eq. (4) are negative, which gives the following conclusion.
Lemma 1
The equilibrium E _{0} is locally asymptotically stable.
For example, a simulation in this case is given in Fig. 6, with the following parameter values: λ _{ N }=0.02, d _{1}=d _{2}=d _{3}=0.01, f _{1}=f _{2}=0.665, a=0.012, b=69, Δ=10^{6} and λ=0.
Next, we investigate the stability of the equilibrium E _{1}. The Jacobi matrix at E _{1} is given as follows:
The characteristic equation is:
Let
and
Then we obtain
As \(D_{1}^{1}=D_{2}^{1},N_{1}^{1}=N_{2}^{1},A_{1}^{1}=A_{2}^{1}\), we have F _{1}=F _{3}. We further have
where
Similarly, we have
where
Therefore, Eq. (5) becomes
Using the RouthHurwitz criterion for Eq. (6), we obtain
and
Namely, M _{1} M _{2}−M _{3}>0, and M _{1} M _{2}−M ^{′} _{3}>0. Because M _{3}>M ^{′} _{3}, we only need the condition M _{1} M _{2}−M _{3}>0 to be satisfied. That is, if \(\frac {{{M}_{3}}}{{{M}_{1}}{{M}_{2}}}<1\), the eigenvalues of F _{1}+F _{2}=0 and F _{1}−F _{2}=0 are negative; thus, E _{1} is locally stable.
Lemma 2
If \(\frac {{{M}_{3}}}{{{M}_{1}}{{M}_{2}}}<1\), the equilibrium E _{1} is locally asymptotically stable; if\(\frac {{{M}_{3}}}{{{M}_{1}}{{M}_{2}}}>1\), it is unstable.
For example, a simulation in this case is given in Fig. 7 with the following parameter values: λ _{ N }=0.027, d _{1}=d _{2}=d _{3}=0.01, f _{1}=f _{2}=0.665, a=0.012, b=69, Δ=10^{6} and λ=35.
Please note that the equilibria and stability analysis above show that the model would reach an equilibrium in either case of λ, which corresponds to that the system would result in a stable pattern in either normal or mutant conditions. If there are small disturbances of parameters (e.g., small environmental noises), the system would always converge to a stable state after a period of time when λ=0 according to Lemma 1, but the system would converge to a stable state for λ>0 when the stability conditions in Lemma 2 are satisfied. Besides, we will also use Lemma 2 to carefully choose appropriate parameter values for preserving the stability of the system. This is expected to be shown in the the boundary formation model of the Drosophila large intestine.
Deterministic simulation results
However, it is impossible to analyze the model above with a large number of cells in theory, e.g., a model with 60 cells would result in 180 dimensional ordinary differential equations. In this section, we will explore the model using numerical simulation.
The simulation starts from a prepattern of Delta expression in normal large intestines, that is, Delta is expressed only in the ventral region [13]. For this, we set the production rate of the Delta level as follows: 0 for cells at the first three rows (dorsal cells) and λ (greater than 0) for cells at the other rows (ventral cells). We set the time span to [0, 1500]. The parameter values are given in Table 1; except λ and λ _{ N }, other parameter values are fixed throughout the paper. We use the same setting for all the following simulations.
Besides, we determine a boundary cell using the following rule. If the concentration of the active Notch proteins (A) in a cell is equal to or greater than 1.0, we regard the cell as a boundary cell.
In the following, we run the model given in System (1) with the parameter values given in Table 1 to obtain deterministic simulation results. We first validate the ability of our model to reproduce the wildtype result observed in the wet lab. We run simulation and obtain the result in the wildtype condition, illustrated in Fig. 8a. We can see that a single strand of boundary cells that are abutting the ventral region is produced, which is consistent with the experimental observation (corresponding to Scenario 1 and Fig. 3a). We also run the model in other conditions such as overexpression of Notch, and obtain simulation results (see Fig. 8bc), which are also consistent with experimental observations (see Fig. 3 b1c2). That is, our model has the ability to reproduce the experimental observations obtained in the wet lab.
Phenotype change due to the change of Delta or Notch expression
During the boundary formation in the Drosophila large intestine, different Delta or Notch expression (normal, misexpression or overexpression) may result in different boundary cell distributions (phenotypes) [1]. In our model, we try to incorporate all these phenotypes reported so far. As shown in System (1), Delta or Notch expression is controlled by the rates λ and λ _{ N }, respectively. Therefore, we need to map different Delta and Notch expression observed in experiments to different values of λ and λ _{ N } in the model.
To do this, we explore the phenotype change due to the change of Delta or Notch expression in their parameter space and then determine what parameter value results in what kind of phenotype. As well, we cannot adopt a theoretical analysis due to the large number of equations involved, and instead, we still use numerical simulation. Besides, each time we tune only one parameter, either λ or λ _{ N }, by fixing the other.
Phenotype change due to the change of Delta expression
By fixing λ _{ N }=0.0005, we gradually increase λ from 0 to 12 with a step 0.1 and then run simulation to obtain the simulation result, given in Table 2.
From the table above, we can see when λ>0.2, the first three rows are stable boundary cells, while when λ≤0.2, all the five rows are boundary cells. This means there are only two phenotypes appearing due to the change of Delta expression.
Phenotype change due to the change of Notch expression
By fixing λ=9, we gradually increase λ _{ N } from 0 to 0.1 with a small step 0.0001 and then run simulation to obtain the simulation result, given in Table 3.
Table 3 shows with the increase of λ _{ N }, the number of boundary cells increases, which results in a couple of different phenotypes. Besides, we should notice that a tiny parameter change would result in a big change of the number of boundary cells, so we can deduce that parameter λ _{ N } is very sensitive.
According to the analysis result above, we finally map different Delta and Notch expression observed in experiments to different values of λ and λ _{ N } in the model. Besides, we can reproduce all known biological phenotypes we have found by combining different values of λ and λ _{ N }.
Simulation results with noises of Notch expression
As described in [1], the boundary formation of the Drosophila large intestine is greatly affected by temperature; the random effect of temperature may cause fluctuations of gene expression levels, resulting in different phenotypes, in which boundary cells may randomly be distributed in the dorsal and ventral domains. In this section, we will discuss how to use our model to obtain these random phenotypes.
From Tables 2 and 3, we know that only λ _{ N } is sensitive to System (1), so we only add noises to λ _{ N }. Then, System (1) with a random noise becomes the following form:
where e is a small noise.
Running this model with different λ _{ N } values, we can obtain different random phenotypes. Figure 9 illustrates some phenotypes. So far, we have illustrated that our model reproduces all the phenotypes corresponding to those four scenarios we give above.
Perturbation influences of binding and inhibition parameters on boundary formation
During the development of the Drosophila large intestine, there are many environmental factors or noises, which may vary and even heavily affect the boundary formation. To make the model more realistic, we need to consider noises into the model, e.g., adding a random term to each parameter of interest, and analyze the model in noisy conditions. In the section above, we have discussed the random effects of Notch expression on the boundary formation. But in this section, we will further analyze the perturbation influences of the binding and inhibition parameters of the DeltaNotch pathway on the boundary formation. That is, we will explore the inhibition parameters, a and b, and binding parameters, f _{1} and f _{2}.
For these parameters, which ones are of interest? To answer this, we may alternatively find those parameters which have significant influences on the dynamics of the model when they vary. For this, we use sensitivity analysis on the mathematical model to study the influence of parameter perturbation [28]. After that, we can design experiments by considering noises in vivo, and carefully add appropriate noise items to key parameters. We will explore this in the case of overexpression of Notch proteins. For the other cases, we can do the same exploration.
Next, with the parameter setting given in Fig. 8b (of course we can use any other setting), we will explore the influence of parameter perturbation on the boundary formation by changing parameters, a, b, f _{1} and f _{2}, to different perturbation intensities.
Perturbation influence of parameter a
First, we explore the perturbation influence of parameter a. All the parameters take the values given in Fig. 8b. In order to do this, we add a perturbation term ε β to parameter a, where ε denotes the disturbance intensity and β is a random number that satisfies the uniform distribution between [1, 1]. Then, System (1) with a random disturbance becomes the following form:
At the beginning of each simulation, a random noise ε β _{ i } for each cell i is generated for a by sampling β _{ i }. By varying ε between [0, 0.001] and performing simulation for 50 times at each ε level, we compute the mean and standard error of the number of boundary cells. For example, Fig. 10 gives the result when ε=10^{−7}, ε=10^{−6}, ε=10^{−5} and ε=10^{−4}, respectively.
From Fig. 10, we can see that a small perturbation of parameter a may cause a big change of the number of boundary cells (mean), i.e., the bigger the perturbation intensity is, the less the number of boundary cells is (from 36 boundary cells without perturbation to about 15 boundary cells with a perturbation intensity ε=10^{−4}). Furthermore, from the error bars, we know that the bigger the perturbation intensity is, the bigger the standard error is, although the smaller the mean is.
Perturbation influence of parameter b
We further explore the perturbation influence of parameter b. Similarly we add a perturbation term ε β to b, i.e., replacing all b with b+ε β in System (1). We also simulate the model in the same way as for parameter a, and obtain the mean and standard error of the number of boundary cells when ε=0.01, ε=0.05, ε=0.1 and ε=0.2, respectively, illustrated in Fig. 11.
Figure 11 shows that a big perturbation causes a small change (e.g., there are 36 boundary cells without perturbation and around 34 boundary cells with the perturbation intensity ε=0.2), which is quite different from the effect of parameter a. Therefore, the perturbation influence of parameter b is much smaller than that of a.
Perturbation influence of parameters f _{1} and f _{2}
We further explore the perturbation of f _{1} and f _{2}. For simplicity, we assume f _{1}=f _{2}=f. Similarly, we add a perturbation term ε β to f _{1} and f _{2}, i.e., replacing all f _{1} and f _{2} with f+ε β in System (1). We also simulate the model in the same way as above, and obtain the mean and standard error of the number of boundary cells when ε=0.001, ε=0.005, ε=0.01 and ε=0.02, respectively, given in Fig. 12.
From Fig. 12, we can see, similar to parameter b, the perturbation of parameter f does not have an obvious influence on the number of boundary cells in the given perturbation intensities.
In conclusion, the analysis results above show that the perturbation of parameter a has an obvious influence on the boundary formation, while the other three parameters not. Therefore, when we consider noises in the model, we only add a random term to parameter a, and neglect others. That is, we finally obtain a model with random noises, which is given in System (8). With this model, we can make reasonable predictions in noisy conditions by setting the appropriate intensity of the noise according to the analysis result above (see, e.g., Fig. 13).
Conclusion
In this paper, we give a mathematical model for the DeltaNotch dependent boundary formation in the Drosophila large intestine. We aim to use this model to better interpret related biological phenomena and therefore we perform not only theoretical but also simulation analysis to achieve the goal. By combining different analysis techniques, we finally confirm that the model is both mathematically and biologically sound and is sufficient to interpret related experimental findings of the boundary formation in the Drosophila large intestine. Moreover, by modulating parameters, our simulation can generate various abberant patterns of boundary cells that have not been reported in biological studies so far. Biologically, these varieties of phenotypes are assumed to be results of variation of gene products in the DeltaNotch pathway. To further verify the validity of our simulation method, we are now trying to develop a biological experimental system to manipulate levels of forcedgene expression. In a next step, we will shift the purpose of the model from interpretation to prediction. That is, we will use this model to make different predications in different conditions and then help biologists to have an idea of detailed mechanisms of complicated biological systems, and to design experiments to validate the predictions obtained by our simulation.
References
 1
Takashima S, Yoshimori H, Yamasaki N, Matsuno K, Murakami R. Cellfate choice and boundary formation by combined action of Notch and engrailed in the Drosophila hindgut. Dev Genes Evol. 2002; 212:534–41.
 2
Murakami R, Shiotsuki Y. Ultrastructure of the hindgut of Drosophila larva, with special reference to the domains identified by specific gene expression patterns. J Morphol. 2001; 248:144–50.
 3
Takashima S, Murakami R. Regulation of pattern formation in the Drosophila hindgut by wg, hh, dpp, and en. Mech Dev. 2001; 101:79–90.
 4
Hamaguchi T, Takashima S, Okamoto A, Imaoka M, Okumura T, Murakami R. Dorsoventral patterning of the Drosophila hindgut is determined by interaction of genes under the control of two independent gene regulatory systems, the dorsal and terminal systems. Mechanisms of Development. 2012; 29:236–43.
 5
Kumichel A, Knust E. Apical localisation of crumbs in the boundary cells of the Drosophila hindgut is independent of its canonical interaction partner stardust. PLoS ONE. 2014; 9(5):94038.
 6
Wang Y, Zhang SX, Chen L. Modelling biological systems from molecules to dynamical networks. BMC Systems Biology. 2012; 6(Suppl 1):1.
 7
Grimm O, Coppey M, Wieschaus E. Modelling the bicoid gradient. Development. 2010; 137:2253–64.
 8
Murray JD. Mathematical Biology. Berlin: Springer; 2003.
 9
Nakamasu A, Takahashi G, Kanbe A, Kondo S. Interactions between zebrafish pigment cells responsible for the generation of Turing patterns. PNAS. 2009; 106(21):8029–34.
 10
Liu F, Blätke M, Heiner M, Yang M. Modelling and simulating reactiondiffusion systems using coloured Petri nets. Comput Biol Med. 2014; 53:297–308.
 11
Morelli LG, Uriu K, Ares S, Oates AC. Computational approaches to developmental patterning. Science. 2012; 36:187–91.
 12
Iwaki D, Lengyel JA. A DeltaNotch signaling border regulated by engrailed/invected repression specifies boundary cells in the Drosophila hindgut. Mech Dev. 2002; 114:71–84.
 13
Matsuno H, Murakami R, Yamane R, Yamasaki N, Fujita S, Yoshimori H, Miyano S. Boundary formation by notch signaling in Drosophila multicellular systems: experimental observations and gene network modeling by genomic object net. In: Proceedings of Pac Symp Biocomput: 37 January 2003; Kauai, Hawaii. World Scientific Press: 2003. p. 152–63.
 14
Hayward P, Kalmar T, Arias AM. Wnt/notch signalling and information processing during development. Development. 2008; 135:411.
 15
Muskavitch MA. DeltaNotch signaling and Drosophila cell fate choice. Dev Biol. 1994; 166(2):415–30.
 16
Bray SJ. Notch signalling: a simple pathway becomes complex. Nat Rev Mol Cell Biol. 2006; 7:678–89.
 17
Kopan R, Ilagan MXG. The canonical Notch signaling pathway: unfolding the activation mechanism. Cell. 2009; 137:216–33.
 18
Baker NE. Notch signaling in the nervous system: pieces still missing from the puzzle. BioEssays. 2000; 22:264–73.
 19
Collier JE, Monk NA, Maini PK, Lewis JH. Pattern formation by lateral inhibition with feedback: a mathematical model of DeltaNotch intercellular signaling. J Theor Biol. 1996; 183(4):429–46.
 20
Boareto M, Jolly MK, Lu M, Onuchic JN, Clementi C, BenJacob E. JaggedDelta asymmetry in Notch signaling can give rise to a sender/receiver hybrid phenotype. PNAS. 2015; 22(5):402–9.
 21
Sprinzak D, Lakhanpal A, LeBon L, Santat LA, Fontes M, Anderson GA, GarciaOjalvo J, Elowitz MB. Cisinteractions between Notch and Delta generate mutually exclusive signalling states. Nature. 2010; 465:86–90.
 22
Boyce WE. Elementary Differential Equations and Boundary Value Problems, Binder Ready Version, (10th Edition). USA: Wiley; 2012.
 23
Hamby DM. A review of techniques for parameter sensitivity analysis of environmental models. Environ Monit Assess. 1994; 32(2):135–54.
 24
Turányi T. Sensitivity analysis of complex kinetic systems: Tools and applications. J Math Chem. 1990; 5(3):203–48.
 25
Seydel R. Practical Bifurcation and Stability Analysis. New York: Springer; 2010.
 26
Lakshmikantham V, Leela S, Martynyuk AA. Stability Analysis of Nonlinear Systems. Switzerland: Springer; 2015.
 27
Khanh NH, Huy NB. Stability analysis of a computer virus propagation model with antidote in vulnerable system. Acta Math Sci. 2016; 36B(1):49–61.
 28
Shivamoggi BK. Perturbation Methods for Differential Equations. New York: Springer; 2003.
Acknowledgment
We would like to thank Prof. Monika Heiner for fruitful discussions and good comments. We also would like to thank the anonymous reviewers for their constructive comments.
Funding
This work was mainly supported by National Natural Science Foundation of China (61273226), which also financed the articleprocessing charge. FL was supported by JSPS Invitation Fellowships for Research in Japan (L15542) for his stay in Japan during this research. None of the funding agencies had any role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.
Availability of data and materials
The data supporting the results of this article are included within the article.
About this supplement
This article has been published as part of BMC Systems Biology Volume 11 Supplement 4, 2017: Selected papers from the 10th International Conference on Systems Biology (ISB 2016). The full contents of the supplement are available online at https://bmcsystbiol.biomedcentral.com/articles/supplements/ volume11supplement4.
Author information
Author notes
Affiliations
Contributions
FL conceived the study and developed the model. DS performed the simulations and analyzed the model. HM and RM provided the biological knowledge. All authors read and approved the final manuscript.
Corresponding author
Correspondence to Fei Liu.
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Published
DOI
Keywords
 Boundary formation
 Drosophila large intestine
 DeltaNotch signaling pathway
 Local stability analysis
 Simulation validation
 Perturbation analysis