Abstract

Convolution representation manifests itself as an important tool in the reduction of partial differential equations. In this study, we consider the convolution representation of traveling pulses in reaction-diffusion systems. Under the adiabatic approximation of inhibitor, a two-component reaction-diffusion system is reduced to a one-component reaction-diffusion equation with a convolution term. To find the traveling speed in a reaction-diffusion system with a global coupling term, the stability of the standing pulse and the relation between traveling speed and bifurcation parameter are examined. Additionally, we consider the traveling pulses in the kernel-based Turing model. The stability of the spatially homogeneous state and most unstable wave number are examined. The practical utilities of the convolution representation of reaction-diffusion systems are discussed.

1. Introduction

Many phenomena are described by mathematical and physical models, including traffic flow, stripe of dune, bubble in the granular medium, and plasma flow [15]. Among these phenomena, spatio-temporal patterning is one of the attractive topics in the field of nonlinear physics. To describe these patterns by a simple mathematical model, many reaction diffusion (RD) systems have been proposed. In pattern formation in RD systems, the Turing instability is a key concept; the difference in diffusion rates of the activator and inhibitor leads to lateral inhibition. The Turing pattern is caused by lateral inhibition around the excited region: short-range activation and long-range inhibition [6]. Recently, the reduction of a two-component RD system to a one-component RD equation with a convolution term has been considered. Considering that the system converges to a stationary state, the inhibitor is expressed in the form of convolution using the activator. The convolution representation using a weight function, called a kernel, has a nonlocal effect on the activator. A Mexican-hat-type function is often used as a kernel [79]. The Mexican-hat-type kernel creates a suppressive region far from the central enhanced region, which imitates the lateral inhibition around the excited region, resulting in pattern formation [10]. In the 1970s, the convolution term was introduced in the time evolution equation of the membrane potential in a neural model [11]. This term caused a nonlocal effect on the membrane potential via the neighbour cells in the network. To choose adequate parameters in a model, the traveling wave was calculated self-consistently. Inspired by this study, traveling waves in a neural network have been considered [12, 13].

Among many Turing patterns in nature, the pigmentation on the skin of zebrafish has been well-reported [1417]. Fish pigmentation shows many patterns such as spot, stripe, and labyrinth. Critical pigment cells are yellow and black, which form interaction networks, resulting in Turing patterns. This pattern formation has been modeled using a multicomponent RD system. Given the variety of patterns found on the skin of zebrafish, the dependence of pattern formation on the kernel function was studied recently [18]. The multicomponent RD system was reduced to a single-component activator equation with a convolution term, called as kernel-based Turing model (KT model). The Mexican-hat-type kernel composed of two Gaussian functions with reverse sign was employed. Generally, patterns are sensitive to the shape of a kernel; amplitude, width, distance between centers of a source, and the peak position of the distribution. Therefore, variations in observed patterns in nature are due to variations in the kernel. Additionally, we can predict spatial patterns in two dimensions by using the KT model. That is, we first obtain the most unstable mode (wave number) and the corresponding eigenvalues through linear stability analysis. Considering this information in the kernel, we can predict what pattern will occur when the spatially homogeneous state becomes unstable [19]. From these perspectives, the KT model is expected to be a promising equation for pattern formation. Generally, the KT model is valid when the system converges to a stationary state; transient patterns cannot be described by the KT model.

Aside from the static stationary patterns, many temporally changing patterns exist: stripe patterns in nematic liquid crystal, traveling band in colloidal suspension, traveling waves in shallow-water, and fluidized granular media [2024]. Although these dynamical patterns were modeled and analyzed by many mathematical methods, not many studies focused on the traveling waves in the RD equation using the kernel function [2528]. The interaction of traveling pulses and fronts in one dimension described by convolution representation was studied [29]. The interaction between two adjacent waves decayed exponentially with an increase in the distance between them. Although the asymptotic behavior of interacting traveling waves was considered, an explicit traveling wave solution using convolution representation was not given. Stationary traveling waves in one and two dimensions have been observed in experiments [3032], traveling waves in the RD equation with the convolution representation should be investigated.

In this study, we consider the procedure to transform traveling pulses in RD systems into convolution representation, employing a mathematically tractable RD system. Accordingly, we consider two types of RD systems. The first one is a one-dimensional RD system with a global coupling term. An experimental setup corresponding to the global coupling term has been considered: the current filaments in gas-discharge and --- diode systems and the chemical reaction on catalytic ribbon surface. In electric circuit systems, although the current filament distributes on the electrode, the total current is feedback-controlled [33, 34]. Similarly, in the chemical reaction on the catalytic ribbon, the spatial average temperature of the surface is kept constant [35]. These feedback controls are modeled by space-averaged terms called global couplings. Although the stability of standing pulse in these systems has been studied, the properties of traveling pulse solutions have not been clarified yet. For the second RD system, we consider a usual two-component RD system with cubic nonlinearity. We derive the KT model from this system and consider localized pulses from the perspective of Turing instability. The remainder of this paper is organized as follows. In Section 2, we introduce an RD system with a global coupling term. We consider the stability of the standing pulse solution and show the relation between the traveling speed and bifurcation parameter. Using these results, we demonstrate the choice of an adequate speed in the kernel representation of traveling pulse. In Section 5, we consider pulse solutions in the KT model. To confirm the validity of numerical simulations, the stability of the spatially homogeneous solution is examined. Finally, the application of our study to experiments and further extensions is discussed in Section 8.

2. RD System with Global Coupling

We first introduce a two-component RD system with a global coupling term in one dimension [36]. Activator and inhibitor satisfy the following time evolution equations: where where is the step function satisfying for and for We choose and such that and , respectively. , , and are positive constants. When the magnitude of global coupling , is a constant, , and the RD system (1) is a local coupling system, called the McKean model. By contrast, when , the RD system (1) is a global coupling system; when is large, the second term of controls such that the width of excited region is in the limit

Let us briefly summarize the solution of the RD system (1) [37]. First, we consider the case . When is large, there is a stationary pulse solution, the standing pulse. Although the pulse has left and right interfaces in the limit , these interfaces become a transition layer (inner layer) with a width for finite . On decreasing , the standing pulse is destabilized through the oscillatory bifurcation; the interfaces oscillate symmetrically with respect to the center of the pulse (breathing motion). On decreasing further, the pulse is secondarily destabilized through the translational bifurcation, resulting in the traveling pulse. The velocity of the traveling pulse increases as decreasing . By contrast, when , on decreasing , the standing pulse is primarily destabilized through the translational bifurcation; the oscillatory bifurcation is suppressed strongly for large . The global coupling term in the definition of suggests the strict restriction of for large . This condition prohibits the breathing motion of the interface, and the translational bifurcation occurs primarily on decreasing .

3. Convolution Representation of the RD System

For smaller , the RD system (1) has traveling pulse solutions. We consider the convolution representation of the stationary traveling pulse solution of the RD system (1). In the moving coordinate , where is a traveling speed, the time evolution equations of and are

In the interval , the stationary solution of in the system (3), subject to periodic boundary conditions, is given in the form of cyclic convolution. The detailed derivation is given in Appendix A, and the final expression is where . Using Equation (4), the time evolution of is where corresponds to a kernel and the cyclic convolution is defined as with where is a sign function, satisfying for and for . For sufficiently large with , kernel Equation (7) is approximately , the kernel for diffusion. When with , we note that Equation (7) recovers the results in ref. [38]. Equation (5) yields a time evolution of for a given initial condition until converges to a stationary solution; it is valid for given , , and , for which the RD system (3) has a stationary traveling pulse. Unlike the Korteweg-de Vries (KdV) equation [39], traveling speed is not arbitrary. To find the adequate value of for given and , we perform stability analysis of the standing pulse in the RD system (1).

4. Stability Analysis Using a Singular Perturbation Method

For the RD system (1), the stability of the standing pulse is examined using the singular perturbation method [37, 40]. The exponential growth of perturbation to the standing pulse solution is determined by the roots of . Here, is the standing pulse width, and is a stationary solution, symmetric with respect to , of the RD system (1) in the limit The brief derivation of the stability formula is given in Appendix B, and the final expression is

In the derivation of Equation (8), we assumed that perturbations were supplied to the stationary position, and the left and right interface positions were shifted as and , respectively. We consider two modes of perturbation: the symmetric perturbation with respect to the center of the pulse and antisymmetric one . and correspond to the symmetric and antisymmetric cases, respectively. We observe that , which corresponds to the translational invariance. When is degenerated at , that is, the standing pulse is destabilized through the translational bifurcation; condition (9) gives as

Meanwhile, when has a pair of the imaginary solutions with real number , the standing pulse is destabilized through the oscillatory bifurcation; interfaces oscillate symmetrically with respect to the center of the pulse. This condition gives , which is given by the solution

Using these formula, we consider the dependence of bifurcation of standing pulse on and in Section 7.

5. Traveling Pulse Solution of the KT Model

In the above discussion, we considered the convolution representation of traveling pulse in an RD system where the inhibitor was expressed by using kernel integral in the stationary state. By contrast, there is another type of RD equation with convolution representation: the KT model. The Mexican-hat-type kernel is employed in this model.

Let us consider the following RD system with a cubic nonlinear function in the activator: where where is a constant in the range . We consider the stationary pulse solutions in the interval under periodic boundary conditions. We use the Mexican-hat-type kernel, which causes the activation in the short range while suppressing it in the long range.

To consider the traveling pulse, we introduce the moving coordinate , where is a traveling speed. Following a similar procedure as in Section 3, the RD equation with a convolution term, the KT model, is proposed where with kernel where and . When , Equation (16) recovers the kernel of a standing pulse in ref. [38]. To make the Mexican-hat-type kernel, we assume that and is under condition . The KT model (14) is the time evolution equation of with an initial condition , which is valid until converges to a stationary solution. In the KT model, parameters , , and are not restricted, except for condition .

In Section 3, Equation (6) represents stationary solution of , and Equation (7) corresponds to its kernel, which is composed of a single term. By contrast, in the KT model, the kernel Equation (16) is composed of two similar terms with the reverse sign. Because , the second term of Equation (16) operates suppressively on far from the center of the pulse. For a larger , a suppressive region of the kernel is wider.

6. Linear Stability of the Homogeneous Solution of the KT Model

In this section, we consider the Turing instability in a two-component RD system (12). There is a spatially homogeneous solution of the RD system (12); . The KT model corresponding to the RD system (12) is proposed as where and is given by Equation (16) with .

We consider the pulse solution of the KT model (17) in the interval under periodic boundary conditions. To clarify the mechanism of the pulse solution occurrence in the KT model (17), we consider the linear stability of the spatially homogeneous solution. The growth rate of instability is obtained through linear stability analysis. A brief derivation is given in Appendix D, and the eigenvalue is where is and represents the wave number that takes the values of . and are given in Section 5. Although the above eigenvalue is derived for the stationary state with , the formulation is valid in the moving coordinate system with finite . The linear stability analysis of the spatially homogeneous solution in the two-component RD system was considered in refs. [6, 7]. The eigenvalue was given by a quadratic function in terms of resulting from the diffusion terms of the activator and inhibitor in the system. By contrast, in the stability analysis of the KT model, the diffusion of the inhibitor is considered in the kernel function. This results in a different expression of an eigenvalue. Using this formula, we consider the validity of the pulse solutions in the KT model in Section 7.

7. Numerical Results

In this section, we consider the numerical results. In the calculation of the RD system (1), the time evolution of using kernel representation Equations (5) and (14), we fix , .275, .3, and and choose differences and . The calculation of spatial direction is performed using an implicit method, Crank-Nicolson method, under periodic boundary conditions. Here, is in the range of .

We first show the profiles of and of the stationary pulses in the RD system (1) in Figure 1. The standing pulse is symmetric with respect to and remains in the same position. Although the standing pulse is stable for large , it is destabilized through the translational bifurcation on decreasing . The traveling pulse moves to the right direction with a constant speed. The traveling speed and pulse width are determined self-consistently in the RD system (1).

To consider the stability of the stationary standing pulse, we show the dependence of the pulse width on [37]. We consider the symmetric standing pulse with respect to , the position of interface is , and the pulse width . There are two cases depending on , as shown in Figure 2. When , , on increasing , changes and converges to . The stability of the standing pulse is calculated using . The result is shown in Figure 3. When , the standing pulse is first destabilized through the oscillatory bifurcation on decreasing . By contrast, for large finite values of , the oscillatory bifurcation is suppressed, and the standing pulse is first destabilized through the translational bifurcation on decreasing . The suppression of the oscillatory bifurcation by global coupling results from the conservation of the pulse width. Under strong global coupling, although the deviation of the pulse width is strictly prohibited, the translational motion with the same width is allowed. We briefly summarize the stability of the standing pulse solution of the RD system (1). When , the RD system (1) is a local coupling system, and standing pulse solutions exist for large . On decreasing , the standing pulse solutions are destabilized through the oscillatory bifurcation, and the interfaces oscillate (breathing motion). On decreasing further, the standing pulses with breathing motion bifurcate secondarily to traveling pulse solutions. By contrast, when , the RD system (1) is a global coupling system. Although standing pulse solutions exist for large , the standing pulse is destabilized through a translational bifurcation on decreasing because the oscillatory bifurcation is suppressed. Thus, the traveling pulses primarily occur on decreasing .

Although the traveling speed of the soliton solution of the KdV equation is an arbitrary parameter, that of the RD system is determined self-consistently. To obtain a stationary traveling pulse solution by using Equation (5), it is necessary to choose an adequate value of for given and . Accordingly, we first consider the dependences of velocity on and , where is defined as the interval of two interfaces in the limit The relation between , , and is calculated self-consistently by using (given by Equation (B.1)), where represents evaluated at interface positions . The brief derivation of is given in Appendix C. The dependence of on is shown in Figure 4(a). For , the traveling pulse appears subcritically at ; the traveling pulse appears suddenly with a finite value of , which increases with decreasing [4144]. By contrast, for with .4, the traveling pulse appears supercritically at with . The dependence of on is shown in Figure 4(b). For , increases with increasing . By contrast, for , increases in small range of , which is almost constant. This property results from the strong global coupling in the RD system (1); the oscillatory instability is suppressed strongly [36, 37]. By the definition of , Equation (2), the pulse width is for large in the limit From these results, the value of for chosen and is restricted such that the standing pulse is destabilized through the translational bifuraction; here, is less than . Additionally, should be less than to find the value of on - relation, as shown in Figure 4(a). For , the value of is chosen on the upper branch with solid curve in Figure 4(a); the upper (lower) branch is stable (unstable) [41, 44]. In the case of , the lower branch disappears, and the value of is chosen on the dashed curve in Figure 4(a). The dependences of on and are confirmed by the simulation of the RD system (1). The data points are in good agreement with theoretical curves, except for a small value of . This discrepancy results from the fact that is chosen small but a finite value in the RD system (1).

Thus, we have considered the stability of the standing pulse and its dependence on , , and in the two-component RD system (1). Consequently, the standing pulse is destabilized by decreasing . Based on this knowledge, we consider the traveling pulse in a one-component RD equation with a convolution term (5). The kernel given by Equation (7) is shown in Figure 5. When , the kernel is symmetric with respect to and is positive for all . We note that the kernel is not a Mexican-hat type. Although there is no suppressive region in the kernel, the convolution term operates suppressively on in Equation (5) owing to its sign: minus sign as . On increasing , the kernel becomes asymmetric with respect to and remains a positive value (see dashed curve in Figure 5). The convolution term maintains the traveling pulse solution for . As examples of the pulse, choosing the adequate values of and , the profiles of and of the stationary pulse in Equation (5) are shown in Figure 6. Initial condition was given by a small positive excitation at . Comparing the profiles of and in Figure 6 with that in Figure 1, we note that these agree well with each other. Thus, both the standing and traveling pulses in the RD system (1) are recovered by using Equation (5) with an adequate value of determined by Figure 4(a). We note that the pulse solutions calculated by using Equation (5) are valid in the stationary state and the profiles in the transient state do not agree with those in the RD system (1).

Let us consider the pulse solutions of the KT model (14). The profiles of the kernel Equation (16) are shown in Figure 7(a). The kernel of the KT model is composed of two similar terms with reverse sign, and its profile is a Mexican-hat type. Although the kernel is positive around the center, it is negative far from the center. Using the KT model (14), the profiles of in the stationary state are calculated, where the initial condition was given by a random number with an amplitude of .1. The stationary profiles of are shown in Figure 7(b). When , the profile represents the standing pulse (solid curve). When is chosen positive, the profile represents the traveling pulse with the right direction in the moving coordinate (dashed curve). For sufficiently large , there remains only one traveling pulse in the system. The pulse number depends on and , as shown in Figure 8(a). As increasing and , asymmetry of the pulse is remarkable, and the interval of pulses is larger due to the strong lateral inhibition; the number of pulses in the system is reduced. To confirm the validity of the occurrence of the pulse solution in the KT model, the stability of a spatially homogeneous solution is considered by using Equation (19). The dependence of eigenvalue on wave number is shown in Figure 8(b). When is positive, the corresponding wave number destabilizes the homogeneous solution. The figure suggests that the homogeneous solution is unstable for . The most unstable mode, by which the maximum eigenvalue is attained, is for both and 0.5. In our simulation of the KT model (14), we obtain three (two) standing (traveling) pulses for , as shown in Figure 7(b). Thus, the number of pulses obtained by the simulation is reasonable, and the occurrence of pulses results from the Turing instability.

8. Conclusions

In this study, we considered the convolution representation of traveling pulses in a two-component RD system and the KT model. For analytical tractability, a McKean-type nonlinearity was adopted in a two-component RD system (1) with Equation (2). Under the assumption that the relaxation of an inhibitor was fast, a two-component RD system was reduced to a one-component RD equation with a convolution term given by Equation (5). To obtain the stationary standing pulse or traveling pulse by using Equation (5), it was necessary to proceed with numerical calculations until an activator converged to a stationary state. When the RD system (1) did not converge to a stationary state, our formulation was not valid. Additionally, the transient state calculated by using Equation (5) did not agree with that calculated by the RD system (1). In the one-component RD equation with a convolution term given by Equation (5), the traveling speed was unknown a priori. To find an adequate value of , the stability analysis of the standing pulse of the RD system (1) and - relation was necessary. The fluctuation of parameter around the bifurcation point in the case of the strong global coupling was much smaller than that in the case of local coupling. This is because that the oscillatory instability was strongly prohibited under a strong global coupling and the standing pulse was primarily destabilized through the translational bifurcation. Using an adequate value of for a chosen , we could obtain the stationary traveling pulse solution.

By contrast, in the previous studies of the KT model [18, 38], it was shown that a stationary standing pulse occurred using a symmetric Mexican-hat-type kernel. As the extension to the traveling pulse, we have shown that the traveling pulse was described by the asymmetric Mexican-hat-type kernel. The spatially homogeneous solution was destabilized through the Turing instability, and the stationary standing or traveling pulses occurred. In a multicomponent RD system, the Turing instability is caused by different magnitudes of the diffusion of activators and inhibitors. In the KT model, the kernel was composed of different magnitudes of a Gaussian-type function with a reverse sign, which formed the Mexican-hat-type kernel. The difference between and created a lateral inhibition, resulting in the localized standing or traveling pulses. The number of pulses in the system corresponded to the most unstable mode of the wave number. Because there is no restriction on parameters , , and in the KT model, it is necessary to choose them by considering other types of information, for example, matching patterns in the model with experimental observations. Another method to determine the parameters is to match the most unstable wave number obtained through the linear stability analysis with observation. The most unstable wave number is determined self-consistently by , , , and under the periodic boundary conditions.

We discuss the practical utilities of the convolution representation of RD systems. First, we discuss the global coupling system in the experiments. The current filament in the gas-discharge system [45, 46] corresponds directly to the theoretical model [36]. The total current in these gas-discharge systems is calculated by the integral of the current density over the electrode. The external resistor in the circuit corresponds to the magnitude of global coupling . On decreasing , the total current is increased, and the sequence of transition of current filament occurs; on decreasing , the spatially homogeneous filament is destabilized to a localized filament, an oscillatory filament (a rocking current filament; antisymmetric oscillation keeping a total current constant), and finally, a traveling filament. The difference between the experimental setup and theoretical model is that the integrand is the sum of the activator and inhibitor in the theoretical model [36]. Although this assumption is unpractical, the physics remains the same. When the traveling pulse is described by using a one-component RD equation with a convolution term, an adequate traveling speed is obtained following our procedures.

Second, we consider the -component RD system. There are many experimental setups that show traveling waves in one and two dimensions, such as nematic liquid crystal, camphor boat, and action potential in nerves [4749]. Reducing the multicomponent RD system to a single-component RD equation with a convolution term (the KT model), it is necessary to determine (), , and in the kernel. The stability condition of the spatially homogeneous solution is obtained by following similar procedures with this study. The most unstable wave number coincides with the number of waves in the system. Thus, , , and can be chosen such that the numerical results by the KT model match with the observed pattern. The multicomponent RD system or complicated cascade reactions can be reduced to an equivalent KT model, and the parameters in the kernel are determined by the above procedures.

Appendix

A. Derivation of the Kernel Function of a Traveling Pulse

In this appendix, we briefly derive the kernel of a traveling pulse in an RD system (1). The equations of an RD system are where and . In the moving coordinate , Equation (A.1) is

In the stationary state, the left-hand side of Equation (A.2) is zero. We derive the stationary solution in the interval under the periodic boundary conditions. The stationary solution is obtained by applying the variation of parameters as where and are coefficients. Imposing the periodic boundary conditions, and , on Equation (A.3), we can determine as

Substituting Equation (A.4) into Equation (A.3), we finally obtain as

We note that Equation (A.5) is in the form of cyclic convolution. Defining the convolution with kernel , is expressed as where is given by Equation (7) and is a sign function, satisfying for and for . In the derivation of Equation (A.3), we assumed that is in the stationary state. This assumption corresponds to the adiabatic approximation; is a fast variable compared to and is regarded as a stationary state.

B. Stability Formula by a Singular Perturbation Method

In this appendix, following the singular perturbation method [40], we derive briefly the stability formula given in Section 4. In a stationary state, we consider symmetric solution with respect to . Here, denotes the stationary standing pulse solution of the RD system (1) in the limit , the interfaces are located at , and the pulse width is . Perturbations are supplied to the standing pulse solution; we assume that the left and right interface positions are shifted as and , respectively.

The velocity of the pulse is given (for derivation, see ref. [42]) by where . Using Equation (B.1), the equations of motion on the left front at and the right front at are given by respectively, and corresponds to evaluated at the interface position . We expand around the interface position, as , where corresponds to at in the stationary state. The velocity is expanded around with as where and and represents evaluated using with . After several calculations, we obtain where is given by

To express Equation (B.2) in a linear form of and , we define the Laplace transformation of as

Using the facts that where sign in Equation (B.7) corresponds to the value at , respectively, we consider two modes of perturbations: the symmetric perturbation with respect to the center of the pulse and antisymmetric one . We apply the Laplace transformation to Equation (B.2) and finally obtain where and correspond to the symmetric and antisymmetric cases, respectively.

C. Traveling Pulse in the Limit of

In this appendix, we show the traveling pulse solution of the RD system (1). In the limit , the RD system (1) is greatly simplified. Putting , the traveling pulse satisfies a time evolution equation

The stationary traveling pulse solution is obtained under the boundary conditions in the limit and . The positions of left and right interfaces are denoted by and , respectively, and the final expression of is where and () are given by

The corresponding stationary solution is given by , and the pulse width is .

D. Stability Analysis of the Homogeneous Solution of the KT Model

In this appendix, following similar procedure used in ref. [38], we derive the linear growth rate of instability of the spatially homogeneous solution in the KT model (17). In our choice of parameters, there is a spatially homogeneous solution of the RD system (12); . We denote as and . We note that, in the moving coordinate , the system has a spatially homogeneous solution .

In the moving coordinate , the KT model (17) is where is given by Equation (18) and in the RD system (12) is expressed by using kernel in the KT model as where is given by Equation (16). We note that satisfies

Using this relation, , for an arbitrary spatially homogeneous solution . To consider the linear stability of Equation (D.1), we put with a small deviation . The corresponding deviation of is given by

We define the Fourier transformation of as where . Applying Fourier transformation to Equation (D.1), and putting , we obtain the eigenvalue as where is where and are given in Section 5. Because a pure imaginary number does not affect the stability of the solution, omitting the second term of Equation (D.6), we finally obtain the growth rate of instability of the spatially homogeneous solution as

Data Availability

All data that support the findings of this study are included in the article.

Conflicts of Interest

The author declares that there are no conflicts of interest.