A continuous-time stochastic model of cell motion in the presence of a chemoattractant

We consider a force-based model for cell motion which models cell forces using Hooke's law and a random outreach from the cell center. In previous work this model was simplified to track the centroid by setting the relaxation time to zero, and a formula for the expected velocity of the centroid was derived. Here we extend that formula to allow for chemotaxis of the cell by allowing the outreach distribution to depend on the spatial location of the centroid.


1.
Introduction. Cell motion is key to most biological systems ranging from bacteria to single cell amoeba to vertebrates; in this work we focus on amoeboidal cell motion. In [4], a differential equation model for cell motion was introduced and in [5] a simplified version of the model was considered where the relaxation time was set to zero and the centroid was tracked. In both models, a "cell center" is coupled with adhesion sites. In the case where the interevent times for the adhesion sites are exponentially distributed, the continuous-time centroid model (CTCM) was shown to be a pure jump-type continuous-time Markov process and a formula for the velocity of the expected value of the centroid was derived. Here we extend those results to allow for chemotaxis of cells by generalizing the formula for the velocity of the expected position.
Cell motion and chemotaxis has been modelled for many years with differential equations based on random walks [13,10]. From those early works, a rich body of work has been developed [8,2,16,15]. Our work differs from this body of literature in that it is not a partial differential equation model of an aggregate of cells. It is a model for an individual cell with stochastic attachments of the adhesion sites more closely related to [3,1].
We begin by reviewing the differential equation model introduced in [4] and the CTCM introduced in [5].
2. Review of models. In [4], the cell was modeled as a nucleus and multiple integrin-based adhesion sites (I-sites) [6,7,14] which exert forces on the nucleus. I-sites attach to an external substrate and, once attached, remain fixed to that substrate location for a duration determined by a given probability distribution. Likewise the duration which an I-site remains unattached is determined by a (potentially different) given probability distribution.
Let n denote the number of I-sites, u i denote the location of the ith I-site in a Euclidean space E (usually R 2 or R 3 ), and let the random variable ψ i ptq denote whether it is attached. For a fixed realization the equation describing x, the location of the cell at time t, is where Here the ζ i are spring constants associated with the forces exerted by the I-sites and a p,i is the last time that ψ i changed from 0 to 1, and b p,i are independent, identically distributed, random vectors with distribution η. In terms of the attachment and detachment times, a p`1,i " κ p,i`dp,i`ap,i where κ p,i is the random variable drawn from an exponential distribution indicating how long I-site i remains attached after time a p,i , in other words κ p,i is the time of the attachment event. Likewise d p,i is the time the I-site remains detached after the detachment event occurs at time a p,i`κp,i . In [5] the above model was modified to the continuous time centroid model. This modification was motivated by informally considering the limit of the above model as the forces become very large (that is the ζ i 's become large). In this limit, the smooth motion of the cell nucleus disappears and instead the motion of the nucleus (which is located at the centroid of the I-sites) jumps from position to position. The motion of the centroid is a right-continuous jump process. We note that the location of a given I-site only changes at the moment it attaches and it's new location is a random perturbation from the centroid the moment before attachment occurs. Whenever the only attached I-site detaches (leaving no I-sites attached), the "centroid" is considered to be the same as it was immediately before this detachment occurred.
We now modify Equation (1) to reflect these changes and, for clarity to distinguish between the two models, we will denote the location of the I-site's by v i . We restrict the number of I-sites to n, and keeping with the notation of [5], allow i to range from 0 to n´1. Again fixing a realization, the equation describing c, the location of the centroid of the cell at time t, is 0 " Here as before a p,i is the last time that ψ i changed from 0 to 1, and b p,i are independent, identically distributed, random vectors with distribution η.
In this manuscript we add chemotaxis to the model. In order to do this, we let the outreach distribution η be a family of distributions parameterized by location denoted η y . Thus, the cell can place a new attachment in a preferential direction based on a chemotactic signal. This necessitates the addition of a regularity assumption to prove Proposition 3.1*. Specifically, assume that η y is weakly continuous with respect to y. (That is, if y k Ñ y 0 and f is a continuous, bounded, and real-valued function, then ş f pxq dη y k pxq Ñ ş f pxq dη y0 pxq.) Additionally, we modify the previous boundedness assumption, which is necessary for Lemma 4.5*. In particular, we assume that η y is supported on a |¨|-ball of radius R centered at the origin for every y where |¨| is the 8-norm on E.
3. Mathematical analysis of the chemotactic model. As was done for the simpler model in [5], we wish to show that the chemotactic model informally described above can be formulated in a mathematically precise way, that this precise mathematical formulation gives rise to a jump-type, continuous-time, Markov process (denoted X t ), and that the mean cell velocity corresponding to this process is the solution of a tractable integral equation. A key part of this task is the construction of a rate kernel.
At a fixed time, each I-site has a location in a Euclidean space E and a status in t0, 1u, with 0 meaning "detached" and 1 meaning "attached". Because the centroid of the attached I-sites is important, we wish to explicitly include it in the state variable. Our state space X can thus be thought of as the disjoint union of a finite number of subspaces of the Euclidean space of dimension pn`1q dimpEq. We denote our state space as Here and below, we adopt the standard convention that if A and B are sets, then A B is the collection of functions from B to A. Furthermore, we introduce the abbreviation n˘for n˘1 and (for convenience) use rns to represent the set of nonnegative integers less than n. The vector vpiq should be interpreted as the location of I-site i if 0 ď i ă n, whereas vpnq is the location of the centroid of the attached I-sites. The quantity ψpiq should be interpreted as the status of I-site i.
In accordance with the usage in [9], the rate kernel α will be the product of a rate function c and a transition probability kernel µ. While, cpxq represents the reciprocal of the expected wait time to the next jump when in state x, and µpx, Bq represents the probability that the next jump lands in set B if it starts in state x.
The rate function is fairly easy to define. It is given by the formula where |ψ| :" ř iPrns ψpiq (the number of attached I-sites), and θ a and θ d are, respectively, the expected time for a detached site to attach and for an attached site  Figure 1. Panel a) and c) depicts the state space x e given as an example. The circle represents the cell centroid. The black "x"s represent attached I-sites (labeled 0 and 1) and the red "x" denotes the detached I-site (labeled 2). The sets A and B represent the example sets used in the measure for case 1 and case 2. Panels a) and b) are for case 1 where I-site 1 dettaches. Panel a) is the initial configuration and panel b) is the new state. Likewise, panels c), and d) are for cases 2 where I-site 2 attaches. The setsĀ and B are the sets A and B transformed by S´1˝F´1. In panel a), A " A andB " B`p.5, .5q. In panel c),Ā " A´p.5, .5q and B " 3B´p1.5, 1.5q and their intersection contains η " p1, 0q. In panel d), I-site 2 has attached at location x 2 " c`η " p.5, .5qp 1, 0q.
to detach. The validity of this formula depends on the behavior of independent exponential random variables.
On the other hand, the correct formula for µ is complicated, so, of necessity, we build it up in stages. Two families of auxiliary functions that will facilitate the technical aspects of this process are F i px, yq :" tpi, xq, pn, yqu and S pa,b,a,bq px, yq :" papx´aq, bpy´bqq.
Observe F i indexes tuples and S pa,b,a,bq is a scale-and-translate function.
Example cont.
Using the example state x e , if I-site 2 was about to change status it is useful to consider F 2`p 0, 1q, p.5, .5q˘" The measures µ pψ,vq ti,nu are the lowest-level objects that deal specifically with our evolutionary model, and also have the most complicated formulas. There are three cases to consider. The first is when the only attached I-site detaches. In this scenario, the centroid and I-sites do not change location, only the attachment state of the I-site changes. The second case is when there are several I-sites attached and one detaches. The only changes are in the centroid location and the state of the I-site that detached. The final case is when an I-site attaches. Both the I-site and the centroid have new locations. The resulting formula for µ pψ,vq ti,nu is given by where I : EˆBpEq Ñ r0, 1s is the inclusion kernel defined by the formula Ipx, Cq " δ x pCq, so Given that the current state is pψ, vq, the measure µ pψ,vq ti,nu represents the conditional probability that the I-site i will attach in a given subset of E at the next jump, conditioned on the event that it is I-site i which will attach or detach at that time.
Thus if p1, 1q P A and p0, 0q P B the probability that if I-site 1 detaches the new state will have I-site 1 located in A and the centroid located in B is 1.
This means that µ xe t2,3u`AˆB˘" pη vp3qˆI q´`A´p.5, .5q˘ˆ`3B´p1.5, 1.5q˘" η p.5,.5q´`A´p .5, .5q˘X`3B´p1.5, 1.5q˘Ī f x 2 is the new location of I-site 2, then the new location of the centroid is c n " 1 3`p 1, 1q`x 2˘. If x 2 P A and c n P B the probability that if I-site 2 attaches to location x 2 then new state will have I-site 2 located in A and the centroid located in B is the measure ofĀ XB (shown in Figure 1) as determined by η p.5,.5q . In the figure η p.5,.5q " δ p1,0q and thus the measure is 1.
Corresponding measures that account for the other ("inactive") I-sites are µ pψ,vq tiu where G i pxq :" tpi, xqu. If we let P i :" ttju : j P rnsztiuu Y tti, nuu (so that tP i : i P rnsu is the partition of rn`s consisting of ti, nu and n´1 singletons), the product ą pPPi µ pψ,vq p measures the probability of the I-sites' and centroid's locations being in a given subset of E rn`s at the next jump, conditioned on I-site i being the one that changes status.
Example cont.
Since I-site 0 is "inactive" in both cases of our example µ pψ,vq tiu is the appropriate measure to use when considering I-site 0. Let C and D be subsets of R 2 and D " CˆAˆBˆD, then If p0, 0q P C then the measure of C is 1. For case 1, when I-site 1 detaches, the partition P 1 " t0u, t1, 3u, t2u ( is needed for the product measure This conditional probability, which deals with physical locations, needs to be combined with a corresponding conditional probability that deals with attachment/detachment statuses. The latter is simply given by δ sipψq , where s i pψq :" pψztpi, ψpiqquq Y tpi, 1´ψpiqqu, so s i pψq is the necessary combined status of ψ that results from I-site i changing status. Given a system in state pψ, vq, δ sipψqˆą pPPi µ pψ,vq p is the conditional probability of a jump landing in a given subset of X, conditioned on I-site i being the one to change status at the time of the jump.
These conditional probabilities need to be multiplied by the corresponding probabilities that I-site i changes status at the next jump. The hypothesis that wait times are exponential leads to the formula for the latter probability.
Example cont.
For case 1, In light of the Law of Total Probability, the formula for µ we finally arrive at is (the restriction to X of) Many results concerning this chemotactic system follow readily from results established in [5] for the special case in which η x is independent of x. To make it easier to keep track of the connection between that paper and this one, we number results here the same way they were numbered there. Ones that change due to the new assumption but have no appreciable changes in the proofs are denoted with an asterisk. Theorem 4.11 which involves substantial changes is denoted with a :.
Proposition 3.1* The map µ is a probability kernel from X to X.
Proof. Replace η with η vpnq in the proof of Proposition 3.1 in [5]. Because η x is weakly continuous with respect to x, µpx, Aq is measurable with respect to x.
Proposition 3.2* For any Borel probability measure ρ on X, there is a discretetime Markov process Y on X with transition kernel µ such that Y 0 (i.e., the value of Y at time 0) is ρ-distributed.
Proposition 3.3* For every Borel probability measure ρ on X, there is a pure jump-type continuous-time Markov process X on X with rate kernel α such that X 0 is ρ-distributed. If Y is as defined in Proposition 3.2*, pγ i q is a sequence of standard exponential random variables, and tY, γ 1 , γ 2 , γ 3 , . . .u is an independent family, then X can be defined by the formula X t " Y k for t P rτ k , τ k`q , where τ k :" The map πpψ, vq :" |ψ| connects X with the much simpler state space rn`s, and there is an evolution law on the latter that corresponds with that derived above on X. We claim that the appropriate rate function for this reduced system iŝ cpiq :" θ d i`θ a pn´iq, with the transition probability kernel beinĝ µpi,¨q :" so the rate kernel isαpi,¨q :"ĉpiqμpi,¨q " θ d iδ i´`θa pn´iqδì .
Proposition 4.2* If X is as in Proposition 3.3*, thenX :" π˝X is a pure jump-type continuous-time Markov process with rate kernelα and initial distribution ρ˝π´1.
Proof. In [5] we derive a version of Dynkin's criterion that demonstrates a natural correspondence between two discrete Markov processes. This criterion is used in the proof of Proposition 4.2 in [5] and similar application gives the proof of this proposition. IfẐ is a pure jump-type continuous-time Markov process with rate kernelα, then the distribution ofẐ t converges to σ as t Ñ 8, regardless of the distribution ofẐ 0 .
Proof. See the proof of Proposition 4.3 in [5].
Lemma 4.4* Let f be a nonnegative measurable function on X, and let x " pψ, vq P X.
Then ż and Proof. Replace η with η vpnq in the proof of Lemma 4.4 in [5].
Now, let |¨| be the 8-norm on E, let f i pψ, vq :" vpiq, and let gpxq :" maxt|f i pxq| : i P rn`su. Lemma 4.5* If Y is as in Proposition 3.2*, and k is a whole number, then gpY k q ď gpY 0 q`kR almost surely.
Proof. Modify the proof of Lemma 4.5 in [5] as follows: Lemma 4.6* Let Y be as in Proposition 3.2*, and let k 1 and k 2 be whole numbers. Then |f n pY k2 q´f n pY k1 | ď 2gpY 0 q`pk 1`k2 qR almost surely.
Proof. See the proof of Lemma 4.6 in [5], with Lemma 4.5* being applied rather than Lemma 4.5.
Proposition 4.8* Let ρ be a distribution on X such that f i is ρ-integrable for every i. Let X be as in Proposition 3.3*. Then for every i P rn`s and t ě 0, the expected value of f i pX t q is well-defined and finite.
Proof. See the proof of Proposition 4.8 in [5], with references to Proposition 3.2, Proposition 3.3, and Lemma 4.5 replaced by references to their starred counterparts.
Proposition 4.10* Let ρ be a distribution on X such that, for every i, f i is ρintegrable. Let X be as in Proposition 3.3*. Then t Þ Ñ Epf n pX t qq is continuous.
Proof. See the proof of Proposition 4.10 in [5], with references to Proposition 3.3, Lemma 4.6, and Proposition 4.8 replaced by references to their starred counterparts.
Theorem 4.11 † Let σ be as in Proposition 4.3* and let ρ be a distribution on X such that σ " ρ˝π´1 and such that f i is ρ-integrable for every i. Let X be as in Proposition 3.3*. Then

B Bt
Epf n pX t qq " for every t ą 0, whereηpxq :" ş R N x dη fnpxq pxq and µpX t p¨qq is the distribution of the random process X at time t.
Proof. Following the proof of Theorem 4.11 in [5], by Proposition 4.10* and the mean value theorem for one-sided derivatives in [12], it suffices to verify that (20*) holds for the right-hand derivative. The time-homogeneity of X means that it suffices to set t " 0. (Proposition 4.8* implies that the hypothesized integrability condition translates to a corresponding integrability condition after a time shift, and Propositions 4.2* and 4.3* imply that the projected distribution σ does not change.) Let Y , pγ k q, and pτ k q be as in Proposition 3.3*, δ ą 0, and hpx, yq " f n pyq´f n pxq. For each pi, jq P rn`sˆrn`s, define Ω i,j :" tπpY 0 q " i, πpY 1 q " ju. Then (It should be noted that for t ‰ 0 the statement would be that pY , Y `q is pµ ˆµq-distributed where µ is the distribution of the discrete-time Markov process Y described in Proposition 3.2* at iteration . Of course the important property for our result is that the projected distribution is the stationary distribution σ, which both ρ and µ for any satisfy.) So PpΩ i,j q " pρˆµqpπ´1ptiuqˆπ´1ptjuqq " " |ψ| θ d θ d |ψ|`θ a pn´|ψ|q 1 tj`u p|ψ|q`pn´|ψ|q θ a θ d |ψ|`θ a pn´|ψ|q 1 tj´u p|ψ|q " θ d pj`q θ d pj`q`θ a pn´j`q 1 π´1ptj`uq pxq`θ a pn´j´q θ d pj´q`θ a pn´j´q 1 π´1ptj´uq pxq.
Plugging this into (22*) and using the fact that ρpπ´1p¨qq " σ and that for t ‰ 0, µ pπ´1p¨qq " σ, we get For the other factor of the summand in (21*), we have From the formula for τ k , the independence of Y and tγ k u, and the fact that pY 0 , Y 1 q is pρˆµq-distributed, we get Epf n pY 1 q´f n pY 0 q | Ω i,j X tt P rτ 1 , τ 2 quq " Eˆf n pY 1 q´f n pY 0 qˇˇˇˇΩ i,j X " γ 1 cpY 0 q ď t ă γ 1 cpY 0 q`γ 2 cpY 1 q

*"
Eˆf n pY 1 q´f n pY 0 qˇˇˇˇΩ i,j X " γ 1 cpiq ď t ă γ 1 cpiq`γ 2 cpjq *" Epf n pY 1 q´f n pY 0 q | Ω i,j q " 1 PpΩ i,j q ż π´1ptiuqˆπ´1ptjuq pf n pyq´f n pxqq dpρˆµqpx, yq if PpΩ i,j q ą 0. The proof proceeds as in the proof of Theorem 4.11 up to (27) since anything involving Lemma 4.4* or µ either uses a fixed x P X, uses the uniform bound R on the support of η y , or gives something that does not involve η y . Equation (27)  If t ‰ 0 the equation is
As a simple example, consider the case where the source for a chemoattractant is located at a fixed location p. If the outreach is deterministic and a unit vector in the direction of the source, it is given by p´f n pxq |p´f n pxq| .
If all the I-sites start at location y, the initial distribution of X 0 is non-zero if and only if the I-sites and the centroid are all at y and compatible with the steadystate distribution of the number of attached I-sites. In this caseηpxq is just ppf n pxqq{|p´f n pxq|, so the formula becomes B Bt Epf n pX 0 qq " ÿ iPrns θ a pn´iq i`ż π´1ptiuq p´f n pxq |p´f n pxq| P X0 pdxq " ÿ iPrns θ a pn´iq i`ż π´1ptiuq p´f n pxq |p´f n pxq| dρpxq p´y |p´y| 1 pθ d`θa q n θ d ppθ d`θa q n´θ d q If the random process X t starts in this configuration, the centroid will always lie on the line defined by y and p. Thus p´f n pxq " p´αpp´yq, so p´f n pxq |p´f n pxq| " p´αpp´yq |p´αpp´yq| .
This is a unit vector on the line so it is˘pp´yq{|p´y|. Thus for any time t

B Bt
Epf n pX t qq "˘p´y |p´y| where the sign is determined by which side of p the centroid, f n pX t q, is located. The cell has a constant speed but the direction changes sign since the cell will usually pass by the source due to the unit vector outreach.
4. Discussion. We have extended the main result previously reported in [5] to include a space dependent distribution for the location of the attachment sites. Theorem 4.11* gives a formula for the velocity of the expected position of the centroid when the attachment sites are governed by a space dependent distribution. This increases the value of the model introduced in [5] by extending its uses to situations where the cell motion is governed by external cues such as chemotaxis, haptotaxis, durotaxis, or other space-based taxis. The formula without space-based taxis (Equation (20) in [5]) is more amenable to analytical solutions. When spacebased taxis is added, the formula given in Equation (20*) is in terms of the distribution of X at time t. This distribution is not easily determined from the initial distribution of the system and thus numerical methods will typically be required to determine the velocity.