# Phase field model for cell spreading dynamics

Mohammad Abu Hamed  Department of Mathematics, Technion - Israel Institute of
Technology, Haifa 32000, Israel  Department of Mathematics, The College of
Sakhnin - Academic College for Teacher Education, Sakhnin 30810, Israel
Alexander A. Nepomnyashchy  Department of Mathematics, Technion - Israel
Institute of Technology, Haifa 32000, Israel

###### Abstract

We suggest a 3D phase field model to describe 3D cell spreading on a flat
substrate. The model is a simplified version of a minimal model that was
developed in Winkler  _et al._ (2019). Our model couples the order parameter
$u$ with 3D polarization (orientation) vector field P of the actin network. We
derive a closed integro-differential equation governing the 3D cell spreading
dynamics on a flat substrate, which includes the normal velocity of the
membrane, curvature, volume relaxation rate, a function determined by the
molecular effects of the subcell level, and the adhesion effect. This equation
is easily solved numerically. The results are in agreement with the early fast
phase observed experimentally in Dobereiner  _et al._ (2004). Also we find
agreement with the universal power law Cuvelier  _et al._ (2007) which suggest
that cell adhesion or contact area versus time behave as $\sim t^{1/2}$ in the
early stage of cell spreading dynamics, and slow down at the next stages.

††preprint: APS/123-QED

##  I Introduction

Understanding the phenomenon of cell spreading has numerous potential
applications, which include designing biomaterials for optimal control of cell
behavior Gong  _et al._ (2018), insight into cell morphology Folkman and
Moscona (1978), and developing efficient methods for gene transfection in
biomaterials Yang  _et al._ (2019).

In the last two decades, several models have been developed that describe cell
spreading on a flat substrate. Those models take into account the elastic Li
_et al._ (2015), Xiong  _et al._ (2010) or visco-elastic properties of the
cell and/or the substrate Gong  _et al._ (2018), Nisenholz  _et al._ (2014).
In addition, those models describe the dynamics of some subcellular components
such as cortical cytoskeleton, cell nuclear, actin filaments, and microtubules
Fang and Lai (2016), Vernerey and Farsad (2014). Also, they consider the
mechanical interactions between cell adhesion molecules like cross-membrane
protein, molecular clutches, and the extracellular property of the substrate.
There exist other models that describe cell spreading on non-flat substrate
such as V or Y-shaped micro-patterned substrates McEvoy  _et al._ (2017). All
of the previous models were compared with and validated by experimental
measurements.

Typically, the implementation of computational models needs hard numerical
simulations based on finite elements methods Vernerey and Farsad (2014),
Odenthal  _et al._ (2013) or minimizing some free energy functionals Fang and
Lai (2016); some models include stochastic effects McEvoy  _et al._ (2017).

Based on experimental data and measurements, some universality property of
cell spreading have been discovered. Usually early spreading is isotropic.
Cell spreading may experience three sequential phases, basal (cell touches the
substrate), fast continuous spreading (generation of lamellipodial sheet), and
periodic local contractile spreading Dobereiner  _et al._ (2004). These phases
obey a power-law area growth with distinct exponents when we plot cell
adhesion area (contact area) versus time. Later the authors in Cuvelier  _et
al._ (2007) succeed to explain these power-law relationships with a relatively
simple physical model. They consider energy balance and assume that actin
cortex is a viscous liquid McGrath (2007).

A minimal computational phase field model of 3D cell crawling on general
substrate topography was developed in Winkler  _et al._ (2019). In the present
paper we consider a simplified version of that model. We choose the substrate
to be a flat surface, $z=0$, in order to model the dynamics of cell spreading
on the plane. Unlike the models mentioned above, our model is simple. We
describe the cell spreading dynamics by a single scalar non-local partial
differential equation of the cell interface (16), which could be solved easily
with Wolfram Mathematica program. Our model is in qualitative agreement with
observations at the early fast phase, and the universal power law at the
earlier stages of cell spreading.

The structure of the paper is as follows. In Sec. II we present the minimal 3D
phase field model. In Sec. III we introduce the proper length and time scales
of the spreading dynamics. We perform asymptotic analysis and find the fields
at the leading order, and then we use the solvability condition to derive a
closed evolutionary nonlocal equation that describes the cell interface
dynamics (16). We solve this equation numerically via the function `NDSolve`
of Wolfram Mathematica. We reveal the agreement with the universal power law.
Finally, in Sec. IV we present the conclusions.

##  II Formulation of the problem

In order to describe the dynamics of cell located in region $z>0$ and
spreading on the flat substrate $z=0$, see Fig. 1, we extend the problem into
the whole space, postulating the reflection symmetry or antisymmetry of our
fields under the transformation $z\to-z$.

Let us consider the following simplified version of the model that was
formulated in Winkler  _et al._ (2019):

$\displaystyle u_{t}=D_{u}\nabla^{2}u-(1-u)(\delta-u)u-$
$\displaystyle\alpha\nabla u\cdot\textbf{P}-k\nabla\Psi_{u}\cdot\nabla u,$
(1a) $\displaystyle\delta=\frac{1}{2}+\mu\delta V-\sigma|\textbf{P}|^{2},\
\delta V(t)=\int u\textrm{d}^{3}r-v_{0},$ (1b)
$\displaystyle\textbf{P}_{t}=D_{p}\nabla^{2}\textbf{P}-\tau^{-1}\textbf{P}-$
$\displaystyle\beta\Psi_{p}(z)\left[(1-\nu)\hat{P}\nabla u+\nu\nabla
u\right],$ (1c) $\displaystyle u(r=0)=1,\quad u(r\rightarrow\infty)=0,$ (1d)
$\displaystyle\textbf{P}(r=0)=\textbf{P}(r\rightarrow\infty)=0;$ (1e)

see Fig. 1; here $u$ is the order parameter that is close to $1$ inside the
cell and $0$ outside, and P is the three-dimensional polarization vector field
representing the actin orientations. In (1c), $\hat{P}=\hat{I}-\hat{n}\hat{n}$
is the projection operator onto the local tangential plane, where
$\hat{n}=\nabla\Psi_{p}/|\nabla\Psi_{p}|$ (in our case $\hat{n}=\pm\hat{z}$).
Therefore,

$\hat{P}\nabla u=\nabla u-\hat{z}\hat{z}\nabla u=u_{x}\hat{x}+u_{y}\hat{y}.$

The parameter $0\leq\nu\leq 1$ model the contribution of actin polarization
from the tangential limit $\nu=0$, and the isotropic limit $\nu=1$.

The constant parameters of the problem are: $D_{u}$ is the stiffness of
diffuse interface, $D_{p}$ is the diffusion coefficient for P, $\alpha$ is the
coefficient characterizing advection of $u$ by P, $\beta$ determines the
creation of P at the interface, $\tau^{-1}$ is the inverse time of the
degradation of P inside the cell, $v_{0}$ is twice the overall initial volume
of the cell due to the reflection symmetry, $\mu$ is the stiffness of the
volume constraint, and $\sigma$ is the contractility of actin filament
bundles. All the parameters listed above are positive.

In addition, we make the basic assumption that the ratio $\epsilon$ of the
thickness of the cell wall (i.e., the width of the transition zone, where $u$
is changed from nearly 1 to nearly 0) to the characteristic size of the cell
is small, $\epsilon\ll 1$, see Fig. 1.

Figure 1: A schematic description of cell spreading dynamics in spherical
coordinate system, and the boundary conditions of the simplified model (1).
The cell is only the shape in the upper region $z\geq 0$ while the model (1)
is formulated in the hole space with reflection symmetry assumption with
respect of the substrate plane $z=0$, (4). The thickness of the cell wall
(i.e., the width of the transition zone, where $u(r,t)$ is changed from nearly
1 to nearly 0) is $O(1)$, and the cell size is large. Therefore the ratio
$\epsilon$ of the thickness of the cell wall to the size of the cell is small.
In addition we emphasize the regions where the fields
$\nabla\Psi_{u}\cdot\nabla u$, and $\Psi_{p}\nabla u$ give their main
contribution and exponentially small otherwise. These are the regions where
protrusions developed.

Motivated by Winkler  _et al._ (2019), we define the static fields as

$\Psi_{u}(z)=\textrm{e}^{-(\epsilon
z)^{2}/D_{u}},\quad\Psi_{p}(z)=\textrm{e}^{-(\epsilon z)^{2}/\l_{p}},\ \
l_{p}=\tau D_{p},$ (2)

see Fig. 2. We choose $l_{p}>D_{u}$ to allow more substantial actin inside the
cell. The appearance of $\epsilon$ in the exponents of (2) allows to avoid
boundary layer problem complications both in time and space. The expression
$k\nabla\Psi_{u}\cdot\nabla u$ in (1a) models the adhesion effect of the
substrate; $k$ is adhesion strength parameter. Notice that
$\nabla\Psi_{u}\cdot\nabla u=O(1)$ only in region nearby the substrate and
also at the cell boundary or membrane i.e., where protrusion holds. Also the
appearance of $\Psi_{p}\nabla u$ in (1c) allows high actin concentration
nearby the flat substrate $z=0$ where protrusions are developed during cell
spreading, and low actin concentration otherwise, see Fig. 1. This scenario is
in agreement with experimental studies, see EL  _et al._ (2011), Li  _et al._
(2015), and the review paper Mattila and Lappalainen (2008).

We apply the spherical coordinate system, see Fig. 1, hence
$u=u(r,\theta,\varphi,t)$,
$\textbf{P}(r,\theta,\varphi,t)=p\hat{r}+q\hat{\theta}+w\hat{\varphi}$. We
define the iso-surface of the interface as $u(r=\rho(\theta,\varphi,t))=1/2$.
As a result of our definition of the spherical coordinate, we have

$\displaystyle(r,\theta,\varphi),\ 0<r<\infty,\
-\frac{\pi}{2}<\theta<\frac{\pi}{2},\ 0<\varphi<2\pi,$ (3a) $\displaystyle
x=r\cos\theta\cos\varphi,\ y=r\cos\theta\sin\varphi,\ z=r\sin\theta,$ (3b)
$\displaystyle\nabla=\hat{r}\partial_{r}+\hat{\theta}\frac{\partial_{\theta}}{r}+\hat{\varphi}\frac{\partial_{\varphi}}{r\cos\theta},$
(3c)
$\displaystyle\nabla^{2}=\partial_{r}^{2}+\frac{2\partial_{r}}{r}+\frac{1}{r^{2}}\Big{(}\partial_{\theta}^{2}-\tan\theta\partial_{\theta}+\frac{\partial_{\varphi}^{2}}{\cos^{2}\theta}\Big{)},$
(3d)
$\displaystyle\nabla^{2}\textbf{P}=\hat{r}\nabla^{2}p+\hat{\theta}\nabla^{2}q+\hat{\varphi}\nabla^{2}w+O\left(\frac{1}{r^{2}}\right).$
(3e)

One can calculate,

$\displaystyle\left[(1-\nu)\hat{P}+\nu\hat{I}\right]\nabla
u=u_{x}\hat{x}+u_{y}\hat{y}+\nu u_{z}\hat{z}=\nabla u+(\nu-1)u_{z}\hat{z}$
$\displaystyle=\left[(1+(\nu-1)\sin^{2}\theta)u_{r}+(\nu-1)\frac{\sin
2\theta}{2r}u_{\theta}\right]\hat{r}+$
$\displaystyle\left[(1+(\nu-1)\cos^{2}\theta)\frac{u_{\theta}}{r}+(\nu-1)\frac{\sin
2\theta}{2}u_{r}\right]\hat{\theta}+\frac{u_{\varphi}}{r\cos\theta}\hat{\varphi}.$

Notice that due to the appearance of $\Psi_{u,p}(z)$, the system (1a)-(1e)
does not have any rotationally spherical symmetric solutions. Therefore, we
have to look for general shape solutions.

Figure 2:  Plot of the density functions $\Psi_{u,p}$ in (2). Notice that
$\Psi_{u,p}=O(1)$ nearby the substrate $z=0$ (red region) while attenuate far
away (blue region) . The width of the layer where $\Psi_{p}=O(1)$ is
$O(\sqrt{\tau D_{p}})$. Therefore it is thicker than layer where
$\Psi_{u}=O(1)$ that is $O(\sqrt{D_{u}})$.

The reflection symmetry assumption relative to the substrate plane $z=0$,
yields the conditions,

$\displaystyle u(-\theta)=u(\theta),\quad\rho(-\theta)=\rho(\theta),$ (4a)
$\displaystyle p(-\theta)=p(\theta),\quad w(-\theta)=w(\theta),\quad
q(-\theta)=-q(\theta).$ (4b)

##  III Dynamics of general shape interface

In order to balance the front dynamics with curvature we impose the following
scaling that describes slow dynamics of a large-size cell Hamed and
Nepomnyashchya (2020),

$\tilde{t}=\epsilon^{2}t,\quad\rho(\theta,\varphi,t)=\epsilon^{-1}R(\theta,\varphi,t),\quad\epsilon\ll
1.$ (5)

The transition zone variable is defined as

$\zeta=r-\rho(\theta,\varphi,t)=O(1).$ (6)

Also we define,

$\displaystyle R(\theta,\varphi,t)=\tilde{R}(\theta,\varphi,\tilde{t}),\
u(r,\theta,\varphi,t)=\tilde{u}(\zeta,\theta,\varphi,\tilde{t}),$ (7a)
$\displaystyle\textbf{P}(r,\theta,\varphi,t)=\tilde{\textbf{P}}(\zeta,\theta,\varphi,\tilde{t}).$
(7b)

The chain rule yields

$\partial_{t}=-\epsilon\tilde{R}_{\tilde{t}}\partial_{\zeta}+\epsilon^{2}\partial_{\tilde{t}},\quad\partial_{r}=\partial_{\zeta}.$
(8)

later on we drop the tildes. It holds that

$\displaystyle\frac{1}{r}=\frac{\epsilon}{R}-\frac{\epsilon^{2}\zeta}{R^{2}}+...,$
(9a) $\displaystyle\Psi_{p}(z)=\textrm{e}^{-(\epsilon
z)^{2}/l_{p}}=\textrm{e}^{-\epsilon^{2}r^{2}\sin^{2}\theta/l_{p}}=$
$\displaystyle\textrm{e}^{-\epsilon^{2}(\zeta+\epsilon^{-1}R)^{2}\sin^{2}\theta/l_{p}}\sim\textrm{e}^{-R^{2}\sin^{2}\theta/l_{p}},$
(9b) $\displaystyle\nabla\Psi_{u}\sim\frac{-\epsilon
R}{D_{u}}\textrm{e}^{-R^{2}\sin^{2}\theta/D_{u}}\left(2\sin^{2}\theta\hat{r}+\sin
2\theta\hat{\theta}\right).$ (9c)

In addition one can calculate,

$\displaystyle\partial_{\theta}=-\epsilon^{-1}R_{\theta}\partial_{\zeta}+\partial_{\theta},\
\partial_{\varphi}=-\epsilon^{-1}R_{\varphi}\partial_{\zeta}+\partial_{\varphi},$
(10a)
$\displaystyle\partial_{\theta}^{2}=\epsilon^{-2}R_{\theta}^{2}\partial_{\zeta}^{2}-\epsilon^{-1}(R_{\theta\theta}\partial_{\zeta}+2R_{\theta}\partial_{\zeta\theta}^{2})+\partial_{\theta}^{2},$
(10b)
$\displaystyle\partial_{\varphi}^{2}=\epsilon^{-2}R_{\varphi}^{2}\partial_{\zeta}^{2}-\epsilon^{-1}(R_{\varphi\varphi}\partial_{\zeta}+2R_{\varphi}\partial_{\zeta\varphi}^{2})+\partial_{\varphi}^{2},$
(10c)
$\displaystyle\frac{\partial_{\theta}}{r}=-\frac{R_{\theta}}{R}\partial_{\zeta}+O(\epsilon),\
\frac{\partial_{\theta}}{r^{2}}=-\epsilon\frac{R_{\theta}}{R^{2}}\partial_{\zeta}+O(\epsilon^{2}),$
(10d)
$\displaystyle\frac{\partial_{\varphi}}{r}=-\frac{R_{\varphi}}{R}\partial_{\zeta}+O(\epsilon),\
\frac{\partial_{\varphi}}{r^{2}}=-\epsilon\frac{R_{\varphi}}{R^{2}}\partial_{\zeta}+O(\epsilon^{2}),$
(10e)
$\displaystyle\nabla^{2}u=\left(1+\frac{R_{\theta}^{2}}{R^{2}}+\frac{R_{\varphi}^{2}}{R^{2}\cos^{2}\theta}\right)u_{\zeta\zeta}+$
$\displaystyle\epsilon\Bigg{[}\left(\frac{2}{R}+\frac{R_{\theta}}{R^{2}}\tan\theta-\frac{R_{\theta\theta}}{R^{2}}-\frac{R_{\varphi\varphi}}{R^{2}\cos^{2}\theta}\right)u_{\zeta}$
$\displaystyle-\frac{2}{R^{2}}\left(R_{\theta}u_{\zeta\theta}+\frac{R_{\varphi}}{\cos^{2}\theta}u_{\zeta\varphi}\right)$
$\displaystyle-\frac{2\zeta}{R^{3}}\left(R_{\theta}^{2}+\frac{R_{\varphi}^{2}}{\cos^{2}\theta}\right)u_{\zeta\zeta}\Bigg{]}+O(\epsilon^{2}).$
(10f)

We can approximate the nonlocality in (1b) as follows,

$\int
u\textrm{d}^{3}r\sim\frac{\epsilon^{-3}}{3}\int_{0}^{2\pi}\textrm{d}\varphi\int_{-\pi/2}^{\pi/2}R^{3}(\theta,\varphi,t)\cos\theta\textrm{d}\theta.$
(11)

Consider the following scaling of the model parameters

$\alpha=\epsilon A,\ \frac{4\pi\mu}{3}\epsilon^{-3}=\epsilon M,\
\sigma=\epsilon S,\ v_{0}=\epsilon^{-3}V_{0},\ k=O(1).$ (12)

Let us introduce the expansions

$u=u_{0}+\epsilon u_{1}+...,\quad p=p_{0}+\epsilon p_{1}+...$ (13)

We define the auxiliary function ,

$\Lambda(\theta,\varphi,t)=\left(1+\frac{R_{\theta}^{2}}{R^{2}}+\frac{R_{\varphi}^{2}}{R^{2}\cos^{2}\theta}\right)^{-1/2},$

and the function,

$\displaystyle\Phi(\tau,D_{u},D_{p},\zeta)=$
$\displaystyle\frac{1}{8}\sqrt{\frac{\tau}{2D_{u}D_{p}}}\int_{-\infty}^{\infty}\textrm{e}^{-|s|/\sqrt{\tau
D_{p}}}\cosh^{-2}\left(\frac{s-\zeta}{\sqrt{8D_{u}}}\right)\textrm{d}s,$

that are basic for our next analysis.

We substitute the length, time (5), and the parameters scaling (12) into
system (1). We write the system (1) in the transition zone variable (6)-(7)
via the expansions and the chain rules (8)-(11). We substitute the asymptotic
expansions (13) and finally we collect terms of the same order.

Consequently we obtain at the leading order the following system,

$\displaystyle
D_{u}\Lambda^{-2}u_{0\zeta\zeta}=(1-u_{0})(\frac{1}{2}-u_{0})u_{0},$
$\displaystyle D_{p}\Lambda^{-2}p_{0\zeta\zeta}-\tau^{-1}p_{0}=$
$\displaystyle\beta\textrm{e}^{-R^{2}\sin^{2}\theta/l_{p}}\left[1+(\nu-1)\sin^{2}\theta-(\nu-1)\frac{\sin
2\theta}{2}\frac{R_{\theta}}{R}\right]u_{0\zeta},$ $\displaystyle
D_{p}\Lambda^{-2}q_{0\zeta\zeta}-\tau^{-1}q_{0}=$
$\displaystyle-\beta\textrm{e}^{-R^{2}\sin^{2}\theta/l_{p}}\left[(1+(\nu-1)\cos^{2}\theta)\frac{R_{\theta}}{R}-(\nu-1)\frac{\sin
2\theta}{2}\right]u_{0\zeta},$ $\displaystyle
D_{p}\Lambda^{-2}w_{0\zeta\zeta}-\tau^{-1}w_{0}=-\beta\textrm{e}^{-R^{2}\sin^{2}\theta/l_{p}}\frac{R_{\varphi}}{R\cos\theta}u_{0\zeta}.$

Following the Ginzburg-Landau theory and Fourier transform method the solution
of this system is given by,

$\displaystyle
u_{0}(\zeta)=\frac{1}{2}\left[1-\tanh\left(\frac{\Lambda\zeta}{\sqrt{8D_{u}}}\right)\right],$
$\displaystyle p_{0}(\zeta)=\beta\lambda_{p}\Lambda\Phi(\Lambda\zeta),$
$\displaystyle\lambda_{q}=\textrm{e}^{-R^{2}\sin^{2}\theta/l_{p}}\left[1+(\nu-1)\sin^{2}\theta-(\nu-1)\frac{\sin
2\theta}{2}\frac{R_{\theta}}{R}\right],$ $\displaystyle
q_{0}(\zeta)=-\beta\lambda_{q}\Lambda\Phi(\Lambda\zeta),$
$\displaystyle\lambda_{q}=\textrm{e}^{-R^{2}\sin^{2}\theta/l_{p}}\left[(1+(\nu-1)\cos^{2}\theta)\frac{R_{\theta}}{R}-(\nu-1)\frac{\sin
2\theta}{2}\right],$ $\displaystyle
w_{0}(\zeta)=-\beta\lambda_{w}\Lambda\Phi(\Lambda\zeta),$
$\displaystyle\lambda_{w}=\textrm{e}^{-R^{2}\sin^{2}\theta/l_{p}}\frac{R_{\varphi}}{R\cos\theta}.$

Notice that these results satisfy the symmetry conditions (4), if
$R(-\theta)=R(\theta)$.

The equation for the correction term $u_{1}$ at the order $O(\epsilon)$ have
the form,

$\displaystyle L[u_{1}]=\textrm{RHS},$ (15) $\displaystyle
L=D_{u}\Lambda^{-2}\partial_{\zeta}^{2}-\left(\frac{1}{2}-3u_{0}+3u_{0}^{2}\right)\hat{I},$
$\displaystyle\textrm{RHS}=-R_{t}u_{0\zeta}-$ $\displaystyle
D_{u}\Bigg{[}\left(\frac{2}{R}+\frac{R_{\theta}}{R^{2}}\tan\theta-\frac{R_{\theta\theta}}{R^{2}}-\frac{R_{\varphi\varphi}}{R^{2}\cos^{2}\theta}\right)u_{0\zeta}$
$\displaystyle-\frac{2}{R^{2}}\left(R_{\theta}u_{0\zeta\theta}+\frac{R_{\varphi}}{\cos^{2}\theta}u_{0\zeta\varphi}\right)$
$\displaystyle-\frac{2\zeta}{R^{3}}\left(R_{\theta}^{2}+\frac{R_{\varphi}^{2}}{\cos^{2}\theta}\right)u_{0\zeta\zeta}\Bigg{]}$
$\displaystyle+A\left(p_{0}-\frac{R_{\theta}}{R}q_{0}-\frac{R_{\varphi}}{R\cos\theta}w_{0}\right)u_{0\zeta}$
$\displaystyle-\frac{k}{D_{u}}R\cdot\textrm{e}^{-R^{2}\sin^{2}\theta/D_{u}}\left(2\sin^{2}\theta-\sin(2\theta)\frac{R_{\theta}}{R}\right)u_{0\zeta}+$
$\displaystyle(1-u_{0})u_{0}\left\\{\tilde{V}(t)-S(p_{0}^{2}+q_{0}^{2}+w_{0}^{2})\right\\},$

where the volume variation have the form,

$\tilde{V}(t)=M\left[\frac{1}{4\pi}\int_{0}^{2\pi}\textrm{d}\varphi\int_{-\pi/2}^{\pi/2}R^{3}(\theta,\varphi,t)\cos\theta\textrm{d}\theta-\frac{3}{4\pi}V_{0}\right].$

We apply the solvability condition, which is the orthogonality of the right-
hand side (RHS) of equation (15) to the solution $u_{0\zeta}$ of the
homogenous equation $L[u]=0$ of (15) i.e.,

$\int_{-\infty}^{\infty}\textrm{RHS}(\zeta)\cdot
u_{0\zeta}(\zeta)\textrm{d}\zeta=0.$

We therefore obtain a closed equation governing the interface dynamics
$R(\theta,\varphi,t)$,

$a\Lambda R_{t}=-2aD_{u}\mathcal{H}-\tilde{V}+\Omega-N,$ (16)

where

$\mathcal{H}=\frac{1}{2}\nabla\cdot\hat{n}=\frac{1}{2}\nabla\cdot\left(\frac{\nabla(r-R)}{|\nabla(r-R)|}\right)$
(17)

is the mean local curvature of the surface $r=R(\theta,\varphi,t)$, see
Appendix A, and

$\displaystyle\Omega(\theta,\varphi,t)=6\beta
A\Omega_{1}\Lambda^{2}\left(\lambda_{p}+\frac{R_{\theta}}{R}\lambda_{q}+\frac{R_{\varphi}}{R\cos\theta}\lambda_{w}\right)+$
$\displaystyle
6\beta^{2}S\Omega_{2}\Lambda^{2}\left(\lambda_{p}^{2}+\lambda_{q}^{2}+\lambda_{w}^{2}\right),$
$\displaystyle\Omega_{1}(\tau,D_{u},D_{p})=\int_{-\infty}^{\infty}\Phi(\xi)\bar{u}_{0\xi}^{2}(\xi)\textrm{d}\xi,$
$\displaystyle\Omega_{2}(\tau,D_{u},D_{p})=$
$\displaystyle\int_{-\infty}^{\infty}\Phi^{2}(\xi)(\bar{u}_{0}(\xi)-1)\bar{u}_{0}(\xi)\bar{u}_{0\xi}(\xi)\textrm{d}\xi>0,$

where

$\bar{u}_{0}(\xi)=\frac{1}{2}\left[1-\tanh\left(\frac{\xi}{\sqrt{8D_{u}}}\right)\right].$

The adhesion effect is implemented by

$N=\frac{ak}{D_{u}}\Lambda
R\cdot\textrm{e}^{-R^{2}\sin^{2}\theta/D_{u}}\left(2\sin^{2}\theta-\sin(2\theta)\frac{R_{\theta}}{R}\right)$

For more details and explanations about arguments for the derivation of the
governing equation (16) see Hamed and Nepomnyashchya (2016), Hamed and
Nepomnyashchya (2020), and Hamed and Nepomnyashchya (2021).

Equation (16) comes in conjunction with the Neumann boundary conditions

$R_{\theta}(\theta=0)=R_{\theta}(\theta=\pi/2)=0,$ (18)

and some initial interface $R(t=0)$.

In order to model cell spreading, we may consider the axi-symmetric case
$R_{\varphi}=0$, since according to experimental observation, the onset of
cell spreading is isotropic Dobereiner  _et al._ (2004). As for the initial
interface shape, we take the truncated sphere with radius $R_{0}$ and center
$(0,0,\eta)$, see Fig. 3(a).

$R(t=0)=\eta\sin\theta+\sqrt{R_{0}^{2}-\eta^{2}\cos^{2}\theta},\ \
0\leq\eta\leq 1.$ (19)

Then the initial normalized volume in this case is given by

$V_{0}=\int_{0}^{\pi/2}R(t=0)^{3}\cos\theta\textrm{d}\theta,$

which is twice the volume of the truncated sphere.

In Fig. 3(a)-(f) we present the sequence of plots that show the results of the
numerical simulation of the interface $R(\theta,\varphi,t)$ according to
equation (16) and boundary condition (18). We use the function `NDSolve` of
Wolfram Mathematica. We plot $R(\theta,\varphi,t)$ only in the upper region
$0\leq\theta\leq\pi/2$, the plot in the lower region is only a mirror
reflection of the upper surface due to our symmetry assumption (4). Following
experimental scenarios where the spherical-like cell almost touch the
substrate we may take the parameters of our initial interface as $\eta=0.95$,
and $R_{0}=1$, see Fig. 3(a). This simulation as we see describes cell
spreading. We begin from almost full sphere and end up with ellipsoid-like
shape which is the steady state solution of the system (16), (18), and (19).
In Fig. 3(g) we display the cell hight $R(\pi/2,t)$ which decreases from
almost 2 to 1.21, while in Fig. 3(h) we display the cell contact area (radius)
$R(0,t)$ which increases to 2.

In addition Fig. 3(g)-(h) display the fast spreading phenomena at the
beginning of cell spreading in agreement with the continuous spreading fast
phase that was observed experimentally in Dobereiner  _et al._ (2004). Also,in
Fig 3(i) we consider the Log-Log plot of the cell radius versus time, also we
plot the piecewise function that connect two function of the form
$b_{1}t^{1/2}$, and $b_{2}t^{1/4}$, for a proper choice of the parameters
$b_{1,2}$, and for the connecting point. We notice the agreement with the
universal power law Cuvelier  _et al._ (2007) that suggest that cell adhesion
or contact area versus time behave as $\sim t^{1/2}$ in the early state of
cell spreading dynamics, and slow down in the next states. The plot of the
slope $\sim t^{1/4}$ in Fig. 3(i) is only to emphasize the slowing down of the
next phase.

In Fig. 4 we perform similar analysis where we choose $D_{u}=0.5$, and
$D_{p}=0.02$ while the other parameters remain as those of Fig. 3. Notice that
in this case we have $\tau D_{p}<D_{u}$, unlike the previous case of Fig. 3

##  IV Conclusion

We utilize a simplified version of minimal 3D phase field model that was
developed in Winkler  _et al._ (2019), in order to model cell spreading
dynamics on a flat substrate. The model (1) couples the order parameter $u$
with 3D polarization (orientation) vector field P of the actin network. The
model is formulated in the whole space but with appropriate symmetry
conditions with respect to transformation $z\rightarrow-z$, (4).

After we introduce the proper time and length scale and perform asymptotic
expansion, we solve equations for the fields at the leading order. As a result
of the solvability condition we derive a closed integro-differential equation
(16) governing the 3D cell spreading dynamics, which includes the normal
velocity $\omega_{n}=\Lambda R_{t}$ of the membrane, curvature $\mathcal{H}$,
volume relaxation rate $\tilde{V}$, a function $\Omega(t)$ determined by the
molecular effects of the subcell level, and the adhesion effect $N$.

Excluding the adhesion effect this result is similar to the 2D case which
describe the onset of 2D cell dynamics on flat substrate Hamed and
Nepomnyashchya (2020) and 3D case that describe the onset of 3D cell motility
immersed in 3D extracellular matrices Hamed and Nepomnyashchya (2021).

The equation governing the interface or membrane dynamics during spreading may
be presented in the form:

$\omega_{n}=-2D_{u}\mathcal{H}-\tilde{V}+\Omega-N,$

after we put the proper scaling transformation, $t\rightarrow a^{2}t$ and
$R(t)\rightarrow aR(t)$ in (16). This equation is easily solved numerically
via the function `NDSolve` of Wolfram Mathematica. The simulation present cell
spreading with significant hight decreasing and radius increasing of the
initial truncated spheres, see Fig 3, Fig. 4

These results are in agreement with the early fast phase that was observed
experimentally in Dobereiner  _et al._ (2004). Surprisingly, the result are in
qualitative agreement with universal power law which suggest that adhesion or
contact area versus time behave as $\sim t^{1/2}$ in the early state of cell
spreading dynamics ,and then it slow down. The appearance of the slope $\sim
t^{1/4}$ in Fig. 3(i) is only to emphasize the slowing down of the later
phase. This is a surprising result since in our phase field model we did not
assume any viscosity property of the cell membrane as it is assumed in
Cuvelier  _et al._ (2007).

## Appendix A

Here we give an explicit expression for the mean curvature of a surface given
in spherical coordinate description $r=R(\theta,\varphi,t)$, see Fig. 1.
Following the definition (17),(3), one can calculate,

$\displaystyle\nabla\cdot\hat{n}=\Lambda\Bigg{\\{}\frac{2}{R}-\frac{R_{\theta\theta}}{R^{2}}+\frac{R_{\theta}\tan\theta}{R^{2}}-\frac{R_{\varphi\varphi}}{R^{2}\cos^{2}\theta}\Bigg{\\}}$
$\displaystyle+\Lambda^{3}\Bigg{\\{}\frac{1}{R^{3}}\left(R_{\theta}^{2}+\frac{R_{\varphi}^{2}}{\cos^{2}\theta}\right)+\frac{R_{\theta}^{2}R_{\theta\theta}}{R^{4}}+\frac{2R_{\theta}R_{\varphi}R_{\theta\varphi}}{R^{4}\cos^{2}\theta}$
$\displaystyle+\frac{R_{\varphi}^{2}}{R^{4}\cos^{4}\theta}\left(\frac{1}{2}R_{\theta}\sin
2\theta+R_{\varphi\varphi}\right)\Bigg{\\}}$ (20)

Figure 3: (a)-(f) Time sequence of the simulation of the axisymmetric case of
equation (16) with initial upwards shifted sphere (19) with radius $R_{0}=1$
and center $(0,0,0.95)$. We employ the following value of parameters
$\beta=5,A=1,\tau=10,D_{u}=1,D_{p}=0.2,M=8,S=2,\nu=0.5,k=15$. (g) is the plot
of the ellipsoid like hight $R(\pi/2,t)$, (h) is the plot of the ellipsoid
like radius $R(0,t)$, both in the time interval $0\leq t\leq 0.1$. Notice the
stationary ellipsoid-like has hight 1.21 and radii 2. (i) Log-Log plot of the
cell radius $R(0,t)$ verses time $t$. Also we plot the piecewise function that
connect two function of the form $b_{1}t^{1/2}$, and $b_{2}t^{1/4}$, for a
proper choice of the parameters $b_{1,2}$, and for the connecting point.
Notice the qualitative agreement with the universal power law in the initial
fast phase and next the slower phase. The appearance of the slope $\sim
t^{1/4}$ in (i), is shown only to emphasize the slowing down of the later
phase.

Figure 4: (a)-(h) Time sequence of the simulation of the axi symmetric case of
equation (16) . We employ following value of parameters
$D_{u}=0.5,D_{p}=0.02$, the other parameters and the initial sphere are as in
Fig 3. (i) is the plot of the ellipsoid like hight $R(\pi/2,t)$, (j) is the
plot of the ellipsoid like radius $R(0,t)$, both in the time interval $0\leq
t\leq 0.1$. Notice the stationary ellipsoid-like has hight 1.3 and radii 2.74.
Notice that here we have $\tau D_{p}<D_{u}$, unlike the previous case of Fig.
3. (k) Log-Log plot as that of Fig. 3.

## References

  * Winkler  _et al._ (2019) B. Winkler, I. S. Aranson, and F. Ziebert, “Confinement and substrate topography control cell migration in a 3d computational model,” Communications Physics 2:82 (2019).
  * Dobereiner  _et al._ (2004) Dobereiner, Dubin-Thaler, Giannone, Xenias, and Sheetz, “Dynamic phase transitions in cell spreading,” Physical review letters 93 (2004).
  * Cuvelier  _et al._ (2007) Cuvelier, The´ry, Chu, Dufour, Thie´ry, Bornens, Nassoy, and Mahadevan, “The universal dynamics of cell spreading,” Current Biology 17, 694–699 (2007).
  * Gong  _et al._ (2018) Z. Gong, S. E. Szczesny, S. R. Caliarie, E. E. Charrier, O. Chaudhuri, X. Cao, Y. Lin, R. L. Mauck, P. A. Janmey, J. A. Burdick, and V. B. Shenoy, “Matching material and cellular timescales maximizes cell spreading on viscoelastic substrates,” PNAS 115 (12), E2686–E2695 (2018).
  * Folkman and Moscona (1978) J. Folkman and A. Moscona, “Role of cell shape in growth control,” Nature 273, 345–349 (1978).
  * Yang  _et al._ (2019) Y. Yang, X. Wang, X. Hu, N. Kawazoe, Y. Yang, and G. Chen, “Influence of cell morphology on mesenchymal stem cell transfection,” ACS Appl. Mater. Interfaces 11, 1932−–1941 (2019).
  * Li  _et al._ (2015) Y. Li, D. Lovett, Q. Zhang, S. Neelam, R. A. Kuchibhotla, R. Zhu, G. G. Gundersen, T. P. Lele, and R. B. Dickinson, “Moving cell boundaries drive nuclear shaping during cell spreading,” Biophysical Journal 109 (2015).
  * Xiong  _et al._ (2010) Y. Xiong, P. Rangamani, M. A. Fardin, A. Lipshtat, B. D. Thaler, O. Rossier, M. P. Sheetz, and R. Iyengar, “Mechanisms controlling cell size and shape during isotropic cell spreading,” Biophysical Journal 98, 2136–2146 (2010).
  * Nisenholz  _et al._ (2014) N. Nisenholz, K. Rajendran, Q. Dang, H. Chen, R. Kemkemer, R. Krishnan, and A. Zemel, “Active mechanics and dynamics of cell spreading on elastic substrates,” Soft Matter 10 (2014).
  * Fang and Lai (2016) Y. Fang and K. W. C. Lai, “Modeling the mechanics of cells in the cell–spreading process driven by traction forces,” Physical Review E 93 (2016).
  * Vernerey and Farsad (2014) F. J. Vernerey and M. Farsad, “A mathematical model of the coupled mechanisms of cell adhesion, contraction and spreading,” J Math Biol 68(4), 989–1022 (2014).
  * McEvoy  _et al._ (2017) E. McEvoy, V. S. Deshpande, and P. McGarry, “Free energy analysis of cell spreading,” Journal of the Mechanical Behavior of Biomedical Materials 74, 283–295 (2017).
  * Odenthal  _et al._ (2013) T. Odenthal, B. Smeets, P. V. Liedekerke, E. Tijskens, H. V. Oosterwyck, and H. Ramon, “Analysis of initial cell spreading using mechanistic contact formulations for a deformable cell model,” PLOS Computational Biology 9, e1003267 (2013).
  * McGrath (2007) J. L. McGrath, “Dispatch: Cell spreading: The power to simplify,” Current Biology 17, R358 (2007).
  * EL  _et al._ (2011) Barnhart EL, Lee K-C, Keren K, Mogilner A, and Theriot JA, “An adhesion–dependent switch between mechanisms that determine motile cell shape,” PLoS Biol 9(5) (2011).
  * Mattila and Lappalainen (2008) P. K. Mattila and P. Lappalainen, “Filopodia: molecular architecture and cellular functions,” Nature Publishing Group 9 (2008).
  * Hamed and Nepomnyashchya (2020) M. Abu Hamed and A.A. Nepomnyashchya, “A simple model of keratocyte membrane dynamics: The case of motionless living cell,” Physica D 408 (2020).
  * Hamed and Nepomnyashchya (2016) M. Abu Hamed and A.A. Nepomnyashchya, “Dynamics of curved fronts in systems with power-law memory,” Physica D 328–329, 1–8 (2016).
  * Hamed and Nepomnyashchya (2021) M. Abu Hamed and A.A. Nepomnyashchya, “Three–dimensional phase field model for actin–based cell membrane dynamics,” submitted to The European Journal of Applied Mathematics (2021).