A continuation method for spatially discretized models with nonlocal interactions conserving size and shape of cells and lattices

In this paper, we introduce a continuation method for the spatially discretized models, while conserving the size and shape of the cells and lattices. This proposed method is realized using the shift operators and nonlocal operators of convolution types. Through this method and using the shift operator, the nonlinear spatially discretized model on the uniform and nonuniform lattices can be systematically converted into a spatially continuous model; this renders both models point-wisely equivalent. Moreover, by the convolution with suitable kernels, we mollify the shift operator and approximate the spatially discretized models using the nonlocal evolution equations, rendering suitable for the application in both experimental and mathematical analyses. We also demonstrate that this approximation is supported by the singular limit analysis, and that the information of the lattice and cells is expressed in the shift and nonlocal operators. The continuous models designed using our method can successfully replicate the patterns corresponding to those of the original spatially discretized models obtained from the numerical simulations. Furthermore, from the observations of the isotropy of the Delta–Notch signaling system in a developing real fly brain, we propose a radially symmetric kernel for averaging the cell shape using our continuation method. We also apply our method for cell division and proliferation to spatially discretized models of the differentiation wave and describe the discrete models on the sphere surface. Finally, we demonstrate an application of our method in the linear stability analysis of the planar cell polarity model.


Introduction
The development of multicellular organisms is regulated by intercellular communication and signaling pathways of various types. These include diffusible proteins acting as ligands and cell membrane proteins communicating with the neighboring cells. In the last fifty years, approaches comprising mathematical models and numerical simulations have been extensively used to understand the mechanisms underlying the biological phenomena. It is a common practice to divide a region of interest either into square or hexagonal elements representing cells, as shown in Fig. 1; this allows for the discrete spatial independent variables to be used. We also assume that the unknown dependent variables of the model are uniform on the lattices. With these preconditions, we model the phenomena on the lattices mathematically. In this paper, we label the mathematical models with the discrete spatial independent variable as discrete models, and the ones with the continuous spatial independent variable as the continuous models. Modeling the phenomena on the divided lattices often demonstrates good reproducibility and presents good agreement with experimental results.
One of the examples in which the cellular interaction is conserved in various organisms is the Delta-Notch signaling. The intercellular communication in this signaling is based on the informational exchange between neighboring cells (Collier et al. 1996;Sato et al. 2013;Yasugi et al. 2008). The function of the Delta-Notch signaling is known as lateral inhibition. During neural development of the fly embryo, binding of the Delta ligand to the Notch receptor suppresses the expression of achaete-scute complex (AS-C) proneural genes. On the contrary, AS-C genes induce Delta expression. Consequently, signal-sending cells demonstrate a high level of the AS-C genes, while signal-receiving cells express low level of the AS-C genes. During embryogenesis, neuroepithelial cells (NEs) that express high Delta and AS-C differentiate into neural progenitor cells. In contrast, the surrounding cells express low levels of AS-C genes and differentiate into non-neuronal cells. In accordance with these interactions between the neighboring cells, the expression patterns of Delta and Notch activation show a salt-and-pepper like pattern that distinguish the neuronal cells from the nonneuronal ones, as shown in Fig. 2. Information regarding discreteness, such as the size and shape of each cell, affects the entire pattern in the developmental process,

Notch mediated lateral inhibition
Salt-and-pepper pattern Delta therefore, modeling in the framework of the discrete model is compatible with the phenomenon described (Collier et al. 1996;Lehotzky and Zupanc 2019;Sato et al. 2016). A discrete model shows good reproducibility of the experimental results for the differentiation propagation in the developing fly brain (Sato et al. 2016;Tanaka et al. 2018). Using the described type of discrete model for the Delta-Notch interaction, the salt-and-pepper pattern appearance and regulated differentiation propagation in the fly brain have been explained (Collier et al. 1996;Jörg et al. 2019;Sato et al. 2016;Tanaka et al. 2018). It is well known that the function of Delta-Notch signaling is diverse and Notch activation shows several different patterns. For example, Notch activation oscillates in the segmentation in vertebrates and progresses unidirectionally in the fly optic lobe development Kageyama et al. (2012).
Another good example of a biological system that utilizes discrete modeling is the planar cell polarity (PCP) (Adler 2002;Goodrich and Strutt 2011). The intercellular proteins and membrane in the cells of the fly wing interact with each other, among the neighboring cells, and they are localized asymmetrically. Owing to this asymmetric localization, the direction of the epithelial hair in the wing of a fly is determined Ayukawa et al. (2014). It is reported that the discrete model proposed by this paper can reproduce the biological experiments of the PCP.
The analytical study of the discrete models was further conducted to specify the function of the intercellular interactions and the discreteness in the dynamics. The analytical results for the discrete type of reaction diffusion systems have been reported in Bates et al. (2001), Chow (2003). The results related to the traveling wave solutions in the system of the discrete models are reported in Guo et al. (2019), Hupkes and Sandstede (2010), Straatman and Hupkes (2019).
Although the discrete models are useful in describing the abovementioned behavior and dynamics phenomenologically, the analysis of discrete models is rather difficult in general, and the technique of analyzing a discrete model is poor compared to the analysis of the continuous models. For example, as the discrete models usually comprise numerous unknown variables, it is usually difficult to compute the eigenvalues in case of higher-dimensional space. Thus, the methods of analysis for the partial differential equation are being reconstructed, such that they are applicable to the discrete models (Chow 2003). Moreover, discrete models are not compatible with the description of regional expansion caused by cell division and the three-dimensional information. In order to overcome the difficulties mentioned, the limit of the cell and lattice size is often set to zero, and the differential operator is derived. There are some examples of the continuation. The discrete models to prion dynamics and coagulationfragmentation processes were investigated mathematically through the continuation method using the piece-wise constant functions (Crampin et al. 1999;Laurençot and Mischler 2002). In these papers, the continuous models were derived on the lattices with sufficiently small length. By taking the limit of the lattice size, the convergence is rigorously shown, and the integro-differential equations were derived.
However, the method taking the limit of the lattice size to zero may cause a problem because the patterns caused by the spatially discretized structures, such as lattice and cell membrane, sometimes disappear in the continuous models. In this paper, we questioned if it is possible to convert a discrete model into a continuous model while retaining the size and shape of cells and lattices. In the light of this question, we propose a novel continuation method for the spatially discretized models, while conserving the size and shape of cells and lattices. In our proposed method, we perform the continuation part of the process by introducing the shift operators, instead of deriving the differential operators. Thus, the nonlinear discrete models can be systematically converted into continuous models with spatially discretized structures. Moreover, by reducing the shift operators using the integral operators of the convolution type with suitable kernels, we propose a nonlocal evolution equation that can approximate the solution of the spatially discrete model, for the application to the biological experiments and mathematical analysis. The approximation of nonlinear discrete models by nonlocal evolution equations in the one-dimensional periodic region is assured by singular limit analysis. In the continuous models, the information of the size and shape of cells and lattice was reflected in the part of shift operator and the kernel in the convolutions. Furthermore, we confirmed the isotropy of the Delta-Notch signaling system for the irregularly shaped cell in the fly brain. According to the biological results, we propose a radially symmetric kernel in the nonlocal evolution equation for discrete models and prove the effectiveness of the kernel by replicating the spatially discretized patterns in the numerical simulations. As a result of the description using nonlocal evolution equation, we model the cell division and proliferation in the discrete model for the wave of the differentiation and extend the model to the sphere surface. Moreover, we show that linear stability analysis can be performed by the continuation method.
This paper is organized as follows: In Sect. 2, we first introduce the concept of our continuation method by modifying the general discrete model. Our continuation method is characterized by the combination of the shift operator and the Friedrichs mollifiers as the convolution kernel based on the lattice shape. In Sect. 2.2 we state the result of the singular limit analysis of the discrete models and nonlocal evolution equations in the general form of the spatial interactions. In Sect. 2.3 we explain that our method can be extended to the discrete models on nonuniform lattice. In Sect. 2.5.3, we introduce the radially symmetric kernel for the averaged shape of the cells, based on the results of the real biological experiments for Delta-Notch signaling in the fly brain. In Sect. 3, we show the results of numerical simulations in biological applications of the continuation method: the proneural wave in the fly brain and planar cell polarity in the fly wing. Our results indicate that the continuation method using shift operators and integral operators can be applicable for diverse multicellular systems.
2 Continuation method with shift and convolution operator for discrete models

Scalar equation in one-dimensional space
In this section, we explain the concept of the continuation method, while retaining the shape and size of cells and lattices. First, we describe the continuation method applied on a typical discrete model containing the intercellular interaction terms and the reaction term. In this paper, we do not distinguish between the spatial and intercellular interactions. Suppose N cells of the uniform length l > 0 are packed in the one-dimensional space, then the following discrete model is considered: where u i = u i (t) denotes the concentration or density of some substances on the ith cell c i at time t > 0, f : R 3 → R is the function corresponding to the intercellular interactions, and g : R → R is the function for the reactions. Setting the one-dimensional space as we impose the periodic boundary condition u 0 (t) = u N (t), and u 1 (t) = u N +1 (t).
The linear intercellular interaction can be defined as: where a i , (i = −1, 0, 1) are constants. The typical examples of the function f are diffusion and lateral inhibition such as the Delta-Notch interaction given by where the denominator of f lat is the total number of neighboring cells referred from Collier et al. (1996), and the sign of the lateral inhibition f lat can be changed in the system. If the dynamics of u i is more influenced by the other cells than by the neighboring cells, f becomes the nonlocal interactions. As introduced in Doumic et al. (2009), Laurençot and Mischler (2002), we will utilize the piece-wise constant functions for our continuation method. For equation (1) with i = 1, . . . , N on each cell c i , we define the characteristic function as and also we define at position x ∈ T and at time t > 0. For the continuous method of the discrete model, we set the following assumption: For any N there exists a unique global solution u(x, t) ∈ C([0, T ], L 1 (T)) of (1).
As in Proposition 2 and "Appendix B", the existence and uniqueness of the global solution of (1) is shown for specified functions f and g. Changing the variable in the ith equation (1) by multiplying the unknown function u i by the characteristic function

As we can compute
for j = 0, ±1, · · · , ±N , we obtain As the spatially independent variable is continuous, the discrete model (1) is successfully converted into a continuous model. The equation of (7) is point-wisely equivalent to the equation (1). Thus, if the initial conditions of equation (1) and (7) are the same, the solutions are equivalent as described by the following remark: Remark 1 Using the initial data of discrete model {u i (0)} N i=1 , and imposing the initial data as u 0 (x) , the solution of continuous model (7) is equivalent to that of the discrete model (1). Furthermore, to apply the continuous model to the experiments and analyze conveniently, we approximate the shift operator using the convolution on the mollifier. We define the shift operator as follows: The shift operator is regarded as the convolution of the shifted Dirac Delta function δ l := τ l δ = δ(x + l), and we can describe the model (7) as follows: Here we suppose that the Dirac Delta function is periodic with Nl, i.e., δ(x) = δ(x + Nl), and we define the convolution k * v with respect to x in T as where T can be replaced with a given region in this paper. Setting the Friedrichs mollifier with a small parameter 0 < ε 1 as with a constant for the normalization of integration C 0 > 0, we assume that ρ ε is also periodic with Nl. We use the symbol ρ ε for the mollifier in higher-dimensional case.
Approximating the Dirac Delta function by the mollifier ρ ε (x), we have where the shifted mollifier is given by ρ ε,l := ρ ε (x + l), and we denote the unknown variable by u ε (x, t) as the solution of this equation depends on ε. If the intercellular interaction f is linear, we derive the typical nonlocal evolution equation by summarizing the kernel of the convolution as follows: where we put the kernel as K = a −1 ρ ε,−l + a 1 ρ ε,l . Consequently, we have the following nonlocal evolution equation: Such type of a nonlocal evolution has been analyzed in numerous papers (Bates et al. 1997;Coville and Dupaigne 2007). Figure 3 shows the profile of the kernel K for f Δ . Figure 4 shows the results of the numerical simulations for both continuous and discrete heat equations fed with the spatially discretized initial data.
As in Fig. 4a, it is observed that the solution of the equation with shift operator is not continuous until the solution attains the steady-state in the numerical simulation. On the contrary, as in Fig. 4b, it is observed that the solution of the equation with the mollifier becomes continuous before the solution attains the steady-state in the numerics.
Remark 2 If f and g are linear, and if u ε = u * ρ ε , u ε becomes the solution of (7). This is owing to the linearity of the convolution operator, which can be described as follows; using the convolution of the mollifier in the equation (7), we have We compute as follows thereby satisfying Eq. (7).
Furthermore, our proposed method is consistent with the continuation method where the cell limit was set to 0 or lattice size was set to l, as we can derive the differential operator by setting the limit l → +0 after applying our continuation method.
Remark 3 If f is equal to f Δ , we see that Even if the intercellular interaction is nonlocal, which means it is affected by not only the neighboring cells but also the other cells, our continuation method is applicable to the discrete model, in a similar way. A discrete model in which intercellular interactions are influenced by the cells other than neighboring cells, is given as follows: where [·] is the Gauss's symbol, and f : R 2[N −1/2]+1 → R is a function corresponding to the interaction here. If f is linear, the function f is generally defined as where i=−[(N −1)/2] are constants. Following the calculation in (6), we derive the equivalent continuous model as follows: By describing the nonlocal interactions using the convolution of the mollifier ρ ε , the kernel is provided by and thus, the nonlocal evolution, which can approximate the solution of (P S ), is given as follows:

Singular limit analysis
In this subsection, we will explain that the solution of (P ε ) is sufficiently close to that of (P S ) in L 2 (T) by singular limit analysis. As the interaction f with the form in (P D ) includes the case of intercellular interaction in the discrete model, we deal with the equations (P S ) and (P ε ) in this analysis. We firstly assume that f is the form of (10). For the condition g, we assume that there exist positive constants g 0 , g 2 , g 4 and a nonnegative constant g 1 , g 3 such that for u, v ∈ R, A typical example of g is g(u) = u(1 − u 2 ), where g 0 = g 2 = g 4 = 1, g 1 = g 3 = 0 and p = 3. First, we calculate the fundamental solution of the (P S ) without the reaction term g(u).

Proposition 1 The fundamental solution of u t
where {q k } N k=1 are real constants determined by the initial data.
The proof is in "Appendix B". Since the equation (P S ) is equivalent to (P D ), and the associated matrix from f of the system (P D ) is the cyclic, the eigenvalue and the eigenvector can be calculated. Next, we have the uniqueness and global existence of the solutions of (P S ) and (P ε ).
Proposition 2 Assume that f is given by (10), and where C 1 is a positive constant.
The proves are in "Appendix B". We note that each global solution of (P S ) and (P ε ) is in L 2 (T) due to L ∞ (T) ⊂ L 2 (T). Thus, we have global boundedness in L 2 (T) as from the estimations (13) and (14). For the solution of the model (P S ) and (P ε ), we have the following error estimate. Setting the error between the solution of (P S ) and (P ε ) as we have the following convergence result.

Theorem 1 Suppose the same assumptions of Propositions 2 and 3. Let u(x, t) and u ε (x, t) be solutions of (P S ) and (P
where C 5 and C 6 are positive constants independent of ε. Thus, we have The calculation of the energy estimate is put in "Appendix B". From this estimation, the solution of the continuous model (8) converges to that of (7) in L 2 (T) space as ε tends to 0. This implies the solution of nonlocal evolution equation can approximate the solution of discrete model.
The typical example of above f is in the model of the PCP (38). By using the inequality (15), the proof is followed by that of Proposition 2, Proposition 3, and Theorem 1.

Nonuniform lattice in one-dimensional space
In this subsection, we introduce that our continuation method is applicable to the discrete models on nonuniform lattices by adding some conditions as a remark. Labeling the ith cell as c i , (i = 1, . . . , N ) with the nonuniform length l i > 0, we suppose N cells are packed in the one-dimensional space. Let u i = u i (t) be the concentration or density of some substances on c i at time t > 0. Imposing the periodic boundary condition u 0 (t) = u N (t), and u 1 (t) = u N +1 (t) with l 0 = l N , and l 1 = l N +1 , we consider the following discrete model in this subsection: where the definitions of f and g are same as those in (1), and the initial data are given by For the length l i , we define the following functions: for any x ∈ c i . Using the characteristic function (5), we define the piecewise constant function for the shift as follows Changing the variable in the ith equation (16) by multiplying the unknown function u i by the characteristic function χ c i (x), and adding u i (t)χ c i (x) with respect to i = 1, · · · , N , we have As we can compute that we obtain If the initial datum is given by , the solution of continuous model (16) is equivalent to that of the discrete model (17).
For the mathematical analysis, the function l(x) and r (x) can be rendered continuous by using mollifier. We see that l(x) and r (x) belong to we propose an approximation model to (17) as (18) Similarly to Sect. 2.1, we approximate the nonuniform shift operators by the convolutions. We assume (A1) for the function u ε in (18). We can rewrite (18) in the convolution form as (2), we obtain the nonlocal evolution equation where the kernel is given by With ε > 0 and η > 0 the kernel K (x, y) is differentiable. As introduced above, the discrete models on nonuniform lattice can be rendered continuous models. The nonlocal evolution equation (19) is expected to approximate the solutions of that of the original discrete models on nonuniform lattices. The error estimations, the analysis and the application to this proposed model is one of our future works.

System in one-dimensional space
In the case of system we can perform the continuation method similarly to Sect. 2.1. We will explain the method by using the typical reaction diffusion system and Delta-Notch signaling system which are often used for the mathematical modeling in the successive sub-subsections.

Reaction diffusion system
First, we explain the typical reaction diffusion system in the one-dimensional space with periodic boundary conditions. Let u i = u i (t) and v i = v i (t) be the concentration of the diffusive substances on the uniform cells or lattices c i , (i = 1, . . . , N ), respectively. The reaction diffusion system in the framework of the discrete model can be described as follows: where d u , d v > 0 are the diffusion coefficients, f Δ is defined by (3), and g 1 , g 2 : R 2 → R are the functions for reactions in this sub-subsection. For this equation, setting the variables at position x ∈ T and at time t > 0 as we derive the reaction diffusion system performed our continuation method as follows: Similarly to the previous subsection, this above equation is point-wisely equivalent to the equation (20). Indeed, if the initial conditions of (20) and (21) are same, the solutions of (20) and (21) are equivalent. Approximating the shift operator by the mollifier, and describing (20) like the form of (9) , we have the following nonlocal evolution equation which is expected to approximate the solution of the original discrete model (20):

Delta-Notch interaction system
Secondary, we consider the continuation method to the general Delta-Notch interaction system. Let D i = D i (t) and N i = N i (t) be the expression of Delta, and Notch signal in the cell c i , (i = 1, . . . , N ), respectively. The simple description for the Delta-Notch signaling in the framework of the discrete model is given by the following system Collier et al. (1996): where f is the function f lat defined in (4) or the function depending f lat , g 1 and g 2 are the functions for reactions in this sub-subsection, and we replace the literature for the cell number N with M ∈ N for the clear description in this sub-subsection. Similarly to Sect. 2.1, the changing the variables in Eq. (22) as yield the following equation Indeed, if the initial conditions of (22) and (23) are the same, the solutions of (22) and (23) are equivalent for any time t > 0. Regarding the shift operator as the convolution with the Dirac Delta function, we approximate it by the convolution with the mollifier.
If f is linear such as f lat , we have the following nonlocal evolution equation similarly to (9): where K = (ρ ε,−l +ρ ε,l )/2. Even if the number of the unknown variables is increased, our method is applicable to make the discrete model continuous.

Two-dimensional space
In this subsection we will explain our continuous method for the discrete model in the two-dimensional case. As the continuation method for the scalar equation can be applied to the system similarly to one-dimensional case, we firstly deal with the scalar discrete model in two-dimensional space. As considering the model in twodimensional case, the number of terms in the intercellular interaction is increased. Accordingly, the number of the terms of shift or convolution with mollifier is increased in the continuation method. The procedure of continuous method for the reaction term is the same as the explanation of the previous subsections. We set the square region, and impose the periodic boundary condition in this subsection.

Square lattice
We perform the continuation method for the discrete model in the square lattice. This mathematical model corresponds to the situation that square cells are packed in a plane in the development of multicellular organs. Dividing the square region into N 2 square parts of lattices, we label each lattice as c i, j , (i, j = 1, . . . , N ), and denote the horizontal and vertical length of a lattice by l x > 0 and l y > 0, respectively. Thus, the region is described by One divided region corresponds to one cell or lattice as shown in Fig. 1.
For the simplicity, we suppose that region is a regular square. Then the scalar discrete model with intercellular interaction can be described as follows: is denoted by the concentration or density of some substances on c i, j at time t > 0, f : R 5 → R and g : R → R are intercellular and reaction functions, respectively, in this sub-subsection. If f is a linear function, it can be generally written by The typical examples of f are diffusion and lateral inhibition as follows: where f Δ× is referred from Chow (2003). For this discrete model we prepare the following characteristic function at position (x, y) ∈ T 2 : For the Eq. (24) for i, j = 1, . . . , N , we change the variables similarly to Sect. 2.1 by using the characteristic function. Here setting the variable on T 2 as we have from same calculation as that on one-dimensional case. We put the specific calculation in "Appendix C". The discrete model is successfully converted into the continuous model. Similarly to the case in one dimension, approximating the shift operator by the convolution with the mollifier yields the nonlocal evolution equation: where we define the shift operator τ l,m as If f is linear, the description with the kernel is given as follows As mentioned above, our method enables us to derive the continuous model and nonlocal evolution equation for the original discrete model. The profile of the kernel for f Δ on the square lattice is shown in Fig. 5a.

Hexagonal lattice
In this sub-subsection, we explain the continuation method on hexagonal lattice. Due to the hexagonal lattice, the direction of the shift operator is different from that in the square lattice. Dividing the region into the regular hexagons, we label each mesh as c j , j = 1, . . . , N as in Fig. 1b. In this sub-subsection we use the index of j for the label of each cell instead of i. We write the neighboring cells around the cell c j as c Λ k j , k = 1, . . . , 6, i.e., k j (k = 1, cdots, 6) is the index of neighboring cells for the jth cell.
The typical discrete model on the hexagonal lattice can be described as follows: where f : R 7 → R and g : R → R are intercellular and reaction functions, respectively, here. The linear intercellular interaction f on the cell c i is generally given by Similar to the previous sub-subsection, the typical examples of f are diffusion and Delta-Notch interaction as follows: Utilizing the characteristic function χ c j (x, y), we define the variable at position (x, y) ∈ T 2 at time t > 0 as As the derivation is similar to those in previous sections, we put detailed calculations in "Appendix C". Changing the variable through the characteristic function as in the previous Sect. 2.5.1, we obtain the continuous model: We define the shift operator as Approximating the shift operator by the convolution with the mollifier, we can derive the kernel corresponding to the intercellular interaction on the hexagonal lattice. Figure 5b shows the profile of the kernel for the diffusion f Δ on the hexagonal lattice.

Isotropy of Delta-Notch signaling and radially symmetric kernel
In the previous sub-subsections we explained the continuation method for the uniform lattice with the uniform shape and size. However, the shape of the cells during development is not always uniform. Assuming the uniform lattice for the mathematical modeling for the phenomena might be artificial. Then we propose a kernel of our continuation method conserving the lattice size implicitly without assuming the uniform lattice based on the biological experiments. As explained in Sect. 1, the fly brain is often used for the study of neurogenesis mediated by the Delta-Notch signaling system. The shape of the NEs and neuroblasts (NBs), neural stem-like cells, in the fly brain looks various. Activation of Notch signaling is induced by binding of the Notch receptor with the Delta ligand expressed in adjacent cells at the cell surface. Therefore, activation of Notch signaling might be affected by the shape of the cell membrane. However, as various stochastic noise and other signaling pathway, other than Delta-Notch signaling are involved during development (Tanaka et al. 2018), we conjecture that the shape of the activated region of Notch may become isotropic and averaged. We asked how Notch signaling is activated when Delta is artificially expressed in a small number of cells. In this condition, the Notch activity was visualized by using the NRE-dVenus transgenic construct (Housden et al. 2012). As shown in Fig. 6, Notch signaling was activated in a group of cells immediately adjacent to the Delta-expressing cells through trans-activation forming a concentric distribution pattern. The inactivation of Notch signaling within the Delta expressing cells is most likely due to the effect of cis-inhibition (del Álamo et al. 2011;Sprinzak et al. 2010). This result suggests that the shape of cells do not affect the spatial activation pattern of Notch in vivo.
Based on this experimental results, we propose the following shape for the kernel The profile of this kernel is shown in Fig. 7. The donut-like pattern of Notch activation is consistent with the concentric shape of the kernel used in (30). By using this shape of kernel, the nonlocal operator becomes radially symmetric. This radially symmetric kernel is also applicable to describe the signaling system with projections of cells such as pigment cells in the skin of fishes Watanabe and Kondo

Applications
In this section we will apply our continuation method for some discrete models of previous studies, and perform the numerical simulations to investigate how patterns are generated.

Continuous model of the proneural wave
The developing fly brain looks like a hemisphere. During early stages of development, undifferentiated NEs proliferate and these NEs differentiate into NBs later stages of development. The transition from NEs to NBs propagates from the medial to the lateral directions. Since the transition is visualized by the transient expression of L'sc, which is one of the AS-C complex member and acts as a proneural factor, the propagation of the differentiation in the fly brain is called the proneural wave (PW) (Yasugi et al. 2008).
where the calculation region is set as From the assumption (A5), using the characteristic function, we change the variable as follows: We propose the following the continuous model: where the profile of the kernel K = K (x, y) is determined by the profile of the lattice. For the simple description, we impose the local term of the EGF in the max function of the fourth equation based on the assumption (A6). We perform the numerical simulation to investigate whether the continuation method for the model (31) is effective or not. First, we perform the numerical simulations by using the kernel corresponding to Delta-Notch interaction on the square a b Fig. 9 The results of the numerical simulations for the state of the differentiation A in the continuous model  Figure 9 shows the numerical results of (32) with the kernel corresponding to square lattice. The uniform propagation of the differentiation corresponding to the PW is replicated with suitable parameters as in Fig. 9a. As the parameter corresponding to the strength of activation of EGF a e is decreased, the strength of the lateral inhibition by Delta-Notch becomes relatively larger. In this situation, the continuous model with the kernel reproduces the nonuniform propagation of the differentiation of AS-C corresponding to the salt-and-pepper pattern as in Fig. 9b. In this numerics, the square salt-and-pepper pattern does not depend on the numerical mesh in the code of numerical simulation. This one square region reflects the size and shape of one cell. We put the numerical results of the original discrete model for the PW reported in Sato et al. (2016), Tanaka et al. (2018) in "Appendix A". Secondary, the numerical results of the continuous model with the kernel corresponding to Delta-Notch interaction on hexagonal lattice are shown in Fig. 10. Similarly to Fig. 9, with the large value of strength of the activation for the EGF a e the continuous model with the kernel corresponding to the hexagonal mesh has replicated the mode of the PW. As decreasing the value of a e , the stripe propagation of the differentiation is reproduced. Furthermore, by decreasing the value of a e more, the nonuniform propagation of the differentiation corresponding to the salt-and-pepper pattern has been reproduced. We can observe that each differentiated region exhibits the hexagonal shape. Even though the prepared mesh in the code of the numerical simulation is square, we can replicate the hexagonal patterns in the continuous model. This numerical results suggest that we can directly introduce the information of the spatial discreteness into the continuous model, and the solution of the continuous model with the suitable kernel can reproduce the solution of discrete model.
Next, we perform the numerics with the radially symmetric kernel (30). As shown in Fig. 11, the mode of uniform PW and stripe propagation have been reproduced a b c Fig. 10 The results of the numerical simulations for the state of the differentiation A in the continuous model (32) with the kernel corresponding to the hexagonal lattice. The parameters are same as that in Fig. (9) and K (x, y) = 6 k=1 (τ Λ k ρ ε ) is used. a PW with a e = 5.0, b Stripe pattern with a e = 2.1, c Salt-and-pepper pattern with a e = 1.0 depending on the value of a e . Moreover, as in Fig. 11c, the continuous model has replicated the propagation of the salt-and-pepper pattern even by using the radially symmetric kernel (30). As shown in Fig. 11c, the profile of the differentiated region is spotted and can be interpreted that it indicates the shape of the averaged cells. Although we also use square mesh in the numerical simulation code in these numerics, the continuous model with radially symmetric kernel can reproduce the uniform and nonuniform propagations. The right side of Fig. 11 shows the section of the numerics of (b) in y = L y /2 at t = 20.0. The blue curve corresponds to the profile of Delta, and we can observe that the Delta is expressed at the wave front. These numerical results are consistent with the observation that Delta expression is localized to the cell membrane in the real fly brain. In the view of mathematical modeling, it is explained that the term of the activation from the AS-C in the front a d A(A 0 − A) is imposed in the third equation of (32). We succeeded in reproducing the realistic pattern through our continuation method.
As mentioned above, the results of our numerics suggest that we can analyze the solution observed in the discrete model in the framework of the continuous model equipping the spatially discretized initial data. In the successive subsections, we perform the numerical simulations of discrete model of the PW on growing domains and expansions of the model on the sphere by using our continuation method to investigate the realistic situation of the developing fly brain.

Modeling of cell division on the discrete model
Owing to our continuation method by the convolution with the kernel, we are able to model the cell division and proliferation in the discrete model. In this subsection, we explain this application by using the model of the PW.
During the process of the PW, the nonuniform cell division occurs on the surface. The fly brain develops via an early NE expansion phase followed by a differentiation phase from NEs to NBs (Egger et al. 2007;Hofbauer and Campos-Ortega 1990). When we try to add the effect of cell division to the discrete model, it is often artificial because we must decide the timing, direction, and number of cell division. However, we can introduce this effect naturally in our continuation method by expressing it as the domain growth. We put the explanations of the basic idea for the domain growth in "Appendix D". Using the method of the domain growth, we add the effect of cell division to (31) as follows: where K (x) = ρ ε (x −l)+ρ ε (x +l),K is the kernel with changed variable of K , Γ y is a derivative bijection function, and η is determined below. The detail ofK and Γ y are explained in "Appendix D". Since AS-C can be regarded as the level of differentiation, we suppose that the cell is the NE if the value of A(x, t) is close to 0, and the cell is the NB if the value of A(x, t) is close to 1. To express the nonuniform cell division of the PW, η = η(y, t, A) is given by the monotone decreasing function with respect to A, because NE is divided on the surface of fly brain. Here, we assume that the point satisfying A(t, x) ≥ A * is not divided on the surface. Now, we set where η 0 is a constant. Figure 12 shows the numerical results of (33) in the cases of fixed domain and the nonuniform cell division, respectively. In the beginning of the numerical simulations, we observed the similar patterns in the both (a) and (b). However, the NBs newly appear between the valleys of regions of NBs in Fig. 12b as the time passes. This numerical result can be explained in the view of mathematical modeling as follows. The differentiation of A at each point is inhibited by N in the max function of the fourth equation in (33). N in a cell is activated by D which is activated at the wave fronts of the regions corresponding to the adjacent cells. Therefore, when the region corresponding to the NEs is close to the wave front, the differentiation is inhibited. Conversely, farther from the wave front the region corresponding to the NEs is, the weaker the strength of the lateral inhibition of Notch becomes. As the EGF diffuses to the region, the differentiation of NBs occurs in the valleys of the regions corresponding to NBs. At present, we consider the continuous model of the PW in one-dimensional space. In the future, we will try to calculate in two-dimensional space and on the sphere surface in order to apply the experiment. Furthermore, in Kawamori et al. (2011), it is reported that the wave front of the PW is dented in the clone of fly brain due to the fast NE's division. Thus, we want to understand the reasons why the profile of the wave front is affected by the speed of the NE's division from the viewpoint of the mathematical model.

Description of discrete model on sphere surface
For another example of applications of our continuation method to mathematical modeling, we explain the description of the discrete model on sphere surface. We show that we can deal with the discrete model on the sphere surface by using the radially symmetric kernel (30).
Various pattern formations often occur in the region of sphere surface in the development of multicellular organisms. In the case of the PW investigated in the previous subsections, the fly brain is the hemisphere-like shape, and the PW sweeps across the surface. It is natural to construct the discrete model for the PW on the sphere surface, but the mathematical studies of the PW have been discussed on the 2D plane due to the technical difficulties of the discreteness in the numerical simulations on the sphere. Here, our continuation method can overcome these difficulties and enables us to deal with the model on the sphere surface as the continuous model equation. In practice, by applying the continuation method with the radially symmetric kernel with a radius r > 0 in Sect. 2.5.3, we can compute the continuous model (32) on the sphere surface by the spectrum method. We put the explanations of the basic idea for the spectrum method on sphere surface in "Appendix E". Using the spectrum method, we can compute the following model of the PW on the sphere surface numerically: where the equation of D and A are same as the equation (31), and r S 2 is a sphere with radius r > 0 and K : [0, 2r ] → R is defined as The Laplace-Beltrami operator Δ r S 2 on the general sphere with the radius r > 0 is given by where the definition of the Laplace-Beltrami operator is in "Appendix E". The convolution operator on the sphere * r S 2 is computed as dΩ r is denoted by the standard measure on r S 2 , and * S 2 is the convolution on the unit sphere. According to this calculation, we can rewrite the equation (35) on the unit sphere as follows: where k in is a positive constant, J = J (x, t) reproduces the profile of the JAK/STAT signaling, which cancels the biological noise in the fly brain reported in Tanaka et al. (2018), and the notation of unknown variables are based on (36). When we calculate the equation (37) by the spectrum method, the spatial noise arises from finite spherical harmonic expansion. Therefore, we need the effect of JAK/STAT as the role of the noise reduction (Tanaka et al. 2018). For simplicity, we assume that the value of J (x, t) is spatially uniform. Furthermore, as the spatial interaction of the equation (37) are the diffusion term in first equation and the convolution term in second equation, we calculate numerically the evolution of E and N by the spectrum method and compute the evolution of D and A by the Euler method. In numerical simulation of Fig. 13 at which parameters are the same as those of Fig. 11, we obtain numerical results of the propagation of AS-C similar to the case of the 2D plane when r = 10.0. We can interpret the reason of this result as follows. When the radius of the sphere r is relatively large compared to the cell size, the curvature of the cell surface becomes small. From this, the cell on the sphere surface can be regarded as the same state as the case of the plane. Therefore, we obtain similar numerical results to those of Fig. 11.
We observe the PW is accelerated as the wave directs from the equator to the north pole as in Fig. 14. We think that this arises from the diffusion of the EGF ligand. Because the space become narrower as the wavefront approaches the pole, the EGF ligand accumulates and induce faster NB differentiation. It is not clear whether the speed of the PW progression is accelerated when the PW reaches close to the pole in vivo. This will be one of the interesting questions to be solved in the future by performing live imaging of the PW and quantitatively measured the speed of the wave progression.

Application to the model of planar cellular polarity
In Ayukawa et al. (2014), a discrete model for planar cellular polarity (in short, PCP) of epithelial hair in the fly wing has been proposed by focusing on the interactions of Fig. 13 The results of the numerical simulations for the state of the differentiation A in the continuous model on the sphere surface (37). The parameters are same as that in Fig. 11, and r = 10.0 and k in J ≡ 1.0×10 −3 . a PW with a e = 2.0. b Stripe pattern with a e = 0.7. c Salt-and-Pepper pattern with a e = 0.4 Fig. 14 The velocity of PW for each position, when d t = 0 and other parameters are the same as that in Fig. 11. φ represents the latitude transmembrane proteins, distal complexes and proximal complexes. The intercellular protein and the cytoplasmic component are asymmetry localized by the intercellular interactions. Due to the asymmetric localization of the proteins, the direction of an epithelial hair in a cell is determined. If the transmembrane receptor, Frizzled (Fz) is localized in a cell, the other transmembrane protein, Strabismus (Stbm), is localized in the opposite side of the cell membrane in the same cell. Each membrane protein interacts distal and proximal complexes in a cell, which causes the localization of the polarization of Fz and Stbm in the neighboring cells. Fz and Stbm in the neighboring cells are localized near sides. Localization of these proteins between adjacent cells leads to the local alignment of PCP among small group of cells.
A simple mathematical model succeeded in describing the mechanism of the PCP Ayukawa et al. (2014). In the modeling of Ayukawa et al. (2014), the region corresponding to the fly wing is divided into hexagonal lattices and each lattice is labeled as c i , i = 1, . . . , N . By denoting the unknown variable for the direction of Fz in the ith cell c i by θ i = θ i (t) the following discrete model is proposed by the authors of Ayukawa et al. (2014): where the number of the term in the linear combination with sin function depends on the spatial dimension and arrangement of the lattice. We perform our continuation method to this discrete model for PCP. Setting as we have the following continuous model by the changing the variables as in Sect. 2.5.2, where the shift operator is defined by (28). This model is equivalent to (38) if the initial data is equal. The form of the shift operator is changed depending on the shape of the set lattice. We performed numerical simulations with the discretized initial datum on square lattice.
As in Fig. 15, we observe that the continuation model (39) with the discretized initial data can replicate the patterns in the discrete model (38). In this simulation, the color corresponding to the direction of an epithelial hair gets gradually uniform. This solution corresponds to that all epithelial hair grow in the same direction. On the other hand, as in Fig. 16, the solution corresponding to a swirl of the epithelial hairs is obtained in the steady-state as the number of cells is increased. This generation of the pole of θ is consistent to the report by Ayukawa et al. (2014).
In the discrete model, the number of the unknown variables is required for the same number of the cells. Accordingly, it is sometimes hard to compute the analytic calculations, for example, eigenvalue problem in the linear instability around the equilibrium solution. However, our method can reduce the discrete model with multiple components into scalar continuous model, and it gives us simpler calculations. We perform the linear instability of the model (38)  one-dimensional space. Suppose that the number of cells is N , and that the length is l. Setting the region as T = [0, Nl], we impose the periodic boundary condition. The model of PCP in one-dimensional space is given by the following form: For this interval T, we found that one equilibrium solution is given by θ Additionally, the constant solution θ * (x) = α, α ∈ [0, 2π ] is also equilibrium solution of (40). The linear stability analysis around the equilibrium solutions is explained as follows. Letting the range of the linear operator be in R, and setting the perturbation as θ(x, t) = θ * + ε(x, t) and substituting the model (40) by it, we have Linearizing this problem around equilibrium solutions, we have the following eigenvalue problem: where ϕ = ϕ(x) is the eigenfunction associated by the eigenvalue λ. Plugging the nth term of the Fourier series expansion ϕ n = a n exp − 2nπi Nl x , a n := 1 Nl to (41), where i is imaginary number, we obtain that We have the eigenvalues λ n = 2 cos 2π N cos 2nπ N − 1 . The calculation of the eigenvalue of f Δ in the matrix form is also written in "Appendix B", and each result is consistent. From this calculation, if the number of cell N is bigger than 3, we see that equilibrium solutions is linearly stable. By the same calculus, it is shown that the constant solution θ * is also linearly stable. Even in the two-dimensional case, our above method enables us to compute the linear stability analysis if we have the equilibrium solutions.

Discussion
In this paper, we proposed a continuation method for discrete models, using shift and convolution operators while conserving the size and shape of cells and lattices. The proposed method enables the conversion of nonlinear discrete models into continuous models in a systematic manner, retaining the discreteness information. As the continuous model applied to our method with the shift operator is point-wisely equivalent to the original discrete model, the solutions are equal if the initial data are the same. The framework of the continuous model provides a few advantages, as per the analysis results. As in Sect. 3.4, we reproduced the pattern for the PCP in epithelial hair corresponding to that of the discrete model. Furthermore, we constructed the equilibrium solutions and performed the linear stability analysis in the continuous PCP model, using the Fourier series expansion. In Sect. 2.3, we showed that our continuation method can be applied to the discrete models on nonuniform lattices. Although the framework of discrete model is technically difficult to express the dynamics on nonuniform lattices, The proposed method enables us to treat the spatial non-uniformity on the continuous models mathematically. In the future, we will extend our work to the analyses and applications to this continuous models from the discrete models on the nonuniform lattices.
As a next step, we reduced the continuous model with the shift operators to the nonlocal evolution equation, using the approximation of the shift operator by convolution of a mollifier. We have also conducted the singular limit analysis of the discrete model and the nonlocal evolution equation in a one-dimensional interval with periodic boundary condition, showing that every solution is sufficiently close in L 2 (T) space. This suggests that nonlocal evolution with suitable kernels is capable of approximating the solution of discrete models, if the initial data are the same. Using the nonlocal evolution equation with the kernel of mollifiers, we have succeeded in replicating the pattern observed in the original discrete PW model. When the intercellular interaction was linear, the profile of kernel was determined by the lattice shape as shown in (26) and (29). Using these kernels in the continuous model for the PW, we reproduced the square and hexagonal shapes of the salt-and-pepper patterns.
Furthermore, we experimentally confirmed the isotropy of the Delta-Notch signaling system in the real fly brain. Based on the biological experiment, we proposed a radially symmetric kernel for the domain comprising averaged cell shapes. Even with the radially symmetric kernel for the Delta-Notch signaling interaction, we could reproduce the various propagation patterns of the PW. The radially symmetric kernel can also be applied to the discrete models on nonuniform lattices if the molecular and cellular system is not affected by the shape of cells and lattice as explained in Sect. 2.5.3. For application to the biological experiments, using a kernel with a small width, such as the Friedrichs mollifier, yields results that are more compatible with experiments than the combination of the shift operators. However, the shift operator can be more convenient for the analysis. In Ei et al. (submitted), by arranging the Dirac Delta function radially as the kernel in the continuation method, the reduction method of the continuous model for PW into 1 or 2 variable system, and its numerical simulations are addressed. Nonlocal evolution equation with certain kernel is sometimes difficult to analyze. However, nonlocal evolution equation provides us with a unifying concept to mathematical modeling and analysis as it is capable of including partial differential equations such as the reaction diffusion systems and discrete models. Moreover, various analytical and numerical methods for partial differential equations are applicable to nonlocal evolution equations as explained in Sects. 3.2, 3.3 and 3.4. In future, we also intend to extend our work to other domains and include higher dimensions of singular limit analysis.
In Sect. 3.2, 3.3 and 3.4, we applied the continuous model on the PW progression and the PCP formation. We demonstrated that the continuous model can easily include the effect of NE cell proliferation and can expand the simulation result from the 2D plane to the 3D spherical surface. In fly brain development, the final numbers of NBs and neurons are dependent on the NE proliferation. Using biological experiments, it has been shown that NE expansion is regulated by several signaling pathways. The PI3K/Akt/TOR pathway promotes NE proliferation in a diet-dependent manner in the early stages of development (Franco and Carmena 2019;Lanet et al. 2013). The Hippo pathway inhibits the overgrowth of NEs and inactivation of Hippo signaling inhibits NB differentiation (Kawamori et al. 2011;Reddy et al. 2010;Richter et al. 2011). By combining biological experiments on these signaling pathways with numerical calculation, it will be possible to understand the in vivo situation of the PW progression in more detail. In this paper, we presented two biological examples, for which our numerical method is applicable. Here, we emphasize that the numerical method is also useful for other biological systems because our method is based on fundamental and conserved intercellular interactions. The continuation method and numerical calculation with continuous models will facilitate our understanding of a wide variety of biological processes, both quantitatively and qualitatively. The propagation of the salt-and-pepper pattern, and stripe pattern are observed in the numerical simulations. We see that the profile of A i, j are uniform in each cell c i, j .

Energy estimate
We estimate error between the solutions of (P S ) and (P ε ) with (11) in the singular limit analysis with linear function f given by (10). We assume the initial condition of (P S ) and (P ε ) is the same, and given by u( First, we compute the fundamental solution of (P S ) without g (u). As the equation (P S ) and the equation (P D ) is equivalent, we solve the equation (P D ). The system (P D ) is described by where V = (u 1 (t), . . . , u N (t)), F : R N → R N and , (i f N is even).
Since the matrix A is cyclic, we obtain the eigenvalue (12) i.e., Associated eigenvector V k of which components are {v n } N n=1 with eigenvalue λ k is computed by v n = ω nk , (k = 1, . . . , N ). Thus, the fundamental solution of (42) is given by Therefore, we obtain the exact solution as Setting a vectorṼ k as {ṽ n } = ω −nk andλ k = [(N −1)/2] j=−[(N −1)/2] a j ω − jk , we find that AṼ k =λ kṼk . Thus, if a j = a − j , ( j = 1, . . . , . From this calculation, setting a vector W k as {w n } = cos(2π nk/N ) + sin(2π nk/N ). Then AW k = λ k W k . Thus the fundamenal solution of (42) is given by where the matrix (W 1 , W 2 , . . . , W N ) t is the transpose of the matrix (W 1 , W 2 , . . . , W N ). Next, we show the existence the global solution of (P S ) and (P ε ) under the above assumptions, respectively.

Proof of Proposition 2
The existence and uniqueness of the mild solution for (P S ) in C([0, T ], L ∞ (T)) is guaranteed by the fixed point theorem. We show the global existence in L ∞ (T) by the argument of the maximum principle. For a contradiction, we assume that there exists a positive finite constant T > 0 such that lim sup τ T max 0≤t≤τ u(·, t) L ∞ (T) = ∞ Then we take a sequence {T n } n∈N , T n T as n → ∞ such that as n → ∞. Indeed, for any R there exists a positive constant n 0 ∈ N such that for all n ≥ n 0 , we see that R n > R, and for any 0 < t < T n , there exist j n ∈ N such that u(·, t) L ∞ (T) = |u j n (t)|. We define the points (x n , t n ) ∈ T × (0, T n ) as satisfying |u(x n , t n )| = |u j n (t n )| = R n . As R n → ∞, we can choose Multiplying the principal equation of (P S ) by u, then a j u(x n + jl, t n )u(x n , t n ) + g(u(x n , t n ))u(x n , t n ), from (A2) and considering the degree of polynomial of R n . This yields a contradiction.
Similar to the proof of Proposition 2, we show the global existence of the equation (P ε ).

Proof of Proposition 3
The existence and uniqueness of the mild solution for (P ε ) in C([0, T ], L ∞ (T)) is guaranteed by the fixed point theorem. By the same argument as the proof of Proposition 2, we obtain the L ∞ (T) estimate for the solution of (P ε ). For a contradiction, we assume that there exists a positive finite constant T > 0 such that Then we take a sequence {T n } n∈N , T n T as n → ∞ such that as n → ∞. Indeed, for any R there exists a positive constant n 0 ∈ N such that for all n ≥ n 0 , we have R n > R, and for any 0 < t < T n , we define the points satisfying |u ε (x n , t n )| = R n by (x n , t n ) ∈ T × (0, T n ). In the case that the candidates of the maximum point is on the discontinuous point, employing the larger value at the point by taking the left-sided and right-sided limits, we define the point as (x n , t n ). As R n → ∞, we can choose Moreover, for any {a j } N j=1 , and positive constants p, g 0 , g 1 and g 2 , there exist sufficient large constant r , we see that K L 1 (T) r 2 + g(r )r ≤ ⎛ ⎜ ⎜ ⎝ N −1 2 j=− N −1 2 |a j | + g 2 ⎞ ⎟ ⎟ ⎠ |r | 2 − g 0 |r | p+1 + g 1 |r | 3 < 0 from (A2) and considering the degree of polynomial of r . Thus, there exists a positive constant n 1 ∈ N such that for n ≥ n 1 0 ≤ 1 2 replacing r with u ε (x n , t n ). This yields a contradiction. Now we show the singular limit analysis between the solutions of (P S ) and (P ε ).
Proof (Proof of Theorem 1) Taking the difference between the solutions of (P S ) and (P ε ), we have Multiplying the above equation by U ε and integrating it with respect to x, we obtain that from (A3) and where θ ∈ (0, 1), and C 7 = (g 3 sup t>0 u + θ u ε L ∞ (T) + g 4 ).
Utilizing the classical Gronwall lemma yields that

Derivation on square and hexagonal lattices
In this appendix we put the detailed calculations for the derivations of the continuous models in square and hexagonal lattices. For (25), we can calculate as N i, j=1 u i− p, j−q χ c i, j (x, y) = u(x − pl x , y − ql y , t), ( p, q ∈ {1, . . . , N }).
In the case of hexagonal lattice, the derivation is given as follows. Here for simple description, introducing the complex variable x = x + yi ∈ C, we identify the two-dimensional Euclid space R 2 with the complex plane C. Then we compute that N j=1 u Λ k j (t)χ c j (x, y) = u x + le i( π 2 − (k−1)π 3 ) , t = u x + cos π 2 − π(k − 1) 3 l, y + sin π 2 − π(k − 1) 3 l, t .
Using the shift operator (28) and approximating the shift operator by the convolution with the mollifier, we can derive the kernel corresponding to the intercellular interaction on the hexagonal lattice.

Reaction, diffusion and nonlocal interaction on growing domain
We explain the notion of the mathematical model on a growing domain by using a reaction diffusion equation with nonlocal interactions. Based on the previous reports (Crampin et al. 1999(Crampin et al. , 2002, general scalar reaction diffusion equation with convolu-tion term on a growing domain in one-dimensional space is given by ∂u ∂t + ∇(a · u) = Δu + K * u + f (u), in (0, L(t)) × {t > 0}, where u = u(t, x) is an unknown function, f : R → R is a nonlinear function, a = a(x, t) is the velocity field of the flow and satisfies dx dt = a(x, t), x ∈ (0, L(t)).