Phase Field Modeling of Dictyostelium Discoideum Chemotaxis

TL;DR

This paper introduces a phase field model for Dictyostelium discoideum chemotaxis, integrating chemical signaling and cell movement to replicate pseudopodia patterns observed during chemotaxis.

Contribution

It develops a novel coupled phase field and reaction-diffusion model that captures cell movement and chemical dynamics simultaneously.

Findings

Successfully mimics pseudopodia patterns during chemotaxis

Efficient spectral method for numerical solution

Extends phase-field formulation for membrane-restricted components

Abstract

A phase field approach is proposed to model the chemotaxis of Dictyostelium discoideum. In this framework, motion is controlled by active forces as determined by the Meinhardt model of chemical dynamics which is used to simulate directional sensing during chemotaxis. Then, the movement of the cell is achieved by the phase field dynamics, while the reaction-diffusion equations of the Meinhardt model are solved on an evolving cell boundary. This task requires the extension of the usual phase-field formulation to allow for components that are restricted to the membrane. The coupled system is numerically solved by an efficient spectral method under periodic boundary conditions. Numerical experiments show that our model system can successfully mimic the typically observed pseudopodia patterns during chemotaxis.

Tables

Table 1. Table 2.1: Parameters in the Meinhardt model coupled with a phase field membrane.

Parameter	Value	Parameter	Value
$D_a$	8e-3	$r_a$	0.2
$D_c$	1.8$D_a$	$r_b$	0.3
${\nabla }_{⟂}$	1.0	$r_c$	0.13
$s_a$	5e-4	$b_a$	0.1
${\tau }_0$	0.01	$b_c$	0.05
${\tau }_0$	0.01	$s_c$	0.2
$dr$	0.02	$C_{\text{chem}}$	0.02

Equations

$$\frac{\partial a}{\partial t}$$

$$\frac{\partial b}{\partial t}$$

$$\frac{\partial c}{\partial t}$$

$$s(\text{r},t)=(1+d_r\xi )(1+C_{\text{chem}}f(\text{r},t)),$$

$$f(\text{r},t)=1-\frac{\text{dist}(\text{r},t)-{\text{dist}}_{\min }(t)}{{\text{dist}}_{\max }(t)-{\text{dist}}_{\min }(t)}$$

$${\tau }_0\frac{\partial (g(\varphi )a)}{\partial t}+\nabla \cdot (g(\varphi )a\mathbf{v})=D_a{\nabla }_{\parallel }\cdot (g(\varphi ){\nabla }_{\parallel }a)+D_{\perp }{\nabla }_{\perp }\cdot (g(\varphi ){\nabla }_{\perp }a)$$

$$+g(\varphi )\left(\frac{s(\text{r},t)(a^2{b}^{-1}+b_a)}{(s_c+c)(1+s_a{a}^2)}-r_{a}a\right),$$

$${\tau }_0\frac{\partial b}{\partial t}=r_b\frac{\int g(\varphi )a\text{d}\text{r}}{\int g(\varphi )\text{d}\text{r}}-r_{b}b,$$

$${\tau }_0\frac{\partial (g(\varphi )c)}{\partial t}+\nabla \cdot (g(\varphi )c\mathbf{v})=D_c{\nabla }_{\parallel }\cdot (g(\varphi ){\nabla }_{\parallel }c)+D_{\perp }{\nabla }_{\perp }\cdot (g(\varphi ){\nabla }_{\perp }c)+g(\varphi )\left(b_{c}a-r_{c}c\right).$$

$${\nabla }_{\parallel }=\left[\begin{array}{cc}[c]n_y^2& -n_x{n}_y\\ -n_x{n}_y& n_x^2\end{array}\right]\left[\begin{array}{c}[c]{\partial }_x\\ {\partial }_y\end{array}\right],{\nabla }_{\perp }=\left[\begin{array}{cc}[c]n_x^2& n_x{n}_y\\ n_x{n}_y& n_y^2\end{array}\right]\left[\begin{array}{c}[c]{\partial }_x\\ {\partial }_y\end{array}\right],$$

$$E_{\text{ten}}=\gamma {\int }_{\Omega }\left(\frac{ϵ}{2}|\nabla \varphi {|}^2+\frac{1}{ϵ}G(\varphi )\right)\text{d}\mathbf{r},$$

$${\mathbf{F}}_{\text{ten}}=\frac{\delta E_{\text{ten}}}{\delta \varphi }\frac{\nabla \varphi }{ϵ|\nabla \varphi {|}^2}=\frac{\gamma }{ϵ}\left(-ϵ{\nabla }^2\varphi +\frac{1}{ϵ}{G}^{\prime }(\varphi )\right)\frac{\nabla \varphi }{|\nabla \varphi {|}^2}.$$

$$E_{\text{bend}}=\frac{\kappa }{2}{\int }_{\Omega }\frac{1}{ϵ}{\left(ϵ{\nabla }^2\varphi -\frac{1}{ϵ}{G}^{\prime }(\varphi )\right)}^2\text{d}\mathbf{r},$$

$${\mathbf{F}}_{\text{bend}}=\frac{\delta E_{\text{bend}}}{\delta \varphi }\frac{\nabla \varphi }{ϵ|\nabla \varphi {|}^2}=\frac{\kappa }{{ϵ}^2}\left(ϵ{\nabla }^2-\frac{1}{ϵ}{G}^{\prime \prime }\right)\left(ϵ{\nabla }^2\varphi -\frac{1}{ϵ}{G}^{\prime }\right)\frac{\nabla \varphi }{|\nabla \varphi {|}^2}.$$

$${\mathbf{F}}_{\text{area}}=M_{\text{area}}\left({\int }_{\Omega }\varphi \text{d}\mathbf{r}-A_0\right)\frac{\nabla \varphi }{|\nabla \varphi |},$$

$${\mathbf{F}}_{\text{chem}}=-\alpha \tilde{a}\frac{\nabla \varphi }{|\nabla \varphi |},$$

$${\mathbf{F}}_{\text{tot}}={\mathbf{F}}_{\text{ten}}+{\mathbf{F}}_{\text{bend}}+{\mathbf{F}}_{\text{area}}+{\mathbf{F}}_{\text{chem}}+{\mathbf{F}}_{\text{fr}}=0$$

$$\mathbf{v}=-\frac{1}{\tau }{\mathbf{F}}_{\text{fr}}=\frac{1}{\tau }({\mathbf{F}}_{\text{ten}}+{\mathbf{F}}_{\text{bend}}+{\mathbf{F}}_{\text{area}}+{\mathbf{F}}_{\text{chem}}).$$

$$\tau \frac{\partial \varphi }{\partial t}=$$

$$-M_{\text{area}}\left(\int \frac{ϵ}{2}|\nabla \varphi {|}^2+\frac{1}{ϵ}G(\varphi )\text{d}\text{r}-P_0\right)|\nabla \varphi |+\alpha \tilde{a}|\nabla \varphi |.$$

$${\varphi }^0(\mathbf{x})=0.5+0.5\tanh \left(\frac{r-\text{dist}(\mathbf{x},0)}{ϵ/3}\right),$$

$$a^0(\mathbf{x})\equiv 0,b^0(\mathbf{x})\equiv 0.01,c^0(\mathbf{x})\equiv 0.$$

$$\text{CI}=\frac{\sum_{n}\Big{\langle}\mathbf{x}_{\text{center}}^{n+1}-\mathbf{x}_{\text{center}}^{n},\frac{\mathbf{r}_{0}-\mathbf{x}_{\text{center}}^{n}}{\|\mathbf{r}_{0}-\mathbf{x}_{\text{center}}^{n}\|}\Big{\rangle}}{\sum_{n}\Big{\langle}\mathbf{x}_{\text{center}}^{n+1}-\mathbf{x}_{\text{center}}^{n},\mathbf{x}_{\text{center}}^{n+1}-\mathbf{x}_{\text{center}}^{n}\Big{\rangle}}.$$

$$\tau \frac{{\varphi }^{n+1}-{\varphi }^n}{\Delta t}=$$

$$-M_{\text{area}}\left(\int \frac{ϵ}{2}|\nabla {\varphi }^n{|}^2+\frac{1}{ϵ}G({\varphi }^n)\text{d}\text{r}-P_0\right)|\nabla {\varphi }^n|+\alpha \tilde{a}|\nabla {\varphi }^n|.$$

$$\left(\frac{\tau }{\Delta t}+\kappa {\nabla }^4-\gamma {\nabla }^2\right){\varphi }^{n+1}=\text{RHS}({\varphi }^n),$$

$${\tau }_0\frac{\partial (\tilde{g}(\varphi )a)}{\partial t}+\nabla \cdot (\tilde{g}(\varphi )a\mathbf{v})=D_a{\nabla }_{\parallel }\cdot (g(\varphi ){\nabla }_{\parallel }a)+D_{\perp }{\nabla }_{\perp }\cdot (g(\varphi ){\nabla }_{\perp }a)$$

$$+\tilde{g}(\varphi )\left(\frac{s(\text{r},t)(a^2{b}^{-1}+b_a)}{(s_c+c)(1+s_a{a}^2)}-r_{a}a\right),$$

$${\tau }_0\frac{\partial b}{\partial t}=r_b\frac{\int \tilde{g}(\varphi )a\text{d}\text{r}}{\int \tilde{g}(\varphi )\text{d}\text{r}}-r_{b}b,$$

$${\tau }_0\frac{\partial (\tilde{g}(\varphi )c)}{\partial t}+\nabla \cdot (\tilde{g}(\varphi )c\mathbf{v})=D_c{\nabla }_{\parallel }\cdot (g(\varphi ){\nabla }_{\parallel }c)+D_{\perp }{\nabla }_{\perp }\cdot (g(\varphi ){\nabla }_{\perp }c)+\tilde{g}(\varphi )\left(b_{c}a-r_{c}c\right),$$

$${\nabla }_{\parallel }\cdot (g(\varphi ){\nabla }_{\parallel }a)$$

$$=\frac{ϵ}{2}\left[\begin{array}{cc}[c]n_y^2& -n_x{n}_y\\ -n_x{n}_y& n_x^2\end{array}\right]\left[\begin{array}{c}[c]{\partial }_x\\ {\partial }_y\end{array}\right]\cdot \left[\begin{array}{c}[c]({\partial }_y\varphi {)}^2{\partial }_{x}a-({\partial }_x\varphi {\partial }_y\varphi ){\partial }_{y}a\\ -({\partial }_x\varphi {\partial }_y\varphi ){\partial }_{x}a+({\partial }_x\varphi {)}^2{\partial }_{y}a\end{array}\right]$$

$$\displaystyle=\frac{\epsilon}{2}\Bigg{[}\Big{(}n_{y}^{2}\partial_{x}-n_{x}n_{y}\partial_{y}\Big{)}\Big{(}(\partial_{y}\phi)^{2}\partial_{x}a-(\partial_{x}\phi\partial_{y}\phi)\partial_{y}a\Big{)}$$

$$\displaystyle\qquad\qquad\qquad+\Big{(}-n_{x}n_{y}\partial_{x}+n_{x}^{2}\partial_{y}\Big{)}\Big{(}-(\partial_{x}\phi\partial_{y}\phi)\partial_{x}a+(\partial_{x}\phi)^{2}\partial_{y}a\Big{)}\Bigg{]},$$

$${\nabla }_{\perp }\cdot (g(\varphi ){\nabla }_{\perp }a)$$

$$=\frac{ϵ}{2}\left[\begin{array}{cc}[c]n_x^2& n_x{n}_y\\ n_x{n}_y& n_y^2\end{array}\right]\left[\begin{array}{c}[c]{\partial }_x\\ {\partial }_y\end{array}\right]\cdot \left[\begin{array}{c}[c]({\partial }_x\varphi {)}^2{\partial }_{x}a+({\partial }_x\varphi {\partial }_y\varphi ){\partial }_{y}a\\ ({\partial }_x\varphi {\partial }_y\varphi ){\partial }_{x}a+({\partial }_y\varphi {)}^2{\partial }_{y}a\end{array}\right]$$

$$\displaystyle=\frac{\epsilon}{2}\Bigg{[}\Big{(}n_{x}^{2}\partial_{x}+n_{x}n_{y}\partial_{y}\Big{)}\Big{(}(\partial_{x}\phi)^{2}\partial_{x}a+(\partial_{x}\phi\partial_{y}\phi)\partial_{y}a\Big{)}$$

$$\displaystyle\qquad\qquad\qquad+\Big{(}n_{x}n_{y}\partial_{x}+n_{y}^{2}\partial_{y}\Big{)}\Big{(}(\partial_{x}\phi\partial_{y}\phi)\partial_{x}a+(\partial_{y}\phi)^{2}\partial_{y}a\Big{)}\Bigg{]}.$$

Full text

††journal: Journal of XXX

1 Introduction

Many cells have an internal “compass”, which enables them to navigate through various environments. This “compass” detects the gradients in chemical concentrations, the rigidity of extra-cellular matrix, cellular adhesion sites, fluidic shear stress, etc. Interestingly, extensive studies on how such a “compass” is realized have led to “taxis-mania, focusing on chemotaxis, durotaxis, mechanotaxis, haptotaxis, plithotaxis and so on.

In this work, we specifically concentrate on chemotaxis, which plays an extensive role in many physiological processes [1]. For example, primordial cells are capable of figuring out their way to proper locations by sensing chemical clues, thus correctly forming the organs. Chemotaxis is also essential for immune responses and wound healing. In addition to these normal physiological processes, the pathology of numerous diseases, such as cancer metastasis and inflammatory disorders, is believed to be related to chemotaxis [3, 2].

Dictyostelium discoideum, a type of amoeboid cell, is a popular model system for the study of chemotaxis. These cells rely on chemotaxis to find nutrition. When suffering from starvation, they are capable of chemotaxing in response to cAMP gradients in order to aggregate and enhance their chances for survival. Dicty chemotactic behavior is very similar to that of human leukocytes. Conceptually, chemotaxis can be divided into motility, directional sensing, and polarity [4, 5, 6]. According to experiments [6, 7, 8], motility involves periodic extensions and retractions of pseudopods – temporary actin-filled protrusions of the cell membrane. Directional sensing refers to the process by which cells sense the chemical gradients and adjust their direction. Polarity refers tp the reorganization of the cell interior to favor moving in a fixed direction. Once polarized, protrusions mainly extend from the cell anterior, regardless of whether a chemical gradient exists or not.

1.1 Directional sensing: LEGI-BEN and Meinhardt models

Detailed experimental investigations have uncovered many important features of chemotaxis [9]. First of all, the actin cytoskeleton in motile cells exhibits many characteristics of an excitable medium, such as the presence of propagating waves in the cell membrane [10, 11, 12, 13].This fact has suggested that an excitable network, which composes a simple activator-inhibitor system, may help explain the spontaneous migration of these cells. Then, to include the cell’s response to external gradients, one can modify the activator-inhibitor system by adding a steering bias: higher concentrations of chemoattractants will lower the threshold of excitability, thus causing more excitation. Over time, cells with this bias will tend to move along the directions towards higher chemoattractant density. Chemotactic cells can display surprising sensitivity, responding to a chemical gradient as small as 1%. Such sensitivity can be captured by a biased excitable network [14].

Although successful in qualitatively explaining spontaneous cell motion in the absence of chemoattractants as well as directed motion in the presence of chemoattractants, the biased excitable network (BEN) approach still misses some features of realistic chemotactic behavior. One such feature is the adaptability of cells to external chemical clues. Since higher concentrations of chemoattractants lower the threshold for excitable behavior according to a biased excitable network, the cells would be predicted to become ”hyper-excited” if they are exposed to a spatially uniform increase in chemoattractant concentration In fact, this does not occur and after transient responses to uniform increments in chemoattractants, the excitability returns to baseline levels. To account for this behavior, a local-excitation, global-inhibition (LEGI) mechanism was proposed [15, 16]. According to this mechanism, chemoattractants give rise to the release of a slowly diffusing activator which is accompanied by a rapidly diffusing inhibitor. Thanks to the regulation by the global inhibitor, the level of excitation returns to a threshold level.

Thus, one natural strategy is to directly combine the LEGI mechanism and the biased excitable network into a hybrid LEGI-BEN model [17]. One can however use a different model which already contains the two essential components of excitation and adaptation. More than a decade ago, Meinhardt proposed a three-component model [18], which shares a similar conceptualization with the LEGI-BEN approach. Here, a biased excitable network, which includes a slowly diffusing activator and a fast inhibitor, is further regulated by an extra even faster global inhibitor. The total quantity of the activator is held approximately constant over time by this global inhibitor, thus preventing a large increase in excitable behavior. For the sake of mathematical simplicity, we choose to use the Meinhardt model in our modeling efforts, whereas the LEGI-BEN model can also in the future be embedded in our phase field model, if needed.

The governing equations of Meinhardt model are:

[TABLE]

Here, $a,b$ and $c$ are respectively the local activator, the global inhibitor, and the local inhibitor. By referring to $c$ as a global inhibitor, we mean $D_{b}\gg D_{c}$, so that the inhibiting effects of $b$ can spread over the whole interface in almost no time, thus regulating the total quantity of the activator. Therefore we can assume the concentration of $b$ to be uniform everywhere on the membrane, which leads to the replacement of partial differential equation by an ordinary differential equation with a nonlocal source. This system exhibits nice bifurcating patterns [18], which have been successfully compared in experimental findings [19]. A local excitation bifurcates into a pair of competing daughter bursts of excitations, which travel in opposite directions. One of these daughter bursts will win out over the other, which vanishes. A new bifurcation will occur once the “loser” dies. The decisive factor for the competition between the pair of bifurcations is the spatial factor $s(\textbf{r},t)$ in the excitation of the activator:

[TABLE]

where $f(\textbf{r},t)$ is a function of the spatial clue related to the concentrations of chemoattractant. The factor $(1+d_{r}\xi)$ represents the effect of stochastic fluctuations, with $\xi$ taken as white noise and $d_{r}$ being the fluctuation strength. In this work, we assume the existence of a chemoattractant source $\textbf{r}_{0}$, and simply set $f(\textbf{r},t)$ as:

[TABLE]

where $\text{dist}(\textbf{r},t)$ represents the distance between any point r and the chemoattractant source $\textbf{r}_{0}$, while $\text{dis}_{\max}(t)$ and $\text{dis}_{\min}(t)$ respectively represent the maximal and minimal distance at time $t$, between the chemoattractant source $\textbf{r}_{0}$ and all points on the cell membrane. It is evident that $f(\textbf{r},t)$ varies between [math] and $1$, monotonically decreasing with the distance from the chemoattractant source $\textbf{r}_{0}$, along the cell membrane. The parameter $C_{\text{chem}}$ regulates the strength of the bias on the threshold of the activator’s excitability. It turns out that the values of $C_{\text{chem}}$ as small as 0.01, can still significantly affect the bifurcation patterns in space and time. This fact is consistent with the chemotactic cells’ sensitivity to chemical gradients as mentioned above. Statistically, excitation bursts propagating toward favorable positions in the chemical gradients are more likely to survive and continue to bifurcate, thus fostering the directed navigation of chemotactic cells.

1.2 Phase field model framework

In the past few decades, phase field models have emerged as one of the most successful methods for studying interfacial problems; see the two review articles [22, 23] and the references therein. In the phase field model framework, a phase field function $\phi$ is introduced, assigning a value (say, 0) for one phase, and another value (say, 1) for the other. In the interfacial region, the phase field $\phi$ rapidly but smoothly transitions from 0 to 1. The interface is tracked by the 1/2-level set during the morphological evolution. The main advantage of the phase field approach is that it can allow for the computation of the temporal evolution of arbitrary morphologies and complex microstructures without explicitly tracking interfaces.

In the realm of biology, cell shape dynamics and cell migration processes have been simulated by using phase field models. In [25], a quantitative model for cell shape and motility dynamics was constructed based on the original phase field concept [26]. An auxiliary field is introduced to distinguish the cell’s interior ($\phi=0$) from the exterior ($\phi=1$). The dynamics of the cell are governed by equations that couple this field to the actual physical degrees of freedom, and the diffuse layer separating the interior from the exterior marks the membrane location.

This cell motility model is part of a larger set of recent theoretical studies that have attempted to model cell migration. For example, some studies have attempted to calculate the ”flow” of the actin cytoskeleton in a one-dimensional [27, 28, 29] or fixed two-dimensional cell geometry [30]. In some works, the cell boundary was allowed to change according to a phenomenological function of protrusion rate [31, 32] while other approaches implemented physical forces along the cell membrane, obtained cell shape and speed, but ignored actin flow and detailed adhesion mechanisms [25, 33]. Yet others examined adhesion dynamics and cell-substrate coupling while ignoring cell deformations [34] or focused on the dynamics of the leading edge [35]. Ziebert, Aranson and their coworkers studied the cell shape dynamics by coupling a vector field model of the actin filament network with the cell shape [36]. Finally, a more comprehensive model for cell migration was presented in [37] which couples actin flow with discrete adhesion sites and deformable cell boundaries. Other interesting patterns such as periodic migration [38] or circular motion [39, 40] have also been studied using phase field framework.

2 Phase field model of Chemotaxis

2.1 Pseudopodia: Chemical dynamics on a phase field membrane

Our goal in this paper is to couple the aforementioned directional sensing system to a computational model of the resultant motion. With the intricate spatio-temporal patterns of the Meinhardt model in hand, our immediate challenge is how to make this happen on the membrane of a dynamically evolving cell described by a phase field. To our best knowledge, there are very few published attempts to study similar problems. For example, Nelson et al. applied a hybrid computational framework to couple the Meinhardt model with cell movement [19]. There, the movement of the cell is achieved using a level set method [24], while the reaction-diffusion equations of the Meinhardt model are approximated on an evolving cell boundary using an arbitrary Lagrangian-Eulerian surface finite element method (ALE-SFEM). In our phase field model, we can achieve the same effect in a much simpler manner. To accomplish this, we modified the approach used to couple bulk reaction-diffusion systems to phase field cells [25, 39]. Specifically, instead of a factor of $\phi$ to limit reaction to the interior, we use $g(\phi(\textbf{r}))=\frac{\epsilon}{2}|\nabla\phi|^{2}$ to restrict concentrations to the membrane (see below). In addition, we found it necessary to insert diffusion in the normal direction of the phase field interface, so that reaction-diffusion processes in different layers synchronize with each other. Our revised equations for the membrane Meinhardt system are:

[TABLE]

where $\mathbf{v}$ is the interface velocity equal to $-\partial_{t}\phi\frac{\mathbf{\nabla}\phi}{|\mathbf{\nabla}\phi|^{2}}$. As already mentioned, the global inhibitor $b$ immediately spreads over the whole membrane, and it satisfies an ordinary differential equation instead of a partial differential equation. The term $g(\phi(\textbf{r}))=\frac{\epsilon}{2}|\nabla\phi|^{2}$ is only nonzero in the interfacial region so that the reaction-diffusion dynamics only occur near the interface, and $D_{\perp}$ refers to the diffusion we add in the normal direction of the membrane. More specifically, $\nabla_{\parallel}$ and $\nabla_{\perp}$ read:

[TABLE]

in which the normal vector $\textbf{n}=[n_{x},n_{y}]^{T}$ can be calculated by $\textbf{n}=-\frac{\nabla\phi}{|\nabla\phi|}$. In practice, $\nabla_{\perp}$ has to be much larger than the other diffusion coefficients in the chemical systems. The detailed parameters used in our model simulations are listed in Table 2.1.

2.2 Chemotaxis dynamics of Dictyostelium discoideum

We model the Dictyostelium discoideum cell as a 2d region with a fixed area $A_{0}$. The evolving shape of the cell membrane is determined by the competition of several forces, including surface tension, bending force, the pressure that constrains the cell area, the chemical protrusive force which is proportional to the density of local activator $a$, and the effective friction due to the interaction between cell membrane and the substrate. All of them are formulated under the phase field framework, as follows.

Given the surface energy in phase field formulation [41]:

[TABLE]

in which $\gamma$ is the surface tension, $\epsilon$ is the phase field parameter controlling the width of the cell membrane (the width of phase field interface), and $G(\phi)=18\phi^{2}(1-\phi)^{2}$ is a double well potential with minima at $\phi=0$ and $\phi=1$, the surface tension is derived by taking the variational derivative of the surface energy [25],

[TABLE]

Similarly, given the bending energy in phase field formulation [41]:

[TABLE]

with $\kappa$ being the bending rigidity, we obtain the bending force [25]:

[TABLE]

The area force is given as a soft penalty on cell area:

[TABLE]

in which $M_{\text{area}}$ is the penalty constant. We assume the protrusion force is simply proportional to the density of the local activator $a$ on the membrane. Given the fact that the activator’s concentration may change over several magnitudes, we added a saturating restriction to the force $|\mathbf{F}_{\text{chem}}|\propto\tilde{\alpha}=\max(10,a)$,

[TABLE]

in which $\alpha$ is the strength of the chemical protrusion force. One could alternatively use a sigmoidal function to achieve the same effect. A friction force due to the interaction between the cell and the substrate (such as adhesion, attachment and detachment of the cell from the substrate) is introduced which is proportional to the local speed: $\mathbf{F}_{\text{fr}}=-\tau\mathbf{v}$. The force balance at quasi-steady state

[TABLE]

yields

[TABLE]

Finally using the transport equation of the phase field $\phi$ along the velocity field $\mathbf{v}$: $\frac{\partial\phi}{\partial t}+\mathbf{v}\cdot\nabla\phi=0$, we obtain the following equation for $\phi$:

[TABLE]

Physically, this chemotaxis dynamics of $\phi$ implies that the friction force on the cell is balanced by the chemical protrusion force, which is transmitted from the substrate onto the cell via adhesion complexes.

3 Numerical Simulations

In this section, we present several numerical experiments of cell chemotaxis dynamics. The numerical method described in the Appendix is adopted to solve the chemotaxis dynamics (2.2) coupled with the Meinhardt dynamics (2.1)-(2.3). In all of the simulations, we take $L_{x}=L_{y}=10,N_{x}=N_{y}=2^{8},N_{t}=60$ and $\Delta t=5e-4$. The interfacial width of $\phi$ is fixed as $\epsilon=10h_{x}$, where $h_{x}=2L_{x}/N_{x}$ is the grid spacing in the $x$ direction.

For initial data, we take $\phi^{0}$ as a disk with center at origin and radius $r=4$:

[TABLE]

in which $\text{dist}(\mathbf{x},\mathbf{0})$ stands for the Euclidean distance between $\mathbf{x}$ and the origin. We further take

[TABLE]

Unless otherwise specified, the chemoattractant source is located at $\mathbf{r}_{0}=(0,-40)^{T}$ and the strength of the bias is $C_{\text{chem}}=0.02$.

3.1 Meinhardt dynamics on the membrane of a fixed cell

In this example, we test the Meinhardt dynamics on the membrane of a fixed phase field cell as given in (3.1). The numerical simulation is presented in Figure 3.1 from which the bifurcating patterns are clearly observed on the cell membrane. In this simulation, several time snapshots are taken at $t=10.0,11.0,13.0,14.0,16.8,18.0$. A bifurcation occurs at $t=10.0$. After a short time period, two branches are formed at $t=11.0$. One branch defeats the other at $t=13.0$. Then the bifurcation repeatedly recurs. In each subfigure, the black dashed circle represents the cell membrane. The yellow dashed lines indicate the small box in which the Meinhardt equations are solved (see Appendix for the details about the small box).

3.2 Meinhardt dynamics on the membrane of a free cell

We now allow the cell to deform and move by set the cell solving the phase field equation (2.2) together with the Meinhardt system (2.1)-(2.3). A set of numerical results is presented in Figure 3.2. In this figure, the top row is the four snapshots at $t=20,30,40,50$, on which the yellow color indicates a high concentration of activator $a$. The bottom row shows the cell trajectories for times up to $t=20,30,40,50$. Our simulations turn out to be generally consistent with observed chemotactic behavior. Pseudopods are randomly generated in the cell membrane, while the biasing effects towards the direction of the chemoattractant source $\textbf{r}_{0}$ accumulate over time, which eventually leads to the cell’s translation along favorable directions. The trajectory of our simulated cell has demonstrated high efficiency in the cell’s navigation.

Moreover, even when the chemoattractant source $\mathbf{r}_{0}$ is suddenly moved, a strong adaptability of our simulated cell is also observed in Figure 3.3. In this simulation, the chemoattractant source $\mathbf{r}_{0}$ is located at $\mathbf{r}_{0}=(0,-40)^{T}$ over time $[0,40]$. During this time period, the cell moves towards $\mathbf{r}_{0}$ similarly to that shown in Figure 3.2. At $t=40$, $\mathbf{r}_{0}$ is suddenly changed to $\mathbf{r}_{0}=(40,0)^{T}$. The cell can quickly adjust the direction and move towards the new location over time $[40,80]$, still with a slight biasing effect. At $t=80$, a new location $\mathbf{r}_{0}=(0,40)^{T}$ is assigned, and the cell changes direction and moves towards the north. For the sake of clear display, the cells are plotted in a window of $[-15,25]\times[-25,15]$, without showing the locations of chemoattractant source.

3.3 Chemotaxis index

To quantify the efficiency of cells’ navigation, we measured the chemotaxis index (CI) of our simulated cells over time. The chemotaxis index is defined by the ratio between the distance the simulated cell has traveled in the direction toward the chemoattractant source $\mathbf{r}_{0}$ and the total distance it has traveled:

[TABLE]

In Figure 3.4, we plot the CI for three simulated cells with various $C_{\text{chem}}=0.01,0.03,0.05$ of the strength of the bias. The top row, from left to right, are the trajectories of three cells up to time $T=40$, with different values of $C_{\text{chem}}=0.01,0.03,0.05$, respectively. The cells are plotted in the box $[-20,20]\times[-30,10]$, with the chemoattrant source located at $\text{r}_{0}=[0,-40]^{T}$. The bottom row, from left to right, are the CI of the three cells up to time $T=40$, with different values of $C_{\text{chem}}=0.01,0.03,0.05$, respectively. In each plot of the CI, CI curve is in blue, with the values indicated by the left $y$-axis. The solid and dashed orange curves (values indicated by the right $y$-axis) are the cumulative distance travelled by the cell in the direction toward $\textbf{r}_{0}$, and the cumulative total distance, respectively.

Our result shows that when the strength of the bias is weaker ($C_{\text{chem}}=0.01$), the cell wanders more randomly along the trajectory toward the chemoattractant source; while when the strength of the bias becomes stronger ($C_{\text{chem}}=0.05$), the cell’s moving direction is more straightforward. On the other hand, even for the case with a weaker strength of the bias $C_{\text{chem}}=0.01$, our simulation still shows a high value for CI (CI becomes close to $0.6$ at $t=40$), indicating efficient cell navigation. For the case with a stronger strength of the bias $C_{\text{chem}}=0.05$, the CI can reach an even higher value, close to $0.9$.

4 Discussion and outlook

Here we have presented a new computational approach to the problem of chemotactically driven motion in Dictyostelium. Our results are consistent with experimental findings and recapitulate results obtained by the level set method [19]. While both the level set method and phase field model can couple the Meinhardt dynamics with membrane evolution, thereby successfully reproducing the pseudopod morphology, our method enjoys the fact that the coupling of the phase field with the Meinhardt dynamics is much more straightforward, by simply incorporating $g(\phi)=\frac{\epsilon}{2}|\nabla\phi|^{2}$ in the Meinhardt equations (2.1-2.3). In addition, although it has been claimed that an evolving cell boundary solved by (ALE-SFEM) plus the level set approach has the advantage of efficiency [19, 20], such a computational strategy suffers from a very complicated implementation. Indeed, the level-set-modeled cell profile has to “communicate” with the evolving cell boundary at every single time step. More specifically, two set of meshes need to be introduced, the finite element mesh for the cell membrane update, and level set mesh for the update of the level set function. In each time step, one needs to project the finite element mesh points onto level set mesh points by using nearest-neighbor point, in order to update level set function; then use the level set mesh points to form a new finite element mesh, on which the Meinhardt is updated. Since the new finite element mesh may fail to be equidistributed, a step of re-gridding the finite element mesh is often. In contrast, in our phase field framework, when the cell movement and membrane reaction-diffusion system are combined using the same implicit tracking language, no further “communication” between mechanical and chemical systems is needed, as they are solved on the same uniform mesh (see numerical details in Appendix).

Finally, in addition to computational efficiency and formulaic complexity, our phase field framework has an obvious advantage when it becomes necessary to couple intra-cellular flow and focal adhesions into the model [37]. In fact, it is definitely worth exploring chemotaxis with a more biophysically complete model, i.e. with all the effects of bulk reaction-diffusion dynamics, cytoplasmic flow, Meinhardt patterns, focal adhesions, and membrane forces. Quite a few other interesting features, apart from bifurcating pseudopods, were observed in experiments more than a decade ago [21]. For example, when Dictyostelium discoideum cells adhere to the substrate, they exert opposing pole forces that are orders of magnitude higher than required to overcome the resistance from their environment. Also, the strain energy exerted by migrating Dicty on the substrate is (almost) quasi-periodic and can be used to identify different stages of the cell motility cycle. Moreover, the period displays an inversely proportional relation with cell velocity. In recent work by Copos et al. [44], a simple mechanochemical model of 2D (in the vertical plane) cell motility was used to study the periodic changes in cell length and the related spatiotemporal dynamics of traction forces. Our phase field model can provide a platform to carefully study these features in future work.

5 Acknowledgements

H. Levine’s work is supported by the NSF, grants Nos. PHY-1935762 and PHY-2019745. Y. Zhao’s work is supported by a grant from the Simons Foundation through Grant No. 357963 and NSF grant DMS2142500.

Appendix

In this appendix, we present in detail the numerical algorithm to solve the coupled system (2.1)-(2.3) and (2.2). We take the computational domain $\Omega=[-L_{x},L_{x})\times[-L_{y},L_{y})$. Periodic boundary conditions are used for the coupled system. A uniform grid $\Omega_{h}$ is generated over $\Omega$ by taking $h_{x}=\frac{2L_{x}}{N_{x}}$ and $h_{y}=\frac{2L_{y}}{N_{y}}$. The grid points are given as $(x_{i},y_{j})=(-L_{x}+(i-1)h_{x},-L_{y}+(j-1)h_{y})$. Given initial data $(\phi^{0},a^{0},b^{0},c^{0})$, we aim to find $(\phi^{n},a^{n},b^{n},c^{n})$ for $n=1,2,\cdots,N_{t}$ with $N_{t}=\frac{T}{\Delta t}$.

For solving the equation of $\phi$ in (2.2), we adopt a semi-implicit Fourier spectral method. More specifically, we discretize the equation as

[TABLE]

This discretization can be rewritten as

[TABLE]

which can be efficiently solved by the Fourier spectral method.

Next, we consider the numerical method for the Meinhardt system (2.1)-(2.3). For the sake of numerical stability, the Meinhardt equations (2.1)-(2.3) are replaced by:

[TABLE]

in which $\tilde{g}(\phi)=\frac{1}{\epsilon}G(\phi)$. The replacement of $g$ by $\tilde{g}$ is reasonable due to the fact that in the Ginzburg-Landau functional (2.5), the term $g(\phi)=\frac{\epsilon}{2}|\nabla\phi|^{2}$ plays identical role as $\tilde{g}(\phi)=\frac{1}{\epsilon}G(\phi)$ at the system equilibrium [42, 43].

Note that we replace all the terms of $g(\phi)$ by $\tilde{g}(\phi)$ except for those in the parallel and perpendicular diffusion terms. We do this is because by taking $g(\phi)=\frac{\epsilon}{2}|\nabla\phi|^{2}$ together with the parallel and perpendicular gradient operators (2.4), the diffusion terms can be significantly simplified. Explicitly, for the parallel and perpendicular diffusion terms in the equation of $a$,

[TABLE]

and

[TABLE]

To numerically discretize the above two terms, firstly we evaluate $(\partial_{x}\phi,\partial_{y}\phi)$ using Fourier spectral method, and calculate $(n_{x},n_{y})$ as

[TABLE]

in which $\epsilon_{0}$ is a sufficiently small constant (say, $\epsilon_{0}=1e-8$) to avoid dividing zero. Secondly, we evaluate

[TABLE]

using Fourier spectral method. Thirdly, the first and second derivatives of $a$ are evaluated by central difference:

[TABLE]

Inserting all evaluations above back, we obtain the numerical approximation of $\nabla_{\parallel}\cdot(g(\phi)\nabla_{\parallel}a)$ and $\nabla_{\perp}\cdot(g(\phi)\nabla_{\perp}a)$.

The advection term $\nabla\cdot(\tilde{g}(\phi)a\mathbf{v})$ in the equation of $a$ is approximated by a central difference scheme:

[TABLE]

in which $\mathbf{v}=[v^{x},v^{y}]^{T}=-\partial_{t}\phi\frac{\nabla\phi}{|\nabla\phi|^{2}}$ and is calculated by taking $\partial_{t}\phi\approx\frac{\phi^{n+1}-\phi^{n}}{\Delta t}$, and $\nabla\phi$ evaluated by Fourier spectral approximation.

The time derivative $\frac{\partial(\tilde{g}(\phi)a)}{\partial t}$ is approximated by forward Euler scheme,

[TABLE]

Finally, with all terms discretized in the equation of $a$, we get an update on $a$: $a^{n}\rightarrow a^{n+1}$. The equation (A.3) can be solved numerically in a similar manner. The equation (A.2) is an ODE, so we can adopt an efficient fourth-order Runge-Kutta method (RK4) to solve it.

Since the phase field cell $\phi$ moves around in the computational domain $\Omega$ and may near the edge, we do not solve the phase field equation and the Meinhardt equations in the entire domain. We only solve these equations in a smaller box of size $1.75L_{x}\times 1.75L_{y}$ near the cell. This box is re-centered if the cell is close to one of its four boundaries: if $\phi\geq 0.5$ within $\frac{N_{x}}{16}$ (or $\frac{N_{y}}{16}$) pixels of the boundary, the box is shifted $\frac{N_{x}}{4}$ (or $\frac{N_{y}}{4}$) pixels away from the boundary. We treat the small box as having periodic boundary conditions, which is appropriate as we keep the cell from too closely approaching the edge of $\Omega$.