Stability analysis and simulations of coupled bulksurface reactiondiffusion systems
Abstract.
In this article we formulate new models for coupled systems of bulksurface reactiondiffusion equations on stationary volumes. The bulk reactiondiffusion equations are coupled to the surface reactiondiffusion equations through linear Robintype boundary conditions. We then state and prove the necessary conditions for diffusiondriven instability for the coupled system. Due 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 reactiondiffusion system can induce patterning on the surface independent of whether the surface reactiondiffusion system produces or not, patterning. On the other hand, the surface reactiondiffusion system can not generate patterns everywhere in the bulk in the absence of patterning from the bulk reactiondiffusion system. For this case, patterns can only be induced in regions close to the surface membrane. Various numerical experiments are presented to support our theoretical findings. Our most revealing numerical result is that, Robintype boundary conditions seem to introduce a boundary layer coupling the bulk and surface dynamics.
Key words.
Bulksurface reactiondiffusion equations, bulksurface finite elements, Turing diffusivelydriven instability, linear stability, pattern formation, Robintype boundary conditions
Mathematics subject classification
35K55, 35K57, 37B25, 37B55, 37C60, 37C75,
1 Introduction
In many fluid dynamics applications and biological processes, coupled bulksurface partial differential equations naturally arise in (). In most of these applications and processes, morphological instabilities occur through symmetry breaking resulting in the formation of heterogeneous distributions of chemical substances [11]. In developmental biology, it is essential the emergence and maintenance of polarised 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 polarisation in biological cells due to responses to external signals through the outer cell membrane [24, 25]. In the context of reactiondiffusion 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 wellknown as the Turing diffusiondriven instability [28]. 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 longrange inhibition and shortrange activation [8, 18].
Recently, there has been a surge in studies on models that coupled bulk dynamics to surface dynamics. For example, Rätz and Röger [25] study symmetry breaking in a bulksurface 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, while on the cellsurface, a system of two membrane reactiondiffusion equations is formulated. The bulk and cellsurface membrane are coupled through Robintype boundary conditions and a flux term for the membrane system [25]. Elliott and Ranner [5] 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 Robintype boundary conditions. Novak et al. [22] present an algorithm for solving a diffusion equation on a curved surface coupled to a diffusion model in the volume. Checkkin et al. [3] study bulkmediated 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. [21] presents a detailed review on surface and bulk characterisation of electrospun membranes of porous and fibrous polymer materials. To explain the longrange proton translocation along biological mombranes, Medvedev and Stuchebrukhov [17] propose a model that takes into account the coupled bulkdiffusion that accompanies the migration of protons on the surface. More recently, Rozada et al., [26] presented singular perturbation theory for the stability of localised 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 surfaceelliptic or surfacediffusion or surfacereactiondiffusion models posed on the surface through Robintype boundary conditions [3, 6, 17, 21, 22, 24, 25]. Here, our focus is to couple systems of reactiondiffusion 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 [3, 6, 17, 21, 22, 24, 25, 30]. We employ the bulksurface finite element method as introduced by Elliott and Ranner in [5] to numerically solve the coupled system of bulksurface reactiondiffusion equations. Details of the surfacefinite element can be found in [4]. The bulk and surface reactiondiffusion systems are coupled through Robintype boundary conditions. The coupled bulksurface finite element algorithm is implemented in deall II [1].
The key contributions of our work to the theory of pattern formation are:

We derive and prove Turing diffusiondriven instability conditions for a coupled system of bulksurface reactiondiffusion equations.

Using a bulksurface 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 surfacereactiondiffusion system has the longrange inhibition, shortrange activation form and the bulkreactiondiffusion system has equal diffusion coefficients, then the surfacereactiondiffusion 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 bulkreactiondiffusion system has the longrange inhibition, shortrange activation form and the surfacereactiondiffusion system has equal diffusion coefficients, then the bulkreactiondiffusion system can induce pattern formation on the surface.

Furthermore, we prove that if the bulk and surface reactiondiffusion systems have equal diffusion coefficients, no patterns form.

These theoretical predictions are supported by numerical simulations.
Hence this article is outlined as follows. In Section 2 we present the coupled bulksurface 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 Section 2.2 where we derive Turing diffusiondriven instability conditions for the coupled system of bulksurface reactiondiffusion equations. To validate our theoretical findings, we present bulksurface finite element numerical solutions in Section 3. In Section 4, we conclude and discuss the implications of our findings.
2 Coupled bulksurface reactiondiffusion systems on stationary volumes
In this section we present a coupled system of bulksurface reactiondiffusion equations (BSRDEs) posed in a threedimensional volume as well as on the boundary surface enclosing the volume. We impose Robintype boundary conditions on the bulk reactiondiffusion system while no boundary conditions are imposed on the surface reactiondiffusion system since the surface is closed.
2.1 A coupled system of bulksurface reactiondiffusion equations (BSRDEs)
Let be a stationary volume (whose interior is denoted the bulk) enclosed by a compact hypersurface which is . Also, let be some time interval. Moreover, let denote the unit outer normal to , and let be any open subset of containing , then for any function which is differentiable in , we define the tangential gradient on by, where denotes the regular dot product and denotes the regular gradient in . The tangential gradient is the projection of the regular gradient onto the tangent plane, thus . The LaplaceBeltrami operator on the surface is then defined to be the tangential divergence of the tangential gradient . For a vector function the tangential divergence is defined by
To proceed, we denote by and two chemical concentrations (species) that react and diffuse in and and and 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 nondimensionalised coupled system of BSRDEs
(2.1) 
with coupling boundary conditions
(2.2) 
In the above, represents the Laplacian operator. and are a positive diffusion coefficients in the bulk and on the surface respectively, representing the ratio between and , and and , respectively. and represent the length scale parameters in the bulk and on the surface respectively. In this formulation, we assume that and are nonlinear reaction kinetics in the bulk and on the surface. and are reactions representing the coupling of the internal dynamics in the bulk to the surface dynamics on the surface . As a first attempt, we will consider a more generalised form of linear coupling of the following nature [12]
(2.3)  
(2.4) 
where , , , , and are constant nondimensionalised parameters. Initial conditions are given by the positive bounded functions , , and .
Activatordepleted reaction kinetics: An illustrative example
From now onwards, we restrict our analysis and simulations to the wellknown activatordepleted substrate reaction model [8, 10, 23, 27, 29] also known as the Brusselator given by
(2.5) 
where and are positive parameters. For analytical simplicity, we postulate the model system (2.1) in a more compact form given by
(2.6) 
with coupling boundary conditions (2.2)(2.4). In the above, we have defined appropriately
(2.7)  
(2.8)  
(2.9)  
(2.10) 
2.2 Linear stability analysis of the coupled system of BSRDEs
Definition 2.1 (Uniform steady state).
Proposition 2.1 (Existence and uniqueness of the uniform steady state).
Proof.
Remark 2.1.
Remark 2.2.
Note that there exists an infinite number of solutions to problem (2.12).
Linear stability analysis in the absence of diffusion
Next, we study Turing diffusiondriven 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 convenience’s sake, let us denote by the vector of the species , , and . Furthermore, defining the vector such that for all , , and , it follows that writing the linearized system of coupled BSRDEs can be posed as
(2.13) 
where represents the Jacobian matrix representing the first linear terms of the linearization process. Its entries are defined by
(2.14) 
where by definition represents a partial derivative of with respect to . We are looking for solutions to the system of linear ordinary differential equations (2.13) which are of the form . Substituting into (2.13), results in the following classical eigenvalue problem
(2.15) 
where is the identity matrix. Making appropriate substitutions and carrying out standard calculations we obtain the following dispersion relation for
where
(2.16)  
(2.17)  
(2.18)  
(2.19) 
For convenience’s sake, let us denote by
(2.20) 
the submatrices of corresponding to the bulk reaction kinetics and the surface reaction kinetics respectively. We can now define
Theorem 2.1 (Necessary and sufficient conditions for ).
The necessary and sufficient conditions such that the zeros of the polynomial have are given by the following conditions
(2.21)  
(2.22)  
(2.23)  
(2.24)  
(2.25)  
(2.26) 
Proof.
The proof enforces that is a Hurwitz polynomial and therefore satisfies the RouthHurwitz criterion in order for . The first condition to be satisfied is that otherwise is a trivial root, thereby reducing the 4th order polynomial to a cubic polynomial. The first four conditions are a result of requiring that each coefficient with , , and of the polynomial are all positive. The rest of the conditions are derived as shown below.
Remark 2.3.
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.
Linear stability analysis in the presence of diffusion
Next we introduce spatial variations and study under what conditions the uniform steady state is linearly unstable. We linearize around the uniform steady state by taking small spatially varying perturbations of the form
(2.30) 
Substituting (2.30) into the coupled system of BSRDEs (2.1)(2.4) with reaction kinetics (2.5) we obtain a linearized system of partial differential equations
(2.31)  
(2.32)  
(2.33)  
(2.34) 
with linearised boundary conditions
(2.35)  
(2.36) 
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
(2.37) 
where each satisfies the boundary conditions (2.35) and (2.36). On the surface the eigenvalue problem is posed as
(2.38) 
Remark 2.4.
For the case of circular and spherical domains, if , then .
Taking , , then writing in polar coordinates , we can define, for all , , , the following power series solutions [24, 25]
(2.39)  
(2.40) 
On the surface, substituting the power series solutions (2.40) into (2.33) and (2.34) we have
(2.41)  
(2.42) 
Similarly, substituting the power series solutions (2.39) into the bulk equations (2.31) and (2.32) we obtain the following system of ordinary differential equations
(2.43)  
(2.44) 
Equations (2.43) and (2.44) are supplemented with boundary conditions
(2.45)  
(2.46) 
where . Writing
and substituting into the system of ordinary differential equations (2.41)(2.44), we obtain the following eigenvalue problem
(2.47) 
where
and
Note that the boundary conditions (2.45) and (2.46) have been applied appropriately to the surface linearised reactiondiffusion equations. Since
it follows that the coefficient matrix must be singular, hence we require that
Straight forward calculations show that the eigenvalue solves the following dispersion relation written in compact form as
(2.48) 
where we have defined conveniently