# Origin of anisotropic diffusion in Turing Patterns

H. Koibuchi<EMAIL_ADDRESS><EMAIL_ADDRESS>F. Kato
National Institute of Technology (KOSEN), Ibaraki College, Hitachinaka, Japan.
G. Diguet  Micro System Integration Center, Tohoku University, Sendai, Japan
T. Uchimoto  Institute of Fluid Science (IFS), Tohoku University, Sendai,
Japan  ELyTMaX, CNRS-Universite de Lyon-Tohoku University, Sendai, Japan

###### Abstract

In this paper, we numerically study Turing patterns by the Finsler geometry
(FG) modeling technique on thermally fluctuating triangular lattices, which
are often used for modeling cell membranes or lipid membranes, focusing on the
origin of diffusion anisotropy. The FG modeling prescription allows us to
assume direction-dependent diffusion described by Laplacian. To implement such
diffusion anisotropy in the FG modeling, we need an internal degree of freedom
(IDF), which depends on direction and position and is controlled by some
external forces or stimuli. For such a direction-dependent IDF, we use
velocity directions corresponding to thermal fluctuations of the lattice
vertices. We find that anisotropic Turing patterns emerge in the direction
along which vertices fluctuate. In the simulations, direction-dependent
diffusion coefficients are unnecessary for input, and instead, the IDF aligns
the direction of vertex fluctuation along a direction implemented by external
stimuli. Our results and techniques provide insight into the origin of
diffusion anisotropy connected to Turing patterns.

##  I Introduction

Turing patterns (Fig. 1(a)) are described by a set of partial differential
equations composed of diffusion and reaction terms [1]. This equation has an
unstable plane wave solution in an unstable environment called Turing
instability [2, 3, 4, 5]. In this situation, the appearance of isotropic
Turing patterns can be checked numerically with suitable pairs of diffusion
constants. Moreover, if the diffusion constants are direction-dependent,
anisotropic patterns such as those on zebra emerge. However, the origin of
this anisotropy in diffusion coefficients is unclear.

Figure 1: (a) An example of a Turing pattern, (b) two-dimensional triangulated
lattice, where small balls represent lipid molecules or lumps of lipids, (c)
omni-directional or isotropic fluctuations of lipids corresponding to
isotropic diffusion, and (d) aligned or anisotropic fluctuations of lipids
corresponding to anisotropic diffusion, where the arrows represent the
fluctuation directions. The fluctuation is in-plane directions only.

We numerically study this problem on two-dimensional triangular lattices of
fixed connectivity [6, 7, 8] and dynamically-triangulated fluid lattices [9]
(Fig. 1(b)). Triangulated surfaces have been used to study the shape
morphology of lipid membranes or cell membranes [10, 11] by assuming a
Helfrich Hamiltonian, which is the curvature energy of membranes in three-
dimensional space [12]. Triangular meshes are an extension of the linear chain
model for polymers [13]. In such a model for membranes, triangle vertices
correspond to lipid molecules. For this reason, fluctuating triangular
lattices are suitable for studying the origin of zebra patterns. We use the
Finsler geometry (FG) modeling technique, which connects anisotropy in vertex
diffusion with anisotropy in Laplacian. To connect these two different
anisotropies, we need to include internal degrees of freedom (IDF) for
distance anisotropy on the triangulated lattices in the FG modeling. Another
possible origin of IDF is bacterial swarming or navigation on surfaces [14,
15]. The bacterial cell size is $0.4(\rm{\mu m})$; therefore, the cells can
thermally fluctuate. However, this movement of bacteria is mainly activated by
some environmental conditions and is not always described by the standard
diffusion equation [14, 15]. Therefore, we study in this paper whether
anisotropic Turing patterns are described by diffusion anisotropy, which is
controlled by some IDF connected to the diffusion itself.

The problem is to find the content or identity of IDF. One possible answer is
that this IDF corresponds to the direction of thermal fluctuations of vertices
of the lattice. The movement of vertices is characterized by velocity, which
has a directional degree of freedom and plays a role in IDF (Fig. 1(c),(d)).
Without a specific condition, thermal fluctuations, including the directions,
are at random leading to isotropic diffusion. Therefore, anisotropic Turing
patterns are expected to emerge only for aligned velocity along a direction
oriented by some external stimuli or forces. This velocity alignment,
representing an actual anisotropic diffusion of lattice vertices, causes
anisotropic Turing patterns on both fixed connectivity and fluid surfaces.

##  II Models

###  II.1 FitzHugh-Nagumo equation and anisotropic Turing patterns

Let $\phi(x,y)$ and $\psi(x,y)$ be scalar functions at two-dimensional space
point $(x,y)$ satisfying a set of partial differential equations called the
FitzHugh-Nagumo equation:

$\displaystyle\begin{split}\frac{\partial\phi}{\partial
t}=D_{\phi}{\mathop{}\\!\mathbin{\bigtriangleup}}\phi+f(\phi,\psi),\quad\frac{\partial\psi}{\partial
t}=D_{\psi}\mathop{}\\!\mathbin{\bigtriangleup}\psi+\gamma g(\phi,\psi),\quad
f=\phi-\phi^{3}-\psi,\;\quad g=\phi-\alpha\psi-\beta,\end{split}$ (1)

where Laplace operator is defined by
$\mathop{}\\!\mathbin{\bigtriangleup}=\frac{\partial^{2}}{\partial
x^{2}}+\frac{\partial^{2}}{\partial y^{2}}$, and $D_{\phi}$ and $D_{\psi}$ are
the diffusion constants. Other parameters $f$ and $g$ are called reaction
terms, and $\alpha$, $\beta$ and $\gamma$ are constants.

Figure 2:  (a) An isotropic configuration obtained under the condition
$(a,b)\\!=\\!(1,1)$ in the modified Laplacian in Eq. (2), and anisotropic
configurations obtained at (b) $(a,b)\\!=\\!(0.99,1)$, (c)
$(a,b)\\!=\\!(0.98,1)$ with (d) the other parameters assumed in the FitzHugh-
Nagumo equation in Eq. (1). The size of the regular square lattice is
$N\\!=\\!100^{2}$.

It is well-known that a time-dependent plane wave solution is close to a
constant solution, and the time-dependent solution possesses an unstable
behavior corresponding to Turing patterns in a specific range of constants. We
are interested in anisotropic Turing patterns, and therefore, we modify the
Laplacian to

$\displaystyle\mathop{}\\!\mathbin{\bigtriangleup}\phi\to
a\frac{\partial^{2}\phi}{\partial x^{2}}+(2-a)\frac{\partial^{2}\phi}{\partial
y^{2}},\;(0<a<2),\quad\mathop{}\\!\mathbin{\bigtriangleup}\psi\to
b\frac{\partial^{2}\psi}{\partial x^{2}}+(2-b)\frac{\partial^{2}\psi}{\partial
y^{2}},\;(0<b<2)$ (2)

both of which restore to the isotropic $\mathop{}\\!\mathbin{\bigtriangleup}$
in the limit of $a\\!\to\\!1$, $b\\!\to\\!1$. Here, we use the regular square
lattice of size $N\\!=\\!100^{2}$ to solve the steady-state solutions
satisfying the conditions ${\partial\phi}/{\partial t}\\!=\\!0$ and
${\partial\psi}/{\partial t}\\!=\\!0$ in Eq. (1). Figures 1(a)–(c) show (a)
isotropic and (b),(c) anisotropic patterns, where the assumed parameters are
shown in Fig. 1(d). We find that the anisotropic pattern in Fig. 1(c) is more
evident than that in Fig. 1(b). The reason for the appearance of such an
anisotropic pattern is that $D_{\phi}$ is effectively direction-dependent such
that
$D_{\phi}\to(D_{\phi}^{x},D_{\phi}^{y})=D_{\phi}(a,2\\!-a\\!)\\!=\\!(0.98D_{\phi},1.02D_{\phi})$
for $a\\!=\\!0.98$. Note that $(a,2\\!-a\\!)$ or $(b,2\\!-b\\!)$ is the reason
for the anisotropy in diffusion constants
$(D_{\phi}^{x},D_{\phi}^{y})\\!=\\!(aD_{\phi},(2-a)D_{\phi})$ or
$(D_{\psi}^{x},D_{\psi}^{y})\\!=\\!(bD_{\psi},(2-b)D_{\psi})$. As a result,
$(D_{\phi}^{x},D_{\phi}^{y})$ changes from isotropic to anisotropic such that
$(D_{\phi}^{x},D_{\phi}^{y})\\!=\\!(D_{\phi},D_{\phi})\to(D_{\phi}^{x},D_{\phi}^{y})\\!=\\!(0.98D_{\phi},1.02D_{\phi})$
if $(a,b)$ change as $(a,b)\\!=\\!(1,1)\\!\to\\!(a,b)\\!=\\!(0.98,1)$. Here,
we assume anisotropy only in the variable $\phi$. However, it is still unclear
the origin or mechanism for $a$ to deviate from $a\\!=\\!1$ to $a\\!=\\!0.98$.

###  II.2 Finsler geometry modeling of anisotropic diffusion

To find the reason for $(D_{\phi}^{x},D_{\phi}^{y})$ being anisotropic, we use
the FG modeling technique. In this FG modeling, a new IDF $\vec{\tau}$ is
introduced for the variables $\phi$ and $\psi$ in Eq. (1), and the diffusion
constant $D_{\psi}$ also becomes isotropic for the variable $\psi$. The
variable $\vec{\tau}$ treated as a dynamical variable in addition to $\phi$
and $\psi$, is updated by Monte Carlo (MC) simulations. For this purpose, we
need a Hamiltonian depending on $\vec{\tau}$.

Thus, we start with the following Hamiltonian:

$\displaystyle S=S_{\phi}+S_{\psi},\quad
S_{\phi}=\frac{1}{2}\int\sqrt{g}d^{2}x\,g^{ab}\frac{\partial\phi}{\partial
x^{a}}\frac{\partial\phi}{\partial x^{b}},\quad
S_{\psi}=\frac{1}{2}\int\sqrt{g}d^{2}x\,g^{ab}\frac{\partial\psi}{\partial
x^{a}}\frac{\partial\psi}{\partial x^{b}},$ (3)

The $S_{\phi}$and $S_{\psi}$ in $S$ are written by using the inverse $g^{ab}$
and the determinant $g$ of the Finsler metric $g_{ab}$, where $(x^{1},x^{2})$
represents a local coordinate of the surface. By the variation principle, we
have a differential equation for $\phi$ and $\psi$ associated with the
operator expression

$\displaystyle\mathop{}\\!\mathbin{\bigtriangleup}=\frac{1}{\sqrt{g}}\frac{\partial}{\partial
x_{a}}\left(\sqrt{g}g^{ab}\frac{\partial}{\partial x^{b}}\right),$ (4)

which is a Laplacian on curved surfaces. In this expression, $g$ and $g^{ab}$
depend on $\vec{\tau}$, which is also dependent on the local coordinate
directions $x^{1}$ and $x^{2}$. Therefore, the operator
$\mathop{}\\!\mathbin{\bigtriangleup}$ in Eq. (4) can be written as
$\mathop{}\\!\mathbin{\bigtriangleup}_{\phi}\\!=\\!a_{x}(\tau)\frac{\partial^{2}}{\partial
x^{2}}\\!+\\!a_{y}(\tau)\frac{\partial^{2}}{\partial y^{2}}\\!+\\!\cdots$ and
$\mathop{}\\!\mathbin{\bigtriangleup}_{\psi}\\!=\\!b_{x}(\tau)\frac{\partial^{2}}{\partial
x^{2}}\\!+\\!b_{y}(\tau)\frac{\partial^{2}}{\partial y^{2}}\\!+\\!\cdots$,
where $a_{x}(\tau)$, $a_{y}(\tau)$ and $b_{x}(\tau)$, $b_{y}(\tau)$ play roles
in $a$, $2\\!-\\!a$ and $b$, $2\\!-\\!b$ in
$\mathop{}\\!\mathbin{\bigtriangleup}_{\phi}$ and
$\mathop{}\\!\mathbin{\bigtriangleup}_{\psi}$ in Eq. (2), respectively. The
detailed information on the discretization will be reported elsewhere.

Figure 3: (a) Two triangles $ijk$ and $ilj$ sharing bond $ij$, and a local
coordinate $(x^{1},x^{2})$ of triangle $ijk$ at vertex $i$, (b) a line
$C(x^{1})$ along the $x^{1}$ axis, one-dimensional energy $E_{\phi}(\vec{r})$
for Finsler length along $C(x^{1})$ and the corresponding IDF $\vec{\tau}$ at
$\vec{r}$, (c) Expressions of $E_{\phi}(\vec{r})$, Finsler function $L_{E}$
for $\phi$ and the corresponding Finsler metric $g$ along $C(x^{1})$. The
definition of $v_{ij}$ for $\phi$ is different from that of $v_{ij}$ for
$\psi$ as shown in Eq. (5).

Using $S_{\phi}$ and $S_{\psi}$, we introduce two different Finsler time
lengths in the local coordinate system on the lattice. First we introduce one-
dimensional discrete energies $E_{\phi}$ and $E_{\psi}$ corresponding to
$S_{\phi}$ and $S_{\psi}$ to define the Finsler length $v_{ij}$ for $\phi$ and
$\psi$

$\displaystyle\begin{split}&E_{\phi}=\sum_{ij\in
C(x^{1})}\left(\phi_{j}-\phi_{i}\right)^{2},\quad
v_{ij}^{\phi}=|\vec{e}_{ij}\cdot\vec{\tau}_{i}|+v_{0}\\\ &E_{\psi}=\sum_{ij\in
C(x^{1})}\left(\psi_{j}-\psi_{i}\right)^{2},\quad
v_{ij}^{\psi}=\sqrt{1-(\vec{e}_{ij}\cdot\vec{\tau}_{i})^{2}}+v_{0}\end{split}$
(5)

where $ij$ denotes a bond along the line $C(x^{1})$ representing a local
coordinates axis (Figs. 3(a),(b)). This energy is used to define the Finsler
time length $s$ along the $x^{1}$ axis by using the Finsler function
$L_{E}\\!=\\!\sqrt{E^{\prime 2}}/|v|$ such that $L_{E}\\!=\\!ds/dx^{1}$, where
$v$ is a velocity or Finsler length along the $x^{1}$ axis (Fig. 3(c)) (see
also Ref.[16] for more detailed information on this part). To define $v$, we
introduce a variable $\vec{\tau}(\in S^{1}/2)$, which will be described in the
following subsection. The Finsler length along the $x^{2}$ axis is defined in
the same way as that along the $x^{1}$ axis and is implemented in the model
through the metric tensor $g_{ab}$ such that
$(g_{11},g_{22})\\!=\\!(v_{12}^{-2},v_{13}^{-2})$. Finsler lengths $v_{ij}$
for $\phi$ and $\psi$ denoted by $v_{ij}^{\phi}$ and $v_{ij}^{\psi}$ are
defined by using the unit tangential vector $\vec{e}_{ij}$ from $i$ to $j$ and
$\vec{\tau}_{i}$ at $i$ as in Eq. (5). Note that $v_{ij}^{\phi}$
($v_{ij}^{\psi}$) is long along the direction parallel (perpendicular) to the
$\vec{e}_{ij}$ direction and is controlled by some external force $\vec{B}$
because the direction of $\vec{\tau}_{i}$ is controlled by $\vec{B}$, which
will be introduced below.

Using those obtained metric tensors corresponding to IDF $\vec{\tau}$
including several additional energies defined on triangular lattice, we obtain
the discrete Hamiltonian [17]

$\displaystyle\begin{split}&S(\vec{r};\vec{\tau})=S_{1}+D_{\phi}S_{\phi}+D_{\psi}S_{\psi}+\lambda_{\phi}S_{0}^{\phi}+\lambda_{\psi}S_{0}^{\psi}+S_{B}^{\phi},\\\
&S_{1}=\sum_{ij}\ell_{ij}^{2},\quad\ell_{ij}^{2}\\!=\\!(\vec{r}_{i}\\!-\\!\vec{r}_{i})^{2},\\\
&S_{\phi}=\sum_{i}\sum_{j(i)}\gamma_{ij}^{\phi}\left(\phi_{j}-\phi_{i}\right)^{2},\quad
S_{\psi}=\sum_{i}\sum_{j(i)}\gamma_{ij}^{\psi}\left(\psi_{j}-\psi_{i}\right)^{2},\quad
S_{B}^{\phi}=-\sum_{i}\left(\vec{\tau}_{i}\cdot{\vec{B}}\right)^{2},\\\
&\gamma_{ij}^{\phi}=\frac{1}{6}\left(\frac{v_{ij}^{\phi}}{v_{ik}^{\phi}}+\frac{v_{ji}^{\phi}}{v_{jk}^{\phi}}+\frac{v_{ij}^{\phi}}{v_{il}^{\phi}}+\frac{v_{jl}^{\phi}}{v_{jl}^{\phi}}\right),\quad\gamma_{ij}^{\psi}=\frac{1}{6}\left(\frac{v_{ij}^{\psi}}{v_{ik}^{\psi}}+\frac{v_{ji}^{\psi}}{v_{jk}^{\psi}}+\frac{v_{ij}^{\psi}}{v_{il}^{\psi}}+\frac{v_{jl}^{\psi}}{v_{jl}^{\psi}}\right),\\\
\end{split}$ (6)

where $ij$ in $\sum_{ij}$ in $S_{1}$ denotes the bond connecting two
neighboring vertices $i$ and $j$ and
$\ell_{ij}^{2}\\!=\\!(\vec{r}_{i}\\!-\\!\vec{r}_{i})^{2}$ is the bond length
squares. $\sum_{i}\sum_{j(i)}$ in $S_{\phi}$ and $S_{\psi}$ denotes the sum
over vertex $i$ and the sum over bonds $j(i)$ connected to $i$. The first term
$S_{1}$ is the bond potential, which makes the mean value of $\ell_{ij}$
uniform. The symbol ${\vec{B}}$ in the term $S_{B}^{\phi}$ denotes some
external stimuli or forces aligning the IDF $\vec{\tau}_{i}$. We should note
that $\vec{\tau}_{i}$ is forced to align along the direction of ${\vec{B}}$
when $|{\vec{B}}|\\!\not=\\!0$, and as a result, thermal fluctuations of
$\vec{r}$ also align along ${\vec{B}}$. More detailed information on the
Hamiltonian, including the discretization of $S_{\phi}$ and $S_{\psi}$, will
be reported elsewhere.

###  II.3 Hybrid numerical technique

Figure 4: (a) Initial configuration of a triangulated lattice with periodic
boundary conditions and of size $N\\!=\\!20^{2}$, which is the total number of
vertices smaller than $N\\!=\\!50^{2}$ assumed in the simulations. (b) An
illustration of the bond flip process, which is a basic process of dynamical
triangulation $\mathcal{T}$ in $Z_{\rm fluid}$ and allows vertices $i$ and $j$
to diffuse over the surface, (c) an illustration of a vertex move from
$\vec{r}({\rm old})$ to $\vec{r}({\rm new})$, which is a basic process of
thermal fluctuations in both $Z_{\rm fixed}$ and $Z_{\rm fluid}$, and (d)
definition of IDF $\vec{\tau}$ by using $\Delta\vec{r}=\vec{r}({\rm
new})\\!-\\!\vec{r}({\rm old})$. Vertices on the upper and left boundaries in
(a) are identified with those on the opposite boundaries. These vertices
fluctuate like the internal vertices in both $Z_{\rm fixed}$ and $Z_{\rm
fluid}$, and the boundary bonds also flip like the internal bonds in $Z_{\rm
fluid}$.

As indicated in $S(\vec{r};\vec{\tau})$, the variable to be updated by MC
simulations is the vertex position $\vec{r}$ of the lattice, and IDF
$\vec{\tau}$ is defined by the deformation of $\vec{r}$. Therefore, the
variables $\\{\vec{r}\\}$ are the source of IDF $\\{\vec{\tau}\\}$, or we
simply consider that the vertex position $\vec{r}$ plays a role in the IDF
$\vec{\tau}$ for the diffusion anisotropy for $\phi$ and $\psi$. Here, we
should note that vertex move in MC does not correspond to the velocity of the
vertex because MC iterations of the variable $\vec{r}$ are not the time
evolution of $\vec{r}$. In this sense, vertex fluctuations realized in the MC
process do not exactly correspond to thermal fluctuations of $\vec{r}$ the
position of lipids, although the mean value of physical quantity obtained by
those fluctuating ensemble configurations of vertices is expected to be a
meaningful quantity. Thus, we assume that vertex fluctuations in the MC
process represent real thermal fluctuations. Hence, we consider that movement
of lipid molecules or lumps of lipids is reflected in the IDF $\vec{\tau}$.

We use two different types of lattices, fixed connectivity and fluid, of
periodic boundary conditions to update the IDFs in MC (Fig.4(a)). The
partition functions are given by

$\displaystyle\begin{split}&Z_{\rm
fixed}=\int\prod_{i}d\vec{r}_{i}\exp(-S),\quad(\textrm{fixed model}),\\\
&Z_{\rm
fluid}=\sum_{\mathcal{T}}\int\prod_{i}d\vec{r}_{i}\exp(-S),\quad(\textrm{fluid
model}).\end{split}$ (7)

The symbol $\sum_{\mathcal{T}}$ in $Z_{\rm fluid}$ denotes the sum over
possible triangulation ${\mathcal{T}}$, which is performed by the bond-flip
technique [9] (Fig. 4(b)). For this bond-flip process of MC updates for
triangulation ${\mathcal{T}}$, the connectivity of triangles changes, and as a
result, vertices freely diffuse over the surface. We should note that the
simulation unit is defined by $k_{B}T\\!=\\!1$, where $k_{B}$ and $T$ are the
Boltzmann constant and the temperature, respectively.

Thus, the numerical procedure is as follows:

  1. (i)

Discrete versions of Eq. (1) are iterated such that
$\phi(t+dt)\leftarrow\phi(t)+dt(\cdots)$ and
$\psi(t+dt)\leftarrow\psi(t)+dt(\cdots)$ with the operators
$\mathop{}\\!\mathbin{\bigtriangleup}_{\phi}$ and
$\mathop{}\\!\mathbin{\bigtriangleup}_{\psi}$ in Eq. (4)

  2. (ii)

In each discrete time step $t\to t+dt$, the variable $\vec{r}$ is updated in
MC simulations, as shown in Fig. 4(c), for sufficiently many sweeps on a fixed
connectivity lattice (fixed model) or dynamically-triangulated fluid lattices
(fluid model), which are defined by the partition functions in Eq. (7). At
each update of $\vec{r}$, the variable $\vec{\tau}$ is also updated by using
$\Delta\vec{r}(=\vec{r}({\rm new})\\!-\\!\vec{r}({\rm old})$ (Fig. 4(d)).

  3. (iii)

Steps (i) and (ii) are repeated until the conditions $\partial\phi/\partial
t\\!=\\!0$, and $\partial\psi/\partial t\\!=\\!0$ are numerically satisfied in
Eq. (1), in which the Laplace operator
$\mathop{}\\!\mathbin{\bigtriangleup}=\frac{\partial^{2}}{\partial
x^{2}}+\frac{\partial^{2}}{\partial y^{2}}$ is replaced by
$\mathop{}\\!\mathbin{\bigtriangleup}_{\phi}$ for $\phi$ and
$\mathop{}\\!\mathbin{\bigtriangleup}_{\psi}$ for $\psi$ defined in Eq. (4).

##  III Results

Figure 5: (a) Isotropic and (b) anisotropic patterns on a fixed connectivity
surface, and (c) isotropic and (d) anisotropic patterns on the fluid surface.
The assumed parameters are shown on the top of the snapshots. The small cones
represent IDF $\tau$ indicating the direction of vertex diffusion. The size of
the lattice is $N\\!=\\!50^{2}$.

The simulations are performed on triangulated lattices of periodic boundary
conditions of size $N\\!=\\!50^{2}$. The results of the fixed connectivity and
fluid models are shown in Figs. 5(a), (b) and Figs. 5(c), (d), respectively.
The assumed parameters $\alpha,\beta,\gamma$, $D_{\phi},D_{\psi}$ in Eqs. (1),
$\lambda_{\phi},\lambda_{\psi},\vec{B}$ in $S$ of Eq. (6), and $v_{0}$ in
$\vec{v}^{\,\phi,\psi}$ in Eq. (5) are shown in the upper part of each
snapshot. We find that isotropic (anisotropic) patterns appear when
$\vec{B}\\!=\\!(B,0)$ with $B\\!=\\!0$ ($B\\!>\\!0$) in both fixed and fluid
surfaces in Figs. 5(a), (c) (5(b),(c)), where $\lambda_{\phi}$ is fixed to be
nonzero for anisotropic patterns. The direction of $\vec{\tau}$ controlled by
$\vec{B}$ is closely related to whether the pattern is isotropic or
anisotropic. This correspondence between the directions of $\vec{\tau}$ and
the patterns explains the anisotropic Turing patterns. Because anisotropy in
$\vec{\tau}$ corresponds to anisotropy in vertex diffusion, and the
anisotropic Turing pattern is a consequence of anisotropy in diffusion
constants as confirmed in Fig. 2, we consider that anisotropy in vertex
diffusion is reflected in the anisotropy in Turing patterns.

##  IV Concluding Remarks

Preliminary simulation results on anisotropic Turing patterns are presented in
this paper. The FitzHugh-Nagumo equation is numerically solved on fixed
connectivity and dynamically triangulated surfaces. The vertex position of
these lattices plays a role in an internal degree of freedom (IDF) for Finsler
geometry (FG) modeling of anisotropy in the Laplacian in the FitzHugh-Nagumo
equation, and this IDF is updated in the Metropolis Monte Carlo simulation
technique. We find that anisotropy in Laplacian or diffusion anisotropy is
connected with anisotropy in vertex diffusion on triangulated lattices via the
IDF, implying that Turing patterns can be controlled by the vertex diffusion.
Detailed information on the induced anisotropy in diffusion constants of
Laplacian should be studied, and further numerical studies are necessary for
more detailed information on the relation between Turing patterns and
anisotropy in IDF.

###### Acknowledgements.

The author H.K. acknowledges Dr. M. Nakayama and Dr. S. Tasaki for comments.
The authors acknowledge Dr. Jean-Paul Rieu for suggestions. Numerical
simulations were performed on the Supercomputer system "AFI-NITY" at the
Advanced Fluid Information Research Center, Institute of Fluid Science, Tohoku
University.

## References

  * [1] A. M. Turing, Phil. Trans. Roy. Soc. London 237, pp.37-72 (1952).
  * [2] Y. Kuramoto, “Dynamics of Reaction and Diffusion”, Seibutu Butsuri (in Japanese) . 19, pp.11-18 (1979).
  * [3] J. W. Chiu and K.-H. Chiam, Phys. Rev. E. 78, 056708 (2008).
  * [4] S. Kondo and R. Asai, Letters to Nature 376, pp.765-768 (1995).
  * [5] P. Dziekan, A. Lemarchand, and B. Nowakowski, J. Chem. Phys. 137, 074107 (2012).
  * [6] Y. Kantor and D.R. Nelson, Phys. Rev. A 36, 4020 (1987)
  * [7] G. Gompper and D. M. Kroll, Phys. Rev. A 46, 7466 (1992).
  * [8] H. Koibuchi and T. Kuwahata, Phys. Rev. E 72, 026124 (2005).
  * [9] J. -S. Ho and A. Baumg${\ddot{\rm a}}$rtner, Europhys. Lett. 12, 295 (1990).
  * [10] L. Peliti and S. Leibler, Phys. Rev. Lett. 54, 1690 (1985).
  * [11] D. Nelson, in Statistical Mechanics of Membranes and Surfaces, Second Edition, eds. D. Nelson, T. Piran, and S. Weinberg (World Scientific, Singapore, 2004), p. 1.
  * [12] W. Z. Helfrich, Naturforsch 28c, 693 (1973).
  * [13] M. Doi and F. Edwards, The Theory of Polymer Dynamics (Oxford University Press, 1986).
  * [14] G. Ariel, A. Shklarsh, O. Kalisman, C. Ingham and E. Ben-Jacob, New J. Phys. 15, 125019 (2013),
  * [15] G. Ariel, A. Rabani, S. Benisty, J. D. Partridge, R. M. Harshey and A.Be’er, Nature Comm. 6: 8936 (2015), DOI: 10.1038/ncomms9396.
  * [16] S. El Hog, F. Kato, S. Hongo, H. Koibuchi, G. Diguet, T. Uchimoto and H. T. Diep, Results in Phys. 38, 105578 (2022), https://doi.org/10.1016/j.rinp.2022.105578.
  * [17] H. Koibuchi, M. Okumura and S. Noro, J. Phys. Conf. Ser. 1730, 012035 (2020), doi:10.1088/1742-6596/1730/1/012035.