Full Finite Element Scheme for Reaction-Diffusion Systems on Embedded Curved Surfaces in R3

The purpose of this article is to study numerically the Turing diffusion-driven instability mechanism for pattern formation on curved surfaces embedded in R3, 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.


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][11][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 ℝ 3 . 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][22][23][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 ℝ 3 ; 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 reactiondiffusion 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.

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 uðx, tÞ and vðx, tÞ, of the two chemical species.
Here, ∇ 2 is the Laplace operator in some bounded domain Ω ∈ ℝ n , D u and D v are the corresponding diffusion coefficients, and the functions f ðu, vÞ and gðu, vÞ give the local reaction kinetics. Usually, concentrations u and v are assumed to be C 1 ð½0, TÞ and C 2 ðΩÞ in time and space, respectively. Moreover, for the system to be well-posed, the kinetics given by f and g 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 u and v and the kinetics f and g must be L 2 ðΩ, ð0, TÞÞ. 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 ðu s , v s Þ. 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: 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 f and g indicate partial derivatives. This implies that there is a set of values k 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: In the case of curved surfaces embedded in ℝ 3 , 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 k 2 = nðn + 1Þ; in such a case, the dispersion relation considered is λðnðn + 1ÞÞ. For the subsequent numerical implementation of the reaction kinetic models, it is important to choose constants whose values are in this range.

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 ℝ 3 , such as a torus or a sphere. To do this, we first consider the general problem of the system of reactiondiffusion equations on a surface that can be stated as follows: where (i) uðx, tÞ and vðx, tÞ represent the local distribution of the two constituents at time t on the surface Σ that diffuse at different rates and react according to the nonlinear functions f and g (ii) D u and D v are the diffusion coefficients of the u and v 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 ∇ 2 Σ becomes the so-called Laplace-Beltrami operator of the surface (iv) f ðu, vÞ and gðu, vÞ are the reaction kinetic functions.
There are several different models for f and g 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 D u ≠ D v , the diffusiondriven instability can cause spatially inhomogeneous patterns to appear in the system [1,2].

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. 3

Advances in Mathematical Physics
Find u and v ∈ H 1 ðΣÞ such that Remark 2. The functions u and v 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)D u and (-1)D v times the Laplace-Beltrami operator for the distributions u and v, respectively. (7). Let us define the time discretization step Δt as Δt = T/N, where T is the final time and N a positive integer. Next, we approximate system (7) by

Time Discretization of the Reaction-Diffusion System
where u n = uðx, nΔtÞ and v n = vðx, nΔtÞ. 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. (7). Now, we will consider the space discretizations of the surface of the torus and the sphere. For the torus, we consider a wellknown parametrization to discretize the space. However, for the case of the sphere, we discretize directly to avoid the singularities due to the usual parameterization.

Full Discretization of the Reaction-Diffusion
System on the Surface of a Torus. For the space discretization of the torus surface, denoted by Σ T , we proceed as in [23,24] using the following parameterization: Þcos ϕ, to map Σ T over the square b Ω = ð0, 2πÞ × ð0, 2πÞ of the plane ðϕ, θÞ (see Figure 1). We should note that it is necessary Advances in Mathematical Physics to consider periodic boundary conditions to take into account the fact that Σ T 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 τ h of b Ω with the following properties: τ h is a finite collection of tri- 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 ℙ 1 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 n = 1, ⋯, N, where u 0h and v 0h are an approximation of u 0 and v 0 in V h , respectively.
Remark 5. The above discrete linear elliptic problems (14) are associated with the same matrix, differing only by their righthand 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 Δt, implying that the matrix associated with the backward Euler scheme is not too badly conditioned, allowing thus the solution of these discrete elliptic problems.

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 ℝ 3 , 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 x 2 + y 2 + z 2 = R 2 , where R 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

Advances in Mathematical Physics
Σ by the polyhedral surface Σ h , as shown in Figure 2. Let us assume that the elements of Σ h are triangular facets, and we denote by T h the set they form, with h the maximum diameter of the elements T ∈ T h . From T h , we approximate H 1 ðΣÞ 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 T ∈ T h is described by the where ðε, ηÞ are the coordinates used to define the usual reference elementT ∈ ℝ 2 , as in Figure 3; , are the physical coordinates; and φ i ðε, ηÞ 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 u, can be written as follows: Let u 0,h ∈ V h and v 0,h ∈ V h be approximations of u 0 and v 0 , respectively. Then, for n = 1, ⋯, N, The second equation of (9), corresponding to distribution v, is approximated in a similar way. On the other hand, the tangential gradient ∇ T u n h can be computed by projecting the usual gradient ∇u n h on the surface of T, that is, where u n h ∈ V h , I is the identity matrix in ℝ 3 , and ℙ T = n T n t T is the orthogonal projection matrix to the normal direction of T. This quantity is well defined, since ∇u n h is a constant vector for any T ∈ T h and n T is the (unique) unitary vector normal to T. Let us simplify the notation used in (17), so that the discrete version of system (9) Advances in Mathematical Physics 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 T ∈ τ h and T ∈ T h , taking advantage of ð

Numerical
Examples on a Torus Embedded on ℝ 3 . Next, we present some numerical examples with different functions for the reaction kinetics. In the first case, we consider a torus with minor radius r and major radius R. For the parameterization in the time, we use Δt = 0:01. Also, a regular uniform mesh on b Ω = ð0, 2πÞ × ð0, 2πÞ is considered for the finite element discretization, with Δϕ = Δθ = 2π/100 (see Figure 1).

Example 7.
We first consider the FitzHugh-Nagumo model, see [12], given by where D u and D v are the diffusion coefficients of u and v, respectively.
The FitzHugh-Nagumo kinetic model was first proposed to study the propagation of electrical signals in neurons [6]. In that case, u and v are interpreted as electric potentials, the fixed value −0:6 corresponds to the current pulse I, 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 D u = 1, D v = 1:75, R = 20, and r = R/3. Figure 4 shows the initial condition u 0 and the approximation of the final state uðx, TÞ, with T = 100.
In Figure 5, we show a different initial condition u 0 and the approximations of the state uðx, TÞ, with T = 100 and T = 300.
Example 8. Now, we choose the reaction-diffusion system proposed by Barrio et al. [7], where D u and D v are the diffusion coefficients of u and v, respectively.
The coefficients show a conservation relation between the chemical products [7,35]. The principal interaction parameters are r 1 and r 2 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 same parameters and conditions of the reactiondiffusion 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 u of the reactiondiffusion system (22) with different radii, corresponding to different surfaces in the curvature landscape.
We observe that the patterns that appear earlier if the surface have a more pronounced curvature at θ = π as seen in

Numerical Examples on a Sphere in ℝ 3
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 D u and D v are the diffusion coefficients of species u and v, 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 u and v, produced at rates a and b, respectively.
For the first case of this example, we fix the parameters at D u = 1/R, D v = 10/R, a = 0:01, b = 1:2, T = 600, and R = 40 (R is the radius of sphere). For the discretization, we used h = 0.0051 and δt = 1/100. In Figure 10, we show the final distribution u for T = 100, 300, 400, and 600.
For the next example, we considered D u = 1/20, D v = 1, a =0.1, b = 0:9, T = 300, and R = 20. 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 D u and D v are the diffusion coefficients of the morphogens u and v, respectively.
In the Gierer and Meinhardt model, species u acts as an activator of both species and v as an inhibitor of activator sources. In agreement with the Turing conditions, this model assumes that v diffuses faster than u, and a and b are again the production rates [4]. Saturation i in this model is introduced by modifying the quadratic production term in u 2 v, by u 2 v/ð1 + ku 2 Þ , where k is a sufficiently large constant, so    11 Advances in Mathematical Physics that when u tends to k, it is said to be saturated. In this case, the production and consumption constants with saturation are ρ u and ρ v , respectively; μ u is a decrease constant and σ v 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.

Advances in Mathematical Physics
Initial solution u 0 ,

Advances in Mathematical Physics
Stripe patterns are shown in Figure 14; to obtain these patterns, we use the following parameters D u = 0:000028, D v = 0:000168, and k = 0:35. Figures 14(a) and 14(b) show the initial distributions, which are different to the above cases.
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 D u = 0:000026, D v = 0:000182, k = 0:45, and σ u = 0:0019. Also, the initial distribution that we considered is the same as in the cases for Figures 14(a) and 14(b).

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 ℝ 3 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 N 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 N 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.