# Numerical study of anisotropic diffusion in Turing patterns based on Finsler
geometry modeling

Gildas Diguet  Micro System Integration Center, Tohoku University, Sendai,
Japan  Madoka Nakayama  Research Center of Mathematics for Social Creativity,
Research Institute for Electronic Science, Hokkaido University, Sapporo, Japan
Sohei Tasaki  Department of Mathematics, Faculty of Science, Hokkaido
University, Sapporo, Japan  Fumitake Kato  Hiroshi Koibuchi
<EMAIL_ADDRESS><EMAIL_ADDRESS>National Institute of
Technology (KOSEN), Ibaraki College, Hitachinaka, Japan.  Tetsuya Uchimoto
Institute of Fluid Science (IFS), Tohoku University, Sendai, Japan  ELyTMaX,
CNRS-Universite de Lyon-Tohoku University, Sendai, Japan

###### Abstract

We numerically study the anisotropic Turing patterns (TPs) of an activator-
inhibitor system, focusing on anisotropic diffusion by using the Finsler
geometry (FG) modeling technique. In the FG modeling prescription, the
diffusion coefficients are dynamically generated to be direction dependent
owing to an internal degree of freedom (IDOF) and its interaction with the
activator and inhibitor under the presence of thermal fluctuations. In this
sense, FG modeling contrasts sharply with the standard numerical technique,
where direction-dependent diffusion coefficients are assumed in the reaction-
diffusion (RD) equations of Turing. To find the solution of the RD equations,
we use a hybrid numerical technique as a combination of the metropolis Monte
Carlo method for IDOF updates and discrete RD equations for steady-state
configurations of activator-inhibitor variables. We find that the newly
introduced IDOF and its interaction are one possible origin of spontaneously
emergent anisotropic patterns on living organisms such as zebra and fishes.
Moreover, the IDOF makes TPs controllable by external conditions if the IDOF
is identified with lipids on cells or cell mobility.

##  I Introduction

Turing patterns (TPs) are described by partial differential equations of
Turing Turing-TRSL1952  using two different variables $u$ and $v$, which are
scalar functions on a domain in ${\bf R}^{2}$ in the case of a two-dimensional
system. These $u$ and $v$ parameters are usually called the activator and
inhibitor, respectively Koch-Meinhardt-RMP1994 ; FitzHugh-BP1961 ; Nagumo-
etal-ProcIRE1962 , owing to their interaction properties as implemented in the
reaction and diffusion (RD) terms in the equation. TPs emerge on scales
ranging from macroscopic Sekimura-etal-PRSL2000 ; Kondo-Miura-Science2010 ;
Bullara-Decker-NatCom2015  to microscopic Tan-etal-Science2018 ; Fuseya-etal-
NatPhys2021 .

These patterns emerge as a result of competition between diffusion and
reaction, and therefore, many studies have been conducted to extend the
Laplace operators of the RD equation to the graph Laplacian to find TPs in
random networks Nakao-etal-NatPhys2010 ; Carletti-Nakao-PRE2020 ; Asllani-
etal-NatCom2014 ; Asllani-etal-PRE2014 ; Petit-etal-PhysA2016 . In those
networks, topological properties rather than geometric properties play a
significant role in forming TPs. Pattern formation depends on the node number,
which is the total number of connections, and TPs can be visualized by using a
“coordinate axis” of node numbers. Another extension is to modify diffusivity
to accommodate non-Gaussian behavior of Brownian particles confined in narrow
plates by including fluctuations in diffusion constants, where position-
dependent and anisotropic diffusion constants are assumed Alexandre-etal-
PRL2023 . Such a non-Gaussian distribution of particle displacements is
considered anomalous diffusion corresponding to anomalous transport phenomena
in crowded biological materials such as cellular membranes Bressloff-Newby-
RMP2013 ; Hofling-Franosch-RPP2013 ; Sokolov-etal-PhysToday2002 ; Metzler-
Klafter-JPhysA2004 . This anomalous transport is characterized by
subdiffusion, which is described by a power law behavior of the mean-square
displacement $\sim t^{\alpha}(0\\!<\\!\alpha\\!<\\!1)$, observed at
intermediate time scales. So-called super diffusion characterized by
$t^{\alpha}(\alpha\\!>\\!1)$ is observed in bacterial swarming Ariel-etal-
NJP2013 ; Ariel-etal-NatCom2015 . This phenomenon is described by the Levy
walk, which is a model of a random walk with a constant speed.

Figure 1:  Anisotropic Turing patterns on (a) zebra and (b), (c) fishes. The
patterns vary depending on the individual, while the direction of patterns
depends only on the species.

Anisotropic TPs observed on zebra and fishes (Fig. 1) are known to emerge as a
consequence of a difference in diffusion constants between the activator and
inhibitor. This difference in diffusion constants corresponding to marine
angelfish was estimated by Kondo and Asai in Ref. Kondo-Nature1995 . In Refs.
Shoji-etal-DevDyn2003 ; Iwamoto-Shoji-RIMS2018 , Shoji et al. reported that
anisotropy in diffusion constants is effective in determining the direction of
stripe patterns in numerical studies. TPs appear on curved surfaces Varea-
etal-PRE1999 . Krause et al. reported that anisotropy in the evolution of TPs
is sensitive to curvature on growing domains, which are two-dimensional curved
surfaces embedded in ${\bf R}^{3}$, assuming the induced metric for describing
the Laplace Beltrami operator in the diffusion terms Krause-etal-
BulMatBiol2019 .

Anisotropic diffusion is ubiquitously expected in many phenomena. Particles in
cubic crystals under a stress field undergo direction-dependent diffusion
Dederichs-Schroeder-PRB1978 , and anisotropic surface diffusion of CO on
Ni(110) has been experimentally observed Xiao-etal-PRL1991 . Light propagates
faster in one direction than in the other directions in strongly scattering
media Bret-Lagendijk-PRE2004 , and skyrmion transport is also direction
dependent, accompanying a shape deformation due to an applied in-plane
magnetic field Kerber-etal-PRAP2021 . All these phenomena, including
anisotropic TPs such as those in Fig. 1, are well described by assuming
direction-dependent diffusion constants because the underlying origins of
anisotropy are mostly well understood. However, such a phenomenological
description is not always satisfactory and can be further improved.

In this paper, we study the origin of anisotropic diffusion in TPs. For this
purpose, we focus on a mechanism of dynamic anisotropy in which diffusion
constants dynamically appear as direction-dependent constants, allowing us to
understand anisotropic diffusion without assuming these constants
ICMsquare2020-JPCconf2022 ; Koibuchi-etal-ICMsquare2022 . Here, we recall that
such dynamic anisotropy is successfully implemented in interaction
coefficients in several statistical mechanical models using the Finsler
geometry (FG) modeling technique Koibuchi-Sekino-PhysA2014 ; Takano-Koibuchi-
PRE2017 ; Proutorov-etal-JPC2018 ; ElHog-etal-PRB2021 ; ElHog-etal-RinP2022 .
This technique simply changes length scales for interactions from the
Euclidean length to position- and direction-dependent lengths along the local
coordinate systems, and detailed information on new interactions is
unnecessary. Thus, we consider that the FG modeling technique takes advantage
of phenomenological modeling. Indeed, FG modeling prescription is applicable
to Laplace operators in the diffusion terms of the RD equation. This
applicability comes from the fact that differential equation models such as RD
equations share the same property in their interactions with those statistical
mechanical models in which FG modeling is successful. In addition, the
stochastic nature in the internal degree of freedom (IDOF) introduced in Monte
Carlo (MC) methods is shared with the convergent configurations of activator-
inhibitor variables, and hence, physical quantities are obtained by
calculating the mean values of many convergent steady-state configurations of
the RD equation.

This paper is organized as follows: In Section II, we review the RD equation
of the FN type for the variables $u$ and $v$ from a numerical point of view
and show numerical data such as snapshots of isotropic and anisotropic TPs
obtained on a regular square lattice of size $N\\!=\\!100^{2}$, where
anisotropic TPs appear as a result of assumed direction-dependent diffusion
coefficients in the RD equations. To quantify the anisotropic TPs and evaluate
the effect of the direction-dependent diffusion coefficients, we introduce
absolute second-order derivatives of $u$ and $v$ as well as squares of the
first-order derivatives. In Section III, we introduce the FG modeling
technique to implement the dynamical anisotropy in the Laplace operator in the
RD equation by including a new IDOF on two types of triangulated lattices:
fixed-connectivity and dynamically triangulated lattices. A numerical
technique, which we call the hybrid technique, is introduced to update $u$ and
$v$ and the new IDOFs. In Section IV, we show numerical data including
snapshots obtained on these two lattices. Finally, we summarize the results in
Section V, where remaining problems are also mentioned.

##  II Standard approach to Turing patterns

In this section, we review TPs described by the FitzHugh-Nagumo (FN) equation
and show snapshots of isotropic and anisotropic TPs and related physical
quantities, which we call absolute second-order partial differentials,
obtained by standard numerical techniques on regular square lattices with
periodic boundary conditions (PBCs). In Fig. 2(a), we show a small lattice of
size $N\\!=\\!10^{2}$ to visualize the lattice structure, and a lattice site
$(i,j)$ and its four neighboring sites are illustrated in Fig. 2(b). Only a
regular square lattice is considered in this section.

Figure 2:  (a) Regular square lattice with PBCs of size $N\\!=\\!10^{2}$,
which is small compared with $N\\!=\\!100^{2}$ assumed in the simulations in
the following subsection, and (b) lattice site $(i,j)$ and its four nearest
neighbor sites, where
$1\\!\leq\\!i\\!\leq\\!n_{x},1\\!\leq\\!j\\!\leq\\!n_{y}$. The lattice spacing
is assumed to be $\Delta x\\!=\\!\Delta y\\!=\\!1$ in the simulations.

###  II.1 FitzHugh-Nagumo equation with diffusion anisotropy

Let $u(x,y)$ and $v(x,y)$ be the variables corresponding to the activator and
inhibitor, respectively, which satisfy the RD equations of FN type

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

on the two-dimensional plane Shoji-etal-DevDyn2003 ; Iwamoto-Shoji-RIMS2018 .
The first and second terms on the right-hand side are called the diffusion and
reaction terms, respectively, where
$\mathop{}\\!\mathbin{\bigtriangleup}\\!=\\!\frac{\partial^{2}}{\partial
x^{2}}\\!+\\!\frac{\partial^{2}}{\partial x^{2}}$ is the Laplace operator. The
symbols $D_{u}$ and $D_{v}$ are the diffusion coefficients, and $\alpha$ and
$\gamma$ are constants.

For suitable ranges of the parameters, the RD equations have certain steady-
state solutions called TPs with a periodicity in spatial directions. When the
periodicity appears almost regularly in one direction, the patterns become
anisotropic, as shown in Fig. 1(a)–(c). These anisotropic patterns can be
reproduced by using the RD equation with diffusion anisotropy introduced with
the parameters $a$ and $b$ such that Iwamoto-Shoji-RIMS2018

$\displaystyle\begin{split}&\mathop{}\\!\mathbin{\bigtriangleup}u\to{\mathop{}\\!\mathbin{\bigtriangleup}}_{a}u=a\frac{\partial^{2}u}{\partial
x^{2}}+(2-a)\frac{\partial^{2}u}{\partial y^{2}},\;(0<a<2),\\\
&\mathop{}\\!\mathbin{\bigtriangleup}v\to{\mathop{}\\!\mathbin{\bigtriangleup}}_{b}v=b\frac{\partial^{2}v}{\partial
x^{2}}+(2-b)\frac{\partial^{2}v}{\partial y^{2}},\;(0<b<2).\end{split}$ (2)

For the ranges $0\\!<\\!a\\!<\\!2$ and $0\\!<\\!b\\!<\\!2$, the diffusion
constants $D_{u}$ and $D_{v}$ in Eq. (LABEL:FN-eq-Eucl) effectively becomes
direction dependent such that

$\displaystyle\begin{split}&D_{u}\to\left(D_{u}^{x},D_{u}^{y}\right)=\left(aD_{u},(2-a)D_{u}\right),\\\
&D_{v}\to\left(D_{v}^{x},D_{v}^{y}\right)=\left(bD_{u},(2-b)D_{v}\right).\end{split}$
(3)

The direction-dependent coefficients can be represented by the ratio

$\displaystyle c^{u}_{v}=\frac{c_{u}}{c_{v}}=\frac{a(2-b)}{b(2-a)},\quad{\rm
where}\quad c_{u}=\frac{D_{u}^{x}}{D_{u}^{y}}=\frac{a}{2-a},\quad
c_{v}=\frac{D_{v}^{x}}{D_{v}^{y}}=\frac{b}{2-b}.$ (4)

We call $c^{u}_{v}$ the anisotropy coefficient or simply anisotropy. Note that
the conditions $a\\!=\\!b\\!=\\!1$ correspond to the isotropic diffusion
represented by the anisotropy $c^{u}_{v}\\!=\\!1$, which implies isotropy. The
condition $c^{u}_{v}\\!=\\!1$ is expected for all $a\\!=\\!b$ even when
$a\\!\not=\\!1$, where $c_{u}\\!\not=\\!1$ and $c_{v}\\!\not=\\!1$.

###  II.2 Numerical solutions on a regular square lattice with periodic
boundary conditions

The time evolution equations in Eq. (LABEL:FN-eq-Eucl) are replaced by the
discrete time evolution equations using the discrete time step $\Delta t$ such
that

$\displaystyle\begin{split}&u_{ij}(t+{\Delta}t)\leftarrow
u_{ij}(t)+{\Delta}t\left[D_{u}{\mathop{}\\!\mathbin{\bigtriangleup}}_{a}u_{ij}(t)+f\left(u_{ij}(t),v_{ij}(t)\right)\right],\\\
&v_{ij}(t+{\Delta}t)\leftarrow
v_{ij}(t)+{\Delta}t\left[D_{v}{\mathop{}\\!\mathbin{\bigtriangleup}}_{b}v_{ij}(t)+g\left(u_{ij}(t),v_{ij}(t)\right)\right],\end{split}$
(5)

where $u_{ij}$ and $v_{ij}$ denote the discrete analogs of $u(x,y)$ and
$v(x,y)$ defined at lattice site $(i,j)$,
$(1\\!\leq\\!i\\!\leq\\!n_{x},1\\!\leq\\!j\\!\leq\\!n_{y})$ (Fig. 2(a)). The
discrete diffusion terms ${\mathop{}\\!\mathbin{\bigtriangleup}}_{a}u_{ij}$
and ${\mathop{}\\!\mathbin{\bigtriangleup}}_{b}v_{ij}$ are given by

$\displaystyle\begin{split}&{\mathop{}\\!\mathbin{\bigtriangleup}}_{a}u_{ij}=\frac{a}{(\Delta
x)^{2}}\left(u_{i+1,j}+u_{i-1,j}-2u_{i,j}\right)+\frac{2-a}{(\Delta
y)^{2}}\left(u_{i,j+1}+u_{i,j-1}-2u_{i,j}\right),\\\
&{\mathop{}\\!\mathbin{\bigtriangleup}}_{b}v_{ij}=\frac{b}{(\Delta
x)^{2}}\left(v_{i+1,j}+v_{i-1,j}-2v_{i,j}\right)+\frac{2-b}{(\Delta
y)^{2}}\left(v_{i,j+1}+v_{i,j-1}-2v_{i,j}\right)\end{split}$ (6)

for all $(i,j)$ with lattice spacing $\Delta x$ and $\Delta y$, where the
positions $(i\\!\pm\\!1,j)$ and $(i,j\\!\pm\\!1)$ are shown in Fig. 2(b). The
convergent criteria of the iterations in Eqs. (5) are given by

$\displaystyle\begin{split}&{\rm Max}\left\\{|u_{ij}(t+\Delta
t)-u_{ij}(t)|\right\\}<1\times 10^{-8},\\\ &{\rm Max}\left\\{|v_{ij}(t+\Delta
t)-v_{ij}(t)|\right\\}<1\times 10^{-8}.\end{split}$ (7)

Figure 3: Snapshots of the standard model on a regular square lattice of size
$N\\!=\\!100^{2}$ with PBCs. The anisotropy parameters $(a,b)$ fixed to (a)
$(1,1)$, (b) $(1.02,1)$, (c) $(0.98,1)$, (d) $(1.02,1.02)$, (e) $(1,1.02)$,
(f) $(1,0.98)$ and (g) $(1.02,0.98)$, and the anisotropies $c_{v}^{u}$ in Eq.
(4) are shown. The pattern is isotropic when $a\\!=\\!b$, as in (a) and (d),
represented by $c_{v}^{u}\\!=\\!1$, while the stripe lies along the $x$
direction when $a\\!>\\!b$, as in (b), (f) and (g), represented by
$c_{v}^{u}\\!>\\!1$, and it lies along the $y$ direction when $a\\!<\\!b$, as
in (c) and (e), represented by $c_{v}^{u}\\!<\\!1$.

Snapshots of the convergent configurations of $u$ on the lattice of size
$N\\!=\\!100^{2}$ ($\Leftrightarrow n_{x}\\!=\\!n_{y}\\!=\\!100$) are plotted
in Fig. 3(a)–(g), where $\Delta t\\!=\\!1\times 10^{-3}$ and $\Delta
x\\!=\\!\Delta y\\!=\\!1$. The value of $u$ in the snapshots is normalized to
$[0,1]$ to visualize the patterns, and the corresponding color code is shown,
as are the assumed parameters. We confirm that the patterns are isotropic when
$a\\!=\\!b$ in (a) $a\\!=\\!b\\!=\\!1$ and (d) $a\\!=\\!b\\!=\\!1.02$, while
they are anisotropic when $a\\!>\\!b$ in (b) and (f) and $a\\!<\\!b$ in (c)
and (e). This change from isotropy to anisotropy is enhanced when the
anisotropy $c_{v}^{u}$ is increased from $c_{v}^{u}\\!=\\!1$ in (a) to
$c_{v}^{u}\\!=\\!1.02/0.98$ in (b) and $c_{v}^{u}\\!=\\!(1.02/0.98)^{2}$ in
(g). We find from (b) and (f) that the increase ($\nearrow$) in $a$ and the
decrease ($\searrow$) in $b$ cause the same effect, enforcing the pattern
anisotropy along the $x$ direction, and from (c) and (e) that the decrease
($\searrow$) of $a$ and the increase ($\nearrow$) of $b$ cause the same
effect, enforcing the pattern along the $y$ direction. These effects are
consistent with the observation that the anisotropy of the pattern in (g) is
stronger than that in (b) and (f), where $c_{v}^{u}\\!=\\!(1.02/0.98)^{2}$ in
(g) is larger than $c_{v}^{u}\\!=\\!1.02/0.98$ in (b) and
$c_{v}^{u}\\!=\\!1.02/0.98$ in (f). Snapshots of $v$ are almost identical to
those of $u$ and are not plotted.

To quantify the observed anisotropy in the TPs caused by anisotropic diffusion
constants in the diffusion terms of Eq. (6), we calculate the mean values of
the absolute of the second-order partial differentials of $u$ and $v$

$\displaystyle\begin{split}&d^{2}_{x}u=\frac{1}{N}\sum_{ij}\left|(u_{i+1,j}+u_{i-1,j}-2u_{i,j}\right|,\quad
d^{2}_{y}u=\frac{1}{N}\sum_{ij}\left|u_{i,j+1}+u_{i,j-1}-2u_{i,j}\right|,\\\
&d^{2}_{x}v=\frac{1}{N}\sum_{ij}\left|v_{i+1,j}+v_{i-1,j}-2v_{i,j}\right|,\quad
d^{2}_{y}v=\frac{1}{N}\sum_{ij}\left|v_{i,j+1}+v_{i,j-1}-2v_{i,j}\right|,\end{split}$
(8)

which are “interaction” parts, while the factors $a$, $2\\!-\\!a$ and $b$,
$2\\!-\\!b$ in Eq. (6) are considered “coefficients”, from which $D_{u}$ and
$D_{v}$ are removed for simplicity. These interactions and coefficients are
closely connected to each other. Therefore, the interactions are expected to
be dependent on the parameters $a$ and $b$ and, as a consequence, on the input
anisotropy $c_{v}^{u}$. The reason the absolute symbol $|*|$ is used is that
the sign of the second-order partial differentials $\partial^{2}_{x,y}u$ and
$\partial^{2}_{x,y}v$ are locally changeable in their sign and cancel out when
they are summed over, and hence, we need $|*|$ in the sum to see how much
$\partial^{2}_{x,y}u$ and $\partial^{2}_{x,y}v$ deviate from zero.

Figure 4: Absolute second-order partial differentials $d^{2}_{\mu}u$ and
$d^{2}_{\mu}v$ obtained with (a) $a$ variation under the conditions
$a\\!\geq\\!b(\\!=\\!1)$ and (b) $b$ variation under $(1\\!=\\!)a\\!\leq\\!b$
on the $N\\!=\\!100^{2}$ lattice. The horizontal lower axes represent $a$ and
$b$. Axes $a$ and $b$ are more intuitive than $c_{v}^{u}$, which is also shown
on the upper axes. The results $d^{2}_{x}u\\!<\\!d^{2}_{y}u$ and
$d^{2}_{x}v\\!<\\!d^{2}_{y}v$ obtained under the conditions
$a\\!\geq\\!b(=\\!1)$ in (a) and those $d^{2}_{x}u\\!>\\!d^{2}_{y}u$ and
$d^{2}_{x}v\\!>\\!d^{2}_{y}v$ obtained under $(1\\!=\\!)a\\!\leq\\!b$ in (b)
are consistent with the pattern anisotropy in the $x$ and $y$ directions in
Fig. 3, respectively, as discussed in the text.

Figure 4(a) shows $d^{2}_{\mu}u$ and $d^{2}_{\mu}v$ ($\mu\\!=\\!x,y$)
calculated by varying $a$ under the conditions $b\\!=\\!1$ and $a\\!\geq\\!1$.
The other input parameters shown in the figure are the same as those assumed
in the snapshots in Fig. 3. We find from Fig. 4(a) that $d^{2}_{x}u$ and
$d^{2}_{x}v$ decrease ($\searrow$) and that $d^{2}_{y}u$ and $d^{2}_{y}v$
increase ($\nearrow$) with increasing $a$ ($\nearrow$). The results plotted in
Fig. 4(b) also show how $d^{2}_{\mu}u$ and $d^{2}_{\mu}v$ vary when $b$
increases ($\nearrow$) under the conditions $a\\!=\\!1$ and $b\\!\geq\\!1$.
These results are summarized as follows:

$\displaystyle\begin{split}{\rm
Fig.\ref{fig-4}(a)(b=1):}\;&d_{x}^{2}u\searrow\;({\rm
if}\;a\nearrow(>1)),\;\;d_{y}^{2}u\nearrow\;({\rm if}\;2-a\searrow(<1)),\\\
\left[\Leftrightarrow\;\right.&\left.d_{x}^{2}u\nearrow\;({\rm
if}\;a\searrow(<1)),\;\;d_{y}^{2}u\searrow\;({\rm
if}\;2-a\nearrow(>1)),\right]\\\
&d_{x}^{2}v\searrow\;(b=1),\;\;d_{y}^{2}v\nearrow\;(2-b=1),\end{split}$ (9)

and

$\displaystyle\begin{split}{\rm
Fig.\ref{fig-4}(b)(a=1):}\;&d_{x}^{2}u\nearrow\;(a=1),\;\;d_{y}^{2}u\searrow\;(2-a=1),\\\
&d_{x}^{2}v\nearrow\;({\rm if}\;b\nearrow(>1)),\;\;d_{y}^{2}v\searrow\;({\rm
if}\;2-b\searrow(<1)).\\\
\left[\Leftrightarrow\;\right.&\left.d_{x}^{2}v\searrow\;({\rm
if}\;b\searrow(<1)),\;\;d_{y}^{2}v\nearrow\;({\rm
if}\;2-b\nearrow(>1)).\right]\end{split}$ (10)

The statements in the second line of Eq. (9) and those in the third line of
Eq. (10) have no corresponding data in Fig. 4(a) and (b) and appear not to be
supported by the obtained data. However, this outcome is expected because
anisotropic patterns depend only on the ratio $a/b$ at least in the region
close to $a/b\\!=\\!1$ under the condition $a\\!=\\!1$ or $b\\!=\\!1$. For
this reason, we remove the data for $a\\!<\\!1,b\\!=\\!1$ and
$b\\!<\\!1,a\\!=\\!1$ from Fig. 4(a) and (b) to simplify the discussion.

The behaviors of $d^{2}_{\mu}u$ and $d^{2}_{\mu}v$ plotted in Fig. 4(a) and
(b) are expected from the snapshots in Fig. 3. Indeed, we expect
$d^{2}_{x}u\\!<\\!d^{2}_{y}u$ in the patterns that are anisotropic along the
$x$ direction, such as those in Fig. 3(b), (f) and (g). On the patterns
showing anisotropy along the $x$ direction, $u$ is almost constant, implying a
small $d^{2}_{x}u$. This behavior $d^{2}_{x}u\\!<\\!d^{2}_{y}u$ is also
expected from the fact that the direction-dependent diffusion constants
$\left(D_{u}^{x},D_{u}^{y}\right)\\!=\\!\left(aD_{u},(2-a)D_{u}\right)$ in Eq.
(3) are $D_{u}^{x}>D_{u}^{y}$ for $a\\!>\\!1$ and $b\\!=\\!1$ because a large
diffusion constant $D_{u}^{x}$ induces a long-distance spatial correlation of
$u$, implying almost constant $u$ along the $x$-axis, which is reflected in a
small $d^{2}_{x}u$. In fact, if the long-range correlation of $u$ is
accompanied by a large $d^{2}_{x}u$, which corresponds to a large value of the
integral of $(\partial_{x}u)^{2}$ over the domain, then the corresponding
energy $\int dxdy((\partial_{x}u)^{2}\\!+\\!(\partial_{y}u)^{2})$ becomes very
large. This contradicts the fact that Turing instability lies close to a
stable state, which is the energy minimum, of solutions in the RD equation.

On the other hand, $d^{2}_{x,y}v$ has a notably similar behavior to
$d^{2}_{x,y}u$ even under the isotropic constants
$D_{v}^{x}\\!=\\!D_{v}^{y}\\!=\\!D_{v}$ due to the condition $b\\!=\\!1$, as
shown in the third line of Eq. (9). This nontrivial behavior in $d^{2}_{x,y}v$
comes from the implemented interaction between $u$ and $v$ in the reaction
terms $f$ and $g$ in Eq. (LABEL:FN-eq-Eucl). We note that this interaction is
present without the diffusion anisotropy $(a,b)\\!=\\!(1,1)$, and therefore,
this interaction is a reason for the deviations
$d^{2}_{x}u\\!\not=\\!d^{2}_{y}u$ at $a\\!=\\!1$ and
$d^{2}_{x}v\\!\not=\\!d^{2}_{y}v$ at $b\\!=\\!1$ in Fig. 4(a) and (b). These
small deviations imply that the patterns are not completely isotropic but
slightly anisotropic and spontaneously generated, even though the snapshot in
Fig. 3(a) looks isotropic.

The result obtained under the conditions $a\\!>\\!1$ and $b\\!=\\!1$ indicates
that the patterns of $u$ anisotropy in the $x$ direction are caused by the
assumed diffusion anisotropy, and moreover, the anisotropic patterns in $v$
are caused by both the anisotropic pattern of $u$ and the interaction between
$u$ and $v$. The data $d^{2}_{\mu}u$ and $d^{2}_{\mu}v$ plotted in Fig. 4(b)
are calculated by varying $b$ with fixed $a\\!=\\!1$. In this case, the
condition $a\\!<\\!b$ is satisfied and opposite to $a\\!>\\!b$ in Fig. 4(a),
and hence, we obtain the results $d^{2}_{x}u\\!>\\!d^{2}_{y}u$ and
$d^{2}_{x}v\\!>\\!d^{2}_{y}v$, which are opposite to those obtained in Fig.
4(a), as summarized in Eq. (10). The result confirms that the anisotropic
patterns in $u$ are caused by the anisotropic patterns of $v$ and the
interaction between $u$ and $v$.

Figure 5: $S^{\mu}_{u},(\mu=x,y)$ and $S^{\mu}_{v},(\mu=x,y)$ obtained with
(a) $a$ variation under the conditions $a\\!\geq\\!b(\\!=\\!1)$ and (b) $b$
variation under $(1\\!=\\!)a\\!\leq\\!b$ on the $N\\!=\\!100^{2}$ lattice. The
assumed parameters are the same as those in Fig. 4.

Anisotropies in diffusion constants in Eq. (3) can also be reflected in the
quantities $\int(\frac{\partial u}{\partial\mu})^{2}dxdy,(\mu\\!=\\!x,y)$ and
$\int(\frac{\partial v}{\partial\mu})^{2}dxdy,(\mu\\!=\\!x,y)$. The discrete
expressions, denoted by $S^{\mu}_{u}$ and $S^{\mu}_{v}$, are given by

$\displaystyle\begin{split}&S^{x}_{u}=\frac{1}{4}\sum_{ij}\left(u_{i+1,j}-u_{i-1,j}\right)^{2},\quad
S^{y}_{u}=\frac{1}{4}\sum_{ij}\left(u_{i,j+1}-u_{i,j-1}\right)^{2},\\\
&S^{x}_{v}=\frac{1}{4}\sum_{ij}\left(v_{i+1,j}-v_{i-1,j}\right)^{2},\quad
S^{y}_{v}=\frac{1}{4}\sum_{ij}\left(v_{i,j+1}-v_{i,j-1}\right)^{2},\end{split}$
(11)

where differentials are replaced by differences; $\frac{\partial u}{\partial
x}\to\frac{1}{2\Delta x}\left(u_{i+1,j}-u_{i-1,j}\right)$ with $\Delta
x\\!=\\!1$. The behaviors of $S^{\mu}_{u}/N$ and $S^{\mu}_{v}/N$ plotted in
Fig. 5(a) and (b) are almost the same as those in Fig. 4(a) and (b).

##  III Finsler geometry modeling of Turing patterns

In the preceding section, we confirm that anisotropic TPs appear when the
diffusion constants are anisotropic in the DR equations under the condition of
suitable input parameters. However, as emphasized in the introduction, the
origin of such anisotropy in the diffusion constants is unclear. In this
section, which is the main part of this paper, we show one possible origin of
the anisotropic diffusion constants.

###  III.1 Internal degree of freedom on the triangulated lattice and the
Monte Carlo update

Figure 6:  (a) Triangulated square lattice of size $N\\!=\\!100$, which is the
total number of vertices given by $N\\!=\\!n_{x}n_{y}$, where the total number
of boundary vertices is $n_{\mu}\\!=\\!10(\mu\\!=\\!x,y)$. The edge length
$L_{\mu}(\mu\\!=\\!x,y)$ is given by $L_{\mu}\\!=\\!(n_{\mu}\\!-\\!1)d$, where
$d$ is the lattice spacing. PBCs are assumed on the boundary vertices. The
dashed square is the boundary frame of area $A\\!=\\!L_{x}L_{y}$ implemented
by the PBCs. (b) A lattice site $i$ and its neighboring sites denoted by
$j(i)$, $1\\!\leq\\!j(i)\\!\leq\\!J_{\rm max}$ with $J_{\rm max}\\!=\\!6$.
This $J_{\rm max}$ depends on the site $i$ and is also called the coordination
number denoted by $q_{i}$.

First, we present detailed information on triangulated lattices, in which the
vertex position plays a role in the IDOF for anisotropic diffusion. This new
IDOF is governed by the Hamiltonian $S$, as introduced in the following
subsection. A triangulated lattice is plotted in Fig. 6, where PBCs are
assumed for all vertex positions because their positions are allowed to move
as shown below. This assumption is slightly different from those associated
the regular square lattice in Fig. 2(a), where PBCs are assumed only on the
boundary vertices. On the initial configuration in Fig. 6(a), vertices
enclosed by the solid rectangles are identified with those enclosed by the
same-sized oblong dashed rectangles, which are plotted at the positions
$n_{\mu}d,(\mu\\!=\\!x,y)$ distant from the original positions, where
$n_{\mu}$ and $d$ stand for the total number of vertices on the lattice edges
and the lattice spacing suitably assumed in the simulations, respectively. The
dashed square represents an implemented frame by PBCs, in which the vertex
position $\vec{r}$ is identified with $\vec{r}\\!\pm\\!(n_{x}d,n_{y}d)$:

$\displaystyle\vec{r}\equiv\vec{r}\pm(n_{x}d,n_{y}d)\quad({\rm PBCs}),$ (12)

Consequently, the frame area is fixed to $A\\!=\\!L_{x}L_{y}$, where
$L_{\mu}(\mu\\!=\\!x,y)$ is the edge length given by
$L_{\mu}\\!=\\!(n_{\mu}\\!-\\!1)d$. A frame tension or surface tension thus
emerges on the surface, as discussed in the following section. Vertices on and
inside the dashed frame in Fig. 6(a) represent the initial configuration, and
those inside are randomly distributed under the constraint that the minimum
distance $r_{\rm min}\\!=\\!0.8d$ from the other vertices under the PBCs in
Eq. (12). This $r_{\rm min}\\!=\\!0.8d$ is assumed only in the initial lattice
construction. The vertices are linked with the Voronoi tessellation technique
Friedberg-Ren-NPB1984 .

Figure 7:  Illustrations of MC updates for the (a) vertex position and (b)
bond flip under the constraint of nonfolding triangles. Lattices in which only
the vertex position is updated are called fixed connectivity lattices, and
those in which both the vertex position and the bond flip are updated are
called fluid lattices. (c) Triangles are allowed to deform in no-folding
configurations ($\bigcirc$) and prohibited from folding ($\times$). (d) The
internal degree of freedom $\vec{\tau}(\in S^{1}/2)$ at $\vec{r}$ is defined
by $\vec{\tau}\\!=\\!{\Delta}\vec{r}/\|{\Delta}\vec{r}\|$ using the
deformation vector ${\Delta}\vec{r}\\!=\\!\vec{r}^{\prime}\\!-\\!\vec{r}$. (e)
Illustrations of vertex fluctuations corresponding to isotropic and
anisotropic configurations of $\vec{\tau}$.

Next, to introduce an IDOF, we briefly explain the MC update Metropolis-
JCP-1953 ; Landau-PRB1976  of the vertex position $\vec{r}(\in{\bf R}^{2})$,
which is illustrated in Fig. 7(a). The new position $\vec{r}^{\,\prime}$ is
accepted with probability ${\rm Min}[1,\exp(\delta S)]$, where $\delta
S\\!=\\!S(\vec{r}^{\,\prime})\\!-\\!S(\vec{r})$ is the energy change before
and after the vertex movement $\vec{r}\\!\to\\!\vec{r}^{\,\prime}$. The new
position $\vec{r}^{\,\prime}$ is randomly fixed in a small circle of radius
$R$ centered at $\vec{r}$, where the radius $R$ is fixed so that the
acceptance rate is approximately equal to $60\%-90\%$ under the constraint
that the triangles are not folded at any bonds (nonfolding and folding
triangles are illustrated in Fig. 7(c)). This MC process for the update of
vertex positions is the same as in MC studies for polymerized membranes
KANTOR-NELSON-PRA1987 ; Gompper-Kroll-PRA1992 ; KOIB-PRE-2005 , which is a
two-dimensional extension of the linear chain model for polymers Doi-
Edwards-1986 . Another MC process is the so-called bond flip procedure in Fig.
7(b) Ho-Baum-EPL1990 . In this process, the bond $ij$ connecting vertices $i$
and $j$ is removed, and vertices $k$ and $l$ are connected by a bond. This
bond flip is also accepted with probability ${\rm Min}[1,\exp(\delta S)]$ with
energy change $\delta S$ under the constraint for the nonfolding triangles.
The IDOF $\vec{\tau}$ is defined by
$\vec{\tau}\\!=\\!{\Delta}\vec{r}/\|{\Delta}\vec{r}\|(\in S^{1}/2)$ (Fig.
7(d)), which has values on the half circle $S^{1}/2$ due to the nonpolar
nature assumed on the variable. Figure 7(e) illustrates isotropic and
anisotropic configurations of $\vec{\tau}$.

Figure 8:  Snapshots of vertex configurations updated by the metropolis MC
procedure of lattice size $N\\!=\\!400$: (a) initial configuration without the
MC update, (b) fixed connectivity lattice obtained after $1\\!\times\\!10^{5}$
MCS for the vertex updates, and (c) fluid lattice obtained after
$1\\!\times\\!10^{5}$ MCS for the updates of vertex and bond flip. Vertices on
the boundaries of the fixed lattice in (b) fluctuate only locally, while those
of the fluid lattice in (c) diffuse over the lattice. The total number of
bonds $N_{B}$ and the total number of triangles $N_{T}$ are given by
$N_{B}\\!=\\!3N$ and $N_{T}\\!=\\!2N$ on triangulated lattices with PBCs
satisfying zero Euler characteristics $N\\!-\\!N_{B}\\!+\\!N_{T}\\!=\\!0$.

Finally, in this subsection, we show snapshots of vertex configurations
updated by MC processes, including the initial configuration (Fig. 8(a)). To
clearly show the vertices, we use a lattice of size $N\\!=\\!400$. The
snapshot denoted by “fixed” in Fig. 8(b) is obtained after
$1\\!\times\\!10^{5}$ MC sweeps (MCSs) for the update of vertex position
$\vec{r}_{i}$ only, where 1 MCS represents $N$ consecutive updates of
$\vec{r}_{i},(i\\!=\\!1,\cdots,N)$. Vertices outside the dashed square, which
is shown only in Fig. 8(a), are plotted at the opposite side due to the PBC.
We find that the vertices fluctuate only locally on the fixed connectivity
lattice. On the other hand, the vertices diffuse almost freely over the
lattice denoted by “fluid” in Fig. 8(c) Ho-Baum-EPL1990 ; Peliti-Leibler-
PRL1985 ; we find that the boundary vertices are randomly mixed on the fluid
lattice, where 1 MCS represents $N$ consecutive updates of vertex positions
followed by $N$ random updates of bond flips.

We emphasize that these two different types of triangulated lattices are not
simply different in the polygon shape from the regular square lattice in
Section II but rather have extra dynamic degrees of freedom, such as the
vertex position. To suitably treat this new IDOF, we introduce the so-called
Gaussian bond potential, usually assumed for a model of cell membranes
HELFRICH-1973 ; Peliti-Leibler-PRL1985 ; NELSON-SMMS2004 , in the following
subsection. Note also that IDOF $\vec{\tau}$ is connected to directions of the
position movements and hence is called “internal”, in contrast to the
variables $u(x,y)$ and $v(x,y)$, which are functions of the positions.

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

Now, we introduce the FG modeling technique for anisotropic TPs based on the
DR equations of the FN type in Eq. (LABEL:FN-eq-Eucl). FG modeling modifies
length scales for interactions to be direction dependent by using an IDOF that
is not included in the original system, and therefore, the FG model is an
extension of the original model owing to the new IDOF. As a consequence, the
interaction coefficients in the original system are dynamically direction
dependent. The original system in this paper is the DR equations in Eq.
(LABEL:FN-eq-Eucl), which contain the variables $u$ and $v$. The time
evolution of these variables is also iterated in the FG modeling by using the
DR equations for obtaining TPs, while the new IDOF $\vec{\tau}$ is updated by
the MC procedure because we have no differential equation for $\vec{\tau}$. To
perform the MC procedure for $\vec{\tau}$, we assume the following discrete
Hamiltonian $S$ composed of a linear combination of several terms such that

$\displaystyle\begin{split}&S(\vec{r};\vec{\tau})=\left\\{\begin{array}[]{@{\,}ll}S_{1}+U_{V}+D_{u}S_{u}+D_{v}S_{v}+\lambda
S_{\tau}+S_{\tau}^{F},\quad({\rm fixed})\\\ S_{1}+U_{V}+U_{\cal
T}+D_{u}S_{u}+D_{v}S_{v}+\lambda S_{\tau}+S_{\tau}^{F},\quad({\rm
fluid})\end{array}\right.,\end{split}$ (13)

where $S(\vec{r};\vec{\tau})$ denotes that the variables $\vec{r}$ and
$\vec{\tau}$ are not independent; $\vec{\tau}$ is defined by the fictitious
time evolution of $\vec{r}$ in the MC update as shown in Fig. 7(c). The
difference in $S$ between the fixed and fluid models in Eq. (13) is that a
potential $U_{\cal T}$ is included in $S$ for the fluid model, where the
symbol ${\cal T}$ denotes a triangulation, which is considered as a variable
only in the fluid model (Fig. 7(b)). The partition functions are written as

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

where $\prod_{i}d\vec{r}_{i}$ denotes $2N$-dimensional multiple integrations
on a domain $\cal{D}$ in ${\bf R}^{2}$ and where $\sum_{\mathcal{T}}$ in
$Z_{\rm flu}$ denotes the sum over all possible triangulations. These are
simulated by the MC updates in Fig. 7(a) and (b).

The terms on the right-hand side of $S$ are defined as follows:

$\displaystyle\begin{split}&S_{1}=\sum_{ij}\ell_{ij}^{2},\quad\ell_{ij}^{2}\\!=\\!(\vec{r}_{i}\\!-\\!\vec{r}_{j})^{2},\\\
&U_{V}=\sum_{ij}U_{ij},\quad
U_{ij}=\left\\{\begin{array}[]{@{\,}ll}0\quad(\ell_{\min}\leq\ell_{ij}\leq\ell_{\rm
max})\\\ \infty\quad({\rm otherwise})\end{array}\right.,\\\
&\quad(\ell_{\min}=0.01d\quad\ell_{\rm max}=3d,\quad d:{\rm latt.\;sp.}),\\\
&U_{\cal T}=\sum_{i}U_{i},\quad
U_{i}=\left\\{\begin{array}[]{@{\,}ll}0\quad(4\leq q_{i}\leq 9)\\\
\infty\quad({\rm otherwise})\end{array}\right.,\\\
&S_{u}=\sum_{ij}\gamma_{ij}^{u}\left(u_{j}-u_{i}\right)^{2},\quad\gamma_{ij}^{u}=\frac{1}{6}\left(\frac{\chi_{ij}^{u}}{\chi_{ik}^{u}}+\frac{\chi_{ji}^{u}}{\chi_{jk}^{u}}+\frac{\chi_{ij}^{u}}{\chi_{il}^{u}}+\frac{\chi_{ji}^{u}}{\chi_{jl}^{u}}\right),\\\
&S_{v}=\sum_{ij}\gamma_{ij}^{v}\left(v_{j}-v_{i}\right)^{2},\quad\gamma_{ij}^{v}=\frac{1}{6}\left(\frac{\chi_{ij}^{v}}{\chi_{ik}^{v}}+\frac{\chi_{ji}^{v}}{\chi_{jk}^{v}}+\frac{\chi_{ij}^{v}}{\chi_{il}^{v}}+\frac{\chi_{ji}^{v}}{\chi_{jl}^{v}}\right),\\\
&S_{\tau}=-\sum_{ij}\left(\vec{\tau}_{i}\cdot\vec{\tau}_{j}\right)^{2},\quad
S_{\tau}^{F}=-\sum_{i}\left(\vec{\tau}_{i}\cdot{\vec{F}}\right)^{2},\quad{\vec{F}}=(F_{x},F_{y})\end{split}$
(15)

The first term $S_{1}$ is the Gaussian bond potential, which is a spring
potential defined by the sum of bond length squares $\ell_{ij}^{2}$, and
$\sum_{ij}$ denotes the sum over bonds $ij$. The meaning of the inclusion of
$S_{1}$ in $S$ is that the domain $\cal{D}$, defined by a fixed frame of side
length $L_{x}$ and $L_{y}$ for PBCs in Fig. 6(a), can be regarded as a
membrane surface with an internal structure rather than a plane in ${\bf
R}^{2}$ KANTOR-NELSON-PRA1987 ; Gompper-Kroll-PRA1992 ; KOIB-PRE-2005 ; Ho-
Baum-EPL1990 ; HELFRICH-1973 ; Peliti-Leibler-PRL1985 ; NELSON-SMMS2004 .
Therefore, we consider that the frame is spanned by a membrane of area
$A\\!=\\!L_{x}L_{y}$ with surface tension $\sigma$, as mentioned in the
preceding subsection (see Appendix A for detailed information on $\sigma$).

The second term $U_{V}$ is a constraint potential that prohibits the bond
length $\ell_{ij}$ from being out of the range
$\ell_{\min}\leq\ell_{ij}\leq\ell_{\rm max}$ with $\ell_{\min}\\!=\\!0.01d$
and $\ell_{\rm max}\\!=\\!3d$. The lattice spacing $d$ is fixed to two
different values in the simulations, as discussed in Section IV, to check the
influence of the surface tension $\sigma$ (Appendix A). The vertex moves are
constrained inside $\cal{D}$ by the fixed frame of side lengths $L_{x}$ and
$L_{y}$ in Fig. 6(a) and by the condition for nonfolding triangles in Fig.
7(c). These constraints can also be written in constraint potentials; however,
we omit these potentials for simplicity. We note that $S_{1}/N\\!=\\!1$ is
satisfied in the case without any constraints for the bond length, such as the
constraint for the side lengths $L_{x}$ and $L_{y}$ (Fig. 6(a)) and the
constraint $U_{V}$, due to the scale-invariant property of the partition
function (see Appendix A for further detail).

The term $U_{\cal T}$ for the fluid model is a constraint potential that
enforces $4\leq q_{i}\leq 9$, where the coordination number $q_{i}$ is the
total number of bonds connected to vertex $i$ (Fig. 6(b)). The terms $S_{u}$
and $S_{v}$ are discrete versions of the continuous Hamiltonians

$\displaystyle S_{u}=\frac{1}{2}\int\sqrt{g^{u}}d^{2}x\,g^{u,ab}\frac{\partial
u}{\partial x^{a}}\frac{\partial u}{\partial x^{b}},\quad
S_{v}=\frac{1}{2}\int\sqrt{g^{v}}d^{2}x\,g^{v,ab}\frac{\partial v}{\partial
x^{a}}\frac{\partial v}{\partial x^{b}}.$ (16)

In these expressions, $g^{u,ab}$ and $g^{v,ab}$ are the inverse and $g^{u,v}$
is the determinant of the Finsler metrics (see Appendix B for further detail
on the Finsler metric)

$\displaystyle g^{u}_{ab}=\begin{pmatrix}1/(\chi_{ij}^{u})^{2}&0\\\
0&1/(\chi_{ik}^{u})^{2}\end{pmatrix},\quad
g^{v}_{ab}=\begin{pmatrix}1/(\chi_{ij}^{v})^{2}&0\\\
0&1/(\chi_{ik}^{v})^{2}\end{pmatrix},$ (17)

where $\chi_{ij}^{u}$ and $\chi_{ij}^{v}$ are given by

$\displaystyle\begin{split}&\chi_{ij}^{u}=|\vec{\tau}_{i}\cdot\vec{e}_{ij}|+\chi_{0},\\\
&\chi_{ij}^{v}=\sqrt{1-|\vec{\tau}_{i}\cdot\vec{e}_{ij}|^{2}}+\chi_{0}.\end{split}$
(18)

These continuous Hamiltonians in Eq. (16) correspond to the diffusion terms in
Eq. (LABEL:FN-eq-Eucl). The position- and direction-dependent coefficients
$\gamma_{ij}^{u}$ and $\gamma_{ij}^{v}$ in Eq. (15) are given in Appendix B.
$S_{\tau}$ represents the interaction energy for nearest neighbor pairs of
$\vec{\tau}_{i}$ and $\vec{\tau}_{j}$, which are assumed to be nonpolar, with
the interaction coefficient $\lambda$. One additional assumption is that some
external force $\vec{F}$ aligns $\vec{\tau}$ along the direction of $\vec{F}$.
This interaction is described by the final term $S_{\tau}^{F}$, where nonpolar
interaction is assumed. If the interaction is connected to electromagnetic
interactions, polar interactions can also be assumed in $S_{\tau}$ and
$S_{\tau}^{F}$.

Here, we comment on an implication implemented in the definitions of
$\chi_{ij}^{u}$ and $\chi_{ij}^{v}$ in Eq. (18). Let us suppose $\vec{e}_{ij}$
is parallel to the $x^{1}$-axis. Then, the definition of $\chi_{ij}^{u}$
implies that the Finsler length $ds_{u}\\!=\\!dx^{1}/\chi_{ij}^{u}$ (Appendix
B) for interactions in $u$ becomes short along the direction of $\vec{e}_{ij}$
when $\vec{\tau}_{i}$ is aligned along $\vec{e}_{ij}$, where $\chi_{ij}^{u}$
plays a role in the unit Finsler length, and hence, a large $\chi_{ij}^{u}$
means a short Finsler length. As a consequence, this alignment of
$\vec{\tau}_{i}$ along $\vec{e}_{ij}$ effectively makes diffusion coefficients
$D_{u}$ direction dependent, which is numerically confirmed in the following
section. This direction dependence of $D_{u}$ is observed in the model of
Section II when the parameter $a$ is increased in Eq. (3). In contrast, the
definition of $\chi_{ij}^{v}$ implies that an alignment of $\vec{\tau}_{i}$
along $\vec{e}_{ij}$ effectively makes the $D_{v}$ direction dependent,
corresponding to the anisotropic TPs in Fig. 3(f), (g) obtained when $b$ is
decreased. Thus, we expect that these $\chi_{ij}^{u}$ and $\chi_{ij}^{v}$
correspond to the condition $a\\!>\\!b(=\\!1)$ ($\Leftrightarrow
c_{v}^{u}\\!>\\!1$) or $(1\\!=)a\\!<\\!b$ ($\Leftrightarrow
c_{v}^{u}\\!<\\!1$) shown in Fig. 4(a), (b) if $\vec{\tau}_{i}$ is aligned
along $\vec{e}_{ij}$. Therefore, we expect that if the alignment of
$\vec{\tau}_{i}$ is controlled to be aligned along the $x$ direction for
almost all $i$ for example, then $\chi_{ij}^{u}$ and $\chi_{ij}^{v}$ in Eq.
(18) effectively plays a role in $(a,b)$ in Eq. (3).

The discrete Laplace operators on triangulated lattices are obtained by the
discrete $S_{u}$ and $S_{v}$ in Eq. (15) such that

$\displaystyle\begin{split}&\mathop{}\\!\mathbin{\bigtriangleup}u_{i}=2\left(\sum_{j(i)}\gamma_{ij}^{u}u_{j}-u_{i}\sum_{j(i)}\gamma_{ij}^{u}\right),\\\
&\mathop{}\\!\mathbin{\bigtriangleup}v_{i}=2\left(\sum_{j(i)}\gamma_{ij}^{v}v_{j}-u_{i}\sum_{j(i)}\gamma_{ij}^{v}\right)\end{split}$
(19)

for $1\\!\leq\\!i\\!\leq\\!N$, where $\sum_{j(i)}$ denotes the sum over
vertices $j$ connected to $i$ (Appendix C). Note that these Laplace operators
are independent of the anisotropy parameters $a$ and $b$ in contrast to those
in Eq. (6).

###  III.3 Hybrid numerical technique

We introduce a numerical technique to find steady state configurations of $u$
and $v$ under the presence of IDOF $\tau$ on two types of triangulated
lattices: fixed and fluid. The variables $u_{i}$ and $v_{i}$ are updated such
that

$\displaystyle\begin{split}&u_{i}(t+{\Delta}t)\leftarrow
u_{i}(t)+{\Delta}t\left(D_{u}{\mathop{}\\!\mathbin{\bigtriangleup}}u_{i}(t)+f(u_{i}(t),v_{i}(t)\right),\\\
&v_{i}(t+{\Delta}t)\leftarrow
v_{ij}(t)+{\Delta}t\left(D_{v}{\mathop{}\\!\mathbin{\bigtriangleup}}v_{i}(t)+g(u_{i}(t),v_{i}(t)\right),\end{split}$
(20)

where the diffusion terms are given by Eq. (19). Direction dependence is not
manually introduced in these diffusion terms, in sharp contrast to those of
Eq. (6). ${\rm Max}\left\\{|u_{i}(t\\!+\\!\Delta
t)\\!-\\!u_{i}(t)|\right\\}\\!<\\!1\\!\times\\!10^{-8}$ and ${\rm
Max}\left\\{|v_{i}(t\\!+\\!\Delta
t)\\!-\\!v_{i}(t)|\right\\}\\!<\\!1\\!\times\\!10^{-8}$,
$(1\\!\leq\\!i\leq\\!N)$, which are the same relations as those in Eq. (7) for
the standard model on regular square lattices.

The hybrid numerical technique is summarized in the following four steps:

  1. (i)

Discrete time evolution of Eq. (20) is iterated with the diffusion terms
$\mathop{}\\!\mathbin{\bigtriangleup}u_{i}$ and
$\mathop{}\\!\mathbin{\bigtriangleup}v_{i}$ in Eq. (19) with the coefficients
$D_{u}$ and $D_{v}$.

  2. (ii)

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

  3. (iii)

Steps (i) and (ii) are repeated $n_{\rm MC}$ times, where $n_{\rm MC}$ is
suitably large.

  4. (iv)

Steps (i) are repeated under the final configurations of
$\\{\vec{\tau}_{i}\\}$ in (iii) until the convergent criteria described above
are satisfied.

##  IV Results

As described in the preceding section, the computational domain is considered
to be spanned by a membrane with a surface tension $\sigma$, which depends on
the lattice spacing $d$ (Fig. 6(a)). For this reason, to check whether the
results are influenced by $\sigma$ or not, we use two different values of $d$
such that

$\displaystyle\begin{split}&d=0.525\quad(\Leftrightarrow\langle\ell^{2}\rangle\\!\simeq\\!\frac{1}{2}\Leftrightarrow\sigma>0),\\\
&d=0.41\quad(\Leftrightarrow\langle\ell^{2}\rangle\\!\simeq\\!\frac{1}{3}\Leftrightarrow\sigma\simeq
0).\end{split}$ (21)

The mean squared bond lengths $\langle\ell^{2}\rangle$ inside the parenthesis
are numerically obtained under the assumed $d$. Therefore, from the expression
of surface tension
$\sigma\\!=\\!({3N}/{A})\left(\langle\ell^{2}\rangle\\!-\\!\frac{1}{3}\right)$
in Eq. (31), we expect that $\sigma\\!>\\!0$ for $d\\!=\\!0.525$ and
$\sigma\\!\simeq\\!0$ for $d\\!=\\!0.41$. We simply expect that the results
are independent of $\sigma$ because $\sigma$ has no direction dependence.

###  IV.1 Snapshots

The final configurations of $u$ and $v$ depend on their initial configurations
and on $\tau$, and the obtained snapshots are almost identical but not always
exactly the same if the initial values of these variables are different from
each other. For this reason, steps (i) to (iv) described in Section III.3 are
iterated once in the simulations to obtain snapshots in this subsection.

Figure 9(a)–(d) show snapshots of the variable $u$ without the IDOF
$\vec{\tau}$ (upper row) and with $\vec{\tau}$ (lower row) of the fixed model
on the $N\\!=\\!10000$ lattice and the lattice spacing $d\\!=\\!0.525$. The
values of $u$ are normalized in the range $[0,1]$ in the snapshots, as shown
in the color code. These snapshots are obtained by varying $\lambda$ and
$\vec{F}$, both of which influence the alignment of $\vec{\tau}$; $\lambda$ is
the strength of the nearest neighbor correlation of the IDOF $\vec{\tau}$, and
$\vec{F}$ aligns $\vec{\tau}$ along $\vec{F}$. The other parameters are fixed
to $(\alpha,\gamma)\\!=\\!(1,8)$, and $(D_{u},D_{v})\\!=\\!(0.4,5)$. The
$\chi_{0}$ in Eq. (18) is assumed to be $\chi_{0}\\!=\\!0.5$. When
$\lambda\\!=\\!0$ and $\vec{F}\\!=\\!(0,0)$, the vertex positions $\vec{r}$
are expected to be random, causing a random fluctuation of IDOF $\vec{\tau}$.
In this case, an isotropic pattern is expected. These expectations are
confirmed in Fig. 9(a), where the variables $\vec{\tau}$ in the snapshots of
the central region are plotted with small cones to clarify their directions,
although $\vec{\tau}$ and $-\vec{\tau}$ are identified in $S_{\tau}$ and
$S_{\tau}^{F}$ of Eq. (15) and in $\nu_{ij}^{u}$ and $\nu_{ij}^{v}$ of Eq.
(18). When $\vec{F}$ is increased to $\vec{F}\\!=\\!(1,0)$ (Fig. 9(b)), the
pattern becomes aisotropic, and $\vec{\tau}$ is also slightly aligned along
the $x$ direction. The alignment of $\vec{\tau}$ becomes clear when $\lambda$
is increased to $\lambda\\!=\\!0.5$ (Fig. 9(c)), and the pattern is also
forced to be more anisotropic by enlarging $\vec{F}\\!=\\!(0,0)$ to
$\vec{F}\\!=\\!(0,2)$, where the direction is changed to the $y$-axis. From
these snapshots, we find that the isotropy/anisotropy in the patterns is
determined by the fluctuation direction $\tau$ of the vertices, which is
controlled by the external force $\vec{F}$. This finding implies that the FG
modeling prescription effectively modifies the diffusion constants $D_{u}$ and
$D_{v}$ in an anisotropic or direction-dependent manner. These effective
diffusion constants are numerically extracted from the simulation data in the
following subsection. Notably, the patterns of $v$ are almost the same as
those of $u$ plotted in Fig. 9.

Figure 9: Snapshots of TPs obtained on the fixed lattice of size
$N\\!=\\!10000$ with lattice spacing $d\\!=\\!0.525$. The assumed parameters
are $(\alpha,\gamma)\\!=\\!(1,8)$ and $(D_{u},D_{v})\\!=\\!(0.2,5)$, with (a)
$\lambda\\!=\\!0,\vec{F}\\!=\\!(0,0)$, (b)
$\lambda\\!=\\!0,\vec{F}\\!=\\!(1,0)$, (c)
$\lambda\\!=\\!0.5,\vec{F}\\!=\\!(1,0)$, and (d)
$\lambda\\!=\\!0.5,\vec{F}\\!=\\!(0,2)$. The upper rows show the normalized
$u(\in[0,1])$ indicated by the color code without the IDOF $\vec{\tau}$, and
the lower row shows enlarged views with $\vec{\tau}$ of the central region
enclosed by the dashed rectangles, in which $\vec{\tau}$ are represented by
small cones.

The patterns obtained with $d\\!=\\!0.41$ corresponding to
$\sigma\\!\simeq\\!0$ in Eq. (21) are almost the same as those in Fig. 9
obtained with $d\\!=\\!0.525$, and therefore, we find that the patterns are
not influenced by $\sigma$, as expected. The fact that the results are
independent of $\sigma$ is natural because $\sigma$ is isotropic (Fig. 19) and
does not change the lattice structure anisotropy. The pattern anisotropy is
caused by direction-dependent alignment of $\vec{\tau}$ reflecting anisotropy
in the vertex fluctuations. Influences of $\sigma$ on TPs are expected in the
case of anisotropic $\sigma$. This expectation is discussed in the following
subsection.

Figure 10: Snapshots of TPs obtained on the fluid lattice of size
$N\\!=\\!10000$ and lattice spacing $d\\!=\\!0.525$. The assumed parameters
are $(\alpha,\gamma)\\!=\\!(1,8)$ and $(D_{u},D_{v})\\!=\\!(0.2,10)$, with (a)
$\lambda\\!=\\!0,\vec{F}\\!=\\!(0,0)$, (b)
$\lambda\\!=\\!0,\vec{F}\\!=\\!(1,0)$, (c)
$\lambda\\!=\\!0.5,\vec{F}\\!=\\!(1,0)$, and (d)
$\lambda\\!=\\!0.5,\vec{F}\\!=\\!(0,2)$. These are the same as those assumed
in Fig. 9 except for $D_{v}\\!=\\!10$, which is twice as large as that of Fig.
9.

Snapshots obtained on the fluid lattice are shown in Fig. 10(a)–(d), where the
parameters are the same as those for the fixed model presented in Fig. 9
except $D_{v}\\!=\\!10$, which is $D_{v}\\!=\\!5$ in Fig. 9. If
$D_{v}\\!=\\!5$ is used for the fluid model, some of the patterns in Fig. 10
do not appear. This difference is considered to come from a difference in
lattice structure between fixed and fluid lattices. The vertex positions
$\vec{r}$ on fluid lattices are expected to be more influenced by the external
force $\vec{F}$ and their nearest neighbors via the correlation energy
$S_{\tau}$ due to the free diffusion of vertices shown in Fig. 8(c). To see
such a difference in the lattice structure, we calculate the mean value
distribution of the bond length $h(\ell)$ as well as the direction-dependent
diffusion constants in the following subsection.

###  IV.2 Direction-dependent diffusion constants and surface tension

In this subsection, we plot direction-dependent diffusion constants (Appendix
D)

$\displaystyle\begin{split}&D_{u}^{x}=\frac{1}{N_{B}}\sum_{ij}\gamma_{ij}^{u}\cos^{2}\theta_{ij},\quad
D_{u}^{y}=\frac{1}{N_{B}}\sum_{ij}\gamma_{ij}^{u}\sin^{2}\theta_{ij},\\\
&D_{v}^{x}=\frac{1}{N_{B}}\sum_{ij}\gamma_{ij}^{v}\cos^{2}\theta_{ij},\quad
D_{v}^{y}=\frac{1}{N_{B}}\sum_{ij}\gamma_{ij}^{v}\sin^{2}\theta_{ij},\end{split}$
(22)

and the corresponding direction-dependent energies, where $N_{B}(=\\!3N)$ is
the total number of bonds and $\theta_{ij}$ is the angle between
$\vec{e}_{ij}$ and the $x$-axis (Appendix D),

$\displaystyle\begin{split}&S_{u}^{x}=\frac{1}{\langle\gamma_{ij}^{u}\cos^{2}\theta_{ij}\rangle}\sum_{ij}\gamma_{ij}^{u}\cos^{2}\theta_{ij}(u_{i}-u_{j})^{2},\;S_{u}^{y}=\frac{1}{\langle\gamma_{ij}^{u}\sin^{2}\theta_{ij}\rangle}\sum_{ij}\gamma_{ij}^{u}\sin^{2}\theta_{ij}(u_{i}-u_{j})^{2},\\\
&S_{v}^{x}=\frac{1}{\langle\gamma_{ij}^{v}\cos^{2}\theta_{ij}\rangle}\sum_{ij}\gamma_{ij}^{v}\cos^{2}\theta_{ij}(v_{i}-v_{j})^{2},\;S_{u}^{y}=\frac{1}{\langle\gamma_{ij}^{v}\sin^{2}\theta_{ij}\rangle}\sum_{ij}\gamma_{ij}^{v}\sin^{2}\theta_{ij}(v_{i}-v_{j})^{2}.\end{split}$
(23)

We calculate these quantities using the final configurations obtained at step
(iv) described in Section III.3. The final configurations of $u$, $v$ and
$\tau$ depend on their initial configurations, and for this reason, we obtain
their mean final values by repeating steps (i) to (iv) $n_{\rm itr}$ with
random initial configurations. The $n_{\rm itr}$, the total number of
iterations $n_{\rm MC}$ of steps (i) and (ii), and the lattice size $N$ for
the simulations in this subsection are assumed as follows:

$\displaystyle n_{\rm itr}=50,\quad n_{\rm MC}=5\times 10^{5},\quad N=6400.$
(24)

Figure 11(a) shows $D_{u}^{\mu}$ vs. $F$ and $S_{u}^{\mu}/N_{B}$ vs. $F$,
where $\vec{F}\\!=\\!(F,0)$ and $N_{B}$ is the total number of bonds. Open and
solid symbols correspond to the data obtained under $d\\!=\\!0.525$ and
$d\\!=\\!0.41$, respectively, and the data are almost independent of the
difference in the lattice spacing $d$. We also find that $D_{u}^{x}$
($D_{u}^{y}$) increases (decreases) with increasing $F$. Note that the IDOF
$\vec{\tau}$ is almost random at $F\to 0$ and aligns along the $x$-axis when
$F$ increases. This behavior of $D_{u}^{x}$ and $D_{u}^{y}$ in Fig. 11(a) is
consistent with those $D_{u}^{x}\\!=\\!aD_{u}$ and
$D_{u}^{y}\\!=\\!(2-a)D_{u}$ expected when $a$ increases from $a\\!=\\!1$ in
Eq. (3) in the preceding section. Moreover, we also find that
$S_{u}^{x}/N_{B}$ ($S_{u}^{y}/N_{B}$) decreases (increases) with increasing
$F$. These results indicate that $S_{u}^{x}/N_{B}$ ($S_{u}^{y}/N_{B}$)
decreases (increases) when $D_{u}^{x}$ increases ($D_{u}^{y}$ decreases).
Here, we note that $S_{u}^{\mu}/N_{B}$ and $D_{u}^{\mu}$ have no relationship
between input and output, but both $S_{u}^{\mu}/N_{B}$ and $D_{u}^{\mu}$ are
outputs to the input $F$. $D_{v}^{\mu}$ vs. $F$ and $S_{v}^{\mu}/N_{B}$ vs.
$F$ are plotted in Fig. 11(b), in which we also find that the variation in
$S_{v}^{\mu}/N_{B}$ is closely related with that of $D_{v}^{\mu}$.

Figure 11: (a) Diffusion constants $D_{u}^{\mu}$ in Eq. (22) and the
direction-dependent energies $S_{u}^{\mu}$ vs. $F$, and (b) $D_{v}^{\mu}$ and
$S_{v}^{\mu}$ vs. $F$ for the fixed model on a lattice of size $N\\!=\\!6400$.
The data obtained under the lattice spacing $d\\!=\\!0.525$ ($d\\!=\\!0.41$)
are denoted by the open (solid) symbols. The results are almost independent of
$d$. The total number of iterations is set to $n_{\rm itr}\\!=\\!50$ at each
$F$ to calculate the mean values of these physical quantities.

To compare the results plotted in Fig. 11(a), (b) with those in Fig. 4(a),
(b), we consider the results in Fig. 11(a) and (b) as the responses
$S_{u}^{\mu}$ and $S_{v}^{\mu}$ to the inputs $D_{u}^{\mu}$ and $D_{v}^{\mu}$
as follows:

$\displaystyle{\rm Fig.\ref{fig-11}(a):}\;S_{u}^{x}\searrow\;({\rm
if}\;D_{u}^{x}\nearrow),\quad S_{u}^{y}\nearrow\;({\rm
if}\;D_{u}^{y}\searrow),$ (25) $\displaystyle{\rm
Fig.\ref{fig-11}(b):}\;S_{v}^{x}\searrow\;({\rm if}\;D_{v}^{x}\searrow),\quad
S_{v}^{y}\nearrow\;({\rm if}\;D_{v}^{y}\nearrow).$ (26)

Thus, we find that the results in Eqs. (25) and (26) are consistent with the
behaviors of $d_{\mu}^{2}u$ and $d_{\mu}^{2}v$ vs. $a$ and $b$ in Fig. 4(a)
and (b). Indeed, the descriptions in Eq. (25) are consistent with the first
line in Eq. (9), and those in Eq. (26) are consistent with the third line in
Eq. (10), although $S_{u}^{\mu}/N_{B}$ ($S_{v}^{\mu}/N_{B}$) does not directly
correspond to $d_{\mu}^{2}u$ ($d_{\mu}^{2}v$). However, it is natural to
expect that variation of $D_{u}^{\mu}$ ($D_{v}^{\mu}$) included in
$S_{u}\\!=\\!D_{u}^{x}S_{u}^{x}\\!+\\!D_{u}^{x}S_{u}^{y}$
($S_{v}\\!=\\!D_{v}^{x}S_{v}^{x}\\!+\\!D_{v}^{x}S_{v}^{y}$) makes the
remaining part $S_{u}^{\mu}$ ($S_{v}^{\mu}$) vary oppositely in MC updates
owing to the expectation that $S_{u}$ ($S_{v}$) is almost at equilibrium and
hence almost stable. For the same reason, the coefficients $a$, $2\\!-\\!a$
and $b$, $2\\!-\\!b$ in the diffusion terms
${\mathop{}\\!\mathbin{\bigtriangleup}}_{a}u_{ij}$ and
${\mathop{}\\!\mathbin{\bigtriangleup}}_{b}v_{ij}$ in Eq. (6) influence the
corresponding second-order partial differentials $\partial^{2}_{x,y}u$ and
$\partial^{2}_{x,y}v$, variations of which are reflected in $d_{\mu}^{2}u$ and
$d_{\mu}^{2}v$. This consistency between Eqs. (25), (26) and Eqs. (9), (10)
implies that the FG modeling with the IDOF $\tau$ suitably implements the
diffusion anisotropy in TPs. The same consistency is obtained with the results
on the fluid lattice plotted in Fig. 12(a) and (b).

Figure 12: (a) Diffusion constants $D_{u}^{\mu}$ in Eq. (22) and direction-
dependent energies $S_{u}^{\mu}$ vs. $F$ and (b) $D_{v}^{\mu}$ and
$S_{v}^{\mu}$ vs. $F$ for the fluid model on a lattice of size $N\\!=\\!6400$.
The data obtained under the lattice spacing $d\\!=\\!0.525$ ($d\\!=\\!0.41$)
are denoted by the open (solid) symbols, indicating no dependence of the
results on $d$.  Figure 13: Surface tension $\sigma$ vs. $F$
($\vec{F}\\!=\\!(F,0)$) on the (a) fixed lattice with (b)–(e) snapshots of $u$
at $F\\!=\\!0$ and $F\\!=\\!2.8$ and on the (f) fluid lattice with (g)–(j)
snapshots of $u$ at $F\\!=\\!0$ and $F\\!=\\!2.8$, where the lattice spacing
is assumed to $d\\!=\\!0.525$ and $d\\!=\\!0.41$. The lattice size
($N\\!=\\!6400$) and the parameters are the same as those assumed in Figs. 11
and 12. We find that $\sigma(d\\!=\\!0.525)>0$ and $\sigma(d\\!=\\!0.41)\simeq
0$, and both of them decrease with increasing $F$.

Next, we show the surface tension $\sigma$ in Eq. (31) for the lattice spacing
$d\\!=\\!0.525$ and $d\\!=\\!0.41$ in Fig. 13(a) and (b). We have checked in
the preceding subsections that the patterns are only minimally influenced by
$\sigma$; however, it is still interesting to see interaction between $\sigma$
and IDOF $\tau$, which is controlled by the external force $\vec{F}$. We find
from Fig. 13(a),(b) that $\sigma(d\\!=\\!0.525)>0$ and
$\sigma(d\\!=\\!0.41)\simeq 0$ as expected, and moreover, $\sigma$ slightly
decreases as $F$ increases in both fixed and fluid lattices. A decrease in
$\sigma$ corresponds to a decrease in the mean bond length-squares
$\langle\ell^{2}\rangle$ for $\sigma(d\\!=\\!0.525)$ and
$\sigma(d\\!=\\!0.41)$ because of the expression of $\sigma\propto
d^{-2}(\langle\ell^{2}\rangle\\!-\\!{1}/{3})$ in Eq. (31) (see Appendix E for
the distribution of $\ell^{2}$ on both fixed and fluid lattices.) This
decrease in $\langle\ell^{2}\rangle$ indicates that the vertex distribution
becomes uniform in the large-$F$ region compared with a nonuniform
distribution, where a larger $\langle\ell^{2}\rangle$ is expected because the
vertices are confined in the fixed area
$A\\!=\\!(n_{x}\\!-\\!1)(n_{y}\\!-\\!1)d^{2}$, where
$n_{x}\\!=\\!n_{y}\\!=\\!80,d\\!=\\!0.525$. This dependence of $\sigma$ on $F$
implies that $\sigma$ depends on IDOF $\tau$ because $F$ directly influences
$\tau$, and therefore, we consider that $\tau$ interacts with $\sigma$. This
interaction between $\tau$ and $\sigma$, which is controlled by a frame or
boundary conditions, allows us to control $\tau$ with $\sigma$. As mentioned
in the first part of this section, the surface tension $\sigma$ imposed by
PBCs in both $x$ and $y$ directions are isotropic and does not deform the
lattice except for the uniform or isotropic expansion/compression, and hence,
a nontrivial effect is not expected in the TPs. However, if $\sigma$ is
anisotropic, $\sigma$ influences the IDOF $\tau$. Such an anisotropic $\sigma$
is expected to be caused by direction-dependent boundary conditions.
Therefore, $\vec{\tau}$ can be controlled by mechanical boundary conditions as
well as the external force $\vec{F}$. This prediction is confirmed in the
following subsection.

Notably, $\sigma$ decreases as the TPs change from isotropic to anisotropic. A
small $\sigma$ comes from the small tensile energy $S_{1}$ in Eq. (15),
implying that anisotropic TPs are mechanically stable. Anisotropic TPs are
caused by an alignment of IDOF $\vec{\tau}$ corresponding to lower energy
configurations with respect to $S_{\tau}$ in Eq. (15). Thus, we find that the
alignment of $\vec{\tau}$ lowers the mechanical energy $S_{1}$ as well as
$S_{\tau}$. This interaction between $\tau$ and $\sigma$ implies that
anisotropy in the fluctuation directions of $\vec{r}$ influences the position
of $\vec{r}$ because $\tau$ represents the fluctuation directions of
$\vec{r}$. Therefore, FG modeling based on IDOF $\vec{\tau}$ is considered
reasonable.

###  IV.3 Control of the pattern direction

Figure 14: Snapshots of TPs, of which the direction is (a) “forced” by
external $\vec{F}$ and (b) spontaneously determined without $\vec{F}$. The
central domains are plotted with small cones representing IDOF $\vec{\tau}$ in
the lower part; (c) “forcing” with the replacement of
$\chi_{ij}^{u}\leftrightarrow\chi_{ij}^{v}$ in Eq. (33), where the pattern
direction is vertical to IDOF $\vec{\tau}$.

First, we show that the pattern direction can be arbitrarily and spontaneously
determined. Snapshots in Fig. 14(a) denoted by “forced” are obtained with
$\vec{F}\\!=\\!(3,3)$ and $\lambda\\!=\\!0$. In this case, the direction of
IDOF $\tau$ is controlled by $\vec{F}$, and consequently, patterns align along
$\tau$. On the other hand, snapshots in Fig. 14(b) denoted by “spontaneous”
are obtained under $\vec{F}\\!=\\!(0,0)$ and $\lambda\\!=\\!3$. In this case,
IDOF $\vec{\tau}$ aligns to a spontaneously determined direction due to a
relatively large $\lambda(=\\!3)$; pattern directions depend on random numbers
for initial random configurations of $\vec{\tau}$, for instance. These
properties on pattern directions, forced and spontaneous, are specific to FG
models. Figure 14(c) is obtained by “forcing” with $\vec{F}\\!=\\!(3,3)$, the
same as in Fig. 14(a); however, the pattern direction is almost vertical to
that in Fig. 14(a). The reason is that $\chi_{ij}^{u}$ in Eq. (33) is replaced
by $\chi_{ij}^{v}$ such that
$\chi_{ij}^{v}\\!=\\!|\vec{\tau}_{i}\cdot\vec{e}_{ij}|\\!+\\!\chi_{0}$ and
$\chi_{ij}^{u}\\!=\\!\sqrt{1\\!-\\!|\vec{\tau}_{i}\cdot\vec{e}_{ij}|^{2}}\\!+\\!\chi_{0}$.
In this case, $\vec{\tau}$ direction is perpendicular to the pattern direction
as we confirm from the snapshot in the lower part of Fig. 14(c).

Figure 15: Initial configurations of deformed and original lattices.
Deformation does not change the lattice size $N\\!=\\!100$, which is very
small for visualization of the lattice structure. The ratio of side lengths is
(a) $R_{XY}\\!=\\!1.3$, (b) $R_{XY}\\!=\\!1$ and (ca) $R_{XY}\\!=\\!0.7$. The
area $L_{x}L_{Y}$ also remains unchanged under deformation.

Now, we deform the side lengths $L_{x}$ and $L_{y}(=\\!L_{x})$ of the lattice
using the parameter $R_{XY}$ in $L_{x}^{\prime}\\!=\\!\sqrt{R_{XY}}L_{x}$,
$L_{y}^{\prime}\\!=\\!L_{y}/\sqrt{R_{XY}}$ so that
$R_{XY}\\!=\\!L_{x}^{\prime}/L_{y}^{\prime}$ (Fig. 15(a)-(c)). Both the area
$L_{x}L_{y}$ and the size $N$ remain unchanged under deformation.

Figure 16: Alignment of anisotropic TPs along the longer direction of lattices
deformed with ratio $R_{XY}\\!=\\!L_{x}/L_{y}$ ranging from (a)
$R_{XY}\\!=\\!1.4$ to (h) $R_{XY}\\!=\\!0.6$. These alignments are caused by
strains applied on the boundary frame of lattices.

Snapshots on deformed lattices with the range
$1.4\\!\geq\\!R_{XY}\\!\geq\\!0.6$ of the fixed model are shown in Fig.
16(a)-(h). The TP direction aligns along the longer direction if the ratio
$R_{XY}$ deviates from $R_{XY}\\!=\\!1$ to a certain extent. The parameters
are shown in the figure. The coefficients $\lambda$ and $\vec{F}$ are fixed to
$\lambda\\!=\\!2$ and $\vec{F}\\!=\\!(0,0)$, and therefore, these alignments
emerge because the spontaneous direction is determined uniquely by the lattice
deformation. As shown in Fig. 15(a), (c), the bond lengths along the longer
direction are longer than those along the shorter direction. For this reason,
vertices move along the longer direction relatively easily compared to the
vertical direction, and therefore, IDOF $\vec{\tau}$ aligns along the longer
direction. This is the alignment mechanism of anisotropic TPs due to the
boundary condition. We note that the lattice deformation is caused by uniaxial
tensile strains or compressions, and therefore, we consider that mechanical
strains applied on the boundary impart TP anisotropy.

Figure 17:  (a) PBC is assumed in (a) the $x$ direction for vertices in red
color ($\bullet$) and blue color ($\bullet$) denoted by “Fixed boundary for
one-dim moves”, where the vertices can move along the boundary and the
boundary bonds are not flipped (Fig. 7(b)). The vertices ($\bullet$) outside
the fixed boundaries are “virtual points”, which are necessary to perform
Voronoi tessellation for triangulation and to define $\gamma_{ij}^{u,v}$ in
Eq. (15) on the boundaries. (b) PBC is assumed in the $y$ direction, where the
directions of the fixed boundary and virtual points are opposite to those in
(a). The lattice size is $N\\!=\\!100$, which is the sum of vertices
($\bullet$) and ($\bullet$). On the virtual vertices, the dynamical variables
are not defined.

On fluid lattices, this mechanism does not work because vertices freely move
due to their fluid nature. Nevertheless, even on fluid lattices, vertices
fluctuate along a direction in which PBCs are assumed, while the fixed
boundary condition is assumed to prohibit vertices from undergoing free
diffusion in the perpendicular direction. In this case, the spontaneous
alignment direction of TPs is also expected to be determined by the boundary
condition (Figs. 17(a),(b)). In this case, the spontaneous alignment direction
of TPs is expected to be determined also by the boundary condition.

Figure 18:  Snapshots obtained on fluid lattices, in which a pair of fixed
boundaries is assumed in (a) the $x$ direction and (b) the $y$ direction. Note
that the boundary conditions assumed here and shown in Figs. 17(a),(b) are
different from those used so far. The lattice size is $N\\!=\\!10000$, and the
parameters are shown in the figure. The alignment of TPs is a spontaneous one
because of $\vec{F}\\!=\\!(0,0)$, and therefore, spontaneously emergent
directions are determined by boundary conditions.

We find from Figs. 18(a),(b) that TPs align along the fixed-boundary
direction.

##  V Concluding remarks

In this paper, we have successfully shown that the Finsler geometry (FG)
modeling technique is applicable to a differential equation model. We
numerically studied anisotropic Turing patterns (TPs) by using the FG modeling
technique with an internal degree of freedom (IDOF) to implement diffusion
anisotropy in the reaction-diffusion (RD) equation of FitzHugh-Nagumo (FN). In
the ordinary numerical techniques used to solve the RD equation, direction-
dependent diffusion coefficients are assumed as an input to reproduce
anisotropic TPs. On the other hand, anisotropy in diffusion coefficients is
dynamically generated in FG modeling. In this sense, the anisotropic TPs are
attributed to the IDOF alignment originating in direction-dependent
fluctuations of vertices. The vertices are biologically interpreted as cells
or lumps of cells, and therefore, our results indicate that one possible
origin of anisotropic TPs is direction-dependent movements of cells compatible
with those observed in developmental processes.

The surface tension due to the boundary frame of membranes depends on the
frame area and interacts with the IDOF. For this reason, the IDOF on fixed-
connectivity lattices is expected to be controlled by uniaxial strains, which
can also be controlled by suitable boundary conditions. We confirm that the TP
direction is controlled by the surface boundary conditions such that the
direction aligns along the tensile or compressive strain direction on fixed-
connectivity lattices. On fluid lattices with a pair of fixed boundaries, we
confirm that anisotropic TPs spontaneously emerge along the direction of fixed
boundaries.

These anisotropic TPs are understood to be anisotropic diffusion in the FG
modeling, and the diffusion rate is calculable on these fluid lattices by
regarding Monte Carlo iteration as time. It is interesting to study the
dependence of anisotropic TPs on the diffusion rate. This remains to be a
future study.

###### Acknowledgements.

This work is supported in part by Collaborative Research Project J20Ly19 of
the Institute of Fluid Science (IFS), Tohoku University. The authors
acknowledge Dr. Jean-Paul Rieu and Dr. B. Ducharne 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.

##  Appendix A Computational domain spanned by a membrane with a nontrivial
surface tension

In this Appendix, we give detailed information on the computational domain
${\cal D}$ bounded by a fixed frame for PBCs in Figs. 6(a) and 8(a). As
mentioned in Section III, this domain ${\cal D}$ has an internal structure
described by the vertex positions $\vec{r}_{i}(\in{\bf
R}^{2},1\\!\leq\\!i\\!\leq\\!N)$, and these $\vec{r}_{i}$ are regarded as
membranes KANTOR-NELSON-PRA1987 ; Gompper-Kroll-PRA1992 ; KOIB-PRE-2005 ; Ho-
Baum-EPL1990 ; HELFRICH-1973 ; Peliti-Leibler-PRL1985 ; NELSON-SMMS2004 ,
which are classical mechanical $N$-particle systems governed by the spring
potential or the Gaussian bond potential $S_{1}(\vec{r})$ and by several
constraints imposed on $\vec{r}_{i}$. Mainly due to $S_{1}$ and the boundary
frame, the membrane ${\cal D}$ is exposed to a tensile stress $\sigma$, which
comes from the so-called scale-invariant property of the partition function,
as mentioned in Section III Wheater-JPA1994 . This property is common to both
fixed and fluid models, and hence, we use $Z_{\rm fix}$ in Eq. (14) here in
this Appendix for simplicity.

Figure 19:  (a) Surface enclosed by a fixed boundary frame of area $A$ under a
constant surface tension $\sigma$. This surface has a mechanical free energy
$F(A)\\!=\\!\sigma A(=\\!\int^{A}\sigma dA)$ because $\sigma$ is constant and
isotropic. (b) Illustrations of the original surface with a boundary frame and
an expanded surface by a scale parameter $\beta(\geq\\!1)$.

First, we note that the variables $\vec{r}_{i}$ are integrated out in $Z_{\rm
fix}$ such that $\int\prod_{i}d\vec{r}_{i}\exp(-S)$, corresponding to
$2N$-dimensional multiple integration. Here, we rewrite $S$ as
$S\\!=\\!S_{1}\\!+\\!U_{\rm fr}(A)\\!+\\!S_{\rm other}$, where $U_{\rm fr}(A)$
is the potential for fixing the frame area $A$ to $A\\!=\\!L_{x}L_{y}$ and not
included in $S$ of Eq. (13), as mentioned in the text, and $S_{\rm other}$ is
the remaining term in $S$. Note that $S_{1}\\!=\\!\sum_{ij}\ell_{ij}^{2}$ and
$U_{\rm fr}(A)$ depend on $\vec{r}$ and that $S_{\rm other}$ is independent of
$\vec{r}$. Now, we replace $\vec{r}$ with $\vec{r}_{\beta}\\!=\\!\beta\vec{r}$
in $Z_{\rm fix}$ with the scale parameter $\beta(\geq 1)$. This replacement
constitutes a simple variable change in the integrations $\int
d\vec{r}_{i}\\!\to\\!\beta^{2}\int d\vec{r}_{i}$,
$S_{1}\\!\to\\!\beta^{2}S_{1}$, $\cdots$, and it does not change $Z_{\rm
fix}$, and therefore, we have $Z_{\rm fix}(\vec{r})\\!=\\!Z_{\rm
fix}(\vec{r}_{\beta})$, where

$\displaystyle Z_{\rm
fix}(\vec{r}_{\beta})=\beta^{2N}\int\prod_{i}d\vec{r}_{i}\exp\left[-\beta^{2}S_{1}\\!-\\!U_{\rm
fr}(A_{\beta})\\!-\\!S_{\rm other}\right].$ (27)

In the expression $U_{\rm fr}(A_{\beta})$, $A_{\beta}$ is given by
$A_{\beta}\\!=\\!\beta^{-2}A$ because the fixed frame means that the area
scales as $A\\!\to\\!A_{\beta}\\!=\\!\beta^{-2}A$ according to the scale
change $\vec{r}\\!\to\\!\beta\vec{r}$. By differentiating both sides of
$Z_{\rm fix}(\vec{r})\\!=\\!Z_{\rm fix}(\vec{r}_{\beta})$ with respect to
$\beta$ and by fixing $\beta\\!\to\\!1$, we have

$\displaystyle\begin{split}0&=\left.\frac{\partial}{\partial\beta}Z_{\rm
fix}(\vec{r}_{\beta})\right|_{\beta=1}\\\
&=\left.\frac{\partial}{\partial\beta}\left(\beta^{2N}\int\prod_{i}d\vec{r}_{i}\exp\left[-\beta^{2}S_{1}\\!-\\!U_{\rm
fr}(A_{\beta})\\!-\\!S_{\rm other}\right]\right)\right|_{\beta=1}\\\
&=\left.2N\beta^{-1}Z_{\rm
fix}(\vec{r}_{\beta})\right|_{\beta=1}+\left.2\beta^{2N+1}\int\prod_{i}d\vec{r}_{i}\,S_{1}\,\exp\left[-S(\vec{r}_{\beta})\right]\right|_{\beta=1}+\left.\frac{\partial
Z_{\rm fix}(\vec{r}_{\beta})}{\partial A_{\beta}}\frac{\partial
A_{\beta}}{\partial\beta}\right|_{\beta=1}\\\ &=2NZ_{\rm
fix}(\vec{r})-2\int\prod_{i}d\vec{r}_{i}\,S_{1}\,\exp\left[-S(\vec{r})\right]-2A\frac{\partial
Z_{\rm fix}(\vec{r})}{\partial A},\end{split}$ (28)

where $\left.\frac{\partial
A_{\beta}}{\partial\beta}\right|_{\beta=1}\\!=\\!\left.\frac{\partial\beta^{-2}A}{\partial\beta}\right|_{\beta=1}\\!=\\!-2A$
is used in the final term of the third line. We note that the boundary of the
expanded surface by $\beta(\geq\\!1)$ is identified with the fixed frame only
when $\beta\\!\to\\!1$ (Fig. 19(b)). Multiplying $2^{-1}Z_{\rm
fix}(\vec{r})^{-1}$ by the final expression in Eq. (28), we obtain

$\displaystyle 0=N-\langle S_{1}\rangle-\frac{A}{Z_{\rm
fix}(\vec{r})}\frac{\partial Z_{\rm fix}(\vec{r})}{\partial A}.$ (29)

Since $\langle S_{1}\rangle\\!=\\!\langle\sum_{ij}\ell_{ij}^{2}\rangle$ can be
replaced by
$\sum_{ij}\langle\ell_{ij}^{2}\rangle\\!=\\!\sum_{ij}\langle\ell^{2}\rangle\\!=\\!\langle\ell^{2}\rangle\sum_{ij}1\\!=\\!N_{B}\langle\ell^{2}\rangle\\!=\\!3N\langle\ell^{2}\rangle$,
we have the mean squared bond length $\langle\ell^{2}\rangle$ such that

$\displaystyle\langle\ell^{2}\rangle=\frac{1}{3}-\frac{1}{3N}\frac{A}{Z_{\rm
fix}(\vec{r})}\frac{\partial Z_{\rm fix}(\vec{r})}{\partial A}.$ (30)

If the boundary frame for PBCs is not assumed, the second term on the right-
hand side is unnecessary, and hence, we have
$\langle\ell^{2}\rangle\\!=\\!{1}/{3}(\Leftrightarrow\langle
S_{1}\rangle/N\\!=\\!1)$, as mentioned in the main text.

The problem is how to evaluate the second term in the case that the frame is
present. The potential $U_{\rm fr}(A)$ can be expressed only as $U_{V}$ in Eq.
(15) and is therefore not differentiable. To evaluate $\frac{\partial Z_{\rm
fix}(\vec{r})}{\partial A}$ for surfaces with a fixed frame, we regard the
membrane of the $N$-particle system bounded by the frame as a sheet of area
$A$ without the internal structure associated with the $N$-particles. In this
simple surface of area $A$, its free energy $F(A)$ is given by
$F(A)\\!=\\!\sigma A$, where $\sigma$ is the surface tension. Therefore, its
partition function can also be expressed as $Z(A)\\!=\\!\exp[-F(A)]$. Thus,
${\partial Z_{\rm fix}(\vec{r})}/{\partial A}$ can be evaluated by replacing
$Z_{\rm fix}(\vec{r})$ with $\exp[-F(A)]\\!=\\!\exp(-\sigma A)$, and we have
$\langle\ell^{2}\rangle\\!=\\!1/3\\!+\\!\sigma A/(3N)$ from Eq. (30), and
therefore

$\displaystyle\sigma=\frac{3N}{A}\left(\langle\ell^{2}\rangle-\frac{1}{3}\right)\simeq\frac{3}{d^{2}}\left(\langle\ell^{2}\rangle-\frac{1}{3}\right),\quad
A=L_{x}L_{y}=(n_{x}-1)(n_{y}-1)d^{2},$ (31)

where $L_{x},L_{y}$ and $n_{x},n_{y}$ are shown in Fig. 6(a). Thus, $\sigma$
depends on the lattice spacing $d$.

We note that the scale-dependent constraint $U_{V}$ in Eq. (15) should be
included in the $\exp$ term of Eq. (27) such that
$\exp\left[-\beta^{2}S_{1}\\!-\\!U_{\rm fr}(A_{\beta})\\!-\\!U_{V}(\ell_{\rm
min}^{\beta},\ell_{\rm max}^{\beta})\\!-\\!S_{\rm other}\right]$ with
$\ell_{\rm min}^{\beta}\\!=\\!\beta^{-1}\ell_{\rm min}$ and $\ell_{\rm
max}^{\beta}\\!=\\!\beta^{-1}\ell_{\rm max}$. Correspondingly, the terms
$-\ell_{\rm min}\frac{\partial Z_{\rm fix}(\vec{r})}{\partial\ell_{\rm
min}}-\ell_{\rm max}\frac{\partial Z_{\rm fix}(\vec{r})}{\partial\ell_{\rm
max}}$ are included on the right-hand side of Eq. (28). However, these terms
can be dropped when $Z_{\rm fix}(\vec{r})$ is identified with
$Z(A)\\!=\\!\exp(-\sigma A)$, which has no internal structure and is
independent of $\ell_{\rm max}$ and $\ell_{\rm min}$. For this reason, we
exclude $U_{V}$ from the $\exp$ of Eq. (27) from the beginning.

##  Appendix B Finsler function and Finsler length

Figure 20:  (a) A local coordinate system $(x^{1},x^{2})$ with the origin at
vertex $i$ on a triangle $\bigtriangleup_{ijk}$, (b) unit Finsler length
$v^{u}_{ij}$ for the diffusion anisotropy of $u$ on the curve $C(x^{1})$ along
the $x^{1}$-axis, (c) unit Finsler length $v^{v}_{ij}$ corresponding to the
variable $v$, and (d) intuitive illustrations of Euclidean length $\ell$ and
Finsler length $s_{u,v}$ between $OP$ along the $x^{1}$\- and $C(x^{1})$-axes,
respectively. The direction-dependent Finsler metrics $g_{u}$ and $g_{v}$
along the $x^{1}$-axis, which is assumed to be parallel to the $x$-axis, from
$i$ to $j$ are obtained by using Finsler functions ${\cal L}_{u}$ and ${\cal
L}_{v}$ such that $g_{u}\\!=\\!1/(\chi_{ij}^{u})^{2}$ and
$g_{v}\\!=\\!1/(\chi_{ij}^{v})^{2}$. Arrows with $x$ and $y$ in (b), (c) and
(d) represent the canonical coordinate axes of ${\bf R}^{2}$.

In this Appendix, we introduce technical details of Finsler metrics for
diffusion anisotropy. The Finsler length can be introduced in a local
coordinate system $(x^{1},x^{2})$ on triangulated lattices in a discrete
manner (Fig. 20(a)). An origin of the local coordinate is the vertex position
$i$ of the triangle $\bigtriangleup_{ijk}$. The $x^{1}$-axis is considered a
curve $C(x^{1})$, on which a Finsler function ${\cal L}$ is introduced using a
coordinate function $E_{u}(\vec{\bf r})\\!=\\!E_{u}(x^{1})$ and its first
derivative $E_{u}^{\prime}(x^{1})\\!=\\!\frac{dE_{u}}{dx^{1}}$ such that

$\displaystyle\begin{split}&{\cal
L}_{u}\left(E_{u},E_{u}^{\prime}\right)=\frac{\sqrt{E_{u}^{\prime\;2}}}{|\chi_{u}|},\\\
&E_{u}=\int_{0}^{x^{1}}\left(\frac{\partial u}{\partial
x^{1}}\right)^{2}dx^{1},\quad
E_{u}^{\prime}=\frac{dE_{u}}{dx^{1}}.\end{split}$ (32)

Here, we use the expression $\left(\frac{\partial u}{\partial
x^{1}}\right)^{2}\\!=\\!(\nabla u)^{2}$ along $C(x^{1})$. $E_{u}(x^{1})$ is an
energy function that is positive and increases with increasing $(\nabla
u)^{2}$ along the $x^{1}$-axis and monotonically increases with respect to
$x^{1}$; hence, $E_{u}(x^{1})$ can be used to define the Finsler length along
$C(x^{1})$. Note that this coordinate function $E_{u}(\vec{\bf
r})\\!=\\!E_{u}(x^{1})$ is different from those in Ref. Krause-etal-
BulMatBiol2019 , where the induced metric is used, and hence, the coordinate
functions are defined on the tangential lines of Euclidean length on curved
surfaces. From the general prescription, we obtain the Finsler length $s$ by
the time integration of the Finsler function ${\cal L}$ along $C(x^{1})$ from
$t\\!=\\!0$ to $t\\!=\\!T$ such that $s\\!=\\!\int_{0}^{T}{\cal
L}(E_{u},E_{u}^{\prime})dt$. Equivalently, we obtain the corresponding Finsler
metric element $g^{u}$ along the $x^{1}$-axis such that
$g^{u}\\!=\\!(1/2)\partial^{2}{\cal L}_{u}^{2}/\partial
E_{u}^{\prime\;2}\\!=\\!1/(\chi_{ij}^{u})^{2}$, which is suitable for our
purpose. Another Finsler metric for anisotropic diffusion of the variable $v$
along $C(x^{1})$ can also be obtained by $g^{v}\\!=\\!(1/2)\partial^{2}{\cal
L}_{v}^{2}/\partial E_{v}^{\prime\;2}\\!=\\!1/(\chi_{ij}^{v})^{2}$. Here, we
assume that $\chi_{ij}^{u}$ and $\chi_{ij}^{v}$ are given by

$\displaystyle\begin{split}&\chi_{ij}^{u}=|\vec{\tau}_{i}\cdot\vec{e}_{ij}|+\chi_{0},\\\
&\chi_{ij}^{v}=\sqrt{1-|\vec{\tau}_{i}\cdot\vec{e}_{ij}|^{2}}+\chi_{0},\end{split}$
(33)

where $\vec{\tau}_{i}$ denotes the internal DOF at vertex $i$, $\vec{e}_{ij}$
is the unit tangential vector from vertices $i$ to $j$, and $\chi_{0}$ is a
small positive number. These parameters $\chi_{ij}^{u}$ and $\chi_{ij}^{v}$ on
$C(x^{1})$ are visualized by the dashed lines in Fig. 20(b) and (c).

Here, we note that the relation between the Euclidean length $dx^{1}$ between
$i$ and $j$ and the Finsler length $ds_{u}$ from $i$ to $j$ is given by using
the Finsler metric $g_{u}$ such that
$(ds_{u})^{2}\\!=\\!g_{u}(dx^{1})^{2}\\!=\\!(dx^{1})^{2}/(\chi_{ij}^{u})^{2}$.
Thus, we have $ds_{u}\\!=\\!dx^{1}/\chi_{ij}^{u}$, and therefore
$ds_{u}\\!=\\!1$ for $dx^{1}\\!=\\!\chi_{ij}^{u}$. This is the reason we call
$\chi_{ij}^{u}$ the unit Finsler length from $i$ to $j$. We sometimes call
$\chi_{ij}^{u}$ a “velocity,” as $\chi_{ij}^{u}$ physically has the units of
velocity because the Finsler length $ds_{u}$ is a time length. Along the
$x^{2}$-axis from vertices $i$ to $k$ in Fig. 20(a), the unit Finsler lengths
$\chi_{ik}^{u}$ and $\chi_{ik}^{v}$ can be defined by replacing $\vec{e}_{ij}$
with $\vec{e}_{ik}$ in Eq. (33). Thus, we obtain the two-dimensional Finsler
metrics $g^{u}$ and $g^{v}$ with respect to the local coordinate
$(x^{1},x^{2})$ such that

$\displaystyle g^{u}=\begin{pmatrix}1/(\chi_{ij}^{u})^{2}&0\\\
0&1/(\chi_{ik}^{u})^{2}\end{pmatrix},\quad
g^{v}=\begin{pmatrix}1/(\chi_{ij}^{v})^{2}&0\\\
0&1/(\chi_{ik}^{v})^{2}\end{pmatrix}.$ (34)

##  Appendix C Discrete Hamiltonian

Figure 21:  (a) Vertices $j(i)$ connected to vertex $i$ and (b) unit
tangential vector $\vec{e}_{ij}$ from vertices $i$ to $j(i)$ and its
decomposition $\vec{e}_{ij}\\!=\\!\vec{e}_{ij}^{\,x}\\!+\\!\vec{e}_{ij}^{\,y}$
into the canonical coordinate axis components.

Using $g_{u}$ and $g_{v}$ in Eq. (34), we have the discrete expression of
$S_{u}$ and $S_{v}$ in Eq. (15). Here, we show the outline of discretization
of $S_{u}$ as follows: The integral is replaced by the sum over triangles with
the determinant $g\\!=\\!\det(g_{u})$ such that
$\int\sqrt{g_{u}}d^{2}x\to\sum_{\bigtriangleup}1/(\chi_{ij}^{u}\chi_{ik}^{u})$,
and the differentials are replaced by differences with the inverse metric
$g_{u}^{ab}\\!=\\!(g_{u})^{-1}$ such that $g_{u}^{ab}\frac{\partial
u}{\partial x^{a}}\frac{\partial u}{\partial
x^{b}}\\!=\\!g_{u}^{11}(\frac{\partial u}{\partial
x^{1}})^{2}\\!+\\!g_{u}^{22}(\frac{\partial u}{\partial x^{2}})^{2}$ $\to$
$(\chi_{ij}^{u})^{2}(u_{j}-u_{i})^{2}\\!+\\!(\chi_{ik}^{u})^{2}(u_{k}-u_{i})^{2}$.
Thus, we have
$S_{u}\\!=\\!\frac{1}{2}\int\sqrt{g_{u}}d^{2}xg_{u}^{ab}\frac{\partial
u}{\partial x^{a}}\frac{\partial u}{\partial x^{b}}$ $\to$
$\frac{1}{2}\sum_{\bigtriangleup}\left(\frac{\chi_{ij}^{u}}{\chi_{ik}^{u}}(u_{j}-u_{i})^{2}\\!+\\!\frac{\chi_{ik}^{u}}{\chi_{ij}^{u}}(u_{k}-u_{i})^{2}\right)$.
Here, we note that there are two different origins of local coordinates at
vertices $j$ and $k$ other than $i$ on triangle $\bigtriangleup_{ijk}$. The
discrete expressions of $S_{u}$ corresponding to these local coordinates on
$\bigtriangleup_{ijk}$ are obtained by replacing the indices $ijk\to jki$ and
$jki\to kij$ and summing over three different expressions with the factor
$1/3$, we obtain

$\displaystyle\begin{split}S_{u}&=\frac{1}{6}\sum_{\bigtriangleup_{ijk}}\left[\frac{\chi_{ij}^{u}}{\chi_{ik}^{u}}(u_{j}-u_{i})^{2}+\frac{\chi_{ik}^{u}}{\chi_{ij}^{u}}(u_{k}-u_{i})^{2}+\frac{\chi_{jk}^{u}}{\chi_{ji}^{u}}(u_{k}-u_{j})^{2}\right.\\\
&\left.+\frac{\chi_{ji}^{u}}{\chi_{jk}^{u}}(u_{i}-u_{j})^{2}+\frac{\chi_{ki}^{u}}{\chi_{kj}^{u}}(u_{i}-u_{k})^{2}+\frac{\chi_{kj}^{u}}{\chi_{ki}^{u}}(u_{j}-u_{k})^{2}\right]\\\
&=\sum_{\bigtriangleup_{ijk}}\left(\bar{\gamma}_{ij}^{u}(u_{i}-u_{j})^{2}+\bar{\gamma}_{jk}^{u}(u_{j}-u_{k})^{2}+\bar{\gamma}_{ki}^{u}(u_{k}-u_{i})^{2}\right),\\\
&\bar{\gamma}_{ij}^{u}=\frac{1}{6}\left(\frac{\chi_{ij}^{u}}{\chi_{ik}^{u}}+\frac{\chi_{ji}^{u}}{\chi_{jk}^{u}}\right),\quad\bar{\gamma}_{jk}^{u}=\frac{1}{6}\left(\frac{\chi_{jk}^{u}}{\chi_{ji}^{u}}+\frac{\chi_{kj}^{u}}{\chi_{ki}^{u}}\right),\quad\bar{\gamma}_{ki}^{u}=\frac{1}{6}\left(\frac{\chi_{ki}^{u}}{\chi_{kj}^{u}}+\frac{\chi_{ik}^{u}}{\chi_{ij}^{u}}\right),\end{split}$
(35)

where the sum over triangles is explicitly denoted by
$\sum_{\bigtriangleup_{ijk}}$. The summation convention can also be changed
from summation over triangles to summation over bonds $\sum_{ij}$. Recalling
that the bond $ij$ is shared by two triangles $\bigtriangleup_{ijk}$ and
$\bigtriangleup_{jil}$ (Fig. 20(a)), we have to include terms
$\frac{1}{6}\left(\frac{\chi_{ij}^{u}}{\chi_{ik}^{u}}+\frac{\chi_{ji}^{u}}{\chi_{jk}^{u}}\right)(u_{i}-u_{j})^{2}$
from $\bigtriangleup_{ijk}$ and
$\frac{1}{6}\left(\frac{\chi_{ji}^{u}}{\chi_{jl}^{u}}+\frac{\chi_{ij}^{u}}{\chi_{il}^{u}}\right)(u_{j}-u_{i})^{2}$
from $\bigtriangleup_{jil}$ in the sum over bonds $\sum_{ij}$. Thus, we have

$\displaystyle\begin{split}S_{u}=\sum_{ij}\gamma_{ij}^{u}(u_{i}-u_{j})^{2},\quad\gamma_{ij}^{u}=\frac{1}{6}\left(\frac{\chi_{ij}^{u}}{\chi_{ik}^{u}}+\frac{\chi_{ji}^{u}}{\chi_{jk}^{u}}+\frac{\chi_{ij}^{u}}{\chi_{il}^{u}}+\frac{\chi_{ji}^{u}}{\chi_{jl}^{u}}\right),\quad\gamma_{ij}^{u}=\gamma_{ji}^{u}.\end{split}$
(36)

Now, we describe the outline of the discretization of the Laplace Beltrami
operator.

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

on triangulated lattices. Since this operator includes second-order
differentials, we adopt an indirect discretization scheme based on the
discrete Hamiltonian $S_{u}$ in Eq. (35). The continuous expression of
$\mathop{}\\!\mathbin{\bigtriangleup}u$ can be obtained by the variational
technique for the continuous $S_{u}$:

$\displaystyle
S_{u}\\!=\\!\frac{1}{2}\int\sqrt{g_{u}}d^{2}xg_{u}^{ab}\frac{\partial
u}{\partial x^{a}}\frac{\partial u}{\partial
x^{b}}\quad\xrightarrow{\rm{variation}}\quad\mathop{}\\!\mathbin{\bigtriangleup}u=-\frac{\partial
S_{u}}{\partial u}=\frac{1}{\sqrt{g}}\frac{\partial}{\partial
x_{a}}\left(\sqrt{g}g^{ab}\frac{\partial u}{\partial x^{b}}\right).$ (38)

Therefore, we obtain the discrete expression
$\mathop{}\\!\mathbin{\bigtriangleup}u_{i}$ from $S_{u}$ in Eq. (35) by the
discrete variational technique

$\displaystyle
S_{u}=\sum_{ij}\gamma_{ij}^{u}(u_{i}-u_{j})^{2}\quad\xrightarrow{\rm{variation}}\quad\mathop{}\\!\mathbin{\bigtriangleup}u_{i}=-\frac{\partial
S_{u}}{\partial u_{i}}=-\frac{\partial}{\partial
u_{i}}\sum_{jk}\gamma_{jk}^{u}(u_{j}-u_{k})^{2}$ (39)

such that

$\displaystyle\mathop{}\\!\mathbin{\bigtriangleup}u_{i}=2\left(\sum_{j(i)}\gamma_{ij}^{u}u_{j}-u_{i}\sum_{j(i)}\gamma_{ij}^{u}\right)=\sum_{j(i)}2\gamma_{ij}^{u}\left(u_{j}-u_{i}\right),$
(40)

where $\sum_{j(i)}$ denotes the sum over vertices $j$ connected to vertex $i$
with bond $ij$ (Figs. 21(a) and 6(b)). To perform this summation, we replace
the sum $\sum_{jk}$ over bonds $jk$ in Eq. (39) with
$\frac{1}{2}\sum_{j}\sum_{k(j)}(=\\!\sum_{jk})$, where $\sum_{j}$ and
$\sum_{k(j)}$ are the sum over vertices $j$ and the sum over vertices $k$,
which is connected to $j$ with bond $jk$, respectively. The factor $1/2$
appears due to the duplicated sums in $\sum_{j}\sum_{k(j)}$; the term with
index $jk\\!=\\!12(\Leftrightarrow j\\!=\\!1,k\\!=\\!2)$ appears twice in
$\sum_{j}\sum_{k(j)}$, for example. Note that
$\sum_{j}\sum_{k(j)}\\!=\\!\sum_{k}\sum_{j(k)}$ because
$\sum_{jk}\\!=\\!\sum_{kj}$. Thus, we have

$\displaystyle\begin{split}&\frac{\partial}{\partial
u_{i}}\sum_{jk}\gamma_{jk}^{u}(u_{j}-u_{k})^{2}=\frac{\partial}{\partial
u_{i}}\sum_{jk}\gamma_{jk}^{u}(u_{j}^{2}+u_{k}^{2}-2u_{j}u_{k})\\\
&=\frac{1}{2}\sum_{j}\sum_{k(j)}\gamma_{jk}^{u}2(u_{j}\delta_{ji}+u_{k}\delta_{ki}-\delta_{ij}u_{k}-u_{j}\delta_{ik})\\\
&=\sum_{j}\sum_{k(j)}\gamma_{jk}^{u}u_{j}\delta_{ji}+\sum_{j}\sum_{k(j)}\gamma_{jk}^{u}u_{k}\delta_{ki}-\sum_{j}\sum_{k(j)}\gamma_{jk}^{u}\delta_{ij}u_{k}-\sum_{j}\sum_{k(j)}\gamma_{jk}^{u}u_{j}\delta_{ik},\end{split}$
(41)

where $\delta_{ij}$ is the Kronecker delta. In this final expression, the
first term is
$\sum_{j}\sum_{k(j)}\gamma_{jk}^{u}u_{j}\delta_{ji}\\!=\\!\sum_{k(i)}\gamma_{ik}^{u}u_{i}\\!=\\!u_{i}\sum_{k(i)}\gamma_{ik}^{u}$,
and the second term can be written as
$\sum_{j}\sum_{k(j)}\gamma_{jk}^{u}u_{k}\delta_{ki}\\!=\\!\sum_{k}\sum_{j(k)}\gamma_{jk}^{u}u_{k}\delta_{ki}\\!=\\!\sum_{j(i)}\gamma_{ji}^{u}u_{i}\\!=\\!u_{i}\sum_{j(i)}\gamma_{ji}^{u}$,
which is identical to the first term. The third term is
$-\\!\sum_{j}\sum_{k(j)}\gamma_{jk}^{u}\delta_{ij}u_{k}\\!=\\!-\\!\sum_{k(i)}\gamma_{ik}^{u}u_{k}$,
which is identical to the fourth term
$-\\!\sum_{j}\sum_{k(j)}\gamma_{jk}^{u}u_{j}\delta_{ik}\\!=\\!-\\!\sum_{k}\sum_{j(k)}\gamma_{jk}^{u}u_{j}\delta_{ik}\\!=\\!-\\!\sum_{j(i)}\gamma_{ji}^{u}u_{j}\\!=\\!-\\!\sum_{j(i)}\gamma_{ij}^{u}u_{j}$.
This proves that Eq. (40).

Here, we comment on the relation between
$\mathop{}\\!\mathbin{\bigtriangleup}$ in Eq. (40) and the network Laplacian
$\sum_{j=1}^{N}L_{ij}u_{j}\\!=\\!\sum_{j=1}^{N}A_{ij}(u_{j}\\!-\\!u_{i})$,
where $N$ is the total number of nodes (or vertices) and
$L_{ij}\\!=\\!A_{ij}-k_{i}\delta_{ij}$ with
$k_{i}\\!=\\!\sum_{j=1}^{N}A_{ij}$. The adjacency (or connectivity) matrix
$A_{ij}$ is defined by $A_{ij}\\!=\\!1$ if $i$ and $j$ are connected and
$A_{ij}\\!=\\!0$ otherwise. Rewriting
$\mathop{}\\!\mathbin{\bigtriangleup}u_{i}$ in Eq.(40) as

$\displaystyle\begin{split}\mathop{}\\!\mathbin{\bigtriangleup}u_{i}=\sum_{ij}2\gamma_{ij}^{u}\left(A_{ij}-k_{i}\delta_{ij}\right),\end{split}$
(42)

we find that $\mathop{}\\!\mathbin{\bigtriangleup}u_{i}$ is considered to be a
weighted network Laplacian with the weight $2\gamma_{ij}^{u}$. Note that
$\gamma_{ij}^{u}$ is well defined for networks on two-dimensional surfaces, as
shown in Fig. 8. In the isotropic case, $\chi_{0}\to\infty$ in Eq. (33), we
have $\chi_{ij}^{u}/\chi_{kl}^{u}\\!\to\\!1$ and $\gamma_{ij}\\!\to\\!2/3$
from Eq. (36), and therefore, $2\gamma_{ij}^{u}\\!\to\\!4/3$. Thus, with a
suitable normalization factor in the definition of $\gamma_{ij}^{u}$ in Eq.
(36), $\mathop{}\\!\mathbin{\bigtriangleup}u_{i}$ in Eq. (42) reduces to the
standard network Laplacian
$\mathop{}\\!\mathbin{\bigtriangleup}u_{i}\to\sum_{ij}\left(A_{ij}-k_{i}\delta_{ij}\right)$.

##  Appendix D Direction-dependent diffusion constants and effective density
energies

In this Appendix, we describe direction-dependent diffusion constants and
corresponding effective direction-dependent energies obtained by the discrete
Laplace operator in Eq. (19). We simply call these energies “density energies”
because the variables $u$ and $v$ represent the density of the activator and
inhibitor, respectively. First, we note that the unit tangential vector
$\vec{e}_{ij}$ from vertices $i$ to $j$ can be decomposed into
$\vec{e}_{ij}\\!=\\!\vec{e}_{ij}^{\,x}\\!+\\!\vec{e}_{ij}^{\,y}$ (Fig. 21(b)).
Then, by using the angle $\theta_{ij}$ and the relations
$1\\!=\\!\vec{e}_{ij}\cdot\vec{e}_{ij}\\!=\\!(e_{ij}^{x})^{2}\\!+\\!(e_{ij}^{y})^{2}\\!=\\!\cos^{2}\theta_{ij}\\!+\\!\sin^{2}\theta_{ij}$
satisfying $\cos^{2}\theta_{ij}\\!=\\!\sin^{2}\theta_{ij}\\!=\\!1/2$ for
$\theta_{ij}\\!=\\!\pi/4$, it is easy to see that $\gamma_{ij}^{u}$ is
decomposed into two different parts, which can be considered the $x$ and $y$
components such that

$\displaystyle\begin{split}\gamma_{ij}^{u}=\gamma_{ij}^{u,x}+\gamma_{ij}^{u,y}=\gamma_{ij}^{u}\cos^{2}\theta_{ij}+\gamma_{ij}^{u}\sin^{2}\theta_{ij}.\end{split}$
(43)

This $\gamma_{ij}^{u}$ is considered to have values on bond $ij$. Here, we
assume an isotropic case in which $\gamma_{ij}^{u}$ is very large on bonds
that are almost parallel to the $x$-axis
($\Leftrightarrow\theta_{ij}\\!\simeq\\!0$ in Fig. 21(b)), and additionally,
$\gamma_{ij}^{u}$ is assumed to be very small on bonds that are almost
parallel to the $y$-axis ($\Leftrightarrow\theta_{ij}\\!\simeq\\!\pi/2$ in
Fig. 21(b)). In this case, it is easy to understand that
$\gamma_{ij}^{u,x}\\!>\\!\gamma_{ij}^{u,y}$ on bonds $ij$ almost parallel to
the $x$-axis because $\cos^{2}\theta_{ij}\\!\simeq\\!1$ and
$\sin^{2}\theta_{ij}\\!\simeq\\!0$, and that
$\gamma_{ij}^{u,x}\\!<\\!\gamma_{ij}^{u,y}$ on bonds $ij$ almost parallel to
the $y$-axis because $\cos^{2}\theta_{ij}\\!\simeq\\!0$ and
$\sin^{2}\theta_{ij}\\!\simeq\\!1$. Therefore, such a direction dependence of
$\gamma_{ij}^{u}$ is considered to be reflected in $\gamma_{ij}^{u,x}$ and
$\gamma_{ij}^{u,y}$. Moreover, the mean values of $\gamma_{ij}^{u,x}$ and
$\gamma_{ij}^{u,y}$, denoted by $\langle\gamma_{ij}^{u,x}\rangle$ and
$\langle\gamma_{ij}^{u,y}\rangle$, become isotropic in the case that
$\gamma_{ij}^{u}$ is constant or randomly distributed over the lattice.
Indeed, in the case of constant $\gamma_{ij}^{u}$, the lattice averages
$\langle\gamma_{ij}^{u,x}\rangle$ and $\langle\gamma_{ij}^{u,y}\rangle$ depend
only on the mean values $\langle\cos^{2}\theta_{ij}\rangle$ and
$\langle\sin^{2}\theta_{ij}\rangle$, which are the same on the regular square
lattice and on the regular hexagonal lattice. On randomly triangulated
lattices,
$\langle\cos^{2}\theta_{ij}\rangle\\!\simeq\\!\langle\sin^{2}\theta_{ij}\rangle$.
In the case of randomly distributed $\gamma_{ij}^{u}$, it is also natural to
expect that
$\langle\gamma_{ij}^{u}\cos^{2}\theta_{ij}\rangle\\!\simeq\\!\langle\gamma_{ij}^{u}\sin^{2}\theta_{ij}\rangle$
on both regular-triangular and random lattices. Thus, we have direction-
dependent diffusion constants defined by

$\displaystyle
D_{u}^{x}=\frac{1}{N_{B}}\sum_{ij}\gamma_{ij}^{u}\cos^{2}\theta_{ij}=\langle\gamma_{ij}^{u}\cos^{2}\theta_{ij}\rangle,\quad
D_{u}^{y}=\frac{1}{N_{B}}\sum_{ij}\gamma_{ij}^{u}\sin^{2}\theta_{ij}=\langle\gamma_{ij}^{u}\sin^{2}\theta_{ij}\rangle,$
(44)

where $N_{B}\\!=\\!\sum_{ij}1(=\\!3N)$ denotes the total number of bonds.

It is also possible to define effective density energies, corresponding to
$D_{u}^{x}$ and $D_{u}^{y}$ in Eq. (44), such that

$\displaystyle\begin{split}&S_{u}^{x}=\frac{1}{\langle\gamma_{ij}^{u}\cos^{2}\theta_{ij}\rangle}\sum_{ij}\gamma_{ij}^{u}\cos^{2}\theta_{ij}(u_{i}-u_{j})^{2},\\\
&S_{u}^{y}=\frac{1}{\langle\gamma_{ij}^{u}\sin^{2}\theta_{ij}\rangle}\sum_{ij}\gamma_{ij}^{u}\sin^{2}\theta_{ij}(u_{i}-u_{j})^{2}.\end{split}$
(45)

Indeed, $S_{u}$ in Eq. (36) can be expressed as

$\displaystyle\begin{split}S_{u}&=\sum_{ij}\gamma_{ij}^{u}(u_{i}-u_{j})^{2}=\sum_{ij}\left(\gamma_{ij}^{u}\cos^{2}\theta_{ij}+\gamma_{ij}^{u}\sin^{2}\theta_{ij}\right)(u_{i}-u_{j})^{2}\\\
&=\langle\gamma_{ij}^{u}\cos^{2}\theta_{ij}\rangle\frac{\sum_{ij}\gamma_{ij}^{u}\cos^{2}\theta_{ij}(u_{i}-u_{j})^{2}}{\langle\gamma_{ij}^{u}\cos^{2}\theta_{ij}\rangle}+\langle\gamma_{ij}^{u}\sin^{2}\theta_{ij}|\rangle\frac{\sum_{ij}\gamma_{ij}^{u}\sin^{2}\theta_{ij}(u_{i}-u_{j})^{2}}{\langle\gamma_{ij}^{u}\sin^{2}\theta_{ij}|\rangle}\\\
&=D_{u}^{x}S_{u}^{x}+D_{u}^{y}S_{u}^{y}.\end{split}$ (46)

##  Appendix E Bond-length distribution

Figure 22: Bond-length distribution $h(\ell)$ vs. $\ell$ of the fixed model
(red line) and fluid model (blue line) with (a) $\lambda\\!=\\!0$,
$\vec{F}\\!=\\!(0,0)$ and (b) $\lambda\\!=\\!0.5$, $\vec{F}\\!=\\!(0,2)$,
which respectively correspond to the isotropic and anisotropic patterns
plotted in Figs. 9(a), 10(a) and Figs. 9(d), 10(d). The diffusion coefficients
are $(D_{u},D_{v})\\!=\\!(0.2,5)$ for the fixed model and
$(D_{u},D_{v})\\!=\\!(0.2,10)$ for the fluid model, and the other parameters
$(\alpha,\gamma)\\!=\\!(1,8)$ are the same for the two different models. The
cutoff is $\ell_{\rm max}\\!=\\!3d\\!=\\!1.575$, where $d(=\\!0.525)$ is the
assumed lattice spacing. The curves $h(\ell)$ vs. $\ell$ for $d\\!=\\!0.41$
are almost the same as those in the figures.

In this Appendix, we plot the bond length distribution $h(\ell)$ on the fixed
and fluid lattices of size $N\\!=\\!6400$ in Fig. 22(a), (b). The lattice size
is the same as that assumed for the calculations of the diffusion constants
$D_{u}^{\mu}$, $D_{v}^{\mu}$ and $\sigma$ as plotted in Figs. 11, 12 and 13,
respectively. The total numbers of iterations $n_{\rm itr}$ and $n_{\rm MC}$
are also the same as those in Eq. (24). The coefficients are
$(D_{u},D_{v})\\!=\\!(0.2,5)$ on the fixed model, while
$(D_{u},D_{v})\\!=\\!(0.2,10)$ on the fluid model. These parameters, including
$(\alpha,\gamma)$, are the same as those assumed in the calculation of
isotropic and anisotropic snapshots in Fig. 9(a) and (d) for the fixed model
and in Fig. 10(a) and (d) for the fluid model. We find from both isotropic and
anisotropic cases in (a) and (b) that the distribution of $h(\ell)$ on the
fluid lattice is slightly wide in both directions of small and large $\ell$.
This slightly broad spectrum of $h(\ell)$ in the fluid model is due to the
free diffusion of vertices, implying that the external force $\vec{F}$
influences the movement of vertices more strongly on the fluid lattice than on
the fixed lattice. In other words, the vertex movement, which involves free
diffusion, on fluid lattices is caused by a relatively weak external force
$\vec{F}$.

## References

  * (1) A. M. Turing, Chemical Basis of Morphogenesis, Phil. Trans. Roy. Soc. London 237, pp.37-72 (1952).
  * (2) A. Koch and H. Meinhardt, Rev. Mod. Phys. 66, No. 4 pp.1491-1507 (1994).
  * (3) D. Fitzhugh, Impulses and Physiological States in Theoretical Models of Nerve membrane, Biophysical J. 1, pp.446-466 (1961).
  * (4) J. Nagumo, S. Arimioto, and S. Yoshizawa, An Active Pulse Transmission Line Simulating Nerve Axon, Proc. of the IRE. pp.2061-2070 (1962).
  * (5) T. Sekimura, A. Madzvamuse, A.J. Wathen and P.K. Maini, A model for colour pattern formation in the butterfly wing of Papilio dardanus, Proc. R. Soc. Lond. B 267, pp. 851-859 (2000), DOI: 10.1098/rspb.2000.1081.
  * (6) S. Kondo and T. Miura, Reaction-Diffusion Model as a Framework for Understanding Biological Pattern Formation Science 329, pp.1616-1620 (2010), DOI: 10.1126/science.1179047.
  * (7) D. Bullara and Y. De Decker, Pigment cell movement is not required for generation of Turing patterns in zebrafish skin Nature Com. 6, 6971 (2015), DOI: 10.1038/ncomms7971.
  * (8) Z. Tan, S. Chen, X. Peng, L. Zhang, C. Gao, Polyamide membranes with nanoscale Turing structures for water purification, Science 360, pp. 518?521 (2018), DOI: 10.1126/science.aar6308.
  * (9) Y. Fuseya, H. Katsuno, K. Behnia and A. Kapitulnik, Nanoscale Turing patterns in a bismuth monolayer, Nature Phys. 17, pp.1031-1036 (2021), https://doi.org/10.1038/s41567-021-01288-y.
  * (10) H. Nakao, and A.S. Mikhailov, Turing patterns in network-organized activator?inhibitor systems, Nature Physics 6, pp.544-550 (2010), DOI: 10.1038/NPHYS1651.
  * (11) T. Carletti and H. Nakao, Turing patterns in a network-reduced FitzHugh-Nagumo model, Phys, Rev. E 101, 022203 (2020), DOI: 10.1103/PhysRevE.101.022203.
  * (12) M. Asllani, J.D. Challenger, F.S. Pavone, L. Sacconi and D. Fanelli, The theory of pattern formation on directed networks, Nature Comm. 5 4517 (2014) DOI: 10.1038/ncomms5517.
  * (13) M. Asllani, D.M. Busiello, T. Carletti, D. Fanelli, and G. Planchon, Turing patterns in multiplex networks, Phys. Rev. E 90 042814 (2014), DOI: 10.1103/PhysRevE.90.042814.
  * (14) J. Petit, M. Asllani, D. Fanelli, B. Lauwens, T. Carletti, Pattern formation in a two-component reaction?diffusion system with delayed processes on a network, Physica A 462 pp. 230-249 (2016), http://dx.doi.org/10.1016/j.physa.2016.06.003.
  * (15) A. Alexandre, M. Lavaud, N. Fares, E. Millan, Y. Louyer, T. Salez, Y. Amarouchene, T. Gu?rin, and D.S. Dean, Non-Gaussian Diffusion Near Surfaces, Phys. Rev. Lett 130 077101 (2023), DOI:https://doi.org/10.1103/PhysRevLett.130.077101.
  * (16) P.C. Bressloff and J.M. Newby, Stochastic models of intracellular transport, Rev. Mod. Phys. 85 pp. 135-196 (2013), DOI: 10.1103/RevModPhys.85.135.
  * (17) F. H${\ddot{\rm o}}$fling and T Franosch, Anomalous transport in the crowded world of biological cells, Rep. Prog. Phys. 76 046602(50pp) (2013), doi:10.1088/0034-4885/76/4/046602.
  * (18) I.M. Sokolov, J. Klafter, and A. Blumen, Fractional Kinetics, Physics Today 55, 11, 48 (2002), https://doi.org/10.1063/1.1535007.
  * (19) R. Metzler and J. Klafter, The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics, J. Phys. A: Math. Gen. 37, R161?R208 (2004), DOI 10.1088/0305-4470/37/31/R01.
  * (20) G. Ariel, A. Shklarsh, O. Kalisman, C. Ingham and E. Ben-Jacob, From organized internal traffic to collective navigation of bacterial swarms, New J. Phys. 15, 125019 (2013), doi:10.1088/1367-2630/15/12/125019.
  * (21) G. Ariel, A. Rabani, S. Benisty, J. D. Partridge, R. M. Harshey and A.Be’er, Swarming bacteria migrate by L$\acute{\rm e}$vy Walks, Nature Comm. 6: 8936 (2015), DOI: 10.1038/ncomms9396.
  * (22) P.H. Dederichs and K. Schroeder, Anisotroyic diffusion in stress fields, Phys. Rev. B 17 pp. 2524-2536 (1978).
  * (23) Xu-D. Xiao, X.B. Zhu, W. Baum, and Y.R. Shen, Anisotropic Surface Diffusion of CO on Ni(110), Phys. Rev. Lett. 66 pp. 2352-2355 (1991).
  * (24) B.P.J. Bret and A. Lagendijk, Anisotropic enhanced backscattering induced by anisotropic diffusion, Phys. Rev. E. 70 036601 (2004).
  * (25) N. Kerber, M. Wei${\rm\beta}$enhofer, Klaus Raab, K. Litzius, J. Z$\acute{\rm a}$zvorka, U. Nowak, and M. Kl$\ddot{\rm a}$ui, Anisotropic Skyrmion Diffusion Controlled by Magnetic-Field-Induced Symmetry Breaking, Phys. Rev. Appl. 15 044029 (2021).
  * (26) S. Kondo and R. Asai, A reaction-diffusion wave on the skin of the marine anglefish Pomacanthus , Nature 376, pp.765-768 (1995).
  * (27) H. Shoji, A. Mochizuki, Y. Iwasa, M. Hirata, T. Watanabe, S. Hioki, and S. Kondo, Origin of Directionality in the Fish Stripe Pattern, Developemental Dynamics 226, pp.627-633 (2003).
  * (28) R. Iwamoto and H. Shoji, Kakusan Ihousei Turing Pattern, RIMS Kokyuroku (in Japanese) 2087, pp.108-117 (2018).
  * (29) C. Varea, J.L. Arag$\check{\rm o}$n, and R.A. Barrio, Turing patterns on a sphere, Phys. Rev. E 60, pp.4588-4592 (1999), DOI: 10.1103/physreve.60.4588.
  * (30) A.L. Krause, M.A. Ellis, R.A. Van Gorder, Influence of Curvature, Growth, and Anisotropy on the Evolution of Turing Patterns on Growing Manifolds, Bulletin of Mathematical Biology 81, pp.759-799 (2019), https://doi.org/10.1007/s11538-018-0535-y.
  * (31) H. Koibuchi, M. Okumura and S. Noro, J. Phys. Conf. Ser. 1730, 012035 (2020), doi:10.1088/1742-6596/1730/1/012035.
  * (32) H. Koibuchi, F. Kato, G. Diguet, B. Ducharne and T. Uchimoto, Origin of anisotropic diffusion in Turing Patterns, http://arxiv.org/abs/2208.09977, To appear in AIP Conf. Ser.
  * (33) H. Koibuchi and H. Sekino, Monte Carlo studies of a Finsler geometric surface model, Physica A 303, pp.37-50 (2014), https://doi.org/10.1016/j.physa.2013.08.006.
  * (34) Y. Takano and H. Koibuchi, J-shaped stress-strain diagramof collagen fibers: Frame tension of triangulated surfaces with fixed boundaries, Phys. Rev. E 95, 042411(1-11) (2017), https://doi.org/10.1103/PhysRevE.95.042411.
  * (35) E. Proutorov, N. Matsuyama and H. Koibuchi, Finsler geometry modeling and Monte Carlo study of liquid crystal elastomers under electric fields, J. Phys. Cond. Mat. 30, 405101 (2018), https://doi.org/10.1088/1361-648X/aadcba.
  * (36) S.El Hog, F. Kato, H. Koibuchi, and H.T. Diep, Finsler geometry modeling and Monte Carlo study of skyrmion shape deformation by uniaxial stress, Phys. Rev. B 105, 024402 (2021), https://doi.org/10.1103/physrevb.104.024402.
  * (37) S.El Hog, F. Kato, S. Hongo, H. Koibuchi, G. Diguet, T. Uchimoto and H.T. Diep, The stability of 3D skyrmions under mechanical stress studied via Monte Carlo calculations, Results in Phys. 38, 105578 (2022), https://doi.org/10.1016/j.rinp.2022.105578.
  * (38) R. Friedberg and H.-C. Ren, Field theory on a computationally constructed random lattice, Nucl. Phys. B 235, pp. 310-320 (1984), https://doi.org/10.1016/0550-3213(84)90501-7.
  * (39) N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, and A.H. Teller, Equation of State Calculations by Fast Computing Machines, J. Chem. Phys. 21, 1087 (1953), https://doi.org/10.1063/1.1699114.
  * (40) D.P. Landau, Finite-size behavior of the simple-cubic Ising lattice, Phys. Rev. B 13, 2997 (1976),] DOI:https://doi.org/10.1103/PhysRevB.14.255.
  * (41) Y. Kantor and D.R. Nelson, Phase transitions in flexible polymeric surfaces, Phys. Rev. A 36, 4020 (1987), DOI:https://doi.org/10.1103/PhysRevA.36.4020.
  * (42) G. Gompper and D. M. Kroll, Shape of inflated vesicles, Phys. Rev. A 46, 7466 (1992), DOI:https://doi.org/10.1103/PhysRevA.46.7466.
  * (43) H. Koibuchi and T. Kuwahata, First-order phase transition in the tethered surface model on a sphere, Phys. Rev. E 72, 026124 (2005), DOI:https://doi.org/10.1103/PhysRevE.72.026124.
  * (44) J. -S. Ho and A. Baumg${\ddot{\rm a}}$rtner, Simulations of Fluid Self-Avoiding Membranes, Europhys. Lett. 12, 295 (1990), DOI 10.1209/0295-5075/12/4/002.
  * (45) W. Z. Helfrich, Elastic Properties of Lipid Bilayers: Theory and Possible Experiment, Naturforsch 28c, 693 (1973), DOI: 10.1515/znc-1973-11-1209.
  * (46) L. Peliti and S. Leibler, Effects of Thermal Fluctuations on Systems with Small Surface Tension, Phys. Rev. Lett. 54, 1690 (1985), DOI:https://doi.org/10.1103/PhysRevLett.54.1690.
  * (47) D. Nelson, The Statistical Mechanics of Membranes and Interfaces, in Statistical Mechanics of Membranes and Surfaces, Second Edition, eds. D. Nelson, T. Piran, and S. Weinberg (World Scientific, Singapore, 2004), p. 1.
  * (48) J F Wheater, Random surfaces: from polymer membranes to strings, J. Phys. A: Math. Gen. 27, 3323 (1994), DOI 10.1088/0305-4470/27/10/009).
  * (49) M. Doi and F. Edwards, The Theory of Polymer Dynamics (Oxford University Press, 1986).