Analytical results on the Beauchemin model of lymphocyte migration

Abstract

The Beauchemin model is a simple particle-based description of stochastic lymphocyte migration in tissue, which has been successfully applied to studying immunological questions. In addition to being easy to implement, the model is also to a large extent mathematically tractable. This article provides a comprehensive overview of both existing and new analytical results on the Beauchemin model within a common mathematical framework. Specifically, we derive the motility coefficient, the mean square displacement, and the confinement ratio, and discuss four different methods for simulating biased migration of pre-defined speed. The results provide new insight into published studies and a reference point for future research based on this simple and popular lymphocyte migration model.

Background

A unique property of the immune system is that it mainly consists of constantly moving cells. Technological progress has facilitated the study of lymphocyte migration at increasingly higher resolutions, culminating in the application of two-photon microscopy to tracking single lymphocytes in the living animal [1–3]. Accompanying this progress, mathematical and computational models of lymphocyte migration have been developed and applied to both qualitative and quantitative questions. These models can be roughly categorized as follows according to their level of detail: Whole-population models, usually formulated using ordinary or partial differential equations, do not distinguish between individual cells, but treat cell subpopulations as continuous quantities. Particle-based models do consider each cell individually, but treat cells as freely moving particles without mass, volume, and shape. Finally, whole-cell models such as the Cellular Potts Model [4] explicitly represent the cell and its interaction with the environment. In the present paper, we analyze a particle-based model that was introduced by Beauchemin, Dixit, and Perelson [5]. Throughout the paper, we refer to this model as the "Beauchemin Model".

Particle-based approaches such as the Beauchemin model are used for studying questions where a whole-population approach does not provide sufficient information, but a whole-cell model is not necessary, not feasible for computational reasons or would require too many assumptions on unknown parameters. For instance, Grigorova et al. [6] and ourselves [7] recently used the Beauchemin model to study the transit of T cells through a lymph node in the absence of antigen, and we used it to determine the amount of directional bias that could be detected in random T cell migration using contemporary two-photon imaging techniques [7].

The Beauchemin model is defined by three parameters - a speed v_free, a time interval t_free, and another time interval t_pause. The model describes a particle moving freely in a three-dimensional space according to the following rules: The particle starts at a fixed position and turns in a random direction (more precisely, a direction is sampled from a uniform distribution on the unit sphere). The turn takes a fixed time t_pause, during which the cell does not move. This reflects the time it takes to "displace the internal structure which brings about motion" [4]. Then, the cell moves in the chosen direction with velocity v_free for a fixed time t_free, after which it stops and restarts the process.

The main distinguishing characteristic of the model is the introduction of the pause phase. In this sense, it can be seen as a generalization of the "ideal chain" model of polymer physics [8], in which a chain of rigid rods of fixed length ℓ is considered that are freely jointed to each other. This is mathematically equivalent to the Beauchemin model if we set v_free · t free = ℓ and t_pause = 0. In particular, analytical results e.g. on the mean square displacement of the ideal chain [8] are directly transferable to this case. As we shall see later on in this paper, for some quantitative properties even the more general case with t_pause ≥ 0 can be reduced to the ideal chain model.

This paper is structured as follows: In the upcoming section, we begin by formulating a mathematical framework for our analysis, which explicates sufficiently but not overly general notions of random walks and Brownian motion. We deliberately include a fairly large amount of detail, including proofs, in order to present our results in a self-contained fashion. Thereafter, we present analytical results about random motion in the Beauchmin model, and discuss four different ways for simulating biased migration with pre-defined speed. We wrap up with concluding remarks.

We point out that around half of the results discussed in this paper originate from our recent study that used the Beauchemin model [7]. We chose to integrate these results into the present paper in rewritten and extended form because they were formerly only published in condensed form as supporting online material. In this manner, we hope to provide a coherent account of all available analytical results, which we anticipate to serve as a useful reference for future users of the model. We mark results that have appeared earlier in our previous study or elsewhere (e.g. Proposion 10) by giving the appropriate reference in proposition or theorem statements, while propositions or theorems stated without such a reference (e.g. Proposion 12) have, to our knowledge, not been published before.

Mathematical foundations

In this section, we outline the basic mathematical framework that we use to describe stochastic cell motion and, later on, the Beauchemin model. A similar, but slightly less general framework (based on random walk on a lattice) is used in Berg's textbook Random Walks in Biology [9], while a much more general account of random walks than given here is usually presented in probability textbooks (e.g., [10]). At the end of this section, we will outline the connection of this random walk to the corresponding macroscopic partial differential equation model, which is the convection-diffusion-equation.

First of all, we introduce some notational conventions. We denote the expectation of a random variable ξ by E[ξ], its variance by Var(ξ) = E [(ξ − E[ξ])²] and its probability density function by f_ξ. Elementary random variables are written in lowercase Greek letters (usually ξ _i ) while derived random variables that are functions of elementary random variables are written in capital Latin letters. All random variables are assumed to be continuous and real-valued. Mean and variance of random variables will also be written as μ and σ², if the referenced variable is clear from context.

Observed tracks and associated measures

The Beauchemin model is usually validated against, or used to generate, cell tracking data as it would be recorded in a two-photon microscopy experiment. As we will see later, there can be subtle differences between the real underlying state and the observed state of a cell. We thus make this difference explicit in our analysis by means of the following definitions for cell tracking terminology.

Definition 1 (Cell track). A d-dimensional cell track is a finite sequence

$$\left(x\left(1\right),t\left(1\right)\right),\dots ,\left(x\left(k\right),t\left(k\right)\right)$$

consisting of positions $x\left(1\right),\dots ,x\left(k\right)\in {ℝ}^d$and increasing time points $t\left(1\right),\dots ,t\left(k\right)\in ℝ$.

The time indices are written in brackets to indicate that the positions are discrete-time observations of an underlying continuous motion process. For population-based cell track analysis, it is common to align a set of tracks to a common starting point.

Definition 2 (Aligned cell track sets). For a cell track T = (x(1), t(1)) , ..., (x(k), t(k)), the track that results from subtracting the first element from all elements in the sequence, i.e.

$$T_0=\left(0,0\right),\left(x\left(2\right)-x\left(1\right),t\left(2\right)-t\left(1\right)\right),....,\left(x\left(k\right)-x\left(1\right),t\left(k\right)-t\left(1\right)\right)$$

is called the zero-aligned version of T .

Throughout the paper, we will assume that all cell track sets are zero-aligned.

Definition 3 (Mean displacement and mean square displacement). Let x₁(t) , ..., x _n (t) be the positions of n zero-aligned cell tracks at some fixed time t. Then the mean displacement D(t) for S is defined by

$$D\left(t\right)=\frac{\left|\right|x_1\left(t\right)\left|\right|+\cdots +\left|\right|x_n\left(t\right)\left|\right|}{n}$$

where || · || is the d-dimensional euclidian norm or vector length, and the mean square displacement D²(t) is defined by

$$D^2\left(t\right)=\frac{\left|\right|x_1\left(t\right)|{|}^2+\cdots +|\left|x_n\left(t\right)\right|{|}^2}{n}$$

Beauchemin et al. [5] estimated the three parameters v_free, t_free and t_pause of their model by fitting it to mean displacement data: For each combination of parameters, 10⁶ simulations of the model were run, and the mean displacement data was recorded and averaged over all simulations. The resulting mean displacement curve was compared to real data from several publications. The parameter triplets were ranked according to their squared sum of errors to the real data (see Table 1). It was found that the values t_free = 2.0 min, t_pause = 0.5 min and v_free = 18.8 μ m/min gave the best agreement to the real data. These values are biologically reasonable and similar to those estimated by Miller et al. [1].

Table 1 Different parameters can lead to similar motility coefficients.

Microscopic random walk with bias

In this section, we will define the notion of a one-dimensional random walk used in this paper and derive its relation to a one-dimensional convection-diffusion process. These concepts can be easily generalized to higher dimensions under the assumption that the dimensions are pairwise uncorrelated, while stochastic independence in the strict sense is not necessary. This assumption holds for the particular random walks we analyze in this paper.

Definition 4 (One-dimensional random walk [9, 11]). Let ξ _i , $i\in \mathbb{N}$, be independent, identically distributed real-valued random variables with mean μ and variance σ². Then the sequence (S _t ), $i\in \mathbb{N}$, defined by

$$S_t={\displaystyle \sum _{i=1}^t}{\xi }_i$$

is called a random walk with step bias μ and step variance σ². If μ =0, (S _t ) is called an unbiased random walk.

Intuitively, the sequence (S _t ) is the random walk trajectory and the ξ _i are the individual steps. The above definition is a slightly generalized version of the commonly used entirely discrete random walk with fixed steps to the left and the right: We assume that each step is sampled from a (possibly continuous) random variable with finite first and second moments. The following corollary provides an important characteristic of this random walk:

Corollary 5 (Mean and variance of a random walk). For a random walk (S _t ) as given by Definition 4,

$$E\left[S_t\right]=t\mu andVar\left(S_t\right)=t{\sigma }^2.$$

Proof. Since the ξ _i are independent and thus uncorrelated, this is a direct consequence of the linearity of variance and expectation.    □

If the random walk is unbiased, Var(S _t ) can be geometrically interpreted as the expected mean square squared displacement of the particles at time t (Definition 3). This leads to an important characteristic of the unbiased random walk: In a particle ensemble, the particles will go nowhere on average, and their expected mean square displacement is linear in time. This characterization is often used in the immunological literature [2] to convey the difference between a random walk versus directed migration, where the root mean square displacement is linear in time (Figure 1).

Figure 1

Directed vs. random motion. Traveled distance (directed movement) and root mean square displacement (random movement) of directed and random motion plotted as functions of time (left) and square root of time (right). On the left hand side, directed motion yields a straight line, with a slope given by its speed. On the right hand side, random motion yields a straight line, with a slope given by its diffusion coefficient (also called motility coefficient in this paper).

For sufficiently long observations, we do not only know the variance and mean of a random walk, but can furthermore obtain a good approximation of the entire particle distribution from the central limit theorem.

Theorem 6 (Central limit theorem [9, 10]). Let (S _t ) be a random walk as given in Definition 4. Then as t → ∞, the distribution of the random variable

$$Z_t=\frac{S_t-t\mu }{\sqrt{t}}$$

converges to a normal distribution with mean 0 and variance σ² .

Thus, after a sufficiently long time, the distribution of the random walk (S _t ) at fixed times is approximately normal. Of course, the terms "sufficiently long" and "approximately normal" are rather vague. In the context of intravital two-photon imaging, which is limited to a finite, often axially thin imaging region, we can usually not safely assume that we are in the range of validity of the normal distribution approximation: We are observing only a few "steps" (prolonged periods of rather persistent motion followed by short pauses) of a random walk. Beltman et al. [3] and ourselves [7] have discussed the implications of this technical limitation for detecting directional bias and estimating the motility coefficient.

We now generalize the random walk to higher dimensions.

Definition 7 (Multi-dimensional random walk [7]). Let $d\in \mathbb{N}$. Suppose $\left({\xi }_i^{\left(1\right)},\dots ,{\xi }_i^{\left(d\right)}\right)$, $i\in \mathbb{N}$are independent, identically distributed ${ℝ}^d$-valued random vectors with ${\xi }_i^{\left(1\right)},\dots ,{\xi }_i^{\left(d\right)}$pairwise uncorrelated, $E\left[{\xi }_i^{\left(1\right)}\right]={\mu }^{\left(1\right)},\dots ,E\left[{\xi }_i^{\left(d\right)}\right]={\mu }^{\left(d\right)}$and $Var\left({\xi }_i^{\left(1\right)}\right)=\dots =Var\left({\xi }_i^{\left(d\right)}\right)={\sigma }^2$. Then the sequence (S _t ), $t\in \mathbb{N}$, defined by

$$S_t=(S_t^{\left(1\right)},\dots ,S_t^{\left(d\right)}={\displaystyle \sum _{i=1}^t}\left({\xi }_i^{\left(1\right)},\dots ,{\xi }_i^{\left(d\right)}\right)$$

is called a d-dimensional random walk with step bias (μ⁽¹⁾, ..., μ^(d)) and step variance σ². If μ⁽¹⁾ = ... = μ^(d) = 0, then (S _t ) is called an unbiased d-dimensional random walk.

The following theorem gives the limit distribution of suitably scaled multi-dimensional random walks within our framework.

Theorem 8 (Multi-dimensional central limit theorem [7]). Let (S _t ) be a d-dimensional random walk as given in Definition 7. Then as t → ∞, the distribution of the random vector

$$Z_t=\frac{\left(S_t^{\left(1\right)}-t{\mu }^{\left(1\right)},\dots ,S_t^{\left(d\right)}-t{\mu }^{\left(d\right)}\right)}{\sqrt{t}}$$

converges to a d-dimensional normal distribution with zero means and covariance matrix

$$\sum =\left(\begin{array}{ccc}\hfill {\sigma }^2\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\sigma }^2\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {\sigma }^2\hfill \end{array}\right).$$

Proof. Let $a_1,\dots ,a_d\in ℝ$. According to the Cramér-Wold Theorem, it suffices to show that $\left(a_1\left(S_t^{\left(1\right)}-t{\mu }^{\left(1\right)}\right)+\dots +a_d\left(S_t^{\left(d\right)}-t{\mu }^{\left(d\right)}\right)\right)/\sqrt{t}$ converges in distribution to a normal random variable with zero mean and variance $\left(a_1^2+\dots +a_d^2\right){\sigma }^2$ (see [12]). Note that the random variables $a_1{\xi }_i^{\left(1\right)}+\dots +a_d{\xi }_i^{\left(d\right)}$, $i\in \mathbb{N}$, are independent, identically distributed and have mean a₁μ⁽¹⁾ + ... + a _d μ^(d). Furthermore, since ${\xi }_i^{\left(1\right)},\dots ,{\xi }_i^{\left(d\right)}$ are pairwise uncorrelated,

$$\mathsf{\text{Var}}\left(a_1{\xi }_i^{\left(1\right)},\dots ,a_d{\xi }_i^{\left(d\right)}\right)=\left(a_1^2+\dots +a_d^2\right){\sigma }^2.$$

Hence, the result follows from the central limit theorem for one-dimensional random walks.    □

The central limit argument shows that even complex cell migration processes might have a simple description when observed for a sufficiently long time: The diffusion coefficient is a function of all "microscopic" cell motility parameters such as persistence time, speed, turning angle distributions and so on. In the scaling limit, this single coefficient describes an unbiased random walk completely. Therefore, it is more common to describe random walk processes on such timescales in terms of a differential equation instead of a stochastic process. This can be done using a scale transition, as we will explain below.

The convection-diffusion equation

The random walk (S _t ) introduced in Definition 4 is a plausible microscopic model for the position of particles. However, if the positions are observed on large time scales (e.g., only every 20th state of the random walk can be observed due to a limited temporal resolution), then a macroscopic continuous model can be used to describe the migration. This continuous model no longer depends on the detailed characteristics of the microscopic model such as the turning angle distribution or persistence length (i.e., the precise shape of the ξ _i ).

Mathematically, we construct the continuous model as follows: For $n\in \mathbb{N}$ let (Z_n(t)) with t ≥ 0 be given by

$$Z_n\left(t\right)=\frac{1}{\sqrt{n}}{\displaystyle \sum _{i=1}^{⌊tn⌋}}\left({\xi }_i-\mu \right)+t\mu .$$

(1)

We write the time index t in brackets in order to emphasize that the processes are given in continuous time. The process (Z _n (t)) is obtained by observing n steps of the random walk (S _t ) in every time interval (t − 1, t] with $t\in \mathbb{N}$. It is scaled such that E [Z _n (t)] = E [S _t ] = tμ and Var(Z _n (t)) = Var(S _t ) = tσ² for every $t\in \mathbb{N}$. Between two successive steps, the process is constant.

Let (Z(t)), t ≥ 0, denote the limiting process obtained as n, the number of steps per time interval, tends to infinity. According to the Functional central limit theorem (sometimes also referred to as Donsker's Theorem, see [12]), (Z(t)) has the following two properties:

• ♦ Any finite-dimensional distribution of (Z(t)) is Gaussian.

• ♦ E[Z(t)] = tμ and Cov(Z(s), Z(t)) = min(s, t)σ² for s, t ≥ 0.

A process with these properties is called Brownian motion with drift μ. If t> 0, then Z(t) has a probability density ϕ(x, t) given by

$$\varphi \left(x,t\right)=\frac{1}{\sqrt{2\pi t\sigma }}\text{exp}\left(-\frac{{\left(x-t\mu \right)}^2}{2t{\sigma }^2}\right)$$

for $x\in ℝ$.

As one can easily verify, this probability density function solves the following well-known differential equation for the initial condition ϕ(x, 0) = δ(x), where δ denotes the Dirac delta distribution.

Theorem 9 (Convection-diffusion equation [11, 13]). Let C = μ and M = σ²/2. Then

$$\frac{\partial }{\partial t}\varphi \left(x,t\right)=\left[-C\frac{\partial }{\partial x}+M\frac{{\partial }^2}{\partial x^2}\right]\varphi \left(x,t\right).$$

Taken the initial condition given above, the solution of the convection-diffusion-equation exists and is unique [13].

As the sum of independent Brownian motions is again a Brownian motion, the convection-diffusion equation describes the dynamics of single diffusing particles as well as of ensembles of particles carrying out motions independently: Assuming the ensemble is large enough such that stochastic perturbations can be ignored, the convection-diffusion-equation describes the evolution of the particle concentration over time. This scaling argument justifies the use of the convection-diffusion equation for simulating large cell populations on timescales that are substantially longer than the rhythmicity of the cell migration process. For example, we used the convection-diffusion equation to simulate the transit of T cells through the lymph node, which takes approximately half a day [7].

Note that in the context of molecular motion, the quantity M is usually called diffusion coefficient. Throughout this paper we will use the term motility coefficient to indicate that we are talking about particles that are substantially larger than molecules; this is common in the immunological literature since Miller et al.'s early work [1].

Similarly as for the discrete model, the continuous model can be generalized to multiple dimensions. In this case, (Z _n (t)) is obtained by observing n steps of a d-dimensional random walk (S _t ) in every time interval of unit length. For n tending to infinity, the Functional central limit theorem yields convergence to a d-dimensional Brownian motion with independent components and convection coefficients C⁽¹⁾ = μ⁽¹⁾ , ..., C^(d)= μ^(d). The motility coefficient is equal to M = dσ² /2, where d is the number of observed dimensions.

Summary

As a microscopic model for stochastic motility of particles, we have first introduced one-dimensional random walks. According to this model, particles move by taking independent and identically distributed steps at equidistant discrete time points. Therefore, just like in the Beauchemin model, our framework does not consider persistence to last longer than the duration of a single step; the model is fully characterized by the absolute time between successive steps and the distribution of steps. Important characteristics of the step distribution are the step bias μ and the step variance σ². The random walk is called unbiased if μ = 0, and biased, otherwise. If we observe many particles carrying out simultaneous unbiased random walks, the particles will go nowhere on the average, but they will spread out increasingly, their distribution being approximately normal. The model can be easily generalized to d dimensions if we assume that the steps are dimension-wise uncorrelated (not necessarily statistically independent). Essentially, we can then treat the d dimensions like d separate one-dimensional random walks.

On the macroscopic scale, (biased) random walks can be approximated by Brownian motion (with drift). In this case the motion of particles is described sufficiently well by only two parameters, namely the motility coefficient M = σ² /2 and the convection coefficient C = μ. The motility coefficient M quantifies the random component of a biased random walk, while the convection coefficient C (which can be interpreted as a velocity) characterizes the directed component, if there is any (Figure 2). Convergence to Brownian motion can also be obtained for much more general random walk models using appropriate central limit theorems. However, the particular random walk chosen here is general enough for our purpose of analyzing the Beauchemin model.

Figure 2

Position distributions for biased random walks. Expected distribution of zero-aligned cell displacements measured in one dimension in the case of an unbiased random walk with motility coefficient M (left) and a biased random walk with convection coefficient C (right) as per theorem 9. In both cases, the distribution is Gaussian with standard deviation $\sqrt{2tM}$. In the biased case, the mean of the distribution moves away from the origin with time.

The Beauchemin model

In this section, we analytically derive three quantitative characteristics of the particular random walk defined by the Beauchemin model, namely the motility coefficient, the mean square displacement, and the confinement ratio.

Motility coefficient

The fact that the model is three-dimensional might seem to complicate our analysis. However, fortunately, it is possible to decouple the dimensions of the model, as will be shown in the following. The upcoming result was stated by ourselves [7] and, more recently, by Donovan and Lythe [14].

Proposition 10 (Motility Coefficient of the Beauchemin model [7, 14]). Let $\left({\xi }_i,v_i,{\zeta }_i\right)\in {ℝ}^3$, $i\in \mathbb{N}$, be independent random vectors distributed uniformly on a sphere with radius r. Set $\left(X_t,Y_t,Z_t\right)={\sum }_{i=1}^t\left({\xi }_i,v_i,{\zeta }_i\right)$. Then X _t , Y _t and Z _t are unbiased random walks as defined in 4, and as t → ∞, the distribution of $\left(X_t,Y_t,Z_t\right)/\sqrt{t}$ converges to a 3-variate normal distribution with zero mean and covariance matrix

$$\sum =\left(\begin{array}{ccc}\hfill {\sigma }^2\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\sigma }^2\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {\sigma }^2\hfill \end{array}\right),{\sigma }^2=r^2/3.$$

Proof. The key ingredients to the proof are

(i) ξ _i , υ _i , ζ _i are uniformly distributed on [−r, r].

(ii) Cov(ξ _i , υ _i ) = Cov(ξ _i , ζ _i ) = Cov(υ _i , ζ _i ) = 0.

In order to establish (i), note that P (ξ _i < h), i.e., the probability that ξ _i lies below the height h in the sphere, is proportional to the part of the sphere's surface area lying below the plane ξ _i = h. The surface of a spherical cap is proportional to its height, and thus ξ _i must be uniformly distributed on [−r, r] (see [15]). By isotropy, the same arguments hold for υ _i and ζ _i . The statement (ii) on pairwise uncorrelatedness is obtained by (ii) is obtained by observing that (ξ _i , υ _i ) is uniformly distributed on a circle with radius r − |z| given that ζ _i = z, and hence Cov(ξ _i , υ _i ) = 0.

Now, (i) implies E [ξ _i ] = E [υ _i ] = E [ζ _i ] = 0, and hence (X _t ), (Y _t ) and (Z _t ) are unbiased random walks. The statement on the asymptotic distribution of (X _t , Y _t , Z _t ) follows by (ii) together with Theorem 8, and the fact that a random variable with uniform distribution on [−r, r] has variance r²/3.    □

Hence, particles in the Beauchemin model behave exactly like the superposition of three uncorrelated random walks with uniformly distributed steps.

Mean square displacement

Putting together our result from the last section with the relationship between the motility coefficient M and the dimension-wise step variance σ² of a random walk, we obtain the following relation between the diffusion coefficient and the microscopic parameters:

$$M=\frac{{\left(v_{\text{free}}\times t_{\text{free}}\right)}^2}{6\left(t_{\text{free}}+t_{\text{pause}}\right)}$$

(2)

The corresponding diffusion coefficients for the 5 best-fitting parameter triplets determined by Beauchemin et al. are given in Table 1. It turns out that all the best-fitting triplets have very similar motility coefficients. On the other hand, the microscopic motility parameters vary quite a lot: For example, the pause time of the 4th best fitting triplet is 5 times longer than that of the 3rd best fitting triplet. The triplets were ranked by Beauchemin et al. according to the sum of squared residuals (SSR) to the experimental data. However, the SSR of the 5th best fitting triplet is only less than 4% higher than that of the best fitting triplet. Our upcoming analysis will provide some insight why this occurs.

The next two propositions will yield a precise characterization of the variance - and with it the mean square displacement - of the position of single particles and ensembles of observed particles in the Beauchemin model. It is mathematically more convenient to analyze the mean square displacement rather than the mean displacement because, as we have shown above, the dimensions can be decoupled. For mean displacement, this does not hold, and we would have to take into account the stochastical dependence between dimensions.

We start out with the variance of the position of one particle that starts its random walk at time t = 0. For simplicity, we will derive the following equations assuming that one dimension is observed. Due to the independence of the dimensions discussed above, multi-dimensional versions of these results are obtained simply by multiplying the derived quantity (i.e., mean square displacement or confinement ratio) with the number of observed dimensions.

Proposition 11 (Expected square displacement of single particles in the Beauchemin model). A particle starts at time t = 0 in position x(0) = 0 to carry out a random walk as defined by the Beauchemin model with fixed parameters t _pause ,t _free ,v _free . Let x(t + τ) denote the position of the particle observed in one dimension at time t + τ where t is an integer multiple of t _free + t _pause and 0 ≤ τ < t _free + t _pause . Then the expected square displacement of the particle is given by:

$$E\left[D^2\left(t+\tau \right)\right]=Var\left(x\left(t+\tau \right)\right)=2Mt+\frac{v_{free}^2}{3}{\left(\text{max}\left(\tau -t_{pause},0\right)\right)}^2$$

with M as defined in equation (2).

Proof. The term 2Mt follows directly from Proposition 10 and the linearity of variance (see also Corollary 5). Thus we continue with t = 0. Then we have from Proposition 10 that the position at time t_free + t_pause is uniformly distributed on the interval [−v_freet_free, v_freet_free], having variance 2M(t_free + t_pause).

For τ ≤ t_pause, the particle's position does not change, hence its variance is equal to 0. For τ > t_pause, the position is uniformly distributed on an interval with a size proportional to τ − t_pause. Hence the variance has quadratic form y = a(τ − t_pause)². Inserting τ = t_free + t_pause, we obtain

$$2M\left(t_{\text{free}}+t_{\text{pause}}\right)=a{t}_{\text{free}}^2$$

Inserting the definition of M from above, this results in

$$a=\frac{{\left(v_{\text{free}}\right)}^2}{3}$$

(3)

Combining the cases τ ≤ t_pause and τ > t_pause, we arrive at the proposition.    □

A plot of the expected square displacement shows a series of iterated pulses that approaches the linear progression described by the diffusion equation from below (see Figure 3). However, we will most likely not see such phasic behaviour if we take the average of several observed particles, since the particles will in general not be synchronized and the pulses that are caused by the free runs between turns will average out. Consequently, the next step in our analysis is to derive the expected mean square displacement for a particle ensemble. We do this by assigning to each observed particle a phase, which reflects the particle's state at the time when we start our observation. Then we average over all possible phases.

Figure 3

The Beauchemin model and the diffusion equation. The expected square displacement (or position variance) of a single particle along with the approximation by the diffusion equation.

Proposition 12 (Empirical mean square displacement of the Beauchemin model). Set t _run = t _free + t _pause . A particle starts at a uniformly distributed random time −t₀, −t _run <−t₀ ≤ 0 in some position to carry out a random walk as defined by the Beauchemin model with fixed parameters t _pause , t _free , v _free . Let $x\left(t\right)\in ℝ$ denote the position of this particle observed in one dimension at time t, where we choose the coordinate system such that x(0) = 0. Then the expected square displacement of the particle at time t is given by

$$\begin{array}{ll}\hfill E\left[D^2\left(t\right)\right]& =\frac{1}{t_{run}}{\int }_{-t_{run}}^{0}Var\left(x\left(t-t_0\right)-x\left(t_0\right)\right)d{t}_0\\ =2Mt-2M{t}_{free}\times \left\{\begin{array}{cc}\frac{1}{3}\hfill & t\ge t_{free}\hfill \\ \frac{1}{3}{\left(\frac{t}{t_{free}}\right)}^3-{\left(\frac{t}{t_{free}}\right)}^2+\left(\frac{t}{t_{free}}\right)\hfill & t<t_{free}\hfill \end{array}\right..\end{array}$$

Proof. See Figure 4 for an illustration of how the parameter δ affects the variance of the observed particle position during the initial observation time. For δ ≤ 0, the cell is resting at the beginning of observation time and remains paused for time t_pause − δ. For δ > 0, the cell is moving at the beginning of observation time and keeps moving for time δ. Let us write Var(S _δ (t)) for the position variance of a particle with phase shift δ.

Figure 4

Effect of the phase parameter on the variance. The plots illustrate how the phase parameter δ affects the variance of a single observed cell in the Beauchemin model. The expected variance for a (large) ensemble of observed particles is obtained by averaging over all values of δ.

Denote by τ₀ the time that the particle was already moving when we started observing it. Hence,

$$\widehat{x}\left(t\right)=x\left(t+{\tau }_0\right)-x\left({\tau }_0\right)$$

where $\widehat{x}\left(t\right)$ denotes our observed position. We are interested in the variance $Var\left[\widehat{x}\left(t\right)\right]$, for which the following holds:

$$\mathsf{\text{Var}}\left[\widehat{x}\left(t\right)\right]$$

(4)

$$=\mathsf{\text{Var}}\left[x\left(t+{\tau }_0\right)-x\left({\tau }_0\right)\right]$$

(5)

$$=\mathsf{\text{Var}}\left[x\right(t+{\tau }_0\left)\right]+\mathsf{\text{Var}}\left[x\left({\tau }_0\right)]-2\mathsf{\text{Cov}}[x(t+{\tau }_0,x\left({\tau }_0\right)\right]$$

(6)

What we are deriving is the expectation of expression 6, which is a function of the random variable τ₀. Due to the markov property of the random walk, only the last observed step matters and we can take τ₀ as uniformly distributed over [0, t_free + t_pause). Thus,

$$\mathsf{\text{E}}\left[\mathsf{\text{Var}}\left[\widehat{x}\left(t\right)\right]\right]=\frac{1}{t_{\mathsf{\text{free}}}+t_{\mathsf{\text{pause}}}}{\int }_{\tau =0}^{t_{\mathsf{\text{free}}}+t_{\mathsf{\text{pause}}}}\mathsf{\text{Var}}\left[x\left(t+{\tau }_0\right)\right]-\mathsf{\text{Var}}\left[x\left({\tau }_0\right)\right]d{\tau }_0$$

We start by integrating the two variance terms. Some algebra yields

$$\begin{array}{l}\frac{1}{t_{\mathsf{\text{free}}}+t_{\mathsf{\text{pause}}}}{\int }_{{\tau }_0=0}^{t_{\mathsf{\text{free}}}+t_{\mathsf{\text{pause}}}}\mathsf{\text{Var}}\left[x\left(t+{\tau }_0\right)\right]+\mathsf{\text{Var}}\left[{\tau }_0\right]d{\tau }_0\\ =\frac{{\left(v_{\mathsf{\text{free}}}{t}_{\mathsf{\text{free}}}\right)}^2}{3\left(t_{\mathsf{\text{free}}}+t_{\mathsf{\text{pause}}}\right)}\left(\left(t+\frac{t_{\mathsf{\text{free}}}}{3}\right)+\frac{t_{\mathsf{\text{free}}}}{3}\right)\\ =2M\left(\left(t+\frac{t_{\mathsf{\text{free}}}}{3}\right)+\frac{t_{\mathsf{\text{free}}}}{3}\right)\end{array}$$

where M is the motility coefficient as defined in equation 2. For the covariance term it can be shown by some case distinctions that

$$\begin{array}{l}\frac{1}{t_{\mathsf{\text{free}}}+t_{\mathsf{\text{pause}}}}{\int }_{{\tau }_0=0}^{t_{\mathsf{\text{free}}}+t_{\mathsf{\text{pause}}}}\mathsf{\text{Cov}}\left[x\left(t+{\tau }_0\right),{\tau }_0\right]d{\tau }_0\\ =M\text{min}\left(t,t_{\mathsf{\text{free}}}\right)+M\underset{=0\mathsf{\text{for}}t\ge t_{\mathsf{\text{free}}}}{\underbrace{\left(\frac{\text{min}{\left(t,t_{\mathsf{\text{free}}}\right)}^2}{t_{\mathsf{\text{free}}}}+\frac{\text{min}\left(t,t_{\mathsf{\text{free}}}\right)3}{3\underset{\mathsf{\text{free}}}{\overset{2}{t}}}+\frac{2t_{\mathsf{\text{free}}}}{3}\right)}}\end{array}$$

Combining the variance and covariance terms as per expression 6, we obtain the claimed identity.    □

As illustrated in Figure 5, the empirical mean square displacement from time t_free onwards is thus perfectly described by the linear equation

Figure 5

Forming the population average. Expected mean square displacement of observed particle ensembles as of Proposition 12. The pulses of the individual variances (see Figure 4) average out, and the exact expression converges quickly to a linear form with slope α = 2M and intercept β = 2Mt_free/3. One degree of freedom from the individual variance is lost due to the averaging process.

$$\mathsf{\text{E}}\left[D^2\left(t\right)\right]=\alpha t+\beta$$

(7)

with slope

$$\alpha =2M=\frac{{\left(v_{\mathsf{\text{free}}}\right)}^2{\left(t_{\mathsf{\text{free}}}\right)}^2}{3\left(t_{\mathsf{\text{free}}}+t_{\mathsf{\text{pause}}}\right)}$$

(8)

and intercept

$$\beta =-\frac{2t_{\mathsf{\text{free}}}}{3}M=-\frac{{\left(v_{\mathsf{\text{free}}}\right)}^2{\left(t_{\mathsf{\text{free}}}\right)}^3}{9\left(t_{\mathsf{\text{free}}}+t_{\mathsf{\text{pause}}}\right)}.$$

(9)

This is an important difference to other models of mean square displacement such as Fürth's equation, which is frequently applied to cell migration data [16, 17]. This equation reads

$$\mathsf{\text{E}}\left[D^2\left(t\right)\right]=2M\left(t-P\left(1-e^{-t/P}\right)\right)$$

(10)

where P is a parameter called the persistence length. Hence, in the Beauchemin Model, the persistence of migration carries on only for the predefined time t_pause + t_free, after which it is immediately and completely lost. On the other hand, persistence in Fürth's equation gradually decays over time following an exponential term. This is why we previously used Fürth's equation to estimate the cell motility coefficient from short-term migration data [7]. However, whether this really leads to more accurate results remains to be determined.

A second notable consequence of the above result is that although the Beauchemin model has three microscopic parameters, its step-wise variance has only two degrees of freedom. Consequently, while the value of t_pause would become apparent when analyzing the track of a single simulated cell at a very high resolution, it does no longer matter when a large population of cells that migrate asynchronously is analyzed, as shown by Proposition 12: Given a fixed t_free, there are infinitely many combinations of t_pause and v_free that yield an identical motility coefficient and thereby an identical expected square displacement of the population. For example, the empirical mean square displacement of the best-fitting triplet (0.5, 2.0, 18.8) is identical (up to some rounding) to that of the triplet (9.32, 2.0, 40). Hence, we cannot expect to get reliable estimates for all three parameters t_free,t_pause and v_free from fitting the model to mean square displacement data. This is consistent with the observation by Beauchemin et al. that the parameter t_pause can always be set to 0 without much influence on the quality of the best fit.

However, Beauchemin et al. fitted their simulations to the mean displacement data, and the above result does not imply that the mean displacement too has just two degrees of freedom. Perhaps surprisingly, simulations show that this is in fact not the case (Figure 6): While mean square displacement plots generated with using the above parameter triplets and our analytical solutions are indeed consistent with each other, the mean displacement plots are significantly different. Hence, the mean displacement plot appears to preserve more information on the underlying migration process than does the mean square displacement plot, at least when applied to cell tracks simulated using the Beauchemin model.

Figure 6

Mean displacement can be more informative than mean square displacement. The pause parameter of the Beauchemin model is not essential for the mean square displacement, but for the mean displacement. Left: Root mean square displacement curves generated for two populations with t_free = 2 min and t_pause = 0.5 min, v_free = 18.8 μ m/min (red squares and error bars) or t_pause = 9.32 min, v_free = 40 μ m/min (blue circles and error bars). Simulations (N = 10000) with these two parameter yield the same results, as predicted by Proposition 12 (black line). Error bars show the standard error of the mean. Right: Although there is no difference with respect to the mean square displacement, the mean displacement curves of the two simulated populations differ.

Confinement ratio

The confinement ratio is a measure of track straightness. In the context of lymphocyte migration analysis, it is also called "chemotactic index" or "meandering index" [2, 3]. The confinement ratio is usually defined as the quotient of a path's displacement from the origin versus the complete length of the path. If the cell changes directions very often, the confinement ratio will tend to zero over time, while for directed migration without a random component, the confinement ratio is equal to 1.

Instead of the confinement ratio, we consider here the squared confinement ratio, which can be easily determined for the Beauchemin model from the previous result.

Definition 13 (Squared confinement ratio). Consider a particle moving according to the Beauchemin model and let D²(t) denote the squared distance between the particle position at time t and its starting point. Then we define the squared confinement ratio C²(t) as follows:

$$C^2\left(t\right)=\frac{D^2\left(t\right)}{t^2{v}_{free}^2}$$

By assuming t_pause = 0 and applying the results from the previous section, we can determine the expected confinement ratio of the Beauchemin model as follows:

$$\begin{array}{ll}\hfill \mathsf{\text{E}}\left[C^2\left(t\right)\right]& =\frac{1}{t^2{v}_{\mathsf{\text{free}}}^2}\left[2Mt-2M{t}_{\mathsf{\text{free}}}\times \left\{\begin{array}{cc}\frac{1}{3}\hfill & t\ge t_{\mathsf{\text{free}}}\hfill \\ \frac{1}{3}\hfill & {\left(\frac{t}{t_{\mathsf{\text{free}}}}\right)}^3-{\left(\frac{t}{t_{\mathsf{\text{free}}}}\right)}^2+\left(\frac{t}{t_{\mathsf{\text{free}}}}\right)\hfill & t<t_{\mathsf{\text{free}}}\hfill \end{array}\right.\right]\\ =\frac{1}{t^2}\left[\frac{t_{\mathsf{\text{free}}}}{3}t-\frac{t_{\mathsf{\text{free}}}^2}{3}\times \left\{\begin{array}{cc}\frac{1}{3}\hfill & t\ge t_{\mathsf{\text{free}}}\hfill \\ \frac{1}{3}\hfill & {\left(\frac{t}{t_{\mathsf{\text{free}}}}\right)}^3-{\left(\frac{t}{t_{\mathsf{\text{free}}}}\right)}^2+\left(\frac{t}{t_{\mathsf{\text{free}}}}\right)\hfill & t<t_{\mathsf{\text{free}}}\hfill \end{array}\right.\right]\\ =\frac{t_{\mathsf{\text{free}}}}{3t}-\frac{t_{\mathsf{\text{free}}}^2}{3t^2}\times \left\{\begin{array}{cc}\frac{1}{3}\hfill & t\ge t_{\mathsf{\text{free}}}\hfill \\ \frac{1}{3}\hfill & {\left(\frac{t}{t_{\mathsf{\text{free}}}}\right)}^3-{\left(\frac{t}{t_{\mathsf{\text{free}}}}\right)}^2+\left(\frac{t}{t_{\mathsf{\text{free}}}}\right)\hfill & t<t_{\mathsf{\text{free}}}\hfill \end{array}\right.\\ =\frac{t_{\mathsf{\text{free}}}}{3t}-\frac{1}{9}\times \left\{\begin{array}{cc}1\hfill & t\ge t_{\mathsf{\text{free}}}\hfill \\ \frac{t}{t_{\mathsf{\text{free}}}}-3+\frac{3t_{\mathsf{\text{free}}}}{t}\hfill & t<t_{\mathsf{\text{free}}}\hfill \end{array}\right.\\ =\frac{t_{\mathsf{\text{free}}}}{3t}-\frac{1}{9}\times \left\{\begin{array}{cc}1\hfill & t\ge t_{\mathsf{\text{free}}}\hfill \\ \frac{3t_{\mathsf{\text{free}}}^2+t^2}{t_{\mathsf{\text{free}}}t}-3\hfill & t<t_{\mathsf{\text{free}}}\hfill \end{array}\right.\\ =\dots \\ =\frac{1}{3\text{max}\left(t,t_{\mathsf{\text{free}}}\right)}\left(t_{\mathsf{\text{free}}}-\frac{t\times t_{\mathsf{\text{free}}}}{3\text{max}{\left(t,t_{\mathsf{\text{free}}}\right)}^2}\right)\end{array}$$

Note that $\mathsf{\text{E}}\left[C^2\left(0\right)\right]=\frac{1}{3}$, because the path length tv_free is measured in three dimensions while the expected observed square distance is projected to one dimension.

Beltman et al. [3] have argued that the confinement ratio is not well suited for comparisons between different experiments because it tends to zero over time, and is therefore sensitive to the duration of the experiment from which it was estimated. They suggested to instead normalize the confinement ratio by multiplying it with the square root of time, such that it converges to a constant value. Applying this normalization to the squared confinement ratio of the Beauchemin model, we obtain

$$t\times \mathsf{\text{E}}\left[C^2\left(t\right)\right]=\frac{t}{3\text{max}\left(t,t_{\mathsf{\text{free}}}\right)}\left(t_{\mathsf{\text{free}}}-\frac{t\times t_{\mathsf{\text{free}}}}{3\text{max}{\left(t,t_{\mathsf{\text{free}}}\right)}^2}\right),$$

(11)

which for t ≥ t_free simplifies to

$$t\times \mathsf{\text{E}}\left[C^2\left(t\right)\right]=\frac{t_{\mathsf{\text{free}}}\left(3t-1\right)}{9t}.$$

(12)

Thus, the normalized squared confinement ratio converges to t_free/3. When there is evidence that the observed migration process really behaves roughly like the Beauchemin model, this could be a good rule of thumb for quick estimations of the persistence time from experimental data.

Simulating biased migration

For our recent study of the detection limits of two-photon imaging [7], we extended the Beauchemin model for simulating three different types of biased migration, for which we used the term taxis modes [18] to indicate that such biased migration in cells usually occurs in response to an external stimulus. These three taxis modes were orthotaxis, where the movement speed is adjusted, klinotaxis, where the duration of the "free runs" is changed, and topotaxis, where the turning angle distribution of the simulated cell is changed. Thus, these modifications can be viewed as semi-mechanistic ones that attempt to describe ways in which a cell could re-adjust its migration behaviour from pure random migration to biased random migration in response to an external stimulus. Before we review these modifications, we start by presenting a fourth method to simulate biased migration which is not mechanistically motivated, but is much easier to analyze.

In all cases, taxis is spatially and temporarily uniform and controlled by two parameters: The vector $\overrightarrow{b}\in {ℝ}^3$, assumed to have unitary length, sets the direction of the taxis, and the parameter p ∈ [0, 1] determines the magnitude of taxis, with p = 0 equivalent to an unbiased random walk and p = 1 to a maximally biased random walk; the precise meaning of the parameter p will be defined for each case below.

Simple phenomenological model

A simple ad-hoc method to extend the Beauchemin model for simulating biased migration is the following: During the pause phase, the cell is no longer kept stationary, but moves into the bias direction $\overrightarrow{b}$ with speed p · v_free. This model is very simple to analyze, because the modification is deterministic and so changes only the mean, but not the variance of each step. Hence, the motility coefficient remains the same.

Proposition 14 (Taxis speed for the simple model). The length of the convection coefficient, ||C||, of the simple biased migration model is given by

$$∥C∥=\frac{p\cdot v_{free}\cdot t_{pause}^2}{t_{free}+t_{pause}}.$$

Proof. During each pause phase, the cell moves a distance p · v_free · t_pause into the bias direction. ||C|| is then simply the product of this distance and the total fraction of time that the cell spends in the pause phase.    □

There might be some use for the simple biased migration model in situations where it is necessary to study the effect of bias independently of the motility coefficient. However, it is difficult to imagine how such a situation could arise biologically. Below we therefore discuss model modifications that have a more realistic motivation. For each of these modifications, it is possible to analytically derive the convection coefficient, i.e., the expected speed of the bias.

Throughout the rest of this section, we let the vector $\overrightarrow{d}$ denote the orientation vector of the Beauchemin model (i.e., a vector randomly sampled from the unit sphere), and $\overrightarrow{b}$ the bias direction. Moreover, we let $\left({\xi }^{\prime },v^{\prime },{\zeta }^{\prime }\right)$ denote the vector of random variables representing a step of the modified model, while (ξ, υ, ζ) denotes a step in the original model.

Orthotaxis model

In the orthotaxis model, the migration speed is no longer constant but depends on the chosen target direction. The more consistent this direction is with the target direction, the faster the simulated cell will move. Such bias could occur for instance in response to some external dragging force, and might be illustrated by imagining a random stroll through a city on a very windy day.

Formally, we define the adjusted speed as follows:

$$v_{\mathsf{\text{free}}}^{\prime }=v_{\mathsf{\text{free}}}\left(1+p〈\overrightarrow{b},\overrightarrow{d}〉\right)$$

Here 〈·,·〉 denotes the scalar product and p the bias strength parameter. Note that the movement perpendicular to $\overrightarrow{b}$is not affected at all, while movement against the direction $\overrightarrow{b}$ is slower and in direction $\overrightarrow{b}$ it is faster. More precisely, the relation between p and the bias speed v_taxis is given by the following proposition.

Proposition 15 (Speed of orthotaxis [7]). The convection coefficient ||C|| of the orthotaxis model is given by

$$∥C∥=v_{free}\cdot \frac{p{t}_{free}}{3\left(t_{free}+t_{pause}\right)}.$$

Proof. We can assume without loss of generality that the bias is along the first dimension, i.e. $\overrightarrow{d}=\left(1,0,0\right)$, and thus $〈\overrightarrow{b},\overrightarrow{d}〉=\xi$. We are interested in the expectation $\mathsf{\text{E}}\left[{\xi }^{\prime }\right]$; it is easy to see that the expectations $\mathsf{\text{E}}\left[v^{\prime }\right]$ and $\mathsf{\text{E}}\left[{\zeta }^{\prime }\right]$ are still 0. Then we get

$$\begin{array}{ll}\hfill {\xi }^{\prime }& ={v^{\prime }}_{\mathsf{\text{free}}}\xi =v_{\mathsf{\text{free}}}\left(1+p\xi \right)\xi \\ =v_{\mathsf{\text{free}}}\left(\xi +p{\xi }^2\right).\end{array}$$

Because ξ is uniformly distributed on [−1, 1] (see proof of Proposition 10), the expected location $\mathsf{\text{E}}\left[{\xi }^{\prime }\right]$ at the end of the step is

$$\begin{array}{ll}\hfill \mathsf{\text{E}}\left[{\xi }^{\prime }\right]& =\mathsf{\text{E}}\left[v_{\mathsf{\text{free}}}\cdot t_{\mathsf{\text{free}}}\left(\xi +p{\xi }^2\right)\right]\\ =v_{\mathsf{\text{free}}}\cdot t_{\mathsf{\text{free}}}\left(\mathsf{\text{E}}\left[\xi \right]+p\mathsf{\text{E}}\left[{\xi }^2\right]\right)\\ =v_{\mathsf{\text{free}}}\cdot t_{\mathsf{\text{free}}}\cdot p/3.\end{array}$$

Dividing this expected location $\mathsf{\text{E}}\left[{\xi }^{\prime }\right]$ by the step duration t_free + t_pause yields the claimed expression for ||C||.    □

Topotaxis model

In the case of topotaxis, the cell adjusts its turning behaviour but keeps the speed constant. This also affects both the variance and the mean of the cell step distribution. Topotaxis could be a likely bias mechanism for cells migrating on some sort of structural cell network in the tissue, like it was suggested to be the case for T cells in secondary lymphoid organs [19].

We implement topotaxis by skewing the distribution from which we pick the cell's orientation $\overrightarrow{d}=\left(\xi ,v,\zeta \right)$. Recall that this is a uniform distribution on a sphere in the unbiased case. Such a distribution can also be generated by sampling ξ uniformly from [−1, 1] and then picking υ, ζ uniformly from a circle with radius $\sqrt{1-{\xi }^2}$. Thus, the probability density function f_ξ(x) is defined as follows:

$$f_{\xi }\left(x\right)=\frac{1}{2}\mathsf{\text{ , }}x\in \left[-1,1\right].$$

For the modified model version, we define a parameterized version of f _ξ (x) as follows:

$$f_{{\xi }^{\prime }}\left(x\right)=\frac{1+px}{2}\mathsf{\text{ , }}x\in \left[-1,1\right].$$

Hence, for p = 0 we obtain the distribution of the unbiased case, while for p = 0.5 we get the slightly skewed distribution and for p = 1 we obtain the maximally skewed distribution .

This modification of the Beauchemin model yields the same convection coefficient as in the previous case.

Proposition 16 (Speed of topotaxis [7]). The speed ||C|| of the directional motion induced by the topotaxis model is given by the equation

$$∥C∥=v_{free}\cdot \frac{p{t}_{free}}{3\left(t_{free}+t_{pause}\right)}.$$

Proof. This is simply shown by calculating $\mathsf{\text{E}}\left[{\xi }^{\prime }\right]$ for the probability density function defined above.    □

Klinotaxis model

Lastly, we can also induce bias by modifying the duration t_free of the free run depending on the picked orientation $\overrightarrow{d}$. A well-known biological example of such bias is the "run-and-tumble" motion by which E. Coli searches for food: Because the bacterium is too small to sense concentration gradients in situ, it is capable of comparing concentration measurements along a persistent path until it figures out whether it is going in the right direction. When it senses that it is moving in the right direction (run), it will try to keep going, while otherwise it will stop its motion soon (tumble) and orient itself in a new random direction.

We implement klinotaxis by modifying the free run duration t_free depending on the angle between $\overrightarrow{d}$ and $\overrightarrow{b}$ as follows:

$$t_{\mathsf{\text{free}}}^{\prime }=t_{\mathsf{\text{free}}}\left(1+p〈\overrightarrow{b},\overrightarrow{d}〉\right)$$

Again, this modification leads to exactly the same bias speed as found for the previous two models.

Proposition 17 (Speed of klinotaxis [7]). The speed ||C|| of the directional motion induced by the klinotactic modification of the Beauchemin model is given by the equation

$$∥C∥=v_{free}\cdot \frac{p{t}_{free}}{3\left(t_{free}+t_{pause}\right)}.$$

Proof. The proof follows along the lines of a similar model of E. Coli run-and-tumble motion discussed by de Gennes [20]. We start again by analyzing a single free run of a cell projected onto the direction of bias, which is assumed without loss of generality to be along the first dimension ξ. Again we need to determine $\mathsf{\text{E}}\left[{\xi }^{\prime }\right]$. We start by writing ${\xi }^{\prime }$ as

$${\xi }^{\prime }=v_{\mathsf{\text{free}}}\cdot \xi \cdot t_{\mathsf{\text{free}}}\cdot \left(1+p\xi \right).$$

Now let us assume for a moment that t_pause = 0. It is easy to verify that $\mathsf{\text{E}}\left[t_{\mathsf{\text{free}}}^{\prime }\right]=t_{\mathsf{\text{free}}}$. Thus, the average speed towards the bias direction over many successive steps is $\mathsf{\text{E}}\left[{\xi }^{\prime }/t_{\mathsf{\text{free}}}\right]$. A simular calculation as in the proof of Proposition 15 yields

$$\mathsf{\text{E}}\left[{\xi }^{\prime }/t_{\mathsf{\text{free}}}\right]=v_{\mathsf{\text{free}}}\cdot p/3.$$

Now we arrive at the claimed expression by multiplying the above with the overall fraction of time that the particle spends in the free run phase, which is

$$t_{\mathsf{\text{free}}}/\left(t_{\mathsf{\text{free}}}+t_{\mathsf{\text{pause}}}\right).$$

   □

Note that the klinotaxis model no longer fits within our mathematical framework because the duration of a step in the model is no longer constant. Thus, the central limit theorem we applied in that section no longer applies, and a different central theorem would be needed to show formally that the model converges to a Brownian motion.

Moreover, we point out that for all modifications discussed in this section except for the simple phenomenological one, the random motility component changes. Therefore, these modifications cause as a side-effect a modified motility coefficient. In general, the random motility component are even affected in a non-isotropic fashion, such that the motility coefficient would need to be described using a 3 × 3 matrix instead of a single scalar value. In Figure 7, we show some simulation results that illustrate this point.

Figure 7

Taxis skews the random motility component. Displacement of in silico cells when applying the four different methods of simulating biased migration along the X axis using the methods discussed in this paper, i.e., the simple phenomenological model (black), orthotaxis (red), topotaxis (blue), and klinotaxis (green). For each panel, N = 1000 cells were simulated using the parameters t_pause = 0.5 min, t_free = 2.0 min, and v_free = 18.8 μ m/min. Cell displacements are shown after one hour of simulation time, and the taxis parameters p are in each case set such that the expected displacement is 300 μ m per hour. The length of the scale bars is proportional to the dimension-wise variances of the cell displacements, demonstrating that these variances change when biased migration is introduced.

Conclusions

In this article, we have compiled both old and new analytical results on the Beauchemin model of lymphocyte migration, which was originally a model of purely random motility, but is now able to accommodate partially biased migration as well. An important practical consequence of these results is that validation of the model against experimentally determined motility coefficient, mean square displacement data, and confinement ratio curves does no longer need to be performed by simulation, because these quantities can be calculated directly from the model parameters. This also allows for the use of standard numerical fitting procedures, such as implemented in popular computer algebra packages, to estimate the model parameters from such data by least-squares-fitting. An important observation was that when simulating a large population of unsynchronized cells, the model is equivalent to the ideal chain model of polymer physics, which is however only true with respect to the quantities named above and not for others like the mean displacement. Interestingly, this shows that we can sometimes infer information from mean displacement data that we cannot infer from mean square displacement data. However, an analytical derivation of the mean displacement of the model has not yet been achieved, and can be expected to be complicated judging from the experiences with similar models [21]. For now, we leave this as an important open problem.

Subtle biased migration of lymphocytes can point to important biological functions or phenomena [22, 23], and the use of mathematical models has helped to understand the quantitative impact of such effects [23, 7]. The modifications discussed above facilitate simulations of such subtle biases with pre-defined speed, again removing the need for fitting the bias strength parameter to data by using simulations. However, all but one of the discussed modifications also affect the random migration component, which is likely also true in biological scenarios where randomly migrating cells respond to an external stimulus. Therefore, another important task that remains to be solved is to analytically derive the motility coefficient matrices for the topotaxis, klinotaxis, and orthotaxis models.