Pattern formation of skin cancers: Effects of cancer proliferation and hydrodynamic interactions

We study pattern formation of skin cancers by means of numerical simulation of a binary system consisting of 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 some of 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.

Figure 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 is infinitely large and we do not consider any out-of-plane deformation of 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 $\varphi$ and $\psi$, respectively. The two types of cell fill all the available space and satisfy the saturation constraint, i.e., $\varphi +\psi =1$. Further, the local velocities are denoted by ${\mathbf{v}}_{\varphi }$ and ${\mathbf{v}}_{\psi }$ for cancer and healthy cells, respectively. We also take into account a mechanical friction between dermis and epidermis that is characterized by the friction coefficient $\zeta$.

Figure 2

Time evolutions of cancer area fraction $\varphi (\mathbf{r},t)$ for four different values of the cancer proliferation rate $\gamma =1,3,4$, and $5\times {10}^{-3}$ (bottom to top) in the presence of full hydrodynamic interactions ($\zeta =0$). The other dimensionless parameters are ${\varphi }_0=0.3,{\varphi }_{\infty }=0.8,\chi =2.5,\kappa =1,\rho =0.3$, and $\eta =1.0$. The system size is $512\times 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 $\gamma =1$ and $5\times {10}^{-3}$, see also SM1.mp4 and SM2.mp4, respectively, in the Supplemental Material [47].

Figure 3

Time evolutions of cancer area fraction $\varphi (\mathbf{r},t)$ for four different values of the friction coefficient $\zeta ={10}^{-3},{10}^{-2},{10}^{-1}$, and $\infty$ (top to bottom) while the cancer proliferation rate is fixed to $\gamma =1\times {10}^{-3}$. Notice that the limit $\zeta \to \infty$ is equivalent to the complete absence of hydrodynamic 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].

Figure 4

Plots of the average area fraction $\langle \varphi \left(t\right)\rangle$, 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 ($\zeta =0$). In the former case, simulations were performed by omitting the advection term in Eq. (6). The cancer proliferation rate is changed as $\gamma =1,2,3,4$, and $5\times {10}^{-3}$ (from bottom to top). The other parameters are the same as those in Fig. 2.

Figure 5

Steady-state diagram of cancer patterns obtained for different cancer proliferation rate $\gamma$ and friction coefficient $\zeta$ (controlling the strength of hydrodynamic interactions). The other parameters are the same as those in Fig. 2. Hydrodynamic interactions are fully present when $\zeta =0$, whereas they are completely absent in the limit of $\zeta \to \infty$ (No HI). The latter situation was simulated by omitting the advection term in Eq. (6). Black circles correspond to cancer-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 $\gamma =4\times {10}^{-3}$ and $t={10}^5$ in Fig. 2).

Figure 6

Plots of the circularly averaged structure factor $S(k,t)$ as a function of the wave number $k$ for different time steps $t$ (from right to left) (a) in the absence of hydrodynamic interactions and (b) in the presence of full hydrodynamic interactions ($\zeta =0$). The cancer proliferation rate is fixed to $\gamma =3\times {10}^{-3}$, while the other parameters are the same as those in Fig. 2. Notice that the real space pattern evolution that corresponds to (b) is presented in Fig. 2.

Figure 7

Log-log plots of the characteristic wave number $\langle k\left(t\right)\rangle$, defined by Eq. (15), as a function of time $t$ (a) in the absence of hydrodynamic interactions and (b) in the presence of full hydrodynamic interactions ($\zeta =0$). The cancer proliferation rate is changed as $\gamma =1,2,3,4$, and $5\times {10}^{-3}$ [from bottom to top at $t={10}^5$ 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).

Figure 8

(a) Log-log plot of the steady-state value of the characteristic wave number $k_{\infty }$ in Fig. 7 as a function of the proliferation rate $\gamma$ in the absence of hydrodynamic interactions. From the slope of the fitted straight line, we find a power-law relation $k_{\infty }\sim {\gamma }^{0.32}$. (b) Log-log plot of $\langle k\left(t\right)\rangle t^{1/3}$ as a function of the dimensionless variable $\gamma t$ using all the data in Fig. 7. The collapse of all the data 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).

Figure 9

Plots of the velocity field $\mathbf{v}(\mathbf{r},t)$ shown by the arrows at (a) $t=7800$ (system size $200\times 200$) and (b) $t=9400$ (system size $150\times 150$) when $\gamma =1\times {10}^{-3}$ in the presence of full hydrodynamic interactions ($\zeta =0$). The other dimensionless parameters are ${\varphi }_0=0.3,{\varphi }_{\infty }=0.8,\chi =2.5,\kappa =1$, $\rho =0.3$, and $\eta =1.0$. Both patterns are the closeups of a larger system size simulation as presented by the bottom panels of Fig. 2. See also SM1.mp4 in the Supplemental Material [47].