Inference of random walk models to describe leukocyte migration

The situation we model consists of N non-interacting point particles in ${{\mathbb{R}}}^{2}$. These points represent the centroid of a cell, which is often used to track cells from live imaging data. The approximation of the cell as a point mass introduces simplicity and for many problems the added complexity of particle shapes or sizes is unnecessary. In this way our model is build on the assumption that the movement results from the cells interaction with its environment. In situations where the cell density is high or the cell's compartment is small, cell–cell-interactions need to be considered. The direction of a particle's movement at any time step t is described by two random variables, which are a step length s_t and a turning angle ${\alpha }_{t}$ taken with respect to an arbitrary fixed reference axis. The step length follows the distribution

$$\begin{eqnarray}&&{s}_{t}\sim \sqrt{{\rm{d}}t}\times {N}^{+}(0,1),\end{eqnarray} \tag{ 1 }$$

where ${N}^{+}(0,1)$ is a truncated normal distribution, truncated at 0, and dt = 0.001. The chosen step size model is just a practical implementation; our statistical analysis is independent of step sizes.

The turning angle ${\alpha }_{t}$ is defined as the angle between a motion vector and a reference axis, here the negative y-axis (see figure 1(a), at time t. ${\alpha }_{t}$ follows the wrapped normal distributions (Breitenberger 1963) with the probability density function

$$\begin{eqnarray}&&{N}_{w}({\alpha }_{t}| \mu ,\sigma )\\ &&=\displaystyle \frac{1}{\sigma \sqrt{2\pi }}\displaystyle \sum _{k=-\infty }^{\infty }\mathrm{exp}(-\displaystyle \frac{{({\alpha }_{t}-\mu +2\pi k)}^{2}}{2{\sigma }^{2}}).\end{eqnarray} \tag{ 2 }$$

We note that we have selected the wrapped normal distribution mainly for its flexibility and relative simplicity; the inference scheme described below is, in principle, applicable to any other choice of distribution defined on a circle (e.g. von Mises distribution).

Zoom In Zoom Out Reset image size

The mean μ and variance σ depend on whether the random walker follows a biased or persistent motion. For the biased motion we define $\mu =\beta$ (the direction of bias) and for the persistent distribution we have $\mu ={\alpha }_{t-1}$, which is the direction in the previous time step (see figure 1(a)). The variances σ for the biased and persistent motion are defined as

$$\begin{eqnarray}&&{\sigma }_{p}=-2\mathrm{log}(p)\end{eqnarray} \tag{ 3 }$$

and

$$\begin{eqnarray}&&{\sigma }_{b}=-2\mathrm{log}(b)\end{eqnarray} \tag{ 4 }$$

respectively, with the persistence and bias parameters, p and b, with $p,b\in [0,1]$. They affect the variance of the distributions such that the closer to 1 they are, the smaller their respective variances will be, and the more likely the particle will be to sample an angle in the direction of the bias or the persistence, as appropriate. If p or b is equal to 0, then the corresponding variance will tend to infinity and thus the distribution is a wrapped uniform distribution.

The decision whether the random walker follows biased or persistent motion is based on a further random variable, which follows a Bernoulli distribution with success probability w, so that w describes the probability of a biased motion and $1-w$ describes the probability of a persistent motion; w is assumed to remain constant over time t.

At each time point, ${\alpha }_{t}$ and s_t are determined and, accordingly, the particle moves a distance of s_t in the direction defined by ${\alpha }_{t}$. This constitutes one step in the particle's trajectory, described by three parameters (p, b and w). Five different parameter combinations are considered as reference points for specific types of random walks as can be seen in table 1. A particle that performs Brownian motion does not show either bias or persistence (p = 0, b = 0). The choice of the weight is then arbitrary and we chose without loss of generality w = 0.5. A persistent random walker has no bias (b = 0), but he has some level of persistence which is defined by p and w. Accordingly, a biased random walker has no persistence (p = 0) and the level of bias is defined by b and w. In the case of a biased persistent random walk (BPRW1 and BPRW2) none of the parameters should be 0. From these example parameter combinations we compute TMs, which we will refer to in the following as query TMs.

To summarize the model output in the approximate inference framework we extract the directional transitions from the simulation and compute from these a transition matrix, first introduced in (Taylor et al 2013), which allows for visualization of the dynamics exhibited by the particle trajectories.

A TM T is the expected joint distribution of successively measured turning angles $({\alpha }_{t},{\alpha }_{t+{\rm{d}}t})$. The components of the matrix T_ij gives the joint probability of measuring the angle ${\alpha }_{i,t}$ at some time t and the angle ${\alpha }_{j,t+{\rm{d}}t}$ at the next measurement time $t+{\rm{d}}t$, where ${\alpha }_{i,t}$ represents any angle that lies in the half-open interval $[2\pi i/{N}_{\mathrm{bins}},2\pi (i+1)/{N}_{\mathrm{bins}})$ for $i=0,1,...{N}_{\mathrm{bins}}$.

T_ij can be estimated from the data by constructing a two-dimensional histogram. First, we bin successive (and overlapping) pairs of angles of motions measured for every particle in terms of the intervals above. Then we define the sample TM as

$$\begin{eqnarray}&&{\widehat{T}}_{{ij}}=\displaystyle \frac{1}{N^{\prime} }\displaystyle \sum _{q=1}^{N}\displaystyle \sum _{r=1}^{{T}_{q}/{\rm{d}}t}{I}_{{qr},{ij}},\end{eqnarray} \tag{ 5 }$$

where

$$\begin{eqnarray}{I}_{{qr},{ij}}=\left\{\begin{array}{ll}1, & \begin{array}{l}\mathrm{if}\;{\alpha }_{i,r*{\rm{d}}t}^{(q)}\in \left[\frac{2\pi i}{{N}_{\mathrm{bins}}},\frac{2\pi (i+1)}{{N}_{\mathrm{bins}}}\right)\\ \mathrm{and}\;{\alpha }_{i,(r+1)\ast {\rm{d}}t}^{(q)}\in \left[\frac{2\pi j}{{N}_{\mathrm{bins}}},\frac{2\pi (j+1)}{{N}_{\mathrm{bins}}}\right)\end{array}\\ 0, & \mathrm{otherwise},\end{array}\right.\end{eqnarray} \tag{ 6 }$$

and $N^{\prime} ={\sum }_{q=1}^{N}{T}_{q}/{\rm{d}}t$ the total number of successive angle-pair observations. The choice of ${N}_{\mathrm{bins}}$ is not pre-determined and can, in principle, be tuned for different datasets. The more transitions are observed in a data set the larger ${N}_{\mathrm{bins}}$ can be chosen.

As can be seen in the key in figure 1(b), we have adopted the convention whereby the columns of the TM indicate the direction intervals at time t while the rows the intervals at time $t+{\rm{d}}t$.

The interval between experimental observations dt_obs is determined by the experimental setup. In the general case, the interval between directional changes d t is an independent random variable; specifically, the number of directional changes that take place between any two observations is itself an unobserved quantity. It is precisely this hidden feature that prevents one from deriving exact closed-form expressions for the parameter likelihoods, thereby necessitating an approximate approach. In this scenario, the random walk model adopted here is an effective model (see (Sim et al 2015) for alternative velocity jump process models); concomitantly, the directional terms ${\alpha }_{t}$ are simply functions of the displacement data rather than the true ballistic paths of the particles.

Nevertheless in the simulation examples below we let

$$\begin{eqnarray}&&{\rm{d}}t={\rm{d}}{t}_{\mathrm{obs}}=\mathrm{constant}.\end{eqnarray} \tag{ 7 }$$

This simplifying assumption is made for the sole purpose of validating the proposed ABC approach as it allows us to determine the exact path probabilities and, hence, parameter posteriors for comparing against the approximate versions.

For the set of model parameters Θ, the likelihood resulting from the observation of N paths is simply

$$\begin{eqnarray}&&L(\Theta )=\displaystyle \prod _{i=1}^{N}p({\pi }_{i}\;| \;\Theta ),\end{eqnarray} \tag{ 8 }$$

where ${\pi }_{i}=({\alpha }_{i,0},...,{\alpha }_{i,T})$ is the ith path consisting of a sequence of $T+1$ observations ${\alpha }_{i,t}\in {{\mathbb{R}}}^{2}$. From equation (7) the individual terms can be factorized as

$$\begin{eqnarray}&&p({\pi }_{i}\;| \;\Theta )=p({\alpha }_{i,0}\;| \;\Theta )\displaystyle \prod _{t=1}^{T}p({\alpha }_{{i}_{t}}\;| \;{\alpha }_{i,t-1},\Theta ),\end{eqnarray} \tag{ 9 }$$

where the individual probability terms are simply given by equation (2).

Even though for the proposed model in this study it is straightforward to define the likelihood function and perform exact inference, the motivation of this study is to provide an inference scheme that can also be applied to far more complex random walk models, for which the likelihood is not tractable. The approach taken for the parameter inference is ABC in a sequential Monte Carlo sampling scheme (ABC–SMC) (figure 1(c)) (Toni et al 2009). In order to perform ABC–SMC, one requires some sufficient statistic, μ, which can be calculated from both the real data and the simulated data; a distance function; a prior distribution for each of the parameters; some distance thresholds ${\epsilon }_{i}$; and a set of perturbation kernels. The overall performance is affected by all of the factors (Filippi et al 2013, Silk et al 2013), but the choice of the summary statistic is crucial; even for good summary statistics, however, the distance function can play a major role in setting the rate and reliability of convergence. The ${\epsilon }_{i}$ are typically a series of decreasing distance thresholds where $i=0,1,\ldots ,n$ and n. Here we choose ${\epsilon }_{0}=1.0$ and decrease it adaptively so that in population i the new threshold ${\epsilon }_{i}$ is based on the 10%-ile of the distances accepted in the $(i-1)\mathrm{th}$ population. As summary statistics μ we compute the TMs; the distance function will be discussed at length later on in the paper. The prior distributions for the parameters w, p and b are all uniform distributions between 0 and 1. The ABC–SMC algorithm proceeds as follows:

We chose P = 500, $q=10\%$. The value for d is dependent on the distance function being used since they all return values of various magnitudes, but it tends towards 0.

In order to perform ABC–SMC as described above, it is mandatory to define a distance function that is able to efficiently discriminate between different random walk TM matrices. There is no standard method of computing the distance between two matrices and each proposed matrix metric focusses on distances between specific matrix characteristics. For this reason a number of potential distance function candidates were tested in order to find the best option for our application. The distances investigated are listed in table 2 (Deza and Deza 2006). The TM from real or simulated data is denoted by R and the query TM is denoted by Q, and $A=Q-R$. For the purposes of the following notation, all matrices are m x n. We chose $m=n=13$ in the case study and $m=n=15$ otherwise. These is a sufficient number of intervals to detect random walk characteristics, but also an appropriately low number of intervals to work with experimental data, as the number of data points required for reliable estimation of the TM increases dramatically with the number of intervals. The entry in the ith row and jth column of a matrix A is denoted by a_ij, while its singular values are written as s_i(A). The singular values of a matrix A are defined as the square roots of the eigenvalues of the matrix ${A}^{*}A$ where ${A}^{*}$ is the conjugate transpose of A and ${s}_{1}(A)\leqslant {s}_{2}(A)\leqslant ...$ (Deza and Deza 2006).

Table 2.  The distances that were investigated.

Distance	Formula (Deza and Deza 2006)
Frobenius norm	$| | A| {| }_{{Fr}}=\sqrt{{\displaystyle \sum }_{i=1}^{m}{\displaystyle \sum }_{j=1}^{n}| {a}_{{ij}}{| }^{2}}$
Hellinger distance	$H(R,Q)=\frac{1}{\sqrt{2}}\sqrt{{\displaystyle \sum }_{i=1}^{k}{(\sqrt{{p}_{i}}-\sqrt{{q}_{i}})}^{2}}$
Infinity norm	$| | A| {| }^{\infty }={\mathrm{max}}_{1\leqslant j\leqslant n}{\displaystyle \sum }_{j=1}^{n}| {a}_{{ij}}|$
Kullback–Leibler div. 1	${D}_{\mathrm{KL}}(R| | Q)={\displaystyle \sum }_{i}{\displaystyle \sum }_{j}\left({r}_{{ij}}\cdot \mathrm{ln}\left(\frac{{r}_{{ij}}}{{q}_{{ij}}}\right)\right)$
Kullback–Leibler div. 2	${D}_{\mathrm{KL}}(Q| | R)={\displaystyle \sum }_{i}{\displaystyle \sum }_{j}\left({q}_{{ij}}\cdot \mathrm{ln}\left(\frac{{q}_{{ij}}}{{r}_{{ij}}}\right)\right)$
One norm	$| | A| {| }^{1}={\mathrm{max}}_{1\leqslant j\leqslant n}{\displaystyle \sum }_{i=1}^{n}| {a}_{{ij}}|$
Spectral norm	$| | A| {| }_{{sp}}={({\rm{max}}\;{\rm{eigenvalue}}\;{\rm{of}}\;{\rm{}}{A}^{*}A)}^{\frac{1}{2}}$
Trace norm	$| | A| {| }_{{tr}}={\displaystyle \sum }_{i=1}^{\mathrm{min}\{m,n\}}{s}_{i}(A)$

The distances include the Frobenius norm, which is also sometimes called the Euclidean norm; the Hellinger distance, which is used to calculate how similar two probability distributions are (Duan et al 2012); two natural norms: the infinity norm and one norm, which are defined as the maximum absolute row sum and the maximum absolute column sum, respectively; and two Ky-Fan k-norms: the spectral norm and the trace norm, which are defined to be the sum of an m x n matrix' first k singular values where k = 1 and $k=$ min $\{m,n\}$, respectively. Furthermore we test the Kullback–Leibler divergence. The latter is not strictly a distance since it is asymmetrical. However, it is often used to quantify differences between distributions. For the purpose of parameter inference in an ABC scheme, the distance function does not need to be symmetric, as long as the non-negativity, the identity of indiscernibles and the triangular inequality are fulfilled. We explore both possibilities for the Kullback–Leibler divergence (${D}_{\mathrm{KL}1}(R| | Q)$ and ${D}_{\mathrm{KL}2}(Q| | R)$).

mpx:GFP/fms:RFP zebrafish embryos (Gray et al 2011) (5 days post fertilization) were anesthetized in 0.6 M MS-222 (Tricaine methanesulfonate, Sigma-Aldrich) and the tail fin was transacted using a sterile scalpel. The fish were then transferred to fresh system water for 2 h 28.5°C before transferral to 0.8% low-melt agarose (Flowgen, Lichfield, UK) for time-lapse imaging experiments. Images were captured using a Zeiss Axiovert 200 inverted microscope (Zeiss, Cambridge, UK) controlled by the C-Imaging Simple-PCI acquisition software (Hamamatsu, Sewickley, PA, USA) for up to 11 h post wounding. The temperature was maintained at 28.5°C throughout the experiment using a full incubation chamber with temperature control. The time gap between two consecutive images was 18 s. This resulted in time-lapse movies of GFP-labelled (green) neutrophils and mCherry-labelled (red) macrophages. Image processing was performed as described in (Taylor et al 2013).

In order to determine which of the seven distances is the best candidate for parameter inference we generated 500 TMs based on simulated random walk trajectories with various values for w, p and b, which we refer to as query TMs (Q). We compute each of the seven distances between these query TMs (Q) and the the five reference TMs (R) shown in table 1. The results are then represented graphically to visualize the distances' performance. The parameter values used to generate the query TMs are chosen on a lattice in parameter space with intervals 0.2 for w and 0.1 for p and b.

The first method of visualizing these results is by using an atlas of 2D heat maps. Each heat map shows the distance from (or to, in the case of ${D}_{\mathrm{KL}}(Q| | R)$) the reference TM to (or from) query TMs with varying levels of b and p, and a fixed value of w. 'Atlases' of these heatmaps are produced for each reference TM, where each row of heatmaps uses a different distance function, and each column uses query TMs produced with a different value of w. For a complete set of these heat atlases, please see supplementary figures 1–5.

Examples from the BM and BPRW1 cases with fixed w = 0.5 are shown in figure 2. A good distance will have a minimum value when both the reference and query TM have been produced using the same parameters, and will increase rapidly as the parameters become less similar. We can quantify this latter property by examining the proportion of the full parameter space that reflects a distance below a threshold relative to the minimum distance value; here, a low value would indicate a good distance measure. The results for several thresholds are shown in figures 3(A) and (B) with the Hellinger distance and the trace norm performing well.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As a second approach to visualize the results we produce rank plots of the calculated distances. This means that we group all of the distances calculated to (or from) each reference TM using a particular distance function, order these numbers from lowest to highest, and plot them on a graph. This way it is possible to see the gradient of distances as the parameters change more clearly. Again, the steeper this gradient the better, because a steep gradient means this distance function will elucidate differences between TMs whose parameters are only slightly different. The results for the trace norm and Hellinger distance can be seen in figures 3(C) and (D). For other distances, please refer to supplementary figure 6.

Visual inspection of the heat atlases clearly indicates that the Hellinger distance and the trace norm outperform the other distance functions. Both the Hellinger distance and the trace norm are minimal between TMs with identical parameters, and maximal when parameters are most different. There is also a clear, strong gradient from low to high values of the distances.

We next explore the performance of the Hellinger distance and the trace norm in an ABC–SMC inference scheme. As data we use the five simulated reference TMs, since their parameters are known. We run ABC–SMC in order to infer for each TM the original parameters of bias, persistence and weight using both, the Hellinger distance and the trace norm, as a distance between the data and simulations. We compare the inference results to the posterior distribution obtained using the Metropolis Hastings algorithm. This comparison will indicate (i) if the applied statistics (TMs) are sufficient to describe the data, and (ii) if the chosen metric is appropriate.

We approximate the posterior parameter distribution of a TM generated from a BM with parameters $(w,p,b)=(0.5,0.0,0.0)$. The marginal posterior distributions based on the Hellinger distance and the trace norm, respectively, are shown in figure 4(a). While b and p are inferred well for both distances, the marginal posterior distribution of w spans the entire prior. The reason for this becomes apparent when visualizing the posterior distribution as a 3D scatterplot. All three parameters are highly correlated in a way that a high bias can reproduce the BM characteristics, as long as the weight is supporting the persistence, which itself needs to be close to 0. Vice versa, a high persistence can reproduce the BM characteristics, as long as the weight is supporting the bias, which itself needs to be close to 0. In order to reduce this dependency, we rescale the parameters by simple multiplication of the parameters. Since the weight is equivalent to the probability of the particle choosing the bias distribution over the persistence distribution, $b^{\prime} ={wb}$ and $p^{\prime} =(1-w)p$ where b' and p' are the rescaled bias and rescaled persistence parameters, respectively. The rescaled parameters (p',b') used in generating the BM are (0.0, 0.0). The rescaled posterior distribution matches these values well for both Hellinger distance and the trace norm.

Zoom In Zoom Out Reset image size

We repeat the inference scheme with the remaining four TMs. The results for the BPRW1 are shown in figure 4(b). The true parameters $(w,p,b)\;=(0.5,0.7,0.3)$ are inferred well by both distances. However, the marginal posterior distributions based on the Hellinger distance are narrower around the true parameter values, which shows a better performance of this distance. Note, the 3D scatterplot does not show any correlations between the three parameters, which indicates that no other parameter combinations can reproduce the characteristics of the BPRW model. Therefore it is not a surprise that also the rescaled parameters $(p^{\prime} ,b^{\prime} )=(0.35,0.15)$ are inferred well.

For all five tested scenarios the approximated posterior distributions are in good agreements with the exact posterior distributions obtained by the Metropolis Hastings algorithm. The main deviation is observed for the weight parameter w in the Brownian motion case 4(a). Furthermore, the exact marginal posterior distributions are slightly narrower than the approximated counterparts, which indicates some loss of information in the TMs and/or the applied metric.

In conclusion, the Hellinger distance is the best distance function to infer bias and persistence parameters for 2D random walks using TMs as summary statistics. The inference results are in agreement with exact methods where these are applicable.

To test the Hellinger distance and evaluate its application on real-world problems we apply this distance to infer the random walk parameters of migrating immune cells. More specifically, we extract macrophage and neutrophil tracks from living zebrafish embryos that had undergone tail transection. Macrophages and neutrophils are the first layer of defence during inflammation, here mimicked by wounding. The extracted cell tracks were grouped according to how long after wounding and how far away from the wound they were detected. This results into spatial-temporal clusters, consisting of two temporal groups (T1: 1–6 h; T2: 6–11 h), each of them consisting of seven spatial groups that are generated by a sliding window approach. Figure 5 provides an overview of how the data are separated. For all 28 groups of trajectories (14 for macrophages and 14 for neutrophils) we compute the TM from the extracted cell tracks and infer the random walk parameters (figures 5 and 6). The ABC–SMC framework using the Hellinger distance successfully reproduces the TM computed from the in vivo data (figure 5(B)). We compare the distances from the data TM to the simulated TM to distances from the same data TMs to a given realization of a random walk with $p=b=0$. As shown in supplementary figure 8, the former distances are significantly smaller for every TM, as is expected.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The posterior distributions of the inferred rescaled parameters $p^{\prime} =(1-w)p$ and $b^{\prime} ={wb}$ are given in figure 6. It becomes apparent that the level of observed bias and persistence during the first 6 h after wounding is dependent on the distance from the wound for neutrohils. Macrophages, on the contrary, have a nearly constant level of persistence and only a modest increase of the bias with increasing distance from the wound. After 6 h the observed bias is decreased close to the wound, while the level of observed persistence is increased. Neutrophils show a comparable behaviour before and after 6 h. In general, we find that while the level of bias is similar in both cell types, the level of persistence is significantly higher in neutrophils than in macrophages. This case study demonstrates the potential of TMs as summary statistics combined with the Hellinger distance in an ABC framework.