Stability analysis and simulations of coupled bulk-surface reaction–diffusion systems

In this article, we formulate new models for coupled systems of bulk-surface reaction–diffusion equations on stationary volumes. The bulk reaction–diffusion equations are coupled to the surface reaction–diffusion equations through linear Robin-type boundary conditions. We then state and prove the necessary conditions for diffusion-driven instability for the coupled system. Owing to the nature of the coupling between bulk and surface dynamics, we are able to decouple the stability analysis of the bulk and surface dynamics. Under a suitable choice of model parameter values, the bulk reaction–diffusion system can induce patterning on the surface independent of whether the surface reaction–diffusion system produces or not, patterning. On the other hand, the surface reaction–diffusion system cannot generate patterns everywhere in the bulk in the absence of patterning from the bulk reaction–diffusion system. For this case, patterns can be induced only in regions close to the surface membrane. Various numerical experiments are presented to support our theoretical findings. Our most revealing numerical result is that, Robin-type boundary conditions seem to introduce a boundary layer coupling the bulk and surface dynamics.


Introduction
In many fluid dynamics applications and biological processes, coupled bulk-surface partial differential equations naturally arise in (2D + 3D) [1][2][3]. In most of these applications and processes, morphological instabilities occur through symmetry breaking resulting in the formation of heterogeneous distributions of chemical substances [4]. In developmental biology, it is essential for the emergence and maintenance of polarized states in the form of heterogeneous distributions of chemical substances such as proteins and lipids. Examples of such processes include (but are not limited to) the formation of buds in yeast cells, and cell polarization in biological cells owing to responses to external signals through the outer cell membrane [5,6]. In the context of reaction-diffusion processes, such symmetry breaking arises when a uniform steady state, stable in the absence of diffusion, is driven unstable when diffusion is present thereby giving rise to the formation of spatially inhomogeneous solutions in a process now well known as the Turing diffusion-driven instability [7]. Classical Turing theory requires that one of the chemical species, typically the inhibitor, diffuses much faster than the other, the activator resulting in what is known as the long-range inhibition and short-range activation [8,9].
Recently, there has been a surge in studies on models that couple bulk dynamics to surface dynamics. For example, Rätz & Röger [6] study symmetry breaking in a bulk-surface reactiondiffusion model for signalling networks. In this work, a single diffusion partial differential equation (the heat equation) is formulated inside the bulk of a cell, whereas on the cell surface, a system of two membrane reaction-diffusion equations is formulated. The bulk and cellsurface membrane are coupled through Robin-type boundary conditions and a flux term for the membrane system [6]. Elliott & Ranner [10] study a finite-element approach to a sample elliptic problem: a single elliptic partial differential equation is posed in the bulk, and another is posed on the surface. These are then coupled through Robin-type boundary conditions. Novak et al. [11] present an algorithm for solving a diffusion equation on a curved surface coupled to a diffusion model in the volume. Chechkin et al. [12] study bulk-mediated diffusion on planar surfaces. Again, diffusion models are posed in the bulk and on the surface coupling them through boundary conditions. In the area of tissue engineering and regenerative medicine, electrospun membrane are useful in applications such as filtration systems and sensors for chemical detection. Understanding of the fibres' surface, bulk and architectural properties is crucial to the successful development of integrative technology. Nisbet et al. [13] present a detailed review on surface and bulk characterization of electrospun membranes of porous and fibrous polymer materials. To explain the long-range proton translocation along biological mombranes, Medvedev & Stuchebrukhov [14] propose a model that takes into account coupled bulk diffusion that accompanies the migration of protons on the surface. More recently, Rozada et al. [15] present singular perturbation theory for the stability of localized spot patterns for the Brusselator model on the sphere.
In most of the work above, either elliptic or diffusion models in the bulk have been coupled to surface-elliptic or surface-diffusion or surface-reaction-diffusion models posed on the surface through Robin-type boundary conditions [5,6,[11][12][13][14]16]. Here, our focus is to couple systems of reaction-diffusion equations posed both in the bulk and on the surface, setting a mathematical and computational framework to study more complex interactions such as those observed in cell biology, tissue engineering and regenerative medicine, developmental biology and biopharmaceuticals [5,6,[11][12][13][14]16,17]. We employ the bulk-surface finite-element method as introduced by Elliott & Ranner in [10] to numerically solve the coupled system of bulk-surface reaction-diffusion equations (BSRDEs). Details of the surface-finite-element can be found in reference [18]. The bulk and surface reaction-diffusion systems are coupled through Robin-type boundary conditions. Details of the coupled bulk-surface finite-element method can be found in [19]; the finite element algorithm is implemented in deall II [20].
The key contributions of our work to the theory of pattern formation are: -we derive and prove Turing diffusion-driven instability conditions for a coupled system of BSRDEs; -using the bulk-surface finite-element method, we approximate the solution to the model system within the bulk and on the boundary surface of a sphere of radius one; . -our results show that if the surface-reaction-diffusion system has the long-range inhibition, short-range activation form and the bulk-reaction-diffusion system has equal diffusion coefficients, then the surface-reaction-diffusion system can induce patterns in the bulk close to the surface and no patterns form in the interior, far away from the surface; -on the other hand, if the bulk-reaction-diffusion system has the long-range inhibition, short-range activation form and the surface-reaction-diffusion system has equal diffusion coefficients, then the bulk-reaction-diffusion system can induce pattern formation on the surface; -furthermore, we prove that if the bulk and surface reaction-diffusion systems have equal diffusion coefficients, no patterns form; and -these theoretical predictions are supported by numerical simulations.
Hence, this article is outlined as follows. In §2, we present the coupled bulk-surface reactiondiffusion system on stationary volumes with appropriate boundary conditions coupling the bulk and surface partial differential equations. The main results of this article are presented in §2b where we derive Turing diffusion-driven instability conditions for the coupled system of BSRDEs. To validate our theoretical findings, we present bulk-surface finite-element numerical solutions in §3. In §4, we conclude and discuss the implications of our findings.

Coupled bulk-surface reaction-diffusion systems on stationary volumes
Here, we present a coupled system of BSRDEs posed in a three-dimensional volume as well as on the boundary surface enclosing the volume. We impose Robin-type boundary conditions on the bulk reaction-diffusion system, whereas no boundary conditions are imposed on the surface reaction-diffusion system since the surface is closed.

(a) A coupled system of bulk-surface reaction-diffusion equations
Let Ω be a stationary volume (whose interior is denoted the bulk) enclosed by a compact hypersurface Γ := ∂Ω which is C 2 . In addition, let I = [0, T](T > 0) be some time interval. Moreover, let ν denote the unit outer normal to Γ , and let U be any open subset of R N+1 containing Γ , then for any function u which is differentiable in U, we define the tangential gradient on Γ by, ∇ Γ u = ∇u − (∇u · ν) ν, where · denotes the regular dot product and ∇ denotes the regular gradient in R N+1 . The tangential gradient is the projection of the regular gradient onto the tangent plane, thus ∇ Γ u · ν = 0. The Laplace-Beltrami operator on the surface Γ is then defined to be the tangential divergence of the tangential gradient Γ u = ∇ Γ · ∇ Γ u. For a vector function u = (u 1 , u 2 , . . . , u N+1 ) ∈ R N+1 , the tangential divergence is defined by To proceed, we denote by u : Ω × I → R and v : Ω × I → R two chemical concentrations (species) that react and diffuse in Ω and r : Γ × I → R and s : Γ × I → R be two chemical species residing only on the surface Γ which react and diffuse on the surface. In the absence of crossdiffusion and assuming that coupling is only through the reaction kinetics, we propose to study the following non-dimensionalized coupled system of BSRDEs In the above, ∇ 2 = ∂ 2 /∂x 2 + ∂ 2 /∂y 2 + ∂ 2 /∂z 2 represents the Laplacian operator. d Ω and d Γ are a positive diffusion coefficients in the bulk and on the surface, respectively, representing the ratio between u and v, and r and s, respectively. γ Ω and γ Γ represent the length-scale parameters in the bulk and on the surface, respectively. In this formulation, we assume that f (·, ·) and g(·, ·) are nonlinear reaction kinetics in the bulk and on the surface. h 1 (u, v, r, s) and h 2 (u, v, r, s) are reactions representing the coupling of the internal dynamics in the bulk Ω to the surface dynamics on the surface Γ . As a first attempt, we consider a more generalized form of linear coupling of the following nature [21] and where α 1 , α 2 , β 1 , β 2 , κ 1 and κ 2 are constant non-dimensionalized parameters. Initial conditions are given by the positive-bounded functions u 0 (x), v 0 (x), r 0 (x) and s 0 (x).
(i) Activator-depleted reaction kinetics: an illustrative example From now onwards, we restrict our analysis and simulations to the well-known activator-depleted substrate reaction model [8,[22][23][24][25] also known as the Brusselator given by where a and b are positive parameters. For analytical simplicity, we postulate the model system (2.1) in a more compact form given by ⎧ ⎨ with coupling boundary conditions (2.2)-(2.4). In the above, we have defined appropriately and provided the following compatibility condition on the coefficients of the coupling is satisfied Proof. The proof follows immediately from the definition of the uniform steady state satisfying reaction kinetics (2.7)-(2.10). It must be noted that in deriving this unique uniform steady state the compatibility condition (2.12) coupling bulk and surface dynamics must be satisfied.

Remark 2.4.
Note that there exists an infinite number of solutions to problem (2.12).

(i) Linear stability analysis in the absence of diffusion
Next, we study Turing diffusion-driven instability for the coupled system of BSRDEs (2.1)-(2.4) with reaction kinetics (2.5). To proceed, we first consider the linear stability of the spatially uniform steady state. For the sake of convenience, let us denote by w = (u, v, r, s) T , the vector of the species u, v, r and s. Furthermore, defining the vector ξ such that |ξ i | 1 for all i = 1, 2, 3 and 4, it follows that writing w = w * + ξ , the linearized system of coupled BSRDEs can be posed as where J F represents the Jacobian matrix representing the first linear terms of the linearization process. Its entries are defined by where by definition f 1 u := ∂f 1 /∂u represents a partial derivative of f 1 (u, v) with respect to u. We are looking for solutions to the system of linear ordinary differential equations (2.13) which are of the form ξ ∝ e λt . Substituting into (2.13), results in the following classical eigenvalue problem where I is the identity matrix. Making appropriate substitutions and carrying out standard calculations, we obtain the following dispersion relation for λ and For the sake of convenience, let us denote by the submatrices of J F corresponding to the bulk reaction kinetics and the surface reaction kinetics, respectively. We can now define Theorem 2.5 (Necessary and sufficient conditions for Re(λ) < 0). The necessary and sufficient conditions such that the zeros of the polynomial p 4 (λ) have Re(λ) < 0 are given by the following conditions and Proof. The proof enforces that p 4 (λ) is a Hurwitz polynomial and therefore satisfies the Routh-Hurwitz criterion in order for Re(λ) < 0. The first condition to be satisfied is that a 4 = 0 otherwise λ = 0 is a trivial root, thereby reducing the fourth-order polynomial to a cubic polynomial. The first four conditions are a result of requiring that each coefficient a i with i = 1, 2, 3 and 4 of the polynomial p 4 (λ) are all positive. The rest of the conditions are derived as shown below.
We require that the determinant of the matrix Substituting a 1 , a 2 and a 3 appropriately, we obtain Exploiting the fact that it then follows that Multiplying throughout by 2, we obtain condition (2.25) in theorem 2.5. The last condition results from imposing the condition that It can be shown that On the other hand, Hence, combining (2.28) and (2.29) and simplifying conveniently, we have resulting in condition (2.26).
Remark 2.6. The characteristic polynomial, can also be written more compactly in the form of thereby coupling the bulk and surface dispersion relations in the absence of spatial variations.
(ii) Linear stability analysis in the presence of diffusion , we obtain a linearized system of partial differential equations with linearized boundary conditions In the above, we have used the original reaction kinetics for the purpose of clarity.
In order to proceed, we restrict our analysis to circular and spherical domains where we can transform the cartesian coordinates into polar coordinates and be able to exploit the method of separation of variables. Without loss of generality, we write the following eigenvalue problem in the bulk Taking x ∈ B, y ∈ Γ , then writing in polar coordinates x = ry, r ∈ (0, 1) we can define, for all l ∈ N 0 , m ∈ Z, |m| ≤ l, the following power-series solutions [5,6]  where where we have defined conveniently Det(M) Γ := d Γ l 2 (l + 1) 2 − γ Γ (d Γ f r + g s )l(l + 1) + γ 2 Γ (f r g s − f s g r ).
The above holds true if and only if either In the presence of diffusion, we require the emergence of spatial growth. In order for the uniform steady state w * to be unstable, we require that either 1. Re(λ l,m (k 2 l,m )) > 0 for some k 2 l,m > 0, or 2. Re(λ l,m (l(l + 1))) > 0 for some l(l + 1) > 0, or 3. Both.
Solving (2.49) (and similarly (2.50)), we obtain the eigenvalues Similarly, on the surface, Re(λ l,m (l(l + 1))) > 0 for some l(l + 1) > 0 if and only the following conditions hold We are in a position to state the theorem 2.8.

56)
then the necessary conditions for Re(λ l,m (k 2 l,m )) > 0 for some k 2 l,m > 0 are given by

Similarly, assuming that
Tr(J F ) Γ = f r + g s < 0 and Det(J F ) Γ = f r g s − f s g r > 0, (2.58) then the necessary conditions for Re(λ l,m (l(l + 1))) > 0 for some l(l + 1) > 0 are given by Proof. The proof is a direct consequence of conditions (2.52)-(2.55). Assuming that conditions (2.56) and (2.58) hold, then one of the conditions in (2.52) and (2.54) is violated, which implies that Re(λ l,m (k 2 l,m )) < 0 for all k 2 l,m > 0 and similarly Re(λ l,m (l(l + 1))) < 0 for all l(l + 1) > 0. This entails that the system can no longer exhibit spatially inhomogeneous solutions. The only two conditions left to hold true are (2.53) and (2.55). This case corresponds to the classical standard two-component reaction-diffusion system which requires that (for details, see for example [9]) (2.60) and similarly and Tr(J F ) Γ = f r + g s < 0 and Det(J F ) Γ = f r g s − f s g r < 0.
Then, it follows that condition (2.25) given by -Similarly, the case when and Tr(J F ) Γ = f r + g s < 0 and Det(J F ) Γ = f r g s − f s g r < 0.
This implies that none of the conditions (2.21)-(2.26) are violated, while condition (2.54) fails to hold. -Finally, the cases when Tr(J F ) Γ = f r + g s < 0 and Det(J F ) Γ = f r g s − f s g r > 0, (2.62) and and Tr(J F ) Γ = f r + g s > 0 and Det(J F ) Γ = f r g s − f s g r > 0, (2.63) result in remark 2.10.
The above cases clearly eliminate conditions (2.52) and (2.54) as necessary for uniform steady state to be driven unstable in the presence of diffusion. We are now in a position to state our main result.

Theorem 2.13 (Necessary conditions for diffusion-driven instability). The necessary conditions for diffusion-driven instability for the coupled system of BSRDEs (2.1)-(2.4) are given by
and

(iii) Theoretical predictions
From the analytical results, we state the following theoretical predictions to be validated through the use of numerical simulations.
1. The bulk and surface diffusion coefficients d Ω and d Γ must be greater than one in order for diffusion-driven instability to occur. Taking d Ω = d Γ = 1 results in a contradiction between conditions (2.64), (2.70) and (2.71). As a result, the BSRDEs does not give rise to the formation of spatial structure. For this case, the uniform steady state is the only stable solution for the coupled system of BSRDEs (2.1)-(2.4). 2. The above imply that taking d Ω > 1 and d Γ = 1, the bulk-reaction-diffusion system has the potential to induce patterning in the bulk for appropriate diffusion-driven instability parameter values, whereas the surface-reaction-diffusion system is not able to generate patterns. Here, all the conditions (

3.
Alternatively, taking d Ω = 1 and d Γ > 1, the bulk-reaction-diffusion system fails to induce patterning in the bulk while the surface-reaction-diffusion system has the potential to induce patterning on the surface. Similarly, all the conditions (2.64)-(2.71) hold except (2.70). 4. On the other hand, taking d Ω > 1 and d Γ > 1 appropriately, then the coupled system of BSRDEs exhibits patterning both in the bulk and on the surface. All the conditions (2.64)-(2.71) hold both in the bulk and on the surface.

Numerical simulations of the coupled system of bulk-surface reaction-diffusion equations
Here, we present bulk-surface finite-element numerical solutions corresponding to the coupled system of BSRDEs (2.1)-(2.5). Here, we omit the details of the bulk-surface finite-element method as these are given elsewhere (see [19] for details). Our method is inspired by the work of Elliott & Ranner [10]. We use the bulk-surface finite-element method to discretize in space with piecewise bilinear elements and an implicit second-order fractional-step θ-scheme to discretize in time using the Newton's method for the linerization [19,26]. For details on the convergence and stability of the fully implicit time-stepping fractional-step θ-scheme, the reader is referred to Madzvamuse et al. [19,26]. In all our numerical experiments, we fix the kinetic model parameter values a = 0.1 and b = 0.9 since these satisfy the Turing diffusion-driven instability conditions (2.64)-(2.71). For these model parameter values, the equilibrium values are (u * , v * , r * , s * ) = (1, 0.9, 1, 0.9). Initial conditions are prescribed as small random perturbations around the equilibrium values. For illustrative purposes, let us take the parameter values describing the boundary conditions as shown in table 1; these are selected such that they satisfy the compatibility condition (2.12).

(a) A note on the bulk-surface triangulation
We briefly outline how the bulk-surface triangulation is generated. For further specific details, please see reference [19]. Let Ω h be a triangulation of the bulk geometry Ω with vertices  figure 1. The bulk triangulation induces the surface triangulation as illustrated.

(b) Numerical experiments
Here, we will present only four cases to validate our theoretical predictions outlined in §2b(iii). In most of our simulations, parameter values are fixed as shown in table 1, except for d Ω and d Γ whose values are varied to demonstrate the patterning mechanism of the coupled system of BSRDEs. We present patterns corresponding only to the chemical species u and r in the bulk and on the surface, respectively. Those corresponding to v and s are 180 degrees out of phase to those of u and r and are therefore omitted. It must be noted however that this is not always the case in general, Robin-type boundary conditions may alter the structure of the solution profiles depending on the model parameter values and the coupling compatibility boundary parameters. The bulk-surface finite-element numerical simulations of the coupled system of BSRDEs with d Ω = 1 in the bulk, d Γ = 1 on the surface are shown in figure 2. We observe that no patterns form in complete agreement with theoretical predictions. Similar to classical reaction-diffusion systems, diffusion coefficients must be greater than one. In particular, the diffusion coefficients must be greater than their corresponding respective critical diffusion coefficients in the bulk and on the surface. An example is shown next.   of the surface membrane. Spots are observed to form on the surface, whereas in the bulk, small balls form inside. Far away from the surface, no patterns form, because the necessary conditions for diffusion-driven instability are not fulfilled in the bulk. These results confirm our theoretical predictions. We note that this particular example describes realistically pattern formation in biological systems. We expect skin patterning to manifest in the epidermis layer as well as on the surface.
(iii) Simulations of the coupled system of BSRDEs with (d Ω , d Γ ) = (20, 1) To generate patterns in the bulk, we take d Ω = 20 > d crit Γ = 8.5 and d Γ = 1 on the surface. Figure 4 exhibits stripe, circular and spot patterns in the bulk as illustrated by the cross sections. On the surface, small amplitude patterns occur consistent with theoretical predictions. Although the patterns for the u species (columns one and two) appear uniform on the surface this is simply owing to the colour scale, with the amplitude of the patterns in the bulk larger than those on the surface. This difference in the amplitude of the pattern of the bulk solution in the bulk and on the surface is due to the Robin-type boundary conditions. Unlike zero-flux (also known as homogeneous Neumann), boundary conditions for standard reaction-diffusion systems which imply that no species enter or leave the domain, here, there is deposition or removal of chemical species through the flux on the surface, resulting in differences in amplitude between the bulk and surface solutions.
(iv) Simulations of the coupled system of BSRDEs with (d Ω , d Γ ) = (20,20) In this example, we illustrate how both bulk and surface dynamics induce patterning by taking d Ω = 20 in the bulk, d Γ = 20 on the surface. Figure 5 shows pattern formation in the bulk and on  the surface. In the bulk, we observe the formation of balls (which can be seen as spots through cross sections) and these translate to spots on the surface. The surface dynamics themselves induce spot pattern formation.

Conclusion, discussion and future research challenges
We have presented a coupled system of BSRDEs whereby the bulk and surface reaction-diffusion systems are coupled through Robin-type boundary conditions. Nonlinear reaction-kinetics are considered in the bulk and on the surface and for illustrative purposes, the activator-depleted model was selected because it has a unique positive steady state. By using linear stability theory close to the bifurcation point, we state and prove a generalization of the necessary conditions for Turing diffusion-driven instability for the coupled system of BSRDEs. Our most revealing result is that the bulk reaction-diffusion system has the capability of inducing patterning (under appropriate model and compatibility parameter values) for the surface reaction-diffusion model. On the other hand, the surface reaction-diffusion is not capable of inducing patterning everywhere in the bulk; patterns can be induced in only regions close to the surface membrane. For skin pattern formation, this example is consistent with the observation that patterns will form on the surface as well as within the epidermis layer close to the surface. We do not expect patterning to form everywhere in the body of the animals. Our studies reveal the following observations and research questions still to be addressed: -our numerical experiments reveal that the Robin-type boundary conditions seem to introduce a boundary layer coupling the bulk and surface dynamics. However, these boundary conditions do not appear explicitly in the conditions for diffusion-driven instability and this makes it difficult to theoretically analyse their role and implications to pattern formation. Further studies are required to understand the role of these boundary conditions as well as the size of the boundary layer; -the compatibility condition (2.12) implies that the uniform steady state in the bulk is identical to the uniform state on the surface. We are currently studying the implications of relaxing the compatibility condition; -finally, in this manuscript, we have not carried out detailed parameter search and estimation to deduce the necessary and sufficient conditions for pattern generation as well as isolating excitable wavenumbers in the bulk and on the surface. Such studies might reveal more interesting properties of the coupled bulk-surface model and this forms part of our current studies.
We have presented a framework that couples bulk dynamics (three-dimension) to surface dynamics (two-dimension) with the potential of numerous applications in cell motility, developmental biology, tissue engineering and regenerative medicine and biopharmaceutical where reaction-diffusion-type models are routinely used [5,6,[11][12][13][14]16,17].
We have restricted our studies to stationary volumes. In most cases, biological surfaces are known to evolve continuously with time. This introduces extra complexities to the modelling, analysis and simulation of coupled systems of BSRDEs. In order to consider evolving bulk-surface partial differential equations, evolution laws (geometrical) should be formulated describing how the bulk and surface evolve. Here, it is important to consider specific experimental settings that allow for detailed knowledge of properties (biomechanical) and processes (biochemical) involved in the bulk-surface evolution. Such a framework will allows us to study three-dimensional cell migration in the area of cell motility [16, [27][28][29]. In future studies, we propose to develop a three-dimensional integrative model that couples bulk and surface dynamics during growth development or movement.
Data accessibility. This manuscript does not contain primary data and as a result has no supporting material associated with the results presented.