PHYSICAL REVIEW E 99, 032416 (2019)
Editors’ Suggestion

Featured in Physics

Pattern formation of skin cancers: Effects of cancer proliferation and hydrodynamic interactions
Takuma Hoshino
Department of Chemistry, Graduate School of Science, Tokyo Metropolitan University, Tokyo 192-0397, Japan

Ming-Wei Liu and Kuo-An Wu
Department of Physics, National Tsing Hua University, Hsinchu 30013, Taiwan

Hsuan-Yi Chen
Department of Physics, National Central University, Jhongli 32001, Taiwan
and Institute of Physics, Academia Sinica, Taipei 11529, Taiwan

Tatsuaki Tsuruyama
Center for Anatomical Studies, Kyoto University Graduate School of Medicine, Kyoto 606-8501, Japan

Shigeyuki Komura*
Department of Chemistry, Graduate School of Science, Tokyo Metropolitan University, Tokyo 192-0397, Japan(Received 3 January 2019; published 21 March 2019)
We study pattern formation of skin cancers by means of numerical simulation of a binary system consistingof cancer and healthy cells. We extend the conventional model H for macrophase separations by considering
a logistic growth of cancer cells and also a mechanical friction between dermis and epidermis. Importantly,
our model exhibits a microphase separation due to the proliferation of cancer cells. By numerically solving
the time evolution equations of the cancer composition and its velocity, we show that the phase separation
kinetics strongly depends on the cell proliferation rate as well as on the strength of hydrodynamic interactions.
A steady-state diagram of cancer patterns is established in terms of these two dynamical parameters and someof the patterns correspond to clinically observed cancer patterns. Furthermore, we examine in detail the time
evolution of the average composition of cancer cells and the characteristic length of the microstructures. Our
results demonstrate that different sequence of cancer patterns can be obtained by changing the proliferation rate
and/or hydrodynamic interactions.
DOI: 10.1103/PhysRevE.99.032416
I. INTRODUCTION

Tissue morphogenesis is a process in which multicellularorganisms are dynamically formed in a coherent manner [1].
Several deterministic and stochastic models to describe tissue
regeneration using such as stem cells have been proposed froma theoretical point of view [2,3]. Recently, various analogies
between viscoelastic fluids and biological tissues have been
pointed out to investigate mechanical response of a biologicaltissue to an applied force [4–7]. Needless to say, studies
on tumor dynamics are directly connected with medical diagnosis and there have been many attempts to understandcancer behaviors across multiple biological scales [8–11].
Although some correlations between cancer patterns and theirmalignancies are realized, it is not well-understood why andhow such malignant patterns appear in tissues. For example,
a skin cancer called melanoma often exhibits characteristic
surface patterns which are diagnosed by medical doctors[12]. However, fundamental mechanisms that underlie sucha pattern formation need to be further investigated.

*komura@tmu.ac.jp

2470-0045/2019/99(3)/032416(13)Recently, some dynamical studies on skin lesions havebeen performed to discuss the morphological changes in earlymelanoma development by using a phase separation model[13–16]. Among these works, Chatelain et al. investigated abinary system composed of cancer and healthy cells. They
demonstrated that not only the cell-cell adhesion but also
the coupling to the diffusion of nutrients (oxygen) leads
to the microstructure (e.g., “dots” and “nests”) formation in
the early stage melanoma [13–15]. These microstructures areanalogous to those in block copolymer systems [17]. In themodel by Chatelain et al., the domain coarsening takes place
due to diffusion process whereas hydrodynamic interactions
are not considered. Hence, their model can be regarded as
an extension of “model B” [18–20] to take into account theformation of microstructures. For bacterial colonies without hydrodynamic interactions, an arrested phase separation
was explained only by considering a local density-dependent
motility and the birth/death of bacteria [21].
In general, a biological tissue can be regarded as a viscoelastic material because it responds like a solid with finite
elasticity at short time scales and behaves like a fluid with
an effective viscosity at long time scales [4–7]. Since the“differential adhesion hypothesis” was proposed by Steinberg

032416-1©2019 American Physical Society

HOSHINO, LIU, WU, CHEN, TSURUYAMA, AND KOMURAPHYSICAL REVIEW E 99, 032416 (2019)[22,23], the similarities between tissues and liquids have beenrecognized for a long time. For example, by using particle tracking velocimetry in gastrulating Drosophila embryos,
it was shown that cytoplasmic redistribution during ventral
furrow formation is described by the presence of hydrodynamic flows [24]. In a recent study of tissue dynamics ofa stratified epithelium, it was shown that a steady hydrodynamic flow of stratified epithelium is controlled by the cellproliferation rate [25,26]. Although these works highlight theimportance of liquid flows in the tissue dynamics, the effectsof hydrodynamic interactions on the skin cancer dynamics
have not been considered so far.
For ordinary fluid mixtures, however, it is well-known thathydrodynamic interactions play crucial roles in their phaseseparation dynamics. This is because the convection of thecomposition field kinetically enhances the phase separationin the presence of flows. The standard model that takes into
account the hydrodynamic effects is called “model H” that has
been extensively studied in the literature [18–20]. For ordinary
3D fluid mixtures, model H predicts that the domain size
increases linearly with time [27,28]. This is much faster thanthe Brownian coagulation process [29] or the Lifshitz-Slyozov
evaporation-condensation process [30].
In this paper, we study the pattern formation of skincancers by means of numerical simulation of a binary systemcomposed of cancer and healthy cells. Our main focus is to investigate the effects of cancer proliferation and hydrodynamic
interactions on the phase separation kinetics. For this purpose,
we shall extend the conventional model H by incorporating
a logistic growth of cancer cells and a mechanical friction
between dermis and epidermis. Similar to chemically reactive binary fluid mixtures [31,32] or block copolymer melts[33–35], our model also exhibits a microphase separation due
to the proliferation of cancer cells.
Performing numerical simulations of the time evolution of
the cancer cell composition and the velocity field, we showthat the phase separation dynamics is strongly affected by
the cell proliferation rate as well as by the strength of hydrodynamic interactions. We shall examine in detail how theaverage composition of cancer cells and the characteristic sizeof microstructures depend on these dynamical parameters.
Our results also demonstrate that different sequence of cancer
patterns can be obtained by changing the cancer proliferation
rate and/or the hydrodynamic effects. Furthermore, our model
can reproduce some of the clinically observed microstructuresin melanoma.
In the next section, we discuss the dynamical equationsof a binary cell system in the presence of hydrodynamic
interactions. In Sec. III, we present our simulation resultsfor different proliferation rates and friction coefficients, andsummarize them in terms of a steady-state diagram as afunction of these parameters. For qualitative arguments, we
further perform structure analysis of the obtained patterns
and give a scaling argument for the observed microphaseseparation. In Sec. IV, we discuss the mechanisms for patternformation in the early and late stages by using the amplitude
equations method and the sharp interface model, respectively.
Finally, the summary of our work and some discussions are
given in Sec. V.

FIG. 1. Schematic illustration of an epidermal tissue on dermis.
The cell layer is assumed to be thin enough so that it can be
regarded as a 2D fluid with hydrodynamic flows. The fluid sheet isinfinitely large and we do not consider any out-of-plane deformationof the epidermal layer. The cell layer is composed of cancer cells(shown in black) and healthy cells (shown in white), and their areal
compositions are defined by φ and ψ, respectively. The two types ofcell fill all the available space and satisfy the saturation constraint,i.e., φ + ψ = 1. Further, the local velocities are denoted by vφ andvψ for cancer and healthy cells, respectively. We also take into
account a mechanical friction between dermis and epidermis that ischaracterized by the friction coefficient ζ .
II. MODEL
A. Continuity equations

Let us consider an epidermal cell layer on dermis as
schematically depicted in Fig. 1. The cell layer is assumed
to be thin enough such that it can be regarded as a twodimensional (2D) system characterized by a 2D vector r =(x, y). Here we do not consider any out-of-plane deformation of the epidermal layer. We assume that the cell layer
is composed of cancer cells and healthy cells whose area
fractions are denoted by φ(r, t ) and ψ (r, t ) (0  φ  1 and0  ψ  1), respectively, which depend on time t. For thehydrodynamic description, we define the corresponding local
velocities by vφ (r, t ) and vψ (r, t ) for cancer cells and healthy
cells, respectively. We further assume that the two types ofcell fill all the available space and always satisfy the saturation constraint φ + ψ = 1 at every point. This saturation
constraint leads to the following incompressibility condition:∇ · v = 0,(1)
where we have introduced the local average velocity,v = φvφ + ψvψ ,(2)which is weighted by the respective area fractions.
To take into account the proliferation of cancer cells andthe death of healthy cells simultaneously, we consider thefollowing continuity equations that are consistent with theabove incompressibility condition:

032416-2∂φ+ ∇ · (φvφ ) = (φ),∂t∂ψ+ ∇ · (ψvψ ) = −(φ),
∂t

(3)
(4)

PATTERN FORMATION OF SKIN CANCERS: EFFECTS OF …

PHYSICAL REVIEW E 99, 032416 (2019)

where the function (φ) represents the compositiondependent cancer proliferation rate of epidermal cells. Among
various possibilities, we choose here the following logistic
growth function:


φ
,
(5)
(φ) = γ φ 1 −φ∞coefficient L is independent of φ [41], because a composition
dependent transport coefficient would not alter the asymptotic
dynamics [42,43]. In the present work, we do not include any
stochastic noise.
The last term −ζ v in Eq. (7) represents the frictional
dissipation between the epidermal layer and dermis. In humantissues, such a friction arises from the adhesion of integrinsthat connect a keratin intracellular network to collagen fibersof basement membranes. With this frictional term, the total
momentum is no longer conserved within the 2D fluid sheet.
Furthermore, the friction coefficient ζ controls the strengthof hydrodynamic interactions. Namely, hydrodynamics does
not play any role when ζ → ∞, whereas hydrodynamic
interactions are fully present when ζ → 0. Later, we shallsystematically change the value of ζ to investigate the effectsof hydrodynamic interactions on the phase separation kinetics.
To further obtain the chemical potential μ and the stresstensor , we introduce the total free energy describing thephase separation of a cell mixture. Following Wise et al. who
discussed a continuum model of multi-species tumor growth[44], we use the following form for a binary cellular system:1
F = dr 2 [φ ln φ + (1 − φ) ln(1 − φ)
a β

κ
2
(8)
+ χ φ(1 − φ)] + (∇φ) .
2

where the coefficient γ > 0 is the cancer proliferation rate in
the epidermal layer. Such a logistic growth was considered
before to describe the effects of birth and death in bacterialcolonies [21]. Starting from an initial average composition,
φ0 , the cancer cell composition evolves toward a higher composition, φ∞ , whose value is roughly given by φ∞ ≈ 0.6-0.8
depending on the cancer cell type [36]. Since the function
(φ) is positive, cancer cells proliferate during the phaseseparation while healthy cells die out due to the invasionof increased cancer cells, as described by Eq. (4). Since thetime-evolution of healthy cells is simply given by ψ (r, t ) =1 − φ(r, t ) due to the saturation condition, we shall only
consider Eq. (3) in the following discussion.
We note here that the above introduced functional form of
the proliferation rate, (φ), is analogous to that consideredin the previous model [13–16] in which they also included
the diffusion of nutrient concentration. One can easily showthat the form of Eq. (5) can be obtained by simply assumingthat the nutrient concentration decreases linearly with thecancer composition φ. For the purpose of clarifying the effectsof cancer proliferation and hydrodynamic interactions, it issufficient to consider the above sigmoidal growth without
introducing any additional field variable.
It should be mentioned that the above logistic growth ofcancer cells can also originate from the mechanical coupling
between the net cell division rates and pressure [4]. In general,the cell division rates depend on mechanical pressure [37–40]
and are characterized by the homeostatic pressure, i.e., thepressure for which cell division and apoptosis balance and nonet growth occurs. Near the homeostatic state, we are allowed
to expand both the pressure and the net cell division rate
to linear order in density difference around the homeostatic
density [4]. Such an effect also leads to the growth term in
Eq. (5).
B. Dynamical equations

Here, a has the dimension of length, β −1 has the dimensionof energy, χ is a dimensionless interaction parameter between
cancer and healthy cells, and κ > 0 is a quantity related to theline tension in the 2D cellular sheet.
Since the above equation has the same form as the FloryHuggins free energy, a phase separation occurs for the condition χ > 2 [45]. Notice that the local terms can be replacedby any other phenomenological description which exhibits aphase separation at sufficiently strong repulsion between the
different cell types. Hence, the exact functional form is notimportant and different forms of free energy were proposed inRefs. [13–16].
The chemical potential μ is obtained from the functional
derivative of the total free energy F with respect to φ [45],


1
δF
φ
= 2 ln+ χ (1 − 2φ) − κ∇ 2 φ. (9)
μ=
δφ
a β
1−φ

Next we consider the time evolution equations for φ and v.
By extending the standard model H for phase separations withhydrodynamic interactions [18–20], the dynamical equations
that are consistent with Eq. (3) can be given by

The stress tensor due to the composition gradient  is given
by [45]
∂φ ∂φ
,
(10)
i j = −κ∂ri ∂r j∂φ= −∇ · (φv) + L∇ 2 μ + (φ),∂t∂vρ= η∇ 2 v − ∇ p + ∇ ·  − ζ v,∂twith i, j = x, y.
The coupled Eqs. (6) and (7) together with the incompressibility condition in Eq. (1) constitute our model for skincancers and provide us with a new type of phase separation
dynamics. In the absence of the cancer proliferation effect,i.e., γ = 0, the above model reduces to conventional models for macrophase separations [18–20]. When γ = 0, ourmodel reduces to model H in the limit of ζ → 0 with fullhydrodynamic interactions, while it corresponds to model B
in the limit of ζ → ∞ for which hydrodynamic interactions
are completely suppressed. The case of γ = 0 showing an arrested phase separation was studied for the pattern formation

(6)
(7)

together with the incompressibility condition in Eq. (1). In
the above equations, L is the transport coefficient, μ is thechemical potential, ρ is the mass density, η is the viscosity, p isthe 2D pressure,  is the stress tensor due to the composition
gradient, and ζ is the friction coefficient. For simplicity, we
assume that both ρ and η are constants and do not dependon φ. Moreover, we consider the case when the transport

032416-3

HOSHINO, LIU, WU, CHEN, TSURUYAMA, AND KOMURAPHYSICAL REVIEW E 99, 032416 (2019)

FIG. 2. Time evolutions of cancer area fraction φ(r, t ) for four different values of the cancer proliferation rate γ = 1, 3, 4, and 5 × 10−3(bottom to top) in the presence of full hydrodynamic interactions (ζ = 0). The other dimensionless parameters are φ0 = 0.3, φ∞ = 0.8, χ =2.5, κ = 1, ρ = 0.3, and η = 1.0. The system size is 512 × 512 and the velocity filed is not shown. In the present grayscale representation,
the values 0 and 1 correspond to white and black, respectively. For γ = 1 and 5 × 10−3 , see also SM1.mp4 and SM2.mp4, respectively, in the
Supplemental Material [47].

of bacterial colonies in the absence of hydrodynamic interactions [21].
C. Simulation method

We numerically solve Eqs. (1), (6), and (7) by using astandard Euler’s method on a 2D square lattice of size 512 ×
512 with periodic boundary conditions. The pressure field p is
calculated with the marker-and-cell method in each time step[46]. It is convenient to use the quantities a, β −1 , and a4 β/L
to scale length, energy, and time, respectively. The numericalestimations for these quantities will be discussed in Sec. V.
Then the dimensionless velocity becomes v ≡ (a3 β/L)v andthe dimensionless model parameters are defined by
ρ≡

L2
L
a4 β
ρ, η ≡ 2 η, γ ≡γ , ζ ≡ Lζ , κ ≡ βκ.
4
a β
a
L
(11)

With the above rescaling, we end up with the following
six dimensionless parameters: χ , κ, ρ, η, γ , and ζ . Amongthese parameters, we have fixed four of them as χ = 2.5,
κ = 1.0, ρ = 0.3, and η = 1.0 in our simulations. Moreover,the initial and the final values of the cancer area fractions
are chosen as φ0 = 0.3 and φ∞ = 0.8 [36], respectively.
In the following, we shall mainly vary the two dynamicalparameters, γ and ζ , to see the effects of cancer proliferation and hydrodynamic interactions on the pattern formationof skin cancers. Physically speaking, the strength of thehydrodynamic interaction should be characterized by a dimensionless number ζ a2 /η = ζ /η that involves both the viscosity and the friction coefficient. Since we set η = 1.0 in
our simulations, the parameter ζ controls the strength of thehydrodynamic interaction. When we present the simulation
results in Sec. III, the above tilde notation is omitted and all
the quantities are treated as dimensionless numbers.
III. SIMULATION RESULTS
A. Pattern formation dynamics

In this section, we present the results of the numerical simulations of the proposed model. We first define the spatiallyaveraged composition of cancer cells as
1
φ(t ) =dr φ(r, t ),(12)A
where A is the total area of the system. Because of the cancer
proliferation, φ(t ) varies from the initial value φ0 = 0.3
toward the stationary value φ∞ = 0.8. Typical time evolutionsof cancer patterns are shown in Fig. 2 when ζ = 0 for four
different values of the cancer proliferation rate γ = 1, 3, 4,
and 5 × 10−3 (bottom to top). Notice that ζ = 0 corresponds
to the case with full hydrodynamic interactions.
Let us first discuss the case of small proliferation rate γ =1 × 10−3 (bottom panels in Fig. 2 and the movie SM1.mp4
in the Supplemental Material [47]). In the initial stage ataround t = 5 × 102 , dots of cancer cells (shown in black) are

032416-4

PATTERN FORMATION OF SKIN CANCERS: EFFECTS OF …

PHYSICAL REVIEW E 99, 032416 (2019)

FIG. 3. Time evolutions of cancer area fraction φ(r, t ) for four different values of the friction coefficient ζ = 10−3 , 10−2 , 10−1 , and ∞ (topto bottom) while the cancer proliferation rate is fixed to γ = 1 × 10−3 . Notice that the limit ζ → ∞ is equivalent to the complete absence ofhydrodynamic interactions (No HI). In practice, such a situation was simulated by omitting the advection term in Eq. (6). The other parameters
are the same as those in Fig. 2. The values 0 and 1 correspond to white and black, respectively. For No HI, see also SM3.mp4 in the Supplemental
Material [47].

formed within a continuous healthy region (shown in white).
We shall call such a structure as a “cancer-in-healthy” (C/H)
pattern. As time evolves, smaller cancer domains collide andmerge to form larger domains at around t = 104 . However,
not all the cancer domains are connected to each other eventhough φ(t ) already exceeds the critical composition φc =0.5. The C/H pattern in the late stage no longer evolves intime and the system attains a steady state without undergoinga macroscopic phase separation. This result shows that ourmodel exhibits a microphase separation.
When the cancer proliferation rate is larger such as whenγ = 5 × 10−3 (top panels in Fig. 2 and the movie SM2.mp4
in the Supplemental Material [47]), healthy regions transformto cancer domains even in the early stage, and the C/Hpattern is already formed at around t = 102 . As the averagecomposition φ(t ) increases, a locally bicontinuous cancer
structure is formed at around t = 103 . However, such a locallybicontinuous structure is destroyed later and smaller healthy
domains emerge. At this stage, black cancer domains arealmost fully connected to form a large continuous domain
at around t = 5 × 103 . In the late stage, circular domainsof healthy cells appear in the network of cancer cells. Sucha structure will be called as a “healthy-in-cancer” (H/C)
pattern. These circular healthy domains do not coarsen any
more in the long time and result in a microphase separation.
When the proliferation rate is intermediate such as whenγ = 3 × 10−3 , healthy domains are elongated and form anarrow continuous network. Moreover, cancer domains in thelate stage at around t = 105 take polygonal shapes rather thancircular shapes. For γ = 4 × 10−3 , a coexistence between the
C/H and H/C patterns is observed as a steady-state structure.
So far we have explained the effects of cancer proliferation
rate γ in the presence of full hydrodynamic interactions,
i.e., ζ = 0. Next we investigate the hydrodynamic effects
by changing the friction coefficient ζ . In Fig. 3, we present
the time evolutions of cancer patterns when the proliferation
rate is fixed to γ = 1 × 10−3 while the friction coefficientis varied as ζ = 10−3 , 10−2 , 10−1 , and ∞ (top to bottom).
Notice that hydrodynamic interactions are completely absent
when ζ → ∞. In practice, this situation is simulated byomitting the advection term in Eq. (6) which is then decoupledfrom the Stokes equation. When the friction coefficient issmall such as when ζ = 10−3 (top panels in Fig. 3), thetime evolution of cancer pattern is similar to that obtainedwith full hydrodynamic interactions (bottom panels in Fig. 2).
However, the steady-state cancer domains at around t = 105
are more elongated. The appearance of elongated domains inthe steady state is more remarkable for ζ = 10−2 .
As the hydrodynamic interactions are further weakened
such as when ζ = 10−1 , cancer domains are more elongatedespecially in the late stage patterns. Here we emphasize againthat the major cancer domains are disconnected while theminor healthy domains form a continuous network structure.
When hydrodynamic interactions are completely absent (bottom panels in Fig. 3 and the movie SM3.mp4 in the Supplemental Material [47]), we eventually obtain an asymmetric032416-5

HOSHINO, LIU, WU, CHEN, TSURUYAMA, AND KOMURAPHYSICAL REVIEW E 99, 032416 (2019)

FIG. 5. Steady-state diagram of cancer patterns obtained for
different cancer proliferation rate γ and friction coefficient ζ (controlling the strength of hydrodynamic interactions). The other parameters are the same as those in Fig. 2. Hydrodynamic interactions
are fully present when ζ = 0, whereas they are completely absent in
the limit of ζ → ∞ (No HI). The latter situation was simulated byomitting the advection term in Eq. (6). Black circles correspond tocancer-in-healthy (C/H) patterns (such as the bottom right pattern
in Fig. 2), red (light gray) circles correspond to healthy-in-cancer(H/C) patterns (such as the top right pattern in Fig. 2), and green(open) circles correspond to (locally) asymmetric bicontinuous (AB)patterns (such as the bottom right pattern in Fig. 3) in the respective
steady states. Black triangles indicate the coexistence between C/H
and H/C patterns (such as γ = 4 × 10−3 and t = 105 in Fig. 2).

FIG. 4. Plots of the average area fraction φ(t ) , defined by
Eq. (12), as a function of time t (a) in the absence of hydrodynamic
interactions and (b) in the presence of full hydrodynamic interactions(ζ = 0). In the former case, simulations were performed by omitting
the advection term in Eq. (6). The cancer proliferation rate is changedas γ = 1, 2, 3, 4, and 5 × 10−3 (from bottom to top). The otherparameters are the same as those in Fig. 2.

bicontinuous (AB) structure at least locally. In this structure,both the wider interconnected cancer domain and the narrowerinterconnected healthy domain are convoluted to each otherfor t  104 .
B. Average cancer composition

In Fig. 4, we have plotted the average cancer composition
φ(t ) , defined by Eq. (12), as a function of time t by changing
the cancer proliferation rate γ . To calculate this quantity,
average over five independent runs (starting from different
initial configurations) has been taken. Figure 4(a) is the casewhen hydrodynamic interactions are completely absent. As
γ is increased, the saturation time becomes smaller and thesaturated value of φ(t ) becomes larger. It is interesting tonote that φ(t ) overshoots before it reaches the stationary
value.
When hydrodynamic interactions are fully present (ζ =0), however, the time evolutions of φ(t ) are different aspresented in Fig. 4(b). Here we notice that the value of φ(t )becomes slightly larger when the hydrodynamic interactions
are present especially for larger γ values. However, the overshooting behavior of φ(t ) is suppressed in Fig. 4(b). These
results indicate that hydrodynamic interactions affect not only
the steady-state behavior but also the transient dynamics ofpattern formation.
C. Steady-state diagram

Next we have systematically varied the proliferation rate γ
and the friction coefficient ζ to see how the steady-state structures depend on these dynamic parameters. We have mentioned before that there are at least three different steady-statepatterns: cancer-in-healthy (C/H), healthy-in-cancer (H/C),
and asymmetric bicontinuous (AB) patterns. The obtained
steady-state patterns are classified into these three cases fordifferent combinations of γ and ζ . In Fig. 5, we summarize
the results in terms of a steady-state diagram in which thethree different cases, C/H, H/C, and AB are distinguished.
The triangle indicates the coexistence between C/H and H/Cpatterns.

032416-6

PATTERN FORMATION OF SKIN CANCERS: EFFECTS OF …

PHYSICAL REVIEW E 99, 032416 (2019)

The C/H pattern clinically corresponds to the globule
pattern of melanoma, and is typically observed when the proliferation rate γ is small and hydrodynamic interactions arestrong (small ζ ). The AB pattern appears when hydrodynamic
interactions are weak or fully suppressed (large ζ ) while theproliferation rate γ is relatively small. The AB pattern may
correspond to the stripe pattern of melanoma mainly found in
human palms or soles. Finally, the H/C pattern typically appears when both γ and ζ are large. When the proliferation rateis as large as γ = 5 × 10−3 , only the H/C pattern is obtained
irrespective of the strength of hydrodynamic interactions. In
contrast to the other two cases, however, the H/C pattern isusually not diagnosed in typical skin cancers because domainsof healthy cells are completely destroyed by invasive cancer
cells.
In the case of an ordinary microphase separation, the latestage structure should be the H/C pattern when φ(t ) > 0.5.
As shown in Fig. 5, however, we obtain either the C/H patternor the AB pattern for different combinations of ζ and γ ,especially when γ is small. This is one of the unique featuresof the proposed model for cancer cells with hydrodynamic
interactions. Since these steady-state patterns are typicallyobtained in the presence of hydrodynamic interactions, weconsider that they appear kinetically and do not correspondto equilibrium microstructures. The different mechanisms forthe pattern formation will be discussed in Sec. IV.
D. Structure analysis

To analyze the time evolutions of the patterns quantitatively, we have calculated their structure factors. Let δφ(r, t )be the deviation of φ(r, t ) from its average value, δφ(r, t ) =φ(r, t ) − φ(t ) , where φ(t ) defined in Eq. (12) depends
on time. First we introduce the spatial Fourier transform ofδφ(r, t ) by(13)
δφk (t ) = dr δφ(r, t )e−ik·r ,where k = (kx , ky ) is a 2D wave vector. Then the structure
factor is defined as
S(k, t ) = δφk (t )δφ−k (t ) ,(14)

where the average is over the ensemble of systems. Using thecircularly averaged structure factor S(k, t ) with k = |k|, we
calculate the following (inverse) characteristic length scale ofpatterns [48],k(t ) =−1

dk k S(k, t ),dk k −2 S(k, t )(15)where we omit k = 0 in the integrals.
In Fig. 6, we plot the time evolutions of the circularly
averaged structure factor S(k, t ) as a function of the wave
number k when γ = 3 × 10−3 . Figure 6(a) corresponds tothe case when hydrodynamic interactions are absent, while
Fig. 6(b) presents the case with full hydrodynamics. By comparing these two cases, we see that the early stage structuresare similar as long as the proliferation rate γ is the same. In
the intermediate stage, however, the microstructure formation
is faster in the presence of hydrodynamic interactions, andthe peak position is shifted to a smaller k value in Fig. 6(b).

FIG. 6. Plots of the circularly averaged structure factor S(k, t )as a function of the wave number k for different time steps t (fromright to left) (a) in the absence of hydrodynamic interactions and(b) in the presence of full hydrodynamic interactions (ζ = 0). The
cancer proliferation rate is fixed to γ = 3 × 10−3 , while the otherparameters are the same as those in Fig. 2. Notice that the real spacepattern evolution that corresponds to (b) is presented in Fig. 2.

We also find that the peak height in the late stage is slightly
smaller in Fig. 6(b) than that in Fig. 6(a).
In Fig. 7, we have plotted the characteristic wave numberk(t ) , defined by Eq. (15), as a function of time. The proliferation rate γ is similarly changed as in Fig. 4 and the averageover five independent runs has been taken as before. As shownin Fig. 7(a) when hydrodynamic interactions are absent, theaverage wave number k(t ) saturates at larger values (smaller
structures) when γ is increased. This means that γ is an important parameter that controls the characteristic length scaleof the steady-state microstructures. Comparing Figs. 4(a) and
7(a), we notice that the saturation times for φ(t ) roughly
correspond to those for k(t ) .
The effects of hydrodynamic interactions on k(t ) can beseen in Fig. 7(b) for which we have set ζ = 0. Here k(t )032416-7

HOSHINO, LIU, WU, CHEN, TSURUYAMA, AND KOMURAPHYSICAL REVIEW E 99, 032416 (2019)

FIG. 7. Log-log plots of the characteristic wave number k(t ) ,defined by Eq. (15), as a function of time t (a) in the absence ofhydrodynamic interactions and (b) in the presence of full hydrodynamic interactions (ζ = 0). The cancer proliferation rate is changedas γ = 1, 2, 3, 4, and 5 × 10−3 [from bottom to top at t = 105 in(a) and from right to left for the intermediate time region in (b)]. The
other parameters are the same as those in Fig. 2. The dashed lines
indicate the power-law behaviors with the respective slopes −1/3 in(a) and −2/3 in (b).

shows a large decrease up to the intermediate stage. This
result indicates that hydrodynamic interactions tend to form
larger domains even though they are only transient structures.
Interestingly, a minimum of k(t ) appears at around t = 104
and k(t ) exhibits an undershooting behavior. Hence, thetransient domain size depends not only on the proliferation
rate γ but also on the friction coefficient ζ . According to
Fig. 7(b), however, the late stage dynamics of k(t ) has not
yet reached the steady state completely. Such a long-lived dynamics is also different from the case without hydrodynamic
interactions.

FIG. 8. (a) Log-log plot of the steady-state value of the characteristic wave number k∞ in Fig. 7(a) as a function of the proliferation
rate γ in the absence of hydrodynamic interactions. From the slopeof the fitted straight line, we find a power-law relation k∞ ∼ γ 0.32 .
(b) Log-log plot of k(t ) t 1/3 as a function of the dimensionlessvariable γ t using all the data in Fig. 7(a). The collapse of all thedata confirms the validity of the scaling assumption in Eq. (16). The
dashed lines indicate the power-law behaviors with the respective
slopes 1/3 both in (a) and (b).

E. Scaling analysis of domain growth

The result in Fig. 7(a) can be further analyzed in termsof the scaling argument. In the case of γ = 0, for which theaverage cancer fraction remains constant (conserved case), thesystem exhibits a macrophase separation because such a casewithout any hydrodynamics corresponds to model B. Let us
denote the steady-state characteristic wave number as k∞ . In
Fig. 8(a), we plot k∞ as a function of γ used in Fig. 7(a). We
find a clear power-law behavior, i.e., k∞ ∼ γ s with s ≈ 0.32.
This result suggests that the characteristic wave number obeys

032416-8

PATTERN FORMATION OF SKIN CANCERS: EFFECTS OF …

the following scaling form:k(t ) ∼ t −α f (γ t ),(16)

where α is the domain growth exponent in the absence of thecancer proliferation effect, and f (z) is a scaling function witha dimensionless variable z = γ t. A similar scaling hypothesis
was successfully used to analyze the phase separation dynamics of chemically reactive binary mixtures [31,32] or that ofblock copolymer melts [33–35].
According to the evaporation-condensation process considered by Lifshitz and Slyozov [30], the growth exponent shouldbe α = 1/3 when hydrodynamic interactions are absent. This
exponent is indeed observed and shown by the dashed line in
Fig. 7(a) before the saturation time. The asymptotic behaviorof the scaling function should be f (z) ∼ const. for z
1,
and f (z) ∼ zα for z
1. The latter power-law behavior isrequired because k(t ) should not depend on time t in thesteady state. Hence, we immediately obtain k∞ ∼ γ α ands = α. In Fig. 8(b), we have replotted the quantity k(t ) t 1/3
as a function of γ t using all the data in Fig. 7(a). The collapseof all the curves demonstrates that our simulation results arein good agreement with the above scaling ansatz as long as γ
is small enough.
In Fig. 7(b) with full hydrodynamic interactions, thegrowth exponent in the intermediate stage is as large as α =2/3 which is much larger than that in Fig. 7(a). However,this result does not obey a simple scaling behavior becauseof the complicated undershooting behaviors. Here we pointout that the value α = 2/3 was discussed by Furukawa who
considered the interplay between the inertia of the fluid andthe surface energy density [20,49]. This growth exponent wasalso confirmed by lattice Boltzmann simulations for a criticalquench of a 2D binary fluid when the viscosity is small andstochastic noise is absent [50,51]. Our result cannot be directlycompared with theirs because the average composition varies
with time and also the system exhibits a microphase separation in the late stage. However, it is evident from Fig. 7(b) thata substantial acceleration of phase separation takes place inthe presence of hydrodynamic flows.

PHYSICAL REVIEW E 99, 032416 (2019)

The amplitude equations can be derived by substituting
Eqs. (17) and (18) into Eqs. (6) and (7):


γ
φ
dφ
−
|δφq |2 ,(19)≈ γ φ 1 −
dt
φ∞
φ∞ q
dδφq
≈ −iq ·
(δφq1 vq2 ) − q2 μq
dt
q1 +q2 =q


2φ
γ
−
+ γ δφq 1 −(δφq1 δφq2 ),φ∞
φ∞ q +q =q
1

(20)
ρ

dvq
= −ηq2 vq − iqpq + iq · q − ζ vq ,dt1

q = κ

(q1 ⊗ q2 )δφq1 δφq2 ,

(23)
(24)

q1 +q2 =q

where q̂ ≡ q/q is the unit vector and ⊗ represents the dyadic
product. Then the previous amplitude equations can be simplified asdδφq≈ λ1 δφq + λ2
δφq1 δφq2
dt
q +q =q2+ λ3

δφq1 δφq2 δφq3
q1 +q2 +q3 =q

In the early stage of phase separation, when φ is mostly
uniform with small perturbations, we are able to analyze thepattern formation by using the amplitude equations methodwith which the variations of φ and v are viewed as a group of
perturbation waves:

(17)

q

vq (t )eiq·r + c.c.,3

pq = q̂ · q · q̂,

A. Early stagev(r, t ) ≈

2

(22)

IV. MECHANISMS FOR PATTERN FORMATIONδφq (t )eiq·r + c.c. ,(21)

where q = |q|. In the above, μq , pq , and q are the qth
component of the Fourier series of μ, p, and , respectively,
and are given by


1
1
2
+
− 2χ + κq δφq
μq ≈φ
1 − φ


1
1
1
−
+
+
δφq1 δφq22
φ 2(1 − φ )2 q +q =q
1
2


1 1
1
+
+
δφq1 δφq2 δφq3 ,3 φ 3(1 − φ )3 q +q +q =q1

φ(r, t ) ≈ φ(t ) +

2

(18)

q

where c.c. denotes the complex conjugate and the summationof q is taken over the principal modes of the pattern of interest.

− iq ·(δφq1 vq2 ),(25)

q1 +q2 =q

ρdvq= −(ηq2 + ζ )vq + iq · q · (I − q̂ ⊗ q̂),
dt

(26)

where I is the unit tensor and the three coefficients in Eq. (25)are given by




1
2φ
1
λ1 = −q2,+− 2χ + κq2 + γ 1 −φ
1 − φ
φ∞
(27)


γ
1
q2
1
λ2 = −
−
−
+
,
(28)
2
φ 2
(1 − φ )2
φ∞


q2 1
1
+
λ3 = −
< 0.
(29)
3
3 φ
(1 − φ )3032416-9

HOSHINO, LIU, WU, CHEN, TSURUYAMA, AND KOMURAPHYSICAL REVIEW E 99, 032416 (2019)

According to the above amplitude equations, it is clearthat the hydrodynamic interaction, described by the last termin Eq. (25), is a higher order contribution which does not
influence the early stage dynamics. The linear term λ1 δφq
in Eq. (25) is independent of hydrodynamic interactions anddominates when δφq is small. Therefore, the early stage
dynamics must be similar regardless of the values of ζ as seenin Fig. 3 for t  103 .
According to the time evolution of the velocity in Eq. (26),the combination ηq2 + ζ controls the decay of the hydrodynamic flow. This implies that hydrodynamic interactions playa significant role for large length scales. Consequently, theflow is suppressed in the early stage when the average wave
number k(t ) is large, whereas it is strengthened when kdecreases as pattern evolves.
B. Late stage

To discuss the late stage dynamics from the viewpoint ofhydrodynamic flows, we show in Fig. 9 the velocity fieldv(r, t ) together with the cancer fraction field φ(r, t ) in thepresence of full hydrodynamic interactions (ζ = 0) at (a) t =7800 and (b) t = 9400 when γ = 1 × 10−3 (see also SM1.mp4
in the Supplemental Material [47]). In Fig. 9(a), a large scalepair of vortices is created; one of them rotates clockwise andthe other moves counterclockwise. Such a flow is triggered by
the coalescence of two smaller domains into a larger domain.
As a result, a strong flow is induced at the neck region of
the two merging domains. Somewhat later in Fig. 9(b), however, a circular flow appears inside a large domain. Anotherimportant feature in this pattern is the existence of a flow
along the domain boundaries. Such a flow sometimes induces
a large velocity field in the narrow channel between larger
domains.
In the late stage of pattern evolution, the domain structuresof the healthy and cancer cells become relatively robust. The
values of φ within healthy-rich and cancer-rich domains are
saturated to φ ≈ 0.145 and φ ≈ 0.855, respectively, whichcorrespond to the two free energy minima of Eq. (8) when
χ = 2.5. Once the microstructure is formed, the subsequentevolution of pattern is determined by the competition between
two different processes; the shape accommodation and thecoalescence process.
The shape accommodation results from the movementof interfaces that tends to minimize total interfacial energy.
Therefore, the system energetically favors circular domains
and the resultant pattern is the C/H pattern composed of
circular cancer domains separated by healthy cells. However, the interface is not static due to non-zero net proliferation rate and the coalescence occurs when two nearby
cancer domains continue to grow and eventually connect
each other. For larger proliferation rates, the coalescencesurpasses the shape accommodation process. Hence, cancer
domains get interconnected and the length scale of patternincreases. This process leads to a breakdown of the sixfold
symmetry of the C/H pattern owing to the random connecting
processes.
Since the pattern is kinetically controlled by these twoprocesses, the steady state should depend on the values ofγ and ζ , as summarized in Fig. 5. The rate of coalescence

FIG. 9. Plots of the velocity field v(r, t ) shown by the arrows
at (a) t = 7800 (system size 200 × 200) and (b) t = 9400 (systemsize 150 × 150) when γ = 1 × 10−3 in the presence of full hydrodynamic interactions (ζ = 0). The other dimensionless parameters are
φ0 = 0.3, φ∞ = 0.8, χ = 2.5, κ = 1, ρ = 0.3, and η = 1.0. Bothpatterns are the closeups of a larger system size simulation aspresented by the bottom panels of Fig. 2. See also SM1.mp4 in the
Supplemental Material [47].

process is influenced by the domain growth rate γ . The shapeaccommodation is realized through the mass transportationand it is enhanced by the additional hydrodynamics flows
across interface, as presented in Fig. 9. Thus, when γ is aslarge as γ ≈ 5 × 10−3 , the dominating coalescence processconnects all domains together and transforms the pattern into a
uniform cancer cells with few healthy spots, corresponding to
the H/C patterns. However, the shape accommodation process
is faster than the coalescence process for small γ so that theC/H pattern is preserved in the late stage. The intermediate
stripe-like pattern (AB pattern) appears in the steady state
when the coalescence and shape accommodation processes
are comparable.

032416-10

PATTERN FORMATION OF SKIN CANCERS: EFFECTS OF …
V. SUMMARY AND DISCUSSION

In this paper, we have performed numerical simulationsof pattern formation of skin cancers. In our phase separation
model for a binary cellular system, we have taken into accountthe effects of cancer proliferation and hydrodynamic interactions to describe the time evolutions of cancer cells. As a resultof the proliferation effect, the emerging patterns drasticallychange their structures depending on the different stages ofthe phase separation dynamics.
By controlling the cancer proliferation rate γ and thefriction coefficient ζ between dermis and epidermis, we haveobtained various types of steady-state cancer pattern such asa cancer-in-healthy pattern (C/H), a healthy-in-cancer pattern(H/C), and an locally asymmetric bicontinuous (AB) structure. As summarized in Fig. 5, we have constructed the steadystate pattern diagram for different combinations of γ and ζ
values. In particular, the C/H patterns obtained for a small
proliferation rate and strong hydrodynamic interactions (small
ζ ) and the AB structures obtained for weak hydrodynamic
interactions (large ζ ) might correspond to the globule and thestripe patterns, respectively, in real melanoma diagnoses.
For a quantitative analysis, we have calculated the spatiallyaveraged composition of cancer cells, φ(t ) , and the characteristic length of the cancer patterns, k(t ) , as a functionof time t (see Figs. 4 and 7) both in the presence and theabsence of hydrodynamic interactions. We have shown thatφ(t ) and k(t ) depend not only on the proliferation rate but
also on the strength of hydrodynamic interactions. Withouthydrodynamic flows, we have confirmed in Fig. 8 that thescaling behavior of the characteristic length is described by
the form of Eq. (16). With hydrodynamic flows, however,the domain growth exponent in the intermediate stage was aslarge as α = 2/3, showing a pronounced acceleration of the
microphase separation.
First, we shall give some numbers for the quantities mentioned in Sec. II C to scale length, energy, and time that are
relevant to skin cancers [see Eq. (11)]. The typical length scaleobserved in skin cancer patterns is in the order of 10−3 m.
According to Fig. 7, the characteristic wave number in thesteady state of our simulation is k a ≈ 0.1 (notice that werecover the dimensions of the physical quantities in thissection). From these values, we set the unit of length asa ≈ 10−5 m which corresponds to the size of an epidermal
cell [13]. Since the interstitial fluid pressure in skin carcinoma
was estimated to be roughly  ≈ 103 Pa [13,36], we obtainthe typical energy scale as β −1 ∼ a3 ≈ 10−12 J, which ismuch larger than the thermal energy. From the data of theinterphase friction [13,52,53], the 3D transport coefficient canbe evaluated as L3D ≈ 10−15 m2 Pa−1 s−1 . With this value,
we estimate the typical time scale in our model as a4 β/L ∼a5 β/L3D ≈ 102 s.
Having discussed various scales for skin cancers, we canconvert the dimensionless parameters in our simulations to
the physical quantities with dimensions. For example, the dimensionless time t/(a4 β/L) ≈ 105 to reach the steady statesin Fig. 7 roughly corresponds to t ≈ 102 days which arereasonable for cancer spreading. The choice η = Lη/a2 = 1
in our simulation corresponds to η3D ∼ η/a ≈ 105 Pa s thatfits within the previously reported viscosity values [4,25]. As

PHYSICAL REVIEW E 99, 032416 (2019)

for the cancer proliferation rate, the value γ ∼ a4 βγ /L =10−3 roughly corresponds to γ ≈ 10−5 s−1 ≈ 1 day−1 . This
proliferation rate is in agreement with that in the previous
reports [13,54]. Finally, the range of the scaled friction coefficient ζ = Lζ = 10−3 − 1 in our simulation predicts ζ ≈107 − 1010 Pa s m−1 and it coincides with the range of thefriction coefficient in Ref. [55].
Next we discuss the role of cancer proliferation effects on
the phase separation dynamics. In the conventional model B
describing ordinary macrophase separations, a typical timescale is set by the transport coefficient L. In the present
model, however, the proliferation rate γ in Eq. (5) provides
us with additional time scale. Generally speaking, the phaseseparation dynamics should be determined by the competition between these two time scales. In our simulation, theinitial cancer composition started from φ0 = 0.3 and L wasmuch larger than γ . More precisely, we have chosen thedimensionless number as a4 βγ /L ≈ 10−3 in the simulations[see Eq. (11)]. Hence, the compositional instability for thephase separation, that is governed by L, takes place before
the average composition φ(t ) increases with the rate γ .
As shown in Fig. 2, the cancer domains appear as a resultof unstable concentration fluctuations, and they form C/Hpatterns for φ(t ) < 0.5 in the early stage. In the late stage,the initial C/H pattern continues to remain for smaller γ
values, while it transforms into the H/C pattern for larger
γ values. When the quantity a4 βγ /L is much larger andbecomes close to unity, the system always exhibits the H/Cpattern because the average composition will be immediately
saturated at a larger value φ(t ) > 0.5 before the system
undergoes a phase separation. Hence, the cancer proliferation
significantly affects the microstructures of cancer patterns.
In the present work, we have considered a 2D systemcomposed of cancer and healthy cells whose compositions
evolve in time due to the cancer proliferation effect. Although a similar model was proposed by Chatelain et al.
[13,14], the main difference in our work is that the effects ofhydrodynamic interactions are explicitly taken into account.
Moreover, the strength of hydrodynamic interactions can becontrolled by changing the friction coefficient ζ . When hydrodynamic interactions are fully present, the C/H patterns
continue to remain even in the late stage when φ(t ) > 0.5(see bottom panels in Fig. 2 and SM1.mp4 in the Supplemental
Material [47]). Such a transient pattern was not observed in
the previous study by Chatelain et al. [13,14].
Alternatively, Chatelain et al. took into account the diffusion of nutrient (oxygen) concentration chosen as an additional variable [13,14]. Accordingly, they employed a diffusion equation for the nutrient concentration with a source
term. In their model, the cell-nutrient interaction defines atypical diffusive length that controls the saturation of growing
domains. In our model, we did not consider such a coupling
to the diffusion of nutrients from an outer environment, but
simply used the logistic growth model to describe the cancer
proliferation [see Eq. (5)]. As mentioned before, this simplification is justified when the cancer composition is proportional
to the nutrient concentration.
We have assumed that dermal and epidermal boundaryis flat and the epidermal layer was modeled as a 2D fluid.

032416-11

HOSHINO, LIU, WU, CHEN, TSURUYAMA, AND KOMURAPHYSICAL REVIEW E 99, 032416 (2019)

However, the structure of dermis and epidermal can affectthe cell differentiation and also the cancer pattern formation.
For example, Balois et al. considered melanin transport inepidermis and showed that it is influenced by the dermal andepidermal shape [16]. Such a geometrical effect of basal layer
will be considered in our future study by taking into accountthe hydrodynamic interaction.
Cates et al. argued that the appearance of an arrested
phase separation in bacterial colonies can be explained only
by considering a local density-dependent motility and thebirth/death of bacteria [21]. In their work, the competition
between the effects of birth/death and diffusion leads to atypical length scale beyond which domain coarsening does notoccur. The obtained patterns of 2D simulation indeed show
droplets of the high-density phase dispersed in a continuous
low-density phase at large times [21]. Such a situation is very
reminiscent to the results of our model in the absence hydrodynamic interactions (either C/H or H/C pattern). However,
we have shown that hydrodynamic interactions affect not only
the steady-state patterns but also the transient patterns.
In Sec. II A, we have mentioned that the logistic growth ofcancer cells in Eq. (5) can stem from the mechanical coupling
effect that is controlled by the homeostatic pressure [4]. Ranftet al. discussed the propagation of an interface between two
different cell populations when the homeostatic pressures of
two cell types are different [56]. Taking into account bothsubstrate friction and hydrodynamic interactions, Podewitz
et al. performed mesoscopic simulations to investigate interface dynamics of competing tissues [57]. They showedthat the propagation velocity of the interface is proportional
to the homeostatic stress difference. Recently, Williamson
and Salbreux studied the stability and roughness of sucha propagating interface [55]. In these studies, however, theformation of microstructures of cancer cells, such as dots or
stripes, has not been investigated.
As mentioned before, our model can reproduce clinicallyobserved globule and stripe patterns in melanoma. The C/Hpatterns tend to appear when the proliferation rate is small

and the hydrodynamic interactions are strong. By contrast, thestripe patterns, which are often found in human palms or soles,
tend to appear when hydrodynamic interactions are absent.
In reality, palms and soles contain a thick stratum corneum
and an unique cell layer called “stratum lucidum” which hasa finite stiffness. Such a stiffness may reduce hydrodynamic
interactions and results in the formation of stripe patterns.
Our model suggests that the proliferation and invasionof cancer cells in superficial spreading melanoma can bepredicted by observing the epidermis using dermoscopy.
Melanoma cells migrate horizontally in the epidermis in theinitial stage of tumor development, during which the clinical staging is described by “Clark’s level” and “Breslow’s
depth” [58]. In its staging, the diffusion range and the cellspreading pattern of melanoma cells are the most importantmeasures for making prognostic predictions, such as the fiveyear patient survival rate [59,60]. The present work presents
objective diagnostic indicators and methodologies for making
prognostic predictions for these patients that can be verified
by dermoscopic image data. We expect that our work will beapplied to the development and evaluation of future clinicaldiagnosis.

[1] Edited by V. Capasso, M. Gromov, and A. Harel-Bellan, Pattern
Formation in Morphogenesis (Springer, Heidelberg, 2013).
[2] N. L. Komarova, Curr. Opin. Oncol. 17, 39 (2005).
[3] P. R. Paździorek, Bull. Math. Biol. 76, 1642 (2014).
[4] M. Basan, T. Risler, J.-F. Joanny, X. Sastre-Garau, and J. Prost,
HFSP J. 3, 265 (2009).
[5] J. Ranft, M. Basan, J. Elgeti, J.-F. Joanny, J. Prost, and F.
Jüicher, Proc. Natl. Acad. Sci. USA 107, 20863 (2010).
[6] M. Basan, J. Prost, J.-F. Joanny, and J. Elgeti, Phys. Biol. 8,026014 (2011).
[7] D. Gonzalez-Rodriguez, K. Guevorkian, S. Douezan, and F.
Brochard-Wyart, Science 338, 910 (2012).
[8] S. Kumar and V. M. Weaver, Cancer Metastasis Rev. 28, 113(2009).
[9] T. S. Deisboeck, Z. Wang, P. Macklin, and V. Cristini, Annu.
Rev. Biomed. Eng. 13, 127 (2011).
[10] I. A. Rodriguez-Brenes, N. L. Komarova, and D. Wodarz,
Trends Ecol. Evol. 28, 597 (2013).

[11] P. Van Liedekerke, M. M. Palm, N. Jagiella, and D. Drasdo,
Comp. Part. Mech. 2, 401 (2015).
[12] C. Chatelain and M. Ben Amar, Eur. Phys. J. Plus 130, 176(2015).
[13] C. Chatelain, T. Balois, P. Ciarletta, and M. Ben Amar, New J.
Phys. 13, 115013 (2011).
[14] C. Chatelain, P. Ciarletta, and M. Ben Amar, J. Theor. Biol. 290,
46 (2011).
[15] T. Balois and M. Ben Amar, Sci. Rep. 4, 3622 (2014).
[16] T. Balois, C. Chatelain, and M. Ben Amar, J. R. Soc. Interface
11, 20140339 (2014).
[17] I. W. Hamley, The Physics of Block Copolymers (Oxford University Press, Oxford, 1998).
[18] P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys. 49, 435(1977).
[19] P. M. Chaikin and T. C. Lubensky, Principles of Condensed Matter Physics (Cambridge University Press, Cambridge, 1995).

ACKNOWLEDGMENTS

We thank R. Okamoto, K. Yasuda, T. Kato, and R. Kurita
for useful discussions. T.H. thanks the hospitality of National
Tsing Hua University and National Central University wherepart of this research was conducted under the Co-Tutorial
Program. T.H. acknowledges the support by Grant-in-Aid for
JSPS Fellows (Grant No. 17J01643) from the Japan Society for the Promotion of Science (JSPS). K.-A.W. and M.W.L. acknowledge the support of the Ministry of Scienceand Technology, Taiwan (Grant No. MOST 105-2112-M007-031-MY3). K.-A.W., M.-W.L., and H.-Y.C. acknowledge
the support from National Center for Theoretical Sciences,
Taiwan. S.K. acknowledges the support by Grant-in-Aid for
Scientific Research (C) (Grant No. 18K03567) from the JSPS.

032416-12

PATTERN FORMATION OF SKIN CANCERS: EFFECTS OF …

PHYSICAL REVIEW E 99, 032416 (2019)[20] A. Onuki, Phase Transition Dynamics (Cambridge University
Press, Cambridge, 2002).
[21] M. E. Cates, D. Marenduzzo, I. Pagonabarraga, and J. Tailleur,Proc. Natl. Acad. Sci. USA 107, 11715 (2010).
[22] M. S. Steinberg, Science 141, 401 (1963).
[23] M. S. Steinberg, Dev. Biol. 180, 377 (1996).
[24] B. He, K. Doubrovinski, O. Polyakov, and E. Wieschaus,
Nature 508, 392 (2014).
[25] W.-T. Yeh and H.-Y. Chen, Phys. Rev. E 93, 052421 (2016).
[26] W.-T. Yeh and H.-Y. Chen, New J. Phys. 20, 053051 (2018).
[27] E. D. Siggia, Phys. Rev. A 20, 595 (1979).
[28] V. M. Kendon, M. E. Cates, I. Pagonabarraga, J.-C. Desplat, and
P. Bladon, J. Fluid Mech. 440, 147 (2001).
[29] K. Binder and D. Stauffer, Phys. Rev. Lett. 33, 1006 (1974).
[30] E. M. Lifshitz and L. P. Pitaevskii, Physical Kinetics (Pergamon
Press, Oxford, 1981).
[31] S. C. Glotzer, D. Stauffer, and N. Jan, Phys. Rev. Lett. 72, 4109
(1994).
[32] J. J. Christensen, K. Elder, and H. C. Fogedby, Phys. Rev. E 54,
R2212 (1996).
[33] Y. Oono and M. Bahiana, Phys. Rev. Lett. 61, 1109 (1988).
[34] F. Liu and N. Goldenfeld, Phys. Rev. A 39, 4805 (1989).
[35] M. Bahiana and Y. Oono, Phys. Rev. A 41, 6763 (1990).
[36] R. K. Jain, Cancer Res. 47, 3039 (1987).
[37] F. Montel, M. Delarue, J. Elgeti, L. Malaquin, M. Basan, T.
Risler, B. Cabane, D. Vignjevic, J. Prost, G. Cappello, and J.-F.
Joanny, Phys. Rev. Lett. 107, 188102 (2011).
[38] F. Montel, M. Delarue, J. Elgeti, D. Vignjevic, G. Cappello, andJ. Prost, New J. Phys. 14, 055008 (2012).
[39] K. Alessandri, B. R. Sarangi, V. V. Gurchenkov, B. Sinha, T. R.
Kießling, L. Fetler, F. Rico, S. Scheuring, C. Lamaze, A. Simon,
S. Geraldo, D. Vignjević, H. Doméjean, L. Rolland, A. Funfak,
J. Bibette, N. Bremond, and P. Nassoy, Proc. Natl. Acad. Sci.
USA 110, 14843 (2013).
[40] M. Delarue, F. Montel, O. Caen, J. Elgeti, J.-M. Siaugue, D.
Vignjevic, J. Prost, J.-F. Joanny, and G. Cappello, Phys. Rev.
Lett. 110, 138103 (2013).

[41] A. Tiribocchi, R. Wittkowski, D. Marenduzzo, and M. E. Cates,
Phys. Rev. Lett. 115, 188302 (2015).
[42] S. Puri, A. J. Bray, and J. L. Lebowitz, Phys. Rev. E 56, 758(1997).
[43] R. Ahluwalia, Phys. Rev. E 59, 263 (1999).
[44] S. M. Wise, J. S. Lowengrub, H. B. Frieboes, and V. Cristini, J.
Theor. Biol. 253, 524 (2008).
[45] M. Doi, Soft Matter Physics (Oxford University, Oxford,2013).
[46] F. H. Harlow and J. E. Welch, Phys. Fluids 8, 2182
(1965).
[47] See Supplemental Material at http://link.aps.org/supplemental/
10.1103/PhysRevE.99.032416 for movies.
[48] A. Shinozaki and Y. Oono, Phys. Rev. E 48, 2622 (1993).
[49] H. Furukawa, Physica A (Amsterdam) 204, 237 (1994).
[50] W. R. Osborn, E. Orlandini, M. R. Swift, J. M. Yeomans, andJ. R. Banavar, Phys. Rev. Lett. 75, 4031 (1995).
[51] G. Gonnella, E. Orlandini, and J. M. Yeomans, Phys. Rev. E 59,
R4741 (1999).
[52] P. M. van Kemenade, J. M. Huyghe, and L. F. A. Douven,
Porous Media 50, 93 (2003).
[53] E. A. Swabb, J. Wei, and P. M. Gullino, Cancer Res. 34, 2814
(1974).
[54] A. A. Creasey, H. S. Smith, A. J. Hackett, K. Fukuyama, W. L.
Epstein, and S. H. Madin, In Vitro 15, 342 (1979).
[55] J. J. Williamson and G. Salbreux, Phys. Rev. Lett. 121, 238102(2018).
[56] J. Ranft, M. Aliee, J. Prost, F. Jülicher, and J.-F. Joanny, New J.
Phys. 16, 035002 (2014).
[57] N. Podewitz, F. Jülicher, G. Gompper, and J. Elgeti, New J.
Phys. 18, 083020 (2016).
[58] W. D. James, T. G. Berger, and D. M. Elston, Andrews’ Diseasesof the Skin: Clinical Dermatology (Saunders Elsevier, Philadelphia, 2006).
[59] J. L. Bolognia, J. L. Jorizzo, and R. P. Papini, Dermatology(Mosby, St. Louis, 2007).
[60] F. Tas and K. Erturk, Mol. Clin. Oncol. 7, 1083 (2017).

032416-13