Classification of parameter spaces for reaction-diffusion systems on stationary domains
Wakil Sarfaraz, Anotida Madzvamuse
TL;DR
This paper classifies parameter spaces for reaction-diffusion systems on stationary domains, analyzing how domain size, reaction, and diffusion rates influence system stability and bifurcations, supported by theoretical and numerical results.
Contribution
It provides a comprehensive classification of parameter spaces for reaction-diffusion systems, linking domain size and parameters to different bifurcation types, with new theoretical bounds and conditions.
Findings
Domain size bounds determine bifurcation types.
Turing instability occurs within specific domain size conditions.
Numerical simulations support theoretical predictions.
Abstract
This paper explores the classification of parameter spaces for reaction-diffusion systems of two chemical species on stationary domains. The dynamics of the system are explored both in the absence and presence of diffusion. The parameter space is fully classified in terms of the types and stability of the uniform steady state. In the absence of diffusion the results on the classification of parameter space are supported by simulations of the corresponding vector-field and some trajectories around the uniform steady state. In the presence of diffusion, the main findings are the quantitative analysis relating the domain-size with the reaction and diffusion rates and their corresponding influence on the dynamics of the reaction-diffusion system when perturbed in the neighbourhood of the uniform steady state. Theoretical predictions are supported by numerical simulations both in the…
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 10
Figure 11
Figure 12
Figure 13
Figure 14
Figure 15
Figure 16
Figure 17
Figure 2
Figure 7
Figure 8
Figure 9
Figure 22
Figure 23
Figure 24
Figure 25
Figure 26
Figure 27
Figure 28
Figure 29
Figure 30
Figure 31
Figure 32
Figure 33
Figure 34
Figure 35
Figure 36
Figure 37
Figure 38
Figure 39
Figure 40Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsNonlinear Dynamics and Pattern Formation · Mathematical and Theoretical Epidemiology and Ecology Models · Quantum chaos and dynamical systems
Classification of parameter spaces for reaction-diffusion systems on stationary domains
Wakil Sarfaraz11footnotemark: 1
University of Sussex, School of Mathematical and Physical Sciences, Department of Mathematics, Pevensey 3, Brighton, BN1 9QH, UK
Anotida Madzvamuse
Abstract
This paper explores the classification of parameter spaces for reaction-diffusion systems of two chemical species on stationary domains. The dynamics of the system are explored both in the absence and presence of diffusion. The parameter space is fully classified in terms of the types and stability of the uniform steady state. In the absence of diffusion the results on the classification of parameter space are supported by simulations of the corresponding vector-field and some trajectories around the uniform steady state. In the presence of diffusion, the main findings are the quantitative analysis relating the domain-size with the reaction and diffusion rates and their corresponding influence on the dynamics of the reaction-diffusion system when perturbed in the neighbourhood of the uniform steady state. Theoretical predictions are supported by numerical simulations both in the presence as well as in the absence of diffusion. Conditions on the domain size with respect to the diffusion and reaction rates are related to the types of diffusion-driven instabilities namely Turing, Hopf and Transcritical types of bifurcations. The first condition is an upper bound on the area of a rectangular domain in terms of the diffusion and reaction rates, which forbids the existence of Hopf and Transcritical types of bifurcations, yet allowing Turing type instability to occur. The second condition (necessary) is a lower bound on the domain size in terms of the reaction and diffusion rates to give rise to Hopf and Transcritical types of bifurcations as well as Turing instability.
keywords:
Reaction-diffusion systems , Dynamical systems, Bifurcation analysis, Stability analysis, Turing diffusion-driven instability , Hopf Bifurcation , Transcritical bifurcation , Parameter spaces
††journal: arxiv.org
1 Introduction
Reaction-diffusion systems (RDSs) attract a significant degree of attention from researchers in applied mathematics [1, 2, 3, 4, 5, 6, 7, 8, 9], mathematical and computational biology [10, 11, 12, 13, 14], chemical engineering [15, 16, 17] and so forth. Alan Turing was one of the first scientists to realise the significance of RDSs as a self-governing dynamical system [18], and showed that RDSs can be responsible for the emergence of spatial patterns in nature. A large number of scientists [1, 2, 11, 8, 19, 20, 21], since the publication of [18], have contributed to investigating RDSs, with various types of reaction kinetics. The most popular models of reaction kinetics explored in the literature are the activator-depleted model (also known as the Schnakenberg reaction kinetics) [20, 22, 23, 2], Meinhardt [12, 13] and Thomas [24] reaction kinetic models. From a research perspective the study of RDSs is conducted through different types of approaches, one of which focuses on the local behaviour of the dynamics of the RDS near a uniform steady state, which in turn relates to the subject of stability analysis [8, 6, 25, 26, 20, 12] of RDSs using the stability matrix. Linear stability analysis offers a great deal of insight regarding the behaviour of RDSs in the neighbourhood of a uniform steady state of a particular system. The other usual approach is numerical computation of the actual solution of RDSs using finite element, finite difference and other numerical methods [1, 2, 4, 5, 7, 9, 27, 28]. Numerical solution of RDSs helps to visualise the evolution of the dynamics in time. Numerical computations of RDSs with non-linear reaction kinetics with a particular choice of parameters partially encapsulates the possible types of dynamics that a certain RDS can exhibit. In order to better understand the behaviour of RDSs, it is necessary to have prior knowledge on the classification of the parameter values. Madzvamuse et al., in [20, 26] found regions of parameter space, corresponding to diffusion-driven instability with and without cross-diffusion respectively using the well-known activator-depleted reaction kinetics. Their approach to finding unstable regions in the parameter space is restricted to Turing instability only. One of the few complementary contributions from the present work is the application of a numerical method (exclusive to this paper) in order to obtain the full classification of parameter space. The numerical method is also employed to solve the equations for the implicit curves forming the partitioning of the classification within the parameter space. Liu et al., in [29] attempted to find constraints on the parameters of RDSs with activator-depleted reaction kinetics that causes the system to exhibit Hopf and transcritical bifurcations. The proofs in [29] are focused on the existence of bifurcation points for some theoretical constraints of parametrised variables of the actual parameters, with no relation between the domain size and reaction-diffusion rates. Comparing their work to the present study, our results are robust in the sense that we explicitly relate domain size to the reaction-diffusion rates. Using this relationship, the parameter space is classified for different types of bifurcations. Moreover, in the present work the parameter constraints are a consequence of the relationship between the domain size and the reaction-diffusion rates. An additional drawback in the analysis of [29] is that their results are produced on parametrised variables of the model and not on the actual parameters of the equations, which makes their results applicable to non-realistic possibilities (negative values) of the actual parameters of the model. This drawback is effectively resolved in the current work as the analysis is conducted on the actual two-dimensional positive real parameter space, which in addition to confirming the existence of different bifurcation regions, it offers concrete quantitative classification of the parameter space that guarantees the dynamics of RDSs to exhibit these bifurcations.
The majority of RDSs in the literature [2, 3, 9, 19, 21, 30] that exhibit spatial or temporal pattern, contain nonlinear terms in their reaction kinetics, which makes the mathematical and numerical analysis of such systems extremely challenging. With no closed analytical solutions, studies of the local behaviour of the systems are generally conducted by use of linear stability theory close to bifurcation points. Here, the behaviour of a system can be theoretically predicted in the neighbourhood of the uniform steady states. Linear stability analysis can help derive certain conditions on the reaction kinetics, which lead to causing instability in the dynamics of RDSs in the presence of diffusion. Numerous papers [6, 8, 23, 26, 20] exist in literature that routinely apply linear stability analysis to RDSs. The majority of the published work [6, 23] focuses on deriving linear stability conditions in terms of the characteristics of the stability matrix for diffusion-driven instability, lacking to explore what the numerical application of these conditions induce on the admissible choice of the parameter space for a certain RDS. Spatio-temporal pattern formation occurs when an RDS undergoes diffusion-driven instability [18, 22, 31]. This occurs when a uniform steady state which is stable in the absence of diffusion becomes unstable upon adding diffusion to the system. Diffusion-driven instability crucially depends on the values of diffusion and reaction rates of each reactant, however, more importantly it depends on the parameter choice of the reaction kinetics. In the current work the existing knowledge on the conditions for diffusion-driven instability in the literature is extended using a series of analytical and numerical techniques, to obtain new insights on the combined effects of diffusion and reaction rates, and in turn relating these to domain size of the evolution of the pattern. The detailed and quantitative analysis on the relationship between the domain size and diffusion-reaction rates in light of the parameter classification is an aspect that has not yet to-date received sufficient attention in the literature. The usual approach in selecting parameters for numerical computations of RDSs [31, 14] is based on the behaviour of RDSs in the absence of diffusion by use of trial and error or is based on previously published work, to observe instability caused by diffusion [1, 2, 6, 9, 19]. The absence of a robust method to fully classify the parameter space for an RDS, creates an arguable platform for the importance of this work. Efficient analytical as well as computational methods are used to demonstrate the quantitative relationship between the domain size and diffusion rate for an activator-depleted RDS. The main findings of the present work, which relate the domain size to the diffusion and reaction rates, are presented in the form of theorems with rigorous mathematical proofs and these theoretical results are supported computationally by finite element numerical solutions corresponding to the activator-depleted RDS on fixed rectangular domains. For each numerical demonstration the relative error plots of the solutions for each successive time-step are presented to visualise the convergence of the numerical approximate solutions.
This article is therefore structured as follows. In Section 2 we carry out detailed theoretical linear stability analysis of the system (1) in the absence of diffusion. Linear stability analysis is conducted by computing the stability matrix through which, the non-dimensional parameter space is derived and classified. Section 2.5.1 presents the methodology to compute the solutions of the partitioning curves for the classification on the parameter space. A combination of analytical and numerical methods using polynomials are applied. In Section 3, the linear stability analysis of the system is conducted in the presence of diffusion and the parameter space is explored to understand the consequences of including diffusion in the analysis. The parameter space is further explored to find the change in regions of the parameter space through varying the non-dimensional diffusion coefficient. Section 4 contains the main findings of this work which are presented in the form of two theorems. Each theorem is supported by finite element solutions of the model system. Section 5 presents conclusions, future directions and possible extensions of the current work.
2 Model equations
The dynamics of two chemical species and in a coupled system of nonlinear parabolic equations is considered. Both species diffuse with independent rates and satisfy the well-known activator-depleted reaction-diffusion system in a closed two-dimensional rectangular domain denoted by with area , where and are the corresponding side lengths in the direction of and axes respectively. The reaction-diffusion system satisfying the activator-depleted model for and in its non-dimensional form [20, 22, 23, 2] reads as
[TABLE]
where , , and are non-dimensional positive constants. In (1) the non-dimensional parameter denotes the quantity , where is the diffusion rate of the variable and is the diffusion rate of the variable . The non-dimensional parameter denotes the reaction rate, which is also known as the scaling parameter for the reaction kinetics. The boundary is subject to zero flux condition, which means on the chemical species and satisfy homogeneous Neumann boundary conditions in the form
[TABLE]
where denotes the outward normal through . For initial conditions the existence of some strictly positive quantity from each of these chemical concentrations in the domain is assumed, which is written as
[TABLE]
2.1 Remark
Despite the fact that all the results obtained in this paper hold for general rectangular geometries of i.e. with possibilities , however for simplicity in the analysis is considered as a square, which means . The results can be readily extended to a rectangular case by taking the area of as , where .
2.2 Stability analysis in the absence of diffusion
In the absence of diffusion the system (1) takes the form of a set of ordinary differential equations of the form
[TABLE]
To analyse the stability of system (4), it is necessary to compute its uniform steady state solution. Let denote the uniform steady state solution of the system (4), then must simultaneously satisfy the system of nonlinear algebraic equations in the form
[TABLE]
The nonlinear algebraic system (5) admits a unique solution in the form
[TABLE]
which enforces a restriction on the parameters of the system such that . Since both of these parameters resemble physical quantities, therefore, strictly positive, an appropriate interpretation of this restriction is that they both simultaneously cannot become zero.
2.3 Stability matrix
The stability of system (4) is analysed by computing the Jacobian matrix [22] of (5) and conducting the stability analysis using the uniform steady state solution , hence
[TABLE]
Substituting the expressions in terms of and for in the matrix becomes
[TABLE]
which is called the stability matrix [22, 32] for system (5).
2.4 Parameter analysis
Let
[TABLE]
denote the trace and determinant of respectively, then the characteristic polynomial for the eigenvalues of in terms of and takes the form
[TABLE]
Hence, the two roots of this characteristic polynomial in terms of and are given by
[TABLE]
Expression (6) for the eigenvalues is studied through investigating the domain of and , which is the positive real cartesian plane . The classification that would cause the uniform steady state to change stability and type due to the selection of the choice of the parameter values is explored by examining the sign of the real part [22, 6, 23, 26] of . For example the parameter space that makes the uniform steady state stable, is the simultaneous combined choice of that ensures the real part of the eigenvalues to be negative, which in turn is related to the discriminant of the roots expressed by (6). The full parameter space is investigated, so that all the possible types of influences due to the choice of parameters and on the stability and types of the uniform steady state are encapsulated. The classification is conducted based on to be a complex conjugate in the first case, then in the second case the parameter space is analysed when are real roots. In each case, the space is further classified into stable and unstable regions. In addition to this the partitioning curves, on which the steady state changes its type, are studied using a numerical technique that is employed exclusively and for the first time in the context of the present topic.
2.5 Analysis for the case of complex eigenvalues
It is clear that the eigenvalues are a complex conjugate pair if and only if the discriminant is negative, which means the parameters must satisfy the inequality
[TABLE]
In order to find what region is satisfied by (7), we must find the critical curves on which the expression on the left of (7) is equal to zero, which means that the discriminant changes sign by moving across these curves in the plane . These curves can be found by solving the equation , which is true for the choice of , satisfying
[TABLE]
Solving (8) implies finding the implicit curves that represent the critical region on the plane, for which the discriminant in the expression for is zero. The solution to (8) provides the boundaries, for the region of the plane that results in to be a complex conjugate pair. The choice of on the curves satisfying (8) enforces the eigenvalues of the system to be repeated real values, therefore, on these curves the steady state becomes a star node, whose stability will be analysed in the Section 2.6 with real eigenvalues. The left hand-side of (8) can be factorised in the form , which provides the equations of the two implicit curves and that determine the boundaries of the region corresponding to complex eigenvalues , where and are respectively given by
[TABLE]
2.5.1 Solving the partitioning curves
A mesh is constructed on a square domain , which is discretised by points in both directions of and , where is a positive integer. This constructs a square mesh of cells, each of size , with points in . To find the implicit solution for (9), at every mesh point in the direction , the roots of the cubic polynomial in namely , denoted by are computed using the Matlab command ‘’. It is worth noting that for every fixed , one obtains , to be a cubic polynomial of degree 3 in of the form
[TABLE]
where the coefficients are given by , , , and Similarly for each in the direction of , there are cubic polynomials in satisfying , with of the form
[TABLE]
where , , , and . Each of the equations and have at most three roots namely for every fixed , which means for every fixed the three points namely for are the three points that lie on the implicit curve given by (9). Since are positive real parameters in the system, therefore, at every fixed , only the positive real roots are considered and any root that is either real negative or complex is ignored.
Solving (8) alone serves to determine the boundaries of the complex region, but does not tell us, which side of these implicit curves correspond to real eigenvalues and which side corresponds to complex eigenvalues. This can be decided, by finding another curve on which the eigenvalues are purely imaginary with zero real parts. Since the real part of the eigenvalues is given by and the determinant is strictly positive, therefore, the choice of that solves the equation
[TABLE]
will ensure that the eigenvalues of the system are purely imaginary. Solving (10) is equivalent to finding the set of such that the equation is true, where is given by
[TABLE]
with coefficients and
The curve satisfying equation (10) forms the boundary for the positive and negative real parts of the eigenvalues when they are complex conjugate pair, because on the real part of the eigenvalues are zero and hence changes the sign by moving across it. The positive real combination of that satisfies (10) is given by the red curve in Figure 1. The location of the curve also indicates that the region between the curves and results in that are complex numbers. If are chosen from values on curves and , then the eigenvalues will be repeated real roots, since on these curves the discriminant is zero. However if are chosen from the values on the curve then the eigenvalues are purely imaginary of the form . This implies that the steady state is a centre with the system exhibiting periodic oscillations around the uniform steady state. It is worth noting that is a sufficient condition for to be a centre, due to the strictly positive expression for in the discriminant. Figure 1, was simulated for and , which means the curve , and are respectively formed from the positive real roots of cubic polynomials of , and for .
The curves and in Figure 1 (a) are the critical boundaries, at which the type of the steady state changes from node to a spiral or vice versa, depending on whether one enters or exists the region between and . The real part of is studied to find the stability classification of the parameter space, when are complex conjugate pair. Knowing that the region bounded by and curves in Figure 1 (a) results in to be complex, the next step is to find which part of this region corresponds to positive real part of the complex numbers , and which part of this region corresponds to negative real part of . In order to determine their stability we must analyse the real part of , simultaneously with the assumption that the roots are a pair of complex conjugate numbers, which is given by
[TABLE]
If the sign of (11) is negative in the region between the curves and in Figure 1 (a), then is a stable spiral, if it is positive then it is an unstable spiral. This means that the steady state is a stable spiral if simultaneously satisfy (7) and the inequality holds, which corresponds to the green region A shown in Figure 1 (b). Similarly, the steady state is an unstable spiral if parameters simultaneously satisfy (7) and the inequality , which is the yellow region B shown in Figure 1 (b). It is worth noting that the trace for positive values of and is bounded by the value of 1 [see section 4. Theorem 2], which means for all .
2.5.2 Numerical demonstration of complex eigenvalues
The numerical package for ordinary differential equations ode45 in MATLAB was employed to simulate and visualise the local behaviour of the system (4) near the uniform steady state . Each simulation was run up to a final time , starting at the initial time , with time step . Parameter values were chosen from the regions , and on the curve from Figure 1 (b) and Figure 1 (a) respectively, which are then plotted with 8 trajectories around the uniform steady state . Figure 2 shows that the behaviour of the trajectories around the uniform steady state in each case is in agreement with the theoretically predicted type of . The summary of the chosen numerical values for these simulations are all summarised in Table 1 and these are also indicated in the relevant captions for each figure. The eigenvalues of the stability matrix for each choice of the parameters and at the corresponding uniform steady state are presented in the corresponding captions. In each simulation the vertical axes corresponds to the variable and the horizontal axes corresponds to the variable . The red trajectories show the interactive behaviour of and in relation to the corresponding steady state.
2.5.3 Interpretation of the dynamics for the case of complex eigenvalues
The physical interpretation of this classification is that, whenever the choice of and is taken from region in Figure 1 then the dynamics of the perturbed system near the constant steady state will spirally move away as time grows, which resembles the shape of an outward spreading spiral. This behaviour is shown in Figure 2 (b). However if the parameter values are chosen from region then the dynamics of the perturbed system will spirally move in towards the constant stable steady state , therefore, any trajectory near the steady state forms the shape of spiral that is uniformly moving towards its centre, which is the constant steady state namely the point . This behaviour is shown in Figure 2 (a). If the values of the parameters are chosen from those on the curve , then the dynamics of the system behaves in such a way that it neither goes away from the steady state, and nor does it reach the steady state, instead the trajectories of the system move on fixed orbits around the uniform steady state , forming either circles or ellipse. This behaviour of and is shown in Figure 2 (c).
2.6 Analysis for the case of real eigenvalues
When both eigenvalues of the system are real numbers at the steady state then the uniform steady state becomes a node. If is a positive repeated root, then the steady state is an unstable star node. If is a negative repeated root, then the steady state is a stable star node. Parameter space resulting in to be a node can be analysed through studying the sign of , when both eigenvalues are real. This consequently means that the discriminant has to satisfy the inequality
[TABLE]
Similarly, to the previous section, treating the equal case of (12) first will provide the boundaries between the real and complex regions for . These are and curves in Figure 1 from the previous section. On and it was concluded that are repeated real roots, which are precisely . In order to determine the stability of , when are on and blue curves in Figure 1, it is required to find the classification of this curve on which the sign of the repeated eigenvalue is positive or negative. The steady state is a star, if is repeated real values, which means that the discriminant must be zero
[TABLE]
The steady state is a stable star if the choice of parameters simultaneously satisfy (13) and the inequality which is the curve in Figure 3 (a) in blue colour. The steady state is unstable star if the choice of parameters simultaneously satisfy (13) and the inequality which is the curve in Figure 3 (a) in red colour.
The remaining region outside and curves in Figure 3 (a) can be classified into two parts, one where satisfy (12) and yet both are distinct negative real values. This means that must satisfy the inequality
[TABLE]
such that both eigenvalues are negative real. This is true if satisfy the inequality
[TABLE]
The region in the plane satisfying (14) is denoted by A in Figure 3 (b) shaded in green colour. The remaining region to analyse is the region outside the curve in Figure 3 (a). Any combination of from this region corresponds to the remaining 2 cases, in which either one of or is positive, causing to be a saddle point which is unstable by definition. The other case is if are both positive real values, which makes to be an unstable node. The region satisfying the criteria for a saddle point or unstable node is indicated by B in Figure 3 (b) in yellow colour. It is worth noting that the steady state with either one or both positive real eigenvalues is by definition a saddle point or an unstable node, respectively. Under the current classification these both types fall under one category namely unstable node. The fact that one of the eigenvalues is real positive, hence as time grows very large, the behaviour of the solution is similar to that of an unstable node.
2.6.1 Numerical demonstration for the case of real eigenvalues
The system (4) is numerically simulated using ode45 package in MATLAB to visualise the behaviour of trajectories in the neighbourhood of , when are real values. The final time and time steps are chosen exactly the same as for complex eigenvalues case in the previous section. Each simulation is tested on phase plane diagram, to observe the trajectories of solution. In each case the steady state is computed and 8 trajectories around the steady state whose behaviour is governed by the relevant vector-field are plotted.
2.6.2 Interpretation of the dynamics for the case of real eigenvalues
If the values of and are selected from those on the curves or then the dynamics of the perturbed system will either radially move away or towards the point of the constant steady state , forming the shape of a star. If parameters are selected from the curve in Figure 3 (b), then the eigenvalues are positive real repeated roots. The dynamics of the perturbed system near the constant steady steady state will radially move outward (see Figure 4 (d)). If the values of and are chosen from those lying on the curve in Figure 3 (b), then the local dynamics of the system near the constant steady state radially move inward towards the point (see Figure 4 (c)). Similarly, if and are chosen from region in Figure 3 (b), it will result in to be a pair of distinct negative real roots. This means for the dynamics of the system to move towards the stable constant steady state. The trajectories near by the steady state in this case move towards a single point namely . If parameters are however chosen from region in Figure 3 (b), this results in the constant steady state to be an unstable node, it means the behaviour of trajectories near the steady state are expected to behaviour similarly to those corresponding to the curve (see Figure 4 (a)).
2.7 Summary of parameter classification in the absence of diffusion
In general, the parameter space that influences the nature of the uniform steady state can be categorised into four different regions, which are separated by three curves. This classification is in principle equivalent to that presented in [22], however the distinction with the current work is that for the specific system given by (4) the full parameter plane is classified subject to the proposed theory given in [22]. Each region is characterised and the behavior of the uniform steady state is established by choosing parameters and from the region. Figure 5 (a) summarises the full classification of the parameter space in the absence of diffusion for the uniform steady state . It can be noted that in terms of stability only and irrespective of the type of the steady state, the full parameter space namely the plane can be classified into two regions namely stable or unstable, this is the partition established by the location of the yellow curve in Figure 5 (a). Stability classification irrespective of the type of steady state is explicitly presented in Figure 5 (b). The purpose of stability partition will prove beneficial in the next section when the diffusion-driven instability is explored. Because diffusion-driven instability as the name suggests, is to explore how a stable steady state becomes unstable for the same choice of parameters , when diffusion is added to the system. Therefore, with the help of stability partition, one is able to observe the change in the location of the yellow curve in Figure 5 (b). Table 1 shows the values of the parameters and the values of the unique steady state with the relative eigenvalues that corresponds to the summary of full classification of parameter space presented in Figure 5 (a), and these parameter values are used in the numerical simulations.
3 Stability analysis in the presence of diffusion
It is intuitive to understand that normally diffusion serves to enhance spatial homogeneity of concentration gradients [22, 31], however, in the case of reaction-diffusion systems, it has the opposite effect in that the uniform steady state which was stable in the absence of diffusion becomes unstable in its presence [30, 33, 22, 31]. In this work, the mathematical rigor of how reaction and diffusion together can give rise to instability is investigated, which is responsible for diffusion-driven pattern formation. One of the explicit aims of the current work is to investigate the consequences of the conditions of diffusion-driven instability on the classification of the parameter space in the presence of diffusion. Furthermore, we want to explore quantitative relationship between the domain size and the types of diffusion-driven instability.
3.1 Linear stability analysis
The uniform steady state as previously mentioned satisfies the system with diffusion (1) and the zero-flux boundary conditions (2), therefore, without loss of generality, is considered a steady state of the system in presence of diffusion as well. System (1) is first linearised by introducing new variables namely and in such a way that they are perturbed slightly from the steady state , so that their relationships to and are . For and , the perturbed variables and are substituted in (1) and expanded using Taylor expansion for functions of two variables to obtain a linearised system which can be written in matrix form as
[TABLE]
In order to complete the linearisation in the presence of diffusion, it is necessary to find the eigenfunctions of the Laplace operator, that satisfies the homogeneous Neumann boundary conditions. Eigenfunctions of the Laplace operator on planar domains are well-studied in the literature [6, 8, 23, 33]. It is presumed that and both are of a similar form. The eigenfunctions for the Laplace operator are found by solving the relevant eigenvalue problems, that satisfy the given boundary conditions of problem (1). Such an eigenvalue problem is of the form
[TABLE]
The method of separation of variables is used to find the solution to problem (16) by presuming that the solution is in the form of a product of two functions namely and , so has the form . This form of solution is substituted in equation (16) to obtain two one-dimensional eigenvalue problems, which are individually solved and the resulting set of eigenfunctions solving (16) are given by \bar{u}_{n,m}(x,y)=C_{n,m}\cos{\big{(}\frac{n\pi x}{L}\big{)}}\cos{\big{(}\frac{m\pi y}{L}\big{)}}, where are the coefficients depending on the mode of the eigenfunctions. The solution to problem (15) can be, therefore, written as the sum of infinite expansion in the form of the product of , and the eigenfunctions for the two dimensional Laplace operator, so we have
[TABLE]
where and are the coefficients of the infinite expansion. Substituting this form of solution in (15), the problem can be written as a two-dimensional discrete set of algebraic equations in the form, with
[TABLE]
which can also be written as a two-dimensional discrete eigenvalue problem
[TABLE]
For the stability of the steady state in the presence of diffusion, the eigenvalues of problem (18) are analysed. The characteristic polynomial for (18) is given by
[TABLE]
which is a quadratic polynomial in . Let and denote the trace and determinant of the stability matrix (18) then the pair of eigenvalues of system (18) satisfy the quadratic polynomial in terms of the trace and determinant in the form
[TABLE]
where and are given by
[TABLE]
The roots of (20) are . We note that, in comparison to the case when no diffusion is present, the trace of the Jacobian matrix becomes more negative when diffusion is introduced since a positive quantity is being subtracted, thereby enhancing the stability of the system. This is equivalent to the default intuition that diffusion may only promote stability, due to its natural property of homogenising concentration gradient. However, the modification to determinant part in the the form of an additional positive term is intuitively expected to dominate the strictly positive term which is subtracted from , namely . These modifications to the trace and the determinant of the stability matrix due to the presence of diffusion to the system are investigated thoroughly to derive possible quantitative relationships between domain size, reaction and diffusion rates in the context of the parameter classification. We also attempt to quantitatively understand how this modification of both and contributes to creating instability. The regions of the parameter space for all types of steady states and their stability is analysed. Once the regions for different types of steady states are classified, then subject to this classification the influence of the domain size is explored on the types and existence of diffusion-driven instabilities.
3.2 Analysis for the case of complex eigenvalues
To determine the types of steady state , with corresponding parameter space, it is important to realise that despite the dependence on and , the type of steady state also depends on the additional parameters namely , and . The variation of the additional parameters are investigated for fixed and to explore how the parameter space varies for different values of , for example. Another important relationship to investigate is between the area of the domain and the reaction rate and how these influence the process of diffusion-driven instability. First, we analyse the case of complex eigenvalues. The eigenvalues are complex values if are chosen such that the discriminant of (20) is negative. It means that must satisfy
[TABLE]
One immediate consequence of this inequality is that must be positive, because if is negative, then for no choice of do the eigenvalues become complex. This immediate observation is automatically in agreement with the conditions for diffusion-driven instability presented in [6, 23, 20]. Furthermore, this condition on one hand will determine the region of the plane on which are a complex conjugate pair. On the other hand it serves to indicate the equations for the partitioning implicit curves on which the eigenvalues change type from real to complex. Bearing in mind, that inequality (21) is a sufficient condition for to have imaginary parts, therefore, it suffices to study (21) for to be a spiral. Before numerical treatment of the critical curve for the complex region, the stability of is studied analytically given that have imaginary parts. Stability of is determined by the sign of the real part of when it is a pair of complex conjugate values. Given that this is the case then the real part of is given by
[TABLE]
If we wish to find the parameter space for which is a stable spiral then we require the real part of both eigenvalues to be negative, which means in addition to the condition (21) we want to be negative or equivalently must satisfy
[TABLE]
In comparison to the absence of diffusion, the condition for stability of when are complex eigenvalues is similar to (22) except that the right-hand side of inequality is a strictly positive real quantity. The quantity on the right-hand side of (22) has in the denominator. The maximum value that in the absence of diffusion could attain on is 1. The expression for in the presence of diffusion is the same as in the absence of diffusion except it is reduced by a strictly positive constant and multiplied by . This relationship indicates that the stability of when are complex, depends on the extra parameters , and . The maximum of the left-hand side of (22) is the non-zero positive constant [see Section 4, Theorem 1], which suggests that the real part of is negative if the term on the right-hand side of (22) satisfies
[TABLE]
In the presence of diffusion, if is a pair of complex conjugate eigenvalues then the real part of is always negative given that the area of satisfies (23), otherwise no choice of , and can change the sign of to positive. This claim is formally proven in Theorem 1 as well as numerically demonstrated that no choice of , and can invalidate inequality (22), given that satisfies (23). For the first eigenvalue of the laplace operator, which corresponds to , inequality (21) takes exactly the same form as the inequality that was necessary for to be a spiral in the absence of diffusion, namely inequality (7). Therefore, to avoid repetition the stability analysis of the steady state with is skipped and the non-trivial case is studied, i.e. . In the absence of diffusion, in the region where is a spiral, there are sub-regions that correspond to a stable spiral as well as sub-regions that correspond to an unstable spiral. One immediate consequence of including diffusion, into the model, is that it puts the condition (23) on the domain size for the real part of to remain negative. Similarly the real part of becomes positive if condition (23) is violated. This means, diffusion-driven instability only serves to destabilise a spiral type of uniform steady state if the parameter violates inequality (23), which means
[TABLE]
Conditions (23) and (24) both have consequences in terms of the stability of the original system. The first consequence is that if have imaginary parts, then condition (23) guarantees that the real part of is negative for all choices of and , allowing the eigenvalues to be positive if and only if they are both real repeated or real distinct (with at least one of them positive). Hence, restricting the system from undergoing Hopf bifurcation. Similarly, if have imaginary parts, then condition (23) must be violated (condition (24) must be satisfied) in order for the real part of to become positive whilst the eigenvalues are still complex conjugate pair. It means that for the original system to undergo Hopf bifurcation, the necessary condition on the domain size is (24). Both of these statements are formally proven in the form of theorems in the Section 4. If the real part of the complex eigenvalues becomes zero, then the roots become purely imaginary, indicating that the system undergoes a time-periodic oscillations whilst experiencing spatial stability. This type of instability is referred to as Transcritical bifurcation and this entails that the real part of must become zero, which means
[TABLE]
Rearranging (25) and solving for , the necessary condition for transcritical bifurcation is obtained in the form
[TABLE]
If the parameter satisfies condition (24) simultaneously with condition (26) and , then the system exhibits a transcritical bifurcation.
3.3 Interpretation of the dynamics for the case of complex eigenvalues
For a domain size satisfying (23) the dynamics of the system with diffusion is expected to undergo spatial patterning only. If the diffusion rate is increased, yet satisfying (23), instability is still expected to invade stable regions in the parameter plane, however the only type of instability one may expect is spatial and not temporal. This type of instability is referred to as Turing Instability [22, 2, 31], where a certain region in the parameter space is stable in the absence of diffusion and yet becomes spatially unstable upon adding diffusion to the system. Findings of the current study suggest that, if the relationship between , and is maintained as given by (23), the type of diffusion-driven instability is always restricted to Turing Instability. It means that the only pattern one can achieve with (23) is spatially periodic pattern, with no temporal periodicity. However, if , and are given to satisfy condition (24) then the system may undergo both temporal and/or spatial periodicity. Upon increasing the value of , yet maintaining condition (24) the system is expected to become temporally unstable. This type of instability is referred to as Hopf Bifurcation, where adding diffusion to a temporally stable system, causes temporal instability. In this case one may expect the system to undergo temporal periodicity (pattern along the time-axis).
3.4 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 both eigenvalues are real values if the relationship between and is such that
[TABLE]
The equal case of (27) is looked at first, where we have
[TABLE]
which means that the discriminant is zero, hence the roots are repeated real values of the form , given by
[TABLE]
When and satisfy condition (28), the stability of the steady state is determined by the sign of the root itself. The expression given by (29) can be easily shown to be negative if the area of the domain satisfies the inequality
[TABLE]
Otherwise, the repeated root is positive provided that satisfies
[TABLE]
Analysing (30) and (31) 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 area of is not compared against a negative quantity, such a restriction is
[TABLE]
It must be noted that (32) is the same restriction on the parameter choice obtained for the case of repeated real eigenvalues in the absence of diffusion. The region where the eigenvalues are real repeated roots are implicit curves in the parameter space satisfying (27), these curves are computed numerically in the last part of Section 4. These curves form the boundary between the regions of complex and real eigenvalues. Varying the diffusion rate causes a shift to the location of the curves indicating clearly regions that are subject to diffusion-driven instability. The remaining case to look at is when both eigenvalues are real distinct. This happens if and are chosen such that the strict inequality case of (27) is satisfied. This case corresponds to the diffusion-driven instability Turing type only, because both eigenvalues are real and distinct.
3.5 Interpretation of the dynamics for the case of real eigenvalues
If both eigenvalues are negative distinct real values, then the system is spatially as well as temporally stable, the dynamics will achieve no patterns, hence the system returns to the uniform constant steady state as time grows, [see Section 5 Figure 11] with no effect from diffusion. If the eigenvalues are both real with different signs, then the type of instability caused by diffusion is spatially periodic or oscillatory in space, because this case corresponds to the steady state becoming a saddle point. If both eigenvalues are positive real distinct values, then the dynamics are expected to exhibit a spatially periodic pattern, in the form of stripes or spots.
4 Main results
This section contains the main results of our current study. The analytical findings are formally presented in the form of two theorems with proofs. Each of which indicates the relationship between the domain size with the diffusion coefficient and reaction rate . The numerical technique for computing the critical boundary curves is also briefly explained and results of parameter space classification are numerically computed. Numerical verification of both theorems is carried out in this section, where the reaction-diffusion system is numerically solved using the finite element method on a unit square domain. The relationship between domain length and parameters and is shown to be in agreement with theoretical predictions presented.
Theorem 1** (Turing type diffusion-driven instability).**
Let and satisfy the non-dimensional reaction-diffusion system with activator-depleted reaction kinetics (1)-(2) on a square domain with area and positive parameters , , and . Given that the area of the square domain satisfies the inequality (23) where then for all in the neighbourhood of the unique steady state (u_{s},v_{s})=\big{(}\alpha+\beta,\frac{\beta}{(\alpha+\beta)^{2}}\big{)} the diffusion driven instability is restricted to Turing type only, forbidding the existence of Hopf and transcritical bifurcation.
Proof.
The strategy of this proof is through detailed analysis of the real part of the eigenvalues of the linearised system, when the eigenvalues are a complex conjugate pair. This can be done through studying the surface , and finding that it has a unique extremum point at . The method of second derivative test and Hessian matrix is used to determine the type of this extremum. Upon finding its type, then the monotonicity of is analysed in the neighbourhood of the extremum point in both directions and . The monotonicity analysis and the type of the extremum leads to proving the claim of the theorem.
The eigenvalues in the presence of diffusion, in terms of trace and determinant are given by , where
[TABLE]
It can be immediately observed that in the neighbourhood of for the system to exhibit Hopf or transcritical bifurcation the discriminant of the characteristic polynomial must satisfy the inequality . Therefore, the stability and type of the steady state in this case is determined by the sign of the real part of . The aim is to investigate and derive from it condition (23) on as a requirement for to be negative for all strictly positive choices of , , and . First derivative test is used on to find all stationary points of on the domain . All stationary points of must satisfy , which is true if and only if
[TABLE]
Similarly all stationary points of must also satisfy , which implies
[TABLE]
The system of nonlinear algebraic equations obtained from (33) and (34) has a unique solution namely and [see Remark 4.0.1]. Therefore, has a unique stationary point at the origin. The type of this stationary point is determined by the second derivative test for which the Hessian matrix must be computed and evaluated at the point .
[TABLE]
It is clear that the entries of upon direct evaluation at the point are undefined. This is treated by using L’Hopital’s rule. L’Hopital’s rule sometimes does not work for functions of two variables defined on cartesian coordinates, therefore a transformation of the entries to polar coordinates might be applied. We will exploit this technique to express the Hessian matrix in polar coordinates and differentiate accordingly. The entries of are transformed to polar coordinates using and , so the rule can be applied by taking the . Using coordinates the entries of take the following form
[TABLE]
L’Hopital’s rule is applied to each entry of separately and the is computed for . Starting with the entry and cancelling , the expression takes the form
[TABLE]
Let and respectively denote the numerator and the denominator of the expression for , then the application of L’Hopital’s rule suggests that
[TABLE]
Applying the same procedure for , and , all the entries of are computed and given by
[TABLE]
Since the , therefore, the second derivative test suggests that is a saddle point of . Since it was previously shown that attains a unique stationary point in the domain , i.e. by solving the equations (33) and (34), therefore, if was a maximum and , this would suggest that, whenever has a non-zero imaginary part then regardless of the choice of , and , however due to fact that is a saddle point, it remains to show that is negative at and its first derivatives in the neighbourhood of of and for positive values of and in both directions are negative and do not change sign. Let and denote the curves for constants and respectively on the plane , then
[TABLE]
The expression for clearly satisfy that it is negative at and its first derivative in the direction of is for all . The expression for however is not trivially negative for all values, since the sign of the constant in the expression is positive, which if computed at , leads to the desired condition
[TABLE]
It has been shown that the condition (23) is necessary for to be negative at the unique stationary point namely , it remains to show that the first derivative of , for all which completes the proof. ∎
4.0.1 Remark
The parameters are not the admissible choices for the original system, because this choice of parameters leads to the trivial steady state , which is globally and unconditionally stable and has no physical relevance. However, the surface defined by , which determines the sign of the real part of attains a unique extremum at , which makes the analysis of the neighbourhood of this point important for studying the sign of . It is reasonable to argue that the point is not permissible to use as the platform of the proof of Theorem 1, in that case, proving the results become a step closer by only showing that is a strictly monotonically decreasing function (with no extrema in ) and it can only attain bounded positive values in the neighbourhood of if and only if the condition (23) on the domain size is violated. Therefore, given that maintains to satisfy (23), the sign of the real part i.e. of is guaranteed to be negative, which is a step shorter to reach the claim of Theorem 1. Therefore, it is brought to the attention of the reader that the use of the point as an extremum of is more of complementary factor to the proof rather than an essential one.
Theorem 2** (Hopf or transcritical bifurcation).**
Let and satisfy the non-dimensional reaction-diffusion system with activator-depleted reaction kinetics (1)-(2) on a square domain with length and positive parameters , , and . For the system to exhibit Hopf or transcritical bifurcation in the neighbourhood of the unique steady state (u_{s},v_{s})=\big{(}\alpha+\beta,\frac{\beta}{(\alpha+\beta)^{2}}\big{)}, the necessary condition on the area of the square domain is (24) where and are any integers.
Proof.
The strategy to this proof is some what different from that of Theorem 1, because in Theorem 1 the aim was to show that condition (23) must be satisfied by for to be negative and to maintain negative sign throughout the plane. For this proof the only important step is to derive the condition for to be positive, because when are complex eigenvalues then the sign of is precisely what determines the stability of . The system undergoes Hopf or transcritical bifurcation in the neighbourhood of , if the sign is positive. Therefore for diffusion-driven instability to influence when are complex conjugate pair it is necessary that
[TABLE]
which can equivalently be written as
[TABLE]
The expression on the left hand-side of (38) is explored in particular to find its upper bound, it can be written it in the form of difference between two rational functions as
[TABLE]
where and . The range for and are independently analysed to find the supremum of the expression on the left of (39). 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 (37) can be written as
[TABLE]
which by rearranging leads to the desired condition (24). ∎
4.1 Equations of the implicit partitioning curves
A similar method to that used in the case of the absence of diffusion is applied to solve for the partitioning curves in the parameter plane. The partitioning curve for the complex region, on the plane must satisfy
[TABLE]
since this is the curve on which the discriminant of the expression for eigenvalues change sign. It means the curve satisfying (40) determines the boundary on one side of which are both real and on the other side of it, are both complex conjugate eigenvalues. Equation (40) is solved by a similar method to that used for solving (9) on a square domain and for each fixed mesh point in the direction of it is found that its solution is equivalent to finding the positive real roots of the polynomial of degree 6 in for fixed values of in the form
[TABLE]
where for each the coefficients are given by
[TABLE]
[TABLE]
[TABLE]
[TABLE]
[TABLE]
[TABLE]
The implicit curve satisfying (41) indicates the choices of for which the eigenvalues are repeated real roots, since this is the curve on which the discriminant is zero. Therefore it forms the boundary region in the plane corresponding to complex eigenvalues. The polynomial of degree 6, given by (41) is solved for fixed on a domain with . For each the positive real roots of are extracted and plotted on the plane. The algorithm is run for five different values of the non-dimensional parameter , associated to both cases, where the area of satisfies condition (23) of Theorem 1 as well as condition (24) of Theorem 2.
The equation for the second curve that partitions the region corresponding to complex eigenvalues of the parameter plane is the curve on which the real parts of the complex roots are zero, hence satisfying the equation
[TABLE]
The solution of (42) by algebraic manipulation can be found to be equivalent to finding the positive real roots of the cubic polynomial for fixed , where is given by
[TABLE]
with C_{0}(\alpha_{i})=-\big{(}L^{2}\gamma+(d+1)(n^{2}+m^{2})\pi^{2}\big{)}-L^{2}\gamma\alpha_{i}^{3}, and .
Variations of the parameter are investigated, whilst the area remains to satisfy condition (23) and as expected it is found that in the region of the parameter plane corresponding to complex eigenvalues, there is no sub-region that corresponds to to have positive real parts. For each value of the region corresponding to complex eigenvalues is tested by looking for a critical curve on which is purely imaginary, i.e. satisfying (43). If such a curve exists, it would correspond to the system undergoing periodic oscillations around , thus the system exhibiting transcritical bifurcation, which also implies the existence of a region in the parameter space that corresponds to real parts of to be positive. Upon investigating this region, it is found that such a curve under condition (23) does not exists, and all roots corresponding to the cubic polynomial given by (43) are either complex eigenvalues or they are negative real values, therefore cannot be the choice of admissible parameters of the system. Another observation is that as increases the area of the region corresponding to complex eigenvalues in the plane gradually decreases. The domain size of length satisfying the condition of Theorem 1 is tested for stability analysis of when are complex eigenvalues. The results are shown for different values of by a distinct colour in Figure 4 (a). Finding the solutions of these curves does not however indicate, which side of them correspond to complex and which side to real . This is decided by some numerical trial and error by evaluating using few values from both sides of each curve. Trial and error indicates that the regions under these curves correspond to to be a complex conjugate pair and hence any combination of from this region ensures that the eigenvalues contain a non-zero imaginary part. When the area of satisfies condition (23), the absence of a region satisfying (43) verifies the statement of Theorem 1, therefore the region corresponding to complex eigenvalues for the choice of satisfying (23) has no sub-partitions, because everywhere in this region the real part of the eigenvalues is negative, hence no choice of parameters could result in the system to exhibit Hopf or transcritical bifurcation. The eigenvalues only become positive when they are both real values, therefore condition (23) restricts the diffusion-driven instability to Turing type only. Figure 6 (a) shows how the partitioning curve changes location, as the value of is varied. The region under those curves corresponds to the complex eigenvalues. Similarly Figure 6 (b) indicates the regions corresponding to complex eigenvalues for the corresponding values of . Each stripe in Figure 6 (b) is denoted by a letter that represents the set of all points corresponding to a distinct colour stripe. Similarly, the regions where the eigenvalues are negative real roots are presented by Figure 7 (a), corresponding to the same values of the parameter . Figure 7 (b) shows regions where at least one or both eigenvalues are positive real roots. The summary of Figures 6 and 7 is presented in Table 2.
The algorithm is also run for the case when is chosen such that it satisfies (24) and we find that the region corresponding to complex eigenvalues is further partitioned by the curve satisfying (42). This is the curve on which the eigenvalues are purely imaginary. This curve also indicates that within the region corresponding to complex eigenvalues there is a sub-region in which the eigenvalues are complex with positive real part, which corresponds to the system exhibiting Hopf bifurcation. For choices of parameter values on the curve the system is expected to undergo transcritical bifurcation. The area of the domain is taken as in order to satisfy the condition given by Theorem 2 with respect to and , so that there is enough space for varying and yet maintaining condition (24). It can be easily observed from the parameter space classification that there is a relatively small region in each case corresponding to diffusion-driven instability. This indicates the importance of making sure to choose parameter choices wisely, in order to expect the dynamics of the system to evolve to certain types of patterns or maybe no patterns at all. Table 3 presents the summary for how the regions of the parameter space change with varying the parameter . It would be reasonable to use the same variation of the parameter in both cases of the domain sizes, however, in the first case where the domain size satisfies (23) the span of varying parameter is relatively small, yet causing significant observable change in the parameter space. However, when the domain size is chosen according to condition (24), small variations in the diffusion coefficient makes insignificant changes to the parameter spaces, therefore, in order to pictorially observe the consequential change in the parameter plane the span of variations for have to be significantly large as seen in Table 2.
In the case when satisfies Theorem 1, it can be observed that the type of diffusion-driven instability that can occur is restricted to Turing type only, which is increased as the non-dimensional diffusion coefficient was increased. Turing diffusion-driven instability also occurs for the case when satisfies condition (24) and it also increases with increased value of as shown in Figure 9 (b). The interpretation is that regions in the parameter space exist, which result in the system to be stable in the absence of diffusion, but when diffusion is added to the system, the choice of the parameters from these particular regions result in the system exhibiting instability. This type of instability is restricted to space and hence leads to spatial patterning only, because the eigenvalues are both real. However, if satisfies (24), then in addition to the existence of regions of the parameter space corresponding to Hopf and transcritical bifurcations, it can also be observed that with increased values of parameter , unlike Turing instability, the regions for Hopf and transcritical bifurcations reduce as shown in Figure 10 (a) and (b). The mathematical intuition for this inverse relation is the direct proportionality of on the lower bound for the domain size i.e. the right-hand side of (24). Therefore, as grows, one gets closer to the violation of the necessary condition (24) for Hopf and transcritical bifurcations as proposed by Theorem 2. However, the region for Hopf and transcritical bifurcations remains to exist as long as satisfies (24), which ultimately means that for the system to exhibit temporal periodicity (patterning in time) the domain size has to be sufficiently large satisfying (24).
4.1.1 Remark
The eigenfunction modes used for the current simulations corresponds to the first non-trivial (non-zero) eigenvalues, namely the case where . The results can be readily obtained for any positive integer values of and . Despite the fact that using larger integer values for and will shift or scale the classification for the parameter space, the conditions given by Theorem 1 and Theorem 2 will remain intact irrespective of the values of and , so long as they are positive integers.
4.2 Numerical experiments
To validate the proposed classification of parameter space, the reaction-diffusion system (1)-(3) is numerically solved using the finite element method [1, 5, 9, 34] on a unit square domain with a uniform triangular mesh. Numerical simulations are performed for various choices of parameter values and , chosen from the appropriate parameter spaces to demonstrate and validate our theoretical findings. In all our simulations, we vary the parameters and and keep fixed the domain length size and this allows us to keep constant the well refined number of degrees of freedom for the mesh. The initial conditions for each simulation are taken as small random perturbations around the neighbourhood of the steady state of the form [19, 35]
[TABLE]
The final time is chosen to ensure that beyond the convergence of solutions in successive time-step differences decay to a threshold of for all the simulations. Therefore, in the case where the area of the domain size satisfies condition (23) proven in Theorem 1, when temporal instability is forbidden by the condition on the domain size, the final time of numerical simulations is relatively shorter compared to that used for other cases. For numerical simulations on the domain size satisfying condition (24) proven in Theorem 2, the final time is experimented for longer periods to capture possible existence of temporal periodicity in the dynamics. Simulations with domain size satisfying (24) are captured at different times, which varies with cases, for which details are included in Table 4. In all our numerical results, we only exhibit numerical solutions corresponding to the component, those of are known to be 180-degrees out of phase to those of .
Figure 11 presents the case where the parameters and are chosen from stable node region of Table 2, with final time . It can be observed that the evolved profile of the concentration uniformly converges to the steady state value forming neither spatial nor temporal patterns (Figure 11 (a)). As predicted in Theorem 1, no choice of parameters and can influence the dynamics to exhibit temporal periodicity, therefore any choice of the parameters outside Turing space, given that the domain size satisfies (23) will uniformly converge to the stable steady state as seen in Figure 11 (a). Figure 11 (b) shows the uniform convergence of the discrete -norm difference between solutions at successive time-steps to the constant steady state , where there is no sign of instability occurring, when parameters are outside the Turing space.
Figure 12 (a) presents the evolved profile of the solution captured at the final time step, for the choice of parameters from Turing space presented in Table 4. The convergence in the discrete -norm of difference in the solutions for the successive time steps from start to the end of the simulations are also plotted against time and presented in Figure 12 (b).
Figure 13 shows the simulation for the choice of parameters from the Turing space presented in Table 3 for the domain size satisfying (24). In this case as predicted in Theorem 2, regions of parameter space exist for which the dynamic can exhibit Hopf type bifurcation, therefore, the possibility of temporal periodicity in the dynamics. Figure 14 presents such periodicity in time for spatial patterns. The relative discrete -norm of the difference in the solutions for successive time-steps is therefore showing time periodicity as illustrated in Figure 14 (d), which indicates the transition of the solution from the initially achieved spatial pattern to a different spatial pattern.
5 Conclusion
In this work the full parameter space for a reaction-diffusion system with activator-depleted reaction kinetics was classified by use of linear stability theory and in each region of the parameter space the dynamics of the reaction-diffusion system was investigated. In the absence of diffusion, theoretical results on the dynamics of the system were supported by use of the phase-plane analysis, where in each case the numerical solution of the system was observed to be in agreement with the theoretically predicted behaviour. In the presence of diffusion, for a two-component reaction-diffusion system, two conditions relating the domain size to the diffusion and reaction rates were derived. The proofs of these conditions were presented in Theorems 1 and 2, respectively. For full classification of the parameter space, a numerical method was used to compute the solutions of the implicit curves in the parameter space forming the partitions of classification. In particular, using condition (23) the numerical method for solving the partitioning curves showed the non-existence of a region in parameter space that (if existed) would lead to Hopf bifurcation or transcritical instability. Similarly, applying the numerical method to compute the partitioning curves under condition (24), it was shown that regions in parameter space exist for both Hopf and transcritical bifurcation. Parameters from Hopf bifurcation region as well as from Turing regions under both conditions on the domain size were shown to be in agreement with the theoretical prediction, when the system was numerically solved by using the finite element method. For each simulation the discrete -norm of the successive time-step difference of the solutions is also given to visualise the temporal dynamics of the behaviour of the solutions during the convergence process.
5.1 Ideas for future work
This work sets the premises to study more complex systems of non-autonomous reaction-diffusion equations during growth development whereby the parameter spaces are continuously evolving due to domain growth. It will be revealing to study how the Turing diffusion-driven parameter space, the Hopf and Transcritical regions evolve with time and how the dynamics of the model system evolve. The motivation for such extension is to find whether the conditions (23) and (24) for the domain size with reaction and diffusion rates continue to hold or whether a threshold for the domain size exists, beyond which, these conditions can be invalidated. Another direction of further investigation of the current work is to try and find similar relationships for different geometries such as spherical, elliptical and cylindrical domains. It may be noted that the eigenvalues and eigenfunctions of the Laplace operator change significantly depending on the boundary conditions and the geometry of the domain, therefore, application of the present method suggests, that finding similar conditions for other geometries require a careful step-by-step procedure to find analogous conditions to those given by (23) and (24). It is also possible to apply the idea of the present work to the problems of pattern formation on bulk-surface geometries, where the relationship of surface area with reaction and diffusion rates can be explored to find the relevant influence of this relationship on the formation of spatial and/or temporal patterns. In [20] it is shown that the known condition of is no longer necessary for pattern formation in the presence of linear cross-diffusion, which also provides a possible platform for the extension of the current work, to see, whether the invalidation of the condition in the presence of linear cross-diffusion has any dependence with the relationship of reaction and diffusion rates with the domain size.
6 Acknowledgments
WS acknowledges support of the EPSRC Doctoral Training studentship. AM acknowledges support from the Leverhulme Trust Research Project Grant (RPG-2014-149), the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 642866, and his work was partially supported by the Engineering and Physical Sciences Research Council, UK grant (EP/J016780/1). The authors (WS, AM) thank the Isaac Newton Institute for Mathematical Sciences for its hospitality during the programme (Coupling Geometric PDEs with Physics for Cell Morphology, Motility and Pattern Formation; EPSRC EP/K032208/1). AM was partially supported by Fellowships from the Simons Foundation. AM is a Royal Society Wolfson Research Award Holder, generously supported by the Wolfson Trust.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] A. Madzvamuse, A. H. W. Chung , The bulk-surface finite element method for reaction-diffusion systems on stationary volumes ., J. FE in Analysis and Design. 108 (2015) 0–21.
- 2[2] A. Madzvamuse, A. H. W. Chung , Fully implicit time-stepping schemes and non-linear solvers for systems of reaction-diffusion equations ., J. App. Math. Comp. 244 (2014) 361–374.
- 3[3] S. Yan., X. Lian, W. Wang, R. K. Upadhyay , Spatiotemporal dynamics in a delayed diffusive predator model ., J. App. Math. Comp. 244 (2013) 524–534.
- 4[4] R. Barreira, C. M. Elliot, A. Madzvamuse , The surface finite element method for pattern formation on evolving biological surfaces ., Online. J. Math. Biol. 29.
- 5[5] R. Zhang. X. Yu, J. Zhu, A. F. D. Loula , Direct discontinuous Galerkin method for nonlinear reaction-diffusion systems in pattern formation ., J. App. Math. Mod. 38 (2014) 1612–1621.
- 6[6] A. Madzvamuse , Stability analysis of reaction-diffusion systems with constant coefficients on growing domains ., Int. J. Dyn. Diff. Eq, (2008).
- 7[7] C. Venkataraman, O. Lakkis, A. Madzvamuse , Global existence for semilinear reaction-diffusion systems on evolving domains ., J. Math. Biol. 64 (2012) 41–67.
- 8[8] J. A. Mackenzie, A. Madzvamuse , Analysis of stability and convergence of finite difference method for a reaction-diffusion problem on a one dimensional growing domain ., J. Num. Anal. 322 (10) (2009) 891–921.