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. 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 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 can not generate patterns everywhere in the bulk in the absence of patterning from the bulk reaction-diffusion 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, Robin-type boundary conditions seem to introduce a boundary layer coupling the bulk and surface dynamics.
Bulk-surface reaction-diffusion equations, bulk-surface finite elements, Turing diffusively-driven instability, linear stability, pattern formation, Robin-type boundary conditions
Mathematics subject classification
35K55, 35K57, 37B25, 37B55, 37C60, 37C75,
In many fluid dynamics applications and biological processes, coupled bulk-surface 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 . 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 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 . 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, 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  study symmetry breaking in a bulk-surface reaction-diffusion 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 cell-surface, a system of two membrane reaction-diffusion equations is formulated. The bulk and cell-surface membrane are coupled through Robin-type boundary conditions and a flux term for the membrane system . Elliott and Ranner  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.  present an algorithm for solving a diffusion equation on a curved surface coupled to a diffusion model in the volume. Checkkin et al.  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.  presents a detailed review on surface and bulk characterisation of electrospun membranes of porous and fibrous polymer materials. To explain the long-range proton translocation along biological mombranes, Medvedev and Stuchebrukhov  propose a model that takes into account the coupled bulk-diffusion that accompanies the migration of protons on the surface. More recently, Rozada et al.,  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 surface-elliptic or surface-diffusion or surface-reaction-diffusion models posed on the surface through Robin-type boundary conditions [3, 6, 17, 21, 22, 24, 25]. 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 [3, 6, 17, 21, 22, 24, 25, 30]. We employ the bulk-surface finite element method as introduced by Elliott and Ranner in  to numerically solve the coupled system of bulk-surface reaction-diffusion equations. Details of the surface-finite element can be found in . The bulk and surface reaction-diffusion systems are coupled through Robin-type boundary conditions. The coupled bulk-surface finite element algorithm is implemented in deall II .
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 bulk-surface reaction-diffusion equations.
Using a 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.
These theoretical predictions are supported by numerical simulations.
Hence this article is outlined as follows. In Section 2 we present the coupled bulk-surface reaction-diffusion 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 diffusion-driven instability conditions for the coupled system of bulk-surface reaction-diffusion equations. To validate our theoretical findings, we present bulk-surface finite element numerical solutions in Section 3. In Section 4, we conclude and discuss the implications of our findings.
2 Coupled bulk-surface reaction-diffusion systems on stationary volumes
In this section we present a coupled system of bulk-surface reaction-diffusion equations (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 while no boundary conditions are imposed on the surface reaction-diffusion system since the surface is closed.
2.1 A coupled system of bulk-surface reaction-diffusion 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 Laplace-Beltrami 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 cross-diffusion and assuming that coupling is only through the reaction kinetics, we propose to study the following non-dimensionalised coupled system of BSRDEs
with coupling boundary conditions
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 
where , , , , and are constant non-dimensionalised parameters. Initial conditions are given by the positive bounded functions , , and .
Activator-depleted reaction kinetics: An illustrative example
where and are positive parameters. For analytical simplicity, we postulate the model system (2.1) in a more compact form given by
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).
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 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 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
where represents the Jacobian matrix representing the first linear terms of the linearization process. Its entries are defined by
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
where is the identity matrix. Making appropriate substitutions and carrying out standard calculations we obtain the following dispersion relation for
For convenience’s sake, let us denote by
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
The proof enforces that is a Hurwitz polynomial and therefore satisfies the Routh-Hurwitz criterion in order for . The first condition to be satisfied is that otherwise is a trivial root, thereby reducing the 4-th 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.
We require that the determinant of the matrix
Substituting , and appropriately we obtain
Exploiting the fact that
it then follows that
if and only if
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
with linearised 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
For the case of circular and spherical domains, if , then .
where . Writing
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
where we have defined conveniently