# Bifurcation of traveling waves in a Keller-Segel type free boundary model of cell motility

## Abstract

We study a two-dimensional free boundary problem that models motility of eukaryotic cells on substrates. This problem consists of an elliptic equation describing the flow of cytoskeleton gel coupled with a convection-diffusion PDE for the density of myosin motors. The two key properties of this problem are (i) presence of the cross diffusion as in the classical Keller-Segel problem in chemotaxis and (ii) nonlinear nonlocal free boundary condition that involves curvature of the boundary. We establish the bifurcation of the traveling waves from a family of radially symmetric steady states. The traveling waves describe persistent motion without external cues or stimuli which is a signature of cell motility. We also prove existence of non-radial steady states. Existence of both traveling waves and non-radial steady states is established via Leray-Schauder degree theory applied to a Liouville-type equation (which is obtained via a reduction of the original system) in a free boundary setting.

figuresection \counterwithoutfiguresubsection \MakePerPagefootnote

## 1 Introduction

For decades, the persistent motion exhibited by keratocytes on flat surfaces has attracted attention from experimentalists and modelers alike. Cells of this type are found, e.g., in the cornea and their movement is of medical relevance as they are involved in wound healing after eye surgery or injuries. Also, keratocytes are perfect for experiments and modeling since they naturally live on flat surfaces, which allows capturing the main features of their motion by spatially two dimensional models. The typical modes of motion of keratocytes are rest (no movement at all) or steady motion with fixed shape, speed, and direction. That is why the most important solutions will be steady state solutions (corresponding to a resting cell) and traveling wave solutions (a steadily moving cell).

Traveling wave solutions for cell motility models have been investigated both analytically and numerically for free boundary problems in one space dimension, e.g. [36, 37, 2], numerically for free boundary models in two dimensions, e.g. [3, 42], as well as for phase field models, analytically in one dimension, e.g. [5], and numerically in two dimensions, e.g. [40, 47, 39], for an overview we refer to [46, 1] and references therein. In this work we consider a two-dimensional model that can be viewed both as an extension of the analytical one-dimensional model from [36, 37] to 2D and as a simplified version of the computational 2D model from [3]. Our objective is to study the existence of traveling wave solutions for this model. These solutions describe steady motion without external cues or stimuli which is a signature of cell motility.

In [36, 37], the authors introduced a one dimensional model capturing actin (more precisely, filamentous actin or F-actin) flow and contraction due to myosin motors. They proposed a model that consists of a system of an elliptic and a parabolic equation of Keller-Segel type in the free boundary setting. It was shown in [36] that trivial steady states bifurcate to traveling wave solutions. The Keller-Segel system in fixed domains was first introduced and analyzed in [21, 22, 23] and studied by many authors due to its fundamental importance in biology most notably for modeling chemotaxis. There is a vast body of literature on Keller-Segel models with prescribed (fixed rather than free) boundary, see, e.g., [35], [9],[43], [44] review [17] and references therein. The key issue in such problems is the blow up of the solutions depending on the initial data.

In [3] a two-dimensional free boundary model consisting of PDEs for actin flow, myosin density and, additionally, a reaction-diffusion equation for the cell-substrate adhesion strength was introduced based on mechanical principles. Simulations of this model reveal steady state and traveling wave type solutions in two-dimensions that are compared to experimental observations of keratocyte motion on the flat surfaces. The steady state solutions are characterized by a high adhesion strength (high traction) whereas the moving cell solutions correspond to a low overall adhesion strength. In both cases, the adhesion strength is spatially almost homogeneous. Therefore in this work we consider a simplified two-dimensional problem with constant adhesion strength parameter similar to the one dimensional model of [36, 37]. We further simplify the model in [3], see also review [38] by considering a reduced rheology of the cytoskeleton based on the high contrast in numerical values for shear and bulk viscosities cited in [3]. Thus following [32] we consider equations from [3] with shear viscosity and bulk viscosity scaled to .

The main building block of the model considered in this work is a coupled Keller-Segel type system of two partial differential equations. The first one (obtained after the above simplification of equation from [3]) in dimension-free variables writes as follows:

(1.1) |

where is the time dependent domain in occupied by the cell, is the velocity of the actin gel, and is the myosin density. This equation represents the force balance between the stress in the actin gel on the left hand side and the friction (proportional to the velocity) between the cell and the substrate on the right hand side. Since the shear viscosity , the stress is a scalar composed of a hydrodynamic (passive) part and the active contribution generated by myosin motors. Identifying with the corresponding scalar matrix (tensor ), equation (1.1) can be rewritten in the standard form . Equation (1.1) is coupled to an advection-diffusion equation for the myosin density :

(1.2) |

Myosin motors are transported with the actin flow if bound to actin and freely diffuse otherwise, reflected by the second and first term on the right hand side of (1.2), respectively. Assuming that the time scale for binding and unbinding is very short compared to those relevant for our problem, the densities of bound and unbound myosin motors can be combined into the effective density (see e.g. [36, 37]).

Following [3], the evolution of the free boundary is described by the kinematic boundary condition for the normal velocity ,

(1.3) |

where is the unit outward normal, stands for the curvature of , and constant defined by enforces area preservation. The kinematic condition (1.3) equates the normal velocity of the boundary to the contributions from the normal component of the actin velocity, the surface tension of the membrane ( being the curvature), and the area preservation term . The latter term is constant along the boundary and is interpreted as actin polymerization at the membrane, it compensates for the difference between velocities of the actin gel and the membrane.

On the boundary, equation (1.1) is supplied with the zero stress condition

(1.4) |

whereas for the equation (1.2), a no-flux condition is assumed:

(1.5) |

Similar parabolic-elliptic free boundary problems frequently occur in modeling of biological and physical phenomena. One type of problem arises in tumor growth models, e.g. [13, 11, 16, 18] (see also reviews [12, 30]), however, these are typically linear problems, and the domain area is not preserved. For these models, steady state solutions have been described, and bifurcations to different steady states or growing/shrinking domain solutions have been investigated. Another type of problem arises in the modeling of wound healing, see, e.g., [19], where a free boundary problem for a reaction diffusion equation is used to model the evolution of complex wound morphologies. These models are often agent based rather than continuum models, see, e.g., [7]. More recently, mechanical tumor models have been devised leading to Hele-Shaw type problems, e.g. [34].

In the above works the focus is on solutions describing motion with constant velocity in domains that expand or contract rather than domains of fixed size and shape moving with constant velocity. Besides this shift of focus, the main novelty of the free boundary problem under consideration is the cross diffusion term in equation (1.2) giving rise to the Keller-Segel structure of the bulk equations. This structure was introduced in one dimensional models of cell motility in [36, 37]. While the boundary in one dimensional models (e.g. [36, 37]) consists of just two points, in two dimensional free boundary models the shape of the domain is unknown. This poses questions that do not arise in one dimensional settings and leads to novel challenges in analysis, for example, bifurcations from radially symmetric to non-radially symmetric shapes.

We are interested in traveling wave solutions of (1.1) - (1.3), i.e. solutions of the form , , . Thus after passing to the moving frame and rewriting system (1.1)-(1.5) in terms of the scalar stress we are led to the following free boundary problem

(1.6) |

(1.7) |

(1.8) |

We now outline the main result of the paper (see, Section 6 for further details) and key ingredients of the proof.

###### Theorem 1.1.

Without loss of generality we assume motion in -direction and, slightly abusing notation, write . Furthermore, for a given all nonnegative solutions of (1.7) ( represents the density of myosin and therefore cannot be negative) are given by with some constant . Indeed, it is straightforward that is a solution of (1.7). The uniqueness up to a multiplicative constant follows from the Krein-Rutman theorem [26], or alternatively using the factorization , considering as a unknown function, and proving that by showing that it satisfies an advection-diffusion equation with zero Neumann condition. This allows us to eliminate from (1.6)-(1.7) and rewrite the problem of finding traveling waves in the following concise form:

(1.9) |

with boundary conditions

(1.10) |

and

(1.11) |

Note that an ODE similar to the PDE (1.9) was obtained in the analysis of the one dimensional free boundary problem for the Keller-Segel type system in [36, 37]. In problem (1.9)-(1.11) , , and are unknowns and the parameter is given. Note that (1.9)-(1.11) is a free boundary problem, that is, the domain is also unknown. For radially symmetric solutions of (1.9)-(1.10) with and being a disk, the constant can always be chosen so that the boundary condition (1.11) is satisfied. This observation allows us to construct a one-parameter family of radially symmetric steady state solutions by solving the nonlinear eigenvalue problem (1.9)-(1.10). Furthermore, the equation (1.9) contains exponential nonlinearity, as in the classical Liouville equation [29] which has explicit radially symmetric solutions, but the additional zero order term in the left hand side of (1.9) complicates the analysis. Note that non-trivial steady states also exist in the one-dimensional case [36, 37] (they are unstable).

We rely on an argument from [10] (see also [25]) based on the Implicit Function Theorem to show existence of an analytic curve of radially symmetric solutions of (1.9)-(1.10). Moreover these solutions are extended to the case of nonzero in (1.9) and small perturbations of the domain from a given disk. Then (1.9)-(1.11) is reduced to selecting solutions of (1.9)-(1.10) that satisfy (1.11). Considering the linear part of perturbations of radially symmetric solutions we (formally) derive the condition (3.8) (Section 3) for a bifurcation from the steady states to genuine traveling waves (with ). We next show that the condition (3.8) is indeed satisfied on a nontrivial radially symmetric steady state solution belonging to , exploiting a subtle bound on the second eigenvalue of the linearized problem for the Liouville equation from [41]. Yet another technically involved part of this work is devoted to recasting (1.9)-(1.11) as a fixed point problem in an appropriate functional setting. Then a topological argument based on Leray-Schauder degree theory rigorously justifies the existence of traveling waves with . Both the recasting and the topological argument require spectral analysis of various linearized operators appearing in these considerations. Next the techniques developed for establishing traveling waves solutions are also used to find steady states with no radial symmetry.

The rest of paper is organized as follows. In Section 2 we find a one parameter family of radially symmetric steady state solutions and establish their properties. In Section 3 we derive a necessary condition (3.8) for the bifurcating from the family of radially symmetric steady states to a family of traveling wave solutions () and show that this condition is satisfied on the analytic curve of radially symmetric solutions. In Section 4 we investigate the spectral properties of the linearized operator of the equation (1.9) around radially symmetric steady states. This operator appears in a number of the subsequent constructions. In section 5 we establish existence of the solutions to the Dirichlet problem (1.9)-(1.10) and study their properties. This is done for small but not necessarily zero velocity in a prescribed domain , which is a perturbation of a disk. Section 6 completes the proof of the main result on the bifurcation of the steady states to traveling waves. To this end we rewrite (1.9)-(1.11) as a fixed point problem, and study the local Leray-Schauder index of the corresponding mapping. We show that this index jumps at the potential bifurcation point (identified in Section 3). This establishes the bifurcation at this point. Finally, in section 7 we prove existence of nonradial steady states. In the Appendix A we construct three terms of the asymptotic expansion of traveling wave solutions in powers of small velocity, which allow us to describe the emergence of non-symmetric shapes both analytically and numerically.

## 2 Family of radially symmetric steady states

Problem (1.9)-(1.11) has a family of steady solutions, with , found in a radially symmetric form. Namely, let be a disk with radius , then we seek radially symmetric solutions , , of the equation

(2.1) |

with boundary conditions

(2.2) |

Note that (2.1)-(2.2) is a nonlinear eigenvalue problem, i.e. both the constant and the function are unknowns in this problem. Every solution of (2.1)-(2.2) also satisfies (1.9)-(1.11) with and some constant , that is we can always choose in this radially symmetric problem, so that the condition (1.11) is satisfied. Equation (2.1) is the classical Liouville equation [29] with an additional zero order term (the second term on the left hand side of (2.1)). Various forms of the Liouville equation arise in many applications ranging from the geometric problem of prescribed Gaussian curvature to the relativistic Chern-Simons-Higgs model [33], the mean field limit of point vortices of Euler flow [8] and the Keller-Segel model of chemotaxis [45]. For a review of the literature on Liouville type equations we address the reader to [28] and references therein. While the above works mostly address the issues related to the blow-up in the Liouville equation, see e.g., [27], in contrast our focus is on the construction of the family of solutions and its properties. Since we are concerned with special solutions of (1.1)-(1.5) such as traveling waves and steady states rather than general properties of this evolution problem, the issue of blow-up does not arise.

The following theorem establishes existence of solutions of problem (2.1)-(2.2), and the subsequent lemma lists some of their properties.

###### Remark 2.1.

It is natural to expect that the set of solutions of (2.1)-(2.2) has the same structure as the explicit solutions of the classical Liouville equation [41] in the disk. However, the presence of the additional term in (2.1) complicates the analysis even in the radially symmetric case, in particular, the standard trick based of Pohozhaev identity no longer can be used to establish non-degeneracy (see condition (2.8)).

###### Theorem 2.2.

Fix , then

(i) There exists a continuum (a closed connected set) of nonnegative solutions , of (2.1)-(2.2), emanating from the trivial solution . There exists a finite positive

in particular, for all . On the other hand is not bounded in , and moreover

(2.3) |

(ii) For every there exists a pointwise minimal solution (solution which takes minimal values at every point among all solutions) of (2.1)-(2.2), and these minimal solutions are pointwise increasing in . They form an analytic curve in which can be extended to an analytic curve . The curve is the connected component of that contains , where

(2.4) |

and denotes the second eigenvalue of the linearized eigenvalue problem

(2.5) |

###### Remark 2.3.

Summarizing part of the theorem we have the following inclusions

continuum of solutions | minimal solutions |

where at most may be disconnected. The theorem establishes existence of the analytic curve of radial solutions along which bifurcations to traveling waves with nonzero velocity occur (see Lemma 3.1).

###### Proof.

(i) By the maximum principle every solution of (2.1)-(2.2) with is positive for . Let denote the first eigenvalue of in with homogeneous Dirichlet boundary condition, and let be the corresponding eigenfunction. Then multiplying (2.1) by and integrating we find

and therefore .

To show the existence of the continuum , we rewrite (2.1) as

(2.6) |

with , , and the new unknown parameter in place of . Then we resolve (2.6) with Dirichlet condition on , considering the right hand side of (2.6) as a given function. This leads to an equivalent reformulation of (2.1)-(2.2) as a fixed point problem of the form

(2.7) |

By standard elliptic estimates is a compact mapping in , moreover is a bounded subset of . Therefore we can apply Leray-Schauder continuation arguments, see, e.g., [31], and find that there is a continuum of solutions of (2.7) emanating from and takes all nonnegative values. Then in view of the boundedness of we conclude that . This in turn implies (2.3) by Corollary 6 of [6].

(ii) According to [20] there is a minimal solution of (2.1)-(2.2) for each with depending monotonically on . Consider now any, not necessarily minimal, solution such that the second eigenvalue of the linearized problem (2.5) is positive. By using well-etablished techniques based on the Implicit Function Theorem, see, e.g. [25], we obtain that all the solutions of (2.1)-(2.2) in a neighborhood of belong to a smooth curve through , provided that either the linearized problem (2.5) has no zero eigenvalue or this eigenvalue is simple and the corresponding eigenfunction satisfies the non-degeneracy condition

(2.8) |

Since by assumption , the zero eigenvalue, if any exists, is the first eigenvalue of (2.5) and therefore has a fixed sign and necessarily (2.8) holds. Thus is indeed a smooth curve, it contains the minimal solutions (those, for which the first eigenvalue of linearized problem (2.5) is nonnegative) but extends beyond these. Finally, since the nonlinearity in (2.1) is analytic the curve is analytic as well, see the proof of Proposition (5.1). ∎

###### Lemma 2.4.

###### Proof.

To show (2.9) we first prove that is decreasing. Assume to the contrary that takes a local minimum at and there is such that . Multiply (2.1) by and integrate from to to get

(2.11) |

On the other hand, the left hand side of (2.11) is

Therefore is constant on , this in turn implies that is constant on , a contradiction. Thus for . Next, assuming that at a point we get . This also implies that is constant on . Finally, by the Hopf Lemma.

## 3 Necessary condition for bifurcation of traveling waves

We seek traveling wave solutions with small velocity, i.e. solutions of (1.9)-(1.11) for small , as perturbations of radially symmetric steady states given by a pair of solutions to (2.1)-(2.2). To this end we plug the ansatz

(3.1) |

into (1.9)-(1.11). The function describes the deviation of from the disc while describes the deviation of the stress from . Note that in this first order approximation the constant is not perturbed (see Appendix A, where it is shown that the first correction ). Equating like powers of , the terms of order in (1.9) yields the linear inhomogeneous equation for :

(3.2) |

Furthermore, equating terms of the order in the boundary conditions (1.10), (1.11) we get

(3.3) |

and

(3.4) |

To get rid of trivial solutions arising from infinitesimal shifts of the disk , we require to satisfy the orthogonality condition

(3.5) |

A solution of (3.2)-(3.3) is sought in the form of the Fourier component . Then, has to satisfy

(3.6) |

and, owing to (3.5) and (3.3), the boundary condition

(3.7) |

Now multiply (3.6) by and integrate from to . Taking into account that differentiating (2.1) yields , we integrate by parts to obtain

(3.8) |

where we have also used (3.7) and (3.4). This is a necessary condition for existence of traveling waves bifurcating from the steady state curve at the point , and it can be equivalently rewritten using (2.1)-(2.2),

(3.9) |

or, using (2.10),

(3.10) |

The following Lemma shows that there exists a pair satisfying (3.8), and subsequent Corollary 3.2 specifies such a pair which is used in the proof of Theorem 1.1.

###### Lemma 3.1.

###### Proof.

Let us consider minimal solutions in corresponding to small , and small . We show that the left hand side of (3.9) is strictly less than its right hand side by considering the leading term of the asymptotic expansion of solutions in the limit . Linearizing (2.1) -(2.2) about we get

(3.13) |

By the maximum principle for , and therefore on the left hand side of (3.9) we have

for some independent of , while on the right hand side of (3.9),

Next we show existence of satisfying (3.12).

Case 1: . According to items (i) and (ii) of Theorem 2.2, the curve satisfies

(3.14) |

or, if this is false, at least

(3.15) |

If (3.14) holds then right hand side (3.9) becomes negative, while the left hand side is positive, and we are done.

Now consider the case that (3.15) holds. By continuity of there is a pair such that the second eigenvalue of (2.5) is less than . In other words, the second eigenvalue of

(3.16) |

is negative. Then, according to Proposition 2 in [41], we have

Assume by contradiction, that the right hand side of (3.9) is bigger than or equal to its left hand side, then in view of the equivalent reformulation (3.10) of (3.9), we find

(3.17) |

which in turn implies that

(3.18) |

On the other hand, multiplying (2.1) by and integrating we find

(3.19) |

Combining (3.19) with (3.17) and the first inequality in (3.18) we get

(3.20) |

Finally, applying the Cauchy-Schwarz inequality and using the second inequality in (3.18) leads to

(3.21) |

Thus, (3.20) and (3.21) yield the lower bound for the radius, , so that the Lemma is proved for .

Case 2: . Observe that the maximal value of admits the lower bound . Indeed, considering the initial value problem

(3.22) |

we find that continuously varies from to as decreases from to . Therefore there exists some such that is a solution of (2.1)-(2.2). Now consider the minimal solution of (2.1)-(2.2) with and introduce the function solving the auxiliary problem

(3.23) |

Since is a positive subsolution of (2.1)-(2.2), we have

Therefore, in order to prove the inequality

(3.24) |

it suffices to show that

(3.25) |

The solution of (3.23) is explicitly given by

where , and is the modified Bessel function of the first kind. Since

is increasing in and

the inequality (3.25) holds for , and so does (3.24). This completes the proof of Lemma 3.1. ∎

###### Corollary 3.2.

###### Proof.

The result follows from Lemma 3.1 thanks to analyticity of the curve and to the fact that is connected. ∎

## 4 Fourier analysis of the linearized operator

To construct solutions of problem (1.9)-(1.10) as perturbations of radially symmetric steady states we need to study the properties of the linearized operator of this problem. Namely, we consider the linearized spectral problem

(4.1) |

###### Proposition 4.1.

For any , the eigenvalue corresponding to the Fourier modes and ,

(4.2) |

is positive, .

###### Proof.

For each and any solution of (2.1)-(2.2), the function is strictly positive and satisfies (by differentiating (2.1))

(4.3) |

or, for any given ,

(4.4) |

Multiplying (4.2) by and integrating from to yields

(4.5) |

where we introduced the abbreviation

We represent as and multiplying (4.4) by , integrate from to . Integrating by parts in the first term we get

(4.6) |

Subtracting (4.6) from (4.5), we find

(4.7) |

Now pass to the limit in this equality as . Observing that the as of the last term in (4.7) is nonnegative we obtain that and if , then , where is a constant. In the latter case , contradiction. Thus . ∎

###### Corollary 4.2.

For each satisfying

the problem

(4.8) |

has a solution. Moreover precisely one such a solution is orthogonal in to all radially symmetric functions .

###### Proof.

Introduce the solution of

(4.9) |

and observe that . Then a solution of the problem

is obtained by separation of variables and applying Proposition 4.1. ∎

## 5 Existence of solutions of the problem (1.9)-(1.10)

For a given we consider a fixed steady state . Using well-established techniques based on the Implicit Function Theorem, see, e.g., Chapter I in [25], we construct a family of solutions of (1.9)-(1.10) in domains given by

(5.1) |

with sufficiently small , , and with small, but not necessarily zero, velocity . Hereafter, slightly abusing the notation, we identify the angle with the corresponding point on the unit circle .

In order to reduce the construction to a fixed domain we introduce the mapping defined in polar coordinates by

(5.2) |

where is such that when and when . Clearly, (5.2) defines a -diffeomorphism whenever is sufficiently small together with its first and second derivatives.

Among all perturbations we single out those satisfying the area preservation condition

(5.3) |

or in linear approximation

The following proposition establishes existence of solutions of problem (1.9)-(1.10). These solutions are obtained as perturbations of the radially symmetric steady states from Section 2.

###### Proposition 5.1.

There exists some such that for all in -neighborhood of the problem (1.9)-(1.10) admits a solution , in the domain (given by (5.1)). Here is an auxiliary real parameter (to be specified in the proof) such that

(5.4) |

defines an analytic parametrization of the curve in a neighborhood of . Moreover, the mappings

belong to and , respectively. The derivatives and at are given by

(5.5) |

where is a unique, as in Corollary 4.2, solution of

(5.6) |

with , and . The derivatives and at satisfy

(5.7) |

for such that , where is a unique, as in Corollary 4.2, solution of of the problem

(5.8) |