Modeling and analysis of the Delta-Notch dependent boundary formation in the Drosophila large intestine

Background The boundary formation in the Drosophila large intestine is widely studied as an important biological problem. It has been shown that the Delta-Notch signaling pathway plays an essential role in the formation of boundary cells. Results In this paper, we propose a mathematical model for the Delta-Notch 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.

In order to understand the mechanism of the Delta-Notch pathway, there have been some mathematical models proposed. For example, Collier et al. proposed a simple ordinary differential equation (ODE) model for the Delta-Notch signaling pathway, and discussed numerical simulation of multiple-cell systems [19]. Boareto et al. devised a theoretical framework that includes a couple of ODEs to explore the effects of Jagged in cell-fate 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 cis-interactions between Notch and Delta generate mutually exclusive signaling states [21]. Specifically, Matsuno et al. analyzed the mechanism of the Notch-dependent boundary formation in the Drosophila large intestine and built a hybrid Petri net model to numerically explore how the Delta-Notch pathway affects the boundary formation in two-dimensional space [13].
In this paper, we propose a mathematical model for the Delta-Notch 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 Delta-Notch 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 two-cell 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 Delta-Notch 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 one-cell-wide strand of boundary cells forms in wild-type embryos.
The Delta-Notch 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 Delta-negative dorsal cell, a Delta-Notch signaling cascade is activated. First the site 2 cleavage of the Notch proteins occurs to generate a transmembrane form (N E ), followed by the Presenilin-dependent 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 Presenilin-mediated site 3 cleavage and thus inhibits Notch signal transduction within Delta-positive 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.  In what follows, we will construct and analyze our model step by step.

Mathematical model
Cell-to-cell interaction mediated by the Delta-Notch 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 Delta-Notch dependent boundary formation in the Drosophila large intestine (see Fig. 4 for a one-cell 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. This pathway shows the interaction of two neighboring cells, which is adapted from [1] 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 The mathematical model is then given as follows: Here, 1 ≤ i ≤ NC. 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 λ 1+ ·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.

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 NG(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 aN i bD 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.  [1,4,13]

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: 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 Delta-Notch 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][26][27] from the theoretical point of view. The mathematical model of two cells is given as follows: (2) 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 mis-expression of Delta, we obtain the equilibrium (2) When λ > 0, which corresponds to normal or overexpression of Delta, we obtain the equilibrium Now, let Then Eq. (3) becomes 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 |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 Next, we investigate the stability of the equilibrium E 1 . The Jacobi matrix at E 1 is given as follows: The characteristic equation is: Then we obtain Similarly, we have Therefore, Eq. (5) becomes Using the Routh-Hurwitz criterion for Eq. (6), we obtain  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 wild-type result observed in the wet lab. We run simulation and obtain the result in the wild-type 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 over-expression of Notch, and obtain simulation results (see Fig. 8b-c), which are also consistent with experimental observations (see Fig. 3b1-c2). That is,

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, mis-expression or over-expression) may result in different boundary cell distributions (phenotypes) [1]. In our model, we try to incorporate all these phenotypes  (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 By fixing λ = 9, and gradually increasing λ N from 0 to 0.1 with a small step 0.0001, we obtain the phenotype change with the change of Notch expression 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 Delta-Notch 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 over-expression 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,  13 Two simulation runs for parameter a with a perturbation intensity ε = 10 −6 . All the parameters take the values given in Fig. 8b ε = 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 Delta-Notch 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 Delta-Notch pathway. To further verify the validity of our simulation method, we are now trying to develop a biological experimental system to manipulate levels of forced-gene 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.