Exploring Mechanisms for Pattern Formation through Coupled Bulk-Surface PDEs in Case of Non-linear Reactions

This work explores mechanisms for pattern formation through coupled bulk-surface partial differential equations of reaction-diffusion type. Reaction-diffusion systems posed both in the bulk and on the surface on stationary volumes are coupled through linear Robin-type boundary conditions. The presented work in this paper studies the case of non-linear reactions in the bulk and surface, respectively. For the investigated system is nondimensionalised and rigorous linear stability analysis is carried out to determine the necessary and sufficient conditions for pattern formation. Appropriate parameter spaces are generated from which model parameters are selected. To exhibit pattern formation, a coupled bulk-surface finite element method is developed and implemented. The numerical algorithm is implemented using an open source software package known as deal.II and show computational results on spherical and cuboid domains. Also, theoretical predictions of the linear stability analysis are verified and supported by numerical simulations. The results show that non-linear reactions in the bulk and surface generate patterns everywhere. Keywords—Bulk-surface; Reaction-diffusion; Finite-ElementMethod (FEM); Partial Differential Equations (PDEs)


I. INTRODUCTION
Most biological and chemical processes that can be explored through reaction and diffusion of chemical species are often modelled by systems of partial differential equations (PDEs) [1]- [3]. A special class of these are reaction-diffusion equations, which are used to analyse and quantify various biological processes such as the natural evolution of pattern formation on animal coats, developmental embryology, immunology, ecological dynamics [4]- [7]. The study of reactiondiffusion systems in general has been and continues to be an interesting topic for research in various branches of scientific studies. In order to quantify the evolution of chemical reaction kinetics associated to biological processes, it is a usual approach to employ a system of partial differential equations describing the chemical reactions, which is investigated through mathematical techniques to reveal the long-term behaviour of the evolving kinetics [8], [9].
Alan Turing was one of the first scientists to suggest in 1952 the use of a system of reaction-diffusion equations to model how two or more chemical substances evolve when they are simultaneously subject to a specific reaction rate and each one of them diffuses independently of the other. Alan Turing suggested that the theory of biological pattern formation can be mathematically formulated by a system of partial differential equations [10]. The study in [11] is a theoretical set-up through coupling reaction-diffusion system to provide insight on trajectory of a particle during the process of bulk excursion, when it unbinds from the surface without a regular occurrence.
Turing's theory suggests that pattern formation occurs, when a system experiences diffusion-driven instability [10], [12], which is a concept that is hypothetically responsible for the emergence of spatial variation in the concentration density of a chemical species. Diffusion-driven instability takes place in the evolution of a system, when a uniform stable steady state is destabilised by including the effects of the diffusion process in the system. It is a non-trivial property of the diffusion operator that it can be responsible to destabilise a stable steady state of a system of partial differential equations, because a diffusion operator by itself has the property to homogenise small spatial perturbations, therefore, intuitively if diffusion is added to a system of reaction kinetics that is stable in the absence of diffusion, then small perturbations near a uniform steady state are expected to ensure that the evolution of the reaction-kinetics converges to the uniform steady state.
Researchers in applied mathematics and computational science also explore bulk-surface reaction-diffusion systems (BSRDSs), which are employed in special kinds of models for biological processes, where species react and diffuse in the bulk of a domain and these are coupled with other species that react and diffuse on the surface of the domain. Bulksurface reaction-diffusion systems are employed as a framework to model the chemical interaction of bulk-surface problems arising in cell biology [13]. In particular the framework proposed by [13] aims to provide improved computational and algorithmic efficiency, which is mainly achieved, through employing the usual diffusion on local tangential planes as an approximation of Laplace-Beltrami operator. The framework proposed by [13] is applied to a realistic cell-like geometry, which produces results that are in agreement with quantitative experimental analysis on fluorescence-loss in photo-bleaching. Another example of a computational approach to solving coupled systems of BSRDEs is the work presented in [14], where they proposed a computational approach to bulk-surface reaction-diffusion systems on time-dependent domains.
In general there are two main aspects to the study of bulk-surface reaction-diffusion equations. The first approach is to solve systems of bulk-surface numerically. Finite element method is the usual choice of the numerical method in the literature, for example there is a detailed study in [15], suggesting some results on the numerical analysis, existence and convergence of finite element approximation when bulksurface reaction-diffusion equations (BSRDEs) are posed with Robin-type boundary conditions. A priori error bounds on the finite element approximate numerical solution are also both derived in certain norms and verified numerically. The work in [15] is concentrated mainly on the numerical analysis side of the particular scheme they present, which lacks to provide any insight on the stability analysis of the proposed system. Though, it is a reasonable decision to exclude stability analysis due to consistency and relevance of contents, however with improvements in computational efficiency of BSRDEs, it is crucial that attention is given to stability analysis of such systems.
Bulk-surface systems with a single PDE posed in the bulk and coupled with another PDE on the surface also play a vital role in the understanding the interaction of receptorligand in the process of a signalling cascade [16]. In [16] the existence of solutions is proven with some computational results associated to the theoretical problem, again lacking to provide insight on the stability behaviour of the dynamics modelled by the coupled system. Even though the results achieved in [16] are mathematically sound from a numerical analysis and computational viewpoint, it would provide a complementary back-up to the work if it is equipped with detailed results of stability analysis.
Stability and bifurcation analysis are two other usual analytical approaches to understanding the dynamical properties of reaction-diffusion system near a uniform steady state [17]- [21]. It is evident from the literature on the subject of stability analysis that a very limited amount of work is done on stability analysis in a coupled bulk-surface set-up. This is mainly due to the extensive complexity associated in deriving the relevant conditions for diffusion-driven instability when equations from the bulk are coupled with equations on the surface. One of the first detailed studies on stability analysis of BSRDEs is conducted in [20], where it is analytically proven that a certain suitable parameter range exists for equations in the bulk that can induce spatial pattern on the surface.
Coupled systems of bulk-surface reaction-diffusion equations (BSRDEs) are one of the several generalisations of reaction-diffusion theory to explore numerous applications in mathematical biology. Processes that involve bulk-surface reaction and/or diffusion are found in various research disciplines such as experimental research in organic chemistry, where a bulk-surface photografting process is used as an efficient tool to create thick grafted layers of hydrophobic polymers in a very short span of time [22], [23].
Bulk-surface reaction kinetics are also used to investigate the behaviour of chemical reactions in the interior of a cell, and to explore how a set of specific reaction kinetics in the interior of a cell evolve to influence the surface of the cell [24]. We also find bulk-surface reaction-diffusion equations that model a particular aspect of cellular functions with relevance to chemical signalling. In [25] a detailed mathematical model is developed for this particular investigation, to explore the dynamics of pattern formation in the consequences of bulk-surface coupling reaction kinetics. Moreover, bulk-surface reaction-diffusion equations help to reveal the mechanism of symmetry breaking which is one of the essential steps before the emergence of polarisation of biological cells or buds in yeast cells, the direction of cell motility [25].
Bulk-surface reaction-diffusion system are also used to model how surface active agents (surfactants) evolve on the surface of a system, in which the chemical concentration is coupled through a given reaction with the substance in the bulk [26]. BSRDSs also arise in mathematical models for the dynamics of lipid raft formation on biological membranes [27], where the formation of the layer on a biological membrane is modelled as the consequence of coupling conditions with species that react and diffuse in the bulk. A further example of biological application employing bulk-surface reactiondiffusion systems is presented in [28], where they model the mediation of cellular metabolism and signalling in part by trans-membrane receptors that undergo the process of diffusion in cell membrane. From the variety of applications that employ BSRDSs, one realises that a robust study of such systems can provide solutions to a great number of important questions in mathematical biology. This in turn requires in-depth and rigorous study of BSRDSs in an attempt to achieve extensive insight on the evolving properties of these models. Most of the published work presented in the current section on the study of BSRDSs either investigate an over-simplified case scenario with the aim of mathematical tractability or a complex model with limitations on the robustness of analytical and numerical findings.
In this paper, the presented study is motivated to explore BSRDSs with a realistic degree of complexity through a fourcomponent reaction-diffusion system, two of which are posed on the surface and the other two are posed the bulk. The equations in the bulk and on the surface also satisfy coupling conditions through the evolution dynamics on the surface is influenced by the reaction-diffusion process inside the bulk. It can prove of great importance to obtain insight on the pattern formation properties of such systems. The tools to achieve this in the current thesis are the combined application of linear stability theory, mode isolation and the finite element method.
The remaining part of the paper is organised as follows. Section 2 a study is conducted through the application of rigorous linear stability theory which is applied to analytically explore and predict the pattern formation properties associated to the adopted bulk-surface reaction-diffusion system. This is done through investigating the necessary conditions for diffusion-driven instability for the system. Section 3 presents deriving a set of sufficient conditions for diffusion-driven instability, which complements the necessary conditions of the previous section in order to insure that spatial pattern is obtained. In Section 4 the theoretical formulation for the finite element method is presented for the investigated system. Section 5 contains the numerical simulations obtained using Deal.II library to verify the analytical predictions associated to the pattern formation properties for the three systems. Finally, Section 6 concludes the presented work in this paper.

II. ANALYSIS OF COUPLED SYSTEM OF BULK-SURFACE REACTION-DIFFUSION EQUATIONS
In this section, we formulate and present the coupled systems of bulk-surface reaction-diffusion equations on stationary volumes, in which two of the equations are posed in the bulk and coupled with two other equations that are posed on the surface bounding the corresponding stationary volume. Reaction-diffusion systems posed both in the bulk and on the surface are coupled through linear Robin-type boundary conditions. For the investigated system we analyse non-linear reaction kinetics both in the bulk and on the surface. The details of the scaling process that makes the system studied in this paper dimensionless is presented. Also, linear stability analysis is carried out both in the absence and presence of diffusion, the necessary and sufficient conditions for steady state to be stable are derived in the absence of diffusion. In the presence of diffusion, the necessary conditions for diffusiondriven instability are derived. The theoretical results for this system show that the bulk dynamics and the surface dynamics drive pattern formation.
Let Ω ⊂ R 3 be a stationary domain with boundary that is a compact hyper-surface denoted by Γ ⊂ R 2 . Let u : Ω × (0, T ] → R and v : Ω × (0, T ] → R denote the concentration of two chemical species which react and diffuse in Ω. Let r : Γ × (0, T ] → R and s : Γ × (0, T ] → R denote two chemical species residing on the surface.
When the species from the bulk and surface are coupled only through the reaction kinetics and there is no crossdiffusion, it means that all four species diffuse independently of each other, which can be written in dimensional form with independent diffusion rates as follows: with coupling boundary conditions where, Ω is a three dimensional fixed domain bounded by a compact surface denoted by Γ, which means that it is a boundary-free connected and closed surface. The strictly positive constants D u > 0, D v > 0, D r > 0 and D s > 0 the independent diffusion rates corresponding to the variables indicated in the respective subscripts of each D.
We assume f (., .) and g(., .) to be non-linear functions. The coupling conditions of the system is represented by h 1 and h 2 which are functions of u, v, r and s. h 1 and h 2 denote reactions of substances through boundary interface, therefore they depend on all four species namely u, v, r and s.
We explicitly define h 1 (u, v, r, s) and h 2 (u, v, r, s) to be The constants α 1 , α 2 , β 1 , β 2 , κ 1 and κ 2 are positive parameters of system (1). We also assume that from all the species we initially have some positive quantity present, which we denote by u 0 , v 0 , r 0 and s 0 , which provides the initial conditions for system (1) written as , and s(x, 0) = s 0 (x).
In this system, we focus on the widely known activatordepleted model also known as the Brusselator model where the reaction kinetics are non-linear, given by with positive parameters k 1 , k 2 , k 3 and k 4 .

A. Non-Dimensionalisation
The system of equations is non-dimensionalised using a specific scale, in space or time, for observing the prospective solution within the specified scale range. In the new system after non-dimensionalisation, the variables and parameters are all unitless and the parameters will be fewer than in system (1). The non-dimensional variables is introduced with a hat and these are written asû,v,r andŝ with the corresponding scaling factors u * , v * , r * and s * , respectively. The process of nondimensionalisation is only represented for the bulk-equations in three spatial dimensions, and the process is identical to nondimensionalise the surface equations where a two-dimensional surface is embedded in three dimensional space. We choose L to denote the scaling factor for length (L b for the bulk and L s for the surface) and t * to denote the scaling factor for time (t * b for the bulk and t * s for the surface), The dimensional and the non-dimensional variables are related through where for the bulk we use the scaling given by τ and for the surface equations we use We substitute for each dimensional variable its corresponding product of non-dimensional variable and the scaling factor leading to whereΩ andΓ, respectively denote unit cube and its six sided surface. The scalingT denotes the final time for the nondimensional system.
Multiplying (6), (7), (8) and (9) by s * , respectively provided that u * , v * , r * and s * are non-zero. We may choose to define t Du and t * s = L 2 s Dr and factoring out some parameters will result in writing the system as: where d Ω = Dv Du and d Γ = Ds Dr express the non-dimensional positive ratios of diffusion parameters. Requiring the terms k3 k2 u * 2 = 1 and k3 k2 r * 2 = 1 to be non-dimensional respectively imply defining u * = k2 k3 and r * = k2 k3 . The scaling factors v * and s * through a similar process may be derived as Substituting (12) in system (10) results in where the new dimensionless parameters γ Ω = k2 , δ 2 = κ1 k2 and δ 3 = κ2 k2 are defined as a consequence of the scaling choice used for u * , v * , r * and s * .
The boundary and initial conditions are nondimensionalised through the same choice of scaling factors for all variables. For notational convenience we drop all the hats from the non-dimensional variables to obtain the full system of BSRDEs given by (1) in its non-dimensional form as in Ω × (0, T ] ∂r ∂t with linear boundary conditions The non-dimensional initial conditions for all equations are given by The parameter γ Ω is known as the reaction scaling parameter in the bulk and γ Γ is the reaction scaling parameter on the surface and both are non-dimensional. B. Linear Stability Analysis in the Absence of Diffusion Definition 2.1: (Uniform steady state): [10], [29] A point (u 0 , v 0 , r 0 , s 0 ) is a uniform steady state of the coupled system of bulk-surface reaction-diffusion equations (13) if it solves the nonlinear algebraic system given by f i (u 0 , v 0 , r 0 , s 0 ) = 0 , for all i = 1, 2, 3, 4 and satisfies the boundary conditions given by (14).
We derive the uniform steady state by solving the algebraic system (18) such that the boundary conditions given by (14) are also satisfied We add (15) and (16) to obtain Upon substituting u 0 into (16), we find Through a similar straightforward algebraic manipulations we also find the steady state expressions for r 0 and s 0 in the form Therefore, the uniform steady state solution satisfying system (13) is of the form (23) Substituting the uniform steady state (23) in (17) and (18) parameters is derived by direct substitution of 23 and algebraic manipulations through the following steps result in Combining (24) and (25) we obtain the required condition on the parameters in the form Therefore, in order for (23) to be a steady state of system (13), a condition on the parameters is required to hold, which is These findings are summarised in the following theorem.
Theorem 2.1: (Existence and uniqueness of the uniform steady state) [20] The coupled system of BSRDEs (13) with conditions (14) admits a unique non-zero steady state given by provided the following compatibility condition on the coefficients of the coupling terms is satisfied Finally, we set out the summary of the necessary and sufficient conditions for Re(λ) < 0 in Theorem 2.2.
Theorem 2.2: (Necessary and sufficient conditions for Re(λ) < 0 ) [10], [29] The necessary and sufficient conditions such that the zeros of the polynomial p 4 (λ) have Re(λ) < 0 are given by the following conditions:

C. Linear Stability Analysis in the Presence of Diffusion
We start by analysing the system by taking the diffusion terms into account and performing the linear stability analysis. We introduce a small perturbation in the neighbourhood of the steady state, namely, (u 0 , v 0 , r 0 , s 0 ). We introduce the small perturbations up to the linear term in the form of If we substitute these small perturbations into the system we obtain and also Similarly we substitute such perturbations in the reaction terms. Since we know that at the steady state f (u 0 , v 0 , r 0 , s 0 ), g(u 0 , v 0 , r 0 , s 0 ), h 1 (u 0 , v 0 , r 0 , s 0 ) and h 2 (u 0 , v 0 , r 0 , s 0 ) all equal zero, therefore we aim to collect terms in such a way to determine the relative expressions for the steady state in each equation. Furthermore, we aim to perform linear stability analysis. Performing the algebra and cancelling the expressions for steady state and ignoring higher order terms will transform the equations into linearised system of equations.
For the remaining of this work, the analysis is restricted to circular and spherical domains, where the Cartesian coordinates are transformed to polar coordinates. The coordinate transformation is done mainly for the convenience of applying the separation variables. A close form solution can be written in the form: which are substituted in the linearised system of equations, to obtain ψ k l,m (x)u l,m (t) = ∆ψ k l,m (x)u l,m (t), For equations on the surface the relations may be written as whereas for the bulk the relations take the form We consider a coordinate transformation in which a vector x may define every point in the bulk by the variables r (radial distance from the origin) and y (a point on the surface), with the relationship x = ry where r ∈ (0, 1), y ∈ Γ.
Since the eigenvalue of the problem on the surface depends on l itself, therefore we may consider positive integers only, and m can be any integer with the restriction |m| ≤ l. This is because the eigenvalues of both problems are equal at r = 1. Note that if r = 1 for the eigenvalue problem in the bulk, then the eigenvalues associated to the usual diffusion operator must coincide with those associated to Laplace-Beltrami operator on the surface, which means the following relation must hold.

−k 2
l,m = −l(l + 1) Now we summarise the results in the following theorem: Theorem 2.3: [10], [29] The necessary conditions for diffusion-driven instability for the coupled system of BSRDEs (13) and (14) are given by and and/or

III. MODE ISOLATION AND PARAMETER SPACE GENERATION
In this section we proceed with the process of deriving sufficient conditions for diffusion-driven instability that shall complement the necessary conditions found in section II to ensure the emergence of spatial patterns. As the standard requirement of this process, we start by extracting excitable wavenumber through the analysis of critical diffusion ratio. Eigenvalues and eigenfunctions of the laplace operator are briefly discussed on the surface. The results for mode isolation for the excitable wavenumber are employed to computationally find Turing parameter spaces on the real positive parameter plane. We also present the process of coordinate transformation from Cartesian to spherical of the usual laplace operator. Finally, we analyse and compare the shift and dependence of Turing spaces for equations in the bulk with those Turing spaces that are derived for equations on the surface.

A. Critical Diffusion Ratio and Excitable Wavenumber
For the bulk, the conditions (34), (35) and (38) are necessary but not sufficient for the emergence of an inhomogeneous spatial structure. The sufficient condition requires the existence of some finite wavenumbers k 2 ∈ (k 2 − , k 2 + ), where k 2 ± are the roots of the equation H 2 (k 2 ) = 0. Also, for the surface the conditions (36), (37) and (39) are necessary but not sufficient for diffusion-driven instability and the sufficient condition requires the existence of some finite wavenumbers l(l + 1) ∈ (l(l + 1) − , l(l + 1) + ) where l(l + 1) ± are the roots of the equation H 1 (l(l + 1)) = 0.
When the minimum H 2 (k 2 ) = 0, we require that For fixed parameters on the kinetics in the bulk, the critical diffusion d c is obtained using the following form: Corresponding to the critical diffusion coefficient d c , there exists a critical wavenumber k 2 c that obtained by: This is the critical wavenumber, the sufficient condition for Turing instability with the necessary conditions (34), (35) and (38) satisfied, which leads the system to evolve into spatial pattern.
Similarly, the critical diffusion coefficient d c on the surface can be obtained from the following equation: The critical wavenumber on the surface is given by which provides the sufficient condition for diffusion-driven instability on the surface.
Since the diffusion coefficient must be greater than 1, then we only take the critical diffusion coefficient ratio as d c = 8.56762745781.  To verify that d < d c does not allow Turing pattern to evolve, the necessary conditions are tested on the parameter space (a 2 , b 2 ) where a 2 and b 2 are the positive constants of the Schnakenberg reaction kinetics. It is found that when d < d c , there is no region in the parameter space that would become unstable due to diffusion in the system. This is shown in Fig.  2(a). Similarly when d > d c , then we see that the unstable region is formed in the parameter space (yellow region in Fig.  2(b), which corresponds to the parameter values that would result in the system to evolve into a Turing pattern.

B. Mode Isolation in the Bulk
With the help of linear stability analysis, certain modes can be isolated to help find the admissible set of parameter values d Ω and γ Ω for diffusion-driven instability. The necessary conditions for diffusion-driven instability are One of the sufficient conditions however, for diffusion-driven instability is that the eigenvalues of the laplace operator should fall in the real interval between the small and the large eigenvalues of the system. It means that: must hold with L and R expressed by and respectively. Therefore, for sufficient condition to exists for diffusion-driven instability, the excitable modes must exist and belong to the interval (51).
A similar approach is applied to the case in two dimensions. The values of d Ω and γ Ω are computed. We are interested in finding combination of d Ω and γ Ω , such that the curve Re(λ(k 2 )) encapsulates only one excitable wavenumber. The algorithm is outlined through the following steps; • Define d Ω = d c + where 0 < 1 and d c = 8.5676.
• If k 2 l,m > k 2 + as shown in Figure 4 then increase the value of γ Ω by 1, till the curve includes the wavenumber by shifting to the right.
• If k 2 l,m < k 2 − then decrease the value of γ Ω by 1, till the curve includes the wavenumber by shiftting to the left.
• If there exist two excitable wavenumbers as shown in Figure 5(a) then we decrease till we obtain a unique excitable wavenumber as shown in Fig. 5(b).  A similar procedure may be employed to extract excitable wavenumbers from the spectrum of Laplace-Beltrami operator as a sufficient condition for Turing pattern to be formed on the surface.

C. Turing (Parameters) Space on the Surface
The Turing (parameters) spaces for equations posed on the surface and the conditions for these are obtained and outlined as The parameter spaces are derived on the actual positive real parameter plane (a 2 , b 2 ), for two choices of diffusion ratios, namely, d Γ = 20 and d Γ = 30 as shown in Fig. 6. Unstable region is shown in the parameter space (yellow region).

D. Turing Spaces in the Bulk and on the Surface
The following sub-figures 7, 8 show diffusion-driven instability spaces for the conditions on diffusion-driven instability given by (34)-(39) in the bulk and on the surface. We combine the Turing spaces (more than one space) in the bulk and on the surface together. We note that if d Ω is chosen the same as d Γ , there is no difference in the region corresponding to Turing space as shown in Sub-figure 9(c). In Sub-figures 9(a) and (b), it can be seen that for larger value of the diffusion coefficient the Turing space is significantly larger than that for the smaller value of the diffusion coefficient. In the context of pattern formation it means that regions corresponding to diffusion driven instability enlarge with an increase in the diffusion coefficient.   This section serves to provide the theoretical formulation required to obtain numerical solutions through the finite element method for the system that were explored in Section II. The Sobolev and Hilbert function spaces are the basis used to obtain the weak formulation. The methods of space and time discretisations were investigated and also the timestepping schemes. We present the weak formulation with the corresponding finite element formulation through a fully implicit treatment by employing the extended form of Newton's method for vector valued functions.

A. Weak Formulation
In order to derive the weak formulation, we multiply (13) by a test function say ϕ ∈ H 1 (Ω) for the bulk and ψ ∈ H 1 (Γ) for the surface and integrated over Ω for the bulk and over Γ for the surface written as Using the Green's formula for the second terms in the above with the boundary conditions (14), we obtain

B. Spatial Discretisation of the Weak Formulation
We discretise the original domain Ω and its boundary Γ to obtain Ω h and Γ h where Ω h ⊂ Ω and Γ h ⊂ Γ with N Ω and N Γ to be the number of vertices on associated to their respective discretisation. Let V Ω h and V Γ h denote the finite element function spaces associated to the discretised domains Ω h and Γ h respectively.

The finite element formulation is then to seek
are true for all test functions ϕ h ∈ V Ω h and ψ h ∈ V Γ h , respectively.
Let {ϕ i } NΩ i=1 and {ψ i } NΓ i=1 be the set of piecewise bilinear basis functions. It is known that the spaces V Ω h and V Γ h are spanned by the basis functions {ϕ i } NΩ i=1 and {ψ i } NΓ i=1 respectively. Thus, u h , v h , r h and s h may be expanded in terms of linear combinations of its corresponding basis functions S i ψ i in the finite element formulations leads to a system of differential equations written in matrix notation as where the matrices with their corresponding entries are given www.ijacsa.thesai.org by and the entries for M 1 , A 1 , B 1 (R, S) and C 1 , are expressed in similar way to those expressed for matrices with subscript 0. The entries of the matrices that are constructed from the combination of function spaces defined in the bulk and on the surface are defined by where M is the mass matrix and A is the stiffness matrix, B is the matrix corresponding to the non-linear terms and C is the column vector.
C. Mesh Generation (using Deal. II) The usual approach to discretising Ω and Γ is such that, Ω is first discretised and denoted by Ω h . The union of those elements from Ω h whose vertices lie on ∂Ω is considered as the discretisation of Γ, which is denoted by Γ h . Bulk is discretised by quadrilateral elements each with uniform structure throughout Ω h . Triangulation Γ h is also a uniform set of 2 dimensional quadrilaterals consisting of the external faces of all the bulk elements that have at least one vertex on Γ h .

D. Time Discretisation
We discretise the time interval [0, T ] into a finite number of uniform subintervals such that 0 = t 0 < t 1 · · · < t J = T . Let τ be the time steps and J be a fixed positive integer, then T = Jτ. We denote the approximate solution at time t n = nτ by u n h = u h (., t n ) where n = 0, 1, · · · , J and similar for the other variables. A fully implicit Euler scheme is used to solve the system in time. The fully implicit scheme applied to the uniform time discretisation. We can obtain the fully discretised system as Algebraic manipulation and rearrangement of each equation, leads to write the system a different form which is F 1 (U n , V n , R n , S n ) = 0, F 2 (U n , V n , R n , S n ) = 0, F 3 (U n , V n , R n , S n ) = 0, F 4 (U n , V n , R n , S n ) = 0, where In order to solve the system of non-linear equation, the employing the extended form of Newton's method for vector valued functions lead to write where, A = (u n k , v n k , r n k , s n k ) the index ı = 1, 2, 3, 4 and

V. CONCLUSION
Bulk-surface reaction-diffusion system is explored through studying non-linear reaction kinetics with linear Robin-type boundary conditions. For non-linear reaction kinetics are posed both in the bulk and on the surface, then with appropriate parameter choices, such system is able to give rise to pattern formation everywhere. Parameters can also be chosen for this system such that pattern emerges in the bulk and extends to the surface, however it forms no pattern on the internal boundary layer. It is worth noting that the emergence of no pattern in the internal boundary layer is a consequence of parameter choice in the system and not the exhaustive results associated to it.
The weak formulation of coupled bulk-surface reactiondiffusion system was obtained to set-up the premises for discretisation in space through employing the standard finite element method. The full coupled system of BSRDEs was simulated using a fully implicit time-stepping scheme through the application of extended form of Newton's method for vector valued functions. Using fully implicit time-stepping scheme, we numerically demonstrate that this system allows patterns to emerge everywhere.
The generality, robustness and applicability of the presented theoretical computational framework for coupled system of bulk-surface reaction-diffusion equations set premises to study experimentally driven models where coupling of bulk and surface chemical species is prevalent. Examples of such applications include cell motility, pattern formation in developmental biology, material science and cancer biology.
ACKNOWLEDGMENT I would like to thank Dr Chandrasekhhar Venkataraman, University of Sussex, for his initial co-supervising. www.ijacsa.thesai.org