Phase field model of cell motility: sharp interface limit in sub-critical case The work of LB was supported by NSF grant DMS-1106666. The work of VR and MP was partially supported by NSF grant DMS-1106666. The authors are grateful to I. Aronson and F. Ziebert for useful discussions on the phase field model of cell motility introduced in their paper.

Phase field model of cell motility: sharp interface limit in sub-critical case thanks: The work of LB was supported by NSF grant DMS-1106666. The work of VR and MP was partially supported by NSF grant DMS-1106666. The authors are grateful to I. Aronson and F. Ziebert for useful discussions on the phase field model of cell motility introduced in their paper.

Leonid Berlyand L. Berlyand Department of Mathematics, The Pennsylvania State University, University Park, PA 16802, USA
22email: berlyand@math.psu.eduV. Rybalko Mathematical Division, B. Verkin Institute for Low Temperature, Physics and Engineering of National Academy of Sciences of Ukraine, 47 Lenin Ave., 61103 Kharkiv, Ukraine
44email: vrybalko@ilt.kharkov.uaM. Potomkin The Pennsylvania State University, University Park, PA 16802, USA,
66email: mup20@ucs.psu.edu
   Volodymyr Rybalko L. Berlyand Department of Mathematics, The Pennsylvania State University, University Park, PA 16802, USA
22email: berlyand@math.psu.eduV. Rybalko Mathematical Division, B. Verkin Institute for Low Temperature, Physics and Engineering of National Academy of Sciences of Ukraine, 47 Lenin Ave., 61103 Kharkiv, Ukraine
44email: vrybalko@ilt.kharkov.uaM. Potomkin The Pennsylvania State University, University Park, PA 16802, USA,
66email: mup20@ucs.psu.edu
   Mykhailo Potomkin L. Berlyand Department of Mathematics, The Pennsylvania State University, University Park, PA 16802, USA
22email: berlyand@math.psu.eduV. Rybalko Mathematical Division, B. Verkin Institute for Low Temperature, Physics and Engineering of National Academy of Sciences of Ukraine, 47 Lenin Ave., 61103 Kharkiv, Ukraine
44email: vrybalko@ilt.kharkov.uaM. Potomkin The Pennsylvania State University, University Park, PA 16802, USA,
66email: mup20@ucs.psu.edu
Abstract

We consider a system of two PDEs arising in modeling of motility of eukariotic cells on substrates. This system consists of the Allen-Cahn equation for the scalar phase field function coupled with another vectorial parabolic equation for the orientation of the actin filament network.

The two key properties of this system are (i) presence of gradients in the coupling terms (gradient coupling) and (ii) mass (volume) preservation constraints. We first prove that the sharp interface property of initial conditions is preserved in time. Next we formally derive the equation of the motion of the interface, which is the mean curvature motion perturbed by a nonlinear term that appears due to the properties (i)-(ii). This novel term leads to surprising features of the the motion of the interface.

Because of these properties maximum principle and classical comparison techniques do not apply to this system. Furthermore, the system can not be written in a form of gradient flow, which is why recently developed -convergence techniques also can not be used for the justification of the formal derivation. Such justification is presented in a one-dimensional model problem and it leads to a stability result in a class of “sharp interface” initial data.

Keywords:
Ginzburg-Landau phase field system with gradient coupling mass conservation cell motility.

1 Introduction

The problem of cell motility has been a classical subject in biology for several centuries. It dates back to the celebrated discovery by van Leeuwenhoek in 17th century who drastically improved microscope to the extent that he was able to observe motion of single celled organisms that moved due contraction and extension. Three centuries later this problem continues to attract the attention of biologists, biophysicists and, more recently, applied mathematicians. A comprehensive review of the mathematical modeling of cell motility can be found in Mog09 ().

This work is motivated by the problem of motility (crawling motion) of eukariotic cells on substrates. The network of actin (protein) filaments (which is a part of cytoskeleton in such cells) plays an important role in cell motility. We are concerned with the cell motility caused by extension of front of the cell due to polymerization of the actin filaments and contraction of the back of the cell due to detachment of these filaments. Modeling of this process in full generality is at present a formidable challenge because several important biological ingredients (e.g., regulatory pathways Mog09 ()) are not yet well understood.

However, in recent biophysical studies several simplified phase field models have been proposed. Simulations performed for these models demonstrated good agreement with experiments (e.g., ShaRapLev10 (); ZieSwaAra11 () and references their in). Recall that phase field models are typically used to describe the evolution of an interface between two phases (e.g., solidification or viscous fingering). The key ingredient of such models is an auxiliary scalar field, which takes two different values in domains describing the two phases (e.g., and ) with a diffuse interface of a small (non zero) width. The advantage of such an approach is that it allows us to consider one set of PDEs in the whole domain occupied by both phases and therefore avoids an issue of coupling two different sets of PDEs in each phase, which is typically quite complicated in both simulations and analysis.

In this work we present rigorous mathematical analysis of the 2D phase field model proposed in ZieSwaAra11 () that consists of a system of two PDEs for the phase field function and orientation vector with an integral mass conservation constraint. This model can be rewritten in a simplified form suitable for asymptotical analysis, so that all key features of the qualitative behavior are preserved, which can be seen from a comparison of simulations from ZieSwaAra11 () with our analytical results. First, in ZieSwaAra11 () the integral mass conservation constraint is introduced in the PDE system by adding a penalization parameter into the double-well potential (formulas (2.2), (2.5)-(2.6) in ZieSwaAra11 ()). We introduce this constraint via a dynamic Lagrange multiplier defined below, which provides the same qualitative behavior of solutions. Second, for technical simplicity we drop two terms in the second equation (for polarization) in the phase field system ZieSwaAra11 (), since our analysis shows that these terms can be incorporated with minor changes in both the results and the techniques. Thirdly, in order to study the long term behavior of the system, we perform the diffusive scaling (, ). Indeed, the crawling motion is very slow and time variable needs to be “accelerated”. Thus, we arrived at the following system of parabolic PDEs for a scalar phase field function and the orientation vector :

(1)
(2)

On the boundary we impose the Neumann and the Dirichlet boundary conditions respectively and , where is a bounded smooth domain.

Equation (1) is a perturbation of the following Allen-Cahn equation AllCah97 ():

(3)

The latter equation is a scalar version of the celebrated Ginzburg-Landau equation and it plays a fundamental role in mathematical modeling of phase transitions. It consists of the standard linear parabolic equation and a nonlinear lower order term, which is the derivative of a smooth double equal well potential

(4)

Equation (3) was introduced to model the motion of phase-antiphase boundary (interface) between two grains in a solid material. Analysis of (3) as led to the asymptotic solution that takes values and in the domains corresponding to two phases separated by an interface of the width of order , the so-called sharp interface. Furthermore, it was shown that this sharp interface exhibits the mean curvature motion. Recall that in this motion the normal component of the velocity of each point of the surface is equal to the mean curvature of the surface at this point. This motion has been extensively studied in the geometrical community (e.g., Ham82 (); Hui84 (); Gra87 (); Bra78 () and references therein). It also received significant attention in PDE literature. Specifically CheGigGot91 () and EvaSpr91 () established existence of global viscosity solutions (weak solutions) for the mean curvature flow. The mean curvature motion of the interface in the limit was formally derived in RubSterKel89 (), Fif88 () and then justified in EvaSonSou92 () by using the viscosity solutions techniques.

Note that equation (3) is closely related to another well-known model of phase separation, the so-called Cahn-Hilliard equation CahHil58 (), which is a forth order reaction diffusion equation that models how two components of a binary fluid spontaneously separate and form domains of two pure fluids.

There are two distinguishing features in the problem (1)-(2): coupling and a nonlocal mass conservation constraint. We first comment on the coupling. Note that another prominent biological FitzHugh-Nagumo model has similar coupling feature but in (1)-(2) the coupling occurs via spatial gradients (gradient coupling) of the unknown functions where as in FitzHugh-Nagumo RocGeoGiu00 (), SorSou96 () the two equations are coupled via lower order terms (unknown functions rather than their derivatives). There are several phase field models for Allen-Cahn (also Cahn-Hilliard) equation coupled with another parabolic equation via lower order terms EvaSonSou92 (); Che94 (). Our analysis shows that the gradient coupling results in novel mathematical features such as the following nonlinear nonlocal equation for the velocity of the interface curve derived below:

(5)

Here stands for the normal velocity of with respect to the inward normal, and denotes the curvature of , is a constant determined by the potential ( for the specific choice (4) of the double-well potential), is the curve length, and function is given by (82).

Next note that the term in (1) is a Lagrange multiplier responsible for the volume constraint (conservation mass in the original physical problem ZieSwaAra11 ()) and it has the following form

(6)

Solutions of stationary Allen-Cahn equation with the volume constraint were studied in Mod86 () by -convergence techniques applied to the stationary variational problem corresponding to(3). It was established that the -limiting functional is the interface perimeter (curve length in 2D or surface area in higher dimensions). Subsequently in the work RubSte92 () an evolutionary reaction-diffusion equation with double-well potential and nonlocal term that describes the volume constraint was studied. The following asymptotic formula for evolution of interface in the form of volume preserving mean curvature flow was formally derived in RubSte92 ()

(7)

Formula (7) was rigorously justified in the radially symmetric case in BroSto97 () and in general case in CheHilLog10 ().

There are three main approaches to the study of asymptotic behavior (sharp interface limit) of solutions of phase field equations and systems.

When comparison principle for solutions applies, a PDE approach based on viscosity solutions techniques was successfully used in EvaSonSou92 (); BarSonSou93 (); Sou02 (); LioKimSle04 (); BarLio03 (); Gol97 (). This approach can not be applied to the system (1)-(2), because of the gradient coupling and nonlocal multiplier . It is an open issue to introduce weak (e.g., viscosity or Brakke type) solutions in problems with constraints. Furthermore, since there is no comparison principle, the only technique available for the justification of the sharp interface limit is energy bounds which become quite difficult due to the coupling and the volume preservation.

Another technique used in such problems is convergence (see Ser10 () and references therein). It also does not work for the system (1)-(2). Standard Allen-Cahn equation (3) is a gradient flow (in metric) with GL energy functional, which is why one can use the convergence approach. However, there is no energy functional such that problem (1)-(2) can be written as a gradient flow.

As explained above the gradient coupling and the volume constraint are the key features of the problem (1), (2), (6) and they led to both novel results and analysis techniques. Specifically, the objectives of our study are three fold:

(i) To show that there is no finite time blow up and the sharp interface property of the initial data propagates in time.

(ii) To investigate how the gradient coupling combined with the nonlocal volume constraints changes the limiting equation of the interface motion.

(iii) To develop novel techniques for the justification of the limiting equation of the interface motion that, in particular, includes rigorous derivation of asymptotic expansion for the solution of the problem (1)-(2).

The paper is organized as follows. In section 2 is devoted to the objective (i). Section 3 the objectives (ii) and (iii) are addressed in the context of a model one-dimensional problem. In Section 4 the equation for the interface motion (5) is formally derived.

2 Existence of the sharp interface solutions that do not blow up in finite time

In this section we consider the boundary value problem (1), (2) with given by (6). Introduce the following functionals

(8)

The system (1)-(2) is supplied with “well prepared” initial data, which means that two conditions hold:

(9)

and

(10)

The first condition (9) is a weakened form of a standard condition for the phase field variable . If , then by the maximum principle implies for . The presence of nontrivial leads to an “extended interval” for . Second condition (10) means that at the function has the structure of “-transition layer”, that is the domain consists of three subdomains: a subdomain where the phase field function is close to (inside the cell) an another subdomain where (outside the cell) separated by a transition region of width . Furthermore, the orientation field has value close to everywhere except the -transition region.

Theorem 2.1

If the initial data , satisfy (9) and (10), then for any the solution , exists on the time interval for sufficiently small , . Moreover, it satisfies and

(11)

where is independent of and .

For the proof of this theorem we refer to BerPotRyb2016 ().

3 1D model problem: rigorous derivation of the sharp interface limit and remarks on stability

In Section 4 we present derivation of the formal asymptotic expansion of the solution of (1)-(2) and use it to obtain the equation of motion (5), which is the principal object of interest in the study of cell motility. There are two main sources of difficulties in rigorous justification of this derivation (i) the possible non-smoothness of limiting velocity field and (ii) dimension greater than one is much harder to handle technically. That is why in this Section we consider a simplified one-dimensional model and impose smallness assumption on that guarantees regularity of asymptotic solutions as well as plays an important technical role in our proof. Recently, we have developed another approach of justification of Sharp Interface Limit for the system (12)-(13) which is valid for arbitrary , see BerPotRyb2016 ().

Specifically, we study the limiting behavior as of the solution of the system

(12)
(13)

assuming that the initial data , has the following (“very well-prepared”) form

and

where . Here functions are solutions of (25) for and (29) for , and are solutions of (28) for and (30) for . The functions which are involved in the definition of are defined by (31) and (33) ( is the expansion the velocity of the cell’s sharp interface). We also assume that there exists a constant , independent of , such that

(14)

We emphasize that in the RHS of (12) is a given function rather than an unknown Lagrange multiplier in (1). The main distinction of 1D case is because in 1D there is no motion due to curvature of the interface since the interface is a point. Therefore if we take initial data such that the domain is a finite interval, then such a “one dimensional mathematical cell” will not move, which corresponds to a well known fact that the motion in the one-dimensional Allen-Cahn problem is exponentially slow (that is very different from AC in higher dimensions). Thus, we choose initial data to be a step like function that is a transition from an unbounded left interval to an unbounded right interval (the “cell”). In the 2D problem Lagrange multiplier appeared due to the mass (volume) conservation constraint in a finite domain occupied by the cell, which has no analog in 1D because the “mathematical cell” must be infinite as explained above. Therefore an analog of the Lagrange multiplier in 1D is chosen to be a given forcing term .

Hereafter denotes the classical standing wave solution of

(15)

with step-like conditions at infinity

(16)

Since the solution of (15)-(16) is uniquely defined up to translations (shifts), we impose the following normalization

(17)

In the particular case of the double-well potential having the form (4) the solution is explicitly given by . Note that solves (12) if the coupling term and are both identically zero.

The main result of this sections is the following theorem:

Theorem 3.1

Let and solve (12) and (13) on with initial conditions satisfying (14). Assume also that is a given smooth function and , where depends on the potential only. Then we have, for sufficiently small ,

(18)

where denotes the location of the interface between and phases and

(19)

Moreover, converges to which solves the interface motion equation similar to (5):

(20)

where , and is defined as solution of (79).

Remark 1

From the definition of function it follows that it depends linearly on the parameter so that function does not depend on .

Remark 2

While Theorem 3.1 describes the leading term of the asymptotic expansion for , in the course of the proof we also construct the leading term of the asymptotic expansion of in the form .

Remark 3

(on stability) Equation (20) is rigorously derived when but it could be formally derived for any real . This equation has the unique smooth solution when , for some . Roughly speaking, if , then is determined by (20) due to the implicit function theorem otherwise multiple solutions may appear. Thus, assumption on smallness of can be viewed as a stability condition. By contrast, for large enough one can observe instability due to the fact that the limiting equation (20) has multiple solutions and, therefore, perturbation of initial data may result in switching between multiple solutions of equation (20). Indeed, to explain we rewrite equation (20)

(21)

where the left hand side of (21) can be resolved in , but not uniquely.

Remark 4

The estimate (19) justifies the asymptotic expansion (18) and this estimate is the principal claim of this Theorem. However, in the course of the proof we actually obtain and justify a more precise asymptotic expansion of the from , which corresponds to and in the the expansions (23) below.

The proof of Theorem 3.1 is divided into two steps, presented in the following two subsections.

In the first step (Subsection 3.1) we formally construct approximate solution of the order and obtain equations for the residuals

(22)

where

(23)

for some . It would be natural to expect that , however, it turned out that due to the gradient coupling and nonlinearity of the problem, we can only prove boundedness of and for some .

The second step (Subsection 3.2) is the central mathematical part of this paper, and we briefly outline its main ideas. The goal there is to obtain bounds on residuals and for appropriate and . The bounds on play central role and they imply bounds on though these bounds are coupled. Therefore the bound (19) is the main claim of the Theorem.

The techniques of asymptotic expansions that include bounds on residuals were first developed for Allen-Cahn PDE in MotSha95 (). The proofs in MotSha95 () are based on the lower bound of the spectrum of linearized self-adjoint stationary Allen-Cahn operator in an unbounded domain. The techniques of this type were subsequently developed and applied in Alikakos, Bates, Chen [1] for the Cahn�Hillard equation, Caginalp and Chen [6] for the phase field system, and CheHilLog10 () for volume preserving Allen-Cahn PDE.

In the system (1)-(2) or its one-dimensional analog (12)-(13) the corresponding linearized operator is not self-adjoint and the previously developed techniques can not be directly applied.

The results of this Section are based on the analysis of a time-dependent linearized problem that corresponds to the entire system. We represent (split) the residual as a some of two parts , where and are new unknown functions. The first part is easer to estimate since it is chosen to be orthogonal to the eigenfunction corresponding to the zero eigenvalue of the linearized stationary Allen-Cahn operator, but is has a more general form than the second one since depends on both spatial and time variables. The second part is of a simpler form because does not depend on , but it contains the eigenfunction . The difficulty of dealing with such an eigenfunction can be explained by analyzing the equations (36) and (38) obtained by rescaling of the spatial variable. Note that the two highest order terms in (36) dominate other terms, and the sum of these two terms is nothing but linearized stationary Allen-Cahn operator. If in (36) one takes , then the terms vanish ( is an eigenfunction) and one has to estimate by analyzing the lower order terms, which is a much harder task. For example, in order to estimate we represent in (50) as a sum of two parts corresponding to the representation of and write down the leading term for the second part (first term in (53)). The justification of expansion with this leading term is a subtle task because it requires bounds on both and , which leads to a condition on smallness of .

3.1 Construction of asymptotic expansions.

First, we seek formal approximations for and in the form:

We also assume admits a power series expansion,

(24)

so that we also have expansion for the velocity ,

Next we expand ,

where

Substitute the series (35) and (24) into (12)-(13), and equate terms of like powers of to obtain that and for satisfy

(25)
(26)
(27)

and

(28)

The equations for have the following form

(29)
(30)
Remark 5

Due to the fact that is an eigenfunction of the linearized Allen-Cahn operator corresponding to the zero eigenvalue, the following solvability conditions for (26),(27) and (29) arise

(31)
(32)
(33)

and uniquely define the functions , such that the solvability conditions (33) are satisfied. Equations (31), (32) and (33) are solvable for , and , respectively, for sufficiently small . Also we note that once the solvability conditions are satisfied then equations (26),(27) and (29) have a family of solutions: , , where is a particular solution. We choose s.t.

Define functions and by

(34)

where

(35)

Substituting the representation for from (34) into (12) we derive the PDE for (note that the differentiation in time and new spatial variable are no longer independent)

(36)

where is of the form

(37)

where , are bounded functions in , and and square integrable with respect to (except ). Moreover, the function is orthogonal to :

The functions , , and are expressed in terms of , and , their exact form, which is not important for the proof, is given in Appendix.

Substituting the representation for from (34) into (13) we derive also the PDE for :

(38)

For more details on derivation of (36) and (38) we refer to Appendix.

3.2 Bounds for residuals and

In this section we obtain bounds for the coupled system of PDEs:

To this end we write the unknown function in the following form,

(39)

Then (36) becomes

(40)
Lemma 1

The following inequality holds

(41)

P r o o f:
Multiply (40) by and integrate to obtain

(42)

In order to derive (41) we simplify equality (42). First, we notice that the integral in the first term can be rewritten as follows

(43)

We used here the definition of , i.e., .

The first term in the right hand side of (42) vanishes. Indeed,

The second term in the left hand side can be split into three terms:

(44)

Rewrite the first term in (44), using (25) and integrating by parts,

Thus,