Travelling wave and Turing patterns in subdiffusive autocatalytic systems

In this work, we analyse autocatalytic reactions in complex and disordered media which are governed by subdiffusion. The mean square displacement of molecules here scale as 𝑡 𝛾 where 0 < 𝛾 < 1 . These systems are governed by fractional partial differential equations. Two systems are analysed i) in the first a logistic growth expression is used to represent the growth kinetics of bacteria. Here the system dynamics is governed by a single variable. ii) the second system is a two variable cubic autocatalytic system in a porous media. Here each reactant is involved in the autocatalytic generation of the other. These systems have multiple steady states. They exhibit traveling waves moving from an unstable steady state to a stable steady state. The minimum wave velocity has been obtained from phase plane analysis analytically for the first system. In addition, the two variable system also shows Turing patterns in selected regions of parameter space. The stability boundary for Turing patterns for subdiffusive system is found to be the same as that for regular diffusive systems obtained by Seshai et al. [1]. System behaviour as predicted by the stability analysis is verified using a robust implicit numerical method based on L1 scheme. trajectories in the phase space on a two-dimensional plane. There is a trajectory connecting 𝑆𝑆3 to 𝑆𝑆2 and 𝑆𝑆1 to 𝑆𝑆2 . These represent travelling wave solutions.


Introduction
Two commonly exhibited spatio-temporal patterns observed in natural systems are: travelling wave solutions and Turing patterns [2]. These patterns arise from the interaction of nonlinear reaction and diffusion and have been intensively studied in the literature [3]. The reaction diffusion model has been extensively used to describe several phenomena observed in nature such as patterns on animal skin [4], pattern formation in ecological systems [5], epidemiological patterns [6], instabilities in predator-prey systems [7]. The reactions in these systems are governed by a positive feedback effect and are typically autocatalytic in nature. These systems are characterised by a multiplicity of steady states. Autocatalysis has been used to describe self-regulated pattern formation in diverse fields such as RNA replication [8], reaction engineering [9], etc. Travelling waves originate from an unstable steady state and terminate at a stable steady state. These waves play a very important role in transmitting biological information in many natural processes such as Ca 2+ waves during fertilization [10] and embryogenesis [11] etc. While travelling waves can be exhibited by a single variable system. Turing patterns require at least a two-variable system. Travelling waves are characterized by a constant speed and shape. The Fisher-Kolmogorov equation has been analysed extensively for travelling wave solutions [3]. Here the reaction is governed by the logistic reaction kinetics. This model has been used to represent several biological and chemical systems [12].
Travelling wave solutions describing bacterial transport and chemotaxis have been studied in the literature [13,14]. The first attempt to understand these migrating bands mathematically was made by Keller and Segal [15]. They established a travelling wave solution using a pair of reaction diffusion equations to describe the bacterial population and the concentration of the attractant. The Keller-Segel (KS) Model ignored the growth kinetics, an essential factor in the bacterial expansion process. As a result, the travelling waves obtained lose stability. Later, researchers modeled bacterial expansion by including growth rate and demonstrated the existence of stable travelling waves. However, these results did not match with experimental observations [16,17]. Recently, Cremer et al. [14] have experimentally shown the travelling bands of the bacterial population in the presence of a nutrient and chemoattractant. They have also numerically confirmed that growth expression is independent of the attractant even when the model includes the chemotactic effect.
Many microbial processes such as bio fertilization and remediation occur in porous media in nature [18,19]. These processes are affected by the migration of microbes through porous media. Despite the importance of bacterial movement in disordered and complex porous media, the process is not well understood. Bacterial movement in porous media gets restricted by cell wall interactions and viscous drag forces caused by confinements [20]. Bhattacharjee et al. experimentally observed that when bacteria are trapped in pores, they reorient their bodies until they escape from the traps [21]. The trap durations are random and widely distributed from a few seconds to minutes. They found that the jump lengths in pore space while being random are narrowly distributed. It has been proposed that bacterial movement governed by entrapment and collision with the solid wall can be viewed as anomalous subdiffusion [21][22][23][24]. The collective movement of bacteria generates travelling waves. These waves are responsible for the formation of biofilm and colonies of different shapes and morphology.
In the context of chemical species, pattern formation also arises from the interaction between diffusion and reaction. Chemical species such as macromolecules or polymers often encounter molecular crowding due to their large size. In porous media, these molecules get trapped in pores and their transport is characterised by anomalous subdiffusion. The transport of molecules and the reaction kinetics are affected by the subdiffusive behaviour caused by confinements [25][26][27].
Most of the earlier work in literature have discussed autocatalytic reactions governed by a single reactant [3,9,14]. However, reactions involving isomers can be autocatalytic in either species. Each species is engaged in the generation of other species. The system of regular diffusion and cubic autocatalysis has been studied by Seshasai et al. [1]. However, when these reactions occur in disordered and complex media where molecular crowding is dominant, their behaviour is governed by anomalous diffusion. Researchers have also found that confinement and collision with the wall can affect reaction and movement of species [22,28,29].
Researchers have studied the interaction of regular diffusion with reactions extensively in the literature. The role of subdiffusive behaviour in the formation of Turing patterns and travelling waves has not been studied in detail, primarily because the modeling of subdiffusive systems involving nonlinear reaction kinetics is still a challenge. Over time researchers have adopted many approaches to model subdiffusive processes [30,31]. The theoretical description of the fractional reaction subdiffusive process is based on continuous time random walk (CTRW) and fractional dynamic approach [25,30]. Henry et al. have provided a detailed description to model the interaction of linear reaction terms and subdiffusion [32]. These models were derived by the CTRW approach for recombination kinetics and instantaneous creation and annihilation processes in subdiffusive media [25,28]. Here anomalous diffusion is modeled by a fractional differential equation [33].
Many biological experiments demonstrate the prevalence of anomalous diffusion. This has led to a rigorous effort to study the dynamics of fractional order systems. Huang and Wang have investigated the general stability conditions of a class of fractional-order chaotic nonlinear systems [34]. Cermak and Nechvatal have discussed Routh-Hurwitz conditions of fractional type and investigated the fractional Lorenz system [35].
Several numerical methods to solve fractional differential equations have been proposed in the literature, such as variational iteration method, Adomian decomposition method, finite difference methods, spectral methods etc. [36]. Yuste and Acedo et al. [37] have developed an explicit finite difference method for the fractional diffusion equation. The method is based on Grunwald-Letnikov definition for fractional derivatives. Later Chen et al. [38] have developed an implicit numerical method based on Grunwald-Letnikov definition to solve fractional reaction subdiffusion equations. Zhuang et al. have presented an implicit method based on integration methods [39]. The Riemann-Liouville definition of fractional derivative has difficulties at = 0 where some terms become unbounded. To avoid this issue Langlands and Henry [40] have presented an implicit L1 scheme for fractional subdiffusion equation. Nikan et al. [41] have introduced a radial basis function-generated finite difference method for solving the fractional differential equation. They have used Grunwald-Letnikov definition with first-order accuracy to discretize the equation. However, Fractional differential equations with nonlinear reaction terms have not been analysed in the literature. We can apply an implicit numerical scheme to discretise these equations. This results in a system of nonlinear algebraic equations which are solved iteratively using standard methods. The solution to the fractional differential equations is very sensitive to the chosen time step. Most of the numerical schemes for discretisation of fractional derivatives involve summation over all previous values. This involves a higher computational effort, particularly for nonlinear systems. These problems become even more severe in higher dimensions and are accompanied by an increase in simulation time.
In this work, our focus is to study the interaction of nonlinear kinetics and subdiffusion affect. The objective is to determine how subdiffusion affects the origin of various spatiotemporal patterns such as travelling waves and Turing patterns. This is of relevance in porous media where transport maybe governed by anomalous diffusion. We analyse two such systems. First, we consider a single variable system of bacterial growth expansion, which exhibits chemotaxis and subdiffusion. The single variable system can exhibit only travelling patterns and our aim to study effect of subdiffusion on traveling wave solutions and obtain insights on wave propagation features. In the second system, we consider a cubic autocatalytic reaction along with subdiffusion in two species in which each species autocatalytically generate the other species. In this system, we investigate the effect of subdiffusion on travelling wave patterns and Turing patterns. We also seek to understand how the Turing space (the parameter space where Turning patterns are observed) gets modified because of subdiffusive affects. Another important contribution of this work is the development of an implicit L1 scheme for non-linear subdiffusive systems which is numerically efficient. The simulations are used to verify the prediction of the linear stability analysis.

System description and governing equations 2.1 Fractional reaction subdiffusion model for bacterial growth
The bacterial growth model of a subdiffusive system exhibiting chemotaxis and logistic kinetics is described by a fractional partial differential equation. The detailed derivation is given in appendix A and B. A continuum approach based on asymptotic analysis leads to the following evolution equation describing bacterial growth and transport in porous media.
Where, , and represents bacterial diffusivity, chemotactic coefficient and growth rate constant, respectively. ( , ) represents the bacterial concentration as a function of space and time. A linear concentration profile of chemoattractant = + 0 is imposed. This transforms the above equation to In the above model, the first term on the right-hand side represents subdiffusive transport of bacteria, the second term represents chemotactic drift and the last term represents bacterial growth.

Fractional reaction subdiffusion model for cubic autocatalytic system
In this section, we present an autocatalytic cubic reaction system in porous media. The system is described by a fractional reaction subdiffusion process.
We assume species and are generated by precursors 1 and 2 with rate constants and respectively. These reactions are represented by 2 → and undergo an autocatalytic reaction governed by Here, and are rate constant of reaction generating and autocatalytically, respectively. Species and also follow decay reactions are decay rate constant for reactants and , respectively.
The fractional-order model equation describing spatio-temporal evolution in porous media sustaining this network of reaction is given by where and are diffusivity of and .
The above system given by (6), can be written in the following form where,

Numerical methods
We have used a robust numerical method based on the L1 scheme [36,40] to numerically solve the fractional differential equation. The scheme is implicit, and at each time step, the discretized equation is considered as a system of nonlinear algebraic equations. We use the 'fsolve' function in MATLAB at every time step to solve this system of equations.
The numerical method based on the L1 scheme is very sensitive to the time step chosen and, it is advisable to use a small-time step for accuracy. The detailed numerical method to analyse fractional differential equation is described in appendix C.
We use a finite difference scheme for discretising the differential equations. The Euler forward difference is used for the first-order time derivative and the second-order spatial derivative is approximated using a second-order central difference scheme. The firstorder spatial derivative is approximated using the central difference scheme or upwind depending on value. We use the implicit L1 scheme to discretize the fractional derivative that is valid for 0 < < 1.
The numerical scheme developed is used to simulate both travelling waves and Turing patterns.

Travelling wave solutions
We seek travelling wave solutions of the fractional differential equation (2) with boundary conditions (−∞) = 1, and (∞) = 0. The boundary conditions are chosen so that each end of the boundary is at a particular steady state.
Operating by −1 −1 on both sides of (2) we obtain, Here, , and are constant.
We used the following fractional complex transform [42,43] to transform the above equation to an ordinary differential equation.
where, is a constant. This reduces (8) into the following ODE Defining = we obtain the following system of ordinary differential equations We analyse the system in the half-plane as cannot be less than zero {( , ): > 0, −∞ < < ∞}. The system has two steady states 1: (0,0) and 2: (1,0). The travelling wave solutions of (8) are phase plane trajectories of the system (11), connecting these steady states.

Phase plane analysis
Phase plane analysis enables us to analyse the qualitative behaviour of dynamical systems. The phase plane portraits are plotted in MATLAB using the function 'pplane' [46]. A linear analysis around 1: (0,0) and 2: (1,0) shows that 1 is a stable node provided is sufficiently high while 2 is a saddle point. The nature of these steady states is tabulated in Table 1.

Saddle
Phase plane trajectories for various are illustrated in Fig. 1. If 0 < < − 2√ , there is a trajectory going from 1 to 2 as shown in Fig. 1(a). These trajectories are physically unrealistic. As is increased further and lies in the range − 2√ ≤ < + 2√ There are trajectories from 1 to 2. These again are physically unrealistic since < 0 as here spirals around 1 as shown in Fig. 1(b) and (c). For all ≥ + 2√ there is a heteroclinic orbit connecting 2 to 1 lying entirely in half plane ≥ 0 as shown in Fig. 1(d). Hence, the minimum value of of the system is = + 2√ . We have numerically solved (2) for different values to obtain the travelling wave solutions. The numerical simulations are subject to the following initial and boundary conditions.
The initial condition Boundary conditions For this combination of boundary conditions and initial conditions described by (14) and (15) From the figures shown in Fig. 2, it can be seen that the travelling wave moves from left to right with a constant speed and a fixed shape. The travelling wave has deviated from the starting steady state = 1 (there is a sudden dip). This deviation or dip decreases as time increases. This arises as the asymptotic behaviour of fractional-order systems follows power law [44]. It is also observed that the dip or deviation reduces as we increase the fractional power as shown in Fig. 2(b).

Travelling wave solutions
The well-mixed system with fractional order corresponding to (6) has three steady states. The phase plane trajectories of this system obtained using L1 scheme are plotted in Fig.  3(a). The phase plane analysis shows that there is one trajectory connecting the unstable steady state SS2 to the stable steady states 1 and 3 as shown in Fig. 3(a). trajectories in the phase space on a two-dimensional plane. There is a trajectory connecting 3 to 2 and 1 to 2. These represent travelling wave solutions.
We are interested in travelling wave solutions of fractional reaction system described by the subdiffusion equation (6). We use the fractional complex transform given by (9) to transform the one-dimensional form of (7). We transform the independent variables ( , ) to as given by (9) The phase plane trajectories are plotted by integrating the above system of equation (18). We see that there are two trajectories connecting 1 to 2 and 3 to 2 as shown in Fig. 3(b). The travelling waves observed correspond to these trajectories. We also observe that the direction of these trajectories is opposite to the trajectories shown in the phase plane of the well-mixed system Fig. 3(a).
We have numerically solved the system of equation given by (6) for different values and obtained the travelling wave solutions shown in Fig. 4. The numerical simulations are subject to the following initial and boundary conditions. The boundary conditions are chosen in such way that waves start from unstable steady states and moves towards a stable steady state.
We use the implicit numerical scheme described in Appendix C with = 50, = 0.1, = 10, = 0.01, to solve (6). The wave starts from the unstable steady states 1 and 3 and moves towards the stable point 2 with a constant wave speed for the chosen parameters. The waves travel from left to right. The waves from 3 to 2 are shown in Fig. 4(a) and (b). The travelling waves from 1 to 2 are shown in Fig. 4(c) and (d).

Turing patterns
It is well known that steady spatial Turing patterns are induced by diffusion under some conditions in nonlinear systems. We now analyse the formation of these patterns in subdiffusive systems. Here the steady states of the well-mixed system are stable and introduction of diffusive transport renders an instability under some conditions.
The stability of the well-mixed system is governed by the eigen values of the following Jacobian matrix evaluated at steady state where, and are the derivatives of and with respect to species . and represent species and , respectively.
The characteristic equation whose roots govern the stability is The spatially homogeneous steady state of fractional order system is stable if and only if [45] If we assume that eigen values are real then the is always zero and stability is decided by the signs of ( ℎ ) and ( ℎ ) [45].
We assume that the trace of Jacobian ℎ is negative and determinant is positive to ensure the stability of the steady state and also assume eigen values are real. This implies ( ℎ ) = + < 0 and − > 0. Hence, we expect the fractional reaction system exhibits Turing pattern for ( ℎ ) < 0. The critical surface across which Turning patterns appear is given by ( ℎ ) = 0.
To study the linear stability of the steady state of (7), we investigate the evolution of disturbances which are periodic in space. For this we substitute = 0 +̃ , = 0 +̃ in (7), where ( 0 , 0 ) is a steady state. Retaining first order terms in the disturbance represented by the ′~′ variable yields The stability of the steady state of fractional reaction subdiffusion system (24) is determined by eigen values of the following Jacobian matrix (25).
Where 1 and 2 are eigen values of the system.
The trace of Jacobian The determinant of Jacobian The ( ) always remains negative as we assume + < 0 for the well-mixed system. This implies that there is no Hopf bifurcation [45]. Hence, the necessary and sufficient condition for stability of fractional reaction subdiffusion system i.e., 0 < ≤ 1 is ( ) > 0. This condition is identical to that used in modelling regular diffusion and reaction systems [1].
From (27) a necessary condition for instability is + > 0. Considering this Along with the Trace condition we conclude that and should be of opposite signs. Without loss of generality, we consider > 0 and < 0. So, can be considered as an activator and as an inhibitor, which implies > − > 1. From this condition, we can conclude that for Turing pattern formation, diffusivity of inhibitor must be greater than diffusivity of activator.
The critical surface across which we have Turing patterns is ( ) = 0. We determine this critical surface, in the space of two parameters and 2 . The critical surface along which ( ) = 0 divides the space into two regions as shown in Fig. 5. In region I, ( ) < 0, the system exhibits Turing patters, when ( ) > 0 as in region II, we get spatially homogeneous stable steady. For a system varying in only one spatial direction, numerical simulations were carried out with no flux boundary conditions. Periodic concentrations of species and were obtained for the parameters in the region I and spatially homogeneous profile for the parameters in region II as shown in Fig. 6. As we increase the diffusivity ratio of to , we observe that there is a reduction in the height of the peak and an increase in its width as shown in Fig. 7(a) and (b). We have performed numerical simulation for different 2 − combination and for different fractional powers, the observed results are similar to those reported in Fig. 6 and Fig. 7. This verifies the stability boundary. Simulations were also carried out in a 2D domain for the parameters in the region I. The obtained profiles are plotted in Fig. 8.

Conclusions
In this work, we have used a numerical method based on the L1 scheme to analyse the fractional reaction subdiffusive system. Our scheme is implicit and robust and has been used to predict travelling wave solutions and Turing patterns.
The single variable system governing the bacterial dynamics was investigated systematically. For this, the stability of steady states was examined by linear stability analysis and the trajectories in the phase plane. We have obtained the expression of minimum wave speed of bacterial expansion analytically. We have numerically simulated the travelling wave pattern and observed stable travelling waves connecting two steady states. The wave starts from the saddle point and terminates at the stable node as predicted by phase plane analysis. The travelling wave solutions provide insight of spatiotemporal patterns in bacterial systems.
The fractional reaction subdiffusion system of two variables for cubic autocatalytic kinetics has also been analysed. The interaction of nonlinear kinetics and subdiffusion results in the travelling wave patterns. These waves are responsible for the rapid transport of species. The fractional cubic autocatalytic system also shows Turing patterns when the ratio of diffusivity of activator to inhibitor is lower than a critical value. The critical surface boundary has also been obtained by linear stability analysis. The boundary of this system is the same as the regular cubic autocatalytic system studied by Seshai et al. [1]. The only difference is that the fractional system takes a much longer time to reach equilibrium. We also plotted the Turing pattern for a different fractional power of the derivative and observed that the steady state behaviour does not depend on the fractional derivative.
The fractional reaction subdiffusion model can be advantageous in studying biological systems where species transport is hindered by traps, macromolecular crowding, or other obstacles. The fractional power of derivative can be considered as a parameter that depends on pore size distribution for a particular system.

Data Availability Statement
trapped for a random time in pores. The detailed derivation is given in reference [30]. The behaviour is captured by the master equation.
( ) = (0)∅( ) where ( ) is a waiting time density, In equation (A1), the number of bacteria at grid point at time is comprised of those bacteria that were at point at time = 0 , which have not yet jumped (the first term on the right-hand side) together with those bacteria that were at an adjacent grid point ± 1 at an earlier time ' but then jumped to point at time t after waiting a time − '.
( , ) and ( , ) are the probabilities of bacteria jumping from to the adjacent grid point to the right and left directions, respectively.
Where, ( , ) = ( , ) is a sensitivity function that depends on the chemoattractant concentration. We consider a heavy-tailed waiting-time density to represent the trap durations in which there is a higher probability of getting a large trap duration.
Where Ґ(. )is the Gamma function.
In (A8), the first term on the right hand of fractional differential equation represents subdiffusive transport and the second term represents chemotactic drift. The system behaviour is governed by two parameters: diffusivity ( ) that depends on the pore size distribution and chemotactic sensitivity ( ) that signifies how strong bacteria are responding to a particular chemical.

Appendix B Model for fractional reactional subdiffusion
If and are equal, (A7) results in subdiffusion. This is governed by Using scaling arguments, Yuste et al. [25] found that the dependency of reaction term and transport of species is affected by a subdiffusive process. Monte Carlo simulations have also shown that the reaction term can be represented as Where, is the rate of reaction.
Consequently, the updated fractional reaction and subdiffusion equation at long times has the following form, Here, the reaction rate has a memory effect that depends on waiting time distribution of reactants [27]. This is valid when the expansion of reaction zone is very slow and reactant concentrations vary over a typical length governed by diffusion. Here there is a flux of species towards the reaction region.

Appendix C
A numerical method to solve the fractional differential equation The Euler forward difference is used for the first-order time derivative as shown below The second-order spatial derivative is approximated using the second-order centred difference scheme evaluated at the next time step 2 ( , ) The first-order spatial derivative is approximated using the centred difference scheme or upwind depending on value, evaluated at the next time step The central difference scheme We use the implicit L1 scheme to discretize the fractional derivative that is valid for 0 ≤ ≤ 1 [36].
L1 scheme for solving differential equation Where, ∆ is the time step and Ґ(. ) is the Gamma function.
We illustrate this with the fractional differential equation given below