 Proceedings
 Open Access
 Published:
Analytical results on the Beauchemin model of lymphocyte migration
BMC Bioinformatics volume 14, Article number: S10 (2013)
Abstract
The Beauchemin model is a simple particlebased 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 predefined 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 twophoton 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: Wholepopulation models, usually formulated using ordinary or partial differential equations, do not distinguish between individual cells, but treat cell subpopulations as continuous quantities. Particlebased models do consider each cell individually, but treat cells as freely moving particles without mass, volume, and shape. Finally, wholecell 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 particlebased model that was introduced by Beauchemin, Dixit, and Perelson [5]. Throughout the paper, we refer to this model as the "Beauchemin Model".
Particlebased approaches such as the Beauchemin model are used for studying questions where a wholepopulation approach does not provide sufficient information, but a wholecell 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 twophoton 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 threedimensional 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 selfcontained fashion. Thereafter, we present analytical results about random motion in the Beauchmin model, and discuss four different ways for simulating biased migration with predefined 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 convectiondiffusionequation.
First of all, we introduce some notational conventions. We denote the expectation of a random variable ξ by E[ξ], its variance by Var(ξ) = E [(ξ − E[ξ])^{2}] 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 realvalued. Mean and variance of random variables will also be written as μ and σ^{2}, 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 twophoton 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 ddimensional cell track is a finite sequence
consisting of positions x\left(1\right)\phantom{\rule{0.3em}{0ex}},\dots ,\phantom{\rule{0.3em}{0ex}}x\left(k\right)\in {\mathbb{R}}^{d}and increasing time points t\left(1\right),\dots ,t\left(k\right)\in \mathbb{R}.
The time indices are written in brackets to indicate that the positions are discretetime observations of an underlying continuous motion process. For populationbased 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.
is called the zeroaligned version of T .
Throughout the paper, we will assume that all cell track sets are zeroaligned.
Definition 3 (Mean displacement and mean square displacement). Let x_{1}(t) , ..., x_{ n } (t) be the positions of n zeroaligned cell tracks at some fixed time t. Then the mean displacement D(t) for S is defined by
where  ·  is the ddimensional euclidian norm or vector length, and the mean square displacement D^{2}(t) is defined by
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^{6} 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].
Microscopic random walk with bias
In this section, we will define the notion of a onedimensional random walk used in this paper and derive its relation to a onedimensional convectiondiffusion 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 (Onedimensional random walk [9, 11]). Let ξ_{ i }, i\in \mathbb{N}, be independent, identically distributed realvalued random variables with mean μ and variance σ^{2}. Then the sequence (S_{ t }), i\in \mathbb{N}, defined by
is called a random walk with step bias μ and step variance σ^{2}. 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,
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).
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
converges to a normal distribution with mean 0 and variance σ^{2} .
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 twophoton 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 (Multidimensional 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 {\mathbb{R}}^{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
is called a ddimensional random walk with step bias (μ^{(1)}, ..., μ^{(d)}) and step variance σ^{2}. If μ^{(1)} = ... = μ^{(d)} = 0, then (S_{ t }) is called an unbiased ddimensional random walk.
The following theorem gives the limit distribution of suitably scaled multidimensional random walks within our framework.
Theorem 8 (Multidimensional central limit theorem [7]). Let (S_{ t }) be a ddimensional random walk as given in Definition 7. Then as t → ∞, the distribution of the random vector
converges to a ddimensional normal distribution with zero means and covariance matrix
Proof. Let {a}_{1},\dots ,\phantom{\rule{0.3em}{0ex}}{a}_{d}\in \mathbb{R}. According to the CramérWold Theorem, it suffices to show that \left({a}_{1}\left({S}_{t}^{\left(1\right)}t{\mu}^{\left(1\right)}\right)+\phantom{\rule{2.77695pt}{0ex}}\dots \phantom{\rule{2.77695pt}{0ex}}+{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)}+\phantom{\rule{2.77695pt}{0ex}}\dots \phantom{\rule{2.77695pt}{0ex}}+{a}_{d}{\xi}_{i}^{\left(d\right)}, i\in \mathbb{N}, are independent, identically distributed and have mean a_{1}μ^{(1)} + ... + a_{ d }μ^{(d)}. Furthermore, since {\xi}_{i}^{\left(1\right)},\dots ,\phantom{\rule{2.77695pt}{0ex}}{\xi}_{i}^{\left(d\right)} are pairwise uncorrelated,
Hence, the result follows from the central limit theorem for onedimensional 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 convectiondiffusion 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
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σ^{2} 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 finitedimensional distribution of (Z(t)) is Gaussian.

♦ E[Z(t)] = tμ and Cov(Z(s), Z(t)) = min(s, t)σ^{2} 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
for x\in \mathbb{R}.
As one can easily verify, this probability density function solves the following wellknown differential equation for the initial condition ϕ(x, 0) = δ(x), where δ denotes the Dirac delta distribution.
Theorem 9 (Convectiondiffusion equation [11, 13]). Let C = μ and M = σ^{2}/2. Then
Taken the initial condition given above, the solution of the convectiondiffusionequation exists and is unique [13].
As the sum of independent Brownian motions is again a Brownian motion, the convectiondiffusion 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 convectiondiffusionequation describes the evolution of the particle concentration over time. This scaling argument justifies the use of the convectiondiffusion 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 convectiondiffusion 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 ddimensional 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 ddimensional Brownian motion with independent components and convection coefficients C^{(1)} = μ^{(1)} , ..., C^{(d)}= μ^{(d)}. The motility coefficient is equal to M = dσ^{2} /2, where d is the number of observed dimensions.
Summary
As a microscopic model for stochastic motility of particles, we have first introduced onedimensional 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 σ^{2}. 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 dimensionwise uncorrelated (not necessarily statistically independent). Essentially, we can then treat the d dimensions like d separate onedimensional 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} /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.
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 threedimensional 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 {\mathbb{R}}^{3}, i\in \mathbb{N}, be independent random vectors distributed uniformly on a sphere with radius r. Set \left({X}_{t},\phantom{\rule{2.77695pt}{0ex}}{Y}_{t},\phantom{\rule{2.77695pt}{0ex}}{Z}_{t}\right)={\sum}_{i=1}^{t}\left({\xi}_{i},{v}_{i},\phantom{\rule{2.77695pt}{0ex}}{\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 3variate normal distribution with zero mean and covariance matrix
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^{2}/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 dimensionwise step variance σ^{2} of a random walk, we obtain the following relation between the diffusion coefficient and the microscopic parameters:
The corresponding diffusion coefficients for the 5 bestfitting parameter triplets determined by Beauchemin et al. are given in Table 1. It turns out that all the bestfitting 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, multidimensional 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:
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_{free}t_{free}, v_{free}t_{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})^{2}. Inserting τ = t_{free} + t_{pause}, we obtain
Inserting the definition of M from above, this results in
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.
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_{0}, −t_{ run } <−t_{0} ≤ 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 \mathbb{R} 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
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 δ.
Denote by τ_{0} the time that the particle was already moving when we started observing it. Hence,
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:
What we are deriving is the expectation of expression 6, which is a function of the random variable τ_{0}. Due to the markov property of the random walk, only the last observed step matters and we can take τ_{0} as uniformly distributed over [0, t_{free} + t_{pause}). Thus,
We start by integrating the two variance terms. Some algebra yields
where M is the motility coefficient as defined in equation 2. For the covariance term it can be shown by some case distinctions that
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
with slope
and intercept
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
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 shortterm 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 stepwise 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 bestfitting 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.
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^{2}(t) denote the squared distance between the particle position at time t and its starting point. Then we define the squared confinement ratio C^{2}(t) as follows:
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:
Note that \mathsf{\text{E}}\phantom{\rule{2.77695pt}{0ex}}\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
which for t ≥ t_{free} simplifies to
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 twophoton 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 semimechanistic ones that attempt to describe ways in which a cell could readjust 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 {\mathbb{R}}^{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 adhoc 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
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},\phantom{\rule{2.77695pt}{0ex}}{v}^{\prime},\phantom{\rule{2.77695pt}{0ex}}{\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:
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
Proof. We can assume without loss of generality that the bias is along the first dimension, i.e. \overrightarrow{d}\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\left(1,\phantom{\rule{2.77695pt}{0ex}}0,\phantom{\rule{2.77695pt}{0ex}}0\right), and thus \u3008\overrightarrow{b},\overrightarrow{d}\u3009\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\xi. We are interested in the expectation \mathsf{\text{E}}\phantom{\rule{2.77695pt}{0ex}}\left[{\xi}^{\prime}\right]; it is easy to see that the expectations \mathsf{\text{E}}\phantom{\rule{2.77695pt}{0ex}}\left[{v}^{\prime}\right] and \mathsf{\text{E}}\phantom{\rule{2.77695pt}{0ex}}\left[{\zeta}^{\prime}\right] are still 0. Then we get
Because ξ is uniformly distributed on [−1, 1] (see proof of Proposition 10), the expected location \mathsf{\text{E}}\phantom{\rule{2.77695pt}{0ex}}\left[{\xi}^{\prime}\right] at the end of the step is
Dividing this expected location \mathsf{\text{E}}\phantom{\rule{2.77695pt}{0ex}}\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}\phantom{\rule{2.77695pt}{0ex}}=\phantom{\rule{2.77695pt}{0ex}}\left(\xi ,\phantom{\rule{2.77695pt}{0ex}}v,\phantom{\rule{2.77695pt}{0ex}}\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:
For the modified model version, we define a parameterized version of f_{ ξ } (x) as follows:
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
Proof. This is simply shown by calculating \mathsf{\text{E}}\phantom{\rule{2.77695pt}{0ex}}\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 wellknown biological example of such bias is the "runandtumble" 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:
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
Proof. The proof follows along the lines of a similar model of E. Coli runandtumble 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}}\phantom{\rule{2.77695pt}{0ex}}\left[{\xi}^{\prime}\right]. We start by writing {\xi}^{\prime} as
Now let us assume for a moment that t_{pause} = 0. It is easy to verify that \mathsf{\text{E}}\phantom{\rule{2.77695pt}{0ex}}\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}}\phantom{\rule{2.77695pt}{0ex}}\left[{\xi}^{\prime}/{t}_{\mathsf{\text{free}}}\right]. A simular calculation as in the proof of Proposition 15 yields
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
□
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 sideeffect a modified motility coefficient. In general, the random motility component are even affected in a nonisotropic 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.
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 leastsquaresfitting. 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 predefined 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.
References
 1.
Miller MJ, Wei SH, Parker I, Cahalan MD: Twophoton imaging of lymphocyte motility and antigen response in intact lymph node. Science. 2002, 296: 18691873. 10.1126/science.1070051.
 2.
Cahalan MD, Parker I: Choreography of cell motility and interaction dynamics imaged by twophoton microscopy in lymphoid organs. Annual Reviews in Immunology. 2008, 26: 585626. 10.1146/annurev.immunol.24.021605.090620.
 3.
Beltman JB, Marée AFM, de Boer RJ: Analysing immune cell migration. Nature Reviews Immunology. 2009, 9: 789798. 10.1038/nri2638.
 4.
Beltman JB, Marée AFM, Lynch JN, Miller MJ, de Boer RJ: Lymph node topology dictates T cell migration behavior. Journal of Experimental Medicine. 2007, 204: 771780. 10.1084/jem.20061278.
 5.
Beauchemin C, Dixit NM, Perelson AS: Characterizing T cell movement within lymph nodes in the absence of antigen. Journal of Immunology. 2007, 178 (9): 55055512.
 6.
Grigorova IL, Panteleev M, Cyster JG: Lymph node cortical sinus organization and relationship to lymphocyte egress dynamics and antigen exposure. Proceedings of the National Academy of Sciences of the United States of America. 2010, 107: 2044720452. 10.1073/pnas.1009968107.
 7.
Textor J, Peixoto A, Henrickson SE, Sinn M, von Andrian UH, Westermann J: Defining the Quantitative Limits of Intravital TwoPhoton Lymphocyte Tracking. Proceedings of the National Academy of Sciences of the United States of America. 2011, 108 (30): 1240112406. 10.1073/pnas.1102288108.
 8.
Flory PJ: Statistical Mechanics of Chain Molecules. 1969, Interscience Publishers
 9.
Berg HC: Random Walks in Biology. 1993, Princeton University Press
 10.
Kallenberg O: Foundations of modern probability. 2002, Springer
 11.
Murray JD: Mathematical Biology: I. An Introduction. 2002, Springer, 3
 12.
Billingsley P: Convergence of Probability Measures. 1999, Wiley, John & Sons, 2
 13.
Evans LD: Partial Differential Equations. 1998, American Mathematical Society
 14.
Donovan GM, Lythe G: Tcell movement on the reticular network. Journal of Theoretical Biology. 2012, 295: 5967.
 15.
Marsaglia G: Choosing a point from the surface of a sphere. Annals of Mathematical Statistics. 1972, 43 (2): 645646. 10.1214/aoms/1177692644.
 16.
Selmeczi D, Mosler S, Hagedorn PH, Larsen NB, Flyvbjerg H: Cell Motility as Persistent Random Motion: Theories from Experiments. Biophysical Journal. 2005, 89: 912931. 10.1529/biophysj.105.061150.
 17.
Grabher C, Cliffe A, Miura K, Hayflick J, Pepperkok R, Rørth P, Wittbrodt J: Birth and life of tissue macrophages and their migration in embryogenesis and inflammation in medaka. Journal of Leukocyte Biology. 2007, 81: 263271.
 18.
Ionides EL, Fang KS, Isseroff RR, Oster G: Stochastic models for cell motion and taxis. Journal of Mathematical Biology. 2003, 48: 2337.
 19.
Bajénoff M, Egen JG, Koo LY, Laugier JP, Brau F, Glaichenhaus N, Germain RN: Stromal cell networks regulate lymphocyte entry, migration, and territoriality in lymph nodes. Immunity. 2006, 25: 9891001. 10.1016/j.immuni.2006.10.011.
 20.
de Gennes PG: Chemotaxis: the role of internal delays. European Biophysics Journal. 2004, 33: 691693. 10.1007/s002490040426z.
 21.
Codling EA, Plank MJ, Benhamou S: Random walk models in biology. Journal of the Royal Society Interface. 2008, 5 (25): 813834. 10.1098/rsif.2008.0014.
 22.
Castellino F, Huang AY, AltanBonnet G, Stoll S, Scheinecker C, Germain RN: Chemokines enhance immunity by guiding naive CD8+ T cells to sites of CD4+ T celldendritic cell interaction. Nature. 2006, 440: 890895. 10.1038/nature04651.
 23.
Beltman JB, Allen CD, Cyster JG, de Boer RJ: B cells within germinal centers migrate preferentially from dark to light zone. Proceedings of the National Academy of Sciences of the United States of America. 2011
 24.
Mempel TR, Henrickson SE, von Andrian UH: Tcell priming by dendritic cells in lymph nodes occurs in three distinct phases. Nature. 2004, 427: 154159. 10.1038/nature02238.
 25.
Miller MJ, Hejazi AS, Wei SH, Cahalan MD, Parker I: T cell repertoire scanning is promoted by dynamic dendritic cell behavior and random T cell motility in the lymph node. Proceedings of the National Academy of Sciences of the United States of America. 2004, 101: 9981003. 10.1073/pnas.0306407101.
 26.
Miller MJ, Wei SH, Cahalan MD, Parker I: Autonomous T cell trafficking examined in vivo with intravital twophoton microscopy. Proceedings of the National Academy of Sciences of the United States of America. 2003, 100 (5): 26042609. 10.1073/pnas.2628040100.
Acknowledgements
We thank Joost Beltman for helpful discussions. JT was supported by the Netherlands Organisation for Scientific Research (NWO) grant 912.10.066 (to RJdB). The research of Mathieu Sinn was conducted while he was a Postdoctoral Research Fellow at the David R. Cheriton School of Computer Science, University of Waterloo, Waterloo, Ontario, Canada N2L 3G1.
Declarations
This article has been published as part of BMC Bioinformatics Volume 14 Supplement 6, 2013: Selected articles from the 10th International Conference on Artificial Immune Systems (ICARIS). The full contents of the supplement are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/14/S6.
Author information
Affiliations
Corresponding author
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
JT conceived the study, participated in the mathematical analysis, and wrote the paper; MS participated in the mathematical analysis and helped write the paper; RJdB participated in the study conception and helped write the paper.
Rights and permissions
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
About this article
Cite this article
Textor, J., Sinn, M. & de Boer, R.J. Analytical results on the Beauchemin model of lymphocyte migration. BMC Bioinformatics 14, S10 (2013). https://doi.org/10.1186/1471210514S6S10
Published:
DOI: https://doi.org/10.1186/1471210514S6S10
Keywords
 Random Walk
 Cell Track
 Convection Coefficient
 Lymphocyte Migration
 Functional Central Limit Theorem