[

# [

## Abstract

In this work an activator-depleted reaction-diffusion system is investigated on polar coordinates with the aim of exploring the relationship and the corresponding influence of domain size on the types of possible diffusion-driven instabilities. Quantitative relationships are found in the form of necessary conditions on the area of a disk-shape domain with respect to the diffusion and reaction rates for certain types of diffusion-driven instabilities to occur. Robust analytical methods are applied to find explicit expressions for the eigenvalues and eigenfunctions of the diffusion operator on a disk-shape domain with homogenous Neumann boundary conditions in polar coordinates. Spectral methods are applied using chebyshev non-periodic grid for the radial variable and Fourier periodic grid on the angular variable to verify the nodal lines and eigen-surfaces subject to the proposed analytical findings. The full classification of the parameter space in light of the bifurcation analysis is obtained and numerically verified by finding the solutions of the partitioning curves inducing such a classification. Furthermore, analytical results are found relating the area of (disk-shape) domain with reaction-diffusion rates in the form of necessary conditions for the different types of bifurcations. These results are on one hand presented in the form of mathematical theorems with rigorous proofs, and, on the other hand using finite element method, each claim of the corresponding theorems are verified by obtaining the theoretically predicted behaviour of the dynamics in the numerical simulations. Spatio-temporal periodic behaviour is demonstrated in the numerical solutions of the system for a proposed choice of parameters and a rigorous proof of the existence of infinitely many such points in the parameter plane is presented under a restriction on the area of the domain, with a lower bound in terms of reaction-diffusion rates.

R
\checkfont

eurm10 \checkfontmsam10 \newdefinitiondefinition[theorem]Definition European Journal of Applied Mathematics]Domain-dependent stability analysis and parameter classification of a reaction-diffusion model on spherical geometries W. Sarfaraz et al.]W.\nsS\lsA\lsR\lsF\lsA\lsR\lsA\lsZ,\ns and A.\nsM\lsA\lsD\lsZ\lsV\lsA\lsM\lsU\lsS\lsE 2000 \volume000 \pagerange[References

eaction-diffusion systems, Dynamical systems, Bifurcation analysis, Stability analysis, Turing diffusion-driven instability , Hopf Bifurcation, Transcritical bifurcation, Parameter spaces, Polar coordinates, Curved boundary

## 1 Introduction

Analysis of reaction-diffusion systems (RDSs) in the context of pattern formation is a widely studied topic [50, 56, 1, 2, 23, 30, 31, 32] in many branches of scientific research. Scholars of mathematical and computational biology [3, 4, 5, 6, 7, 11, 13, 17] devoted a great deal of attention to fully explore the theory of the dynamics governed by reaction-diffusion systems for a variety of reaction kinetics. Few of the routinely used and popular reaction kinetics in the theory of reaction-diffusion systems are activator-depleted [29, 28, 26, 27, 22, 20, 63], Gierer-Meinhardt [24, 23] and Thomas reaction kinetics [25, 37]. It is crucial to realise that the ultimate and complete knowledge encapsulating all aspects of reaction-diffusion systems for all types of possible reaction-kinetics is beyond the scope and feasibility of a single research paper or even a single book. Therefore, literature to date [37]-[50] on the analysis of reaction-diffusion systems falls under a natural classification of approaches based on the application of methods corresponding to the speciality of a particular scholar. It is in this spirit that the current work is a natural extension of [38], in which a full classification of the admissible parameter space associated to the dynamical behavior of activator-depleted reaction-diffusion system was explored on stationary rectangular domains. It is a common hypothesis in the theory of biological pattern formation [39, 13], that the emergence of spatially periodic pattern in the computational simulations of reaction-diffusion systems is connected to the properties of the associated eigenfunctions of the diffusion operator in the corresponding domain. The associated drawback with computational approaches [10, 12, 15, 16, 19, 21, 37, 53, 54, 55] for reaction-diffusion systems is that it lacks to provide rigorous insight of the role of the eigenfunctions in the emergence of spatial pattern, thus forming a natural platform to explore the evolution of spatial patterns from a perspective of dynamical systems. Numerous researchers have contributed to exploring the computational aspect of reaction-diffusion systems [4, 12, 14, 15, 26, 27], in which a restricted choice of parameter values are used with a partial reliance on trial and error to obtain the emergence of an evolving pattern either in space or in time. A large variety of scholars have also employed the analytical approach [1, 9, 13, 17, 18, 20, 28], to explore reaction-diffusion systems from a perspective of dynamical systems. Majority of cases with analytical approaches concluded with certain conditions [9, 17, 18, 20, 28, 31] on the characteristics of the stability matrix that theoretically predicts reaction-diffusion systems to exhibit spatial patterns, lacking to extend the analysis to computational consequences of these conditions on the admissible parameter spaces in terms of diffusion-driven instability. Furthermore, from the literature to date [28, 30, 39, 45, 46, 47, 48], it is evident that insufficient attention is given to rigorously explore the influence of reaction-diffusion rates on the dynamical behaviour of such systems in the context of domain size. The motivation of the current work originates from the hypothesis of the spatial dependence of the eigenfunctions of the diffusion operator with the domain size itself and extending this idea to further investigate how this spatial dependence induces an influence on the linearised stability-matrix of a reaction-diffusion system through the properties of the eigenfunctions of the diffusion operator. With an attempt to further contribute to the current knowledge of reaction-diffusion systems, the present paper employs a set of rigorous findings obtained from bifurcation analysis of the activator-depleted reaction-diffusion model to explore the corresponding influence induced on the admissible parameter spaces and diffusion-driven instability. Despite the restriction of the current work to a particular type of reaction kinetics namely activator-depleted, the methodology can however serve as a framework for developing a standard independent method by combining the known aspects of exploring the topic namely bifurcation analysis, parameter spaces and computational methods, that could be utilised for general reaction-diffusion systems. The lack of a self-contained complete methodology that combines all the known distinct aspects of exploring this topic creates a credible argument for the importance of this work. The contents of the current paper offer a natural extension of the work presented in [38], where the domain of solution was restricted to rectangular geometries. Equivalent results to those presented in [38] are found in the present work, except that the domain of solution in this work is a two dimensional disk-shape geometry bounded by a circle. The bulk of the current work consists of rigorously proven statements that quantitatively relate the radius of a disk-shape solution domain to reaction-diffusion rates in light of bifurcation analysis, which are computationally verified using the finite element method.

The contents of this work are structured such that in Section 2 the model equations for activator-depleted reaction-diffusion system are stated in cartesian coordinates and transformed to polar coordinates in its non-dimensional form with the corresponding initial and boundary conditions. Section 3 consists of a detailed analytical method to explicitly derive eigenfunctions and the corresponding eigenvalues satisfying the boundary conditions prescribed for activator-depleted reaction-diffusion system in Section 2. Furthermore, in Section 3, an application of spectral approach on polar coordinates using chebyshev and Fourier grid method is presented [57] to demonstrate and verify the analytical derivation of the eigenfunctions as well as the stability matrix and the corresponding characteristic polynomial for a linearised approximation of the original system. Section 4 presents analytical results on the relationship of domain-size (radius of disk) with different types of bifurcations in the dynamics. Section 5 is devoted to the parameter space classification in light of bifurcation analysis of the uniform steady-state of the system. A numerical method for computing the partitioning curves for such classification is presented in detail. The shift of parameter spaces as a consequence of changing reaction-diffusion rates is investigated and theorems are proven that rigorously establish the relation of domain size (radius of disk) with the corresponding types of diffusion-driven instabilities. Furthermore, a full classification of the admissible parameter space in terms of stability and types of the uniform steady state of the system is also presented in Section 4, where it is shown that if the radius of a disk-shape domain satisfies certain inequalities in terms of reaction-diffusion rates, then the admissible parameter space allows or forbids certain types of bifurcations in the dynamics of the reaction-diffusion system. Section 6 presents numerical simulations of an activator-depleted reaction-diffusion system using finite element method to verify the proposed classification of parameter spaces and the theoretically predicted behaviour in the dynamics. Due to the curved boundary of the domain a non-standard technique called distmesh [40, 41, 44] is used to obtain discretisation for simulating the finite element method. The technicality of the algorithm for distmesh is briefly explained. Section 7 presents the conclusion and possible extensions of the current work.

## 2 Model equations

Two chemical species and are modelled by the well-known activator-depleted reaction-diffusion system, where both species are coupled through nonlinear reaction terms. The system assumes independent diffusion rates for both of the species. The RDS is considered on a two dimensional circular domain denoted by , where is defined by

 Ω={(x,y)∈R2:x2+y2<ρ2}.

This forms a two dimensional disk-shape domain with a circle forming its boundary denoted by , which contains the points given by

 ∂Ω={(x,y)∈R2:x2+y2=ρ2}.

The RDS with activator-depleted reaction kinetics in its non-dimensional form on cartesian coordinates has the form

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩{∂u∂t=△u+γ(α−u+u2v),(x,y)∈Ω,t>0,∂v∂t=d△v+γ(β−u2v),∂u∂n=∂v∂n=0,on(x,y)∈∂Ω,t≥0,u(x,y,0)=u0(x,y),v(x,y,0)=v0(x,y),(x,y)∈Ω,t=0, (1)

where , , and are strictly positive real constants. In system (1) the positive parameter denotes the non-dimensional ratio of diffusion rates given by , where and are the independent diffusion rates of and respectively. The non-dimensional parameter is known as the scaling parameter of (1), which quantifies the reaction rate. The boundary of under the current study assumes homogeneous Neumann boundary conditions (also known as zero flux boundary conditions), which means that the chemical species and can neither escape nor enter through . denotes the normal to in the outward direction. Initial conditions for (1) are prescribed as positive bounded continuous functions and . System (1) can be spatially transformed to polar coordinates using and to obtain

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩{∂u∂t=△pu+γf(u,v),∂v∂t=d△pv+γg(u,v),∂u∂r∣∣r=ρ=∂v∂r∣∣r=ρ=0,(r,θ)∈∂Ω,t≥0,u(r,θ,0)=u0(r,θ),v(r,θ,0)=v0(r,θ),(r,θ)∈Ω,t=0, (2)

where denotes the Laplace operator in polar coordinates written as

 △pu(r,θ)=1r∂∂r(r∂u∂r)+1r2∂2u∂θ2. (3)

In (2), and depend on coordinates and the functions and are defined by and . Initial and boundary conditions are transformed in a similar fashion. In (2) the strictly positive constants namely , , and remain to satisfy exactly the same definitions as in (1).

## 3 Stability analysis of the reaction-diffusion system

### 3.1 Stability analysis of the reaction-diffusion system in the presence of diffusion

Let and denote the uniform steady state solution satisfying system (1) (and equivalently system (2)) with activator-depleted reaction kinetics in the absence of diffusion and these are given by [50, 56, 38, 17]. For linear stability analysis, system (2) is perturbed in the neighbourhood of the uniform steady state and the dynamics of the perturbed system are explored i.e. , where and are assumed small. In system (2) the variables and are substituted by the expression in terms of and and expanded using Taylor expansion for functions of two variables up to and including the linear terms with the higher order terms discarded. This leads to the linearised version of system (2) written in matrix form as

 ∂∂t[¯u¯v]=[100d][△p¯u△p¯v]+⎡⎢⎣∂f∂u(us,vs)∂f∂v(us,vs)∂g∂u(us,vs)∂g∂u(us,vs)⎤⎥⎦[¯u¯v]. (4)

The next step to complete the linearisation of system (2) is to compute the eigenfunctions of the diffusion operator namely , which will require to find the solution to an elliptic 2 dimensional eigenvalue problem on a disk that satisfies homogeneous Neumann boundary conditions prescribed for system (2).

The solution of eigenvalue problems on spherical domains is a well studied area [58, 33, 59], with the majority of research focused on problems with boundary-free manifolds such as circle, torus and/or sphere. Considering the restriction imposed from the boundary conditions prescribed for the current problem, entails that the case requires explicit detailed treatment to rigorously find the eigenfunctions satisfying these boundary conditions. Therefore, it is important to present a step-by-step demonstration of the process, starting with writing out the relevant eigenvalue problem all the way through to finding the closed form solution in the form of an infinite set of eigenfunctions satisfying such an eigenvalue problem. The eigenvalue problem we want to solve is of the form

 {△pw=−η2w,η∈R,∂w∂r∣∣r=ρ=0,ρ∈R+∖{0}, (5)

where is the diffusion operator in polar coordinates defined by (3), on a disk with radius . We use separation of variables to obtain the solutions of (5). Application of separation of variables to problem (5) requires a solution of the form , which is substituted into problem (5) to obtain

 d2Rdr2Θ+1rdRdrΘ+1r2Rd2Θdθ=−η2RΘ. (6)

Dividing both sides of (6) by , multiplying throughout by and rearranging, yields

 r21Rd2Rdr2+r1RdRdr+η2r2=−1Θd2Θdθ2. (7)

Using anzats of the form and a change of variable is applied to (7). Taking the usual steps of Frobenius method it can be shown that the general solution satisfying (7) takes the form , where and are respectively given by

 R1(x)=xn∞∑j=0(−1)jC0x2j4j×j!×(n+j)×(n+j−1)×⋅⋅⋅×(n+1) (8)

and

 R2(x)=x−n∞∑j=0(−1)jC0x2j4j×j!×(−n+j)×(−n+j−1)×⋅⋅⋅×(−n+1). (9)

Series solutions (8) and (9) are referred to as the Bessel functions of the first kind [60, 59]. Before writing the general solution to problem (5), it must be noted that we require to be a function of and not of , bearing in mind that is related to under the linear transformation given by . Therefore, the general solution to problem (5) can be written in the form

 w(r,θ)=R(x(r))Θ(θ), (10)

where and , with and defined by (8) and (9) respectively. The homogeneous Neumann boundary conditions are imposed on the set of eigenfunctions (10), and noting that the flux is independent of the variable , we have which implies that . A straightforward application of chain rule yields .

Let and denote the coefficients corresponding the term in the infinite series for and respectively, then for every the expressions for and take the forms

 aj=(−1)jC04j×j!×(n+j)×(n+j−1)×⋅⋅⋅×(n+1), (11)
 bj=(−1)jC04j×j!×(−n+j)×(−n+j−1)×⋅⋅⋅×(−n+1). (12)

This entails that can be written in the form . Differentiating with respect to and equating it to zero, on substituting and using the chain rule we obtain the equation

 0=dRdr∣∣r=ρ=η∞∑j=0[aj(n+2j)xn+2j−1+bj(−n+2j)x−n+2j−1]∣∣x=ηρ=η∞∑j=0[aj(n+2j)(ηρ)n+2j−1+bj(−n+2j)(ηρ)−n+2j−1]=η[(ηρ)n∞∑j=0aj(n+2j)(ηρ)2j−1+(ηρ)−n∞∑j=0bj(−n+2j)(ηρ)2j−1], (13)

which holds true if and only if both of the summations in (13) are independently zero. Investigating the first summation in (13) and expanding it for a few successive terms, it can be shown that the successive terms carry alternating signs (due to the expression for ), therefore, the only way the infinite series can become zero is if each of the successive terms cancel one another. Let and denote the terms of the first and second summations in (13) respectively, then for (13) to hold true, and must be true for all . The full expressions for and can be written as

 Fj=(ηρ)n(−1)jC0(n+2j)(ηρ)2j−14j×j!×(n+j)×(n+j−1)×⋅⋅⋅×(n+1) (14)

and

 Fj+1=(ηρ)n(−1)j+1C0(n+2j+2)(ηρ)2j+14j+1×(j+1)!×(n+j+1)×(n+j)×⋅⋅⋅×(n+1). (15)

For the first , second and third successive pairs, independently equating the sum of the expressions (14) and (15) to zero yields that can be written as

 η2n,k=4(2k+1)(n+2k+1)(n+4k)ρ2(n+4k+2). (16)

#### Remark

The restriction on the order of the corresponding Bessel’s equation namely in Theorem 3.1 can be relaxed by employing Bessel’s function of the second kind [58, 59, 60] and imposing on it the homogeneous Neumann boundary conditions. Therefore, for the purpose of the current study the order of the corresponding Bessel’s equation is precluded from becoming a full or half integer. Bear in mind that the proposed choice of makes the set of eigenvalues a semi discrete infinite set. The spectrum is discrete with respect to positive integers and continuous (uncountable) with respect to .

###### Theorem 3.1

Let satisfy problem (5) with homogeneous Neumann boundary conditions. Given that the order of the associated Bessel’s function belongs to the set , [see Remark 3.1.1] then for a fixed there exists an infinite set of eigenfunctions of the diffusion operator as defined in (3), which is given by

 wn,k=[R1n,k(r)+R2n,k(r)]Θn(θ) (17)

with the explicit expressions for , and as

 R1n,k(r)=∞∑j=0(−1)jC0(ηkr)2j+n4jj!(n+j)(n+j−1)⋅⋅⋅(n+1), (18)
 R2n,k(r)=∞∑j=0(−1)jC0(ηkr)2j−n4jj!(−n+j)(−n+j−1)⋅⋅⋅(−n+1) (19)

and

 Θn(θ)=exp(inθ), (20)

where for the successive pair satisfies (16), when ever . {proof} The proof consists of all the steps from (6) to (16).

It is worth noting that for a set containing terms from the series solutions defining the eigenfunctions (17), one can only obtain eigenvalues that satisfy problem (5). Therefore, it is important to realise that unlike the eigenfunctions of a diffusion operator [38, 63] on a rectangular geometry, here only pairwise terms can qualify in the series solutions to be claimed as the eigenfunction satisfying problem (5) and not individual terms associated discretely to the summation index .

### 3.2 Numerical experiments using the spectral method

Let denote the expression , then the full set of eigenfunctions to problem (5), can be written as . For numerical demonstration of Theorem 3.1, the spectral method [57] is employed on a unit disk centred at the origin of the coordinates. A spectral mesh in polar coordinates is constructed on a unit disk , where a periodic Fourier grid is applied to the angular axis and a non-periodic chebyshev grid is applied to the radial axis [57]. In order to tackle the singularity at , chebyshev discretisation is applied to the whole of the diameter of , which means that is used instead of . The interval is discretised in the form , for , where , is a positive odd integer. The odd number of points on the chebyshev grid serves to preclude complications that could rise from the singular point by locating it radially in the middle of two successive chebyshev grid points. The second coordinate is discretised using , for , where , is a positive even integer. For the full details on the implementation of the spectral method in polar coordinates the reader is referred to [57]. Figure 3 (a) shows a coarse mesh structure on a unit disk using and , with angular step-size of . Similarly, Figure 3 (b) shows the refined mesh using and with angular step-size of , on which all the simulations for eigenfunctions are performed to visualise the corresponding nodal lines and surfaces. Colour encoded plots, corresponding to nine different modes are simulated on the mesh given in Figure 3 (b) using a technique presented in [35, 36]. Due to the fact that eigenfunctions are complex valued functions, therefore, trivial methods of visualising a real valued function of two variables do not suffice to give a meaningful representation to a complex valued function. Note that in the expression for the eigenfunctions , only the function contains the non-zero imaginary part, therefore, the variable is encoded by a polar colour scheme (Hue, Saturation Value) and displayed directly on . For full details on the implementation of this process the reader is referred to [35, 36]. Figure 13 shows a colour encoded representation of nine different modes, each of which corresponds to one of the nine nodal line depictions of in Figure 14.

Note that in Figure 14 different modes namely modes 2 and 3 result in the same eigenvalues that correspond to a different orientation of the nodal line depiction for and . Similarly, one may find that there are a lot of such pairs of positive integers and , that correspond to different orientations of the same nodal line depiction. This occurs, when the multiplicity of an eigenvalue is 2, and in fact there are numerous such pairs of and , that correspond to the same , but for different orientation of eigenmodes. However applying the colour encoded representation using the (HSV) colour scheme indicates a distinction between the plots corresponding to each integer value for as shown in Figure 13.

### 3.3 Stability matrix and the characteristic polynomial

The solution to system (2) using separation of variables, can be written (with bars omitted from and ) as the product of the eigenfunctions of the diffusion operator and in the form

 u(r,θ,t)=∞∑k=0Un,kexp(σn,kt)Rn,k(r)Θn(θ),v(r,θ,t)=∞∑k=0Vn,kexp(σn,kt)Rn,k(r)Θn(θ),

where and are the coefficients associated with the mode of the eigenfunctions in the infinite expansion. Substituting this form of solution in (4) and the steady state values in terms of the parameters and for , one obtains the fully linearised form of (2) as a system of two algebraic equations of the form

 Unknown environment '%' (21)

where . In turn, this can be written as a two dimensional discrete eigenvalue problem of the form

 ⎡⎢⎣γβ−αβ+α−η2n,kγ(β+α)2−γ2ββ+α−γ(β+α)2−dη2n,k⎤⎥⎦[uv]=σ[uv]. (22)

The left-hand matrix in (22) is referred to as the stability matrix [50, 56, 12] for system (2), with given by (16). In order to investigate the stability of the uniform steady state , it is required to analyse the eigenvalues satisfying (22), for which the characteristic polynomial takes the form of a quadratic equation in , written as

 ∣∣ ∣ ∣∣γβ−αβ+α−η2n,k−σγ(β+α)2−γ2ββ+α−γ(β+α)2−dη2n,k−σ∣∣ ∣ ∣∣=0. (23)

Let and respectively denote the trace and determinant of the stability matrix given by (22), then the quadratic polynomial (23) can be written in terms of and as

 σ2−T(α,β)σ+D(α,β)=0, (24)

with and expressed by

 ⎧⎪⎨⎪⎩T(α,β)=γβ−α−(β+α)3β+α−(d+1)η2n,k,D(α,β)=(γβ−αβ+α−η2n,k)(−γ(β+α)2−(d+1)η2n,k)+2γ2β(β+α). (25)

The roots of equation (24) in terms of and are . Stability of the uniform steady state is determined by the signs of the two roots namely if they are real. However, if are a complex conjugate pair (with non-zero imaginary parts), then the sign of the real part is sufficient to predict the stability of the uniform steady state . Since there are always exactly two roots on a complex plane for a quadratic equation, therefore, it is impossible with the same choice of parameters , , , and , for to be real and to be complex or vice versa. Therefore, a reasonable approach to encapsulate all the possibilities for the stability and types of the uniform steady state in light of parameters and is to consider the cases when and . In each case the parameter space is rigorously analysed and the classification of the parameter plane in relation to the diffusion parameter is studied. In light of such classification, the analysis is further extended to explore the effects of domain size on the existence of regions in parameter space that correspond to spatial and/or temporal bifurcations.

## 4 Parameter spaces and bifurcation analysis

Bifurcation analysis of system (2) is better conducted when the parameter plane is appropriately partitioned for both of the cases when as well as when . To obtain such a partition on the parameter plane, it is required to find the equations of the partitioning curves and these can be found through a detailed analysis of the expression for , which in turn requires to explore the domains of and .

### 4.1 Equations of the partitioning curves

Starting with the curve on the parameter plane that forms a boundary for the region that corresponds to eigenvalues containing non-zero imaginary part. It must be noted that the only possibility through which can have a non-zero imaginary part is if the inequality is true. It means that those parameter values and satisfying the equation must be lying on a partitioning curve that determines the boundary between the region on the parameter plane that corresponds to eigenvalues with non-zero imaginary part and that which corresponds to a pair of real eigenvalues. We therefore, state that the set of points on the parameter plane satisfying the implicit equation

 Missing or unrecognized delimiter for \Big (26)

forms the partitioning curve between the region that corresponds to a real pair of and that corresponding to a complex conjugate pair of . For the solution of (26) refer to Section 5, where a numerical method is employed to find combinations of on the plane satisfying (26). It is worth noting that a combination of and satisfying (26), entails that the expression for possesses a repeated real root of the form , whose sign will require to be explored in Section 5 for real eigenvalues. Another important point to be made about the curve satisfying (26) is that it also partitions the parameter space for regions of spatial and temporal bifurcations. Because on one side of curve (26), both of the eigenvalues are always real, which can never excite temporal instability of the uniform steady state , whereas on the other side, the eigenvalues are always a complex conjugate pair, which can never excite spatial instability. Therefore, on the side where are a pair of real values, any instability that occurs in the dynamics of system (2) will be strictly relevant to spatial variation, hence any pattern that system (2) can evolve to will be strictly spatial with stable and invariant evolution in time. However, on the side where have non-zero imaginary part, every possible instability in the dynamics of system (2) will be strictly concerned with temporal periodicity, therefore any pattern that emerges from the dynamics of system (2) is expected to be periodic along the time axis. This analysis raises the interesting question whether it is possible for the dynamics of (2) to cause a spatially periodic pattern to undergo an unstable temporal periodicity as well? To answer this question, on the first sight it sounds absurd to claim on one hand that a specific choice of admissible can either cause spatial or temporal instability, and not both at the same time, which is an intuitive claim to make. Because, any admissible choice of either yields a pair of real or a complex conjugate pair of , therefore, a fixed choice of is not expected to yield a real and a complex or vice versa. Therefore, it is intuitive to presume that system (2) should not admit the evolution of such dynamics, in which a combination of spatial and temporal instabilities can occur. However, the current study finds that this counter-intuitive behaviour is possible for system (2) to exhibit, which is related to the existence and shift in the location of a partitioning curve on which the real part of becomes zero during the course of the temporal evolution. Such a curve is referred to as the transcritical curve, for which the implicit equation in terms of admissible and is of the form

 γβ−α−(β+α)3β+α=(d+1)η2n,k, (27)

under the assumption that . It must be noted, that equations (26) and (27) are the only two equations that fully partition the parameter plane . However, one may expect the number of partitions subject to the types and stability of the uniform steady state , on the parameter plane to be four regions separated by three curves. These regions are obtained by sub-partitioning the region corresponding to a pair of real , into two sub-regions where is a pair of real negative values and another where at least or is positive. The region on the parameter plane , that corresponds to complex eigenvalues can also be sub-divided into two regions, one where is a complex conjugate pair with negative real part, and the other where is a complex conjugate pair, but with a positive real part. The numerical solution of equation (26) will reveal that the region corresponding to complex eigenvalues is in fact bounded by two curves in the parameter plane , each of which satisfies (26). Therefore, (26) is implicitly the equation of two partitioning curves instead of one, which including (27) makes a total of three curves dividing the admissible parameter space into four different regions. It is worth bearing in mind, that subject to the types and stability of the uniform steady state, one obtains at most four regions separated by three curves, with the possibility that parameters such as , and may induce equations (26) and (27) in such a way that the number of partitioning curves may become less than three in total. Consequently this entails that a region corresponding to a certain type of bifurcation may completely disappear from the admissible parameter space. This kind of influence on the location and existence of the partitioning curves is quantitatively investigated, in particular, the effect of embedded in the expression for in (16) is explored to analyse its influence on the bifurcation of the uniform steady state.

#### Analysis for the case of complex eigenvalues

Before determining the region with complex eigenvalues on the admissible parameter plane using a numerical treatment of (26), the real part of is investigated analytically, when it is a complex conjugate pair. It can be noted that can only become a pair of complex roots, if satisfies the inequality

 T2(α,β)−4D(α,β)<0. (28)

Given that (28) is satisfied, then the stability of the uniform steady steady state is decided purely by the sign of the real part of , which is the expression

 Re(σ1,2)=12(γβ−α−(β+α)3β+α−(d+1)η2n,k). (29)

If the sign of the expression given by (29) is negative, simultaneously with assumption (28) satisfied, then it can be predicted that no choice of parameters can cause temporal instability in the dynamics of system (2). Therefore, under the assumption (28), if the dynamics of system (2) do exhibit diffusion-driven instability, it will be restricted to spatially periodic behaviour only, which uniformly converges to a temporal steady state, consequently one obtains spatial pattern that is invariant in time. The sign of the expression given in (29) is further investigated to derive from it, relations between the parameter controlling the domain size and reaction-diffusion rates denoted by and respectively. Given that assumption (28) is satisfied then the sign of expression (29) is negative if parameters , , and satisfy the inequality

 β−α−(β+α)3β+α<(d+1)η2n,kγ, (30)

with defined by (16). Note that the expression on the left hand-side of (30) is a bounded quantity by the constant value of 1 [38], for all the admissible choices of , therefore, substituting for in expression (16) and rearranging, it can be obtained that for inequality (30) to remain true, it induces a restriction on the value of , which is of the form

 ρ2<4(d+1)(2k+1)(n+2k+1)(n+4k)γ(n+4k+2). (31)

Inequality (31) conversely implies that so long as the radius of the disk shape domain satisfies (31), then the dynamics of system (2) is guaranteed to exhibit global temporal stability, which also means that any possible instability in the dynamics must be restricted to spatial periodicity or spatial pattern. The formal proof of this claim is presented in Theorem 33. This type of instability concerning space and not time is referred to as Turing instability [50, 56]. On the other hand assuming that (28) is satisfied and using the upper bound of the quantity on the left hand-side of (30) with as defined in (16), it can be shown that a necessary condition for the sign of expression (29) to become positive is for to satisfy the inequality

 ρ2≥4(d+1)(2k+1)(n+2k+1)(n+4k)γ(n+4k+2). (32)

Conditions (31) and (32) both have quantitative influence on the location and topology of the partitioning curves obtained from the numerical solutions of (26) and (27) in the admissible parameter plane, namely . System (2) is restricted from any type of temporal bifurcation if the radius of the disk shape domain is related to the reaction-diffusion parameters and , through inequality (31). It also means that on a disk shape domain with radius satisfying (31), the dynamics of system (2), either exhibit spatially periodic pattern or no pattern at all, both of which are globally stable in time. If the dynamics of a system become unstable along the time axis and exhibits temporal periodicity, then the system is said to undergo Hopf bifurcation [61, 50, 39, 2]. If the real part of a pair of complex eigenvalues become zero, then the system is expected to exhibit oscillations with orbital periodicity. This behaviour is known as transcritical bifurcation [50, 2, 39]. The consequences of the restriction on in the sense of bifurcation analysis means that, whenever the radius of a disk-shape domain is bounded by (31) in terms of and , then the dynamics of system (2) is guaranteed to forbid Hopf and transcritical bifurcations, only allowing for Turing instability to occur. If system (2) allows only Turing type instability to occur under condition (31), it indicates that the eigenvalues only become positive, when they are a pair of real values with zero imaginary parts. Therefore, it is also an indication that diffusion-driven instability is still possible, but it just becomes strictly spatial with (31) satisfied. If the values of parameters and are chosen such that inequality (32) is satisfied, then the possibility of all three types of diffusion-driven instabilities exist on the admissible parameter plane , namely Turing, Hopf and transcritical types of bifurcations. The influence from the area of through the relationship (32) of with the reaction-diffusion rates namely and is summarised in Theorem 33 with a detailed sketch of the proof.

###### Theorem 4.1 (Hopf or transcritical bifurcation)

Let and satisfy the non-dimensional reaction-diffusion system with activator-depleted reaction kinetics (2) on a disk-shape domain with radius and positive real parameters , , and . For the system to exhibit Hopf or transcritical bifurcation in the neighbourhood of the unique steady state , the necessary condition on the radius of the disk-shape domain is that it must be sufficiently large satisfying

 ρ≥2√(d+1)(2k+1)(n+2k+1)(n+4k)γ(n+4k+2), (33)

where is the associated order of the Bessel’s equations and is any positive integer.

{proof}

[Proof:] For system (2) to exhibit Hopf or transcritical bifurcations the eigenvalues of the stability matrix (22) must have non-zero imaginary part with non-negative real part. Consider the real part of , which is precisely given by under the assumption that the admissible choice of parameters satisfies (28). When , then the stability of the uniform steady state is precisely determined by the sign of , which is given by

 T(α,β)=γβ−α−(β+α)3β+α−(d+1)η2n,k. (34)

System (2) undergoes Hopf or transcritical bifurcation if , given that the strict inequality (28) is satisfied, which can only hold true if . In (34) is given by

 η2n,k=4(2k+1)(n+2k+1)(n+4k)ρ2(n+4k+2), (35)

where is the order of the associated Bessel’s equation and is any positive integer. To show the condition on for Hopf or transcritical bifurcation, one may substitute (35) into (34) and requiring the resulting quantity to be non-negative, which yields the inequality

 γβ−α−(β+α)3β+α≥4(d+1)(2k+1)(n+2k+1)(n+4k)ρ2(n+4k+2). (36)

Noting that the left hand-side of (36) can be written as the difference between two non-negative functions and in the form , where and are given by

 f1(α,β)=βα+β,f2(α,β)=α+(α+β)3α+β. (37)

Note also that resides in the denominator of the right hand-side of (36) and parameter is multiplied by the expression on the left hand-side. In order to find what this inequality induces on the relationship between parameters , and , it is essential to analyse the supremum and infemum of and within their respective domains which is . The range for and are independently analysed to find the supremum of the expression on the left of (36). Starting with , which is bounded below and above in the domain , we have , and the for all Similarly considering the expression for , we have and the , for all Since the ranges of both and are non-negative within their respective domains, therefore the supremum of their difference is determined by the supremum of the function with positive sign, which is . Therefore, inequality (36) takes the form

 4(d+1)(2k+1)(n+2k+1)(n+4k)ρ2(n+4k+2)≤γβ−α−(β+α)3β+α≤γsupα,β∈R+(f1(α,β)−f2(α,β))=γsupα,β∈R+f1(α,β)=γ,

which by rearranging and writing the inequality for in terms of everything else, yields the desired statement of Theorem 33, which is condition (33). The claim of Theorem 33, is also numerically verified by showing that a region in the admissible parameter plane that corresponds to Hopf or transcritical bifurcations emerges only if radius of a disk shape domain is sufficiently large satisfying inequality (33). Otherwise, no choice of parameters exist in the admissible parameter plane , allowing the dynamics of (2) to exhibit Hopf or transcritical bifurcation.

#### Analysis for the case of real eigenvalues

The eigenvalues are both real if the discriminant of the roots is either zero or positive, which in turn means that the relationship between and is such that

 T2(α,β)≥4D(α,β). (38)

The equal case of (38) is looked at first, where we have

 T2(α,β)=4D(α,β), (39)

which implies that the discriminant is zero, hence the roots are repeated real values of the form , given by

 σ1=σ2=12(γβ−α−(β+α)3β+α−4(d+1)(2k+1)(n+2k+1)(n+4k)ρ2(n+4k+2)). (40)

When and satisfy condition (39), the stability of the steady state is determined by the sign of the root itself. The expression given by (40) can be easily shown to be negative if the radius of the disk-shape domain satisfies the inequality

 ρ<2√(α+β)(d+1)(2k+1)(n+2k+1)(n+4k)γ(β−α−(α+β)3)(n+4k+2). (41)

Otherwise, the repeated root is positive provided that satisfies

 ρ>2√(α+β)(d+1)(2k+1)(n+2k+1)(n+4k)γ(β−α−(α+β)3)(n+4k+2). (42)

Analysing (41) and (42) carefully, it can be observed that the only terms that can possibly invalidate the inequalities are in the denominator of the right hand-side, namely the expression . Therefore, a restriction is required to be stated on this term to ensure that the radius of is not compared against an imaginary number, such a restriction is

 β>α+(β+α)3. (43)

It must be noted that (43) is the same restriction on the parameter choice obtained for the case of repeated real eigenvalues in the absence of diffusion [38]. By further comparing with (42), it can be noted that it is very similar to condition (33) of Theorem 33, except that (33) is free from any dependence of the parameters and . This makes (33) a sharper version of (42) in the sense, that the curve satisfying (39) subject to condition (42) must be the one forming the boundary of the region on the admissible parameter plane, that corresponds to complex eigenvalues with positive real parts, which is the region for Hopf bifurcation. Therefore, the region of the admissible parameter plane that corresponds to Hopf bifurcation is on one side bounded by curve (39) and on the other side it is bounded by the curve satisfying (27) under their respective assumptions. In Section 5 it is verified to be the case by the numerical computation of the partitioning curve (26), which is the same curve (39). This analysis motivates to explore the possibility of similar comparison between the conditions (41) and (31). A reasonable intuition behind this comparison is that the sub-region on the admissible parameter plane that corresponds to complex eigenvalues with negative real parts must be bounded by curve (39) subject to condition (41), outside of which every possible choice of parameters and will guarantee the eigenvalues to be a pair of distinct real values, which promotes the necessity to state and prove Theorem 44.

###### Theorem 4.2 (Turing type diffusion-driven instability)

Let and satisfy the non-dimensional reaction-diffusion system with activator-depleted reaction kinetics (2) on a disk-shape domain with radius and positive real parameters , , and . Given that the radius of domain satisfies the inequality

 ρ<2√(d+1)(2k+1)(n+2k+1)(n+4k)γ(n+4k+2), (44)

where is the associated order of the Bessel’s equations and