# Thermodynamically Consistent Diffuse Interface Model for Cell Adhesion and
Aggregation

Lingyue Shen  Ping Lin  Zhiliang Xu  Shixin Xu<EMAIL_ADDRESS>Department of Mathematics, University of Dundee, Dundee DD1 4HN, Scotland,
UnitedKingdom.  Department of Applied and Computational Mathematics and
Statistics, University of Notre Dame,102G Crowley Hall, Notre Dame, IN 46556
Duke Kunshan University, 8 Duke Ave, Kunshan, Jiangsu, China.

###### Abstract

A thermodynamically consistent phase-field model is introduced for simulating
multicellular deformation, aggreggation under flow conditions. In particular,
a Lennard-Jones type potential is proposed under phase field framework for
cell-cell, cell-wall interactions. A second-order accurate in both space and
time $C^{0}$ finite element method is proposed to solved the model governing
equations. Various numerical tests confirm the convergence, energy stability,
and nonlinear mechanical properties of cells of the proposed scheme. Vesicles
with different adhesion are also used to explain the pathological risk for
patient with sickle cell disease.

###### keywords:

Vesicle; Aggreggation; Energy stable scheme; Cell-wall interaction.

##  1 Introduction

Studying interaction between structures (cell-cell and cell-blood-vessel) is
an important subject for understanding hemodynamics, because structural
interactions at the cellular level unambiguously appear in a broad spectrum of
blood flow related problems ranging from red blood cell (RBC) distribution [1]
in blood vessel, the growth of blood clot [2], blood cell aggregation [3],
sickle cell disease [4], tumor cell dynamics [5] and diabetes [6].

Under the framework of sharp interface description in which interfaces
separating different components of matter are idealized as hypersurfaces with
zero thickness, there are a number of studies focusing on modeling cell-cell
and cell-blood-vessel interactions under blood flow conditions. Following the
seminal work of Peskin [7, 8], the immersed boundary method was used to
develop a model for platelet aggregation in [9]. Cell-wall interaction model
is introduced in [5] to simulate the adhesion and deformation of tumor cells
at the vessel wall. Local and non-local models are described in [10] to
investigate the invasion and growth of tumor cells. Cell-cell interaction
modeling at the micro-scale is carried out in [11, 12, 13, 14] to study the
RBC aggregation problem. Multi-scale models are introduced in [15, 16] in
which dissipative particle dynamics (DPD) is used in [15] to establish a blood
cell model in blood flow and a stochastic cellular potts model (CPM) is
introduced in [16] for studying blood clot growth. In [17, 18, 19], to account
for cell-cell or cell-substrate adhesion, a Lennard-Jones type potential is
introduced as a one dimensional function of the distance between the points on
different cell membranes and substrates. The potential is a combination of a
repulsive part and an attractive part which shows repulsion when particles get
too close and shows attraction when the distance increase. In [18], this
potential in combination with DPD method as mentioned before, is utilized to
study cell deformation and doublet suspension.

In this work, we focus on utilizing diffuse interface method, also commonly
known as phase-field approach [20] to model cell-cell and cell-blood-vessel
interactions under blood flow conditions. The phase-field models replace the
sharp interface description with a thin transition region in which microscopic
mixing of the macroscopically distinct components of matter is permitted. The
phase-field approach not only yields systems of governing equations that are
better amenable to further analysis, but also handles topological changes of
the interface naturally [21, 20, 22, 23]. By treating the cell membrane as a
diffuse interface, several phase-field vesicle interaction models have been
reported recently. Fusion of cellular aggregates in biofabrication is modelled
with a single phase-field function in [24]. Vesicle-substrate adhesion model
is established in [25, 26, 27] where two phase-field functions are used to
indicate wall and vesicle, respectively. Vesicle-vesicle adhesion model is
introduced in [28] in which different phase-field functions are used to
represent different vesicles. Further, multiple cell aggregation model and
simulations are reported in [29, 30].

Numerically simulating cell models with phase-field method is challenging.
finite difference method is used in [31] to study the elastic bending energy
of the vesicle. Finite element method is introduced in [32, 26, 33, 29],
spectral methods with an adequate number of Fourier modes is used in [28,
34].However, few works developing energy stable numerical schemes are reported
for vesicle or cell models, despite the fact that most phase-field models obey
energy dissipation law at the continuous level. It is well-known that
preserving a discrete counterpart of the energy dissipation property of the
system in numerical solution is crucial for long term stability of numerical
computation [35, 36]. In this paper, we also try to develop a numerical scheme
retaining discrete energy dissipation law that is consistent with the
continuous one.

The purposes of this paper is twofold. The first goal of this paper is to
derive a thermodynamically consistent phase-field model for vesicles’ motion,
shape transformation and their structural interactions under flow conditions
in a closed spatial domain using an energy variational method [37, 38]. All
the physics that ones are interested in are taken into consideration through
definitions of the energy functional and the dissipation functional, together
with the kinematic relations of dynamic evolution of model state variables.

In this work, all terms of the energy and dissipation functionals are defined
on the bulk region of the domain including cell-cell adhesion and cell-wall
interactions. A Lennard-Jones typ potential within the framework of phase-
field method is derived for the vesicle-vesicle interaction, in which an
interacting potential including both repulsive and attractive parts is
established with respect to the phase-field order parameter. To consider the
cell-wall interaction, the wall is just treated as a new phase-field order
parameter or function $\phi_{w}(\bm{x})$ in the proposed Lennard-Jones type
potential. Then, the cell phase-field function $\phi$ defined on the bulk
region of the domain which would be used in the derivation of model governing
equations by the energy variational approach [39]. The advantage of doing so
is that the variational procedure only needs to consider the bulk unlike
previous works [37, 39] where the boundary needs to be considered since the
wall effect are included in boundary conditions of the system, which
simplifies the model derivation considerably.

Performing variation of the energy and dissipation functionals yields a Cahn-
Hilliard-Navier-Stokes (CH-NS) system with Allen-Cahn (AC) general Navier
boundary conditions (GNBC) [40, 39]. This is in contrast to most previous
works [41, 33, 42] in which boundary effects were rarely derived during the
course of model derivation. Dirichlet or Neumann type boundary conditions were
simply added to these models at the end to close the governing equations [43,
22, 23, 35, 44].

Moreover, our model accounts for the incompressibility of the fluid, and the
local and global inextensibility of the vesicle membrane by introducing two
Lagrangian multipliers, hydrostatic pressure $p$ and surface pressure
$\lambda$ [45] and penalty terms, respectively.

The second goal of this paper is to design an energy stable finite element
scheme for solving model governing equations. This scheme preserves a discrete
counterpart of the energy dissipation property of the model in numerical
solution, which is crucial for long term stability of numerical computation.
We note that although this scheme is designed for solving the model equations,
it can be easily applied to other CH-NS systems. We apply our model to
simulate multiple vesicles interacting with each other and the domain
boundaries under flow conditions in silico, such as red blood cells passing a
branched blood vessel. elaborate

Rest of the paper is organized as follows. In Section 2, the thermodynamically
consistent model considering cell-cell and cell-wall interaction is derived,
the energy dissipation law of the model is given. Then the discrete scheme of
the model and the discrete energy law that is consistent with the continuous
condition is given in Section 3. Section 4 is used to present results of
numerical simulations and compare them with the data collected in laboratory
experiments. The conclusion is drawn in Section 5.

##  2  Model Derivation

Derivation of the model in this paper is based on the energy dissipation law
which holds ubiquitously in physical systems involving irreversible processes
[46, 47, 48, 49, 44]. This law states that for an isothermal and closed system
the rate of change of the energy of the system is equal to the dissipation of
the energy as follows:

$\frac{d}{dt}E^{\rm total}=-\triangle\leq 0~{},$ (1)

where $E^{\rm total}$ is the total energy of the system, which is the sum of
the kinetic energy $\mathcal{K}$ and the Helmholtz free energy $\mathcal{F}$
of the system, and $\triangle$ is the rate of energy dissipation, which in
fact is entropy production. Eq. (1) can be easily derived via the combination
of the First and Second Laws of Thermodynamics. The choices of the total
energy functional and the dissipation functional, together with the kinematic
(transport) relations of the variables employed in the system, determine all
the physical and mechanical considerations and assumptions for the problem
[50]. A general technique to determine these relations of variables may be
found in [38, 44].

###  2.1  Multi-vesicle interaction System

Energy variational method [37] is adopted for model derivation to ensure the
thermodynamics consistency of the derived model. Specifically, the model
respects the energy dissipation law given by (1). The model derivation begins
with defining functionals of the total energy and dissipation of the system,
respectively, and making the kinematic assumptions of the state variables of
the system based on physics laws of conservation. We refer the readers to [37,
39, 51] for detailed discussions of this method.

Figure 1 is a schematic showing the setup of the model. Let the problem domain
be $\Omega$, and its wall boundary be $\partial\Omega_{w}$. For a multi-cell
flow system, the dynamic evolution of the $i^{th}$ cell (or vesicle) under
flow conditions within $\Omega$ is tracked by the phase-field function
$\phi_{i}(\bm{x},t)$. Notice that $\phi_{i}(\bm{x},t)\in[-1,1]$ with
$\phi_{i}=1$ for intracellular space and $\phi_{i}=-1$ for extracellular
space. The membrane of the cell is identified by $\phi_{i}(\bm{x},t)=0$. We
also introduce a phase-field function $\phi_{w}(\bm{x})\in[-1,1]$ to represent
the wall boundary of the domain as shown in Fig. 1. This is for considering
cell-wall interaction described below.

Figure 1:  A schematic showing the model setup. Phase-field functions are used
to represent the domain wall boundary and the cell interface.

We assume that the dynamical evolution of the phase-field function $\phi_{i}$
is a generalized gradient flow, see Eq. (2a). We also utilize the laws of
conservation for describing the dynamics of the momentum and the total mass of
the system, see Eqs. (2b)-(2c). The cell membrane (interface) is assumed to be
inextensible. The equation for local inextensibility of the interface is given
by Eq. (2d).

$\displaystyle\frac{\partial\phi_{i}}{\partial
t}+\nabla\cdot(\mathbf{u}\phi_{i})=\nabla\cdot{\mathbf{q}}_{\phi_{i}}~{},$
(2a) $\displaystyle\rho(\frac{\partial\mathbf{u}}{\partial
t}+(\mathbf{u}\cdot\nabla)\mathbf{u})=\nabla\cdot\bm{\sigma}_{\eta}+\mathbf{F}_{\phi_{1},\phi_{2},...,\phi_{N}}~{},$
(2b) $\displaystyle\nabla\cdot\mathbf{u}=0~{},$ (2c)
$\displaystyle(\mathcal{P}_{i}:{\nabla}\mathbf{u})\delta_{i}=0~{},$ (2d)

where $\rho$ is the averaged density of the system. In this work, $\rho$ is
assumed to be a constant. $\mathbf{u}$ is the macroscopic velocity of the
system, and $\bm{\sigma}_{\eta}$ is the system’s visco-elastic stress to be
determined.

The undetermined flux $\mathbf{q}_{\phi_{i}}$ in Eq. (2a) will be determined
by postulating that $\phi_{i}$ is driven by gradients in a chemical potential.
This leads to the Cahn-Hilliard equation which ensures conservation of the
volume of a cell during its dynamic evolution.
$\bm{F}_{\phi_{1},\phi_{2},...,\phi_{n}}$ in Eq. (2b) is the body force
induced by vesicle-fluid interaction and is to be determined as well.

Equation (2d) is the diffuse interface approximation of the local
inextensibility of the membrane of the $i^{th}$ cell with relaxation [39, 43].
$\delta_{i}=\frac{1}{2}\gamma^{2}|\nabla\phi_{i}|^{2}$ is the surface delta
function with the diffuse interface thickness $\gamma$. $\mathcal{P}_{i}$ is
the projection operator, and is defined to be
$(I-\mathbf{n}_{i,m}\otimes\mathbf{n}_{i,m})$.
$\mathbf{n}_{i,m}=\frac{\nabla\phi_{i}}{|\nabla\phi_{i}|}$ is the unit outward
normal vector of the interface when it is defined as an implicit surface by
the phase-field function. This equation is equivalent to
$\mathcal{P}_{i}:{\nabla}\mathbf{u}=0$ under a sharp interface condition. In
phase field frame, this local inextensibility is extended to the whole domain
by multiplying a scalar function $\delta_{i}$ [37] for the convenience of
computation. Following the idea in [43, 39], here we add a relaxation term
$\xi\gamma^{2}{\nabla}\dot{(}\phi_{i}^{2}{\nabla}\lambda_{i})$ in Eq. (2d)
where $\xi$ is a parameter independent of $\gamma$, $\lambda_{i}$ is a
function that measures the interface “pressure” induced by the inextensibility
of the membrane of the $i^{th}$ cell. Thus Eq. (2d) now takes the form:

$\displaystyle(\mathcal{P}_{i}:{\nabla}\mathbf{u})\delta_{i}+\xi\gamma^{2}{\nabla}\dot{(}\phi_{i}^{2}{\nabla}\lambda_{i})=0$
(3)

Eq.(2a)-(2c) and (3) together constitute the governing equations of the
system.

The boundary conditions on the top and bottom of the domain, denoted by
$\partial\Omega_{w}$, are given as follows:

$\displaystyle\left\\{\begin{array}[]{l}\mathbf{u}\cdot\mathbf{n}=0~{},\\\
\mathbf{u}_{\tau}\cdot\bm{\tau}_{k}=f_{\bm{\tau}_{k}}~{},\\\
\frac{\partial\phi_{i}}{\partial
t}+\mathbf{u}\cdot\nabla_{\Gamma}\phi_{i}=J_{\Gamma_{i}}~{},\\\
\mathbf{q}_{\phi_{i}}\cdot\mathbf{n}=0~{}.\end{array}\right.$ (8)

On the boundary $\partial\Omega_{w}$, an Allen-Cahn
$(\ref{assumption_bd})_{3}$ type boundary condition is employed for
$\phi_{i}$. $\nabla_{\Gamma}=\nabla-\mathbf{n}(\mathbf{n}\cdot\nabla)$ is the
surface gradient operator, and
$\mathbf{u}_{\tau}=\mathbf{u}-(\mathbf{u}\cdot\mathbf{n})\mathbf{n}$ is the
fluid slip velocity with respect to the wall where $\bm{\tau}_{i},i=1,2$ are
the tangential directions of the wall surface. $\mathbf{n}$ is the unit
outward normal vector of the wall. $f_{\bm{\tau}_{i}}$ is the slip velocity of
the fluid on the wall along the $\bm{\tau}_{i}$ direction. $J_{\Gamma}$
represents the Allan-Cahn type of relaxation on the wall by using the phase-
field method, and is to be determined, see [52, 37].

The total energy functional $E_{total}$ for the multi-cell system is assumed
to be the sum of the kinetic energy $E_{kin}$ in the macroscale, elastic
energy $E_{cell}$ of cell membrane, cell-cell interaction energy $E_{int}$ and
cell-wall adhesion energy $E_{w}$ in the microscale. Therefore

$\displaystyle E_{total}$ $\displaystyle=$
$\displaystyle\underbrace{E_{kin}}_{Macroscale}+\underbrace{E_{cell}+E_{int}+E_{w}}_{Microscale}~{}$
$\displaystyle=$
$\displaystyle\int_{\Omega}\left(\frac{1}{2}\rho|\mathbf{u}|^{2}\right)d\bm{x}+\sum_{i=1}^{N}\int_{\Omega}\frac{\hat{\kappa}_{B}}{2\gamma}\left|\frac{f(\phi_{i})}{\gamma}\right|^{2}d\bm{x}+\sum_{i=1}^{N}\frac{\mathcal{M}_{s}}{2}\frac{\left(S(\phi_{i}(t))-S(\phi_{i}(t=0))\right)^{2}}{S(\phi_{i}(t=0))}$
$\displaystyle+$
$\displaystyle\int_{\Omega}Hd\bm{x}+\sum_{i=1}^{N}\int_{\Omega}f_{w}(\phi_{i})ds~{}.$
(10)

Here $\hat{\kappa}_{B}$ is the bending modulus of the membrane, and $\gamma$
is the thickness of the diffuse interface. The membrane energy density is
given by

$\displaystyle f(\phi_{i})$ $\displaystyle=$ $\displaystyle\frac{\delta
G}{\delta\phi_{i}}=-\gamma^{2}\Delta\phi_{i}+(\phi_{i}^{2}-1)\phi_{i}~{},$
(11)

with

$\displaystyle G(\phi_{i})$ $\displaystyle=$
$\displaystyle\frac{\gamma^{2}|\nabla\phi_{i}|^{2}}{2}+\frac{(1-\phi_{i}^{2})^{2}}{4}~{}.$
(12)

The function $S(\phi_{i})=\int_{\Omega}\frac{G(\phi_{i})}{\gamma}dx$ is used
to measure the surface area of the cell [43, 22, 37]. $\mathcal{M}_{s}$ is the
penalty constant.

The term $H$ denotes the interaction energy density induced by the interaction
of cells. There are many different ways to define $H$ [28, 29, 24]. Here we
begin with considering mechanical interaction between two cells identified by
phase-field functions $\phi_{1}$ and $\phi_{2}$, respectively. Recall that
$\phi_{i}=1$ represents the intracellular space, and $\phi_{i}=-1$ represents
the extracellular space of the $i^{\rm th}$ cell, respectively. Whether there
exists mechanical interaction between the two cells can be determined by
measuring the overlapping (i.e., occupying the same physical space) of the
membranes of these two cells. The following two functions are introduced to
measure the overlapping

$d_{1}(\bm{x})=(\phi_{1}(\bm{x})-1)(\phi_{2}(\bm{x})-1)~{},~{}~{}\bm{x}\in\Omega~{},$

and

$d_{2}(\bm{x})=(\phi_{1}(\bm{x})+1)(\phi_{2}(\bm{x})+1)~{},~{}~{}\bm{x}\in\Omega~{}.$

It is obvious that $d_{1},d_{2}\in[0,4]$. From physical point of view, we can
regard $d_{1}$ as the “distance” in the space of phase-field functions since
it reaches maximum when both $\phi_{1}$ and $\phi_{2}$ equal to -1 which means
there is no overlap between the two cells. $d_{1}$ reaches minimum when two
cells fully overlap, i.e., $\phi_{1}=\phi_{2}=1$. Similarly, $d_{2}$ could be
regarded as the extent of overlapping since it reaches maximum when cells
fully overlap, and reaches minimum when there is no overlap. Then we propose
the interaction potential $H$ by a polynomial function of $d_{2}$ and $d_{1}$:

$\displaystyle
H=Q(d_{2}^{\alpha}+C)d_{1}^{\beta}=QC(\phi_{1}+1)^{\beta}(\phi_{2}+1)^{\beta}-Q(\phi_{1}^{2}-1)^{\alpha}(\phi_{2}-1)^{\alpha}(\phi_{1}+1)^{\beta-\alpha}(\phi_{2}+1)^{\beta-\alpha}~{},$
(13)

where $Q$ and $C$ are constant that controls the minimum of the interaction
energy. The first term of Eq. (13) accounts for the repulsion between the two
cells which increase sharply to prevent fully overlapping of the cells, and
the second terms is for the attraction between cells which only appears when
the cells contact each other, which means the they start to overlap.

In this work, we take $\alpha=\beta=2$. Thus the interaction potential (13)
becomes:

$\displaystyle
H=Q_{1}(\phi_{1}+1)^{2}(\phi_{2}+1)^{2}-Q_{2}(\phi_{1}^{2}-1)^{2}(\phi_{2}^{2}-1)^{2}~{}.$
(14)

Now we may consider $H$ as a 2D function of $\phi_{1}$ and $\phi_{2}$, where
both $\phi_{1}$ and $\phi_{2}$ take values between -1 and 1 and their values
represent relative position of vesicle 1 and vesicle 2 when they overlap. Also
we need to point out that it is because $H$ is regarded as a 2D function, we
need two variables ($d_{1},d_{2}$) or ($\phi_{1}$, $\phi_{2}$) to formulate a
single-concave potential energy with only one minima in the domain. Figure 2
plots the energy landscape of the interaction potential due to presence of the
phases $\phi_{1}$ and $\phi_{2}$. The energy is equal to 0 when two phases do
not touch or overlap i.e., $\phi_{1}=\phi_{2}=-1$. When they start to overlap,
the energy firstly decreases which means that the attraction force between
these two phases dominates. Then the energy goes up which indicates the
repulsive force dominates. This prevents the two phases from occupying the
same physical space. So the interacting potential energy behaves conceptually
similar to a 2D kind of Lennard-Jones potential. We should point out that Eq.
(13) or (14) is just constructed to mimic Lennard-Jones type of repulsive-
attractive feature. Other formula with similar behavior may be good too. We
remark that in [25] a potential of form $(\phi_{1}^{2}-1)^{2}(\phi_{2}^{2}-1)$
is used, where only attractive feature is included.

Figure 2: Interaction energy with respect to $\phi_{1}$ and $\phi_{2}$ with
$Q_{1}=50$, $Q_{2}=400$.

More general, the interaction energy in multiple phases case goes:

$\displaystyle
H=\sum_{i=1}^{N}\sum_{i<j}H_{ij}=\sum_{i=1}^{N}\sum_{i<j}\left[Q_{1}(\phi_{i}+1)^{2}(\phi_{j}+1)^{2}-Q_{2}(\phi_{i}^{2}-1)^{2}(\phi_{j}^{2}-1)^{2}\right].$
(15)

For the cell-wall adhesion energy, $\phi_{w}$ is defined along the wall with a
thickness [25] as shown in Figure 4. Following above interaction potential
definition, then the cell-wall interaction energy is defined by setting
$\phi_{2}=\phi_{w}$ in Eq. (14),

$\displaystyle
f_{w}(\phi_{i})=Q_{w_{1}}(\phi_{i}+1)^{2}(\phi_{w}+1)^{2}-Q_{w_{2}}(\phi_{i}^{2}-1)^{2}(\phi_{w}^{2}-1)^{2}~{},$
(16)

where $Q_{w_{1}}$ is repulsive energy density and $Q_{w_{2}}$ is adhesion
energy density. Notice that the cell-wall energy $E_{w}$ is defined in the
bulk not on the boundary. However, $f_{w}$ is non-zero only when the two
phases overlap which could be regarded as the attraction force is only induced
when the vesicle is contacting (or close enough to) the wall.

Then the chemical potential $\mu_{i}$ for each phase $\phi_{i}$ is defined as
follows

$\displaystyle\mu_{i}=\frac{\delta
E_{cell_{i}}}{\delta\phi_{i}}=\frac{\hat{\kappa}_{B}}{\gamma^{3}}g(\phi_{i})+\frac{\partial
H}{\partial\phi_{i}}+\frac{\mathcal{M}_{s}}{\gamma}\frac{S(\phi_{i})-S(\phi_{i,0})}{S(\phi_{i,0})}f(\phi_{i})+\frac{\partial
f_{w}(\phi_{i})}{\partial\phi_{i}}~{},$ (17)

where $g(\phi_{i})=-\gamma^{2}\Delta
f(\phi_{i})+(3\phi_{i}^{2}-1)f(\phi_{i})$.

The dissipation functional of the system consists of the dissipation
introduced by fluid viscosity, friction near the wall, and interfacial mixing
due to the diffuse interface representation [39]:

$\displaystyle\Delta$ $\displaystyle=$
$\displaystyle\int_{\Omega}2\eta|\bm{D}_{\eta}|^{2}d\bm{x}+\sum_{i=1}^{N}\int_{\Omega}\frac{1}{M_{\phi}}|q_{\phi_{i}}|^{2}d\bm{x}+\int_{\partial\Omega_{w}}\beta_{s}|\mathbf{u}_{\bm{\tau}}|^{2}ds$
(18)
$\displaystyle+\sum_{i=1}^{N}\left(\int_{\partial\Omega_{w}}\kappa_{\Gamma}|J_{\Gamma_{i}}|^{2}ds+\int_{\partial\Omega_{w}}\xi|\gamma\phi_{i}\nabla\lambda_{i}|^{2}ds\right)~{},$

where $\bm{D}_{\eta}=\frac{\nabla\mathbf{u}+(\nabla\mathbf{u})^{T}}{2}$,
$\eta$ is the viscosity of fluid mixture, $\beta_{s}$ is wall fraction,
$M_{\phi}$ and $\kappa_{\Gamma}$ are the mobility of phase in the bulk and
boundary. We note that in general, the viscosity $\eta$ could be a function of
all phases $\phi_{i}$

$\displaystyle\eta=\eta(\phi_{w},\phi_{1},\cdots,\phi_{n})~{}.$ (19)

The specific forms of the flux and stress functions in the kinematic equations
(2)-(8) are obtained by taking the time derivative of the total energy
functional and comparing with the defined dissipation functional. Then the
time derivative of the total energy goes:

$\displaystyle\frac{dE_{total}}{dt}$ $\displaystyle=$
$\displaystyle\frac{d}{dt}E_{kin}+\frac{d}{dt}E_{cell}+\frac{d}{dt}E_{w}$ (20)
$\displaystyle\equiv$ $\displaystyle I_{1}+I_{2}+I_{3}~{}.$ (21)

Taking the time derivative of $E_{kin}$, together with the conservation law of
momentum $(\ref{assumption})_{2}$, incompressibility of the fluid
$(\ref{assumption})_{3}$ and inextensibility of the membrane
$(\ref{assumption})_{4}$ yields

$\displaystyle I_{1}$ $\displaystyle=$
$\displaystyle\frac{d}{dt}\int_{\Omega}\frac{\rho|\mathbf{u}|^{2}}{2}d\bm{x}$
(22) $\displaystyle=$
$\displaystyle\int_{\Omega}\frac{1}{2}\frac{\partial\rho}{\partial
t}|\mathbf{u}|^{2}d\bm{x}+\int_{\Omega}\rho\frac{\partial\mathbf{u}}{\partial
t}\cdot\mathbf{u}d\bm{x}$ $\displaystyle=$
$\displaystyle\int_{\Omega}\frac{1}{2}\frac{\partial\rho}{\partial
t}|\mathbf{u}|^{2}d\bm{x}+\int_{\Omega}\rho\frac{d\mathbf{u}}{dt}\cdot\mathbf{u}d\bm{x}-\int_{\Omega}\left(\rho\mathbf{u}\cdot\nabla\mathbf{u}\right)\cdot\mathbf{u}d\bm{x}$
$\displaystyle=$
$\displaystyle\int_{\Omega}\frac{1}{2}\frac{\partial\rho}{\partial
t}|\mathbf{u}|^{2}d\bm{x}+\int_{\Omega}\rho\frac{d\mathbf{u}}{dt}\cdot\mathbf{u}d\bm{x}+\int_{\Omega}\nabla\cdot(\rho\mathbf{u})\frac{|\mathbf{u}|^{2}}{2}d\bm{x}$
$\displaystyle=$
$\displaystyle\int_{\Omega}(\nabla\cdot\bm{\sigma}_{\eta})\cdot\mathbf{u}d\bm{x}+\int_{\Omega}\mathbf{F}\cdot\mathbf{u}d\bm{x}+\sum_{i=1}^{N}\int_{\Omega}\lambda_{i}\delta_{i}\mathcal{P}_{i}:{\nabla}\mathbf{u}d\bm{x}+\sum_{i=1}^{N}\int_{\Omega}\xi\gamma^{2}\lambda_{i}{\nabla}\cdot(\phi_{i}^{2}{\nabla}\lambda_{i})d\bm{x}-\int_{\Omega}pI:\nabla\mathbf{u}d\bm{x}$
$\displaystyle=$
$\displaystyle-\int_{\Omega}((\bm{\sigma}_{\eta}+pI):\nabla\mathbf{u})d\bm{x}+\int_{\Omega}\mathbf{F}\cdot\mathbf{u}dx-\sum_{i=1}^{N}\int_{\Omega}{\nabla}\cdot(\lambda_{i}\delta_{i}\mathcal{P}_{i})\cdot\mathbf{u}d\bm{x}$
$\displaystyle-\sum_{i=1}^{N}\int_{\Omega}\xi\gamma^{2}\phi_{i}^{2}({\nabla}\lambda_{i})^{2}d\bm{x}+\int_{\partial\Omega_{w}}((\bm{\sigma}_{\eta}+\sum_{i=1}^{N}\lambda_{i}\delta_{i}\mathcal{P}_{i})\cdot\mathbf{n})\cdot\mathbf{u}_{\tau}dS~{},$

where the slip boundary condition is used. Here the pressure $p$ is a
Lagrangian multiplier and is introduced to ensure the incompressibility of the
fluid.

Taking the time derivative of $E_{cell}$ together with the conservation of
each phase $(\ref{assumption})_{1}$ yields

$\displaystyle I_{2}$ $\displaystyle=$
$\displaystyle\frac{d}{dt}\int_{\Omega}\frac{\hat{\kappa}_{B}}{2\gamma}\sum_{i=1}^{N}\left|\frac{f(\phi_{i})}{\gamma}\right|^{2}d\bm{x}+\frac{d}{dt}\int_{\Omega}Hd\bm{x}+\sum_{i=1}^{N}\frac{d}{dt}\frac{\mathcal{M}_{s}}{2}\frac{(S(\phi_{i})-S(\phi_{i,0}))^{2}}{S(\phi_{i,0})}$
(23) $\displaystyle=$
$\displaystyle\int_{\Omega}\frac{\hat{\kappa}_{B}}{\gamma}\sum_{i=1}^{N}\frac{f_{i}}{\gamma^{2}}\frac{\partial
f_{i}}{\partial t}d\bm{x}+\int_{\Omega}\sum_{i=1}^{N}\frac{\partial
H}{\partial\phi_{i}}\frac{\partial\phi_{i}}{\partial
t}d\bm{x}+\sum_{i=1}^{N}\frac{d}{dt}\frac{\mathcal{M}_{s}}{2}\frac{(S(\phi_{i})-S(\phi_{i,0}))^{2}}{S(\phi_{i,0})}$
$\displaystyle=$
$\displaystyle\int_{\Omega}\frac{\hat{\kappa}_{B}}{\gamma}\sum_{i=1}^{N}\frac{f_{i}}{\gamma^{2}}\left(-\gamma^{2}\Delta\left(\frac{\partial\phi_{i}}{\partial
t}\right)+(3\phi_{i}^{2}-1)\frac{\partial\phi_{i}}{\partial
t}\right)d\bm{x}+\int_{\Omega}\sum_{i=1}^{N}\frac{\partial
H}{\partial\phi_{i}}\frac{\partial\phi_{i}}{\partial t}d\bm{x}$
$\displaystyle+\sum_{i=1}^{N}\frac{d}{dt}\frac{\mathcal{M}_{s}}{2}\frac{(S(\phi_{i})-S(\phi_{i,0}))^{2}}{S(\phi_{i,0})}$
$\displaystyle=$
$\displaystyle\int_{\Omega}\frac{\hat{\kappa}_{B}}{\gamma^{3}}\sum_{i=1}^{N}\left(-\gamma^{2}\Delta
f_{i}+(3\phi_{i}^{2}-1)f_{i}+\frac{\partial
H}{\partial\phi_{i}}+\frac{\mathcal{M}_{s}}{\gamma}\frac{S(\phi_{i})-S(\phi_{i,0})}{S(\phi_{i,0})}f(\phi_{i})\right)\frac{\partial\phi_{i}}{\partial
t}d\bm{x}$
$\displaystyle-\int_{\partial\Omega_{w}}\frac{\hat{\kappa}_{B}}{\gamma}\sum_{i=1}^{N}f_{i}\frac{\partial}{\partial
t}(\partial_{n}\phi_{i})ds+\int_{\partial\Omega_{w}}\frac{\hat{\kappa}_{B}}{\gamma}\sum_{i=1}^{N}\partial_{n}f_{i}\frac{\partial\phi_{i}}{\partial
t}ds+\sum_{i=1}^{N}\int_{\partial\Omega_{w}}\mathcal{M}_{s}\frac{S(\phi_{i})-S(\phi_{i,0})}{S(\phi_{i,0})}\gamma\partial_{n}\phi_{i}\frac{\partial\phi_{i}}{\partial
t}ds$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\left(\int_{\Omega}\mu_{i}\frac{\partial\phi_{i}}{\partial
t}d\bm{x}+\int_{\partial\Omega_{w}}\frac{\hat{\kappa}_{B}}{\gamma}f_{i}\frac{\partial}{\partial
t}(\partial_{n}\phi_{i})ds+\int_{\partial\Omega_{w}}\mathcal{M}_{s}\frac{S(\phi_{i})-S(\phi_{i,0})}{S(\phi_{i,0})}\gamma\partial_{n}\phi_{i}\frac{\partial\phi_{i}}{\partial
t}ds\right)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\left(\int_{\Omega}\mu_{i}\nabla\cdot
q_{\phi_{i}}d\bm{x}-\int_{\Omega}\mu_{i}\mathbf{u}\cdot\nabla\phi_{i}d\bm{x}+\int_{\partial\Omega_{w}}\frac{\hat{\kappa}_{B}}{\gamma}\partial_{n}f_{i}\frac{\partial\phi_{i}}{\partial
t}ds+\int_{\partial\Omega_{w}}\mathcal{M}_{s}\frac{S(\phi_{i})-S(\phi_{i,0})}{S(\phi_{i,0})}\gamma\partial_{n}\phi_{i}\frac{\partial\phi_{i}}{\partial
t}ds\right)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\left(-\int_{\Omega}q_{\phi_{i}}\cdot\nabla\mu_{i}d\bm{x}-\int_{\Omega}\mu_{i}\mathbf{u}\cdot\nabla\phi_{i}d\bm{x}+\int_{\partial\Omega_{w}}\frac{\hat{\kappa}_{B}}{\gamma}\partial_{n}f_{i}\frac{\partial\phi_{i}}{\partial
t}ds\right.$
$\displaystyle\left.+\int_{\partial\Omega_{w}}\mathcal{M}_{s}\frac{S(\phi_{i})-S(\phi_{i,0})}{S(\phi_{i,0})}\gamma\partial_{n}\phi_{i}\frac{\partial\phi_{i}}{\partial
t}ds\right)~{},$

where the Allan-Cahn boundary condition $(\ref{assumption_bd})_{3}$ for each
phase is used.

Computing $\frac{d}{dt}E_{w}$ yields

$\displaystyle I_{3}=\sum_{i=1}^{N}\left(\int_{\partial\Omega}\frac{\partial
f_{w}(\phi_{i})}{\partial\phi_{i}}\frac{\partial\phi_{i}}{\partial
t}ds\right)~{}.$ (24)

By combining Eqs. (22), (23) and (24), we have:

$\displaystyle\frac{d}{dt}E_{total}$ $\displaystyle=$
$\displaystyle-\int_{\Omega}((\bm{\sigma}_{\eta}+pI):\nabla\mathbf{u})d\bm{x}+\int_{\Omega}(\mathbf{F}-\sum_{i=1}^{N}\mu_{i}\nabla\phi_{i}-\sum_{i=1}^{N}{\nabla}\cdot(\lambda_{i}\delta_{i}\mathcal{P}_{i}))\cdot\mathbf{u}d\bm{x}$
(25)
$\displaystyle-\sum_{i=1}^{N}\int_{\Omega}q_{\phi_{i}}\cdot\nabla\mu_{i}d\bm{x}+\sum_{i=1}^{N}\int_{\Omega}\xi\gamma^{2}\phi_{i}^{2}({\nabla}\lambda_{i})^{2}d\bm{x}+\int_{\partial\Omega_{w}}((\bm{\sigma}_{\eta}+\sum_{i=1}^{N}\lambda_{i}\delta_{i}\mathcal{P}_{i})\cdot\mathbf{n})\cdot\mathbf{u}_{\tau}ds+\sum_{i=1}^{N}\int_{\partial\Omega_{s}}\hat{L}_{i}\frac{\partial\phi_{i}}{\partial
t}ds$ $\displaystyle=$
$\displaystyle-\int_{\Omega}((\bm{\sigma}_{\eta}+pI):\nabla\mathbf{u})d\bm{x}+\int_{\Omega}(\mathbf{F}-\sum_{i=1}^{N}\mu_{i}\nabla\phi_{i}-\sum_{i=1}^{N}{\nabla}\cdot(\lambda_{i}\delta_{i}\mathcal{P}_{i}))\cdot\mathbf{u}d\bm{x}$
$\displaystyle-\sum_{i=1}^{N}\int_{\Omega}q_{\phi_{i}}\cdot\nabla\mu_{i}d\bm{x}+\sum_{i=1}^{N}\int_{\Omega}\xi\gamma^{2}\phi_{i}^{2}({\nabla}\lambda_{i})^{2}d\bm{x}+\int_{\partial\Omega_{w}}((\bm{\sigma}_{\eta}+\sum_{i=1}^{N}\lambda_{i}\delta_{i}\mathcal{P}_{i})\cdot\mathbf{n})\cdot\mathbf{u}_{\tau}ds$
$\displaystyle+\sum_{i=1}^{N}\int_{\partial\Omega_{s}}\hat{L}_{i}(-\mathbf{u}\cdot\nabla_{\Gamma}\phi_{i}+J_{\Gamma_{i}})ds$
$\displaystyle=$
$\displaystyle-\int_{\Omega}((\bm{\sigma}_{\eta}+pI):\nabla\mathbf{u})d\bm{x}+\int_{\Omega}(\mathbf{F}-\sum_{i=1}^{N}\mu_{i}\nabla\phi_{i}-\sum_{i=1}^{N}{\nabla}\cdot(\lambda_{i}\delta_{i}\mathcal{P}_{i}))\cdot\mathbf{u}d\bm{x}$
$\displaystyle-\sum_{i=1}^{N}\int_{\Omega}q_{\phi_{i}}\cdot\nabla\mu_{i}d\bm{x}+\sum_{i=1}^{N}\int_{\Omega}\xi\gamma^{2}\phi_{i}^{2}({\nabla}\lambda_{i})^{2}d\bm{x}+\int_{\partial\Omega_{w}}((\bm{\sigma}_{\eta}+\sum_{i=1}^{N}\lambda_{i}\delta_{i}\mathcal{P}_{i})\cdot\mathbf{n}-\sum_{i=1}^{N}\hat{L}_{i}\nabla_{\Gamma}\phi_{i})\cdot\mathbf{u}_{\tau}ds$
$\displaystyle+\sum_{i=1}^{N}\int_{\partial\Omega_{w}}\hat{L}_{i}J_{\Gamma_{i}}ds~{},$

where
$\hat{L}_{i}=\displaystyle\frac{\hat{\kappa}_{B}}{\gamma}\partial_{n}f_{i}+\mathcal{M}_{s}\frac{S(\phi_{i})-S(\phi_{i,0})}{S(\phi_{i,0})}\gamma\partial_{n}\phi_{i}$.

To close the system, we use the energy dissipation law [49, 37, 53]

$\displaystyle\frac{dE}{dt}=-\Delta~{},$ (26)

which states that the rate of changing total energy of the system is induced
by the dissipation. Comparing Eq. (25) with predefined dissipation functional
in Eq. (18) gives rise to

$\displaystyle\left\\{\begin{array}[]{ll}\bm{\sigma}_{\eta}=2\eta\bm{D}_{\eta}-pI~{},&\mbox{in~{}}\Omega~{},\\\
q_{\phi_{j}}=M_{\phi_{j}}\nabla\mu_{j}~{},&\mbox{in~{}}\Omega~{},\\\
\mathbf{F}=\sum_{i=1}^{N}(\mu_{i}\nabla\phi_{i}+\nabla\cdot(\lambda_{i}\delta_{i}\mathcal{P}_{i}))~{},&\mbox{in~{}}\Omega~{},\\\
J_{\Gamma_{i}}=-\kappa_{\Gamma_{i}}^{-1}\hat{L}_{i}~{},&\mbox{on~{}}\partial\Omega_{w}~{},\\\
u_{\tau_{j}}=\beta_{s}^{-1}(-(\mathbf{n}\cdot(\bm{\sigma}_{\eta}+\sum_{i=1}^{N}\lambda_{i}\delta_{i}\mathcal{P}_{i})\cdot\bm{\tau}_{j})+\sum_{i=1}^{N}\hat{L}_{i}\partial_{\tau_{j}}\phi_{i})~{},~{}j=1,2,&\mbox{on~{}}\partial\Omega_{w}~{}.\end{array}\right.$
(32)

To this end, the proposed multi-cellular interaction model is composed of the
following equations:

$\displaystyle\left\\{\begin{array}[]{lll}\frac{\partial\phi_{i}}{\partial
t}+\nabla\cdot(\mathbf{u}\phi_{i})=M_{\phi_{i}}\Delta\mu_{i}~{},\\\
\mu_{i}=\frac{\hat{\kappa}_{B}}{\gamma^{3}}g(\phi_{i})+\frac{\partial
H}{\partial\phi_{i}}+\frac{M_{s}}{\gamma}\frac{S(\phi_{i})-S(\phi_{i,0})}{S(\phi_{i,0})}f(\phi_{i})+\frac{\partial
f_{w}(\phi_{i})}{\partial\phi_{i}}~{},\\\ g(\phi_{i})=-\gamma^{2}\Delta
f_{i}+(3\phi^{2}-1)f(\phi_{i}),\\\
f(\phi_{i})=-\gamma^{2}\Delta\phi_{i}+(\phi_{i}^{2}-1)\phi_{i}~{},\\\
\rho(\frac{\partial\mathbf{u}}{\partial
t}+(\mathbf{u}\cdot\nabla)\mathbf{u})+\nabla
p=\nabla\cdot(2\eta\bm{D}_{\eta})+\sum_{j=1}^{i}(\mu_{j}\nabla\phi_{j}+\nabla\cdot(\lambda_{j}\delta_{j}\mathcal{P}_{j}))~{},\\\
\nabla\cdot\mathbf{u}=0~{},\\\
\delta_{i}(\mathcal{P}_{i}:{\nabla}\mathbf{u})+\xi\gamma^{2}\nabla\cdot(\phi_{i}^{2}\nabla\lambda_{i})=0~{},\end{array}\right.$
(40)

with the boundary conditions on $\partial\Omega_{w}$

$\displaystyle\left\\{\begin{array}[]{l}\mathbf{u}\cdot\mathbf{n}=0~{},\\\
-\beta_{s}u_{\tau_{j}}=(\mathbf{n}\cdot(\bm{\sigma}_{\eta}+\sum_{i=1}^{N}\lambda_{i}\delta_{i}\mathcal{P}_{i})\cdot\bm{\tau}_{j})-\sum_{i=1}^{N}\hat{L}_{i}\partial_{\tau_{j}}\phi_{i}~{},~{}j=1,2,\\\
f_{i}=0~{},\\\ \kappa_{\Gamma_{i}}\left(\frac{\partial\phi_{i}}{\partial
t}+\mathbf{u}\cdot\nabla_{\Gamma}\phi_{i}\right)=-\hat{L}_{i}~{},\\\
\hat{L}_{i}=\frac{\hat{\kappa}_{B}}{\gamma}\partial_{n}f_{i}+\mathcal{M}_{s}\frac{s(\phi_{i})-S(\phi_{i,0})}{S(\phi_{i,0})}\gamma\partial_{n}\phi_{i}~{},\\\
\partial_{n}\mu_{i}=0~{}.\end{array}\right.$ (47)

###  2.2 Dimensionless governing equations

The viscosity, length, velocity, time, bulk and boundary chemical potentials
in the equations are scaled by their corresponding characteristic values
$\eta_{0}$, $L$, $U$, $\frac{L}{U}$, $\frac{\eta_{0}U}{L}$ and $\eta_{0}U$,
respectively. Write $Q_{w_{1}},Q_{w_{2}},Q_{1},Q_{2}$ into
$Q_{0}q_{w_{1}},Q_{0}q_{w_{2}}.Q_{0}q_{1},Q_{0}q_{2}$, where $Q_{0}$ is the
character energy density. The governing equation of the system can be
rewritten as

$\displaystyle\left\\{\begin{array}[]{ll}Re(\frac{\partial\mathbf{u}}{\partial
t}+(\mathbf{u}\cdot{\nabla})\mathbf{u})+{\nabla}P={\nabla}\cdot(2\eta\bm{D})+\sum_{i}\mu_{i}{\nabla}\phi_{i}+\sum_{i}{\nabla}\cdot(\lambda_{i}\delta_{\epsilon_{i}}\mathcal{P}_{i})~{},&\mbox{in~{}}\Omega~{},\\\\[8.53581pt]
{\nabla}\cdot\mathbf{u}=0~{},&\mbox{in~{}}\Omega~{},\\\\[8.53581pt]
\frac{\partial\phi_{i}}{\partial
t}+\mathbf{u}\cdot{\nabla}\phi_{i}=-\mathcal{M}\Delta\mu_{i}~{},&\mbox{in~{}}\Omega~{},\\\\[8.53581pt]
\mu_{i}=\kappa_{B}g(\phi_{i})+\mathcal{M}_{s}\frac{(S(\phi_{i})-S(\phi_{i,0}))}{S(\phi_{i,0})}f_{i}+\alpha\frac{\partial
H}{\partial\phi_{i}}+\alpha\frac{\partial
f_{w}(\phi_{i})}{\partial\phi_{i}}~{},&\mbox{in~{}}\Omega~{},\\\\[8.53581pt]
f_{i}=-\epsilon\Delta\phi_{i}+\frac{(\phi_{i}^{2}-1)}{\epsilon}\phi_{i},~{}g(\phi_{i})=-\Delta
f_{i}+\frac{1}{\epsilon^{2}}(3\phi_{i}^{2}-1)f_{i}~{},&\mbox{in~{}}\Omega~{},\\\\[8.53581pt]
\delta_{\epsilon_{i}}(\mathcal{P}_{i}:{\nabla}\mathbf{u})+\xi\varepsilon^{2}\nabla\cdot(\phi_{i}^{2}\nabla\lambda_{i})=0~{},&\mbox{in~{}}\Omega~{},\end{array}\right.$
(54)

with the boundary conditions

$\displaystyle\displaystyle\left\\{\begin{array}[]{ll}\kappa\dot{\phi_{i}}+L(\phi_{i})=0~{},&\mbox{on~{}}\partial\Omega_{w}~{},\\\
L(\phi_{i})=\kappa_{B}\partial_{n}f(\phi_{i})+\epsilon\mathcal{M}_{s}\frac{S(\phi_{i})-S(\phi_{i,0})}{S(\phi_{i,0})}\partial_{n}\phi_{i}~{},&\mbox{on~{}}\partial\Omega_{w}~{},\\\
-l_{s}^{-1}u_{\tau_{i}}=\bm{\tau_{i}}\cdot(2\eta\mathbf{D}_{\eta}+\sum_{i}\lambda_{i}\delta_{\epsilon_{i}}\mathcal{P}_{i})\cdot\bm{n}-\sum_{i}L(\phi_{i})\partial_{\tau_{i}}\phi_{i}~{},~{}i=1,2,&\mbox{on~{}}\partial\Omega_{w}~{},\\\
f_{i}=0~{},&\mbox{on~{}}\partial\Omega_{w}~{},\\\
\partial_{n}\mu_{i}=0~{},&\mbox{on~{}}\partial\Omega_{w}~{},\end{array}\right.$
(60)

where
$S(\phi_{i})=\displaystyle\int_{\Omega}\frac{\epsilon}{2}|\nabla\phi_{i}|^{2}+\frac{1}{4\epsilon}(\phi_{i}^{2}-1)^{2}d\bm{x}$
and $\delta_{\epsilon_{i}}=\frac{1}{2}\epsilon^{2}|{\nabla}\phi_{i}|^{2}$.

The dimensionless constants appeared in Eqs. (54)-(60) are given by
$\epsilon=\frac{\gamma}{L}$, $Re=\frac{\rho_{0}UL}{\eta_{0}}$,
$\mathcal{M}=\frac{M_{\phi}\eta_{0}}{L^{2}}$,
$\kappa_{B}=\frac{\hat{\kappa}_{B}}{L^{2}\eta_{0}U}$,
$k=\frac{\hat{\kappa}_{B}}{\eta_{0}L}$, $l_{s}=\frac{\eta_{0}}{\beta_{s}L}$,
$\alpha=\frac{Q_{0}}{\eta_{0}U}$, $\mathcal{M}_{s}=\frac{M_{s}}{\eta_{0}U}$.

If we define the Sobolev spaces as follows [54, 37]

$\displaystyle\bm{W}^{1,3}=(W^{1,3})^{2}~{},$ (61)
$\displaystyle\bm{W}_{N}^{1,\frac{3}{2}}=(W^{1,\frac{3}{2}})^{N}~{},$ (62)
$\displaystyle\bm{W}_{N}^{1,3}=(W^{1,3})^{N}~{},$ (63)
$\displaystyle\bm{W}_{N}^{1,\frac{3}{2}}(\Omega)=\left\\{\Lambda=(\lambda_{1},\lambda_{2},...\lambda_{N})^{T}\right\\}~{},$
(64)
$\displaystyle\bm{W}^{1,3}(\Omega)=\left\\{\mathbf{u}=(u_{x},u_{y})^{T}\in\bm{W}^{1,3}|\mathbf{u}\cdot\bm{n}=0,~{}\mbox{on~{}}\partial\Omega_{w}\right\\}~{},$
(65)
$\displaystyle\bm{W}_{\Phi}^{1,3}(\Omega)=\left\\{\Phi=(\phi_{1},\phi_{2},...\phi_{N})^{T}\in\bm{W}_{N}^{1,3}|-1\leq\phi_{i}\leq
1,i=1,2,...N,~{}\mbox{in~{}}\Omega\right\\}~{},$ (66)
$\displaystyle\bm{W}_{\bm{U}}^{1,3}(\Omega)=\left\\{\bm{U}=(\mu_{1},\mu_{2},...\mu_{N})^{T}\in\bm{W}_{N}^{1,3}|\partial_{n}\mu_{i}=0,i=1,2,...N,~{}\mbox{on~{}}\partial\Omega_{w}\right\\}~{},$
(67)
$\displaystyle\bm{W}_{F}^{1,3}(\Omega)=\left\\{F=(f_{1},f_{2},...f_{N})^{T}\in\bm{W}_{N}^{1,3}|f_{i}=0,i=1,2,...N,~{}\mbox{on~{}}\partial\Omega_{w}\right\\}~{},$
(68)
$\displaystyle\bm{W}_{b}=\bm{W}_{\Phi}^{1,3}(\Omega)\times\bm{W}_{F}^{1,3}(\Omega)\times\bm{W}_{\bm{U}}^{1,3}(\Omega)\times\bm{W}_{N}^{1,\frac{3}{2}}(\Omega)\times
W^{1,\frac{3}{2}}(\Omega)\times\bm{W}^{1,3}(\Omega)~{},$ (69)

and let $\|\cdot\|=\left(\int_{\Omega}|\cdot|^{2}d\bm{x}\right)^{\frac{1}{2}}$
and
$\|\cdot\|_{w}=\left(\int_{\partial\Omega_{w}}|\cdot|^{2}ds\right)^{\frac{1}{2}}$
denote the $L^{2}$ norm defined in the domain and on the domain boundary
respectively, then the system (54)-(60) satisfies the following energy law.

######  Theorem 2.1.

If $(\Phi,~{}F,~{}\bm{U},~{}\lambda,~{}P,~{}\mathbf{u})\in\bm{W}_{b}$ are
smooth solutions of the above system (54)-(60), then the following energy law
is satisfied:

$\displaystyle\frac{d}{dt}\mathcal{E}_{total}=\frac{d}{dt}(\mathcal{E}_{kin}+\mathcal{E}_{cell}+\mathcal{E}_{w})$
(70) $\displaystyle=$
$\displaystyle\frac{1}{Re}\left(-2\|\eta^{1/2}\mathbf{D}_{\eta}\|^{2}-\mathcal{M}\sum_{i}\|{\nabla}\mu_{i}\|^{2}-\xi\sum_{i}\|\epsilon\phi_{i}\nabla\lambda_{i}\|^{2}-\kappa\sum_{i}\|\dot{\phi_{i}}\|_{w}^{2}-\|l_{s}^{-1/2}\mathbf{u}_{\tau}\|^{2}_{w}\right)~{},$

where
$\mathcal{E}_{total}=\mathcal{E}_{kin}+\mathcal{E}_{cell}+\mathcal{E}_{w}$,
$\mathcal{E}_{kin}=\displaystyle\frac{1}{2}\int_{\Omega}|\mathbf{u}|^{2}d\bm{x}$,
$\mathcal{E}_{cell}=\displaystyle\frac{\kappa_{B}}{2Re\epsilon}\sum_{i}\int_{\Omega}|f_{i}|^{2}d\bm{x}+\mathcal{M}_{s}\sum_{i}\frac{(S(\phi_{i})-S(\phi_{i,0}))^{2}}{2ReS(\phi_{i,0})}+\frac{\alpha}{Re}\int_{\Omega}Hd\bm{x}$
and
$\mathcal{E}_{w}\\!\\!=\\!\\!\frac{\alpha}{Re}\\!\\!\sum_{i}\displaystyle\int_{\Omega}f_{w}(\phi_{i})d\bm{x}$.

Proof: Multiplying the first equation in Eq. (54) with $\mathbf{u}$ and
integration by parts yield

$\displaystyle\frac{d}{dt}\mathcal{E}_{kin}$ $\displaystyle=$
$\displaystyle\frac{1}{Re}\left\\{-\int_{\Omega}2\eta|\bm{D}_{\eta}|^{2}d\bm{x}+\int_{\partial\Omega_{w}}(\bm{\sigma}_{\eta}\cdot\mathbf{n})\cdot\mathbf{u}_{\tau}ds+\sum_{i}\int_{\Omega}\mu_{i}\nabla\phi_{i}\cdot\mathbf{u}d\bm{x}-\sum_{i}\int_{\Omega}\lambda_{i}\delta_{\epsilon_{i}}\mathcal{P}_{i}:\nabla\mathbf{u}d\bm{x}\right.$
(71)
$\displaystyle\left.+\sum_{i}\int_{\partial\Omega_{w}}(\lambda_{i}\delta_{\epsilon_{i}}\mathcal{P}_{i}\cdot\mathbf{n})\cdot\mathbf{u}_{\tau}ds\right\\}$
$\displaystyle=$
$\displaystyle\frac{1}{Re}\left\\{-\int_{\Omega}2\eta|\bm{D}_{\eta}|^{2}d\bm{x}-\sum_{i}\int_{\Omega}\lambda_{i}\delta_{\epsilon_{i}}\mathcal{P}_{i}:\nabla\mathbf{u}d\bm{x}-l_{s}^{-1}\int_{\partial\Omega_{w}}|\mathbf{u}_{\tau}|^{2}ds\right.$
$\displaystyle\left.+\sum_{i}\int_{\partial\Omega_{w}}L(\phi_{i})\partial_{\tau}\phi\cdot\mathbf{u}_{\tau}ds+\sum_{i}\int_{\Omega}\mu_{i}\nabla\phi_{i}\cdot\mathbf{u}d\bm{x}~{},\right\\}$

where the slip boundary condition in Eq. (60) is applied.

Taking the inner product of the third equation in Eq. (54) with
$\frac{\mu_{i}}{Re}$ and summing up with respect to $i$ result in

$\displaystyle\frac{1}{Re}\sum_{i}\int_{\Omega}\frac{\partial\phi_{i}}{\partial
t}\mu_{i}d\bm{x}+\frac{1}{Re}\sum_{i}\int_{\Omega}\mathbf{u}\cdot\nabla\phi_{i}\mu_{i}d\bm{x}=-\frac{1}{Re}\mathcal{M}\sum_{i}\int_{\Omega}|\nabla\mu_{i}|^{2}d\bm{x}~{},$
(72)

where $\partial_{n}\mu_{i}=0$ is considered here.

Multiplying the fourth equation in Eq. (54) with
$\frac{1}{Re}\frac{\partial\phi_{i}}{\partial t}$ and integration by parts
give rise to

$\displaystyle\quad\frac{1}{Re}\sum_{i}\int_{\Omega}\mu\frac{\partial\phi_{i}}{\partial
t}d\bm{x}$ $\displaystyle=$
$\displaystyle\frac{1}{Re}\sum_{i}\left\\{\kappa_{B}\int_{\Omega}g_{i}\frac{\partial\phi_{i}}{\partial
t}d\bm{x}+\mathcal{M}_{s}\frac{S(\phi_{i})-S(\phi_{i,0})}{S(\phi_{i,0})}\int_{\Omega}f_{i}\frac{\partial\phi_{i}}{\partial
t}d\bm{x}+\alpha\int_{\Omega}\frac{\partial
H}{\partial\phi_{i}}\frac{\partial\phi_{i}}{\partial
t}d\bm{x}+\alpha\int_{\Omega}\frac{f_{w}(\phi_{i})}{\partial\phi_{i}}\frac{\partial\phi_{i}}{\partial
t}d\bm{x}\right\\}$ $\displaystyle=$
$\displaystyle\frac{\kappa_{B}}{Re}\sum_{i}\int_{\Omega}f_{i}\frac{\partial}{\partial
t}\left(-\Delta\phi_{i}+\frac{1}{\epsilon^{2}}(\phi_{i}^{3}-\phi_{i})\right)d\bm{x}-\frac{\kappa_{B}}{Re}\sum_{i}\int_{\partial\Omega_{w}}\partial_{n}f_{i}\frac{\partial\phi_{i}}{\partial
t}ds+\mathcal{M}_{s}\sum_{i}\frac{d}{dt}\left(\frac{(S(\phi_{i})-S(\phi_{i,0}))^{2}}{2ReS(\phi_{i,0})}\right)$
$\displaystyle-\mathcal{M}_{s}\sum_{i}\left(\frac{S(\phi_{i})-S(\phi_{i,0})}{ReS(\phi_{i,0})}\right)\int_{\partial\Omega_{w}}\epsilon\partial_{n}\phi_{i}\frac{\partial\phi_{i}}{\partial
t}ds+\frac{\alpha}{Re}\sum_{i}\int_{\Omega}\frac{\partial
f_{w}(\phi_{i})}{\partial\phi_{i}}\frac{\partial\phi_{i}}{\partial
t}d\bm{x}+\frac{\alpha}{Re}\sum_{i}\int_{\Omega}\frac{\partial
H}{\partial\phi_{i}}\frac{\partial\phi_{i}}{\partial t}d\bm{x}$
$\displaystyle=$
$\displaystyle\frac{d}{dt}\left(\kappa_{B}\sum_{i}\int_{\Omega}\frac{|f_{i}|^{2}}{2Re\epsilon}d\bm{x}\right)+\mathcal{M}_{s}\frac{d}{dt}\left(\sum_{i}\frac{(S(\phi_{i})-S(\phi_{i,0}))^{2}}{2ReS(\phi_{i,0})}\right)+\frac{\alpha}{Re}\frac{d}{dt}\sum_{i}\int_{\Omega}f_{w}(\phi_{i})d\bm{x}+\frac{\alpha}{Re}\frac{d}{dt}\int_{\Omega}Hd\bm{x}$
$\displaystyle-\sum_{i}\int_{\partial\Omega_{w}}\frac{L(\phi_{i})}{Re}\frac{\partial\phi_{i}}{\partial
t}ds$ $\displaystyle=$
$\displaystyle\frac{d}{dt}(\mathcal{E}_{cell}+\mathcal{E}_{w})-\int_{\partial\Omega_{w}}\frac{L(\phi)}{Re}\frac{\partial\phi}{\partial
t}ds~{},$

where the definitions of $f(\phi)$, $g(\phi)$ and the boundary conditions of
$\phi$ and $f$ are utilized.

Multiplying the last equations with $\frac{\lambda_{i}}{Re}$ and integration
by parts and sum up by $i$ leads to

$\displaystyle\frac{1}{Re}\sum_{i}\int_{\Omega}(\lambda_{i}\delta_{\epsilon_{i}}\mathcal{P}_{i}):\nabla\mathbf{u}d\bm{x}-\frac{1}{Re}\sum_{i}\int_{\Omega}\xi\varepsilon^{2}\phi_{i}^{2}({\nabla}\lambda_{i})^{2}d\bm{x}=0~{}.$
(74)

Finally, the energy dissipation law (70) is obtained by combining Eqs. (71),
(72), (2.2) and (74) considering the boundary conditions in (60).
$\hbox{}\hfill\blacksquare$

##  3 Numerical Scheme and Discrete Energy Law

In this section, an energy stable temporal discretization scheme is first
proposed for the multi-cellular system (54)-(60). Then a $C^{0}$ finite
element method is used for spacial discretization to obtain the fully discrete
scheme.

###  3.1 Time-discrete primitive method

The mid-point method is used for the temporal discretization of Eqs.
(54)-(60). Let $\Delta t$ denote the time step size, $()^{n+1}$ and $()^{n}$
denote the values of the variables at times $(n+1)\Delta t$ and $n\Delta t$,
respectively. The semi-discrete in time scheme to solve Eq. (54) is as
follows:

$\displaystyle\left\\{\begin{array}[]{ll}\frac{\mathbf{u}^{n+1}-\mathbf{u}^{n}}{\Delta
t}+(\mathbf{u}^{n+\frac{1}{2}}\cdot{\nabla})\mathbf{u}^{n+\frac{1}{2}}+\frac{1}{Re}{\nabla}P^{n+\frac{1}{2}}=\frac{1}{Re}{\nabla}\cdot(\eta^{n}({\nabla}\mathbf{u}^{n+\frac{1}{2}}+({\nabla}\mathbf{u}^{n+\frac{1}{2}})^{T}))\\\
~{}~{}~{}~{}~{}~{}~{}~{}+\frac{1}{Re}\sum_{i}\mu_{i}^{n+\frac{1}{2}}{\nabla}\phi_{i}^{n+\frac{1}{2}}+\sum_{i}\frac{1}{Re}{\nabla}\cdot\left(\lambda_{i}^{n+\frac{1}{2}}\mathcal{P}_{i}^{n}\delta_{\epsilon_{i}}\right)~{},\\\
{\nabla}\cdot\mathbf{u}^{n+\frac{1}{2}}=0~{},\\\
\frac{\phi_{i}^{n+1}-\phi_{i}^{n}}{\Delta
t}+(\mathbf{u}^{n+\frac{1}{2}}\cdot{\nabla})\phi_{i}^{n+\frac{1}{2}}=-\mathcal{M}\Delta\mu_{i}^{n+\frac{1}{2}}~{},\\\
\mu_{i}^{n+\frac{1}{2}}=\kappa_{B}g(\phi_{i}^{n+1},\phi_{i}^{n})+\mathcal{M}_{s}\frac{(S(\phi_{i}^{n+\frac{1}{2}})-S(\phi_{0_{i}}))}{S(\phi_{0_{i}})}f(\phi_{i}^{n+1},\phi_{i}^{n})\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}+\alpha\frac{H_{i}^{n+1}-H_{i}^{n}}{\phi_{i}^{n+1}-\phi_{i}^{n}}+\alpha\frac{f_{w}(\phi_{i}^{n+1})-f_{w}(\phi_{i}^{n})}{\phi_{i}^{n+1}-\phi_{i}^{n}}~{},\\\
f_{i}^{n+\frac{1}{2}}=-\epsilon\Delta\phi_{i}^{n+\frac{1}{2}}+\frac{1}{\epsilon}((\phi_{i}^{n+\frac{1}{2}})^{2}-1)\phi_{i}^{n+\frac{1}{2}}~{},\\\
\delta_{\epsilon_{i}}\mathcal{P}_{i}^{n}:{\nabla}\mathbf{u}^{n+\frac{1}{2}}+\xi\epsilon^{2}\nabla\cdot((\phi_{i}^{n})^{2}\nabla\lambda_{i}^{n+\frac{1}{2}})=0~{},\end{array}\right.$
(83)

with boundary conditions on $\partial\Omega_{w}$ ,

$\displaystyle\left\\{\begin{array}[]{ll}\kappa\dot{\phi}_{i}^{n+\frac{1}{2}}=-L_{i}^{n+\frac{1}{2}}~{},\\\
L_{i}^{n+\frac{1}{2}}=\kappa_{B}\partial_{n}f_{i}^{n+\frac{1}{2}}+\mathcal{M}_{s}\epsilon\frac{S(\phi_{i}^{n+\frac{1}{2}})-S_{0_{i}}}{S_{0_{i}}}\partial_{n}\phi_{i}^{n+\frac{1}{2}}~{},\\\
-l_{s}^{-1}u_{\tau_{j}}^{n+\frac{1}{2}}=\bm{\tau_{j}}\cdot(\eta^{n}({\nabla}\mathbf{u}^{n+\frac{1}{2}}+({\nabla}\mathbf{u}^{n+\frac{1}{2}})^{T})+\sum_{i}\lambda_{i}^{n+\frac{1}{2}}\delta_{\epsilon_{i}}\mathcal{P}_{i}^{n})\cdot\bm{n}\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\sum_{i}L_{i}^{n+\frac{1}{2}}\partial_{\tau_{j}}\phi_{i}^{n+\frac{1}{2}}~{},~{}~{}j=1,2,\\\
f_{i}^{n+\frac{1}{2}}=0~{},\\\
\partial_{n}\lambda_{i}^{n+\frac{1}{2}}=0~{},\end{array}\right.$ (90)

with $(\cdot)^{n+\frac{1}{2}}=\frac{(\cdot)^{n}+(\cdot)^{n+1}}{2}$ and
$\mathcal{P}_{i}^{n}=I-\mathbf{n}^{n}_{m}\otimes\mathbf{n}^{n}_{m}$ with
$\mathbf{n}_{m}^{n}=\frac{\nabla\phi_{i}^{n}}{|\nabla\phi_{i}^{n}|}$ and

$\displaystyle\left\\{\begin{array}[]{l}f(\phi_{i}^{n+1},\phi_{i}^{n})=-\epsilon\Delta\phi_{i}^{n+\frac{1}{2}}+\frac{1}{4\epsilon}((\phi_{i}^{n+1})^{2}+(\phi_{i}^{n})^{2}-2)(\phi_{i}^{n+1}+\phi_{i}^{n})~{},\\\
g(\phi_{i}^{n+1},\phi_{i}^{n})=\left(-\Delta
f_{i}^{n+\frac{1}{2}}+\frac{1}{\epsilon^{2}}\left((\phi_{i}^{n+1})^{2}+(\phi_{i}^{n})^{2}+\phi_{i}^{n+1}\phi_{i}^{n}-1\right)f_{i}^{n+\frac{1}{2}}\right)~{},\\\
H_{i}^{n}=q_{1}(\phi_{i}^{n}+1)^{2}\sum_{j\neq
i}\left[(\phi_{j}^{n+\frac{1}{2}}+1)^{2}\right]-q_{2}((\phi_{i}^{n})^{2}-1)^{2}\sum_{j\neq
i}\left[((\phi_{j}^{n+\frac{1}{2}})^{2}-1)^{2}\right]~{},\\\
H_{i}^{n+1}=q_{1}(\phi_{i}^{n+1}+1)^{2}\sum_{j\neq
i}\left[(\phi_{j}^{n+\frac{1}{2}}+1)^{2}\right]-q_{2}((\phi_{i}^{n+1})^{2}-1)^{2}\sum_{j\neq
i}\left[((\phi_{j}^{n+\frac{1}{2}})^{2}-1)^{2}\right]~{},\\\
f_{w}(\phi^{n})=q_{w_{1}}(\phi_{i}^{n}+1)^{2}(\phi_{w}+1)^{2}-q_{w_{2}}((\phi_{i}^{n})^{2}-1)((\phi_{w})^{2}-1)~{}.\end{array}\right.$
(96)

Thus we have

$\displaystyle\frac{H_{i}^{n+1}-H_{i}^{n}}{\phi_{i}^{n+1}-\phi_{i}^{n}}$
$\displaystyle=$
$\displaystyle\frac{1}{\phi_{i}^{n+1}-\phi_{i}^{n}}\left({q_{1}}(\phi_{i}^{n+1}+1)^{2}\sum_{j\neq
i}\left[(\phi_{j}^{n+\frac{1}{2}}+1)^{2}\right]-q_{2}((\phi_{i}^{n+1})^{2}-1)^{2}\sum_{j\neq
i}\left[((\phi_{j}^{n+\frac{1}{2}})^{2}-1)^{2}\right]\right.$ (97)
$\displaystyle-\left.q_{1}(\phi_{i}^{n}+1)^{2}\sum_{j\neq
i}\left[(\phi_{j}^{n+\frac{1}{2}}+1)^{2}\right]+q_{2}((\phi_{i}^{n})^{2}-1)^{2}\sum_{j\neq
i}\left[(\phi_{j}^{n+\frac{1}{2}})^{2}-1\right]\right)$ $\displaystyle=$
$\displaystyle q_{1}(\phi_{i}^{n+1}+\phi_{i}^{n}+2)\sum_{j\neq
i}\left[(\phi_{j}^{n+\frac{1}{2}}+1)^{2}\right]$ $\displaystyle-
q_{2}(\phi_{i}^{n+1}+\phi_{i}^{n})((\phi_{i}^{n+1})^{2}+(\phi_{i}^{n})^{2}-2)\sum_{j\neq
i}\left[((\phi_{j}^{n+\frac{1}{2}})^{2}-1)^{2}\right]$

Similarly,

$\displaystyle\frac{f_{w}(\phi_{i}^{n+1})-f_{w}(\phi_{i}^{n})}{\phi_{i}^{n+1}-\phi_{i}^{n}}$
$\displaystyle=$ $\displaystyle
q_{w_{1}}(\phi_{i}^{n+1}+\phi_{i}^{n}+2)(\phi_{w}+1)^{2}$ (98) $\displaystyle-
q_{w_{2}}(\phi_{i}^{n+1}+\phi_{i}^{n})((\phi_{i}^{n+1})^{2}+(\phi_{i}^{n})^{2}-2)(\phi_{w}^{2}-1)^{2}$

Later on, we keep the form
$\frac{f_{w}(\phi_{i}^{n+1})-f_{w}(\phi_{i}^{n})}{\phi_{i}^{n+1}-\phi_{i}^{n}}$
and $\frac{H(\phi_{i}^{n+1})-H(\phi_{i}^{n})}{\phi_{i}^{n+1}-\phi_{i}^{n}}$for
convenience in later derivation.

The above scheme obeys the following theorem of energy stability.

######  Theorem 3.1.

If $(\phi_{i}^{n},\mu_{i}^{n},\mathbf{u}^{n},P^{n})$ are smooth solutions of
the above system (83)-(90), then the following energy law is satisfied:

$\displaystyle\mathcal{E}_{total}^{n+1}-\mathcal{E}_{total}^{n}$
$\displaystyle\\!\\!\\!=\\!\\!\\!$
$\displaystyle(\mathcal{E}_{kin}^{n+1}+\sum_{i}^{N}\left[\mathcal{E}_{cell_{i}}^{n+1}+\mathcal{E}_{i,int}^{n+1}+\mathcal{E}_{w_{i}}^{n+1}\right])-(\mathcal{E}_{kin}^{n}+\sum_{i}^{N}\left[\mathcal{E}_{cell_{i}}^{n}+\mathcal{E}_{i,int}^{n}+\mathcal{E}_{w_{i}}^{n}\right])$
(99) $\displaystyle=$ $\displaystyle\frac{\triangle
t}{Re}\left(-2\|(\eta^{n})^{1/2}\mathbf{D}^{n+\frac{1}{2}}_{\eta}\|^{2}-\mathcal{M}\sum_{i}^{N}\|{\nabla}\mu_{i}^{n+\frac{1}{2}}\|^{2}-\xi\sum_{i}^{N}\|\
\epsilon\phi_{i}^{n}\nabla\lambda_{i}^{n+\frac{1}{2}}\|^{2}\right.$
$\displaystyle\left.-\frac{1}{\kappa}\|\sum_{i}^{N}L(\phi_{i}^{n+\frac{1}{2}})\|_{w}^{2}-\|l_{s}^{-1/2}\mathbf{u}_{\tau}^{n+\frac{1}{2}}\|^{2}_{w}\right)~{},$

where
$\mathcal{E}_{total}^{n}=\mathcal{E}_{kin}^{n}+\sum_{i}^{N}\left[\mathcal{E}_{cell_{i}}^{n}+\mathcal{E}_{i,int}^{n}+\mathcal{E}_{w_{i}}^{n}\right]$
with $\mathcal{E}_{kin}^{n}=\frac{1}{2}\|\mathbf{u}^{n}\|^{2}$,
$\mathcal{E}_{cell_{i}}^{n}=\frac{\kappa_{B}\|f_{i}^{n}\|^{2}}{2Re\epsilon}+\mathcal{M}_{s}\frac{(S(\phi_{i}^{n})-S(\phi_{i,0}))^{2}}{2ReS(\phi_{i,0})}+\frac{\alpha}{Re}H_{i}^{n}$
and
$\mathcal{E}_{i,w}^{n}\\!\\!=\\!\\!\frac{\alpha}{Re}\\!\\!\displaystyle\int_{\Omega}f_{i,w}^{n}d\bm{x}$.

The following two lemmas are needed for proving Theorem 3.1. Proof of these
two lemmas can be found in [39].

######  Lemma 3.2.

Let

$\displaystyle
f(\phi^{n+1},\phi^{n})=-\epsilon\Delta\phi^{n+\frac{1}{2}}+\frac{1}{4\epsilon}((\phi^{n+1})^{2}+(\phi^{n})^{2}-2)(\phi^{n+1}+\phi^{n})~{}.$
(100)

Then $f(\phi^{n+1},\phi^{n})$ satisfies

$\displaystyle\int_{\Omega}f(\phi^{n+1},\phi^{n})(\phi^{n+1}-\phi^{n})d\bm{x}=S^{n+1}-S^{n}-\int_{\partial\Omega_{w}}\epsilon\partial_{n}\phi^{n+\frac{1}{2}}(\phi^{n+1}-\phi^{n})ds~{},$
(101)

where
$S^{n+1}=\int_{\Omega}G(\phi^{n+1})d\bm{x},S^{n}=\int_{\Omega}G(\phi^{n})d\bm{x}$.

######  Lemma 3.3.

Let $g(\phi^{n+1},\phi^{n})=-\Delta
f^{n+\frac{1}{2}}+\frac{1}{\epsilon^{2}}((\phi^{n+1})^{2}+(\phi^{n})^{2}+\phi^{n+1}\phi^{n}-1)f^{n+\frac{1}{2}}~{}.$
Then $g(\phi^{n+1},\phi^{n})$ satisfies

$\displaystyle\int_{\Omega}g(\phi^{n+1},\phi^{n})(\phi^{n+1}-\phi^{n})d\bm{x}$
(102) $\displaystyle=$
$\displaystyle\int_{\Omega}\frac{1}{2\epsilon}((f^{n+1})^{2}-(f^{n})^{2})d\bm{x}-\int_{\partial\Omega_{w}}\partial_{n}f^{n+\frac{1}{2}}(\phi^{n+1}-\phi^{n})ds~{},$

where
$f^{n+1}=-\epsilon\Delta\phi^{n+1}+\frac{1}{\epsilon}((\phi^{n+1})^{2}-1)\phi^{n+1},f^{n}=-\epsilon\Delta\phi^{n}+\frac{1}{\epsilon}((\phi^{n})^{2}-1)\phi^{n}$.

Proof of Theorem 3.1: Multiplying the first equation in system (83) by $\Delta
t\mathbf{u}^{n+\frac{1}{2}}$ gives

$\displaystyle\int_{\Omega}\frac{1}{2}((\mathbf{u}^{n+1})^{2}-(\mathbf{u}^{n})^{2})d\bm{x}+\int_{\Omega}\Delta
t\mathbf{u}^{n+\frac{1}{2}}\cdot((\mathbf{u}^{n+\frac{1}{2}}{\nabla})\cdot\mathbf{u}^{n+\frac{1}{2}})d\bm{x}$
(103) $\displaystyle-\frac{\Delta
t}{Re}\int_{\Omega}P^{n+\frac{1}{2}}{\nabla}\cdot\mathbf{u}^{n+\frac{1}{2}}d\bm{x}$
$\displaystyle=$ $\displaystyle-\frac{\Delta
t}{Re}\int_{\Omega}{\nabla}\mathbf{u}^{n+\frac{1}{2}}:\eta^{n}({\nabla}\mathbf{u}^{n+\frac{1}{2}}+({\nabla}\mathbf{u}^{n+\frac{1}{2}})^{T})d\bm{x}+\frac{\Delta
t}{Re}\sum_{i}\int_{\Omega}\mathbf{u}^{n+\frac{1}{2}}\cdot{\nabla}\phi_{i}^{n+1}\mu_{i}^{n+1}d\bm{x}$
$\displaystyle-\frac{\Delta
t}{Re}\sum_{i}\int_{\Omega}\lambda_{i}\delta_{\epsilon_{i}}\mathcal{P}_{i}^{n}:\nabla\mathbf{u}^{n+\frac{1}{2}}d\bm{x}+\frac{\Delta
t}{Re}\sum_{i}\int_{\partial\Omega_{w}}\lambda_{i}^{n+\frac{1}{2}}(\delta_{\epsilon_{i}}\mathcal{P}_{i}^{n}\cdot\mathbf{n})\cdot\mathbf{u}_{\tau}^{n+\frac{1}{2}}ds$
$\displaystyle+\frac{\Delta
t}{Re}\int_{\partial\Omega_{w}}\mathbf{u}^{n+\frac{1}{2}}\cdot\eta^{n}(({\nabla}\mathbf{u}^{n+\frac{1}{2}}+({\nabla}\mathbf{u}^{n+\frac{1}{2}})^{T})\cdot\bm{n})ds~{}.$

Multiplying the fourth equation in system (83) by
$\frac{\phi_{i}^{n+1}-\phi_{i}^{n}}{Re}$ and integration by parts lead to

$\displaystyle\frac{1}{Re}\sum_{i}\int_{\Omega}\mu_{i}^{n+1/2}(\phi_{i}^{n+1}-\phi_{i}^{n})d\bm{x}=\frac{\kappa_{B}}{Re}\sum_{i}\int_{\Omega}\frac{1}{2\epsilon}((f_{i}^{n+1})^{2}-(f_{i}^{n})^{2})d\bm{x}$
$\displaystyle+\frac{\mathcal{M}_{s}}{Re}\frac{(S(\phi_{i}^{n+1})-S_{i,0})^{2}-(S(\phi_{i}^{n})-S_{i,0})^{2}}{2S_{i,0}}+\frac{\alpha}{Re}\sum_{i}\int_{\Omega}(H_{i}^{n+1}-H_{i}^{n})d\bm{x}+\frac{\alpha}{Re}\sum_{i}\int_{\Omega}(f_{w}(\phi_{i}^{n+1})-f_{w}(\phi_{i}^{n}))d\bm{x}$
(104)
$\displaystyle-\frac{\kappa_{B}}{Re}\sum_{i}\int_{\partial\Omega_{w}}\partial_{n}f_{i}^{n+\frac{1}{2}}(\phi_{i}^{n+1}-\phi_{i}^{n})ds-\frac{\mathcal{M}_{s}}{Re}\sum_{i}\int_{\partial\Omega_{w}}\frac{S(\phi_{i}^{n+\frac{1}{2}})-S_{i,0}}{S_{i,0}}\epsilon\partial_{n}\phi_{i}^{n+\frac{1}{2}}(\phi_{i}^{n+1}-\phi_{i}^{n})ds~{}.$

Multiplying the third equation in system (83) by
$\frac{\mu_{i}^{n+\frac{1}{2}}\Delta t}{Re}$ yield

$\displaystyle\frac{1}{Re}\sum_{i}\int_{\Omega}\mu_{i}^{n+\frac{1}{2}}(\phi_{i}^{n+1}-\phi_{i}^{n})d\bm{x}+\frac{\Delta
t}{Re}\sum_{i}\int_{\Omega}\mu_{i}^{n+\frac{1}{2}}(\mathbf{u}^{n+\frac{1}{2}}\cdot{\nabla})\phi_{i}^{n+\frac{1}{2}}d\bm{x}$
(105) $\displaystyle=$ $\displaystyle-\frac{\mathcal{M}\Delta
t}{Re}\sum_{i}\int_{\Omega}(\nabla\mu_{i}^{n+\frac{1}{2}})^{2}d\bm{x}~{}.$

Multiplying the last equation in system (83) by
$\frac{\lambda^{n+\frac{1}{2}}\Delta t}{Re}$ and integration by parts then sum
by $i$ give

$\displaystyle\frac{\Delta
t}{Re}\sum_{i}\int_{\Omega}(\lambda_{i}^{n+\frac{1}{2}}\delta_{\epsilon_{i}}\mathcal{P}_{i}^{n}):\nabla\mathbf{u}^{n+\frac{1}{2}}d\bm{x}-\frac{\Delta
t}{Re}\sum_{i}\int_{\Omega}\xi\epsilon^{2}(\phi_{i}^{n})^{2}\left|\nabla\lambda_{i}^{n+\frac{1}{2}}\right|^{2}d\bm{x}=0~{}.$
(106)

The discretized energy dissipation law (99) is obtained by combining Eqs.
(103)-(106) and organizing the terms according to the boundary conditions
$L(\phi_{i})$ as shown in (90). $\hbox{}\hfill\blacksquare$

###  3.2 Fully discrete $C^{0}$ finite element scheme

The spatial discretization using $C^{0}$ finite element is straight forward.
Let $\Omega$ be the domain of interest with a Lipschitz-continuous boundary
$\partial\Omega$. Let ${\mathbf{W}_{b}}^{h}\subset\mathbf{W}_{b}$ be a finite
element space with respect to the triangulation of the domain $\Omega$. The
fully discrete scheme of the system is to find

$\left(\\{\Phi_{h}\\}^{n+1},\\{\bm{U}\\}_{h}^{n+1},\\{F_{h}\\}^{n+1},\\{\Lambda_{h}\\}^{n+1},\\{p_{h}\\}^{n+1},\\{\mathbf{u}_{h}\\}^{n+1}\right)\in\mathbf{W}_{b}^{h},$

such that for any
$(\psi_{1,h},...,\psi_{N,h},\chi_{1,h},...,\chi_{N,h},\zeta_{1,h},...,\zeta_{N,h},\Theta_{1,h},...,\Theta_{N,h},q_{h},\mathbf{v}_{h})\in\mathbf{W}_{b}^{h}$,
the following scheme holds.

$\displaystyle\left\\{\begin{array}[]{ll}\int_{\Omega}\left(\frac{\mathbf{u}_{h}^{n+1}-\mathbf{u}_{h}^{n}}{\Delta
t}+(\mathbf{u}_{h}^{n+\frac{1}{2}}\cdot{\nabla})\mathbf{u}_{h}^{n+\frac{1}{2}}+\frac{1}{Re}{\nabla}P_{h}^{n+\frac{1}{2}}\right)\cdot\mathbf{v}_{h}d\bm{x}\\\
=-\int_{\Omega}\frac{1}{Re}(\eta_{h}^{n}({\nabla}\mathbf{u}_{h}^{n+\frac{1}{2}}+({\nabla}\mathbf{u}_{h}^{n+\frac{1}{2}})^{T})):{\nabla}\mathbf{v}_{h}d\bm{x}\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\sum_{i}\int_{\Omega}\frac{1}{Re}\mu_{i,h}^{n+\frac{1}{2}}{\nabla}\phi_{h}^{n+\frac{1}{2}}\cdot\mathbf{v}_{h}d\bm{x}-\sum_{i}\int_{\Omega}\frac{1}{Re}\lambda_{i,h}^{n+\frac{1}{2}}\mathcal{P}_{i,h}^{n}\delta_{i,h,\epsilon}:\mathbf{v}_{h}d\bm{x}\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\int_{\partial\Omega_{w}}\frac{1}{Re}\mathbf{n}\cdot(\eta_{h}^{n}({\nabla}\mathbf{u}_{h}^{n+\frac{1}{2}}+({\nabla}\mathbf{u}_{h}^{n+\frac{1}{2}})^{T})+\sum_{i}\lambda_{i,h}^{n+\frac{1}{2}}\mathcal{P}_{i,h}^{n}\delta_{i,\epsilon})\cdot\mathbf{v}_{h}d\bm{x}~{},\\\
\int_{\Omega}({\nabla}\cdot\mathbf{u}_{h}^{n+\frac{1}{2}})q_{h}d\bm{x}=0~{},\\\\[8.53581pt]
\int_{\Omega}(\frac{\phi_{i,h}^{n+1}-\phi_{i,h}^{n}}{\Delta
t}+(\mathbf{u}_{h}^{n+\frac{1}{2}}\cdot{\nabla})\phi_{i,h}^{n+\frac{1}{2}})\psi_{i,h}d\bm{x}=-\int_{\Omega}\mathcal{M}{\nabla}\mu_{i,h}^{n+\frac{1}{2}}{\nabla}\psi_{i,h}d\bm{x}~{},\\\\[8.53581pt]
\int_{\Omega}\mu_{i,h}^{n+\frac{1}{2}}\chi_{i,h}d\bm{x}=\int_{\Omega}\left(\kappa_{B}\frac{1}{\epsilon^{2}}((\phi_{i,h}^{n+1})^{2}+(\phi_{i,h}^{n})^{2}+\phi_{i,h}^{n+1}\phi_{i,h}^{n}-1)f_{i,h}^{n+\frac{1}{2}}\right.\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\left.+\mathcal{M}_{s}\frac{(S(\phi_{i,h}^{n+\frac{1}{2}})-S(\phi_{i,h,0}))}{S(\phi_{i,h,0})}(\frac{1}{4\epsilon}((\phi_{i,h}^{n+1})^{2}+(\phi_{i,h}^{n})^{2}-2)(\phi_{i,h}^{n+1}+\phi_{i,h}^{n}))\right)\chi_{i,h}d\bm{x}\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\int_{\Omega}(\kappa_{B}{\nabla}f_{i,h}^{n+\frac{1}{2}}+\mathcal{M}_{s}\epsilon\frac{(S(\phi_{i,h}^{n+\frac{1}{2}})-S(\phi_{i,h,0}))}{S(\phi_{i,h,0})}{\nabla}\phi_{i,h}^{n+\frac{1}{2}})\cdot{\nabla}\chi_{i,h}d\bm{x}\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\int_{\Omega}\alpha\frac{f_{w}(\phi_{i,h}^{n+1})-f_{w}(\phi_{i,h}^{n})}{\phi_{i,h}^{n+1}-\phi_{i,h}^{n}}\chi_{i,h}d\bm{x}+\int_{\Omega}\alpha\frac{(H_{i,h}^{n+1}-H_{i,h}^{n})}{\phi_{i,h}^{n+1}-\phi_{i,h}^{n}}\chi_{i,h}d\bm{x}\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\int_{\partial\Omega_{w}}(\kappa_{B}\partial_{\mathbf{n}}f_{i,h}^{n+\frac{1}{2}}+\mathcal{M}_{s}\epsilon\frac{(S(\phi_{i,h}^{n+\frac{1}{2}})-S(\phi_{i,h,0}))}{S(\phi_{i,h})}\partial_{\mathbf{n}}\phi_{i,h}^{n+\frac{1}{2}})\chi_{i,h}ds~{},\\\
\int_{\Omega}f_{i,h}^{n+\frac{1}{2}}\zeta_{i,h}=\int_{\Omega}\epsilon{\nabla}\phi_{i,h}^{n+\frac{1}{2}}\cdot{\nabla}\zeta_{i,h}+\int_{\Omega}\frac{1}{\epsilon}((\phi_{i,h}^{n+\frac{1}{2}})^{2}-1)\phi_{i,h}^{n+\frac{1}{2}}\zeta_{i,h}d\bm{x}\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\int_{\partial\Omega_{w}}\varepsilon\partial_{\mathbf{n}}\phi_{i,h}^{n+\frac{1}{2}}\zeta_{i,h}d\bm{x}~{},\\\
\int_{\Omega}\delta_{i,h,\epsilon}\mathcal{P}_{i,h}^{n}:{\nabla}\mathbf{u}_{i,h}^{n+\frac{1}{2}}\Theta_{i,h}d\bm{x}-\int_{\Omega}\xi\epsilon^{2}((\phi_{i,h}^{n})^{2}\nabla\lambda_{i,h}^{n+\frac{1}{2}})\cdot{\nabla}\Theta_{i,h}d\bm{x}+\int_{\partial\Omega_{w}}\xi\epsilon^{2}((\phi_{i,h}^{n})^{2}\partial_{\mathbf{n}}\lambda_{i,h}^{n+\frac{1}{2}})\Theta_{i,h}d\bm{x}=0~{}.\end{array}\right.$
(121)

Newton’s iteration method is applied to solve the above nonlinear system. The
unique solubility of (121) can be proved by following the idea introduced in
[39].

##  4 Numerical Results

In this section, we first calibrate the model parameters and validate the
model by comparing with the experimental data on RBC deformation under
different stretching forces. Then cell-wall interaction simulations are used
to study the effects of adhesion force and local insensibility. The energy
stability of the numerical scheme is also illustrated. Then the calibrated
model is used to study the cell-wall interaction and cell-cell aggregation.

###  4.1 BenchMark: vesicle deformation under stretch

Laboratory experiments have tested the non-linear elasticity and deformation
of RBC [55] in which optical tweezers are used to provide stretching force to
the cells. We set up a numerical simulation mimicking RBC deformation in the
experiment. New phases $\phi_{tw_{1}},\phi_{tw_{2}}$ are introduced to
represent optical tweezers. The energy corresponding to optical tweezers
stretching RBC is defined as:

$\displaystyle
H_{tw_{i}}=Q_{tw1}(\phi_{1}+1)^{2}(\phi_{tw_{i}}+1)^{2}-Q_{tw2}(\phi_{1}^{2}-1)^{2}(\phi_{tw_{i}}^{2}-1)^{2}~{},~{}~{}i=1,2~{},$
(122)

where $Q_{tw}$ is a coefficient that controls the attraction (???) force that
the optical tweezers provide. The value of other parameters are listed as
follows: $Re=2\times 10^{-4},\mathcal{M}=0.25,\kappa_{B}=2\times
10^{-3},k=2\times 10^{-12},l_{s}=5\times 10^{-3},\mathcal{M}_{s}=2.$ The force
applied on the cell is calculated by the following equation:

$\displaystyle\mathbf{F}=\int_{\Omega}\frac{\partial
H_{tw_{i}}}{\partial\phi}{\nabla}\phi d\bm{x}~{}.$ (123)

The curves of axial and transverse diameter versus stretching forces are shown
in Figure 3 together with the experimental data (diamonds with bars) from
[55].

Figure 3: Nonlinear elastic deformation of red blood cell. The curve shows the
change of the diameter versus stretching force. The diamond is the experiment
data published in [55]. In the experiment schematic diagram, the centered
phase is set to be the cell and the circulate area on the two side stands for
the light tweezers. The force applied on the membrane is kept a constant when
the tweezers moving. The equilibrium is reached when the membrane no longer
extends under the certain stretch force.

The results show that our model fits the experimental data very well.

###  4.2 Cell-wall attraction

Cell-wall interactions under blood flow conditions play important roles during
blood clotting [56] and cancer cell invasion [10]. This example is used to
investigate the effect of adhesion force on the vesicle-wall interaction. We
first consider the case without external flow. As shown in Figure 4, a vesicle
is initially placed at a location with a point-wise contact with the wall
phase. The parameter values of this simulation are listed as follows:
$Re=2\times 10^{-4},\mathcal{M}=5\times 10^{-4},\kappa_{B}=2\times
10^{-2},\varepsilon=2\times 10^{-3},\mathcal{M}_{s}=10^{2},k=4\times
10^{-11},l_{s}=5\times 10^{-6},\alpha=1000,q_{w_{1}}=2,q_{w_{2}}=1$.

Figure 4: The initial state of the case. Figure 5: The top three pictures
shows the deformation of the vesicle with local inextensibility. The bottom
three pictures shows the deformation of the vesicle without local
inextensibility.

Figure 5 shows the scenarios of the the vesicle with and without local
inextensibility at different time. As shown in Fig. 6, the local
inextensiblity constrains slows down the deformation of the vesicle. For
vesicle without local constrain, it allowed the membrane attached to the wall
to be stretched (red color in Fig.6) to achieve the equilibrium faster. This
also could be seen in the Fig. 7. It confirms that in both cases, the energy
monotonously decays to the same value while the vesicle without local
inextensibility evolves faster.

Figure 6: The top two pictures shows the value of surface divergence
$\mathcal{P}_{i}:{\nabla}\mathbf{u}$ with local inextensibility. The bottom
pictures shows the surface divergence of the vesicle without local
inextensibility. Figure 7: Total energy of the two system versus time.

Then the effect of the strength of the adhesion force on the equilibrium
profiles of the vesicles are illustrated in Figure 8. As expected, when the
strength increases from 400 to 2500, the contact surface area increases.

Figure 8: Equilibrium of the cells at different strength of the adhesion
force. (Left: $\alpha=400$; Mid: $\alpha=1000$; Right: $\alpha=2500$.)

During the blood clotting, the platelets interact with several coagulation
factors and get activated. The full activated plated has strong adhesion force
with vessel wall in order to stop bleeding more efficiently. In Figure 9, we
simulate two platelets under the shear blood flow, one fully activated with
strong adhesion and one partial activated with weak adhesion with vessel wall.
It shows that when the fully activated one is captured by the wall; while the
partial activated one is washed away by the shear flow.

Figure 9: The top three picture shows the motion of the cell in strong
adhesion case. The bottom three pictures shows the cell motion in weak
adhesion case.

###  4.3 Cell aggregation and offset

Aggregation of red blood cell is a phenomenon observed in experiment [1]. In
this part, we set up a simulation with four red blood cells contacted with
each other at a small point. The parameters are shown below. $Re=2\times
10^{-5},\mathcal{M}=5\times 10^{-4},\kappa_{B}=4\times 10^{-2},k=4\times
10^{-12},l_{s}=5\times 10^{-3},\mathcal{M}_{s}=10^{3},\alpha=3\times
10^{3},q_{1}=1,q_{2}=0.5$ The evolution of the system with time is shown in
10.

Figure 10: Aggregating of four red blood cells.

From the result we can see that they are creeping together with time under the
attractive force and form rouleaux which is consistent with the experimental
result shown in [1] and [13].

Also, the deformation of the cells at equilibrium is related to the value of
the attractive force [13, 11, 17]. Fig. 11 shows the status of moderately and
strongly aggregate. The parameters are shown below. $Re=2\times
10^{-4},\mathcal{M}=5\times 10^{-4},\kappa_{B}=2\times 10^{-2},k=2\times
10^{-11},l_{s}=5\times 10^{-3},\mathcal{M}_{s}=10.$

Figure 11: Left: Strong aggregate $\alpha=15\times 10^{3},q_{1}=1,q_{2}=1$.
Right: moderate aggregate $\alpha=100,q_{1}=1,q_{2}=1$.

From the figure we can see that under strong aggregation force, obvious
terminal hemispherical caps is shown and the offset between each adjacent cell
is smaller as well. This result is consistent with the experimental phenomenon
observed well [1]. In the following test we set the cells on a Couette flow
with shear rate equal to be $20s^{-1}$ with dimension as in [14]. The motion
of the cells are shown in figure 12. Under shear flow, the rouleaux with
strong adhesion still aggregate together; while the one with weak adhesion is
broken up by the shear flow, which is consistent with the results reported in
[14].

Figure 12: The top three figures shows the motion of the cells with strong
aggregation. The bottom three figures shows the motion under moderate
aggregation.

###  4.4 Red blood cell motion at vessel bifurcation

In this part, we simulate the motion of RBCs in branched vessels. 4 cells are
initially placed in the Y shaped vessel. Width of the main channel of the Y
shaped vessel and the bottom branch of the vessel is set to be $1\times
10^{-7}$ meter. The top branch of the vessel is $0.7\times 10^{-7}$ meter,
which is close to the size of a red blood cell. A pressure drop boundary
condition is used to introduce a shear flow in the vessel with the velocity
around $5\times 10^{-4}m/s$, which is close to the blood flow in capillaries.
Other model parameters in this simulation are as follows: $Re=2\times
10^{-4},\mathcal{M}=5\times 10^{-4},\kappa_{B}=4\times 10^{-2},k=4\times
10^{-11},l_{s}=2=5\times 10^{-6},\mathcal{M}_{s}=20$.

Motion of the RBCs and the velocity field of the flow are shown in Figure 14,
15, and 13, respectively. We firstly simulate motion of the RBC group with a
moderate aggregation force when they pass a Y-shaped channel with the same
width of the both branches as a baseline. The cells divide equally at the
vessel bifurcation. Then one of the channel is widen and a non-equal deviation
is observed (3 move upwards,1 moves downward) due to the lower resistance in
the upper branch. Then, under the same geometric setup, a strong aggregation
force is applied to the cells. In this case, all of the cells goes into the
wider channel. The simulation result explains the experimental result in [1],
which a great number of red blood cells are observed to be absent in the
branched vessel under strong aggregate case compared with moderate aggregate.

Figure 13: Cells are set in a cluster initially under a moderate aggregating
force with $\alpha=25,q_{1}=1,q_{2}=1$. The main channel width is 1 and the
width of the two branches is 0.7. RBCs divide equally at the vessel
bifurcation. The velocity field is indicated by the vector field. Figure 14:
Cells are set in a cluster initially under a moderate aggregating force with
$\alpha=25,q_{1}=1,q_{2}=1$. The top branch is set to be 1 with the other
stays 0.7. One of the cells is going into the branched vessel. The velocity
field is shown as well. Figure 15: Cells are set in a cluster initially under
a strong aggregating force with $\alpha=1.5\times 10^{3},1_{1}=1,Q_{2}=1$. The
set up is the same as Figure 14. None of the cells is going into the branched
vessel. The velocity field is shown as well.

##  5 Conclusions

In this paper, a thermodynamically consistent diffuse interface (phase-field)
model for describing cell-wall, cell-cell interaction and aggregation is
derived based on the energy variation method. The interactions of cellular and
wall structures are modeled by introducing a new interaction energy with the
help of the phase-field function.

Then an efficient scheme using $C^{0}$ finite element spatial discretization
and the mid-point temporal discretization is proposed to solve the obtained
model equations. Thanks to the mid-point temporal discretization, the obtained
numerical scheme is unconditionally energy stable. The model and parameters
are calibrated with experimental data on cell deformation under different
stretch force. Then effects of the adhesion strength on cell-wall and cell-
cell interaction are studied. In the end, the calibrated model is used to
model red blood cell motion near the vessel bifurcation. The results are
consistent with the experimental observation.

## References

  * [1] T. Kirschkamp, H. Schmid-Schönbein, A. Weinberger, R. Smeets, Effects of fibrinogen and $\alpha$2-macroglobulin and their apheretic elimination on general blood rheology and rheological characteristics of red blood cell aggregates, Therapeutic Apheresis and Dialysis 12 (5) (2008) 360–367.
  * [2] A. L. Fogelson, K. B. Neeves, Fluid mechanics of blood clot formation, Annual review of fluid mechanics 47 (2015) 377–403.
  * [3] K. Lee, M. Kinnunen, M. D. Khokhlova, E. V. Lyubin, A. V. Priezzhev, I. Meglinski, A. A. Fedyanin, Optical tweezers study of red blood cell aggregation and disaggregation in plasma and protein solutions, Journal of biomedical optics 21 (3) (2016) 035001.
  * [4] Y. Alapan, J. A. Little, U. A. Gurkan, Heterogeneous red blood cell adhesion and deformability in sickle cell disease, Scientific reports 4 (1) (2014) 1–8.
  * [5] Y. Connor, Y. Tekleab, S. Tekleab, S. Nandakumar, D. Bharat, S. Sengupta, A mathematical model of tumor-endothelial interactions in a 3d co-culture, Scientific reports 9 (1) (2019) 1–14.
  * [6] Y. Deng, D. P. Papageorgiou, X. Li, N. Perakakis, C. S. Mantzoros, M. Dao, G. E. Karniadakis, Quantifying fibrinogen-dependent aggregation of red blood cells in type 2 diabetes mellitus, Biophysical journal 119 (5) (2020) 900–912.
  * [7] Flow patterns around heart valves: A numerical method, Journal of Computational Physics 10 (2) (1972) 252–271.  doi:https://doi.org/10.1016/0021-9991(72)90065-4.
URL https://www.sciencedirect.com/science/article/pii/0021999172900654

  * [8] C. S. Peskin, The immersed boundary method, Acta Numerica 11 (2002) 479–517.  doi:10.1017/S0962492902000077.
  * [9] A. L. Fogelson, R. D. Guy, Immersed-boundary-type models of intravascular platelet aggregation, Computer methods in applied mechanics and engineering 197 (25-28) (2008) 2087–2104.
  * [10] A. Gerisch, M. A. Chaplain, Mathematical modelling of cancer cell invasion of tissue: local and non-local models and the effect of adhesion, Journal of Theoretical Biology 250 (4) (2008) 684–704.
  * [11] P. Ziherl, S. Svetina, Flat and sigmoidally curved contact zones in vesicle–vesicle adhesion, Proceedings of the National Academy of Sciences 104 (3) (2007) 761–765.
  * [12] Y. Liu, W. K. Liu, Rheology of red blood cell aggregation by computer simulation, Journal of Computational Physics 220 (1) (2006) 139–154.
  * [13] P. Ziherl, Aggregates of two-dimensional vesicles: Rouleaux, sheets, and convergent extension, Physical review letters 99 (12) (2007) 128102.
  * [14] J. Zhang, P. C. Johnson, A. S. Popel, Red blood cell aggregation and dissociation in shear flows simulated by lattice boltzmann method, Journal of biomechanics 41 (1) (2008) 47–55.
  * [15] D. A. Fedosov, H. Noguchi, G. Gompper, Multiscale modeling of blood flow: from single cells to blood rheology, Biomechanics and modeling in mechanobiology 13 (2) (2014) 239–258.
  * [16] Z. Xu, N. Chen, S. C. Shadden, J. E. Marsden, M. M. Kamocka, E. D. Rosen, M. Alber, Study of blood flow impact on growth of thrombi using a multiscale model, Soft Matter 5 (4) (2009) 769–779.
  * [17] D. Flormann, O. Aouane, L. Kaestner, C. Ruloff, C. Misbah, T. Podgorski, C. Wagner, The buckling instability of aggregating red blood cells, Scientific reports 7 (1) (2017) 1–10.
  * [18] B. Quaife, S. Veerapaneni, Y.-N. Young, Hydrodynamics and rheology of a vesicle doublet suspension, Physical Review Fluids 4 (10) (2019) 103601.
  * [19] B. Quaife, S. Veerapaneni, Y.-N. Young, Hydrodynamics and rheology of a vesicle doublet suspension, Physical Review Fluids 4 (10) (2019) 103601.
  * [20] D. M. Anderson, G. B. McFadden, A. A. Wheeler, Diffuse-interface methods in fluid mechanics, Annual Review of Fluid Mechanics 30 (1) (1998) 139–165.  arXiv:https://doi.org/10.1146/annurev.fluid.30.1.139, doi:10.1146/annurev.fluid.30.1.139.
URL https://doi.org/10.1146/annurev.fluid.30.1.139

  * [21] J. Lowengrub, L. Truskinovsky, Quasi-incompressible cah-hilliard fluids and topological transitions, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 454 (1978) (1998) 2617–2654.  doi:10.1098/rspa.1998.0273.
  * [22] Q. Du, C. Liu, R. Ryham, X. Wang, Energetic variational approaches in modeling vesicle and fluid interactions, Physica D: Nonlinear Phenomena 238 (9-10) (2009) 923–930.
  * [23] Q. Du, C. Liu, R. Ryham, X. Wang, Modeling the spontaneous curvature effects in static cell membrane deformations by a phase field formulation, Communications on Pure & Applied Analysis 4 (3) (2005) 537.
  * [24] X. Yang, V. Mironov, Q. Wang, Modeling fusion of cellular aggregates in biofabrication using phase field theories, Journal of theoretical biology 303 (2012) 110–118.
  * [25] R. Gu, X. Wang, M. Gunzburger, Simulating vesicle–substrate adhesion using two phase field functions, Journal of Computational Physics 275 (2014) 626–641.
  * [26] J. Zhang, S. Das, Q. Du, A phase field model for vesicle–substrate adhesion, Journal of Computational Physics 228 (20) (2009) 7837–7849.
  * [27] S. Das, Q. Du, Adhesion of vesicles to curved substrates, Physical Review E 77 (1) (2008) 011907.
  * [28] R. Gu, X. Wang, M. Gunzburger, A two phase field model for tracking vesicle–vesicle adhesion, Journal of mathematical biology 73 (5) (2016) 1293–1319.
  * [29] J. Jiang, K. Garikipati, S. Rudraraju, A diffuse interface framework for modeling the evolution of multi-cell aggregates as a soft packing problem driven by the growth and division of cells, Bulletin of mathematical biology 81 (8) (2019) 3282–3300.
  * [30] R. Gu, X. Wang, M. Gunzburger, A two phase field model for tracking vesicle–vesicle adhesion, Journal of mathematical biology 73 (5) (2016) 1293–1319.
  * [31] Q. Du, C. Liu, R. Ryham, X. Wang, A phase field formulation of the willmore problem, Nonlinearity 18 (3) (2005) 1249.
  * [32] F. Guillén-González, G. Tierra, Unconditionally energy stable numerical schemes for phase-field vesicle membrane model, Journal of computational physics 354 (2018) 67–85.
  * [33] Q. Du, J. Zhang, Adaptive finite element method for a phase field bending elasticity model of vesicle membrane deformations, SIAM Journal on Scientific Computing 30 (3) (2008) 1634–1657.
  * [34] Q. Cheng, J. Shen, Multiple scalar auxiliary variable (msav) approach and its application to the phase-field vesicle membrane model, SIAM Journal on Scientific Computing 40 (6) (2018) A3982–A4006.
  * [35] J. Hua, P. Lin, C. Liu, Q. Wang, Energy law preserving c0 finite element schemes for phase field models in two-phase flow computations, Journal of Computational Physics 230 (19) (2011) 7115–7131.
  * [36] P. Lin, C. Liu, Simulations of singularity dynamics in liquid crystal flows: A c0 finite element approach, Journal of Computational Physics 215 (1) (2006) 348–362.
  * [37]
  * [38] Z. Guo, P. Lin, A thermodynamically consistent phase-field model for two-phase flows with thermocapillary effects, Journal of Fluid Mechanics 766 (2015) 226–271.
  * [39] L. Shen, Z. Xu, P. Lin, H. Huang, S. Xu, An energy stable c^0 finite element scheme for a phase-field model of vesicle motion and deformation, SIAM Journal on Scientific Computing 44 (1) (2022) B122–B145.
  * [40] T. Qian, X.-P. Wang, P. Sheng, A variational approach to the moving contact line hydrodynamics, arXiv preprint cond-mat/0602293 (2006).
  * [41] Q. Du, C. Liu, X. Wang, A phase field approach in the numerical study of the elastic bending energy for vesicle membranes, Journal of Computational Physics 198 (2) (2004) 450–468.
  * [42] R. Chen, G. Ji, X. Yang, H. Zhang, Decoupled energy stable schemes for phase-field vesicle membrane model, Journal of Computational Physics 302 (2015) 509–523.
  * [43] S. Aland, S. Egerer, J. Lowengrub, A. Voigt, Diffuse interface models of locally inextensible vesicles in a viscous fluid, Journal of computational physics 277 (2014) 32–47.
  * [44] Z. Guo, F. Yu, P. Lin, S. Wise, J. Lowengrub, A diffuse domain method for two-phase flows with large density ratio in complex geometries, Journal of Fluid Mechanics 907 (2021).
  * [45] K. C. Ong, M.-C. Lai, An immersed boundary projection method for simulating the inextensible vesicle dynamics, Journal of Computational Physics 408 (2020) 109277\.
  * [46] S. Xu, M. Alber, Z. Xu, Three-phase model of visco-elastic incompressible fluid flow and its computational implementation, Communications in computational physics 25 (2) (2019) 586.
  * [47] B. Eisenberg, Y. Hyon, C. Liu, Energy variational analysis of ions in water and channels: Field theory for primitive models of complex ionic fluids, The Journal of Chemical Physics 133 (10) (2010) 104104.
  * [48] Y. Hyon, C. Liu, et al., Energetic variational approach in complex fluids: maximum dissipation principle, Discrete & Continuous Dynamical Systems 26 (4) (2010) 1291.
  * [49] S. Xu, P. Sheng, C. Liu, An energetic variational approach for ion transport, Communications in Mathematical Sciences 12 (4) (2014) 779–789.
  * [50] M.-H. Giga, A. Kirshtein, C. Liu, Variational modeling and complex fluids, 2017\.
  * [51] S. Xu, B. Eisenberg, Z. Song, H. Huang, Osmosis through a semi-permeable membrane: a consistent approach to interactions, arXiv preprint arXiv:1806.00646 (2018).
  * [52] T. Qian, X.-P. Wang, P. Sheng, A variational approach to moving contact line hydrodynamics, Journal of Fluid Mechanics 564 (2006) 333.  doi:10.1017/S0022112006001935.
  * [53] C. Liu, H. Wu, An energetic variational approach for the cahn–hilliard equation with dynamic boundary condition: model derivation and mathematical analysis, Archive for Rational Mechanics and Analysis (2019) 1–81.
  * [54] Z. Guo, P. Lin, J. S. Lowengrub, A numerical method for the quasi-incompressible cahn–hilliard–navier–stokes equations for variable density flows with a discrete energy law, Journal of Computational Physics 276 (2014) 486–507.
  * [55] J. Mills, L. Qie, M. Dao, C. Lim, S. Suresh, Nonlinear elastic and viscoelastic deformation of the human red blood cell with optical tweezers, Molecular & Cellular Biomechanics 1 (3) (2004) 169.
  * [56] Z. Wu, Z. Xu, O. Kim, M. Alber, Three-dimensional multi-scale model of deformable platelets adhesion to vessel wall in blood flow, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 372 (2021) (2014) 20130380.