# Full Finite Element Scheme for Reaction-Diffusion Systems on Embedded Curved Surfaces inℝ3

Authors: D. Assaely León-Velasco; Guillermo Chacón-Acosta
Journal: Advances in Mathematical Physics
 (2021)

---

## Abstract

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


---

## Body

## 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ℝ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–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 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:(1)∂u∂t=Du∇2u+fu,v,∂v∂t=Dv∇2v+gu,v,that describes the evolution and spatial behavior of the concentrations ux,t and vx,t, of the two chemical species. Here, ∇2 is the Laplace operator in some bounded domain Ω∈ℝn, Du and Dv are the corresponding diffusion coefficients, and the functions fu,v and gu,v give the local reaction kinetics. Usually, concentrations u and v are assumed to be C10,T and C2Ω 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 L2Ω,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 stateus,vs. 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:
(2)∂f∂u+∂g∂v<0,∂f∂u∂g∂v−∂f∂v∂g∂u>0,evaluated in the steady state us,vs. 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:
(3)Dv∂f∂u+Du∂g∂v>0,Dv∂f∂u+Du∂g∂v>2DuDv∂f∂u∂g∂v−∂f∂v∂g∂u.The dependence ofλ in wavenumbers, often called the dispersion relation for the linearized system, is as follows:
(4)λk2=12fu+gv−k2Du+Dv±k2Du+Dv−fu+gv2−4DuDvk4−k2Dvfu+Dugv+fugv−fvgu,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: k−2,k+2, where
(5)k±2=12DuDvDv∂f∂u+Du∂g∂v±Dv∂f∂u+Du∂g∂v2−4DuDv∂f∂u∂g∂v−∂f∂v∂g∂u.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 k2=nn+1; in such a case, the dispersion relation considered is λnn+1. 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ℝ3, 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:
(6)∂u∂t=Du∇Σ2u+fu,v,onΣ,∂v∂t=Dv∇Σ2v+gu,v,onΣ,initialconditions:ux,0=u0x,vx,0=v0x,where
(i)
ux,t and vx,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)
Du and Dv 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)
fu,v and gu,v are the reaction kinetic functions. There are several different models for f and g that generate different patterns depending on the system under studyRemark 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, ifDu≠Dv, 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.Findu and v∈H1Σ such that
(7)∫Σ∂u∂tzdΣ+Du∫Σ∇Σu·∇ΣzdΣ=∫Σfu,vzdΣ,∀z∈H1Σ,a.e.on0,T,ux,0=u0x,∫Σ∂v∂tzdΣ+Dv∫Σ∇Σv·∇ΣzdΣ=∫Σgu,vzdΣ,∀z∈H1Σ,a.e.on0,T,vx,0=v0x,where (i) H1Σ=z∣z∈L2Σ,∫Σ∇Σz2dΣ<+∞; (ii) ∇Σ is the tangential gradient on Σ; (iii) dΣ is the infinitesimal surface measure; (iv) u0x,v0x∈L2Σ; and (v) Du and Dv are the diffusion coefficients.Remark 2.
The functionsu 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)Du and (-1)Dv times the Laplace-Beltrami operator for the distributions u and v, respectively.


## 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.Findu and v∈H1Σ such that
(7)∫Σ∂u∂tzdΣ+Du∫Σ∇Σu·∇ΣzdΣ=∫Σfu,vzdΣ,∀z∈H1Σ,a.e.on0,T,ux,0=u0x,∫Σ∂v∂tzdΣ+Dv∫Σ∇Σv·∇ΣzdΣ=∫Σgu,vzdΣ,∀z∈H1Σ,a.e.on0,T,vx,0=v0x,where (i) H1Σ=z∣z∈L2Σ,∫Σ∇Σz2dΣ<+∞; (ii) ∇Σ is the tangential gradient on Σ; (iii) dΣ is the infinitesimal surface measure; (iv) u0x,v0x∈L2Σ; and (v) Du and Dv are the diffusion coefficients.Remark 2.
The functionsu 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)Du and (-1)Dv times the Laplace-Beltrami operator for the distributions u and v, respectively.


## 4. Discretization of the Reaction-Diffusion System

### 4.1. Time Discretization of the Reaction-Diffusion System (7)

Let us define thetime discretization stepΔt as Δt=T/N, where T is the final time and N a positive integer. Next, we approximate system (7) by
(8)u0=u0,v0=v0,forn=1,⋯,N,
(9)un,vn∈H1Σ,∫Σun−un−1ΔtzdΣ+Du∫Σ∇Σun·∇ΣzdΣ=∫Σfun,vnzdΣ,∀z∈H1Σ,∫Σvn−vn−1ΔtzdΣ+Dv∫Σ∇Σvn·∇ΣzdΣ=∫Σgun,vnzdΣ,∀z∈H1Σ,where un=ux,nΔt and vn=vx,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.


### 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ΣT, we proceed as in [23, 24] using the following parameterization:
(10)x1=R+rcosθcosϕ,x2=R+rcosθsinϕ,x3=rsinθ,to mapΣT over the square Ω^=0,2π×0,2π 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 ΣT has no boundary.Figure 1
(a) A surface of a torus. (b) Triangulation ofΩ^¯.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 Ω^ with the following properties: τh is a finite collection of triangles T in Ω^¯, with h denoting the length of the largest edge of the triangulation τh, and Ω^¯=⋃T∈τhT. Additionally, if T,T′∈τh, with T≠T′, we have T∩T′=∅ or T and T′ have only one vertex in common or one full edge in common.Using the parameterization (10), we obtain
(11)dΣT=rR+rcosθdϕdθ,∇ΣT=1R+rcosθ∂y∂ϕ,1r∂y∂θ,and we approximate the space of doubly periodic functions
(12)Hp1Ω^=z∣z∈H1Ω^,zϕ,0=zϕ,2πa.e.on0,2π,z0,θ=z2π,θ,a.e.on0,2π,by the set
(13)Vh=zz∈C0Ω^¯,zT∈ℙ1,∀T∈τhz0,θ=z2π,θ,zϕ,0=zϕ,2π,∀ϕ,θ∈0,2π×0,2π,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,
(14)uhn,vhn∈Vh,r∫Ω^uhn−uhn−1ΔtzR+rcosθdϕdθ+Du∫Ω^rR+rcosθ∂uhn∂ϕ∂z∂ϕ+R+rcosθr∂uhn∂θ∂z∂θdϕdθ=r∫Ω^fuhn,vhnzR+rcosθdϕdθ,∀z∈Vh,r∫Ω^vhn−vhn−1ΔtzR+rcosθdϕdθ+Dv∫Ω^rR+rcosθ∂vhn∂ϕ∂z∂ϕ+R+rcosθr∂vhn∂θ∂z∂θdϕdθ=r∫Ω^guhn,vhnzR+rcosθdϕdθ,∀z∈Vh,where u0h and v0h are an approximation of u0 and v0 in Vh, 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 Δ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.


#### 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ℝ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 x2+y2+z2=R2, 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 Σ 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 Th the set they form, with h the maximum diameter of the elements T∈Th. From Th, we approximate H1Σ through
(15)Vh=zz∈H1Σh,zT∈ℙ1,∀T∈Th.Figure 2
A polyhedral approximation of a sphere of radiusR=10.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∈Th is described by the
(16)xhε,η=∑i=1nxiφiε,η,where ε,η are the coordinates used to define the usual reference element T^∈ℝ2, as in Figure 3; xi=x1i,x2i,x3i, with i=1,2,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:Figure 3
A usual reference elementT^∈ℝ2.Letu0,h∈Vh and v0,h∈Vh be approximations of u0 and v0, respectively. Then, for n=1,⋯,N,
(17)uhn,vhn∈Vh,∑T∈Th∫Tuhn−uhn−1ΔtzdT+Du∑T∈Th∫T∇Tuhn·∇TzdT=∑T∈Th∫Tfuhn,vhnzdT,∀z∈Vh.The second equation of (9), corresponding to distribution v, is approximated in a similar way. On the other hand, the tangential gradient ∇Tuhn can be computed by projecting the usual gradient ∇uhn on the surface of T, that is,
(18)∇Tuhn=I−ℙT∇uhn,where uhn∈Vh, I is the identity matrix in ℝ3, and ℙT=nTnTt is the orthogonal projection matrix to the normal direction of T. This quantity is well defined, since ∇uhn is a constant vector for any T∈Th and nT is the (unique) unitary vector normal to T.Let us simplify the notation used in (17), so that the discrete version of system (9) is written as follows: for n=1,⋯,N,
(19)uhn,vhn∈Vh,∫Σhuhn−uhn−1ΔtzdΣh+Du∫Σh∇Σhuhn·∇ΣhzdΣh=∫Σhfuhn,vhnzdΣh,∀z∈Vh,∫Σhvhn−vhn−1ΔtzdΣh+Dv∫Σh∇Σhvhn·∇ΣhzdΣh=∫Σhguhn,vhnzdΣh,∀z∈Vh.Remark 6.
Note that Remark5 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∈Th, taking advantage of
(20)∫Ω^=∑T∈τh∫T,∫Σh=∑T∈Th∫T.


## 4.1. Time Discretization of the Reaction-Diffusion System (7)

Let us define thetime discretization stepΔt as Δt=T/N, where T is the final time and N a positive integer. Next, we approximate system (7) by
(8)u0=u0,v0=v0,forn=1,⋯,N,
(9)un,vn∈H1Σ,∫Σun−un−1ΔtzdΣ+Du∫Σ∇Σun·∇ΣzdΣ=∫Σfun,vnzdΣ,∀z∈H1Σ,∫Σvn−vn−1ΔtzdΣ+Dv∫Σ∇Σvn·∇ΣzdΣ=∫Σgun,vnzdΣ,∀z∈H1Σ,where un=ux,nΔt and vn=vx,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.


## 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ΣT, we proceed as in [23, 24] using the following parameterization:
(10)x1=R+rcosθcosϕ,x2=R+rcosθsinϕ,x3=rsinθ,to mapΣT over the square Ω^=0,2π×0,2π 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 ΣT has no boundary.Figure 1
(a) A surface of a torus. (b) Triangulation ofΩ^¯.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 Ω^ with the following properties: τh is a finite collection of triangles T in Ω^¯, with h denoting the length of the largest edge of the triangulation τh, and Ω^¯=⋃T∈τhT. Additionally, if T,T′∈τh, with T≠T′, we have T∩T′=∅ or T and T′ have only one vertex in common or one full edge in common.Using the parameterization (10), we obtain
(11)dΣT=rR+rcosθdϕdθ,∇ΣT=1R+rcosθ∂y∂ϕ,1r∂y∂θ,and we approximate the space of doubly periodic functions
(12)Hp1Ω^=z∣z∈H1Ω^,zϕ,0=zϕ,2πa.e.on0,2π,z0,θ=z2π,θ,a.e.on0,2π,by the set
(13)Vh=zz∈C0Ω^¯,zT∈ℙ1,∀T∈τhz0,θ=z2π,θ,zϕ,0=zϕ,2π,∀ϕ,θ∈0,2π×0,2π,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,
(14)uhn,vhn∈Vh,r∫Ω^uhn−uhn−1ΔtzR+rcosθdϕdθ+Du∫Ω^rR+rcosθ∂uhn∂ϕ∂z∂ϕ+R+rcosθr∂uhn∂θ∂z∂θdϕdθ=r∫Ω^fuhn,vhnzR+rcosθdϕdθ,∀z∈Vh,r∫Ω^vhn−vhn−1ΔtzR+rcosθdϕdθ+Dv∫Ω^rR+rcosθ∂vhn∂ϕ∂z∂ϕ+R+rcosθr∂vhn∂θ∂z∂θdϕdθ=r∫Ω^guhn,vhnzR+rcosθdϕdθ,∀z∈Vh,where u0h and v0h are an approximation of u0 and v0 in Vh, 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 Δ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.


### 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ℝ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 x2+y2+z2=R2, 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 Σ 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 Th the set they form, with h the maximum diameter of the elements T∈Th. From Th, we approximate H1Σ through
(15)Vh=zz∈H1Σh,zT∈ℙ1,∀T∈Th.Figure 2
A polyhedral approximation of a sphere of radiusR=10.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∈Th is described by the
(16)xhε,η=∑i=1nxiφiε,η,where ε,η are the coordinates used to define the usual reference element T^∈ℝ2, as in Figure 3; xi=x1i,x2i,x3i, with i=1,2,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:Figure 3
A usual reference elementT^∈ℝ2.Letu0,h∈Vh and v0,h∈Vh be approximations of u0 and v0, respectively. Then, for n=1,⋯,N,
(17)uhn,vhn∈Vh,∑T∈Th∫Tuhn−uhn−1ΔtzdT+Du∑T∈Th∫T∇Tuhn·∇TzdT=∑T∈Th∫Tfuhn,vhnzdT,∀z∈Vh.The second equation of (9), corresponding to distribution v, is approximated in a similar way. On the other hand, the tangential gradient ∇Tuhn can be computed by projecting the usual gradient ∇uhn on the surface of T, that is,
(18)∇Tuhn=I−ℙT∇uhn,where uhn∈Vh, I is the identity matrix in ℝ3, and ℙT=nTnTt is the orthogonal projection matrix to the normal direction of T. This quantity is well defined, since ∇uhn is a constant vector for any T∈Th and nT is the (unique) unitary vector normal to T.Let us simplify the notation used in (17), so that the discrete version of system (9) is written as follows: for n=1,⋯,N,
(19)uhn,vhn∈Vh,∫Σhuhn−uhn−1ΔtzdΣh+Du∫Σh∇Σhuhn·∇ΣhzdΣh=∫Σhfuhn,vhnzdΣh,∀z∈Vh,∫Σhvhn−vhn−1ΔtzdΣh+Dv∫Σh∇Σhvhn·∇ΣhzdΣh=∫Σhguhn,vhnzdΣh,∀z∈Vh.Remark 6.
Note that Remark5 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∈Th, taking advantage of
(20)∫Ω^=∑T∈τh∫T,∫Σh=∑T∈Th∫T.


## 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ΣT, we proceed as in [23, 24] using the following parameterization:
(10)x1=R+rcosθcosϕ,x2=R+rcosθsinϕ,x3=rsinθ,to mapΣT over the square Ω^=0,2π×0,2π 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 ΣT has no boundary.Figure 1
(a) A surface of a torus. (b) Triangulation ofΩ^¯.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 Ω^ with the following properties: τh is a finite collection of triangles T in Ω^¯, with h denoting the length of the largest edge of the triangulation τh, and Ω^¯=⋃T∈τhT. Additionally, if T,T′∈τh, with T≠T′, we have T∩T′=∅ or T and T′ have only one vertex in common or one full edge in common.Using the parameterization (10), we obtain
(11)dΣT=rR+rcosθdϕdθ,∇ΣT=1R+rcosθ∂y∂ϕ,1r∂y∂θ,and we approximate the space of doubly periodic functions
(12)Hp1Ω^=z∣z∈H1Ω^,zϕ,0=zϕ,2πa.e.on0,2π,z0,θ=z2π,θ,a.e.on0,2π,by the set
(13)Vh=zz∈C0Ω^¯,zT∈ℙ1,∀T∈τhz0,θ=z2π,θ,zϕ,0=zϕ,2π,∀ϕ,θ∈0,2π×0,2π,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,
(14)uhn,vhn∈Vh,r∫Ω^uhn−uhn−1ΔtzR+rcosθdϕdθ+Du∫Ω^rR+rcosθ∂uhn∂ϕ∂z∂ϕ+R+rcosθr∂uhn∂θ∂z∂θdϕdθ=r∫Ω^fuhn,vhnzR+rcosθdϕdθ,∀z∈Vh,r∫Ω^vhn−vhn−1ΔtzR+rcosθdϕdθ+Dv∫Ω^rR+rcosθ∂vhn∂ϕ∂z∂ϕ+R+rcosθr∂vhn∂θ∂z∂θdϕdθ=r∫Ω^guhn,vhnzR+rcosθdϕdθ,∀z∈Vh,where u0h and v0h are an approximation of u0 and v0 in Vh, 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 Δ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.


## 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ℝ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 x2+y2+z2=R2, 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 Σ 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 Th the set they form, with h the maximum diameter of the elements T∈Th. From Th, we approximate H1Σ through
(15)Vh=zz∈H1Σh,zT∈ℙ1,∀T∈Th.Figure 2
A polyhedral approximation of a sphere of radiusR=10.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∈Th is described by the
(16)xhε,η=∑i=1nxiφiε,η,where ε,η are the coordinates used to define the usual reference element T^∈ℝ2, as in Figure 3; xi=x1i,x2i,x3i, with i=1,2,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:Figure 3
A usual reference elementT^∈ℝ2.Letu0,h∈Vh and v0,h∈Vh be approximations of u0 and v0, respectively. Then, for n=1,⋯,N,
(17)uhn,vhn∈Vh,∑T∈Th∫Tuhn−uhn−1ΔtzdT+Du∑T∈Th∫T∇Tuhn·∇TzdT=∑T∈Th∫Tfuhn,vhnzdT,∀z∈Vh.The second equation of (9), corresponding to distribution v, is approximated in a similar way. On the other hand, the tangential gradient ∇Tuhn can be computed by projecting the usual gradient ∇uhn on the surface of T, that is,
(18)∇Tuhn=I−ℙT∇uhn,where uhn∈Vh, I is the identity matrix in ℝ3, and ℙT=nTnTt is the orthogonal projection matrix to the normal direction of T. This quantity is well defined, since ∇uhn is a constant vector for any T∈Th and nT is the (unique) unitary vector normal to T.Let us simplify the notation used in (17), so that the discrete version of system (9) is written as follows: for n=1,⋯,N,
(19)uhn,vhn∈Vh,∫Σhuhn−uhn−1ΔtzdΣh+Du∫Σh∇Σhuhn·∇ΣhzdΣh=∫Σhfuhn,vhnzdΣh,∀z∈Vh,∫Σhvhn−vhn−1ΔtzdΣh+Dv∫Σh∇Σhvhn·∇ΣhzdΣh=∫Σhguhn,vhnzdΣh,∀z∈Vh.Remark 6.
Note that Remark5 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∈Th, taking advantage of
(20)∫Ω^=∑T∈τh∫T,∫Σh=∑T∈Th∫T.


## 5. Numerical Examples

### 5.1. 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 radiusr and major radius R. For the parameterization in the time, we use Δt=0.01. Also, a regular uniform mesh on Ω^¯=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
(21)∂u∂t=Du∇2u+u−u33+v−0.6,∂v∂t=Dv∇2v−u−0.6+9.99v,where Du and Dv 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 atDu=1, Dv=1.75, R=20, and r=R/3. Figure 4 shows the initial condition u0 and the approximation of the final state ux,T, with T=100.Figure 4
(a) Initial distributionu0. (b) Final distribution u. Case: equation (21), Du=1, Dv=1.75, T=100, R=20, and r=R/3.
(a)(b)In Figure5, we show a different initial condition u0 and the approximations of the state ux,T, with T=100 and T=300.Figure 5
(a) Initial distributionu0. (b) Distribution u at T=100. (c) Final distribution u at T=300. Case: equation (21), Du=1, Dv=1.75, R=20, and r=R/3.
(a)(b)(c)Example 8.
Now, we choose the reaction-diffusion system proposed by Barrio et al. [7],
(22)∂u∂t=Du∇2u+αu1−r1v2+v1−r2u,∂v∂t=Dv∇2v+βv1+αr1βuv+uγ+r2v,where Du and Dv 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 r1 and r2 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 atDu=2.322×10−3, Dv=0.0045, α=0.899, β=−0.91, γ=−α, r1=3.5, and r2=0. The major and minor radii that we consider are R=1 and r=0.3, respectively. Figure 6 shows the initial conditions u0 and the approximations of the state ux,T, with T=200.Figure 6
(a) Initial distributionu0. (b) Final distribution u. Case: equation (22), Du=2.322×10−3, Dv=0.0045, T=200, R=1 and r=0.3.
(a)(b)Figure7 shows another initial condition u0 and the approximations of the state ux,T, with T=100, 200, and 300.Figure 7
(a) Initial distributionu0. (b) Final distributions u. Case: equation (22); Du=2.322×10−3; Dv=0.0045; T=100,200,300; R=1; and r=0.3.
(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:
(23)K=cosθRr+r2cosθ.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 Figure8, two different views of each experiment are shown for the final distributions u of the reaction-diffusion system (22) with different radii, corresponding to different surfaces in the curvature landscape.Figure 8
Final distributionux,T with (a) R=0.75, r=0.45, and T=30; (b) R=0.9, r=0.35, and T=100; (c) R=0.95, r=0.25, and T=200. Case: equation (22), Du=2.322×10−3, Dv=0.0045.
(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].Figure 9
(a) Curvature as a function ofθ for different values of radii. For certain values of R decreasing and r increasing, the minimum of K in π becomes more pronounced. (b) Curvature landscape for different radii at θ=π. The tori of different radii in the right figure are shown as points with the indicated color code.
(a)(b)


### 5.2. 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
(24)∂u∂t=Du∇2u+a−u+u2v,∂v∂t=Dv∇2v+b−u2v,where Du and Dv 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 chemicalsu and v, produced at rates a and b, respectively.For the first case of this example, we fix the parameters atDu=1/R, Dv=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.Figure 10
Final distributionu: (a) T=100, (b) T=300, (c) T=400, and (d) T=600. Case: equation (24) R =40, Du=1/40, Dv=1/4, a=0.01, and b=1.2.
(a)(b)(c)(d)For the next example, we consideredDu=1/20, Dv=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.Figure 11
Distributionux,T. Case: R=20, Du=1/20, Dv=1, a=0.1, and b=0.9.Example 10.
Next, we will consider the Gierer-Meinhardt model with saturation, given as follows:(25)∂u∂t=Du∇2u+ρuu2v1+ku2−μuu,∂v∂t=Dv∇2v−ρvu2v1+ku2+σv,where Du and Dv are the diffusion coefficients of the morphogens u and v, respectively.In the Gierer and Meinhardt model, speciesu 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 u2v, by u2v/1+ku2, where k is a sufficiently large constant, so 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.For the following example, we fix the parameters atρu=0.18, ρv=0.36, μu=0.08, σv=0.1, and T=100. For the discretization, we used h=0.0051 and δt=1/100. Figure 12 shows the initial conditions u0 and v0 and the approximations of the final state ux,100 and vx,100 for the case Du=0.0005, Dv=0.025, and k=0.Figure 12
(a) Initial distributionu0. (b) Initial condition v0. (c) Final distribution u. (d) Final state v. Case: equation (25), Du=0.0005, Dv=0.025, T=100, and k=0.
(a)(b)(c)(d)Figure13 shows the final distributions ux,100 and vx,100 for the case Du=0.0005, Dv=0.025, and k=0.35; the initial distribution is the same initial distribution as the above case.Figure 13
(a) Final distributionu. (b) Final distribution v. Case: equation (25), Du=0.0005, Dv=0.025, T=100, and k=0.35.
(a)(b)Stripe patterns are shown in Figure14; to obtain these patterns, we use the following parameters Du=0.000028, Dv=0.000168, and k=0.35. Figures 14(a) and 14(b) show the initial distributions, which are different to the above cases.Figure 14
(a) Initial distributionu0. (b) Final condition v0. (c) Final distribution u. (d) Final state v. Case: equation (25), Du=0.000028, Dv=0.000168, T=100, and k=0.35.
(a)(b)(c)(d)A variant for this example given by(26)∂u∂t=Du∇2u+ρuu2v1+ku2−μuu,∂v∂t=Dv∇2v−ρvu2v1+ku2+σv+σu,The numerical results are shown in Figure15; in order to obtained these results, we fixed the parameters at Du=0.000026, Dv=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).Figure 15
(a) Final distributionu. (b) Final distribution v. Case: equation (26), Du=0.000026, Dv=0.000182, T=100, and k=0.45.
(a)(b)


## 5.1. 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 radiusr and major radius R. For the parameterization in the time, we use Δt=0.01. Also, a regular uniform mesh on Ω^¯=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
(21)∂u∂t=Du∇2u+u−u33+v−0.6,∂v∂t=Dv∇2v−u−0.6+9.99v,where Du and Dv 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 atDu=1, Dv=1.75, R=20, and r=R/3. Figure 4 shows the initial condition u0 and the approximation of the final state ux,T, with T=100.Figure 4
(a) Initial distributionu0. (b) Final distribution u. Case: equation (21), Du=1, Dv=1.75, T=100, R=20, and r=R/3.
(a)(b)In Figure5, we show a different initial condition u0 and the approximations of the state ux,T, with T=100 and T=300.Figure 5
(a) Initial distributionu0. (b) Distribution u at T=100. (c) Final distribution u at T=300. Case: equation (21), Du=1, Dv=1.75, R=20, and r=R/3.
(a)(b)(c)Example 8.
Now, we choose the reaction-diffusion system proposed by Barrio et al. [7],
(22)∂u∂t=Du∇2u+αu1−r1v2+v1−r2u,∂v∂t=Dv∇2v+βv1+αr1βuv+uγ+r2v,where Du and Dv 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 r1 and r2 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 atDu=2.322×10−3, Dv=0.0045, α=0.899, β=−0.91, γ=−α, r1=3.5, and r2=0. The major and minor radii that we consider are R=1 and r=0.3, respectively. Figure 6 shows the initial conditions u0 and the approximations of the state ux,T, with T=200.Figure 6
(a) Initial distributionu0. (b) Final distribution u. Case: equation (22), Du=2.322×10−3, Dv=0.0045, T=200, R=1 and r=0.3.
(a)(b)Figure7 shows another initial condition u0 and the approximations of the state ux,T, with T=100, 200, and 300.Figure 7
(a) Initial distributionu0. (b) Final distributions u. Case: equation (22); Du=2.322×10−3; Dv=0.0045; T=100,200,300; R=1; and r=0.3.
(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:
(23)K=cosθRr+r2cosθ.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 Figure8, two different views of each experiment are shown for the final distributions u of the reaction-diffusion system (22) with different radii, corresponding to different surfaces in the curvature landscape.Figure 8
Final distributionux,T with (a) R=0.75, r=0.45, and T=30; (b) R=0.9, r=0.35, and T=100; (c) R=0.95, r=0.25, and T=200. Case: equation (22), Du=2.322×10−3, Dv=0.0045.
(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].Figure 9
(a) Curvature as a function ofθ for different values of radii. For certain values of R decreasing and r increasing, the minimum of K in π becomes more pronounced. (b) Curvature landscape for different radii at θ=π. The tori of different radii in the right figure are shown as points with the indicated color code.
(a)(b)


## 5.2. 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
(24)∂u∂t=Du∇2u+a−u+u2v,∂v∂t=Dv∇2v+b−u2v,where Du and Dv 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 chemicalsu and v, produced at rates a and b, respectively.For the first case of this example, we fix the parameters atDu=1/R, Dv=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.Figure 10
Final distributionu: (a) T=100, (b) T=300, (c) T=400, and (d) T=600. Case: equation (24) R =40, Du=1/40, Dv=1/4, a=0.01, and b=1.2.
(a)(b)(c)(d)For the next example, we consideredDu=1/20, Dv=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.Figure 11
Distributionux,T. Case: R=20, Du=1/20, Dv=1, a=0.1, and b=0.9.Example 10.
Next, we will consider the Gierer-Meinhardt model with saturation, given as follows:(25)∂u∂t=Du∇2u+ρuu2v1+ku2−μuu,∂v∂t=Dv∇2v−ρvu2v1+ku2+σv,where Du and Dv are the diffusion coefficients of the morphogens u and v, respectively.In the Gierer and Meinhardt model, speciesu 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 u2v, by u2v/1+ku2, where k is a sufficiently large constant, so 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.For the following example, we fix the parameters atρu=0.18, ρv=0.36, μu=0.08, σv=0.1, and T=100. For the discretization, we used h=0.0051 and δt=1/100. Figure 12 shows the initial conditions u0 and v0 and the approximations of the final state ux,100 and vx,100 for the case Du=0.0005, Dv=0.025, and k=0.Figure 12
(a) Initial distributionu0. (b) Initial condition v0. (c) Final distribution u. (d) Final state v. Case: equation (25), Du=0.0005, Dv=0.025, T=100, and k=0.
(a)(b)(c)(d)Figure13 shows the final distributions ux,100 and vx,100 for the case Du=0.0005, Dv=0.025, and k=0.35; the initial distribution is the same initial distribution as the above case.Figure 13
(a) Final distributionu. (b) Final distribution v. Case: equation (25), Du=0.0005, Dv=0.025, T=100, and k=0.35.
(a)(b)Stripe patterns are shown in Figure14; to obtain these patterns, we use the following parameters Du=0.000028, Dv=0.000168, and k=0.35. Figures 14(a) and 14(b) show the initial distributions, which are different to the above cases.Figure 14
(a) Initial distributionu0. (b) Final condition v0. (c) Final distribution u. (d) Final state v. Case: equation (25), Du=0.000028, Dv=0.000168, T=100, and k=0.35.
(a)(b)(c)(d)A variant for this example given by(26)∂u∂t=Du∇2u+ρuu2v1+ku2−μuu,∂v∂t=Dv∇2v−ρvu2v1+ku2+σv+σu,The numerical results are shown in Figure15; in order to obtained these results, we fixed the parameters at Du=0.000026, Dv=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).Figure 15
(a) Final distributionu. (b) Final distribution v. Case: equation (26), Du=0.000026, Dv=0.000182, T=100, and k=0.45.
(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ℝ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 inN 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 Figure8. 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.



---

*Source: 8898484-2021-04-12.xml*