Modelling adhesion-independent cell migration

TL;DR

This paper introduces a mathematical model for cell migration that occurs without adhesion, using an elastic cortex and cortical flow, and demonstrates its ability to replicate experimental behaviors.

Contribution

It presents a novel nonlinear degenerate parabolic model for adhesion-independent cell motility and provides an existence analysis and numerical simulations.

Findings

Model captures key features of adhesion-independent migration

Simulations align with experimental observations

Suggests a mechanical mechanism for cell motility without adhesion

Abstract

A two-dimensional mathematical model for cells migrating without adhesion capabilities is presented and analyzed. Cells are represented by their cortex, which is modelled as an elastic curve, subject to an internal pressure force. Net polymerization or depolymerization in the cortex is modelled via local addition or removal of material, driving a cortical flow. The model takes the form of a fully nonlinear degenerate parabolic system. An existence analysis is carried out by adapting ideas from the theory of gradient flows. Numerical simulations show that these simple rules can account for the behavior observed in experiments, suggesting a possible mechanical mechanism for adhesion-independent motility.

Tables

Table 1. Table 1: Numerical and model parameters for the simulations of the paper.

Parameters
Numerical parameters
Symbol	Numerical Value		Description
$\Delta t$	${4.10}^{-2}$		Time step
$\Delta s$	${5.10}^{-3}$		Discretization step
$N$	200		Number of nodes
$ϵ$	0.1		Distance to obstacle for projection
$a$	adapted		Length of the polymerization zone
Model parameters
Symbol	Numerical Value	Biological value	Description
$L$	1	$48\mu m$	Reference membrane length
$L_{\text{eff}}$	2	94 $\mu m$	Effective membrane length
${\kappa }_S$	1	7 $pN.\mu m^{-1}$	Stiffness constant
$p$	$3.2$	22 $pN.\mu m^{-1}$	Pressure force
$v$	$2$	12 $\mu m.mi{n}^{-1}$	Polymerization speed
$\mu$	1	55 $pN.\mu m^{-1}.min$	Internal friction coefficient

Equations

$$\Gamma (t)=\left\{X(s,t):s\in {\mathbb{T}}^1\right\}\subset {\mathbb{R}}^2.$$

$$\tau =\frac{{\partial }_{s}X}{|{\partial }_{s}X|},n=-{\tau }^{\perp },$$

$$E_{\mathrm{el}}(X)=\frac{1}{2}{\int }_{{\mathbb{T}}^1}{\left(|{\partial }_{s}X|-1\right)}_{+}^{2}ds.$$

$$E_{\text{p}}(X)=-p\left|\Omega (t)\right|=\frac{p}{2}{\int }_{{\mathbb{T}}^1}X\cdot {\partial }_s{X}^{\perp }ds.$$

$$\Omega (t)\subset {\Omega }_c,$$

$$E_{\text{obst},\delta }(X)={\int }_{{\mathbb{T}}^1}{W}_{\delta }(X)ds,W_{\delta }={\left(\delta +{\rho }_{\delta }\ast \text{mathds}{1}_{{\Omega }_c}\right)}^{-1},$$

$$l(s,t)={\int }_0^s|{\partial }_{s}X(\sigma ,t)|d\sigma .$$

$$s(l,t)={\int }_0^l\rho (\lambda ,t)d\lambda ,$$

$${\partial }_t\rho =f.$$

$${\int }_{\Gamma }fdl={\int }_{{\mathbb{T}}^1}f(l(s,t))|{\partial }_{s}X(s,t)|ds=0,t\ge 0.$$

$$\frac{D}{Dt}={\partial }_t+\left({\int }_0^{l(s,t)}f(\lambda ,t)d\lambda \right){\partial }_s,$$

$$g(s,t)=f(l(s,t),t)|{\partial }_{s}X(s,t)|,$$

$$\frac{D}{Dt}={\partial }_t+\left({\int }_0^{s}g(\sigma ,t)d\sigma \right){\partial }_s,$$

$$\frac{DX}{Dt}$$

$$\frac{d}{dt}{\int }_{{\mathbb{T}}^1}Xds=-{\int }_{{\mathbb{T}}^1}\left({\partial }_{s}X{\int }_0^{s}g(\sigma ,t)d\sigma \right)ds={\int }_{{\mathbb{T}}^1}Xgds.$$

$${\partial }_{t}X+{\partial }_{s}X{\int }_0^{s}gd\sigma ={\partial }_s\left((|{\partial }_{s}X|-1{)}_{+}\frac{{\partial }_{s}X}{|{\partial }_{s}X|}\right)-p{\partial }_s{X}^{\perp }-\nabla W_{\delta }(X)+{\mathcal{F}}_{\mathrm{comp}}.$$

$${\int }_{{\mathbb{T}}^1}{\mathcal{F}}_{\mathrm{comp}}ds=-{\int }_{{\mathbb{T}}^1}Xgds.$$

$${\partial }_s{S}_{\mathrm{flow}}[X]={\partial }_{s}X{\int }_0^{s}gd\sigma -{\mathcal{F}}_{\mathrm{comp}}.$$

$${\partial }_{t}X+{\partial }_s(S_{\mathrm{flow}}[X])={\partial }_s\left((|{\partial }_{s}X|-1{)}_{+}\frac{{\partial }_{s}X}{|{\partial }_{s}X|}\right)-p{\partial }_s{X}^{\perp }-\nabla W_{\delta }(X),$$

$$X_{\text{equ}}(s)=\frac{1}{2\pi-p}\bigl{(}\cos(2\pi s),\sin(2\pi s)\bigr{)}\,.$$

$$X_{\delta }\in H_{\mathrm{loc}}^1({\mathbb{R}}_{+};\mathcal{L})\cap L_{\mathrm{loc}}^{\infty }({\mathbb{R}}_{+};\mathcal{H}),$$

$$X\in H_{\mathrm{loc}}^1({\mathbb{R}}_{+};\mathcal{L})\cap L_{\mathrm{loc}}^{\infty }({\mathbb{R}}_{+};\mathcal{H}),$$

$${\partial }_{t}X+{\partial }_s(S_{\mathrm{flow}}[X])={\partial }_s\left((|{\partial }_{s}X|-1{)}_{+}\frac{{\partial }_{s}X}{|{\partial }_{s}X|}\right)-p{\partial }_s{X}^{\perp },$$

$${⟨{\partial }_{t}X+{\partial }_s(S_{\mathrm{flow}}[X])-{\partial }_s\left((|{\partial }_{s}X|-1{)}_{+}\frac{{\partial }_{s}X}{|{\partial }_{s}X|}\right)+p{\partial }_s{X}^{\perp },Y-X⟩}_{\mathcal{L}}\ge 0,$$

$${\partial }_{t}X(s,t)={\mathbf{P}}_c\left({\partial }_s\left((|{\partial }_{s}X|-1{)}_{+}\frac{{\partial }_{s}X}{|{\partial }_{s}X|}\right)-p{\partial }_s{X}^{\perp }-{\partial }_s(S_{\mathrm{flow}}[X])\right),$$

$${\mathbf{P}}_c\left(F(s,t)\right)=\{\begin{array}{ll}(F(s,t)\cdot {\tau }_c){\tau }_c& \text{ for }X\in \partial {\Omega }_c,F(s,t)\cdot n_c>0,\\ F(s,t)& \text{ otherwise},\end{array}$$

$$E_{\delta }(X)$$

$$E_{\delta }(X)\ge \frac{1}{2}\left(1-\frac{p}{2\pi }\right)\Vert {\partial }_{s}X{\Vert }_2^2-\Vert {\partial }_{s}X{\Vert }_2\ge -\frac{\pi }{2\pi -p}.$$

$$E_{\mathrm{p}}(X)$$

$$E_{\mathrm{el}}(X)$$

Full text

1-RICAM-Österreichische Akademie der Wissenschaften, Postgasse 7-9, 1010, Wien, Austria

2-INRIA, LJLL, UMPC, 4, place Jussieu, Couloir 16-26, 3e étage 75252 Paris Cedex 05, France

3-Institute of Science and Technology Austria, Am Campus 1, 3400 Klosterneuburg, Austria

4-Faculty of Mathematics, University of Vienna, Oskar-Morgenstern Platz 1, 1090 Vienna, Austria

Keywords: Variational methods; weak solutions; cell motility modelling; cellular cortex; actin polymerization

AMS Subject Classification: 35K40, 35K51, 35Q92, 35A15, 2C17

1 Introduction

One of the most important cellular behaviors is crawling migration. It is observed in many cellular systems both in culture and in vivo [13, 14], and involved in many essential physiological or pathological processes (wound healing, embryonic development, cancer metastasis etc. [22]). Since the works of Abercrombie in 1970 [1], which described the multistep model of lamellipodia-based cell migration, numerous authors have studied the mechanisms of actin-based cell migration. As a result, despite some remaining open questions, lamellipodial migration is now well understood and described[12]. This migration mode implies that specific adhesion points transmit intracellular pulling forces from the cytoskeleton to the substrate [15, 20]. Actin filaments polymerize below the leading plasma membrane generating pushing forces, and plasma membrane tension resists actin network expansion, pushing back the actin filaments into the cell body. Through adhesion complexes linking the cytoskeleton to the substrate, these retrograde forces are translated into forward locomotion of the cell body [15].

Yet, recent studies indicate that cell migration can be achieved without adhesion in confining three-dimensional environments [2]. Increasing levels of confinement seem to favor adhesion-independent migration in many cell types [5], and can trigger transitions from adhesion-based towards low-adhesive migration modes. However, too strong confinement decreases and even prevents migration [23], due to cell stiffness and nucleus volume. If in the last decade adhesion-independent migration has emerged as a possibly common migration mode, the mechanisms of cell propulsion in this case are still poorly understood. So far in the literature, there is only one known alternative to lamellipodial migration: membrane blebs [3, 4]. These blebs are cellular extensions free of actin filaments, and generated by intracellular hydrostatic pressure [18]. Once generated, these blebs grow until a new actin cortex is reassembled inside, which eventually contract and allow the cell to move [24]. However, numerous studies display cells migrating in an adhesion-free manner without bleb formation.

Several physical mechanisms have been proposed for force transmission between cell and substrate during migration without focal adhesions: (i) cell migration by swimming (by creation of blebs, [24]), (ii) force transmission based on cell-substrate intercalations of lateral protrusions into gaps in the matrix, (iii) chimneying force transmission where cells push against the obstacles [7, 10] or again (iv) flow-friction driven force transmission. In this last case, mechanisms based on non-specific friction between the cell and the substrate have been investigated to account for adhesion-independent migration [7]. Here, intracellular forces generated by the cytoskeleton are transmitted to the substrate via non-specific friction which has been experimentally measured in [2, 9]. The molecular origin of nonspecific friction has not been experimentally investigated. Friction could result from interactions between molecules at the cell surface and the substrate, and unveiling the microscopic origin of nonspecific friction will be an important question for future studies.

In this paper, we propose a simplified 2D model for focal adhesion-independent cell migration, based on the mechanisms (iv). We aim to develop a simplified framework to study whether adhesion-free migration could be primarily driven by simple mechanical features.

Our model is focused on the cell’s cortex, which is assumed to be an elastic material confined in the horizontal plane. Because of the cytoplasmic pressure, it is subject to outwards pressure forces. Mathematically, this takes the form of a system of parabolic equations, where the polymerization process leads to an advection-type term. This continuous formulation also allows to properly define the reaction forces compensating mass displacement due to membrane renewal. With this very simple model, we are able to trigger cell migration, with speed depending on the geometrical characteristics of the obstacles. Our results seem to be in qualitative agreement with the biological observations. In Section 2, we present the biological observation of leukocyte migrating cells and the experimental setting. Section 3 is concerned with the mathematical model which is analyzed in Section 4. Section 5 presents our numerical results.

2 Leukocyte migration in artificial microchannels

To address the question of adhesion-independent migration of leukocytes, we took advantage of a specific line of lymphocytes that allows for genetic modifications. This enabled us to block the synthesis of talin, an adaptor protein essential for adhesion functionality, using a genetic engineering technique known as CRISPR/Cas9 [17]. In addition, we used well-established microfabricated channels [19] to mimic the confined in vivo environment, coupled to a home-made microfluidic set-up [16]. The top and bottom walls of these channels are flat, and the side walls can have various structures. We observe that cells, in which talin has been knocked out, are completely unable to adhere and migrate in channels with flat side walls [16]. Strikingly, their motility is restored in channels with structured side walls with wave length of the wall structure on the order of magnitude of the cell diameter (see Figure 1). This indicates that adhesion-free motility relies on a structured confinement.

3 The mathematical model

Since the mechanisms producing the behavior described in the previous section are not known, we propose a rather simple model for the essential components. The essential idea is that, guided by a chemotactic signal, the cell polarizes with increased actin polymerization near the front end. This is assumed to induce a flow of the cell cortex from front to rear, where depolymerization dominates. The cortex is assumed to be an elastic material with a tendency to equidistribute actin along the cell periphery. This mechanism, together with a constant cytoplasmic excess pressure (actually the pressure difference between cytoplasmic and extracellular pressure), determines the cell shape.

The most critical model ingredient are the forces between the cell and its environment. We assume an unspecific friction with the extracellular liquid, assumed at rest, counteracted by a compensating force, which can be seen as a consequence of the intracellular transport of actin from rear to front. This compensation is chosen such that the cell does not move in an unconstrained environment.

Motivated by the experimental setup we choose a two-dimensional model, which seems reasonable for the experimental situation, where the cell is confined between two flat surfaces. Figure 2 illustrates a discrete version of the model, where the cell cortex is described by mass points. We return to this description in Section 5 for simulation purposes, but here we shall formulate a continuous version of the model, where at time $t\geq 0$ the cortex is represented by a Jordan curve

[TABLE]

The one-dimensional torus $\mathbb{T}^{1}$ will be represented by the unit interval, and the variable $s$ corresponds to the amount of actin material along the cortex, i.e. for some non empty interval $[s_{1},s_{2}]$, the amount of actin on the corresponding piece of cortex is $s_{2}-s_{1}$, this will be formalized later on. The total amount is normalized to $1$. The interior of $\Gamma(t)$ is denoted by $\Omega(t)$, such that $\Gamma(t)=\partial\Omega(t)$. Assuming that $\Gamma(t)$ is smooth enough, we denote by $\tau(t,s)$ and $n(t,s)$ the unit tangent and unit outward normal vectors. Assuming positive orientation of the parametrization, we have

[TABLE]

with the convention $(a,b)^{\bot}=(-b,a)$. The notation is illustrated in Figure 3.

The cortex is assumed to be elastic and in equilibrium if $|\partial_{s}X|=1$, such that $1$ represents the scaled total equilibrium length of the cortex. An elastic resistance against stretching, but not against compression, is described by the potential energy functional

[TABLE]

Neglecting resistance against compression can be seen as a convexification of the elastic energy, which facilitates the analysis of Section 4. Actually, we expect the cortex to be always under tension, such that this assumption should not be relevant from a modelling point of view.

This expectation relies on another model ingredient, a cytoplasmic pressure exceeding the extracellular pressure by a constant amount $p>0$. The associated potential energy contribution is given by

[TABLE]

The assumption of a prescribed constant value of $p$ can be seen as a model simplification. The volume of the cell is mainly dictated by the amount of water it contains, which is subject to osmosis, which can be neglected on the time scales at play here. It would then be closer to reality to assume a fixed prescribed cell volume, measured by the area $|\Omega(t)|$. In this case $p$ would become a time dependent unknown with the mathematical interpretation of a Lagrange multiplier.

The space constraints, i.e. the channel walls, are modelled by requiring

[TABLE]

where $\Omega_{c}\subset\mathbb{R}^{2}$ denotes the inside of the channel. For analysis purposes the obstacles will be softened by introducing the energy contribution

[TABLE]

where $\mathds{1}_{\Omega_{c}}$ denotes the indicator function of the admissible region, the small parameter $\delta>0$ measures the softness of the obstacle, and $\rho_{\delta}$ is a positive regularization kernel, approximating the Delta-distribution as $\delta\to 0$. The formal limit of $E_{\text{obst},\delta}$ as $\delta\to 0$ takes the value $1$ when the constraint (3) is satisfied, and the value infinity when it is violated.

If only the three previous energy contributions are considered, we would expect a relaxation to an equilibrium. Movement requires an active component, coming from actin polymerization and depolymerization in the cortex. We assume that cell polarization manifests itself by local imbalances of this process producing a net increase of actin close to the cell front and a decrease close to the rear of the cell. For simplicity we make the equilibrium assumption that the total amount of actin in the cortex does not change.

We introduce the arclength $l$, which is given by

[TABLE]

This relation between the arc length $l$ and the Lagrangian variable $s$ can be inverted in terms of the actin density $\rho(l,t)$ per arc length:

[TABLE]

implying $|\partial_{s}X|=\rho^{-1}$. We denote by $f(l,t)$ the rate of actin increase ($f>0$) or decrease ($f<0$) per time unit and per arc length, so that $\rho$ satisfies

[TABLE]

The above mentioned equilibrium assumption translates to

[TABLE]

We then obtain the material derivative for functions of $s$:

[TABLE]

which has to be understood relative to the arc length along $\Gamma$, measured from the point $X(0,t)$. In particular, the velocity of the cortex relative to the laboratory coordinates is given by $DX/Dt$.

Whereas $f$ is the actin growth density w.r.t. arclength, we need a description in terms of the variable $s$.

[TABLE]

This roughly gives growth and decay rates per actin filament. This leads to

[TABLE]

Friction between the cell surface and the surrounding fluid, which is assumed nonmoving, is modelled by the $L^{2}$-gradient flow for the total energy with the contributions (1), (2), and (4), where the friction force is given in terms of the material derivative:

[TABLE]

Note that the friction coefficient has been eliminated by an appropriate choice of the time scale. This is however not yet the final form of the model. It would predict movement of freely floating cells in the absence of any confinement, as can easily be seen by integrating the equation in the absence of the last term with respect to $s$:

[TABLE]

Since the right hand side will in general be different from zero, the center of mass will move. An explanation is the force, used to move depolymerized G-actin from regions with $g<0$ to regions with $g>0$. We do not describe the corresponding mechanism further, which could include diffusion or transport by molecular motors, for example. By the action-reaction principle it seems reasonable to introduce a compensating counterforce with density $\mathcal{F}_{\mathrm{comp}}(s,t)$, acting on the cortex:

[TABLE]

It can be chosen arbitrarily, except that the total force is fixed:

[TABLE]

We do not make a choice for $\mathcal{F}_{\mathrm{comp}}$ at this point, the precise choice made for our numerical experiments will be discussed in Section 5. One can note that $\mathcal{F}_{\mathrm{comp}}$ is roughly directed from the front towards to back of the cell, since $g$ is negative (respectively positive) where the depolymerization (respectively polymerization).

4 Existence results

4.1 Formulation of the main results

We shall prove a global existence result for the initial value problem for (10) and a convergence result for the limit $\delta\to 0$ of hard channel walls, i.e. strict enforcement of the constraint (3).

Since the terms describing the cortical flow are chosen to have a vanishing integral with respect to $s$, we introduce $S_{\mathrm{flow}}$, which is chosen such that

[TABLE]

The initial value problem can then be written in the form

[TABLE]

For the function spaces, we use the abbreviations $\mathcal{L}:=L^{2}(\mathbb{T}^{1})^{2}$ and $\mathcal{H}:=H^{1}(\mathbb{T}^{1})^{2}$ with the norms $\|\cdot\|_{2}$ and $\|\cdot\|_{1,2}$, respectively. The torus $\mathbb{T}^{1}$ is represented by the $s$-interval $(0,1)$. For the problem parameters we shall use the following assumptions:

1. (A1)

$p<2\pi$, 2. (A2)

$\|S_{\mathrm{flow}}[X]\|_{2}\leq c_{\mathrm{flow}}\|X\|_{2}$,  $\|\partial_{s}(S_{\mathrm{flow}}[X])\|_{2}\leq c_{\mathrm{flow}}\|X\|_{1,2}$. 3. (A3)

$W_{\delta}$ is given by (4), where the domain $\Omega_{c}\subset\mathbb{R}^{2}$ has a smooth boundary, and $x\cdot\nabla W_{\delta}(x)\geq 0$ for $x\in\mathbb{R}^{2}$ and for $\delta$ small enough. 4. (A4)

$X^{0}\in\mathcal{H}$ and $\{X^{0}(s):\,s\in(0,1)\}\subset\Omega_{c}$.

Assumption (A1) can be motivated by looking at the simplified problem without cortical flow and without obstacles, i.e. $\nabla W_{\delta}=S_{\mathrm{flow}}=0$. In this case the expected circular equilibrium state $X_{equ}$ exists only under Assumption (A1) and is given by

[TABLE]

Assumption (A2) is somewhat restrictive compared to the models discussed in Section 3. It can be satisfied for bounded actin growth rate $g$. The inequality in Assumption (A3) means roughly that the allowed domain $\Omega_{c}$ is star shaped (and that its center has been taken as the origin). It is a technical assumption which in general precludes the type of channels used experimentally. Finally, Assumption (A4) on the initial cortex shape implies that the elastic and pressure energies are finite and that the cell lies in the admissible region. It might be noted that we do not assume that $X_{0}$ parametrizes a simple curve. We do not refer to this property since our results do not guarantee that it is preserved globally in time.

Theorem 1**.**

*(Global existence for the penalized problem)

Let the Assumptions (A1)–(A4) hold and let $\delta>0$ be small enough. Then there exists a solution $X_{\delta}$ of (13), (13), such that*

[TABLE]

uniformly with respect to $\delta\to 0$.

Remark 1**.**

By the Morrey inequality proved in the appendix, $X_{\delta}$ is Hölder continuous with exponent $\frac{1}{2}$ in terms of $s$ and $\frac{1}{4}$ in terms of $t$, again uniformly in $\delta$.

The proof will be carried out in the following three sections. It relies on methods for gradient flows, although they cannot be applied in a straightforward way. The terms on the right hand side of (13) are the $L^{2}$-gradients of the energy functionals $E_{\mathrm{el}}$, $E_{\text{p}}$, and $E_{\text{obst},\delta}$ (see Section 3), the first of which is convex. The second and third are treated as continuous perturbations. The cortical flow term on the left hand side is nonvariational.

Our approach is based on a semi-implicit time discretization, where the cortical flow term is evaluated at the old time step. This allows to solve the discrete problem by energy minimization. A priori estimates for the discrete solution allow to pass to the continuous limit.

These estimates are also uniform in the penalization parameter $\delta$ for the potential, so we can also carry out the high penalization limit:

Theorem 2**.**

*(Limit of hard channel walls)

With the assumptions of Theorem 1, the family $\{X_{\delta}:\,\delta>0\}$ of solutions of (13), (13) contains a sequence, converging (as $\delta\to 0$) uniformly on bounded time intervals to*

[TABLE]

which satisfies $X\in\overline{\Omega_{c}}$ on $\mathbb{R}_{+}\times\mathbb{T}^{1}$, (13), and

[TABLE]

for all $(s,t)$ such that $X(s,t)\in\Omega_{c}$.

The properties of $X$ stated in the theorem are not a complete formulation of the obstacle problem. Information on the behavior at the edges of contact regions is missing. Since the equation is degenerately parabolic, this is not so obvious. Under the additional assumption of convexity of the permissible set $\Omega_{c}$, one could expect that $X$ solves the variational inequality

[TABLE]

for all $Y\in\mathcal{H}$ such that $Y(s)\in\overline{\Omega_{c}}$, $s\in\mathbb{T}^{1}$.

For the numerical experiments, the obstacle is not modelled by a potential. For any $X(s,t)\in\Gamma(t)$, the total resulting force is rather restricted to the tangent cone to $\Omega_{c}$ at $X(s,t)$. In other words, our numerical simulations will be based on the formulation

[TABLE]

where the potential $W_{\delta}$ does not appear. The $\mathbf{P}_{c}$ is the projection on the tangent cone:

[TABLE]

where $\tau_{c}$ and $n_{c}$ are a normalized tangent vector and the normalized outward normal along $\partial\Omega_{c}$.

4.2 The energy functional and its properties

We introduce the total energy functional

[TABLE]

and prove its coercivity:

Lemma 1** (Coercivity).**

For $p<2\pi$ and for all $X\in\mathcal{L}$

[TABLE]

Proof.

W.l.o.g. we assume $X\in\mathcal{H}$ and denote by $\Omega$ the enclosed domain and by $L$ the length of its boundary. Then the isoperimetric and the Cauchy-Schwarz inequalities imply

[TABLE]

For the elastic energy again the Cauchy-Schwarz inequality is used:

[TABLE]

which completes the proof, since $E_{\text{obst},\delta}\geq 0$. ∎

Remark 2**.**

This result shows that the definition of $E_{\delta}$ has been appropriate, since finiteness of $E_{\delta}(X)$ implies $X\in\mathcal{H}$.

Lemma 2**.**

With respect to the weak topology in $\mathcal{H}$, the functional $E_{\mathrm{el}}$ is convex and lower semicontinuous, and the functionals $E_{\mathrm{p}}$ and $E_{\text{obst},\delta}$ are continuous.

Proof.

The integrand of $E_{\mathrm{el}}$ is positive, and a convex function of $\partial_{s}X$. Lower semicontinuity follows from Giaquinta[6], Theorem 2.5. The other properties are straightforward. ∎

4.3 Time discretization

Choose $\tau>0$ and $X^{n-1}\in\mathcal{H}$. Then, by the results of the preceding section, the functional

[TABLE]

is weakly lower semicontinuous on $\mathcal{H}$. It is also bounded from below since $E_{\delta}$ is, and since

[TABLE]

Furthermore, sublevel sets are bounded in $\mathcal{H}$. This is sufficient for the minimum of $\Phi(\tau,X^{n-1};\cdot)$ to be assumed in $\mathcal{H}$, and we choose

[TABLE]

We define $X_{\tau}\in C^{0,1}_{\mathrm{loc}}(\mathbb{R}_{+},\mathcal{H})$ as the piecewise linear interpolation of the $X^{n}$, $n\geq 0$. More precisely, we have

[TABLE]

where $k_{t}\in\mathbb{N}$ is such that $k_{t}\tau\leq t<(k_{t}+1)\tau$.

Lemma 3** (Uniform estimates for the interpolation).**

Let Assumptions (A1)–(A4) hold. Then

[TABLE]

uniformly with respect to small enough $\tau$ and $\delta$.

Proof.

Denoting the duality bracket with $\langle\cdot,\cdot\rangle$, the variation of $\Phi$ in the direction $Y\in\mathcal{H}$ gives

[TABLE]

with the formal gradient of the energy determined from

[TABLE]

With $Y=X^{n}$ and with Assumptions (A1) and (A3) we get, similarly to the proof of Lemma 1,

[TABLE]

and we observe

[TABLE]

Finally, we use

[TABLE]

with $0<\gamma<1-p/(2\pi)$. This implies the existence of positive constants $A_{1},A_{2},A_{3}$ (independent from $n$, $\tau$ and $\delta$) such that

[TABLE]

A discrete Gronwall estimate now gives

[TABLE]

and, as a consequence,

[TABLE]

From the last two bounds we get for any $t_{1},t_{2}>0$

[TABLE]

and

[TABLE]

where $k_{2}\tau\leq t_{2}<(k_{2}+1)\tau$. In other words, we have shown $X_{\tau}\in L_{\mathrm{loc}}^{\infty}(\mathbb{R}_{+},\mathcal{L})\cap L_{\mathrm{loc}}^{2}(\mathbb{R}_{+},\mathcal{H})$ uniformly in $\tau$ and $\delta$.

From the minimization we get

[TABLE]

With the velocity $\partial_{t}X^{n}_{\tau}=\frac{X^{n}-X^{n-1}}{\tau}$ of the linear interpolant this reads

[TABLE]

With the Young inequality and Assumption (A2) we obtain

[TABLE]

Using (20), summation over $n$ and an application of Lemma 1 completes the proof.

∎

4.4 The continuous and high penalization limits

Considering the time discrete solution $(X^{n})_{n\geq 0}$, we recall the definition (18) of $X_{\tau}$ and define two additional piecewise constant continuous-time approximations:

[TABLE]

Then (19) implies

[TABLE]

where $\partial E_{\mathrm{el}}(X)$ is the subdifferential of the elastic energy, given by the distributional derivative with respect to $s$ of $-(|\partial_{s}X|-1)_{+}\frac{\partial_{s}X}{|\partial_{s}X|}$. The equality only holds in the dual space of $\mathcal{H}$ a priori, but thanks to Lemma 3, it also holds in $\mathcal{L}$. It is our goal to pass to the limit $\tau\to 0$ in this equation. By the results of the previous section, there exists a sequence $\tau_{k}\to 0$ such that

[TABLE]

We also have that all the terms in (21) except $\partial E_{\mathrm{el}}(X_{\tau}^{\mathrm{new}})$ are bounded in $L_{\mathrm{loc}}^{2}(\mathbb{R}_{+},\mathcal{L})$ uniformly in $\tau$ and $\delta$. In all these terms we can pass to the limit in the sense of distributions, but also weakly in $L_{\mathrm{loc}}^{2}(\mathbb{R}_{+},\mathcal{L})$.

This also implies weak convergence of $\partial E_{\mathrm{el}}(X_{\tau_{k}}^{\mathrm{new}})$ to some $\eta\in L_{\mathrm{loc}}^{2}(\mathbb{R}_{+},\mathcal{L})$, and we obtain

[TABLE]

It remains to identify $\eta$.

Since the subdifferential is a maximal monotone operator, we shall apply the Minty trick [8, 11, Theorem 2.2]. Testing (21) against $X_{\tau}^{\mathrm{new}}$, its strong convergence immediately implies

[TABLE]

which is sufficient for $\eta\in\partial E_{\mathrm{el}}(X)$ and, thus,

[TABLE]

to be understood as the weak derivative of an $L^{2}$-function, since $X\in L_{\mathrm{loc}}^{\infty}(\mathbb{R}_{+},\mathcal{H})$. This completes the proof of Theorem 1.

For the proof of Theorem 2, we note that by the results of Theorem 1 and following Remark 1 the uniform convergence of a subsequence to $X$ is immediate. The uniform boundedness of $E_{\text{obst},\delta}$ implies $X\in\overline{\Omega_{c}}$ by its continuity and by the uniform convergence, since $W_{\delta}(x)\to\infty$ for $x\notin\overline{\Omega_{c}}$.

Similarly, since $\nabla W_{\delta}(x)\to 0$ for $x\in\Omega_{c}$, the weak formulation of (14) can be derived by using test functions vanishing away from $\{(s,t):\,X(s,t)\in\Omega_{c}\}$. The convergence of the various terms is handled exactly as for the continuous limit.

5 Numerical results

5.1 Discretization

We start with the situation without obstacles and introduce $N\in\mathbb{N}$ grid-points for the discretization of $s$ such that $s_{i}=i\Delta s$, $\Delta s=\frac{1}{N}$ for $i$ considered on a discrete torus, meaning that $i$ is identified with $i+N$. The time step is $\Delta t>0$ and $t^{n}=n\Delta t$, $n\geq 0$. The numerical approximation for $X(s_{i},t^{n})$ is denoted by $X_{i}^{n}$.

We assume a cell which is polarized in the fixed direction $\omega\in\mathbb{R}^{2}$, $|\omega|=1$ and define $i_{0}^{n}$ and $i_{1}^{n}$ such that

[TABLE]

Actin is added to the cortex at the leading end $X_{i_{0}^{n}}^{n}$ and removed at the trailing end $X_{i_{1}^{n}}^{n}$, corresponding to $s_{0}(t^{n})=i_{0}^{n}\Delta s$, $s_{1}(t^{n})=i_{1}^{n}\Delta s$ in the notation of Section 3.

We use the explicit Euler scheme for the time discretization and symmetric finite differences for the discretization in the $s$-direction:

[TABLE]

where the elasticity forces are given by

[TABLE]

and the compensating force by

[TABLE]

Here $\tilde{\delta}_{a}(i)$ is a mollification of the Dirac Delta in $s_{i}$, with smoothing parameter $a$ such that $\sum_{i}\tilde{\delta}_{a}(i)\Delta s=1$. The approximation for the compensating force has been chosen such that the discrete center of mass $\sum_{i}X_{i}^{n}\Delta s$ is independent of time (in the absence of obstacles).

The restriction to the admissible domain $\Omega_{c}$ is enforced by an approximation of the formulation (15), (16):

[TABLE]

where

[TABLE]

with the tangent and normal vectors evaluated at the orthogonal projection of $X_{i}$ to $\partial\Omega_{c}$ (see Figure 4). The width $\epsilon$ of the tube, where the projection is applied, is kept sufficiently large with respect to the time step such that the constraint (3) is met.

5.2 Dimensionalisation of the model parameters

We denote by $\kappa_{S}$ the elasticity force and by $\mu$ the internal friction coefficient in front of the time derivative. Note that all these quantities have been set to 1 so far in the model, without loss of generality. We recover them here for the sake of computing their dimensions, Eq. (15) reads:

[TABLE]

We denote by $x_{0}$ and $t_{0}$ the space and time units of the model. The space unit is chosen to be half the cortex length of leucocytes $L=2x_{0}$, i.e $x_{0}\approx 47.6\mu m$. The polymerization speed in leukocytes is around $12\mu m.min^{-1}$ and the dimensionless polymerization speed writes $v=2\frac{x_{0}}{t_{0}}$ therefore 1 time unit of the model corresponds to $t_{0}\approx 8min$. The elastic properties of the human red blood cell have been studied via micropipette aspiration experiments [21] where the authors show that the stiffness constant of leukocytes membrane is of order $k_{S}=7pN.\mu m^{-1}$. In our model $\frac{k_{S}}{\mu}=\frac{1}{t_{0}}=0.1265min^{-1}$ therefore we can deduce that the internal friction coefficient $\mu\approx 55pN.\mu m^{-1}.min$. Finally, we choose the pressure constant $p=3.2$. We deduce $\frac{p}{\mu}=3.2\frac{1}{t_{0}}=0.4min^{-1}$ and therefore $p\approx 22pN.\mu m^{-1}$ which is in biological range. All the model parameters are summarized in Table 1.

5.3 Simulations without Obstacle

In Figure 5 we show computed equilibrium shapes of the cortex in the absence of obstacles. The polarization of the cell is always to the left (i.e. $\omega=(-1,0)$ in (22)).

Figure 5 (A) shows a simulation without compensating force. As expected, the cell migrates in the direction of the leading end with a speed slightly less than the polymerization/depolymerization speed. The equilibrium shape is a circle.

Figures 5 (B,C,D) show the equilibrium shapes obtained for different values of $h=\frac{a}{L}$, where we recall that $a$ measures the spread of the compensating force and $L$ is the equilibrium circumference. The most important observation is that migration has been turned off successfully by including the compensating force. The deviation of the equilibrium shape from a circle is stronger for more concentrated compensating forces (smaller values of $a$ therefore $h$, compare Figures (D) to (B)).

5.4 Simulations of Migration in Channels

In order to reproduce the biological experiments, our protocol for simulations with channels is started by putting the cell with an initially circular shape close to the opening of a channel. Then we let it evolve to an equilibrium shape as in the preceding section, after which we turn on a force pushing it into the channel. This is achieved by turning off the pressure force in a region around the trailing end. The pushing is removed as soon as the cell is entirely inside the channel.

In channels with flat walls with varying widths and spreads of the compensating force (parameter $h$), no migration is observed in the simulations in spite of polarization and corresponding cortex flow (Figure 6). These results are in agreement with the experiments described in Section 2.

Ratchet channels are described by four parameters: (i) a wave length $L_{0}$ of the width variations, (ii) a minimal width $2w_{0}$, (iii) an amplitude of the width variations $d_{0}$, and (iv) an asymmetry parameter $\alpha$. For the length $x$ along the channel, the walls of the channel are given by $\pm f(x)$, with the function

[TABLE]

where $g(x)$ solves the fixed point equation:

[TABLE]

In all our numerical simulations we choose $\alpha=0.4$. Note that for $\alpha=0$, the walls would be sinusoidal functions of $x$, i.e. symmetric with respect to the $y$ axis. As in the experiments, in each simulated channel we combine three wavelengths, $L_{0}=3.9,7.6,11.7\mu m$, increasing in the direction of migration. Figures 7 (A,B) show two examples with $2w_{0}=1.36\mu m$, $d_{0}=0.76\mu m$ (A) and $2w_{0}=3.7\mu m$, $d_{0}=3\mu m$ (B). As further illustration, Figure 7 (C) shows a cell at equilibrium before being pushed into a ratchet with $L_{0}=3.3\mu m$ and $7.6\mu m$, $2w_{0}=3.7\mu m$, and $d_{0}=1.54\mu m$.

Results of a typical simulation are shown in Figure 8. In this situation the cell is able to migrate in all three different wavelengths. It seems that the average speed is highest in the part with the intermediate wavelength. Within each wavelength region, the behavior seems to be periodic related to the periodicity of the channel walls. The bigger the wavelength, the more pronounced are the peaks in the cell speed.

A study of the dependence of the mean cell speed on the various parameters has been carried out (Figure 9). Not surprisingly, increasing the channel width decreases the cell speed and can prevent cell migration. For ratchets of the smallest wavelength 3.9 $\mu$ m, the cell is unable to migrate for large amplitudes. In this situation the cell is unable to enter the ratchets, thus reducing its contact to the wall. The channel acts in these cases as if it had flat walls. For period $7.6\mu m$ ratchets, the optimal speed is obtained for large enough amplitude. This is expected, since for small amplitude the ratchet again acts as a flat walled channel. This fact is reinforced in channels of period 11.6 $\mu m$, where only ratchets with large amplitude can induce cell migration.

To sum up, cell speed is directly linked to the cell compression induced by the channel geometry and to the possibility for the cell to push against parts of the channel wall facing towards the migration direction.

6 Conclusions

This study has been motivated by the experimental results described in Section 2, where adhesion-free migration of leukocytes in artificial, structured micro-channels has been observed.

A mathematical model for this process has been formulated. Migration is assumed to be due to cortex flow, driven by a local imbalance of polymerization and depolymerization in a polarized cell. Cell shape is stabilized by cytoplasmic pressure and elastic behavior of the cortex. An undesired tangential friction, caused by the cortical flow, is balanced by a compensating force, which can be attributed to the internal transport of depolymerized actin.

The model has the form of an obstacle problem for a strongly nonlinear degenerate parabolic system. Global existence of solutions has been proven under natural assumptions on the data. The analysis relies on ideas from the theory of gradient flows, employing the structure of the dominating elastic and pressure terms. The results are complete for an approximate system with “softened” obstacles. The limit for hard obstacles can be carried out, however with an incomplete characterization of the limiting problem.

For numerical simulations, a conservative, explicit-in-time discretization has been introduced. Under appropriate time step restrictions, simulations are stable and (at least qualitatively) reproduce the behavior observed in the experiments. In particular, migration needs both confinement and sufficiently structured channel walls. A parametric study shows the expected dependencies on geometric properties of the channel.

From a modelling point of view, this study has to be seen as a first step. Reliable experimental information on cell cortex structure and dynamics is still scarce. The fact that in our model migration strongly depends on the force compensating excess polymerization and depolymerization, is rather questionable. In ongoing work, the model is extended by a viscous resistance against cortex bending. This effect seems to be a reasonable alternative providing the necessary pushing force against the channel walls. However, its inclusion poses new challenges both from an analytic and from a numerical point of view.

Appendix A Appendix: A modified Morrey inequality

Lemma 4**.**

Let $u\in L^{\infty}((0,T);H^{1}(0,1))\cap H^{1}((0,T);L^{2}(0,1))=:L_{t}^{\infty}H_{s}^{1}\cap H_{t}^{1}L_{s}^{2}$. Then for any $(s_{0},t_{0}),(s_{1},t_{1})\in(0,1)\times(0,T)$ we have:

[TABLE]

Proof.

We introduce $\Delta s:=|s_{1}-s_{0}|$ and $\Delta t:=|t_{1}-t_{0}|$ and consider a rectangle $W\subset(0,1)\times(0,T)$, containing the points $(s_{0},t_{0}),(s_{1},t_{1})$ with sides parallel to the $t$-axis of lengths $\Delta t$ and parallel to the $s$-axis of lengths $\Delta s+\sqrt{\Delta t}$, such that $|W|=\Delta t(\Delta s+\sqrt{\Delta t})$. We have

[TABLE]

For estimating the second term, we introduce the curve $\{(s_{0}+(\sigma-s_{0})\sqrt{p},t_{0}+(\tau-t_{0})p):\,0\leq p\leq 1\}$, whence it can be estimated by

[TABLE]

where the derivatives of $u$ are evaluated along the curve. Employing the Cauchy-Schwarz inequality, this can be estimated further by

[TABLE]

In the second integral over $W$ we introduce the new coordinates $(s,t)=(s_{0}+(\sigma-s_{0})\sqrt{p},t_{0}+(\tau-t_{0})p)$; in the first one we first estimate the integrand by its supremum with respect to $t$ and then make the coordinate transformation only in $s$:

[TABLE]

An analogous treatment of the first term on the right hand side of (25) completes the proof. ∎

Acknowledgments

This work has been supported by the Vienna Science and Technology Fund, Grant no. LS13-029. G.J. and C.S. also acknowledge support by the Austrian Science Fund, Grants no. W1245, F 65, and W1261, as well as by the Fondation Sciences Mathématiques de Paris, and by Paris-Sciences-et-Lettres.