Bifurcation and Periodic Triangular Pattern Formation in Reaction-Diffusion with Anisotropic Diffusion

Abstract

Turing demonstrated a coupled reaction-diffusion equation with two components produced steady-state heterogeneous spatial patterns, under certain conditions. The instability found by Turing is now called a diffusion-driven instability or Turing instability. Systems in two dimensions produce spot and stripe patterns, and these systems have been applied as models to explain patterns observed in biological and chemical fields and to develop image information processing tools. Previously, we developed a method that utilizes a reaction-diffusion system with anisotropic diffusion that exhibits triangular patterns, thereby introducing a certain anisotropic strength. In this chapter, we discuss the effects of anisotropic diffusion on the generation of triangular patterns. By defining the statistical index characterizing the spatial patterns, we investigated the parameter range over which the triangular patterns were obtained. We determined the explanatory variable based on the relative distance of the pitchfork bifurcation point between the maximum and minimum anisotropic diffusion functions. Its relevance to diffusion instability is also discussed.

1. Introduction

The formation of organisms and developmental phenomena are natural self-constructed phenomena. Mathematical and physical methods play an important role in the analysis of self-organizing phenomena, such as the collapse of sandy mountains [1], wind ripples on sand dunes, and snow crystal patterns [1, 2, 3].

For a number of pattern-formation phenomena, the so-called reaction-diffusion (RD) equation has been applied to model self-organizing patterns [1, 3]. In 1952, Alan Turing, a British mathematician well-known for proposing the basic principles of cryptanalysis and computers, proposed one of the most famous explanations for self-organizing phenomena by postulating a mechanism based on the reaction and diffusion of the materials [4]. This mechanism is now known as diffusion-driven instability or Turing instability [1, 2, 3].

Some researchers observed that Turing patterns can be realized in chemical reaction systems [1, 2, 3] and animal skin patterns [5, 6]. These patterns have two dimensions, and stripe, spot, and labyrinthine patterns can be realized in the RD equations [1].

Furthermore, several researchers have developed methods for applying the RD model to image information processing. Admatzky et al. [7] and Ito et al. [8] developed an algorithm for segmenting 2D periodic stripes or spotted images using an RD model. Applying such a mechanism has various advantages; for example, they do not require for noise removal mechanisms.

We applied the image-processing algorithm for pattern formation of outer hair cells in the ear [9, 10] that arranges periodic triangles on the surface of the organ of Corti in the inner ear [11, 12]. To quantify the arrangement of the periodic triangules, we developed an RD model with anisotropic diffusion. In Ref. [10], we explored the conditions for emerging periodic patterns. We found the following explanatory variables: on the heuristic argument of the relative difference in the bifurcation points between the maximum and minimum values of the anisotropic diffusion coefficient. This chapter presents the recent research on periodic triangles as Turing patterns.

2. Model

2.1 Turing instability

Turing [4] proposed the following concept as one of the mathematical mechanisms. The system is expressed as follows:

$$\frac{\partial u}{\partial t}=D_u{\nabla }^{2}u+f\left(u,v\right),\mathit{\text{and}}\frac{\partial v}{\partial t}=D_v{\nabla }^{2}v+g\left(u,v\right),$$E1

where u and v denote the densities of diffusive and reactive materials, respectively. The first terms on the right-hand side of Eq. (1) are diffusion terms. Constants $D_u$ and $D_v$ are the diffusion coefficients and are all positive. The polynomials f (u, v) and g (u, v) for u and v are called reaction terms. This form of the equation is called the reaction-diffusion (RD) equation. Turing [4] showed that Eqs. (1) and (2) can stably have spatially non-uniform periodic solutions if the ratio of the diffusion coefficients is sufficiently large and the reaction terms satisfy the following conditions [1, 4].

FitzHugh-Nagumo [13] chose the following reaction terms:

$$f\left(u,v\right)=u-u^3-v,\mathit{\text{and}}g\left(u,v\right)=\gamma \left(u-\mathit{\alpha v}-\beta \right).$$E2

where the constants $\alpha ,\beta$, and $\gamma$ are positive. u has a nonlinear reaction term $u^3$, where the presence of u increases v, and v decreases u. Therefore, u is called the activator and v is called the inhibitor.

There are several other choices for reaction terms, such as the Schnaclenberg model [14]. In the following sections, we study cases using the FitzHugh-Nagumo equation [13] given by Eq. (2). However, later, we discuss another choice of the reaction term.

First, we consider the case in which Eq. (1) has only one stable time-independent uniform solution $\left(u_0,v_0\right)$, which satisfies $f\left(u_0,v_0\right)=0$ and $g\left(u_0,v_0\right)=0$. As the uniform solution is stable, we have $f_u+g_v<0$ and $f_u{g}_v-f_v{g}_u>0$, where $f_u,f_v,g_u$, and $g_v$, the partial derivatives, are evaluated at $\left(u_0,v_0\right)$. By introducing the $u=u_0+c_{1}exp\left(\mathit{\lambda t}+\mathit{iqx}\right)$, and $v=v_0+c_{2}exp\left(\mathit{\lambda t}+\mathit{iqx}\right)$ into Eq. (1), and leaving the first-order with respect to $c_1$ and $c_2$, we obtain

$$\lambda \left(\begin{array}{c}{c}_1\\ c_2\end{array}\right)=\left(\begin{array}{cc}-q^2{D}_u+f_u& f_v\\ g_u& -q^2{D}_v+g_v\end{array}\right)\left(\begin{array}{c}{c}_1\\ c_2\end{array}\right)$$E3

where the wavenumber q represents the space modulation. To obtain a solution for $c_1\ne 0,c_2\ne 0$, we obtain

$$\begin{array}{l}{\lambda }^2+\lambda \left\{q^2\left(D_u+D_v\right)-f_u-g_v\right\}+N\left(q\right)=0,\\ N\left(q\right)=D_u{D}_v{q}^4-q^2\left(D_u{g}_v+D_v{f}_u\right)+f_u{g}_v-f_v{g}_u.\end{array}$$E4

The condition for one of the eigenvalues, $\lambda$, in Eq. (4) is positive. In other words, the sufficient condition for instability is $N\left(q\right)<0.N\left(q\right)$ becomes negative at finite q as $D_u$ decreases. Differentiation with respect to $q^2$ shows that

$$N_{\mathit{min}}=-\frac{{\left(D_u{g}_v+D_v{f}_u\right)}^2}{4D_u{D}_v}+f_u{g}_v-f_v{g}_u,\mathit{\text{where}}{q}^2=q_c^2=\frac{D_u{g}_v+D_v{f}_u}{2D_u{D}_v}.$$E5

This indicates that $N_{\mathit{min}}$ is negative when $D_u$ is sufficiently small. Therefore, when $D_u\ll D_v$, the uniform solution is unstable for modulations with finite wave numbers. This phenomenon is known as the diffusion-induced instability or Turing instability [1, 2, 3, 4].

At bifurcation, when $N_{\mathit{min}}=0$ and with fixed parameters except $D_u$, this defines a critical coefficients $D_u^c$, which is the pitchfork bifurcation point. The pitchfork bifurcation point $D_u^c$ is obtained as the appropriate root of

$$g_v^2{\left(D_u^c\right)}^2+2\left(D_v{f}_u{g}_v-2D_v{f}_u{g}_v+2D_v{f}_v{g}_u\right)D_u^c+D_v^2{f}_u^2=0.$$E6

Therefore, in the case of $D_u<D_u^c$, Turing instability occurs [1, 2, 10].

2.2 Turing patterns

The linear analysis performed in the previous section indicated that instability occured. However, an analysis that incorporates nonlinearities is necessary to investigate the time evolution in the region of $D_u<D_u^c$ [1, 2, 15, 16]. We performed numerical simulations for the coupled set of Eq. (1) in two dimensions. Numerical simulations were mainly performed for $D_u=5.00\times {10}^{-4},D_v=5.00\times {10}^{-2},\alpha =0.500,$ and $\gamma =26.0$, by varying $\beta$. Periodic boundary conditions were imposed at the system boundary in a square domain of size: $1.28\times 1.28$ (grid: $128\times 128$). We applied a simple explicit scheme and chose the parameters to satisfy the conditions for Turing instability. The initial conditions of u and v were given as the uniform steady state $\left(u_0,v_0\right)$ plus small random deviations. We performed a numerical simulation for $t=1000$, considering the sufficiently long time to reach the final patterns.

Starting from the random initial distributions shown in Figure 1(a) and (e), the time evolutions of the patterns in a density plot of u over time are shown in Figure 1(b-d,f-h). Figure 1(j-l) displays the spatial variations of the concentrations u and v for the distributions along the arrow in Figure 1(e-h). As shown in Figure 1(i-k), the exponential growth of a specific unstable critical mode revealed by the linear analysis emerges from random initial states in the early stages of pattern formation as explained in Ref. [1]. However, as shown in Figure 1(b-d)or(f-h), the static spot or stripe patterns were eventually obtained. This is because the pattern converges asymptotically through nonlinear responses in the next stage of pattern formation [2, 3, 16]. Here, we note that the distributions of v follow the same periodic patterns with the same periodicities but with different amplitudes, as shown in Figure 1(j-i).

Figure 1.

Time evolutions of $u\left(\overrightarrow{r},t\right)$ and their spatial variation of u and v. (a)-(d) were obtained numerically from Eq. (1) with $\beta =-0.08$; (e)-(h) were obtained numerically from Eq. (1) with $\beta =0.01$; (i)-(l) were the distributions of u (red line) and v (blue line) along the solid arrow in (e)-(h), respectively.

Various discussions have been conducted on pattern selection [1, 2, 3, 15, 16]. Some researchers have discussed this from the viewpoint of variational methods (energy) [2, 3, 16]. Others have discussed the transformation of the solution structure on the central manifold [2, 17].

Notably, Eq. (1) is one of the partial differential equations. This system is not immune to the dependence of the obtained pattern on the initial distribution. In other words, a slight modulation in the initial distribution can affect the directionality of the stripes in the resulting stripe pattern or the arrangement of the spots [1, 2, 3].

3. Triangular patterns obtained by RD model with anisotropic diffusion

3.1 RD model with anisotropic diffusion

Kondo and Miura [6] revealed that some of developmental processes are affected by RD systems, as in Eq. (1). It is also important to determine the direction in which organs and tissues self-organize during development [18, 19].

The development of V-shaped bundles for sensorineural hearing in the inner ear periodic triangular pattern area has self-organized in a specific direction based on the surrounding conditions [11, 12]. For the morphometric measurements of V-shaped bundles, we developed a system that applied a self-organizing periodic triangular system using an RD model with the anisotropic diffusion [9]. As the simple diffusion function model, we applied the periodic function depending on the diffusive directions with adjacent reaction cells [9].

In the study by Shoji and Iwamoto [9], we focused on the case of strong enough anisotropy to obtain periodic triangular patterns. Here, we considered the general case as follows

$$\frac{\partial u}{\partial t}=\nabla \cdot \left(D_u\left(\theta \right)\nabla u\right)+u-u^3-v,\mathit{\text{and}}\frac{\partial v}{\partial t}=D_v{\nabla }^{2}v+\gamma \left(u-\mathit{\alpha v}-\beta \right).$$E7

The diffusion coefficient of u is expressed as

$$D_u\left(\theta \right)=\frac{1}{\sqrt{1-\delta cos3\theta }},$$E8

where $\theta$ indicates the angle of the gradient vector of u. This is expressed as:

$$\theta =arctan\left(\frac{\partial u}{\partial y}/\frac{\partial u}{\partial x}\right).$$E9

We would like to note that the flux is proportional to the gradient vector, however the multiplication factor is determined by the angle of the vector. Eq. (8) implies that the diffusivity of u is the highest for $\theta =0,\frac{2\pi }{3}$, and $\frac{4\pi }{3}.\delta$ denotes the magnitude of the anisotropy of u. This condition satisfies $0\le \delta \le 1$. The case $\delta =0$ implies the isotropic diffusion.

This form of anisotropic diffusion was applied by Kobayashi [20] for self-organizing dendritic crystal. In previous studies [19, 21], we also applied this form of anisotropy to explain the directionality of stripes on fish skin.

Figure 2 shows the time evolution of the two-dimensional distribution of u. The model given by Eqs. (7)–(9) was numerically calculated by changing anisotropy $\delta$. The following parameters were selected: $\alpha =0.50,\beta =0.01,\gamma =26.0$.

Figure 2.

Time evolutions of $u\left(\overrightarrow{r},t\right)$. (a)-(d): A stripe pattern was obtained numerically from Eqs. (7)–(9) with $\delta =0.00$ (without anisotropy); (e)-(h): A stripe pattern emerged from an initially periodic triangular pattern was obtained with $\delta =0.50$; (i)-(l) broken triangles from an initially periodic triangle pattern were obtained with $\delta =0.75$, and (m)-(p) a periodic triangular pattern was obtained with $\delta =0.90$.

We start with the random initial distribution shown in Figure 2(a),(e),(i), and (m). The time evolutions of the patterns in a density plot of u over time are shown in Figure 2(b-d),(f-h),(j-l), and (n-p). If we introduce anisotropy $\left(\delta >0\right)$, the modes for the triangular pattern initially emerge (see Figure 2(f),(j), and (n)). However, the triangular mode could not be maintained in the final pattern if anisotropy was not sufficiently strong (Figure 2(g-h), and (k-l)). Eventually, the domain converges to an ordered periodic triangular pattern, as shown in Figure 2(p), using model (7)-(9) with strong anisotropy $\left(\delta =0.90\right)$.

3.2 Image process

To examine the structure of the triangular pattern shown in Figure 3(a), the numerical data of u were Fourier transformed as follows: Its Fourier transform is

Figure 3.

(a) The triangular pattern obtained by the numerical simulation. (b) Relative peak positions of Bragg spots. (c), (d) peak position of Bragg spots for triangular and spot patterns. (e) the intensity of Bragg spots of (c) and (d). The blue and red lines show the case of the spot and triangular patterns, respectively.

$$\widehat{u}\left(\overrightarrow{q}\right)=\int d\overrightarrow{r}u\left(\overrightarrow{r}\right)exp\left(i\overrightarrow{q}\cdot \overrightarrow{r}\right).$$E10

By displaying the intensity ${\left|\widehat{u}\left(\overrightarrow{q}\right)\right|}^2$ in wavenumber space, the characteristics of the structure emerged [22].

A Fourier analysis are performed for the obtained distribution of $u\left(\overrightarrow{r}\right)$ in the triangular pattern shown in Figure 3(a). For comparison, we also performed a Fourie analysis for the spot pattern shown in Figure 3(b), which was obtained from Eqs. (7)–(9) using the same parameters as those for the obtained stripes $\left(\delta =0.00\right)$, except for $\beta =0.08$. The peak positions centered at the origin were plotted, and the peak intensity was expressed in terms of the relative magnitude of the balls. The patterns in Figure 3 (c) and (d) show that they consist of 12 and 3 Fourier modes, respectively.

Figure 3(e) shows the Bragg peak intensities for $q=\left|\overrightarrow{q}\right|$. The blue and red lines in Figure 3(e) represent the case of the spot and triangular patterns, respectively. The intensities of the higher peaks $\left(q>1\right)$ in the triangular pattern (red line) were larger than those of the spot pattern (blue line). However, the triangular pattern’s first peak $\left(q\approx 0.6\right)$ was almost identical to that of the spot pattern.

Using these features, we introduced the triangular clearness TC to determine the periodic triangular patterns as follows:

$$\mathit{TC}=\left(p_2^2+p_3^2\right)/p_1^2,$$E11

where $p_1,p_2$, and $p_3$ are the intensity of the first peak $(q\approx 0.6$ in the case of Figure 3(e)), and the second peak $(q\approx 1.05$ in the case of Figure 3(e)), and the third peak $(q\approx 1.25$ in the case of Figure 3(e)), respectively. We therefore considered the pattern to be composed of a distinct triangle when the TC was larger. When the TC was small, the pattern was considered to be striped or to have broken triangles and stripes.

4. Conditions for triangular patterns obtained

4.1 Comparison between TC and $\delta$

Using the index TC proposed in the previous section, we examined the obtained patterns. We used the same parameter set used as in the computer simulation shown in Figure 1, except $\delta ,D_u$, and $D_v$. For each set of $D_u$ and $D_v$, we performed a numerical simulation of Eqs. (7)–(9) by varying $\delta$. We calculated TC for the obtained patterns in $\delta$. Figure 4(a-c) show the changes in the index TC depending on anisotropy $\delta$.

Figure 4.

Triangular clearness: TC versus the anisotropy $\delta$. (a)-(c), and the index A in Eq. (12) (d)-(f). The parameters were as follows: $\alpha$=0.50, $\beta$=0.01, and $\gamma$=26.0, which were also used in Ref. [9], as well as (a) and (d) $D_u=5.00\times {10}^{-4},D_v=5.00\times {10}^{-2}$, (b) and (e) $D_u=6.00\times {10}^{-4},D_v=5.00\times {10}^{-2}$, and (c) and (f) $D_u=5.00\times {10}^{-4},D_v=7.50\times {10}^{-2}$.

As $\delta$ is increased, a gradually cleaner triangular pattern is obtained. However, the changes in TC are different for the parameter sets shown in Figure 4(a-c). It is difficult to determine the range of the obtained periodic triangular patterns.

4.2 Heuristic index for triangular

Based on a heuristic mathematical argument, we developed the mathematical relations of the obtained periodic triangular patterns. The relations was provided by focusing on the range of the values of the diffusion coefficient function $D_u\left(\theta \right)$, which ranges from the maximum $D_u^{\mathit{max}}=D_u/\sqrt{1-\delta }$ and minimum $D_u^{\mathit{min}}=D_u/\sqrt{1+\delta }$.

We show the results when the parameters $D_u=5.00\times {10}^{-4},D_v=5.00\times {10}^{-2},$ and the pitchfork bifurcation $D_u^c=6.60\times {10}^{-4}$ is obtained using Eq. (6), which was the case shown in Figure 2. The relative ranges of the values of the coefficient function (8) are shown in Figure 5.

Figure 5.

Schematic image for the ranges of Eq. (8) depending on anistropy $\delta$.

When diffusion anisotropy is not introduced $\left(\delta =0.00\right)$, the diffusion coefficient function $D_u\left(\theta \right)$ is always same, indicated by a single point in the area where Turing instability occurs, as shown in Figure 5. In contrast, when diffusion anisotropy is introduced, the diffusion coefficient function $D_u\left(\theta \right)$ takes between $D_u^{\mathit{max}}$ and $D_u^{\mathit{min}}$. It is also seen that the range was in the region where the equilibrium was stable and the region where Turing instability occurs beyond the pitchfork bifurcation, $D_u^c$, as given by Eq. (6). Moreover, as the anisotropy $\delta$ was increased, the range of $D_u$ spreads widely over the area of equilibrium stable. Here, we noted that Turing instability occured when $D_u$ decreased in fixed $D_v$. We quantified the these difference by the ratio of the distance between $D_u^{\mathit{max}}$ and $D_u^c$ and between $D_u^c$ and $D_u^{\mathit{min}}$ as follows:

$$A=\left(D_u^c-D_u^{\mathit{min}}\right)/\left(D_u^{\mathit{max}}-D_u^c\right).$$E12

Note that A satisfies $A\ge -1$. When $\delta =0.00$ (no anisotropy), A is −1.0, whereas A increased with anisotropy $\delta$. In the case of $\delta =0.50,D_u=5.00\times {10}^{-4},$ and $D_v=5.00\times {10}^{-2}$, a triangular pattern appeared first, but eventually a stripe pattern formed as shown in Figure 2(e-h) for $A=0.197$. When $\delta =0.75$, a triangular pattern replaced a broken triangles, as shown in Figure 2(i-l) for A = 1.22. When $\delta =0.90$, a triangular pattern was sustained, as shown in Figure 2(m-p) for $A=3.13$.

Figure 4(d-f) summarizes the results for TC vs. A with varying anisotropy $\delta$ for each set of $D_u$, and $D_v$.

4.3 Regression analysis

Figure 4(d-f) shows the changes in TC with respect to A. A rapid change from a low to high value of TC is observed. Although the positions of the change points are different in Figure 4(d-f), the high and low values of TC are identical, and the distributions are all S-shaped. These functions known as logistic curves are often applied [23, 24].

Regression analysis was performed using the following function:

$$\mathit{TC}\sim \frac{a}{\left(1+{\mathit{be}}^{\mathit{rA}}\right)},$$E13

where a, b, and r are estimated parameters. In Figure 4(d-f), the regression function is represented by a solid line. The fitted parameters are listed in Table 1.

Figure	Parameters in (13)			$Q_0$	$Q_2$
	a	b	r	$\left(A_0,{\mathit{TC}}_0\right)$	$\left(A_2,{\mathit{TC}}_2\right)$
Figure 4(d)	$2.89\times {10}^{-2}$	7.92	−1.45	$\left(\mathrm{1.25,1.45}\times {10}^{-2}\right)$	$\left(\mathbf{2.34},2.28\times {10}^{-2}\right)$
Figure 4(e)	$2.95\times {10}^{-2}$	11.0	−1.48	$\left(\mathrm{1.61,1.48}\times {10}^{-2}\right)$	$\left(\mathbf{2.49},2.23\times {10}^{-2}\right)$
Figure 4(f)	$2.81\times {10}^{-2}$	2.05	−0.93	$\left(\mathrm{0.72,1.47}\times {10}^{-2}\right)$	$\left(\mathbf{2.19},2.32\times {10}^{-2}\right)$

Table 1.

Parameters fitted to the Eqs. (13)–(15) for the obtained patterns obtained using Eqs. (7)–(9).

The coordinates of the inflection point of the function (13), $Q_0\left(A_0,{\mathit{TC}}_0\right)$, are $A_0=\frac{lnb}{r}$, and ${\mathit{TC}}_0=\frac{a}{2}$. At this point, this function reaches 50% of its maximal value a. Moreover, the inflection point coordinates of the derivative functions of (13), $Q_1\left(A_1,{\mathit{TC}}_1\right)$, and $Q_2\left(A_2,{\mathit{TC}}_2\right)$ can be expressed as

$$A_1=A_0-1.318/r,\mathit{\text{and}}{A}_2=A_0+1.318/r.$$E14

The TC values at these points were calculated as follows:

$${\mathit{TC}}_1=\frac{a\left(1-\sqrt{1/3}\right)}{2}\sim 0.211a,\mathit{\text{and}}{\mathit{TC}}_2=\frac{a\left(1+\sqrt{1/3}\right)}{2}\sim 0.789a.$$E15

Therefore, these points were subjected to changes until they reached approximately 21% and 79% for the final TC values of $Q_1$ and $Q_2$ [23, 24]. The coordinates of $Q_0$ and$Q_2$ are shown in Table 1. It was found that the parameters for a and $Q_2$ were similar in all the three cases. However, $Q_0$ and the parameters in Eq. (13) were different.

We also performed a trial in which the initial distribution was given in a different manner; however, as mentioned above, we obtained the same values for the parameters a and $A_2$ (see Ref [10]).

4.4 Range of $\delta$

As seen in the previous section, TC values can be clearly understood by examining by the characteristic values of the logistic function. Therefore, $A_2$ was defined as the bifurcation point to obtain the triangular patterns, and the maximum value of $A_2$, 2.5 defined as the threshold for the explanatory variable A. In this way, considering the range of $\delta$ that can be obtained by triangular patterns according to Eq. (13), we predicted the bifurcation point as a function of $\delta$

$$D_u/\sqrt{1-\delta }-D_u^c\le A_2\left(D_u^c-D_u/\sqrt{1+\delta }\right).$$E16

Transforming this into the $\delta$ condition, we obtained

$$\delta \le 1-1/B^2$$E17

where $B=\left(A_2+1\right)D_u^c/D_u-A_2/\sqrt{2}$. Once $D_u$ and $D_v$ are determined, $D_u^c$ can be calculated using Eq. (6). Subsequently, the bifurcation point of $\delta$ was determined.

In the case of $\alpha =0.50,\beta =0.01,\gamma =26.0,D_v=5.00\times {10}^{-2},$ and $D_u=5.00\times {10}^{-4}$ which are the same parameters as those in Figure 2, the bifurcation point of $\delta$ was 0.856.

5. Discussion

In this section, we explain how the triangular patterns were obtained by numerical simulation of Eqs. (7)-(9) as a two-dimensional Turing pattern. As shown in Figure 3, a triangular pattern was obtained using a complex combination of Fourier modes, including high wavelengths.

In this section, we discuss when the reaction term in (2). We also examined whether the same conclusion holds for other choice of reaction terms. We examined another type of reaction term introduced by Schnackenberg [13]. The model described a simple chemical reaction for glycolysis given by the following reaction terms

$$f\left(u,v\right)=\alpha -u+u^{2}v,\mathit{\text{and}}g\left(u,v\right)=\beta -u^{2}v,$$E18

where $\alpha$ and $\beta$ are positive constants. The anisotropic diffusion coefficient is of the same type as in Eqs. (8) and (9). We obtained results identical to those of the FitzHugh-Nagumo model for the reaction terms in Eq. (2) [10]. Thus, we conclude that the choice of the functional form of the reaction does not strongly affect the determination of periodic triangular patterns.

Moreover, we examined the substance independence of the diffusion anisotropy. We also examined not only the case in which diffusion anisotropy is introduced only into the inhibitor v, a substance with the faster diffusion, but also the case in which diffusion anisotropies are introduced into both the activator and the inhibitor. In both cases, we obtained very similar results as in the case with the result explained above (results not shown). From these results, we conclude that these properties are not related to the substance independence of the diffusion anisotropy.

Thus far, our understanding has not progressed to the point where we can fully describe the mechanism of triangular pattern formation for Turing patterns. However, the conditions under which the triangular pattern is formed and maintained as a motionless fixed pattern, as represented in this chapter, were derived mainly through numerical calculations. We hope that mathematical progress will be achieved in this area.

In contrast, this model was also created to automatically extract triangle images. Therefore, we also hope that this model will be relevant and developed for applications in image analysis application.

Acknowledgments

We thank Prof. M. Mimura and Prof. K. Osaki for helpful comments, and K. Yamada for technical helps.