A stochastic model for cell adhesion to the vascular wall

TL;DR

This paper models the stochastic process of cell adhesion to blood vessel walls, analyzing how blood flow velocity influences cell dynamics and adhesion probability through nonlinear stochastic models and differential equations.

Contribution

It introduces a novel stochastic model incorporating blood flow effects on cell adhesion, with analytical results on velocity thresholds and adhesion times.

Findings

Identification of shear-velocity threshold for adhesion transition

Derivation of cell mean stopping time as a function of flow and adhesion dynamics

Analysis of cell velocity dynamics under different scaling regimes

Abstract

This paper deals with the adhesive interaction arising between a cell circulating in the blood flow and the vascular wall. The purpose of this work is to investigate the effect of the blood flow velocity on the cell dynamics, and in particular on its possible adhesion to the vascular wall. We formulate a model that takes into account the stochastic variability of the formation of bonds, and the influence of the cell velocity on the binding dynamics: the faster the cell goes, the more likely existing bonds are to disassemble. The model is based on a nonlinear birth-and-death-like dynamics, in the spirit of Joffe and Metivier (1986); Ethier and Kurtz (2009). We prove that, under different scaling regimes, the cell velocity follows either an ordinary differential equation or a stochastic differential equation, that we both analyse. We obtain both the identification of a shear-velocity threshold associated with the transition from cell sliding and its firm adhesion, and the expression of the cell mean stopping time as a function of its adhesive dynamics.

Equations

$$V_t=u-\gamma F_t,$$

$$F_t=f{N}_t,$$

$$V_t=u-N_t.$$

$${\tau }_u:=\underset{t\ge 0}{\inf }\{N_t\ge u\}.$$

$$n\mapsto\left\{\begin{array}[]{ll}n+1&\text{ at rate }\lambda(n)=c(u)+rn\,,\\ n-1&\text{ at rate }\mu(n)=d(n)n\,,\\ \end{array}\right.$$

$$\mathbb{E}\left[\underset{t\in [0,T\wedge {\tau }_u]}{\sup }{N}_t^p\right]<+\infty \forall T>0.$$

$$\mathbb{E}\left[N_{t\wedge {\tau }_u}\right]=\mathbb{E}\left[N_0\right]+c(u)\mathbb{E}[t\wedge {\tau }_u]+\mathbb{E}\left[{\int }_0^{t\wedge {\tau }_u}\left(r-d{e}^{\alpha (u-N_s)}\right)N_s\mathrm{d}s\right]$$

$$\mathbb{E}[N_{t}]=\left\{\begin{array}[]{ll}\mathbb{E}[N_{0}]+c(u)t&\text{ if }r=d\,,\\ \mathbb{E}[N_{0}]e^{(r-d)t}+\frac{c(u)}{r-d}\left(e^{(r-d)t}-1\right)&\text{ otherwise.}\end{array}\right.$$

$$\mathbb{E}[N]_{\infty}:=\left\{\begin{array}[]{ll}\frac{c(u)}{d-r}&\text{ for }r<d\,,\\ +\infty&\text{ otherwise,}\end{array}\right.\Leftrightarrow\mathbb{E}[V]_{\infty}:=\left\{\begin{array}[]{ll}u-\frac{c(u)}{d-r}&\text{ for }r<d\,,\\ -\infty&\text{ otherwise,}\end{array}\right.$$

$$\text{Var}(N_{\infty })=\frac{c(u)}{d{\left(1-\frac{r}{d}\right)}^2}.$$

$$Z_t^K=\frac{1}{K}{N}_t^K\in \frac{1}{K}\mathbb{N}.$$

$$c^K(u)=Kc(u),r^K=r,\text{ and }{d}^K(K{Z}_t^K)=d(Z_t^K).$$

$$\underset{K>0}{\sup }\mathbb{E}\left[(Z_0^K{)}^2\right]<+\infty ,$$

$$n(t)=n_0+{\int }_0^{t}c(u)+(r-d(n(s)))n(s)\mathrm{d}s.$$

$$\underset{t\in [0,T]}{\sup }n(t)\le (n_0+cT)e^{rT}<+\infty ,$$

$$n^{\prime }(t)=c+\left(r-d{e}^{\alpha (u-n(t))}\right)n(t)=F(n(t)).$$

$$F^{\prime }(n)$$

$$F^{\prime \prime }(n)$$

$$U_{\alpha }:=\frac{1}{\alpha }\ln (r/d),{\overline{U}}_{\alpha }:=\frac{1}{\alpha }(1-r/d),U_c:=\frac{c}{d-r},\text{ and }\overline{C}:=\frac{1}{\alpha d}(r-d{)}^2.$$

$$c^K(u)=Kc(u),r^K=r+Ka,\text{ and }{d}^K(K{Z}_t^K)=d(Z_t^K)+Ka,$$

$$\underset{K>0}{\sup }\mathbb{E}\left[(Z_0^K{)}^2\right]<+\infty ,$$

$$\mathrm{d}{N}_t=b(N_t)\mathrm{d}t+\sigma (N_t)\mathrm{d}{B}_t,$$

$$N_{k+1}=\mid N_k+b(N_k)\Delta t+\sqrt{2a\Delta t{N}_k}W\mid ,$$

$$\frac{{\sigma }^2}{8}{\left(\frac{2b(0)}{{\sigma }^2}-1\right)}^2>3P\vee 4{\sigma }^2,$$

$$\frac{a}{4}{\left(\frac{c}{a}-1\right)}^2>(3P\vee 8a).$$

$$\mathrm{d}{N}_t=(c+(r-d)N_t)\mathrm{d}t+\sqrt{2a{N}_t}\mathrm{d}{B}_t,$$

$$f_{T_{x\to y}}(t)=\sum _{n=1}^{+\infty }{o}_n{\lambda }_n{e}^{-{\lambda }_{n}t},$$

$${\lambda }_n=(r-d)s_n,$$

$$o_n=-\frac{\Phi (s_n;c/a;\overline{x})}{s_n{\partial }_s(\Phi (s_n;c/a;\overline{y}))},$$

$${\lambda }_n\underset{n\to +\infty }{\sim }\frac{(d-r){\pi }^2}{4\overline{y}}{\left(n+\frac{c}{2a}-\frac{3}{4}\right)}^2-\frac{(r-d)c}{2a},$$

$$\begin{array}{r}{o}_n\underset{n\to +\infty }{\sim }\frac{(-1{)}^{n+1}2\pi (n+c/(2a)-3/4)}{{\pi }^2(n+c/(2a)-3/4{)}^2-\frac{2c}{a}\overline{y}}\times e^{\frac{1}{2}(\overline{x}-\overline{y})}{\left(\frac{\overline{x}}{\overline{y}}\right)}^{\frac{1}{4}-\frac{c}{2a}}\\ \cos \left(\pi \left(n+\frac{c}{2a}-\frac{3}{4}\right)\sqrt{\frac{\overline{x}}{\overline{y}}}-\frac{\pi c}{2a}+\frac{\pi }{4}\right).\end{array}$$

Full text

[TABLE]

1 Introduction

1.1 Biological context

Cell adhesion to the vascular wall is a major process involved e.g in inflammation and metastasis invasion (Granger and Senchenkova, 2010). The adhesive interaction between a cell circulating in the blood flow and endothelial cells forming the wall occurs in the presence of hemodynamic forces. Adhesion bonds can form between cell transmembrane proteins, called ligands, and adhesion receptors at the vascular wall surface. The first step of interaction happens when some bonds form and are stabilized, therefore slowing the cell down: this is the so called capture phase. Then, the cell rolls or slides along the endothelium, as new bonds form in the direction of motion and bonds at the back disassemble. This step is mediated by the selectin receptor molecules. During rolling, endothelial cells may also be stimulated, leading to the recruitment of integrins, another family of cell-receptor molecules. They mediate the cell firm adhesion which allows it to cross the vascular. The cell extravasation out of the vessel leads to further development of the process involved (immune response, invasion of tissues by metastatic cells e.g, see for example Ley et al (2007)). This process is depicted in Figure 1 for leukocytes.

The fate of a rolling cell is not clear, since it can either adhere to the wall or be released in the blood flow. Many experiments have been dedicated to the study, in vivo and in vitro, of isolated cells rolling on either monolayers of cultured endothelial cells or surfaces coated with selectin or other molecules. This has helped putting to light the adhesive interaction between rolling cells and the substrate (Springer et al, 1990), and also the stochastic variability of the rolling motion (Schmid-Schöounbein et al, 1987; Goetz et al, 1994). In particular, for cells rolling on a substrate bearing a uniform concentration of adhesion molecules (Springer et al, 1990), such a variability suggests that the fluctuations in the cell dynamics near the vascular wall reflects the stochastic nature of the binding dynamics. Our purpose is to both build and study a model of cell dynamics in a blood vessel to investigate the role of the cell’s adhesion activity.

1.2 Existing Models

There already exist several models and approaches that aim at studying cell adhesion to the vascular wall. The simplest approach consists in ignoring the cell spatial structure and in assuming a step-wise, stop-and-go motion. In the work by Zhao et al (1995), the dynamics of the center of mass of the cell follows a stepping process: its trajectory is approximated by a series of steps made when a cluster of bonds dissociates. Two random variables are used to describe the average stepsize and lifetime of bond clusters resisting the applied fluid force. Performing some mean field approximation makes it possible to heuristically derive a Fokker-Planck equation that governs the cell velocity evolution. This approach allowed to explain the dispersion of rolling velocity data acquired under different experimental conditions. In the same spirit, in the absence of blood flow, macroscopic models have been developed for cell adhesion force (Preziosi and Vitale, 2011), where bonds are described as a distribution function. Its dynamics follows a maturation-rupture equation, also called renewal equation. In the limit of large bonds turnover, a macroscopic friction coefficient can be computed (Milišić and Oelz, 2011, 2015).

In the work by Grec et al (2018), the cell is described as a point carried by the blood flow and interacting with the endothelium, seen as a straight line. At the level of the individual receptor molecule, ligand binding and dissociation are stochastic Poisson processes. The bond forces are described by linear elastic forces. For constant binding and dissociation rates, the averaged problem writes as a deterministic linear Volterra integro-differential equation similar to the one considered by Preziosi and Vitale (2011); Milišić and Oelz (2011), and provides some information on the cell location. Linear continuous models are not satisfactory as they can not describe the strong dependency of the cell arrest on the shear flow. In the work by Grec et al (2018), force-dependent binding rates and/or nonlinear elastic laws are also considered. For these nonlinear models a threshold on the blood velocity, under which the cell velocity vanishes, was numerically observed. Furthermore the link between the nonlinear stochastic processes and the deterministic equations was only highlighted numerically.

In the works by Hammer and Lauffenburger (1987); Hammer and Apte (1992), the authors used a numerical model to describe the interplay between hydrodynamic transport and adhesion. In these works, the cell is modeled as a receptor-bearing hard sphere that interacts with a planar ligand-bearing wall. The ligand-receptor binding follows a chemical kinetic dynamics according to Bell’s law, (Bell, 1978). Bonds then exert elastic forces on the sphere while the linear shear flow exerts both hydrodynamic force and torque. In the PhD thesis of Korn (2007), the Brownian motion of the sphere is taken into account in order to model the spatial receptor-ligand encountering more precisely. These models allow the numerical study of the influence of weak bonds on the motion of the sphere. More recently, other mechanical models were built to describe more precisely the hydrodynamic forces exerted on the cell (Reboux et al, 2007; Efremov and Cao, 2011; Li et al, 2018). These works justify a bistable behaviour between either cell release and rolling, or between cell arrest and rolling.

1.3 Results of the paper

The main goal of this article is to study the balance between shear forces and resistive adhesion forces, and its consequences on the dynamics of a cell located near the arterial wall. To this aim, we build a stochastic model of birth-death type for a particle cell that develops adhesive interaction with the vascular wall. Note that in this setting, no information on the repartition of bonds can be obtained, so that cell rolling can not be described and will thus be confounded with cell sliding in the following. In order to model the effect of the cell motion on the adhesion dynamics, we consider cell velocity-dependent binding rates. In such a case the model is nonlinear. We want to capture the transition from cell rolling to either its stopping or its release in the blood flow. For this purpose, we perform and justify some scaling limits to quantify for the effect of the cell adhesion activity on the dynamics of the cell.

Following the approach of Grec et al (2018), we consider a minimal discrete stochastic model where the cell is a point particle carried by the blood flow in a one-dimensional setting. The adhesion dynamics is modeled by a Markovian Jump process for the formation and disassembly of bonds exerting units of resistive force. The choice of a stochastic model follows biological observations of e.g Springer et al (1990). This model has some similarities with the one heuristically derived by Zhao et al (1995) but describes more precisely the interaction with the endothelium, in the spirit of Hammer and Apte (1992). Moreover, by considering constant binding forces instead of elastic ones, the present model differs from the one given by Grec et al (2018). Additionally, we take into account both the adhesions growth, modelled by the reproduction of bonds, and the effect of the cell velocity on the binding dynamics, thus leading to a nonlinear stochastic process.

The timescale of the binding dynamics being actually fast compared to the timescale of cell motion, it is natural to use a time-continuous description of the adhesion dynamics. Therefore, in a second step, we rescale the process in the spirit of Joffe and Metivier (1986); Ethier and Kurtz (2009), and we rigorously derive continuous limiting models for the cellular adhesion dynamics. Depending on the renormalization assumptions, we obtain either a deterministic or a stochastic model, which we both study.

The deterministic model predicts that when the shear rate is high or the vascular wall is in a lowly inflamed state (with a low density of adhesion proteins), the cell develops no bonds with the wall, while it can slow down and slide on the wall otherwise, until it eventually stops and adheres to the wall. More generally, the model successfully predicts bistable behaviours in some parameter spaces, between either cell release and sliding or cell sliding and arrest. Consequently, the model is indeed able to describe the two-step process involved in the cell adhesion to the vascular wall. The analysis of the continuous stochastic model is made in two steps. First, when there is no feedback from the cell velocity on the adhesion activity, the model writes as a Cox-Ingersoll-Ross (CIR) process, for which the probability density of the arrest time is explicitly given and numerically shows a transition between the cell stopping and its release in the bloodstream. Then, in the nonlinear case, when the blood flow exerts a feedback on the binding dynamics, we derive the cell mean stopping time using a Fokker-Planck approach.

We believe that this work can have strong implications for the immune response, drug delivery systems, as well as tumor invasion. More precisely, our model could be used as a first step in the construction of a permeability law for the vascular wall. In a future work, it will be compared with experimental in vivo data obtained by Follain et al (2017) in the context of metastatic invasion.

The plan of this article is the following. In Section 2, we detail the construction of the discrete stochastic model of the individual bond dynamics, and we perform its analysis together with numerical simulations. In Section 3, we proceed to rigorous derivations and study of models of the continuous time-evolution of the cell velocity, which are either deterministic or stochastic. In Section 4 we discuss our results and how they allow the characterization of the long term cell motion.

2 A Markovian Jump process for the cell adhesion dynamics

We first present the stochastic model used to investigate the adhesion process outlined in the introduction. We use a classical birth-and-death process for bonds turnover. As in most of the literature, see Bell (1978), it is assumed that the bond lifetime has an exponential distribution. Then, we study its mathematical properties and perform some numerical simulations.

2.1 Modelling approach

Let us consider a cell carried by the blood flow. We suppose that the distance between the cell and the blood vessel wall is small enough so that bonds may always form. Since the cell is in the vicinity of the wall, we assume that the blood shear flow is 1D, parallel to the vascular wall and with a constant velocity, denoted by $u\in\mathbb{R}_{+}$.

Deformability may play a role in the cell dynamics inside the blood flow and in its interaction with the vessel wall. However, some large cells like Circulating Tumor Cells can be very stiff (Follain et al, 2017), and still follow the behaviour studied in this paper. Therefore, we hereafter choose to neglect cell deformability and to focus on the interaction between the intracellular adhesion dynamics and the vascular wall.

In previous studies, see Chang and Hammer (1999); Jadhav et al (2005) e.g., it was shown that approximating the contact surface by a simple geometrical figure (a circle or a rectangle) and neglecting the increase of the contact surface with the flow shear rate due to cell deformability do not qualitatively change the analysis. Moreover, as suggested in Bhatia et al (2003), the cell adhesion is primarily determined by physicochemical properties of adhesion proteins. In a first approximation, we thus assume the cell to be a point particle whose position at time $t\geq 0$ is denoted by $X_{t}$.

Velocity model

To describe the cell motion, we use a non-inertial approximation. Indeed, in a regime of low Reynolds number, viscous forces outweight inertial forces and the momentum equation reduces to the force balance principle:

[TABLE]

where $V_{t}\in\mathbb{R}$ is the cell velocity, $u$ is the blood shear flow and the cell is subjected to a macroscopic resistive force, denoted by $F_{t}\in\mathbb{R}_{+}$, and induced by the bonds that contribute to decelerating the cell, see Figure 1. The parameter $\gamma$ is such that $\gamma^{-1}$ is a friction coefficient, following a linear force-velocity relation. The previous equation is valid only for $\gamma F_{t}\leq u$, as for a maximal force the cell stops and the model is no longer valid. The resistive force arises from the strength of the cell adhesion to the vessel wall. Cellular adhesion is the macroscopic readout of the forces exerted by the wall on the cell through each bond (Boettiger, 2007). As a result, we assume that

[TABLE]

where $f\geq 0$ is the typical force generated by a stabilized bond, and $N_{t}$ is the number of stabilized bonds at time $t$. Note that $f$ ranges in $\SIunits{\mathrm{p}\mathrm{N}}$ (Boettiger, 2007), but its precise value depends on the experimental conditions.

Non-dimensionalization

We now introduce typical quantities for our problem: the typical force considered is the one generated by a stabilized bond: $\overline{F}=f$ ; the associated typical velocity writes $\overline{V}=\gamma\overline{F}=\gamma f$ ; and the typical timescale of the model is taken to be the typical lifetime of a bond (Alon et al, 1995): $\overline{T}=1\SIunits{\mathrm{s}}$. Therefore, keeping the same notations for simplicity, the nondimensionalized problem writes

[TABLE]

We now construct the process $(N_{t})_{t}$ of the number of stabilized bonds over time.

Stochastic model for the adhesive force

Let us now introduce the simple discrete model that we use to describe the individual bonds dynamics. We write $(N_{t})_{t}$ the Markovian processes for the number of stabilized bonds at time $t$, which follows a classical birth-and-death-like dynamics in the state space $\mathbb{N}$.

• $\bullet$

New bonds form spontaneously at rate $c(u)=c\mathds{1}_{u\leq u^{*}}$, for a velocity threshold $u^{*}$ above which no new bonds can be created, due to the high blood velocity.

• $\bullet$

Each existing bond can reproduce at constant rate $r$. This phenomenon reflects the local reinforcement of the connection to the vessel wall by the involvement of integrins in adhesion growth, which can be imputed to cytoskeletal forces or external stresses (Nicolas et al, 2004). Moreover, intuitively, if an adhesion complex is composed of a large number of bonds, the unbounded molecules can find an attachment more easily than a less stable adhesion formed of fewer molecules.

• $\bullet$

Each bond dissociates at the velocity-dependent rate $d(V_{t})=de^{\alpha V_{t}}=de^{\alpha(u-N_{t})}$. We choose here an exponential relation, where $d$ is the unstressed-bonds dissociation rate, and $\alpha$ is a sensitivity parameter. This choice accounts for the fact that the average lifetime of an adhesion site changes with the applied tension exerted by the blood flow: the faster the cell, the shorter the bonds’ lifetime. Note that since the cell velocity is bounded by $u$, the dissociation rate is also bounded. In the following, we will write indifferently $d(V_{t})$ or $d(N_{t})$.

Remark 1

The rate of a single bond formation is mainly determined by the time spent by the two proteins near one another. Therefore, the rate $c$ should depend on the cell velocity when it is non zero. A more realistic choice for $c$ would be a decreasing function of the instantaneous cell velocity $V_{t}$ such as $c(v)=\left(u_{*}-v\right)_{+}$, where $\left(\cdot\right)_{+}$ denotes the positive part function. Since $V_{t}$ depends on $N_{t}$, the rate rewrites $c(n)=\left(u_{*}-u+n\right)_{+}$ with $n$ the number of stabilized bonds. For such a choice, and assuming that $v\leq u_{*}$, the dependency on $n$ then only provides an additional contribution to the reproduction rate.

Note that these rates are also representative of the adhesive properties of the endothelial cells forming the vessel wall. The key point here is that there is a feedback loop between the instantaneous cell velocity and the bonds dynamics. More elaborate dependency could be considered (see e.g Milišić and Oelz (2011, 2015)), in particular involving age dependencies to model the bond elasticity (see Grec et al (2018)), but we choose to keep a minimal set of parameters, for both simplicity and the sake of clarity.

The balance between the adhesion force on the one hand, and the load and torque created by the blood flow on the other hand then determines the outcome of the dynamics: either the cell rolling and arrest or its release in the blood flow.

2.2 Mathematical properties of the discrete model

In this section we derive some mathematical properties based on stochastic analysis tools, and we perform numerical simulations of the process $(N_{t})_{t}$. The dissociation rate being nonlinear, classical tools do not apply. Since we are interested in the dynamics while $V_{t}\geq 0$, that is, while $N_{t}\leq u$, we define the stopping time

[TABLE]

We are interested in the Markovian jump process $(N_{t})_{t\in[0,\tau_{u}]}$ defined by the following transitions:

[TABLE]

where $\lambda$ and $\mu$ are defined on $\mathbb{N}$, and are bounded. Such a process is classically well-defined (see e.g Ethier and Kurtz (2009)). We also control the mean number of bonds in finite time.

Proposition 1 (Moments propagation)

Assume that there exists $p\geq 1$ such that $\mathbb{E}\left[N_{0}^{p}\right]<+\infty$. Then,

[TABLE]

Proof

This Proposition is proved in a more general framework in Appendix A.

The mean path of this process can not be fully studied in the general case. Indeed, for $\mathbb{E}\left[N_{0}\right]<+\infty$, we classically write the mean equation

[TABLE]

where, even when assuming $t\leq\tau_{u}$, the nonlinearity prevents any analysis. In a simpler case, namely when there is no feedback from the cell velocity on the adhesion dynamics, and without considering the stopping time $\tau_{u}$, we obtain a classical immigration-birth-death process, that was already studied by Bansaye and Méléard (2015); Dessalles et al (2018). More precisely, $(N_{t})_{t}$ then follows a negative binomial distribution of parameters $\left(\frac{c(u)}{r},\frac{r}{d}\right)$. It follows that

[TABLE]

At steady state, one finds

[TABLE]

with

[TABLE]

In the case of a circulating cell, we are only interested in the situation where $v\leq u$. As a consequence, assuming $\mathbb{E}[N_{0}]=0$ and $u>0$, we obtain the following mean asymptotic behaviours:

[TABLE]

These results show that the birth-and-death dynamics without feedback intrinsiquely carries a dichotomic asymptotic behaviour. However, such a model is not fully satisfactory, since the only shear-treshold effect that appears concerns the initiation of the adhesive interaction and not the cell fate.

2.3 Numerical Simulations

The process $(N_{t})_{t}$ being Markovian, it can be simulated directly events after events. Consider the population size $N_{T_{k}}$ at time $T_{k}$. Then,

• $\bullet$

the global jump rate writes $\varsigma_{k}=\lambda(N_{T_{k}})+\mu(N_{T_{k}})$, which means that the time before the next event is a random variable following an exponential law of parameter $\varsigma_{k}$. A realization of this law provides $T_{k+1}$.

• $\bullet$

a new bond is created with probability $\lambda(N_{T_{k}})/\varsigma_{k}$, while with probability $\mu(N_{T_{k}})/\varsigma_{k}$ a randomly chosen bond disassembles, and $N_{T_{k+1}}$ follows.

The above procedure can be iterated to give the time evolution of the process. Numerical simulations of the process are displayed in Figure 2. One may observe that the velocity may either shrink to zero or remain close to $u$ for the same parameter values. Note also that sliding phases are observed in both cases.

3 Continuous limiting models and characterization of the dynamics

In this section we separate the scale of the adhesion dynamics from the one of the cell motion. Such scale separation is justified by the large number of bonds, and by the very fast binding dynamics, as compared to the cell displacement. As an illustration, a bond lifetime ranges around $1\SIunits{\mathrm{s}}$, a typical binding rate ranges aroung $10^{3}$ $\SIunits{\reciprocal\mathrm{s}}$ whereas the cell rolling velocity ranges around $30\SIunits{\mathrm{\SIUnitSymbolMicro}\mathrm{m}\reciprocal\mathrm{s}}$ (see Grec et al (2018) and the references therein). As will be seen next, this assumption allows the use of a scaling approach to derive two continuous limiting models, for which deeper analysis can be pursued.

More precisely, let $K\geq 1$ be a parameter such that $1/K$ scales the force generated by a fraction of bonds. Moreover, we assume that the bond dynamics gets faster and faster as the fraction of bonds gets smaller: we consider now $K$-dependent rates $c^{K},\;r^{K},$ and $d^{K}$, related to the process $(N_{t}^{K})_{t}$. We then define the renormalized process $(Z_{t}^{K})_{t}$ by

[TABLE]

3.1 Deterministic continuous limiting model

We consider the following rates:

[TABLE]

In other words, in considering an approximation for a continuous description of bonds, the formation rate is larger, while the self-enhancement of the adhesion dynamics and the bonds typical lifetime stay unchanged. From the modelling viewpoint, such an assumption amounts to considering adhesion clusters that involve only a small number of proteins on each side (wall and cell). Note that since the reproduction is not accelerated, the clustering that leads to adhesion growth is not large enough to induce stochastic fluctuations at the cell level. In this context, we obtain the following convergence result.

Theorem 3.1

Consider the sequence of processes $(Z^{K})_{K}$ for $(Z^{K}_{t})_{t\geq 0}$ defined by (4) and rates defined by (5). If $Z_{0}^{K}\underset{K\rightarrow+\infty}{\longrightarrow}n_{0}\in\mathbb{R}_{+}$ in probability, and if

[TABLE]

then, for $T>0$, $(Z^{K})_{K>0}$ converges in law in $\mathbb{D}\left([0,T],\mathbb{R}_{+}\right)$ to the unique continuous function $n\in\mathcal{C}([0,T],\mathbb{R}_{+})$ solution to

[TABLE]

Remark 2

By the Gronwall lemma, one has for $T<\infty$,

[TABLE]

showing that the global density stays finite in finite time.

Proof

The proof is displayed in Appendix B.

We now perform the analysis of the limiting problem. Let us define the function $F$ by $F(n)=c\mathds{1}_{u\leq u_{*}}+\left(r-de^{\alpha(u-n)}\right)n$. We prove the following result.

Proposition 2

Assume that the rates are given by (5). Then the stationary state(s) $n^{\infty}$ of (6) are as follows.

If $u>u_{*}$, then the system admits two stationary states $n^{\infty}_{1}=0$ and $n^{\infty}_{2}=u-\frac{1}{\alpha}\ln(\frac{r}{d})$. The smallest is stable and the largest is unstable. 2. 2.

If $u\leq u_{*}$,

1. (a)

for $u\leq\frac{1}{\alpha}\ln\left(\frac{r}{d}\right)$, then $n^{\infty}=+\infty$. 2. (b)

for $u>\frac{1}{\alpha}\ln\left(\frac{r}{d}\right)$, then there exists a unique $0<\mathrm{bar}n<\frac{1}{\alpha}$ such that $F^{\prime}(\mathrm{bar}{n})=0$ and

1. i.

If $F(\mathrm{bar}{n})>0$, then $n^{\infty}=+\infty$. 2. ii.

If $F(\mathrm{bar}{n})=0$, then $\mathrm{bar}{n}$ is the unique stationary solution. 3. iii.

If $F(\mathrm{bar}{n})<0$, then there exists two stationary solutions $n^{\infty}_{1}$ and $n^{\infty}_{2}$, such that $0<n^{\infty}_{1}<\mathrm{bar}{n}<n^{\infty}_{2}<+\infty$, the smallest being stable and the largest unstable.

Proof

The case $u>u_{*}$ follows from a direct computation. Consider the case where $u\leq u_{*}$, then one has

[TABLE]

To understand the dynamics of the velocity of the cell, we have to understand that of $n$. We are therefore interested in the values of $n$ for which $F(n)=0$, which are the steady states, and in the sign of $F$ required to obtain the stability of the steady states. A quick computation shows that

[TABLE]

We can study the sign of $F^{\prime\prime}(n)$ to obtain the variations of $F^{\prime}$. This yields the following variation table:

{variations}

n & \z 2α \pI

\filet

F”(n) + \z -

\filet

F’(n) r -de^αu ¸ \hF’(2α)>0 . r

Since $F^{\prime}(2/\alpha)>0$ and $r>0$, the existence of a stationary state $n^{\infty}$ such that $F(n^{\infty})=0$ depends on the sign of $F^{\prime}(0)=r-de^{\alpha u}$. As a consequence,

• •

If $u\leq\frac{1}{\alpha}\ln\left(\frac{r}{d}\right)$, then $\forall n\in\mathbb{R}_{+}$, $F(n)\geq c>0$, hence $n^{\infty}=+\infty$.

• •

If $u>\frac{1}{\alpha}\ln\left(\frac{r}{d}\right)$, then $F^{\prime}(0)<0$, and there exists a unique $0<\mathrm{bar}n<\frac{2}{\alpha}$ such that $F^{\prime}(\mathrm{bar}{n})=0$. This provides the sign of $F^{\prime}$, from which we obtain the variation table for $F$ that leads to the result:

{variations}

n & \z bar n <2α \pI

\filet

F’(n) - \z +

\filet

F(n) \hc>0 . F( bar n) ¸ \h\pI

Notice that since $u>\frac{1}{\alpha}\ln\left(\frac{r}{d}\right)\Leftrightarrow de^{\alpha u}>r$, both behaviours arise according to the comparison between the reproduction and death rates. Note also that since $F^{\prime}$ is strictly increasing on $\left(0,\frac{2}{\alpha}\right)$ and since $F^{\prime}(1/\alpha)=r>0$, we obtain that $\overline{n}<1/\alpha$.

Remark 3

Note that if $u_{*}\geq u>\frac{1}{\alpha}\ln\left(\frac{r}{d}\right)$ the three cases described in the proposition above may occur. Indeed, consider the particular case where $d>r$ and $u=\frac{1-\frac{r}{d}}{\alpha}$, then $\mathrm{bar}n=u$ and $F(\mathrm{bar}n)=c-\frac{(d-r)^{2}}{\alpha d}$ whose sign depends on the value of $c$. The dynamics is then dependent on the ability of the cell to form bonds at primary contact.

We are now interested in the biological interpretation of this study. Since our model is only valid up until the cell adhesion to the vessel wall, we consider Equation (6) up until the adhesion density reaches $u$. The following Corollary locates $u$ with respect to the stationary state(s) of the system. This allows us to assess the cell fate depending on the parameter values. For $\alpha>0$, denote the key parameters

[TABLE]

Corollary 1

Let $n_{0}\in[0,u]$ with $u>0$. Then, Equation (6) admits either one stationary state denoted by $n^{\infty}$, or two stationary states $n_{1,2}^{\infty}$ such that $0<n_{1}^{\infty}<n_{2}^{\infty}$, the smaller one being stable and the larger unstable.

If $u>u_{*}$:

1. (a)

for $r\leq d$, then $n^{\infty}_{1}=0$ and $n^{\infty}_{2}=u-U_{\alpha}\geq u$. 2. (b)

for $r>d$, if $u\leq U_{\alpha}$, $n^{\infty}=+\infty$ ; if $u>U_{\alpha}$, one has $n^{\infty}_{1}=0<n^{\infty}_{2}=u-U_{\alpha}<u$. 2. 2.

If $u\leq u_{*}$, then

1. (a)

for $r<d$, $\exists!0<\mathrm{bar}n<\frac{1}{\alpha}$ such that $F^{\prime}(\mathrm{bar}{n})=0$.

1. i.

For $u>\overline{U}_{\alpha}$, $F^{\prime}(u)>0$, so that $u>\overline{n}$.

1. A.

If $u>U_{c}$, then we have $0<n_{1}^{\infty}<\overline{n}<u<n_{2}^{\infty}$. 2. B.

If $u=U_{c}$, then $0<n_{1}^{\infty}<n_{2}^{\infty}=u$. 3. C.

If $u<U_{c}$, then if $F(\mathrm{bar}n)<0$, one has $0<n_{1}^{\infty}<\overline{n}<n_{2}^{\infty}<u$ ; if $F(\mathrm{bar}n)=0$, then $n^{\infty}=\overline{n}<u$ ; if $F(\mathrm{bar}n)>0$, then $n^{\infty}=+\infty$. 2. ii.

For $u=\overline{U}_{\alpha}$, then $u=\overline{n}$.

1. A.

If $c>\overline{C}$, then $n^{\infty}=+\infty$. 2. B.

If $c=\overline{C}$, then $n^{\infty}=u$. 3. C.

If $c<\overline{C}$, then $0<n_{1}^{\infty}<u<n_{2}^{\infty}$. 3. iii.

For $u<\overline{U}_{\alpha}$, then $0<u<\overline{n}$.

1. A.

If $u=U_{c}$, then $n_{1}^{\infty}=u<n_{2}^{\infty}$. 2. B.

Otherwise, if $F(\mathrm{bar}n)=0$, then $u<n^{\infty}=\overline{n}$ ; if $F(\mathrm{bar}n)>0$, then $n^{\infty}=+\infty$. Finally, if $F(\mathrm{bar}n)<0$, then $u<n_{1}^{\infty}<n_{2}^{\infty}$ when $u<U_{c}$, and $n_{1}^{\infty}<u<\overline{n}$ when $u>U_{c}$. 2. (b)

For $r=d$, there exists a unique $0<\mathrm{bar}n<\frac{1}{\alpha}$ such that $F^{\prime}(\mathrm{bar}{n})=0$ and $u>\overline{n}$. If $F(\mathrm{bar}n)<0$, we have $0<n_{1}^{\infty}<\overline{n}<n_{2}^{\infty}<u$. If $F(\mathrm{bar}n)=0$, $n^{\infty}=\overline{n}<u$. If $F(\mathrm{bar}n)>0$, $n^{\infty}=+\infty$. 3. (c)

For $r>d$, $F(u)>0$.

1. i.

if $u\leq U_{\alpha}$, then $n^{\infty}=+\infty$. 2. ii.

if $u>U_{\alpha}$, then $\exists!0<\mathrm{bar}n<\frac{1}{\alpha}$ such that $F^{\prime}(\mathrm{bar}{n})=0$ and $u>\overline{n}$. Then, if $F(\mathrm{bar}n)<0$, we have $0<n_{1}^{\infty}<\overline{n}<n_{2}^{\infty}<u$. If $F(\mathrm{bar}n)=0$, $n^{\infty}=\overline{n}<u$. If $F(\mathrm{bar}n)>0$, $n^{\infty}=+\infty$.

Proof

The Corollary results from Proposition 2 combined with the sign analysis of $F(u)$ and $F^{\prime}(u)$ that depends on parameter values.

Let us give an interpretation of the above results. First, note that when there are two stationary states, the largest one denoted by $n_{2}^{\infty}$ is unstable. Therefore, if $n_{2}^{\infty}<u$, it follows that for $n(0)>n_{2}^{\infty}$, the adhesion density will increase until the cell adheres to the wall. This can be interpreted as another stable stationary state for the model, which therefore captures the bistability previously described in experimental and theoretical studies. This property will be further discussed in Section 4.

Table 3a gives the outcome of the cell dynamics in the case where $u>u_{*}$. In this case, the blood flow is so fast that the cell does not initiate any new adhesion. However, the dynamics may still be interesting in the case where some bonds already exist and may be stabilized by the growth of adhesions. Then, if $r\leq d$, dissociation is always more frequent than strenghtening of existing adhesions, so that starting with $n_{0}\leq u$, the cell is released into the bloodstream. When $r>d$, and depending on the sensitivity $\alpha$ of the dissociation rate on the cell velocity, the cell either goes back in the bloodstream or adheres firmly to the wall, while sliding is an unstable state.

Table 3b shows the possible outcomes in the case where $u\leq u_{*}$. In this situation, the whole adhesion dynamics is active. Then, the conditions discriminating between different cell fates are based on the balance between the formation of adhesion bonds (related to $r$ and $c$), and their dissociation (related to $d$ and $\alpha$). This study shows how Equation (6) has more complex behaviours in comparison with the mean linear ODE (3).

Let us comment on these results. First of all, not surprisingly, our model successfully predicts the shear-stress threshold above which nor capture nor sliding occurs when the cell does not initially interact with the wall. This is explained by the regulation of bonds dissociation by the cell velocity. The model also predicts the firm adhesion to the wall that is observed experimentally. Finally, the model provides conditions on the parameter space for both stable sliding and bistable behaviour between cell sliding and cell adhesion. This results come from the competition between bonds formation and rupture that occurs in the cell-wall contact area.

3.2 Stochastic continuous limiting model

In this section, we consider the following rates:

[TABLE]

with $K>0$ and $a>0$. The whole adhesion dynamics is therefore accelerated. Note that using the same acceleration for reproduction and death permits to keep the same bounded individual growth rate $r^{K}-d^{K}=r-d$. This way, even if each adhesion bond reproduces and dies infinitely faster, its contribution to the global adhesion growth remains the same.

Theorem 3.2

Consider the sequence of processes $(Z^{K})_{K}$ for $(Z^{K}_{t})_{t\geq 0}$ defined by (4) and rates defined by (7). If for $K\rightarrow+\infty$ the initial value $Z_{0}^{K}$ converges in law to a $\mathbb{R}_{+}$-valued random variable $N_{0}$, with

[TABLE]

then $(Z^{K})_{K>0}$ converges in law in $\mathbb{D}\left([0,T],\mathbb{R}_{+}\right)$ to the continuous process $N=(N_{t})_{t\in[0,T]}\in\mathcal{C}([0,T],\mathbb{R}_{+})$ solution of

[TABLE]

with $B_{t}$ a Brownian Motion, $b(N_{t})=c(u)+(r-d(N_{t}))N_{t}$ and $\sigma(N_{t})=\sqrt{2aN_{t}}$.

Proof

The proof is displayed in Appendix C.

Remark 4

The solution to the SDE (8) is almost surely positive if $b(n)\geq 0$ for all $n\geq 0$, and for a positive initial state (see the 1D comparison principle in Karatzas and Shreve (1998) e.g.).

Numerical simulations

In order to preserve the positivity of the process, we perform numerical simulations of the SDE (8) using a symmetrized Euler scheme, which consists in taking the absolute value of the classical Euler scheme (see e.g Berkaoui, Abdel et al (2008)). More precisely, the scheme is the following: write $(N_{k})_{k}$ for the discretization of $(N_{t})_{t}$, where $N_{k}$ corresponds to the time $t_{k}=k\Delta t$. Then, define $N_{0}=n_{0}$ and for $k\geq 0$,

[TABLE]

with $W\sim\mathcal{N}(0,1)$. The strong $L^{1}$ convergence of this scheme is proved by Berkaoui, Abdel et al (2008) provided that

[TABLE]

for $P\geq|r-d|$ and $\Delta t\leq\frac{1}{2P}$. This condition allows to deal with the non-Lipschitz diffusion coefficient, and rewrites, in our case,

[TABLE]

As an example, $\frac{a}{4}\left(\frac{c}{a}-1\right)^{2}>8a$ is equivalent to $\left(\frac{c}{a}-1\right)^{2}>32$, which is verified for $c>7a$. The numerical simulations are displayed in Figure 4, and show two typical arrest and release situations.

The rolling motion of individual cells has been observed to fluctuate randomly both in vivo and in vitro. The model (8) obtained in Theorem 3.2 is therefore suited to investigate the effects of the stochastic fluctuations associated with the adhesion activity on the cell dynamics. In particular, we are interested in the expected time until the cell firmly adheres to the vessel wall.

The case without feedback: the CIR process.

When we assume that the cell velocity exerts no feedback on the bonds dissociation, the solution of (8) reduces to a CIR process, see Göing-Jaeschke and Yor (1999); Shreve (2004); Ikeda and Watanabe (1989):

[TABLE]

with $c>0$, $a>0$ and $r-d\in\mathbb{R}$. It is known that such processes arise from the diffusion limit of discrete branching processes with immigration (Athreya and Ney, 2004), and show a dichotomy behaviour. Depending on the parameters, the density is almost surely either close to [math] or large, leading to the almost sure cell arrest in our model. Simulations of this process are displayed in Figure 5. Some general properties of the CIR process are displayed in Appendix D. In particular, its stationary probability density is represented in Figure 5 and shows the transition between both behaviours.

Time to reach $u$:

It is also possible to obtain information on the time needed to reach a given value. More precisely, one can use the Laplace transform of the first hitting time of any value, starting at a given point (Göing-Jaeschke and Yor, 1999; Leblanc and Scaillet, 1998). However, it is not possible to proceed to its inversion analytically. While numerical inversions procedures exist (some of them are compared by Leblanc and Scaillet (1998)), they do not always provide satisfactory results: the integral of the output may not be equal to one, and negative values may appear. The procedure proposed by Abate and Whitt (1992) seems satisfactory in this viewpoint.

In this paper, we follow the work of Linetsky (2004a, b) to numerically compute the first hitting time density using an eigenfunction decomposition, an approach used for diffusions in the litterature (Davydov and Linetsky, 2003; Itô and McKean, 1965; McKean, 1956). For the CIR process, it is established in Linetsky (2004c, b) that the same type of decomposition holds. We now recall the result of Linetsky (2004a) that provides a series expansion for the density $f_{T_{x\rightarrow y}}$ of the first hitting time of $y$ starting from $x$.

Proposition 3

Linetsky (2004a)**

1. i)

For $0<x<y\in I$, and $t>0$ we have

[TABLE]

with uniform convergence on $[t_{0},+\infty)$, $t_{0}>0$, and $(\lambda_{n})_{n}$ a strictly positive and strictly increasing sequence with $\lambda_{n}$ growing to $+\infty$ as $n$ goes to infinity. More precisely, we have that

[TABLE]

with $(s_{n})_{n}$ the strictly decreasing sequence of strictly negative roots of $\Phi(\cdot;c/a;\overline{y})=0$, with $\Phi(w_{1};w_{2};w_{3})$ denotes the Kummer confluent hypergeometric function. The sequence $(o_{n})_{n}$ is defined by

[TABLE]

for $\overline{y}:=-\frac{r-d}{a}y$ and $\overline{x}:=-\frac{r-d}{a}x$. 2. ii)

Moreover, the following asymptotics hold:

[TABLE]

as well as

[TABLE]

Therefore, the proposed numerical method requires the computation of the set of negative roots of $\Phi$ to get approximations of the families $\{\lambda_{n}\}_{n}$ and $\{o_{n}\}_{n}$. The choice of the level of truncation for the approximation of (10) may be made using the following estimate:

[TABLE]

for

[TABLE]

Linetsky also notices that using (13)-(14) instead of computing zeros of the Kummer function provides satisfactory results, in particular when $c/a$ is small. For better accuracy, one can also use the exact expressions (11)-(12) of the first term of the decomposition, then the estimates (13)- (14) for the others. The following numerical simulation was performed using the asymptotic expansions of $\lambda_{n}$ and $o_{n}$ only, even for $n$ small, since it is observed that this approximation does not change qualitatively the profile (see Figure 6). Obviously, in this case, the obtained function is not a probability density, and an inconsistency tends to appear near $t=0$ due to the approximation, but the overall shape is preserved.

The general case

Let us now focus on the general case of Equation (8). As a first approach one can use the 1D comparison principle (see e.g Revuz and Yor (2005)) to compare the process with CIR processes. In this work, we follow another method and derive from (8) a Fokker-Planck equation on $p(n,t):=p(n,t|n_{0},t_{0})$ the probability density of $(N_{t})_{t}$ conditionally to its initial condition. We obtain the following equation:

[TABLE]

where we recall that $b(n)=c+(r(n)-d(n))n$, while $\sigma(n)=\sqrt{2an}$ and $\mathbf{J}(n,t)$ is the associated probability current. The natural boundary conditions are the following:

[TABLE]

We are interested in the mean time necessary for the process to reach the value $u$ starting from $n_{0}\in(0,u)$, that we denote by $\tau_{u}(n_{0})$. This question can be adressed by considering the Fokker-Planck equation on $(0,u)$ with [math] a reflecting and $u$ an absorbing barrier. We show the following proposition.

Proposition 4

The mean time $\tau_{u}(n_{0})$ necessary for the process to reach the value $u$ starting from $n_{0}\in(0,u)$ writes

[TABLE]

Proof

Write $G(n_{0},t)$ the probability that a particle starting at $n_{0}$ is still in $(0,u)$ at time $t$. Then,

[TABLE]

Since the dynamics is homogeneous in time, we deduce that $p(n,t|n_{0},0)=p(n,0|n_{0},-t)$ and $n_{0}\mapsto p(n,t|n_{0},0)$ satisfy the backward Fokker-Planck equation:

[TABLE]

and $(n_{0},t)\mapsto G(n_{0},t)$ satisfies

[TABLE]

The initial and boundary conditions are the following:

[TABLE]

Take $f\in\mathcal{C}^{1}(\mathbb{R},\mathbb{R}_{+})$ non-decreasing. Then, classically, $\mathbb{E}[f(\tau_{u})]=\int_{0}^{+\infty}f^{\prime}(t)\mathbb{P}(\tau_{u}>t)\,{\mathrm{d}}t=\int_{0}^{+\infty}f^{\prime}(t)G(n_{0},t)\,{\mathrm{d}}t$. Hence, we get for $k>1$,

[TABLE]

Integration of (16) in time leads to the following ODEs on the family $(\tau_{k})_{k\geq 1}$:

[TABLE]

and for $k>1$,

[TABLE]

By direct integration, we can solve (17), allowing to solve successively the problems (18). Write

[TABLE]

Then, we have

[TABLE]

In practice, denoting $\epsilon>0$ the lower bound in the integral instead of zero, we find that

[TABLE]

Explicit computations lead to the result.

We perform numerical simulations of $\tau_{u}(0)$ for different values of the blood flow velocity $u$ (see Figure 7). Below a shear-velocity treshold, the particle stops very quickly, while it becomes extremely slow above the treshold. The numerical phase plane shows a natural dependency on the adhesion creation rate $c$. The value of $a$, quantifying the noise intensity, has a direct effect on the range of the arrest time, but does not qualitatively change the phase plane.

4 Discussion

In this work, we have presented a discrete model of cell adhesion on a vessel wall. The model is based on the stochastic formation of both weak bonds between the cell and the wall, and stronger ones arising by self-reinforcement. This phenomenon is modelled by a stochastic birth-and-death-type process. The cell velocity exerts a feedback on the breaking rate of bonds: the faster the cell, the shorter-lived the bonds.

Our purpose is to explore the events that follow the first cell contact with the wall and that may determine whether the cell stops and adheres or goes back into the blood stream. More precisely, we are interested in reproducing the bistable behaviour observed experimentally between cell rolling on the wall and its stationary adhesion. Moreover, we look for an expression of the cell stopping time, which is equivalently the first hitting time of a treshold value for the density of bonds. To achieve these goals, we perform some scaling limits to derive continuous deterministic and stochastic models that allow for deeper theoretical analysis.

4.1 Implication of our work on the understanding of the arrest and extravasation of circulating cells

In vivo, the firm adhesion of circulating cells on the wall is a first step towards its extravasation out of the vessel, which has different consequences depending on the cell type. Extravasation of Circulating Tumor Cells (CTCs) takes part in tumor invasion since it allows the formation of secondary tumors (Follain et al, 2017). During the immune response, leukocytes carried by the bloodstream firmly adhere to the vessel wall at inflammed sites. Their extravasation then allows them to pursue their immune function. Let us also mention that cell adhesion to vessel walls is involved in the study of drug delivery systems.

Overall, cell adhesion to vessel walls is a phenomenon showing major applications in biology and medecine. It is now clear that understanding the processes involved in the determination of the location of cell arrest is of prime importance. This justifies the development of new mathematical models able to explain the experimental observations. Our model focuses on the cell behaviour within the bloodstream, and is consequently no longer valid after the cell velocity has reached zero. A different model is then required to capture the fundamental processes involved in extravasation.

4.2 The linear mean ODE (3) does not agree with known observations

Our first finding concerns the mean linear Ordinary Differential Equation (3) in which the cell velocity has no effect on the adhesion dynamics: ligand binding and dissociation rates are constant. As expected, in such a case, the stationary mean number of bonds is independent of the shear rate $u$. Such a result does not agree with in vitro and in vivo experiments, (Follain et al, 2017; Goetz et al, 1994), which show that cells preferably stop in regions with low shear rate. This shows that a nonlinear model is required to capture this phenomenon.

4.3 The nonlinear model qualitatively agrees with known observations: hemodynamic forces affect the adhesion dynamics of circulating cells (Follain et al, 2017)

In the nonlinear case, we perform a renormalization procedure based on the biological observation that the orders of magnitude between the adhesion dynamics and the cell motion are different. This allows us to rigorously build continuous models, either deterministic or stochastic, of the dynamics of the density of bonds (or equivalently of the cell velocity). Then, the analysis shows that these nonlinear models predict the cell adhesion bistability as a result of the competition taking place in the cell-wall contact area between bond formation and rupture.

The deterministic limiting model provides a parameter space for cell adhesion bistability.

The deterministic limiting model writes as a nonlinear ODE, for which we give the stationary states and their stability. This allows us to outline the existence of a dichotomy behaviour in which the stability of the cell arrest is explicitly related to the balance between the blood flow velocity and the adhesion dynamics. It leads to the identification of a shear stress value that separates the cell arrest from a moving state, in agreement with experimental observations (Follain et al, 2017).

Our study also provides a parameter space associated with the existence of two stationary states, with a stable release state and an unstable sliding one. Moreover, since the model is no longer valid after the bond density reaches $u$, the cell arrest is an additional «stable» stationary state. Therefore, the model is indeed able to reproduce the bistability observed in experimental works and already reproduced in more elaborate models that take into account the types of proteins involved in the adhesion and the role of hydrodynamic forces in the cell fate.

In Finger et al (1996), the authors investigate the specificity of short-lived L-selectin proteins in the adhesion process. More precisely, they put to light a shear-flow treshold above which L-selectin-based binding occurs. The authors emphasize that the first step in cell adhesion may require the recruitment of these short-lived proteins to sufficiently slow the cell down to promote the formation of longer-lived types of bonds. Let us note that our model can explain this observation. Indeed, assuming that the spontaneous formation of bonds is related to L-selectins, the observed treshold is described using a bond formation rate of the form $c(u,v)=c\mathds{1}_{[0,u_{*}]}(u)\mathds{1}_{[\underline{v},+\infty)}(v)$, for a given $\underline{v}>0$. Reproduction then accounts for the integrins dynamics. Moreover, our assumption of velocity-dependent dissociation rate is consistent with observations of Finger et al (1996), since L-selectin bonds form when the cell velocity is high and, consequently, so is the dissociation rate. Integrin-based bonds form for a smaller cell velocity, leading to longer lifetimes. Initially, one can assume that $n(0)=0$, $c\neq 0$, $r<<d$, and that $u>\mathrm{bar}{U}_{\alpha}$. Table 3b then shows that in the long-time limit, the density of bonds converges towards a stable sliding state. Then, if the cell velocity gets low enough, the cell dynamics is ruled by Table 3a with $r>d$, where the bistable behaviour between cell release and cell arrest is given. Our model is therefore able to succesfully explain the two main steps involved in the experimental observations.

More recent works have provided mechanical models taking into account the hydrodynamics forces involved in cell adhesion. In Reboux et al (2007), the cell is a rolling cylinder, and interestingly the model captures both sliding and rolling situations. The model shows a bistable behaviour between cell rolling (or weak resistance to sliding) and its release in the blood flow. In Efremov and Cao (2011); Li et al (2018), the cell is either a rectangle with a given zone for bond rupture, or a 2D circular cell in adhesive interaction with an elastic substrate. Both models predict bistability between the cell arrest and a near-release rolling state for intermediate shear rates. Furthermore, plots of the stationary velocity depending on the blood flow velocity show that, at stationary state, the cell adheres for a low shear rate and is in a rolling state above some threshold. The explanation relies on the competition between formation and rupture of bonds, that vary with the proteins involved in the adhesion.

The model we build in this paper is simpler, since it does not describe the cell geometry, so that the hydrodynamic description is minimal. As a consequence, cell rolling cannot be distinguished from sliding, and we are not able to provide a plot relating the blood and cell velocities. However, remarkably, the model still reproduces a bistable behaviour, suggesting that cell adhesion dynamics is crucial for this feature.

The stochastic limiting model provides an expression for the cell mean arrest time and shows the effect of fluctuations on the dynamics.

The stochastic limiting model writes as a diffusive Stochastic Differential Equation on the adhesion density, that carries a nontrivial noise term. It arises in a regime where the fluctuations in the local reinforcement of adhesions have nontrivial effects on the cell dynamics. When there is no effect of the cell velocity on its adhesion activity, the model writes as a linear CIR process, that already carries a dichotomic behaviour, and for which some first hitting time properties are known. We use a spectral method to perform numerical simulations of the corresponding probability density. Finally, in the full model describing the interaction between the cell velocity and its adhesion activity, we use a Fokker-Planck approach to derive an integral formulation of the mean arrest time of a cell. Numerical simulations of this quantity further confirm the dichotomy arising between the cell arrest and its unstopped displacement in the bloodstream, and show that the range of the arrest times is mostly sensitive to the stochastic fluctuations intensity, in agreement with experimental observations.

4.4 Perspectives

The modelling approach we developped gathers key components involved in the behaviour of cells circulating in blood vessels. It allows the study how the blood flow affects the cells ability to initiate and maintain an adhesive interaction with endothelial cells. These phenomena mediate either the stable cell adhesion to the wall, its rolling or sliding, or its release in the blood flow.

Further improvements would consist in extending the modelling framework in several directions. First, the model could take into account time-dependent rates, in order to consider a variable blood-flow velocity, illustrating the effect of the heart cycle on the adhesion dynamics. Moreover, it is natural to extend the model to a 2D setting, where the vessel wall is a surface and the cell geometry is taken into account, as well as the spatial repartition of adhesive proteins at the cell surface. This situation could be handled by adding a spatial structure to the population of bonds, following the framework of Fournier and Méléard (2004); Champagnat and Méléard (2007) and further works on measure-valued stochastic processes. This framework would allow the precise description of hydrodynamic forces exerted on the cell, so that both rolling and sliding phenomena could occur.

Finally, further works should focus on the comparison with experimental measures to to shed light on cell arrest conditions. In particular, in the work of Follain et al (2017), the authors study Circulating Tumor Cells (CTCs) in vivo, and keep track of the cell arrest sites with respect to a tuned hemodynamic flow velocity. They show that cell arrest in blood vessels occurs at sites with permissive flow profiles. Since CTCs are large, fairly rigid cells, our point-particle cells assumption is particularly relevant in this context. In addition, the density of adhesion molecules on the vessel wall is experimentally controlled, so that each feature of the model is characterized by experimental data. In this perspective, confronting the model to these data is of great interest, and could be supplemented with the theoretical study of the distribution of arrest positions.

Acknowledgements.

The authors are very grateful to V.C. Tran and R. Voituriez for very helpful discussions and suggestions. They also wish to thank the reviewers for their relevant suggestions, and F. Gidel for proofreading the manuscript.

Appendix A Mathematical properties of the rescaled process

In this part, we study the rescaled process $(Z_{t}^{K})_{t}$ for a constant $K$, in a general framework that includes the situations studied in this paper. More precisely, we consider rates that satisfy the following hypothesis:

Hypothesis 1

For all $z\in\mathbb{R}_{+}$,

[TABLE]

where $a>0$. In addition the disassembly rate $z\mapsto d^{K}(z)$ is continuous. The total formation and dissociation rates write

[TABLE]

By construction, $(Z_{t}^{K})_{t\geq 0}$ is also a Markov process, and for $\Phi:\mathbb{R}_{+}\rightarrow\mathbb{R}$ measurable bounded, its infinitesimal generator writes

[TABLE]

We show the following Proposition.

Proposition 5 (Moment and martingale properties)

Under Hypothesis 1, if there exists $p\geq 2$ such that $\mathbb{E}\left[(Z_{0}^{K})^{p}\right]<+\infty$, then

$\forall\;T>0,$**

[TABLE] 2. 2.

for all measurable $\Phi:\mathbb{R}_{+}\rightarrow\mathbb{R}$ for which there exists $C$ such that $\forall z\in\mathbb{R}_{+},\;|\Phi(z)|+|L^{K}\Phi(z)|\leq C(1+z^{p})$,

[TABLE]

is a càdlàg $(\mathcal{F}_{t})_{t\geq 0}$-martingale starting from [math]. 3. 3.

The process defined by

[TABLE]

is a càdlàg square-integrable martingale starting from [math] and of quadratic variation

[TABLE]

Proof

Although the proof is similar to Proposition 2.7 of (Bansaye and Méléard, 2015) we recall it here for clarity. In the following $C$ will denote a positive constant which value can change from line to line. Let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space, $Z_{0}$ an integer-valued random variable, and $M(\,{\mathrm{d}}s,\,{\mathrm{d}}w)$ an independent Poisson Point Measure on $\mathbb{R}^{2}_{+}$, of intensity measure $\,{\mathrm{d}}s\,{\mathrm{d}}w$. Denote by $(\mathcal{F}_{t})_{t\geq 0}$ the canonical filtration generated by these objects. Then, we define the $(\mathcal{F}_{t})_{t\geq 0}$-adapted càdlàg process $(Z^{K}_{t})_{t\geq 0}$ as the solution of the following SDE: $\forall t\geq 0$,

[TABLE]

This representation is classical (see e.g Fournier and Méléard (2004); Champagnat and Méléard (2007)). The Poisson jumps related to the measure are accepted or rejected thanks to the indicator functions. The variable $w$ is then used as an acceptance parameter in order to obtain the desired rates for each event. Now, writing that $\mathbb{P}-a.s$, for a positive and measurable test function $\Phi$,

[TABLE]

so that for $\Phi(Z^{K}_{t})=\left(Z^{K}_{t}\right)^{p}$ and neglecting the negative death term, we obtain

[TABLE]

Taking expectation, and using that for $p\in\mathbb{N}$, $(1+z)^{p}-z^{p}\leq C(p)(1+z^{p-1})$, we obtain that

[TABLE]

We conclude with the Gronwall Lemma together with the assumption on the initial condition. The second item of Proposition 5 is a classical property of Markov processes, and (21) is obtained as a consequence. Finally, to obtain the quadratic variation (22), one has to write the semi-martingale decomposition of $\left(Z_{t}^{K}\right)^{2}$, one the one hand by applying (20) to $\Phi(z)=z^{2}$, and on the other hand by applying (20) to $\Phi(z)=z$ before using the Itō formula. The conclusion follows from the uniqueness of the decomposition.

Appendix B Proof of Theorem 3.1

The proof is similar to the ones of (Joffe and Metivier, 1986; Ethier and Kurtz, 2009), and is based on a compactness-uniqueness argument. First, note that since the rates under study have $K$-independent upper bounds, reproducing the proof detailed in Appendix A lead to

[TABLE]

The uniform tightness of the sequence of laws $(Q^{K})_{K}$ of $(Z^{K})_{K}$ follows. Then, from the Prokhorov theorem we deduce the relative compactness of the family of laws $(Q^{K})_{K}$ on $\mathbb{D}([0,T],\mathbb{R}_{+})$. Consider now a convergent subsequence of limit $Q$, and a corresponding sequence of processes converging in distribution to some $n\in\mathbb{D}([0,T],\mathbb{R}_{+})$ of law $Q$. Our aim is thus to identify this limit. Firstly, since the jumps of $(Z_{t}^{K})_{t}$ are of the form $1/K$, we know that any process of law $Q$ is almost surely strongly continuous. Now, for $t\leq T\wedge\tau_{u}$, denote

[TABLE]

We can easily prove that for all $t\leq T$, $\mathbb{E}_{Q}\left[\left|\Psi_{t}(n)\right|\right]=0$. Finally, the convergence follows from the uniqueness of solutions to (6) in $\mathcal{C}([0,T],\mathbb{R}_{+})$, which comes from the Lipschitz-continuity of the disassembly rate.

Appendix C Proof of Theorem 3.2

C.1 Uniform estimates

We first prove two propositions that provide uniform estimates on the process. They will be used to show the tightness of any sequence of laws associated with $(Z_{t}^{K})_{K}$.

Proposition 6

Under the assumptions of Theorem 3.2, if

[TABLE]

then for $T<+\infty$,

[TABLE]

Proof

First, the infinitesimal generator (19) associated with $\Phi(x)=x^{2}$ writes

[TABLE]

Therefore, we obtain that

[TABLE]

and since the Gronwall Lemma yields $\mathbb{E}[Z^{K}_{s}]\leq C(1+\mathbb{E}[(Z_{t}^{K})^{2}])$, we deduce that there exists $C(T)$ such that $\mathbb{E}[(X_{t}^{K})^{2}]\leq C(T)$, and the Proposition is proved.

Proposition 7

Under the assumptions of Theorem 3.2, if $\sup_{K>0}\mathbb{E}[(Z_{0}^{K})^{2}]<+\infty$, then for $T<+\infty$,

[TABLE]

Proof

We first deduce from (21) that

[TABLE]

We want to take the expectation in this inequality. For that purpose, we first use the Burkholder-Davis-Gundy inequality to write that

[TABLE]

Now, since $\mathbb{E}[Z^{K}_{0}]<+\infty$, we obtain from (23) that

[TABLE]

We then use (22) to get

[TABLE]

thanks to Proposition 6. We conclude using the Gronwall lemma.

C.2 Proof of convergence

The proof follows the same outline as in the deterministic case. First, we prove similarly that the sequence of laws $(Q^{K})_{K}$ of the processes $(Z^{K})_{K}$ is uniformly tight in $\mathcal{L}(\mathbb{D}([0,T],\mathbb{R}_{+}))$. Indeed, denote $(A_{t}^{K})_{t\geq 0}$ the finite variation process associated to $(Z_{t}^{K})_{t}$. The Aldous and Rebolledo criterion (Aldous, 1989) states that we need to prove that for all $T>0$ the following inequalities hold true:

1. a)

$\displaystyle\sup_{K>0}\mathbb{E}\left[\sup_{t\in[0,T]}\left|Z_{t}^{K}\right|\right]<+\infty$. 2. b)

$\forall\varepsilon>0,\;\forall\eta>0,\;\exists\delta>0,\;K_{0}\in\mathbb{N}^{*}$ such that for all sequence $(\sigma_{K},\tau_{K})_{K}$ of stopping times with $\sigma_{K}\leq\tau_{K}\leq T,$

1. (i)

[TABLE] 2. (ii)

[TABLE]

These estimates can easily by obtained by direct computations, using Proposition 7 and the Markov inequality. Now, we aim at identifying the limiting values. For $Y=(Y_{t})_{t\geq 0}\in\mathbb{D}([0,T],\mathbb{R}_{+})$, let us define

[TABLE]

We need to show that $\widetilde{M}_{t}(N)$ is a twice-integrable continuous martingale with quadratic variation process defined by

[TABLE]

First, we show that $\widetilde{M}(N)$ is a martingale. Take $0\leq s_{1}<...<s_{n}<s<t$, and $\Phi_{1},...,\Phi_{n}$ continuous bounded functions from $\mathbb{R}$ to $\mathbb{R}$. Define $\Psi$ on $\mathbb{D}([0,T],\mathbb{R})$ by

[TABLE]

As for the proof of Theorem 3.1, we can show that $\mathbb{E}[\Psi(N)]=0$. The main novelty in this proof consists in showing that the bracket of $\widetilde{M}$ is given by (25). We proceed in two steps.

First, we consider $K$-dependent semimartingale obtained from (20) with $\Phi(Z^{K})\leavevmode\nobreak\ =\leavevmode\nobreak\ (Z^{K})^{2}$, that is related to the infinitesimal generator given by

[TABLE]

We can show that, at the limit, we obtain the following martingale:

[TABLE] 2. 2.

Then, the Itō formula applied to (24) show that

[TABLE]

is a martingale. We conclude by the uniqueness of the semimartingale decomposition of $N_{t}^{2}$.

The equivalence between the martingale problem (24)-(25) and the SDE (8) follows from the classical martingale identification

[TABLE]

see for example Karatzas and Shreve (1998). Finally, the pathwise uniqueness of the solution to (8) is classical in dimension 1 since the drift term is Lipschitz-continuous and the diffusion coefficient is $1/2$-Hölder (see e.g Ikeda and Watanabe (1989)).

Remark 5

The solution is strong and has the strong Markov property.

Appendix D The CIR process

In this part, we give some classical results about the CIR process. For more details we refer to Göing-Jaeschke and Yor (1999); Shreve (2004).

Reaching zero

First, it is known that for$c>0$ and a positive initial state, $\{N=0\}$ is a reflective barrier for the process (9). The properties of the probability $\mathbb{P}^{0}$ to hit zero are gathered in the following table.

[TABLE]

These results can be intuitively understood when $r-d\leq 0$, using the correspondance between CIR processes and Orstein-Uhlenbeck processes. Indeed, consider $D$ such processes $(X^{1},\,\cdots,\,X^{D})$ such that $\forall 1\leq i\leq D$,

[TABLE]

with $(B^{i})_{i}$ independent Brownian motions and $\beta>0$. Each process follows a stochastic dynamics that is drifted to zero. The Itō formula allows the derivation of the SDE satisfied by $R=(X^{1})^{2}+\cdots+(X^{D})^{2}$, the squared euclidean norm of $(X^{1},\,\cdots,\,X^{D})$. We obtain

[TABLE]

with $B$ a Brownian motion. The CIR Equation (9) is therefore obtained for $r-d=-\beta\leq 0$, $\sigma^{2}=2a$ and $D=4c/\sigma^{2}>0$. As a consequence, when $4c/\sigma^{2}$ is an integer, the CIR process writes as the squared norm of $D$ Ornstein-Uhlenbeck processes. Therefore, using the properties of the Brownian motin, the CIR process almost surely hits zero infinitely many times when $D=1$, while it has a null probability of reaching zero when $D\geq 2$. As a remark, note also that the CIR process can be rigorously related to the Squared Radial Ornstein-Uhlenbeck process.

Distribution

When $\delta:=\frac{2c}{a}\in\mathbb{N}$ and $r-d<0$, the CIR process writes as the squared norm of a $\delta$-dimensional Ornstein-Uhlenbeck process. Consequently, we have that

[TABLE]

where $Y_{t}$ follows a non-central Chi-square distribution with $\delta$ degrees of freedom, and a non-centrality parameter $\xi_{t}=n_{0}\frac{2(d-r)}{(1-e^{(r-d)t})a}e^{(r-d)t}$. Moreover, the mean solution and its variance can be computed from (9), leading to

[TABLE]

For $n\geq 0$, $n_{0}>0$ and $\kappa_{t}=\frac{d-r}{(1-e^{(r-d)t})a}$, the probability density $p_{n_{0}}$ of the CIR process writes

[TABLE]

where

[TABLE]

is the modified Bessel function of the first kind, with $\Gamma$ the Gamma function defined by $\Gamma(t)=\int_{0}^{\infty}x^{t-1}e^{-x}\,{\mathrm{d}}x$. A stationary distribution exists if and only if $r-d<0$ (Ikeda and Watanabe, 1989). In this case, following a Fokker-Planck approach, one can check that the stationary density writes

[TABLE]

with $\phi(n)=\left(1-\frac{c}{a}\right)\ln(n)+\frac{d-r}{a}n$ the corresponding potential, and $\mathcal{N}$ a normalization constant.