Abstract
The purpose of this article is to study numerically the Turing diffusion-driven instability mechanism for pattern formation on curved surfaces embedded in , specifically the surface of the sphere and the torus with some well-known kinetics. To do this, we use Euler’s backward scheme for discretizing time. For spatial discretization, we parameterize the surface of the torus in the standard way, while for the sphere, we do not use any parameterization to avoid singularities. For both surfaces, we use finite element approximations with first-order polynomials.
1. Introduction
Reaction-diffusion systems of partial differential equations (PDEs) have been used in the modeling of various processes and systems in physics and financial mathematics and especially in chemistry and biology. An attractive property of these systems is that they give rise to spatiotemporal patterns through Turing’s diffusion-driven instability mechanism [1, 2]. For very particular values of reaction rates, diffusion constants, and boundary and initial conditions, spatiotemporal patterns appear as a system solution. This mechanism has been proposed to model complex patterns that appear in various phenomena, for instance, patterns in the skin of some animals and morphogenesis. These systems have been so well studied, both analytically and numerically, that it is possible to know the type of pattern that will appear depending on the kind of reaction that is carried out, how they change according to the initial and boundary conditions, and other system parameters. The well-known models of Schnakenberg [3] and Gierer and Meinhardt [4] are used to explain the activator-inhibitor relationship between chemical substances such as population cycles and metabolic regulation [5]. The FitzHugh-Nagumo model has been used in neurophysiology and in the theory of the nuclear reactor among others [6]. An interesting system is the BVAM model that has been used to study transitions between fish patterns that go from stripes to spots [7]. In this way, we can talk about many more different types of reactions.
In recent works, it has been seen that Turing conditions can strongly depend on the substrate shape where the chemicals diffuse. Modification in the range of unstable modes appears for curved surfaces [8] (and references therein) and systems with advection [9], on growing surfaces [10–12]; there is also a difference if the membranes are thin [13] or elastic [14]; the stability of traveling waves have been found to also depend on surface curvature [15], etc. Patterns have also been found to tend to accommodate at certain positions on the surface, according to the curvature value [16, 17]. This could be useful for the description of processes in the cell membrane, in some tissues, or in general development processes in organisms. Indeed, for phyllotaxis, it is observed that the chemical patterns associated with the plant growth change due to stress modifications on the surface [18]. All these characteristics are typical of biological systems, pointing out the importance of domain geometry in PDE models.
In many of the previous references, the system of equations is solved numerically using the finite difference method; in some cases, the finite element is used but with automatic meshers. Recently, there are some proposals for extensions of finite element methods to solve this kind of systems [19, 20].
In this work, we solve a reaction-diffusion system on a curved surface embedded in three-dimensional space using its variational formulation and the finite elements. In order to solve this kind of problems, we present the spatial and temporal discretization of the surfaces in . As particular cases, we study discretization of the torus and the sphere separately. First, we chose the surface of the circular torus since it is easier to parameterize and triangulate, see [21–24], although it is not as common to find it in industrial and scientific applications as on the sphere. On the other hand, it presents interesting characteristics since its curvature depends not only on its radius, such as that of the sphere, but on one of its angles, which even allows changing the sign in some areas (see [25]). Finally, we study the numerical solution of a reaction-diffusion system for PDE (Turing mechanism for pattern formation) that takes place in the spherical surfaces in ; for this, it is necessary to consider Cartesian coordinates (unlike the case of the torus) to represent a sphere in order to avoid the singularities introduced by the parameterization with spherical coordinates (see [26, 27]). We study several typical reactions on both surfaces and verify the formation of Turing patterns for the corresponding parameter values.
Variational methods provide a handy analytical tool to study PDEs’ properties, for instance, to study the existence, uniqueness, and stability of solutions, differentiability properties, and spectral theory theorems. For a given differential equation, it is possible to construct a functional so that its minimum becomes a solution to the original problem [28]. In [29], Costa shows how variational methods give exact solutions to some differential equations. Recent works show that variational tools provide nontrivial smooth solutions to nonlinear problems with reactions and different boundary conditions, such as the parametric Robin problem with arbitrary potential [30] or the elliptic Dirichlet problem driven by an anisotropic Laplacian [31]. Moreover, the study of nonlinear dissipative models has also been done using group classification techniques that allow analyzing symmetries and reducing equations to obtain exact solutions to such differential equations [32, 33]. Variational methods are also used to find approximate solutions to partial differential equations in some numerical methods [21]. As already mentioned above, we will use the corresponding variational formulation of the reaction-diffusion system to apply the finite element method and obtain numerical solutions.
The paper is structured as follows. Section 2 briefly summarizes the pattern formation mechanism in reaction-diffusion systems. In the Section 3, we establish the corresponding system of equations on the embedded surface and propose its variational formulation. In Section 4, we present in detail the spatial and temporal discretization of the reaction-diffusion system. This provides the implementation of the numerical scheme tested in Section 5, where some examples of pattern-forming reactions on the torus and sphere are also presented. Finally, the results are discussed in Section 6.
2. Turing Mechanism for Pattern Formation
The mechanism proposed by Turing starts when a transport component is added to a system of two interacting chemical species in a steady state. The diffusion helps to destabilize the system from a spatially homogeneous state to a nonequilibrium steady state with distinguishable structures. Generally, patterns formed in two-dimensional domains with simple boundaries are studied. The particular structure of the spatial patterns depends on the parameters of the kinetics and the diffusion coefficients of each species that must be different.
Let the reaction-diffusion PDE system be as follows: that describes the evolution and spatial behavior of the concentrations and , of the two chemical species. Here, is the Laplace operator in some bounded domain , and are the corresponding diffusion coefficients, and the functions and give the local reaction kinetics. Usually, concentrations and are assumed to be and in time and space, respectively. Moreover, for the system to be well-posed, the kinetics given by and must be sufficiently well-behaved since these functions typically contain nonlinear dissipative terms. As we shall see in the next section, to guarantee that the integrals are finite in the variational formulation, both the concentrations and and the kinetics and must be . For diffusion-driven instability to occur and patterns to form, we must guarantee that the homogeneous steady state is stable under small spatial perturbations when the diffusion coefficients vanish and unstable when they are present. We look for perturbations that decrease exponentially in time and are combinations of spatially oscillatory modes [1, 7]. Furthermore, they must be subject to zero-flux or periodic boundary conditions, which precisely allow self-organizing solutions.
The absence of diffusive terms causes the reactive terms to vanish when evaluated at the spatially homogeneous steady state . Around this solution, it is possible to make a perturbative analysis to consider small instabilities. We can arrange the functions in a vector and, in this case, replace it with a matrix of its derivatives evaluated in the steady state. This imposes certain conditions on the parameters that define the reaction kinetics and that can be written as follows: evaluated in the steady state . These conditions are obtained assuming an exponentially stable behavior of the perturbations. When diffusion is introduced, the system becomes unstable. This can be understood if the perturbations decompose in modes labeled by the wavenumber. These are nothing but the eigenvalues of the corresponding Laplacian operator. The corresponding stability analysis now involves the diffusion coefficients and these wavenumbers and induces the following instability conditions:
The dependence of in wavenumbers, often called the dispersion relation for the linearized system, is as follows: where the subscripts in and indicate partial derivatives. This implies that there is a set of values for which can be positive; this is the so-called range of unstable modes, which depends on the reaction and diffusion constants, and is as follows: , where
In the case of curved surfaces embedded in , the eigenvalues of the tangent Laplacian change in each case. For instance, for the sphere, the eigenvalue equation can be solved in terms of spherical harmonics and the eigenvalues would be ; in such a case, the dispersion relation considered is . For the subsequent numerical implementation of the reaction kinetic models, it is important to choose constants whose values are in this range.
3. Reaction-Diffusion System on a Curved Surface
The goal of this article is to analyze the numerical solution of the reaction-diffusion processes that form the so-called Turing patterns but that takes place on a curved surface embedded in , such as a torus or a sphere. To do this, we first consider the general problem of the system of reaction-diffusion equations on a surface that can be stated as follows: where (i) and represent the local distribution of the two constituents at time on the surface that diffuse at different rates and react according to the nonlinear functions and (ii) and are the diffusion coefficients of the and components, respectively(iii) is the tangential gradient on . We will choose to be the sphere or the torus. Let us realize that the corresponding operator becomes the so-called Laplace-Beltrami operator of the surface(iv) and are the reaction kinetic functions. There are several different models for and that generate different patterns depending on the system under study
Remark 1. Thanks to Turing, we now that, under certain conditions, the system tends to a linearly stable uniform steady state, if there is no diffusion; thus, if , the diffusion-driven instability can cause spatially inhomogeneous patterns to appear in the system [1, 2].
3.1. Variational Formulation of the Reaction-Diffusion System
We will solve problem (6) using the methodology of finite elements; therefore, we need to consider first its variational formulation (see [22, 34]), which is as follows.
Find and such that where (i) ; (ii) is the tangential gradient on ; (iii) is the infinitesimal surface measure; (iv) ; and (v) and are the diffusion coefficients.
Remark 2. The functions and that satisfy the weak formulation (7) are called weak solutions of the state equations (6). Since we use the finite element method to numerically solve system (6), the variational formulation is quite useful since it requires solving an integral problem instead of a differential problem.
Remark 3. The elliptic operators associated with system (7) is clearly (-1) and (-1) times the Laplace-Beltrami operator for the distributions and , respectively.
4. Discretization of the Reaction-Diffusion System
4.1. Time Discretization of the Reaction-Diffusion System (7)
Let us define the time discretization step as , where is the final time and a positive integer. Next, we approximate system (7) by
for , where and . Systems (9) are well-posed elliptic problems; they are not associated with any boundary conditions since is a surface without boundary.
Remark 4. For time discretization of the parabolic problems in (7), we have used the backward Euler scheme, as in [24, 26, 27], getting thus (9). Although this implicit scheme is only fist-order accurate, it is robust and stiff A-stable and preserves the maximum principle if combined with appropriate space approximations.
4.2. Full Discretization of the Reaction-Diffusion System (7)
Now, we will consider the space discretizations of the surface of the torus and the sphere. For the torus, we consider a well-known parametrization to discretize the space. However, for the case of the sphere, we discretize directly to avoid the singularities due to the usual parameterization.
4.2.1. Full Discretization of the Reaction-Diffusion System on the Surface of a Torus
For the space discretization of the torus surface, denoted by , we proceed as in [23, 24] using the following parameterization:
to map over the square of the plane (see Figure 1). We should note that it is necessary to consider periodic boundary conditions to take into account the fact that has no boundary.
The problem being formulated on a planar domain can be approximated using finite element methods discussed in [21]. We first consider a finite element triangulation of with the following properties: is a finite collection of triangles in , with denoting the length of the largest edge of the triangulation , and . Additionally, if , with , we have or and have only one vertex in common or one full edge in common.
Using the parameterization (10), we obtain and we approximate the space of doubly periodic functions by the set where is the space of the two variable first degree polynomials. The time-discrete parabolic problem (9) for a torus can be reformulated as follows: for , where and are an approximation of and in , respectively.
Remark 5. The above discrete linear elliptic problems (14) are associated with the same matrix, differing only by their right-hand sides. This matrix is symmetric positive definite, and sparse, and we solve the associated linear systems by a sparse Cholesky solver. On the other hand, we recommend using a small time step , implying that the matrix associated with the backward Euler scheme is not too badly conditioned, allowing thus the solution of these discrete elliptic problems.
4.2.2. Full Discretization of the Reaction-Diffusion System on the Surface of a Sphere
For the reaction-diffusion processes on the surface of a sphere in , we will consider Cartesian coordinates to represent this surface, to avoid the singularities that the standard parameterization usually introduces. We denote by the surface given by , where is the sphere radius. Now, to obtain the space discretization to numerically solve the elliptic problems in (9), we approximate the sphere by a polyhedral surface and proceed as in [25]. However, to calculate the surface gradient , we will proceed in a slightly different way, e.g., [27]. We approximate by the polyhedral surface , as shown in Figure 2. Let us assume that the elements of are triangular facets, and we denote by the set they form, with the maximum diameter of the elements . From , we approximate through
Unlike [26], we consider an isoparametric parameterization of the surface (the same idea can be used for arbitrary parameterizations), i.e., each triangular facet is described by the where are the coordinates used to define the usual reference element , as in Figure 3; , with , are the physical coordinates; and are the corresponding reference shape functions. Thus, we can discretizate the problems in (9); for instance, the discrete version of the first equation of (9), corresponding to distribution , can be written as follows:
Let and be approximations of and , respectively. Then, for ,
The second equation of (9), corresponding to distribution , is approximated in a similar way. On the other hand, the tangential gradient can be computed by projecting the usual gradient on the surface of , that is, where , is the identity matrix in , and is the orthogonal projection matrix to the normal direction of . This quantity is well defined, since is a constant vector for any and is the (unique) unitary vector normal to .
Let us simplify the notation used in (17), so that the discrete version of system (9) is written as follows: for ,
Remark 6. Note that Remark 5 is also true for the matrix associated to system (19).
To solve numerically the integrals encountered in (14) and (19), we found different numerical methods, such as the trapezoid, Simpson, midpoint, and Monte Carlo. Here, we have used the trapezoidal rule on each triangle and , taking advantage of
5. Numerical Examples
5.1. Numerical Examples on a Torus Embedded on
Next, we present some numerical examples with different functions for the reaction kinetics. In the first case, we consider a torus with minor radius and major radius . For the parameterization in the time, we use . Also, a regular uniform mesh on is considered for the finite element discretization, with (see Figure 1).
Example 7. We first consider the FitzHugh-Nagumo model, see [12], given by where and are the diffusion coefficients of and , respectively.
The FitzHugh-Nagumo kinetic model was first proposed to study the propagation of electrical signals in neurons [6]. In that case, and are interpreted as electric potentials, the fixed value corresponds to the current pulse , and the other parameters are scale constants. However, since it is a simple dynamical system with a variety of structures, the system has been applied to model other systems [12].
For the following examples, we fix the parameters at , , , and . Figure 4 shows the initial condition and the approximation of the final state , with .
(a)
(b)
In Figure 5, we show a different initial condition and the approximations of the state , with and .
(a)
(b)
(c)
Example 8. Now, we choose the reaction-diffusion system proposed by Barrio et al. [7], where and are the diffusion coefficients of and , respectively.
The coefficients show a conservation relation between the chemical products [7, 35]. The principal interaction parameters are and as each favors stripes and dots, respectively. The parameters , , and are related to production and depletion of chemicals. This model has been used mainly to study pattern selection in reaction-diffusion models. Also, it has been used to study transitions between fish patterns that go from stripes to spots.
For the numerical results of these examples, we fix the parameters at , , , , , , and . The major and minor radii that we consider are and , respectively. Figure 6 shows the initial conditions and the approximations of the state , with .
(a)
(b)
Figure 7 shows another initial condition and the approximations of the state , with , 200, and 300.
(a)
(b)
For the same parameters and conditions of the reaction-diffusion system of system (22), we now explore how the formation of patterns changes for different surfaces by changing the radii of the torus; this will be related to the curvature of the surface. Let us recall that the Gaussian curvature of the torus is given as a function of the radii and the angle as follows:
Hence, for each set of fixed radii, the torus has a negative curvature on the inside and positive on the outside. However, if the radii vary, we would have a landscape of possible surfaces. In Figure 8, two different views of each experiment are shown for the final distributions of the reaction-diffusion system (22) with different radii, corresponding to different surfaces in the curvature landscape.
(a)
(b)
(c)
We observe that the patterns that appear earlier if the surface have a more pronounced curvature at as seen in Figure 9. Although we cannot conclude a general behavior, our examples serve as evidence indicating that the appearance of patterns will depend on the curvature and other geometric factors of the surface on which the diffusion process is taking place, as reported in the literature [8, 12, 35].
(a)
(b)
5.2. Numerical Examples on a Sphere in
Example 9. First, we will numerically study the effects of a reaction-diffusion system on spheres of different radii. We consider the Schnakenberg model [3, 8], given by where and are the diffusion coefficients of species and , respectively. The numerical results show spot patterns for different cases.
With a minimal number of reactions and reactants, the Schnakenberg model simulates the chemical reactions that occur in some biological systems, such as population cycles and metabolic regulation. Therefore, it has been extensively studied. It is similar to the Gierer-Meinhardt model used to study the activator-inhibitor relation between chemicals. The model describes the evolution of concentrations of chemicals and , produced at rates and , respectively.
For the first case of this example, we fix the parameters at , , , , , and ( is the radius of sphere). For the discretization, we used 0.0051 and . In Figure 10, we show the final distribution for , 300, 400, and 600.
(a)
(b)
(c)
(d)
For the next example, we considered , , 0.1, , , and . The results are shown in Figure 11. Also in this figure, we are able to see that three different views correspond to different times of distribution.
Example 10. Next, we will consider the Gierer-Meinhardt model with saturation, given as follows: where and are the diffusion coefficients of the morphogens and , respectively.
In the Gierer and Meinhardt model, species acts as an activator of both species and as an inhibitor of activator sources. In agreement with the Turing conditions, this model assumes that diffuses faster than , and and are again the production rates [4]. Saturation in this model is introduced by modifying the quadratic production term in , by , where is a sufficiently large constant, so that when tends to , it is said to be saturated. In this case, the production and consumption constants with saturation are and , respectively; is a decrease constant and a constant source. This model has been used to study the patterns that appear in the ladybeetles, using certain tuned values for the parameters [36]. One of the features of that study was that they worked with finite differences, and therefore, only a very specific section of the sphere could be studied. Since the analysis used here presents no problems at the poles, we can extend the domain of the solution to the entire sphere.
For the following example, we fix the parameters at , , , , and . For the discretization, we used and . Figure 12 shows the initial conditions and and the approximations of the final state and for the case , , and .
(a)
(b)
(c)
(d)
Figure 13 shows the final distributions and for the case , , and ; the initial distribution is the same initial distribution as the above case.
(a)
(b)
Stripe patterns are shown in Figure 14; to obtain these patterns, we use the following parameters , , and . Figures 14(a) and 14(b) show the initial distributions, which are different to the above cases.
(a)
(b)
(c)
(d)
A variant for this example given by
The numerical results are shown in Figure 15; in order to obtained these results, we fixed the parameters at , , , and . Also, the initial distribution that we considered is the same as in the cases for Figures 14(a) and 14(b).
(a)
(b)
6. Discussion
In this manuscript, the numerical solution of different systems of reaction-diffusion equations that use the Turing mechanism for pattern formation has been studied. In biological applications, the systems on two-dimensional surfaces embedded in are of particular interest, so here, we study these processes on the surface of the torus and the sphere. To solve this kind of problem, we propose to use the finite element approximation with first-order polynomials.
To discretize the temporal part, we have used the backward Euler scheme, resulting in well-posed elliptic problems. For spatial discretization, we use two different approaches. For the torus, we used the standard parameterization, and for the sphere, we used direct discretization to avoid the singularities that arise in the standard parameterization. To solve the elliptic problems, we have used finite elements with first-order polynomials. For the torus, we have considered periodic boundary conditions due to the parameterization used, while on the sphere, we have not considered any boundary since we are solving the system directly on the surface.
For both surfaces, it was possible to obtain the concentration patterns after at least 1000 iterations, for several different reactions. In addition to the known influence of the values of the reaction parameters on the formation of the patterns, here, an influence of the curvature of the surfaces through the variation of their radii is also appreciated, specially for the reaction kinetics of the BVAM model, on tori of different radii as seen in Figure 8. In particular, for the last studied reaction kinetics on the sphere, the Gierer-Meinhardt model with saturation, solutions can be obtained on the entire sphere since this method allows evading the singularity at the pole, so that patterns can be seen in the whole sphere, instead of just over a hemisphere as in [36].
Given that the variational methods studied for nonlinear dissipative models give exact solutions [30, 33], it would be interesting to address further the reaction-diffusion systems for some kinetics of interest, such as the saturation one, with these methods.
Data Availability
All results have been obtained by conducting the numerical procedure, and the ideas can be shared with the researchers.
Conflicts of Interest
The authors declare no conflict of interest.
Acknowledgments
The authors are thankful for the support from the Convocatoria 2020 para el Fortalecimiento de Cuerpos Académicos of the Division of Natural Sciences and Engineering (DCNI) of UAM Cuajimalpa, which sponsored this work.