Phase field model of cell motility: sharp interface limit in sub-critical case

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 Gamma-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.


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 [22].
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 [22]) 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., [10,13] 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., 1 and 0) 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 [13] 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 [13] with our analytical results.First, in [13] 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 [13]).We introduce this constraint via a dynamic Lagrange multiplier λ (t) 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 [13], 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 (t → ε 2 t, x → εx).
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 P ε : On the boundary ∂ Ω we impose the Neumann and the Dirichlet boundary conditions respectively ∂ ν ρ ε = 0 and P ε = 0, where Ω ⊂ R 2 is a bounded smooth domain.
Equation ( 1) is a perturbation of the following Allen-Cahn equation [1]: 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 Equation ( 3) was introduced to model the motion of phase-antiphase boundary (interface) between two grains in a solid material.Analysis of (3) as ε → 0 led to the asymptotic solution that takes values 0 and 1 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., [17,18,16,4] and references therein).It also received significant attention in PDE literature.Specifically [9] and [11] established existence of global viscosity solutions (weak solutions) for the mean curvature flow.The mean curvature motion of the interface in the limit ε → 0 was formally derived in [19], [14] and then justified in [12] 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 [6], 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 [24], [27] 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 [12,7].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 Γ (t) derived below: Here V stands for the normal velocity of Γ (t) with respect to the inward normal, and κ denotes the curvature of Γ (t), c 0 is a constant determined by the potential W (c 0 = 3/2 for the specific choice (4) of the double-well potential), |Γ (t)| is the curve length, and function Φ(V ) is given by (91).
Next note that the term λ ε (t) in ( 1) is a Lagrange multiplier responsible for the volume constraint (conservation mass in the original physical problem [13]) and it has the following form Solutions of stationary Allen-Cahn equation with the volume constraint were studied in [21] 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 [25] an evolutionary reaction-diffusion equation with doublewell 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 [25] κ ds (7) Formula ( 7) was rigorously justified in the radially symmetric case in [5] and in general case in [8].
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 [12,3,28,20,2,15].This approach can not be applied to the system (1)-( 2), because of the gradient coupling and nonlocal multiplier λ ε (t).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 [26] and references therein).It also does not work for the system (1)-( 2).Standard Allen-Cahn equation ( 3) is a gradient flow (in L 2 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 onedimensional problem.In Section 4 the equation for the interface motion ( 5) is formally derived.

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: and The first condition (9) is a weakened form of a standard condition for the phase field variable 0 ≤ ρ ε (x, 0) ≤ 1.If λ ε ≡ 0, then by the maximum principle 0 ≤ ρ ε (x, 0) ≤ 1 implies 0 ≤ ρ ε (x,t) ≤ 1 for t > 0. The presence of nontrivial λ ε leads to an "extended interval" for ρ ε .Second condition (10) means that at t = 0 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 1 (inside the cell) an another subdomain where ρ ε ∼ 0 (outside the cell) separated by a transition region of width ε.Furthermore, the orientation field P ε has value close to 0 everywhere except the ε-transition region.
Remark 1 This theorem implies that if the initial data are well-prepared in the sense of ( 9)-( 10) then for 0 < t < T the solution exists and has the structure of ε-transition layer.Moreover, the bound on initial data (9) remains true for t > 0, while it relies on maximum principle argument, it also requires additional estimates on λ ε as seen from ( 13) below.
Proof STEP 1.First multiply (1) by ∂ t ρ ε and integrate over Ω : (12) Here we used that due to (6) the integral of ∂ t ρ ε over Ω is zero and thus Next, using the maximum principle in (1) we get Let T ε > 0 be the maximal time such that and from now on we assume that t ≤ T ε .Using ( 12), ( 14) and integration by parts we obtain We proceed by deriving an upper bound for the integral in the right hand side of (15).By (1) we have The following bounds are obtained by routine applying the Cauchy-Schwarz and Young's inequalities.For the the sum of the first three terms in (16) we get, 3 Finally, in order to bound I 5 we first derive, then and using this inequality, ( 16) and ( 14) in ( 15), then substituting the resulting bound in ( 12) we obtain, for sufficiently small ε, STEP 2. Now we obtain a bound for the last two terms in (18).Taking the scalar product of (2) with 2kP ε + 4|P ε | 2 P ε , k > 0, and integrating over Ω we get We chose ) dx, by ( 18) and ( 19) we have the differential inequality, with a constant C > 0 independent of ε.Considering the bounds on the initial data and assuming that ε is sufficiently small, one can easily construct a bounded supersolution G of (20 for sufficienly small ε.By ( 13) and ( 17) we then conclude that T ε in ( 14) coincides actually with T when ε is small.Theorem is proved.
Remark 2 We briefly summarize the key points of the proof.Equation ( 1) has the form of Allen-Cahn PDE perturbed by a quadratic term P ε • ∇ρ ε and term λ ε which is constant in x.Therefore the blow up issue is determined by the competition between these two parts.The first part is a gradient flow corresponding to scalar Ginzburg-Landau energy.This energy admits well known estimates that rule out the blow up (by nonnegatiovity of energy).On the other hand even ODEs with quadratic nonlinearity are known to exhibit a finite time blow up.In presence of the second term blow up is ruled out by establishing differential inequality (20) (cf.standard Gronwall).The key difficulty here is due to the "fourth power" term 1  2 (12).Here a straightforward way would be to get L ∞ estimates on P ε directly, which is, however, not possible.Thus we first lower the power using integration by parts and maximum principle ( 14) to obtain a "third power term" in the RHS of (15).A subtle point here is that maximum principle gives (13) which contains λ ε and one needs estimate (17) on λ ε to get bound ( 14) on an extended interval.Combining ( 15) with ( 1)-( 2) leads to (20), which in turn rules out the blow up.

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 V 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.
Specifically, we study the limiting behavior as ε → 0 of the solution of the system assuming that the initial data ρ ε (x, 0), P ε (x, 0) has the following ("very well-prepared") form and Here functions θ i are solutions of (34) for i = 0 and (38) for i ≥ 1, and Ψ i are solutions of (37) for i = 0 and (39) for i ≥ 1.The functions V i which are involved in the definition of Ψ i are defined by ( 40) and (42) (V = V 0 + εV 1 + ... is the expansion the velocity of the cell's sharp interface).We also assume that there exists a constant C, independent of ε, such that We emphasize that F(t) in the RHS of ( 21) 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 ρ = 1 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 ρ = 0 to an unbounded right interval ρ = 1 (the "cell").In the 2D problem Lagrange multiplier λ ε (t) 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 F(t).
Hereafter θ 0 denotes the classical standing wave solution of with step-like conditions at infinity Since the solution of ( 24)-( 25) is uniquely defined up to translations (shifts), we impose the following normalization In the particular case of the double-well potential W having the form (4) the solution θ 0 is explicitly given by θ 0 = (1−tanh(x/ √ 8))/2.Note that ρ ε (x,t) = θ 0 (x/ε) solves (21) if the coupling term and F(t) are both identically zero.
The main result of this sections is the following theorem: Theorem 2 Let ρ ε and P ε solve (21) and ( 22) on [0, T ] with initial conditions satisfying (23).Assume also that F(t) is a given smooth function and |β | < β 0 , where β 0 > 0 depends on the potential W only. Then we have, for sufficiently small ε, where x ε (t) denotes the location of the interface between 0 and 1 phases and Moreover, x ε (t) converges to x 0 (t) which solves the interface motion equation similar to (5): where c 0 = (θ ′ 0 ) 2 dy, and Ψ 0 (y,V ) is defined as solution of (88).
Remark 3 From the definition of function Ψ 0 it follows that it depends linearly on the parameter β so that function Φ(V ) does not depend on β .
Remark 4 While Theorem 2 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 P ε in the form Ψ 0 x−x 0 (t) ε , ẋ0 (t) .
Remark 5 (on stability) Equation ( 29) is rigorously derived when |β | < β 0 but it could be formally derived for any real β .This equation has the unique smooth solution x 0 (t) when |β | < β * , for some β * > β 0 > 0. Roughly speaking, if β < β * , then x 0 (t) is determined by (29) due to the implicit function theorem otherwise multiple solutions x 0 (t) 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 ( 29) has multiple solutions and, therefore, perturbation of initial data may result in switching between multiple solutions of equation (29).Indeed, to explain we rewrite equation ( 29) where the left hand side of (30) can be resolved in V , but not uniquely.

Remark 6
The estimate (28) justifies the asymptotic expansion (27) 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 ρ ε = θ 0 ), which corresponds to N = 3 and α = 3 in the the expansions (32) below.
The proof of Theorem 2 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 N and obtain equations for the residuals where for some α ≥ 1.It would be natural to expect that α = N + 1, however, it turned out that due to the gradient coupling and nonlinearity of the problem, we can only prove boundedness of u ε and Q ε for some 1 < α < N + 1.
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 u ε and Q ε for appropriate α and N. The bounds on u ε play central role and they imply bounds on Q ε though these bounds are coupled.Therefore the bound (28) 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 [23].The proofs in [23] 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 CahnHillard equation, Caginalp and Chen [6] for the phase field system, and [8] for volume preserving Allen-Cahn PDE.
In the system (1)-( 2) or its one-dimensional analog ( 21)-( 22) 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 u ε as a some of two parts , where v ε (x/ε,t) and ξ ε (t) are new unknown functions.The first part is easer to estimate since it is chosen to be orthogonal to the eigenfunction θ ′ 0 (x/ε) corresponding to the zero eigenvalue of the linearized stationary Allen-Cahn operator, but is has a more general form than the second one since v ε (x/ε,t) depends on both spatial and time variables.The second part is of a simpler form because ξ ε (t) does not depend on x, but it contains the eigenfunction θ ′ 0 (x/ε).The difficulty of dealing with such an eigenfunction can be explained by analyzing the equations ( 45) and (47) obtained by rescaling of the spatial variable.Note that the two highest order ε −2 terms in (45) dominate other terms, and the sum of these two terms is nothing but linearized stationary Allen-Cahn operator.If in (45) one takes u ε (x,t) = θ ′ 0 (x/ε)ξ ε (t), then the ε −2 terms vanish (θ ′ 0 is an eigenfunction) and one has to estimate ξ ε (t) by analyzing the lower order terms, which is a much harder task.For example, in order to estimate ξ ε (t) we reperesent Q ε in (59) as a sum of two parts corresponding to the representation of u ε and write down the leading term for the second part (first term in (62)).The justification of expansion with this leading term is a subtle task because it requires bounds on both ξ ε (t) and ξε (t), which leads to a condition on smallness of β .

Construction of asymptotic expansions.
First, we seek formal approximations for ρ ε and P ε in the form: We also assume x ε (t) admits a power series expansion, so that we also have expansion for the velocity V = − ẋε , Substitute the series (44) and ( 33) into ( 21)-( 22), and equate terms of like powers of ε to obtain that θ i and Ψ i for i = 0, 1, 2 satisfy and The equations for i > 2 have the following form Remark 7 Due to the fact that θ ′ is an eigenfunction of the linearized Allen-Cahn operator corresponding to the zero eigenvalue, the following solvability conditions for (35),(36) and (38) arise and uniquely define the functions V i (t), i = 0, ..., N − 1 such that the solvability conditions (42) are satisfied.Equations (40), (41) and (42) are solvable for V 0 , V 1 and V i , respectively, for sufficiently small β .Also we note that once the solvability conditions are satisfied then equations ( 35),(36) and (38) have a family of solutions: , where θ i is a particular solution.We choose γ s.t.
The functions a ε , b i,ε , e ε and g ε are expressed in terms of θ i , V i and Ψ i , their exact form, which is not important for the proof, is given in Appendix.
Substituting the representation for P ε from (43) into (22) we derive also the PDE for For more details on derivation of (45) and (47) we refer to Appendix.

Bounds for residuals u ε and Q ε
In this section we obtain bounds for the coupled system of PDEs: To this end we write the unknown function u ε in the following form, Then (45) becomes

Lemma 1
The following inequality holds P r o o f: Multiply (49) by u ε = θ ′ 0 (v ε + ξ ε ) and integrate to obtain In order to derive (50) we simplify equality (51).First, we notice that the integral in the first term can be rewritten as follows We used here the definition of v ε , i.e., (θ ′ 0 ) 2 v ε dy = 0.The first term in the right hand side of (51) vanishes.Indeed, The second term in the left hand side can be split into three terms: Rewrite the first term in (53), using (34) and integrating by parts, Next we rewrite the second and the third terms in (53) using the fact that θ ′′′ 0 = W ′′ (θ 0 )θ ′ 0 (this latter equality is obtained by differentiating (34) with respect to y), Also, we make use the following equality which follows from (36), Finally, we use equalities ( 52),( 54),( 57),( 55) and (56) to rewrite (51) as follows Then ( 50) is obtained by applying the standard Cauchy-Schwarz inequality and the Poincare inequality (98) from Appendix.
It follows from the previous lemma that in order show the boundedness of u ε we need to find an appropriate upper bound on the term To this end, we use the equation ( 47) for Q ε .Note that Ψ 0 = Ψ 0 (y;V 0 ) depends on time t through V 0 (t).For simplicity of the presentation we suppress the second argument Ψ 0 (y) := Ψ 0 (y;V 0 ), if it equals to V 0 = − ẋ0 (t).
It follows from (37) functions Ψ ′ 0 and Ψ 0,V 0 := ∂Ψ 0 /∂V 0 solve the following equations Rewrite (47) substituting Thus, Q ε can be written as where A ε and B ε are solutions of the following problems, Next we note that the function A ε can be written as , where D ε solves the following problem, This representation (62) allows us to rewrite the term (58) as follows, The bounds for these terms are collected in the next lemma.

Lemma 2
The following inequalities hold (i) and the energy relation for Then ( 64) and ( 65) are obtained by applying the Cauchy-Schwarz inequality.Note that Ψ 0 depends linearly on β , this allows us to bound the right hand side of (68) by ε dy and this eventually leads to (65) .Item (iii) easily follows from representation (62).

Corrolary 1 The following inequality holds
Thus, we reduced the estimation of (58) to estimation of ξ 2 ε .The key observation leading to the desired bound on ξ ε can be explained now as follows.Observe that the presence of ξ ε in the RHS of (50) might result in the exponential growth of ξ ε (the best one can guarantee from ξ 2 ε ≤ C ε ξ 2 ε type bound).Fortunately, the term Ψ ′ 0 (θ ′ 0 ) 2 dyξ 2 ε in (50) cancels with the leading term appearing after substitution of expansion (62) for Q ε .However, this results in the appearance of lower order terms depending on ξε .Lemma 3 below provides the control of | ξε |.

Lemma 3 The following inequality holds
P r o o f: Multiply equation ( 49) by θ ′ 0 ξε and integrate in y to obtain Using (48) we simplify the left hand side of (71), In the next four steps inequality (70) will be derived by estimating the right hand side of (71) term by term.
STEP 1.The first term in the right hand side of (71) is estimated by using integration by parts and the Cauchy-Schwarz inequality, Here we introduced small parameter δ which does not depend on ε and will be chosen later.STEP 2. In this step we show that the sum of the second and the third terms in the right hand side of (71) gives zero.Indeed, use integration by parts and the definition of θ ′ 0 (θ STEP 3. In this step we estimate the sum of the fourth, the fifth and the sixth terms in the right hand side of (71), To this end we first show that Indeed, differentiating (35) in y one obtains the equation θ ′′′ , whose solvability condition reads The latter equality contains five terms.The second and the fifth term vanish since they are integrals of derivatives (( V 0 2 (θ ′ 0 ) 2 ) ′ and (F ′ (t)θ ′ 0 ) ′ , respectively).Integrating by parts the third term we get, This immediately implies (75).Now, taking into account the equality (75) and representation (62) for Q ε written in the form where we have also used integration by parts.
To prove (ii) we note that definition (46) of R ε , the Cauchy-Schwarz inequality, Poincare inequalities (98) and (101) yield Now it is convinient to introduce the folloowing notation, In terms of E ε and D ε we can rewrite (79) for α = 3 and N = 4 in the following form, Substituting the above inequality and (79) into (77) we obtain Assume that E ε (t) ≤ c for all t ∈ [0,t ⋆ ).
4 Formal derivation of (5) In this section we formally derive equation (5) for 2D system (1-2) with gradient coupling.Analogous derivation for the single Allen-Cahn equation has been done in [19], and [8].Consider a subdomain ω t ⊂ Ω ⊂ R 2 (ω t is occupied by the cell) so that Γ (0) = ∂ ω 0 (boundary of the cell at t = 0).For all t ∈ [0, T ] consider a closed curve Γ (t) s.t.∂ ω t = Γ (t) and ω t ⊂ Ω .Let X 0 (s,t) be a parametrization of Γ (t).In a vicinity of Γ (t) the parameters s and the signed distance r to Γ (t) will be used as local coordinates, so that x = X 0 (s,t) + rn(s,t) = X(r, s,t), where n is an inward normal.
The inverse mapping to x = X(r, s,t) is given by where in the formula for r we choose − if x ∈ ω t and +, if x / ∈ ω t .Recall that Γ (t) is the limiting location of interface as ε → 0. For fixed ε we are looking for interface in the form of ε-perturbation of Γ (t): Introduce the limiting velocity V 0 := −∂ t r and the distance to Γ ε (t) rescaled by ε: Next we define a rule that for all positive t transforms any given function w of original variable x into the corresponding function w in local coordinates (z, s): w(x,t) = w r(x,t) − εh ε (S(x,t),t) ε , S(x,t),t .

A Appendix
A.1 Derivation of equations for u ε and Q ε .
In the end of this subsection we write expression for functions in (46).

A.3 Uniqueness of original problem (1)-(2)
In this appendix we prove uniqueness of the solution to the problem (1)-( 2) in the following class: ε ∈ L 2 ((0,T ) × Ω ).ε ,P ε and ρ (2) ε ,P ε are solutions satisfying (104) and (105) for some positive T > 0. Take ρ ε = ρ 1 ε − ρ 2 ε and P ε = P Estimate the right had side using integration by parts: Thus, using ) P r o o f.Items (i) and (ii) are proved by means of energy relations that are obtained after multiplying (60) and (61) by B ε and A ε , respectively, and integrating in y.The resulting energy relation for the function B ε is