A stochastic model for protrusion activity

In this work we approach cell migration under a large-scale assumption, so that the system reduces to a particle in motion. Unlike classical particle models, the cell displacement results from its internal activity: the cell velocity is a function of the (discrete) protrusive forces exerted by filopodia on the substrate. Cell polarisation ability is modeled in the feedback that the cell motion exerts on the protrusion rates: faster cells form preferentially protrusions in the direction of motion. By using the mathematical framework of structured population processes previously developed to study population dynamics [Fournier and M{\'e}l{\'e}ard, 2004], we introduce rigorously the mathematical model and we derive some of its fundamental properties. We perform numerical simulations on this model showing that different types of trajectories may be obtained: Brownian-like, persistent, or intermittent when the cell switches between both previous regimes. We find back the trajectories usually described in the literature for cell migration.


Introduction
Cell migration is a fundamental process involved in physiological and pathological phenomena such as the immune response, morphogenesis, but also the development of metastasis from a tumor [1,5]. To ensure these functions, cells have a highly complex out-of-equilibrium internal organization where multiscale reactions occur among polymers and molecules, leading to unpredictable macroscopic behaviours.
In the case of cell crawling, cells spread on an adhesive substrate, and form extensions also called protrusions. Then, molecular adhesion complexes grow and ensure a mechanical connection between protrusions and the substrate, by which forces are transmitted and lead to a displacement. Protrusions of a crawling cell can be divided in two types: lamellipodia are wide and flat and fluctuate continuously, while filopodia are long finger-like extensions able to grow further and probe the substrate.
It has been observed that cells protrusive activity fluctuates a lot, and that these fluctuations are responsible for the long-term characteristics of trajectories [2]. Cell trajectories can be very different even for a single cell type: some do not explore the environment, while others have a much more efficient displacement. It is of interest to try to capture this diversity in a mathematical model.
Existing stochastic models for cell trajectories are either Random Walks, Lévy flights or Active Brownian Particle models [9]. In these models, key-features of the motion are quantified, such as the mean persistence time, or the stationary distribution of the particle's velocity [9]. However, the dynamics is macroscopic and Figure 1. Scheme of a polarised crawling cell: protrusive structures at the front (lamellipodium, filopodia), and contractile fibers at the back. The asymmetric activity combines with an asymmetric repartition of molecular regulators organized in feedback loops. Source: [7] processes such as polarisation are taken into account by an arbitrary positive feedback at the scale of the trajectory. In this work, we present a 2D stochastic particle model for cell trajectories based on the filopodial activity, able to reproduce the diversity of trajectories observed.

Model construction
We choose space and time scales large enough so that the cell is described as an active particle (its center of mass). Cell shape and intracellular dynamics are therefore considered at the cell scale. We define now the velocity model that is based on force equilibrium.

Velocity model
At each time, the cell velocity writes V t , and its polar coordinates (v t , θ t ). The crawling of cells on an adhesive substrate occurs at very small scales. Indeed, cells sizes are of the order of 1 − 10 µm, while their speeds ranges at the scale of µm s −1 . Therefore, inertia is negligible for this system (see more precise justifications of the low Reynolds number setting in e.g [10]), and Newton's second law of motion reduces to instantaneous force equilibrium: at all time t ≥ 0, The cell being an active system, macroscopic forces that apply can be either passive or active. In this case appear → a passive force: the friction force exerted by the substrate on the cell due to motion, that writes f = −γ V t , with γ the global friction coefficient, → active forces related to the protrusion process. Indeed, a complex internal activity gives rise to forces in the body of the cell. Filopodial protrusions can be considered as good readouts [2]. As a consequence, in the following, only filopodial forces will be considered. Note that at this scale, the formation of filopodia is discontinuous in time.
Combining these information, we get: where N t is the number of filopodia adhering on the substrate at time t, and ( F i (t)) i the filopodial forces. The cell motion is then entirely described by the protrusions. Now, denoting θ i = arg( F i ), one can write Modelling cell motion then accounts to modelling the time evolution of the filopodial population in terms of individual orientation.

Protrusion model
Each filopodium is characterized by a quantitative parameter, its orientation θ ∈ [0, 2π). Therefore, we use a measure valued process for the stochastic evolution of the set of filopodia, as in ecological population models [4]. Let us denote M F (χ) the set of positive finite measures on χ = [0, 2π], equipped with the weak topology. Notice that as χ is compact, weak and vague topologies on M F (χ) coincide. Write M for the subset of M F (χ) composed of all finite point measures. Then, a filopodium of orientation θ is described by a Dirac measure δ θ on χ, and the whole population by For any measurable function f on χ and any µ ∈ M F (χ), we have < µ, f >= χ f (θ)µ(dθ) . In particular, , and the population size corresponds to N t =< ν t , 1 >. A simple way to express the velocity equation (1) together with hypothesis (2.1) is to write The cell motion is entirely described by a measure-valued markovian jump process (ν t ) t , as in adaptive stuctured population models. We describe now the different events arising.
• The basic dynamics arising is the isotropic appearance of filopodia. It is responsible for the spontaneous activity that is observed experimentally. We write c for the creation rate. • Each filopodium ends up disappearing: the disappearance or death rate is denoted by d constant.
• Polarisation is characterized by a morphological and functional asymmetry visible both on cell shape and at the microscopic scale [3,8]. Here, we use the mesoscopic scale of the model to account for polarisation by its feedback on the protrusive activity. Two phenomena have to be distinguished: → The formation of a protrusion is induced by several microscopic regulators and generates a local positive feedback on the protrusive machinery. Following that, we assume that each filopodium is able to reproduce. Denote r(θ i , ν t ) the individual reproduction rate of a filopodium of orientation θ i . → Polarisation is also reinforced by intracellular actin flows, (see [6]). In particular, faster actin flows favor the formation of protrusions in a single stable configuration. Denote u for the spaceaveraged actin flow velocity over the cell. As actin flows are inwardly directed, − u characterizes the reinforced direction for protrusions. Moreover, we know that V = − 1 α u, where 1 α depends on the cell type and the experimental setting. Therefore, we consider a positive coupling between the reproduction rate and − u = α V , imposing a global feedback.
• Reproduction of spatially localized filopodia questions the localization of the new protrusions. An individual can reproduce to form a filopodium with the same orientation, or it can have a slightly different location. This phenomenon accounts for the stochastic fluctuations arising in the cell signalling pathways involved in protrusions. In our model, we describe this using the notion of heredity and subsequent mutation event for the orientation of the "offspring". For simplicity, we assume that at each reproduction event, the mutation probability µ is constant. In the case of a mutant new protrusion, its orientation is determined following a probability distribution g(z; θ i ) assumed centered in the parent's orientation θ i , with a constant variance. The possible events are summed up in the following graph: Let us comment on the mathematical features of the model. In the case of no interaction between individuals, or global feedback, the process (N t ) t simply follows an immigration, birth and death dynamics, and mathematical information can be derived. In particular, the branching property still holds. This is no longer the case when adding interactions. For example, if the interaction relies on V t , then knowing only N t is not sufficient and one has to know about the structured quantities (N θ1 , N θ2 , ...) at all time.
A choice of reproduction rate and mutation law As protrusions are located on [0, 2π), it is natural to consider reproduction rates as circular functions. We chose a reproduction function that is positively correlated to α V t to account for polarisation. The idea is to have a function centered in θ t the direction of motion, that gets sharper with increasing v t .
We choose a reproduction rate as a multiple of a circular normal distribution density written with κ(|| V t ||) ≥ 0 a non decreasing shape parameter, and I 0 the 0-order modified Bessel function of the first kind. Therefore, we denote r(θ, ν s ) = r * f (θ; θ t , κ).
The mutation law is a circular normal distribution on [0, 2π) of density g(z; θ i ), centered in θ i and with a constant shape parameter (resp. variance) κ (resp. σ 2 ).

Mathematical properties
From now on, we will use the notation C for any constant, that will change from line to line. We introduce here a stochastic differential equation for (ν t ) t driven by Point Poisson Measures. We will show existence and uniqueness of a solution, and prove that it follows the dynamics previously described.
In order to pick a specific individual in the population, we have to be able to order them or their trait. Indeed, from the n-uplet (θ 1 , ..., θ N ), one can recover ν = N i=1 δ θi , but from ν it is only possible to know {θ 1 , ..., θ N }. Definition 3.1. Let us define the function with θ σ(1) θ σ (2) ... θ σ(n) , for an arbitrary order . Now, an individual can be picked by its label i, and the corresponding trait writes H i (ν) = θ σ(i) .
Let (Ω, F, P) be a probability space, and n(di) the counting measure on N * . We introduce the following objects: • ν 0 ∈ M the finite point measure describing the initial population, eventually equal to the null measure.
In this equation, each term describes a different event. The Poisson Point Measures generate atoms homogeneously in time. However, the dynamics we want to describe follows state-dependent rates. Hence, we use indicator functions to keep only some of the events in order to get the wanted rates. Then, the Dirac measures correspond to the individuals added to or removed from the population.

Existence and uniqueness
In this part let us prove existence and uniqueness of a solution for equation (2). Recall that N t =< ν t , 1 >. Proposition 3.3. Assume the boundedness of the reproduction rate (hypothesis 3.2), and that E[N 0 ] < +∞. Then, the two following properties hold.
Proof of 3.3. The proof is similar to prop. 2.2.5 and 2.2.6 in [4].
(1) Let T 0 = 0, and t ∈ R + . Then, the global jump rate of ν t is smaller than c + (r + d)N t . Hence one can P − a.s define the sequence (T k ) k∈N * of jumping times, as well as T ∞ := lim k→+∞ T k . Now, by construction, it is P − a.s possible to build "step-by-step" a solution of equation (2) on [0, T ∞ [. Showing existence of a solution (ν t ) t∈R+ ∈ D(R + , M(χ)) amounts to showing that P − a.s, T ∞ = +∞. That is equivalent to saying that there cannot be an infinite number of jumps in a finite time interval.
• First, we show the control property (3). For n > 0 define the sequence of stopping times (τ n ) n by is indeed a sequence of stopping times.
→ Now, we prove that for all T < +∞, the quantity E sup t∈[0,T ∧τn] N t is bounded ∀n ≥ 0.
For t ∈ R + , using equation (2) and dropping the non-positive term, one has N t∧τn = < ν t∧τn , 1 >≤ N 0 + As each integrand is positive, bounded, and integrable with respect to the intensity measure, taking the expectation and using the Fubini theorem, we can write , which is in contradiction with equation (4). → Property (3) is proved by the Fatou lemma: • Now, let us show that P − a.s, T ∞ = +∞. If this is not the case, then there exists M < +∞ and a set A M ⊂ Ω such that P(A M ) > 0 and ∀w ∈ A M , T ∞ (ω) < M . Moreover, if the assertion is true, then we would have which contradicts lim n→+∞ τ n = +∞. As a consequence, if we prove (4), the proposition is proved. If (4) is not true, there would exist N > 0 and a set B ⊂ A M such that P(B) > 0 and Then, ∀ω ∈ B, (T k (ω)) k can be seen as the subsequence of a sequence of jumping times (T 1 k (ω)) k of a Point Poisson Process of intensity c + (r + d)N . The only accumulation point of (T 1 k (ω)) k being P − a.s + ∞, it contradicts the definition of B, and proves (4).
(2) The sequence of jumping times (T k ) k∈N being already defined, we only have to show that (T k , ν T k ) k∈N are uniquely determined by D = (ν 0 , M 0 , M 1 , M 2 , M 3 ) defined above. But this is clear by construction of the process.

Markov property
Now, we can show that the solution (ν t ) t of equation (2) is a Markov process in the Skorohod space D(R + , M F (χ)) of càdlàg finite measure-valued processes on χ. For that purpose, we introduce ∀ν ∈ M, Φ : M → R measurable and bounded, the operator L defined by Proposition 3.4. Take (ν t ) t≥0 the solution of equation (2) with E[< ν 0 , 1 >] < +∞. Then, (ν t ) t≥0 is a Markovian process of infinitesimal generator L.
In particular, this proposition ensures that the law of (ν t ) t≥0 is independent of the order involved in (3.1).
Proof of proposition 3.4. The process (ν t ) t≥0 ∈ D(R + , M(χ)) is markovian by construction. Now, let N 0 < N < +∞, and consider again the stopping time τ N . Let Φ : M → R be measurable and bounded. As P − a.s we can write we have Again, as all integrands are bounded, we can take expectations to get On the one hand, ∀t ∈ [0, T ], On the other hand, t → ψ(t ∧ τ N , ν) is derivable in t = 0 P − a.s (as ν ∈ D(R + , M(χ))), and for a given ν 0 , we have

or equivalently
Lφ

Numerical simulations
The construction of the process (ν t ) t furnishes directly an algorithm for simulations. We proceed as follows: start with the population measure ν k at time t k , for a particle located at X k .
Time of next event: let τ = c+ < ν k , r + d > denote the global jump rate of the process. Then, the time of the next event writes t k+1 := t k + ∆t, where ∆t ∼ Exp(τ ) .
Nature of the event: what happens at time t k+1 is determined as follows: • creation of a protrusion occurs with probability c τ . Its orientation is chosen uniformly on [0, 2π).
• reproduction of the protrusion number i occurs with probability r(H i (ν k ),ν k ) τ . Then, → with probability (1 − µ), the new protrusion has orientation H i (ν k ), → with probability µ, its orientation is chosen with the realization of a random variable having a probability density g(·; H i (ν k ), ν k ). • protrusion number i disappears with probability d τ . The measure ν k+1 is then obtained from ν k and the information of the event occuring at time t k+1 . Updates: the particle's new position is One only has to start again to get a trajectory over time.

Results
Let us now present the numerical trajectories we obtained, that are displayed in figures 3 and 4. Recall that polarisation is quantified by −αv, so that the larger α is, the more concentrated in the direction of motion the protrusions are formed. We observe indeed different types of trajectories for varying α, from Brownian-like to persistent. In figure 3, the mutation probability is µ = 0.2, whereas it is µ = 0.8 in figure 4. We observe that the territory exploration is significantly lower for a higher mutation probability. This shows that the mutation events can not be neglected.