Table Of Content

A generalized Cahn-Hilliard equation for biological applications Evgeniy Khain1 and Leonard M. Sander2 1 Department of Physics, Oakland University, Rochester, Michigan 48309 and 2 Department of Physics and Michigan Center for Theoretical Physics, The University of Michigan, Ann Arbor, Michigan 48109 Recentlyweconsideredastochasticdiscretemodelwhichdescribesfrontsofcellsinvadingawound [1]. Inthemodelcellscanmove,proliferate,andexperiencecell-celladhesion. Inthisworkwefocus 8 on a continuum description of this phenomenon by means of a generalized Cahn-Hilliard equation 0 (GCH) with a proliferation term. As in the discrete model, there are two interesting regimes. For 0 subcriticaladhesion,therearepropagating”pulled”fronts,similarly tothoseofFisher-Kolmogorov 2 equation. Theproblemoffrontvelocityselectionisexamined,andourtheoreticalpredictionsarein n agood agreementwithanumericalsolutionoftheGCHequation. Forsupercriticaladhesion,there a is a nontrivial transient behavior, where density profile exhibits a secondary peak. To analyze this J regime, we investigated relaxation dynamics for the Cahn-Hilliard equation without proliferation. 6 We found that the relaxation process exhibits self-similar behavior. The results of continuum and 1 discrete models are in a good agreement with each other for thedifferent regimes we analyzed. ] PACSnumbers: 05.70.Ln,05.45.-a,87.18.Ed,87.18.Hf,05.50.+q,02.30.Jr h c e I. INTRODUCTION nearest-neighbor interactions between particles. Above m the criticalstrengthofcell-celladhesion(orbelowacrit- t- Inthispaperweproposeacontinuummethodfordeal- ical temperature in the Ising model), the homogeneous a ing with cells that move, proliferate and interact via ad- state becomes unstable, which leads to phase separation t s hesion. This problem arisesin models for wound healing between high density clusters and a dilute gas of par- t. [2] and tumor growth [3]. It is easy to formulate a dis- ticles. The dynamics of phase separation and coarsen- a m crete model for these processes [1]. However,proceeding ing (where largerclusters grow at the expense of smaller to the continuum limit is non-trivial [4]. ones) is usually described by the Cahn-Hilliard equation d- Considerfirstasimplediscretemodelfordiffusionand [13], a version of this equation can be derived directly n proliferation. Each lattice site can be empty or once oc- from the microscopic model [14]. o cupied. Ateachtimestep,aparticleispickedatrandom. Wecaneasilyaddproliferationtothislatticemodel[1], c Then it can either jump to a neighboring empty site, or sothatwehavediffusion,proliferation,andcell-celladhe- [ proliferate there (a new particle is born). We can ask sion. In this paper we suggestthat the proper candidate 1 for is the continuum analog of this model? It was shown for a continuum description is a Cahn-Hilliard equation v [5, 6, 7] that for small proliferation rates the propagat- with a proliferation term added, the GCH equation. 4 ingfrontsinthisdiscretesystemcanbedescribedbythe The rest of the paper is organized as follows. In Sec. 7 Fisher-Kolmogorovequation [8] (FK): 2 we present both discrete and continuum models which 5 includediffusion,proliferation,andcell-celladhesionand .2 ∂u =D¯ ∂2u +αu(1 u), (1) presenta generalphase diagram. Section3 describes the 1 ∂t ∂x2 − front propagation problem for subcritical adhesion. Sec- 0 where u is the (local) cell density, D¯ is the cell diffusion tion 4 focuses on a supercritical adhesion both for zero 8 and nonzero proliferation. Section 5 includes a brief dis- 0 coefficient, and α is the rate of proliferation [9]. Eq. (1) cussion and summary of our results. : admits solutions in a form of propagating fronts, but v the velocity selection is a nontrivial problem. There is i X a range of possible velocities; initially sufficiently local- II. DISCRETE AND CONTINUUM MODELS r izeddensityprofilesdevelopintopropagatingfrontswith a the “critical” velocity, v =2√D¯α [10]. Next consider a discrete model which includes diffu- A. Discrete model sion and cell-cell adhesion, but not proliferation. When a particle is picked at random, the probability to jump We review the discrete model for diffusion, prolifera- decreases with the number of nearest neighbors, to take tion, and cell-cell adhesion [1]. Consider a square two- into account cell-cell adhesion [1]. This scheme can be dimensional lattice in a channel geometry. The lattice mapped into the Ising model [11, 12], by identifying an distance is assumed to be equal to cell diameter, taking empty site with spin down, and an occupied site with into account hard-core exclusion. Initially, we put cells spin up. There is a simple relation between the aver- into the left part of the channel. We take x to measure age density, u, and the averagemagnetization, m, in the distance along the channel. A cell is picked at random, Ising model: u = (m+1)/2. The number of particles and one of the four neighboring sites is also picked at is fixed because there is no proliferation, and there are random. If this site is empty, the cell can proliferate 2 to this site (so that a new cell is born there), or mi- rameter q as a function of averagedensity u: c grate there. We denote the probability for proliferation 2 1/8 by α. Cell-cell adhesion is represented by a probability 1 1 16(1 qc) u= 1 − . (2) for migration that decreases with the number of nearest 2 ± 2 − q4 (cid:20) c (cid:21) neighbors: p = (1 α)(1 q)n, where 0 q < 1 is the adhesionmpigarramete−r, and −1 n 4 is th≤e number The unstable region corresponds to q > qc, so for su- ≤ ≤ percriticaladhesionthereisaphaseseparationandlarge of nearest neighbors. The case q =0 means no adhesion clusters are formed. Interestingly, even if we start with andreducestothemodelofRefs. [5,6,7]. Fornonzeroq, q < q , the initially homogeneous state can become un- itismuchhardertoacelltodiffuseifithasmanyneigh- c stable when one turns on proliferation, which leads to bors. After eachsteptime is advancedby 1/N,where N phase separation and clustering. is the current number of cells. Without proliferation the model can be mapped into B. Continuum approach the Ising model, as we pointed out in [1]. In this mapping the adhesion parameter q is identified with 1 exp( J/k T) where T is the temperature, k is Todescribethecoarse-graineddynamicsofthediscrete B B − − Boltzmann’s constant, and J is the coupling strength in model, we take a continuum approach. A model equa- the magnetic model, and the average density u is iden- tion, which describes the dynamics of phase separation tified with (m+1)/2, where m is the average magneti- with conserved order parameter (without proliferation) zation. The mapping is possible because our dynamical istheCahn-Hilliardequation[13,14]. Wenowformulate rules satisfy detailed balance. Therefore, the statics of the GCH equation by adding a proliferation term: our model is the same as in the Ising model. By statics, ∂u ∂2 ∂2u df we mean a phase diagram (m,T) (or (u,q) in our case) = ln(1 q) + +αu(1 u) (3) ∂t ∂x2 − ∂x2 du − which has stable and unstable regions. In the stable re- (cid:18) (cid:19) gion, a homogeneous state (with uniformly distributed whereuisthelocaldensity,f isthelocalfreeenergy,and cells) remains homogeneous; in contrast, in the unsta- α is the rate of proliferation. The gradient term in the ble regionphase separationoccurs and large clusters are total free energy functional is given by (1/2)J(∂c/∂x)2, formed. whereJ representsinteratomicinteractions(forexample, the coupling strength in the Ising model). This leads (in dimensionless form) to the ∂2/∂x2[ (J/kT)(∂2u/∂x2)] − 1 term in Eq. (3). The mapping q = 1 exp( J/k T) B − − explains the ln(1 q) coefficient in Eq. (3). Usually, the − mean field form of the local free energy is assumed: 0.95 unstable region 2 4 f(u)=0.5a(u 0.5) +0.25b(u 0.5) . (4) (phase separation) − − 0.9 Figure 2 shows the local free energy both for subcriti- q calandsupercriticaladhesion. Theonlyextremum(minimum) of f(u) for subcritical adhesion is at u = 1/2, so 0.85 the homogeneous state is stable. For supercritical adhesion, the extremum at u = 1/2 becomes a maximum, 0.8 stable region the homogeneous state is no longer stable, and two new (homogenous state) minima appear, with the densities given by Eq. (2). The constants a and b in Eq. (4) are chosen in such a 0.75 0 0.2 0.4 u 0.6 0.8 1 way that the free energy functional satisfies several important conditions. First, we demand that the phase transition threshold is exact (and not its mean field ap- FIG. 1: Phase diagram without proliferation. The critical proximation)asgivenbyEq.(2). Second,astheadhesion adhesionparameterasafunctionofdensityasgivenbyEq.(2) parameter q tends to zero, b should go to zero. Then isshown bysolid line. This curveseparates twoqualitatively different regions. In the stable region, a homogeneous state Eq. (3) transforms into the FK equation. In addition, (with uniformly distributed cells) remains homogeneous; in the diffusion coefficientinEq.(1)shouldbe D¯ =1/4(as contrast, in the unstable region phase separation occurs and canbe derivedfromthe discretemodelwithoutadhesion large clusters are formed. [7]), so that limq−→0a = 1/4. We chose the following expressions: q q c(q) cr a = − , sagTehre[1t1w],o-adnidmtehnesicounravleIsmin(gT)m,owdheilchwsaespasoralvteeds tbhyeOstna-- −|q−qcr|3/4 4qc1r/4 ble and unstable regions, is known (see Fig. 1). In our q qcr 1/4 c(q) b = − , (5) case the threshold is given by the critical adhesion pa- 1 16(1 q)2/q4 q1/4 (cid:20) − − (cid:21) cr 3 x 10−3 The behavior of the density front in the tail region de- 2 pends on the sign of the determinant 1.5 q = qcr P2 4r 3 P3 4Pr Q2 2 q < q D(v)= + + + , cr − 9 3 −27 3 − 2 1 (cid:18) (cid:19) (cid:18) (cid:19) where f 0.5 d2f 1 P = du2|u=0 ln(1 q), 0 (cid:18) (cid:19) − α v r = , and Q= . ln(1 q) ln(1 q) −0.5 − − q > q cr AsintheFKequation,thereisanintervalofpossibleve- −1 locities,v >v . Wecheckednumericallythatforsmall 0 0.2 0.4 u 0.6 0.8 1 min enough α (or small enough q, see below), velocity selection is determined by the condition D = 0. This gives FIG. 2: The local free energy for subcritical (solid line, q = the minimum velocity of front propagation (the critical 0.81), critical (dashed line), and supercritical (dotted line, velocity): q=0.85) adhesion. The minimal adhesion threshold is given dbyen1o6t(e1−twqocrs)t2a/bqlc4erp=ha1s,ews,hsicehegEiqve.s(2qc)r. =0.8284...Twocircles vc2r = 83α dd2uf2|u=0 − 27ln(21 q) dd2uf2|u=0 3 × (cid:18) (cid:19) − (cid:18) (cid:19) d2f −2 3/2 swhhoeurledttheendontloyurensittyricwtihoennoqngtoheesftuonczteiroon. cT(hqi)sifsutnhcattioint 1− 1+12αln(1−q)(cid:18)du2|u=0(cid:19) ! .(8) will be used to fit the continuum results with results of   Asexpected,v tendstozerowhenαgoestozero(fronts discrete simulations. Note that the theoretical analysis cr donotpropagateforzeroproliferation),andv tendsto isperformedforthegeneralformoflocalfreeenergyand cr the value ofthe FK equationwith D¯ =1/4[see Eq.(1)], the specific relations (5) are used only when comparing v √α,whentheadhesionparameterq goestozero. theoretical predictions with numerical simulations. cr −→ The selection rule is illustrated in Fig. 3. If the ve- Inthenextsectionweconsidertheregimeofsubcritical locity is slightly larger than v , Eq. (7) has four real adhesion and focus on front propagation. cr roots: three negative and one positive. As the velocity approaches v , two negative roots approach each other cr and coincide exactly when D = 0, see Fig. 3. It is rea- III. SUBCRITICAL ADHESION: FRONT sonable to suppose, in view of the results for the FK PROPAGATION equation, that this is the selected velocity of front propagation; see [15]. For smaller velocity these two roots A. Theory become complex, which is not allowed as it results in an oscillatory behavior of the density in the tail region. Initially, we put cells into the left part of the channel; Totestthesepredictions,weperformednumericalsim- x measuresthe distance along the channel. In the initial ulations of the Eq. (3). We used the third order Runge- stateallsiteswithx<0areoccupiedandtherestempty. Kuttamethod,withameshsizeδx=1.0,andatimestep For t > 0 cells diffuse and proliferate along the channel δt=0.03. Initial conditions were localized: u=0 ahead and form an advancing front. To analyze those fronts, of the front, u = 1 behind the front, and the interface we look for the solutions in the form u = u(ξ = x vt) had the form u = exp( x)/[1+exp( x)]. Simulations − − − in Eq. (3), where the front velocity v is unknown. This confirm that the front velocity tends to the value given gives by Eq. (8). As in the FK equation, v approaches vcr quite slowly, see Fig. 4. ′′′′ d2f ′′ d3f ′2 ′ Note that Eq. (8) becomes invalid when α (or q) are ln(1−q)u +du2u +du3u +vu +αu(1−u)=0. (6) large enough, so that 1+12αln(1 q)(d2f/du2 u=0)−2 − | becomesnegative. Inthiscasethecharacteristicequation In order to understand velocity selection, we linearize (7)hastworealroots[onenegative(λ1)andonepositive Eq. (6) in the tail region u = 0, in a similar way to the (λ4)]andtwocomplexconjugaterootswithnegativereal analysisofpulledfrontsintheFKequation. Substituting part (λ2,3 = γ iω). Therefore, the density in the tail ± u exp(λξ), we find region is oscillatory and becomes negative, which is for- ∝ bidden. Thisisaninherentfeatureofthetime-dependent 4 d2f 2 Eq.(3). It is knownthat solutions offourth-orderdiffer- E =ln(1−q)λ + du2|u=0 λ +vλ+α=0. (7) ential equations generally do not remainpositive [16]. A (cid:18) (cid:19) 4 u 0whensolvingEq.(3). Inthiscase,thedensitypro- ≥ 0.04 fileisnotanalytic,asinproblemswithcompactsupport: v < v the density becomes zero at some ξcrit and remains zero cr for ξ >ξ . This profile propagates with a well defined crit 0.02 velocity, see Fig. 5. We believe that this is the only type of solution which can be chosen in this regime. E 0 1 2,3 4 B. Comparison with discrete simulations −0.02 Now we compare the results of deterministic contin- v > v uum approach with stochastic discrete modeling. Fig- cr ure 5 shows the front velocity as a function of the adhe- −0.5 0 λ 0.5 sion parameter q for different values of proliferation α. The theoretical predictions are given by Eq. (8). The front velocity in the discrete system is obtained by aver- FIG. 3: The characteristic equation E [Eq. (7)] for different values of front velocity. The solid line corresponds to the aging over many realizations. One can see an excellent critical velocity [given by Eq. (8)], the two dashed lines cor- agreement over wide range of parameters. (Front veloc- respond to larger and smaller velocities. The intersection of ity computed numerically from Eq. (3) also approaches thehorizontaldottedlineE =0withthesolid linegivesfour thesamevalues,seeFig.4). The theoreticalcurvecorre- rootsof thecharacteristic equation (circles). Theparameters sponding to largeα becomes invalid for large q, which is are q=0.4,α=0.003. relatedtotheoscillatorybehaviorofdensitytails. Never- theless,the numericallycalculatedfrontvelocitiesinthis regionarewell-definedandagreewiththosefromdiscrete 0.046 simulations. Thisshowsthatourtheoreticalunderstand- ing in this case is incomplete. Figure 6 shows an example of the corresponding den- 0.044 sity profiles from discrete and continuum simulations. The form of the fronts is very similar, and discrete and v continuumfrontspropagatewiththesamevelocity. Note 0.042 however, that the transient regime for discrete front is longer. 0.04 IV. SUPERCRITICAL ADHESION 0.038 0.5 1 1.5 2 2.5 3 3.5 4 t The situation for q > q is more complicated. As be- x 104 c fore, we start with a sharp front: u = 1 for x < 0 and u = 0 for x > 0. It turns out there is a nontrivial and FIG. 4: Front velocity from numerical solution of Eq. (3) as a function of time. The velocity slowly approaches the long-livedtransientbehavior [1], which we analyze using theoretical value, similarly to the situation in FK equation. discrete and continuum approaches. The parameters are q=0.6,α=0.003. A. Nonzero proliferation similar effect occurs in the extended Fisher-Kolmogorov equation (EFK) [15]. It was shown that in some region We first consider α = 0, q > q . Figure 7 shows a c 6 ofparameterslocalizedinitialconditionscannotdevelop time series of density profile. To obtain the profiles, we into a uniformly translating front solution. Instead, for averaged the density over the channel width and over sufficiently sharp initial conditions one has an envelope many realizations. One realization is shown in Fig. 8. A front: a moving front creates a periodic array of kinks long-lived transient occurs before the propagating front behind it [15]. These oscillations between u = 1 and isformed. Aninsetshowsamagnifiedpictureoftheden- u = 1 occur due to the cubic nonlinearity in the reac- sity profile for early time. An interesting feature is the − tion term of the EFK equation. In our case, u = 1 is secondary density peak. This peak occurs due to phase − not a fixed point, so negative perturbation are not sta- separationandclusterformationinthe low-densityinva- bilized (in fact, numerical simulations of Eq. (3) show sive region. At later times the main front builds up and a finite-time explosion in this region of parameters.) In catchestheisolatedclusters. Thesamefeatureispresent addition, u < 0 has no physical meaning in our case. inthe continuumapproach: Fig.9 showsatime seriesof A possible way to overcome this problem is to demand density profiles from numerical solutions of Eq. (3). 5 1 α = 0.003 0.05 0.8 0.04 0.6 v u 0.03 α = 0.0007 0.4 0.02 α = 0.0003 0.2 0.01 0 0 0.2 0.4 q 0.6 0.8 400 600 x 800 FIG. 5: Front velocity as a function of the adhesion param- FIG.6: Timeseriesofdensityprofilesindiscrete(solidlines) eter q for different values of proliferation α. The theoretical and continuum (dashed lines) models. The parameters are: predictions are given by Eq. (8) (dashed lines). The front q = 0.4, α = 0.003. The discrete fronts correspond to times velocity in the discrete system is obtained by averaging over td1 = 10000 and td2 = 15000, the density profiles computed mforancy(qr)e=ali1za−ti0o.n7.5qT1/h2eicnaltchuelaetxiponresssinioEnqfo.r(8fr)eweeerneeprgeyr,fowrmhiechd f6r0o0m0 annudmtecr2ic=al1s1o0lu0t0i.on of Eq. (3) correspond to times tc1 = gives the best agreement with discrete simulations. One can see that the theoretical curve corresponding to large α becomes invalid for large q, which is related to the oscillatory 1 behaviorofdensitytails. Nevertheless,thenumericallycalcu- 1 lated frontvelocities inthisregion arewell-definedandagree u with those from discrete simulations. Front velocities com- 0.8 puted from the time-dependent Eq. (3) are shown by dia- 0.5 monds. The values of proliferation are α = 0.0003 (circles), t=t 0.6 1 α=0.0007 (asterisks), and α=0.003 (squares). u 0 0 200 x 400 0.4 The transient behavior with a secondary density peak occursduetotheslow(nonexponential)decayofthetail, 0.2 which occurs only for supercritical adhesion. For α> 0, t=t1 t=t2 t=t3 there are two competing processes: the (slow) propagation of the front and relatively fast formation of the sec- 0 0 500 1000 x 1500 2000 ondary peak on the slowly decaying tail. This slowly decaying tail also exists for zero proliferation. In order togainsomeinsightsintoitsformation,wenowstudythe FIG.7: Timeseriesofdensityprofilesindiscretesimulations. relaxation dynamics in Cahn-Hilliard equation [Eq. (3)] An inset shows a magnified picture of the density profile for without proliferation. early time t = t1. The parameters are: q = 0.9, α =0.0003, t1=2.2 104, t2=6.6 104, t3=11 104. × × × B. Zero proliferation thesamedensitytailscoincidewhenmeasuredasafunc- tionof η =x/√t, as expected for purely diffusive behav- Inthissectionwefocusontherelaxationdynamics(for ior. zero proliferation) both in the discrete and continuum Thesamerelaxationdynamicsoccursinthecontinuum models. model. This behavior canbe easily explained. Consider- We start with the discrete lattice model. Initially, all ingthesmall-densityregion(thetail),wecanneglectthe channel sites with x < 0 are occupied (u = 1), and sites ln(1 q)(∂2u/∂x2) term in Eq. (3). Then, substituting withx>0areempty(u=0). However,sinceq >q ,the − c u=u(η) into Eq. (3), we have: final state consists of two phases: a high density phase u = u1(q) for x < 0 and a low density phase u = u2(q) d2f ′′ 1 ′ for x>0, see Eq. (2). du2 + 2ηu =0, (9) Figure 10 shows the tails of two density profiles calcu- (cid:18) (cid:19) ′ lated from the discrete model. As before, we averaged where is the derivative with respect to η. To solve this over the channel width and over many realizations. The equation we need to specify two boundary conditions at relaxation dynamics is self-similar. An inset shows that η =0. Thefirstoneisjustu(η =0)=u2(thelowdensity 6 0 0.06 u 40 200 400 600 800 1000 1200 0.05 0.1 0.04 0 0.05 u 0.03 40 200 400 600 800 1000 1200 0 0 0.5 1 1.5 2 2.5 0.02 η 0 0.01 40 200 400 600 800 1000 1200 0 100 200 300 x 400 500 600 FIG.8: Timeseriesofsnapshots,correspondingtothedensity profiles in Fig. 7. The upperpanel corresponds to t1=2.2 104, the middle panel to t2 =6.6 104, and the lower pan×el FIG. 10: Density profiles (tails) for different times as com- to t3 =11 104. The parameters×are thesame as in Fig. 7. puted from discrete model. An inset shows the same density × tailsasafunctionofη=x/√t. Theparametersareq=0.85, t1=2 104, t2 =105. × 1 1 0.8 u 0.8 0.1 0.6 0.6 u 0.05 u 0.4 0.4 0 0.2 0 1 2 η 3 0.2 0 0 100 200 x 300 400 0 −200 0 x 200 400 FIG.9: Timeseriesofdensityprofilesfromnumericalsolution FIG. 11: Density profiles for different times (blue curve at o(dfaEshq.ed(3l)in.eT),hte1p=ar1a.m8ete1r0s4a,rte2: =q =2.40.9,1α04=, t03.0=0033.0, t01=040, t1 =1.5×104, red curve at t2 =6×104) as computed from and t4 =4.4 104 (soli×d lines). × × Eq.(3). Aninsetshowsthatthetailsofdensityprofilescorre- × spondingtodifferenttimescoincidewhenmeasuredasafunc- tion of η = x/√t. Blue circles correspond to t1, red squares correspond to t2. The asymptotics computed from Eq. (9) ′ is shown by the solid black line. The adhesion parameter is stable phase). Then we find u(η = 0) by a shooting q = 0.85. The simulations were performed for c = 4q1/4 in procedure, demanding u(η ) 0. c −→∞ −→ theexpression for free energy. Figure 11 shows this self-similar relaxation dynamics. Two density profiles depicted in Fig. 11, are calculated from the time-dependent equation (3). An inset shows V. SUMMARY AND DISCUSSION that these profiles (the tails) corresponding to different timescoincidewhenmeasuredasafunctionofη =x/√t. Inthisworkweformulatedandexaminedacontinuum TheasymptoticscomputedfromEq.(9)isinanexcellent model for motile and proliferative cells, which experi- agreement with numerical simulations. ence cell-cell adhesion. We describe collective behavior This self-similarbehaviormeansthatthe decayrate is of the cells by using a modified Cahn-Hilliard equation inversely proportional to t1/2, so it becomes lower with with an additional proliferation term. We identified and time. This is the base of the formation of secondary analyzed different parameter regimes: front propagation densitypeak inthe transientregimefornonzeroprolifer- for subcritical adhesion, nontrivial transient regime for ation. supercriticaladhesionandnonzeroproliferation. There- 7 sults of the continuum descriptionin variousregimesare effect into account in the continuum description. For compared with the results of a discrete model [1]. example, one can consider the following expression for The continuum approach we used is phenomenologi- mobility: M = 1 qu2, or similar forms where M is a − cal. It is presently unknown how to proceed from the decreasing function of density [17]. microscopic lattice model to a macroscopic continuum The modified Cahn-Hilliardequation admits solutions descriptionevenwithoutproliferation. Onewayistouse in the form of propagating fronts. We postulated that the mean field approximation [14], which neglects fluc- the velocity selection procedure is similar to that of the tuations. However, even the statics, namely the phase FK equation and found a critical velocity, Eq. (8). Nu- diagram, Fig. 1, in the mean field approximation is very merical simulations of the time-dependent Eq. (3) con- differentfromthe exactsolution. Wehavetakenanother firmed this result. However, the expression for critical approach, choosing free energy functional that includes velocity becomes invalid in some region of parameters. exact statics in it. This allows us to compare the re- Nevertheless, demanding u 0 we can still solve Eq. (3) ≥ sultsofcontinuumsimulationswiththeresultsofdiscrete numerically and observe propagating fronts with a well- model. defined velocity, see Fig. 5. This problem still needs to The results of continuum theory are in a qualitative be clarified. agreement with discrete simulations in a wide region of An interesting avenue of future work is applying the parameters, in particular for the velocity of front prop- modifiedCahn-Hilliardequationtotheproblemofcluster agation in subcritical regime. The continuum approach nucleation and growth in a two-dimensional system to reproduced a secondary density peak formation in the model clustering of malignant brain tumor cells. transient regime. However, there are important quan- titative differences. Discrete simulations show that the highdensitypartoftheprofilediffusesmuchmoreslowly Acknowledgments than the low density part. At later times this leads to muchsmootherfrontsthaninthecontinuumsimulations. There is no symmetry between particles and holes (oc- We thankP.SmerekaandJ.Lowengrubforuseful dis- cupied and empty sites) in the discrete model [14]. We cussions and the NIH Bioengineering Research Partner- couldintroduceadensity-dependentmobilitytotakethis ship grant, R01 CA085139-01A2for support. [1] E. Khain, L. M. Sander, and C. M. Schneider-Mizell, J. [7] T.Callaghan, E.Khain,L.M.Sander,andR.M.Ziff,J. Stat.Phys. 128, 209 (2007). Stat. Phys. 122, 909 (2006). [2] D. Drasdo, R. Kree, and J. S. McCaskill, Phys. Rev. [8] R. A. Fisher, Annual Eugenics 7, 355 (1937); A. Kol- E 52, 6635 (1995); J. A. Sherratt and J. C. Dallon, mogorov, I.Petrovsky,and N.Piscounov, Moscow Univ. ComptesRendusBiologies 325,557(2002),P.K.Maini, Bull. Math. 1, 1 (1937). D. L. S. McElwain, and D. Leavesley, Appl. Math. Lett. [9] Note however that the front velocity in the discrete sys- 17, 575 (2004). temissmallerthanitscontinuumanalogandapproaches [3] J.A.SherrattandM.A.J.Chaplain, J.Math.Biol. 43, it extremely slow [6, 7]. 291 (2001); L. M. Sander and T. S. Deisboeck, Phys. [10] D.G.AronsonandH.F.Weinberger,Adv.Math.30,33 Rev. E. 66, 051901 (2002); K. R. Swanson, C. Bridge, (1978). J. D. Murray and E. C. Alvord, J. Neurol. Sci. 216, 1 [11] L. Onsager, Phys. Rev. 65, 117 (1944). (2003);H.Hatzikirou,A.Deutsch,C.Schaller,M.Simon, [12] K. Huang, Statistical Mechanics, Wiley, New York and K. Swanson, Math. Mod. Meth. Appl.Sci. 15, 1779 (1987). (2005); C. Athale, Y. Mansury and T. S. Deisboeck, J. [13] J. S. Langer, in Solids far from equilibrium (ed. C. Go- Theor.Biol.233,469(2005);E.KhainandL.M.Sander, dreche) 297–364 Cambridge, NewYork,Cambridge Uni- Phys. Rev. Lett. 96, 188103 (2006); P. Macklin and J. versity Press (1992). Lowengrub, J. Math. Biol. 245, 677 (2007). [14] J. F. Gouyet, M. Plapp, W. Dieterich, and P. Maass, [4] S.Turner,J.A.Sherratt,K.J.Painter,andN.J.Savill, Adv. Phys. 52, 523 (2003), M. Plapp and J. F. Gouyet, Phys. Rev. E 69, 021910 (2004), S. Turner, Phys. Rev. Phys. Rev.E 55, 5321 (1997). E 71, 041903 (2005). [15] W. van Saarloos, Phys. Rev. Lett. 58, 2571 (1987); W. [5] M. Bramson, P. Calderoni, A. Demasi, P. Ferrari, van Saarloos, Phys. Rev. A 37, 211 (1988); G. T. Dee J. Lebowitz, and R. H. Schonmann, J. Stat. Phys. 45, and W. van Saarloos, Phys. Rev. Lett. 60, 2641 (1988); 905(1986),A.R.Kerstein,J.Stat.Phys.45,921(1986). W. van Saarloos, Phys. Rev.A 39, 6367 (1989). [6] E.BrunetandB.Derrida,Phys.Rev.E56,2597(1997); [16] P. Smereka, privatecommunication. D.A.Kessler,Z.Ner,andL.M.Sander,Phys.Rev.E58, [17] J. S. Langer, M. Bar-on, and H. D. Miller, Phys. Rev. 107(1998);E.BrunetandB.Derrida,J.Stat.Phys.103, A. 11, 1417 (1975); K. Kitahara and M. Imada, Prog. 269(2001);E.Moro,Phys.Rev.Lett.87,238303(2001); Theor. Phys. Suppl. 64, 65 (1978); S. Puri, A. J. Bray, D.Panja, G. Tripathy,andW. vanSaarloos, Phys.Rev. and J. L. Lebowitz, Phys. Rev.E 56, 758 (1997). E 67, 046206 (2003); E. Moro, Phys. Rev. E 68, 025102 (R) (2003); D.Panja, Phys. Rep.393, 87 (2004).