Modeling three-dimensional T-cell motility using clustering and hidden Markov models

Abstract

Recent advances in imaging technologies now allow for real-time tracking of fast-moving immune cells as they search for targets such as pathogens and tumor cells through complex three-dimensional tissues. Cytotoxic T cells are specialized immune cells that continually scan tissues for such targets to engage and kill, and have emerged as the principle mediators of breakthrough immunotherapies against cancers. Modeling the way these T cells move is of great value in furthering our understanding of their collective search efficiency. T-cell motility is characterized by heterogeneity at two levels: (a) Individual cells display different distributions of translational speeds and turning angles, and (b) each cell can during a given track, its motility, switch between local search and directional motion. Despite a likely considerable influence on a motile population’s search performance, statistical models that accurately capture both such heterogeneities in a distinguishing manner are lacking. Here, we model three-dimensional T-cell trajectories through a spherical representation of their incremental steps and compare model outputs to real-world motility data from primary T cells navigating physiological environments. T cells in a population are clustered based on their directional persistence and characteristic “step lengths” therein capturing between-cell heterogeneity. The motility dynamics of cells within each cluster are individually modeled through hidden Markov model to capture within-cell transitions between local and more extensive search patterns. We explore the importance of explicitly capturing altered motility patterns when cells lie in close proximity to one another, through a non-homogenous hidden Markov model.

Keywords

Three-dimensional spatio-temporal trajectories clustering heavy tail distribution homogenous hidden Markov models andnon-homogenous hidden Markov models persistency T cells trajectory metrics

Introduction

Our immune systems are vital to our physiological wellbeing, protecting us from pathogens and eradicating cancerous cells, for example.^1–4 The immune response emerges as a system-level phenomenon mediated by local interactions among numerous cell types. Many immune cells are highly motile, traversing barriers and tissues of varying biochemical and mechanical compositions, “scanning” other cells and their environment, and monitoring for patterns and signals that indicate any damages, infection, or cancer, for example.⁵ The delivery of immune cell function is dependent on proximity to other cells, and motility is therefore a critical process in effecting an immune response. Advances in imaging technology have enabled the live three-dimensional (3D) tracking of fast-moving cells, such as T cells, as they traverse their environment and interact with one another; see the work of Galeano et al.⁶ This has spurred the detailed study of both individual and collective immune cell motility and the system-level impacts of local-level processes. Through in silico reproduction, we can model cellular motility, and consider how manipulations will impact the wider immune response. The field is still in its infancy, with the ongoing discourse on how best modeling practice and which metrics fully capture and specify a process as complex as population-wide 3D motility. The emerging picture is one where cells within a population are heterogeneous in their motility characteristics, and that motility at a population level differs between tissues and pathological contexts.

A number of competing models for T-cell motility have been developed. Miller et al.¹ and Textor et al.⁷ consider random walk models for lymphocytes observed in living tissue. In the latter work it was concluded using Hotelling’s test that the displacement of T cells over time is consistent with random motions. These models can be framed more generally within populations of T cells moving according to advection/diffusion equations, which can be altered to model T-cell motions with heterogeneous tracks, such as T cells slowing down and rapidly changing direction for several time steps; see, for example, the literary works of Textor et al.,⁷ Schienbein and Gruler,⁸ Niculescu et al.,⁹ Manolopoulou et al.,¹⁰ Beauchemin et al.,¹¹ and Beauchemin et al.¹²

Competing motility models are often evaluated by comparing various summary statistics describing T-cell trajectories, termed motility metrics, due to the complexity of 3D spatio-temporal trajectories. These include, among many others, asphericity, straightness, outreach ratio, displacement ratio, Hurst index, mean squared displacement, and root mean squared displacement, see the literary works of Banigan et al.,¹³ Mokhtari et al.,¹⁴ each of which we discuss below. For instance, Banigan et al.¹³ obtained the root mean squared displacement and normalized correlation between the displacement from in-vitro lymphocyte data and compared them with those obtained from their proposed model. Their model, a mixture model, also suggests heterogeneities between T-cell displacement.

Although the above-referenced models are sometimes able to reproduce trajectories with similar motility metrics as the observed trajectories, there are often sizable differences. One potential explanation for this is that T cells are naturally heterogenous, i.e. there appear to be different types of movement between T-cell trajectories. Some T cells have predominantly persistent motion and move rapidly through a given area of space while others have less persistent motion, and tend to linger locally in a given area of space for long periods of time. They also exhibit lower or higher speeds, i.e. short or long step lengths, over a fixed time; see the work of Mokhtari et al.¹⁴ This raises the potential of using clustering as a step to first categorize the data to homogenous groups and accordingly fit a cluster-specific statistical model that explains the underlying motility within each group.

There exists an established body of work on clustering general trajectories. Ester et al.¹⁵ clustered two-dimensional (2D) or 3D trajectories using the density defined on the neighborhood around each index of a trajectory. Kalnis et al.¹⁶ considered moving trajectories; the proposed methodology was based on clustering points at a fixed time. Lee et al.¹⁷ considered the parallel, perpendicular, and angular distances between subtrajectories for the purpose of clustering. Jebara et al.¹⁸ clustered trajectories by fitting homogenous hidden Markov models (HHMMs) to each trajectory and then measuring the distance between the fitted models. Zhang et al.¹⁹ introduced a new clustering scheme of spatio-temporal trajectories based on kernel density estimation. Others, see the literary works of Ester et al.,²⁰ Wang et al.,²¹ Wu et al.,²² and Pan et al.,²³ discussed clustering of trajectories in alternative settings.

In the context of lymphocyte motility, Mokhtari et al.¹⁴ clustered T cells’ trajectories using a hierarchical clustering method. A combination of staggered measures, i.e. volume asphericity, volume prolateness, confinement ratio, displacement ratio, and outreach ratio, were used to cluster trajectories into three subpopulations: Fairly straight tracks, strongly confined tracks, and purely random tracks. Notably, step lengths were not considered in the analysis as a tool for clustering. An alternative approach for clustering of 3D spatio-temporal trajectories is to distinguish Brownian type motion from Levy type motion; see the work of Plank et al.,²⁴ as well as a persistent versus non-persistent motion which we present in this paper.

Individual T cells also appear to exhibit a mixture of movements within a given area of space. Stop and go models, proposed in Dustin et al.,²⁵ Beauchemin et al.,¹² Textor et al.,¹¹ are partial differential equation models which attempt to explain this heterogeneity of movements, while Méhes and Vicsek²⁶ used Bayesian analysis of autoregressive processes of the first order to study the change of parameters observed in a given area of space, similar to Metzner et al.²⁷ where they have the parameters updated at each step for the vectorial displacements in 2D. Banigan et al.¹³ considered a mixture model of the aggregated increments under the assumption of spatial isotropy, where their model is a mixture of populations for the trimmed data. Read et al.²⁸ considered random effect models for individual T-cell trajectories. Their results indicate a correlated heterogeneous random walk where the turn speed and translational speed are negatively related while the model provides the best fit in terms of Kolmogorov–Smirnov statistic. Bartolucci and Farcomeni²⁹ considered computer simulations of a random walk process to model the movement in 2D.

An approach to modeling heterogeneity within T-cell trajectories might be to use hidden Markov models (HMMs). Raffa and Dubin³⁰ proposed to model unobserved longitudinal smoking behavior in a smoking cessation clinical trial using HMM’s. HMM’s have been used in the past to model the latent switching behavior observed in animal movements; see the literary works of Patterson et al.,³¹ Langrock³² and Patterson et al.³³ Importantly, these models to date are formulated for 2D motion. To the best of our knowledge, HMM’s have not yet been considered for populations of agents searching in 3D.

This problem shares some similarities to that of modeling the trajectories of marine animals moving in three dimensions, which has been conducted in several cases using HMM’s. Pedersen et al.³⁴ applied a direct Fokker–Planck-based methodology using HMMs to data of fish positions in 2D. Thygesen et al.³⁵ tried to accommodate the measurement error issue for the 2D movement also for fish motility. Vermard et al.³⁶ in a research project considered a Bayesian Hierarchical Model using hidden Markov process in 2D for the three behavioral states of boat activity, specifically stopping, steaming, and fishing; see also Mitani et al.,³⁷ Liu et al.,³⁸ Laplanche et al.,³⁹ Orellana et al.,⁴⁰ Sims et al.,⁴¹ and Humphries et al.⁴²

The idea of HHMM and non-homogenous hidden Markov models (NHMM’s) statistically have been also studied for forecasting univariate time series, among many, see the literary works of Spezia,⁴³ Koki et al.,⁴⁴ Meligkotsidou and Dellaportas,⁴⁵ Holsclaw et al.,⁴⁶ and Maruotti and Rocci.⁴⁷ Spezia⁴³ analyzed the natural logarithm of the daily maximum hourly mean concentrations of ozone in Venice. Koki et al.⁴⁴ considered analyzing the natural logarithm of the realized volatility series; see the literary works of Meligkotsidou and Dellaportas,⁴⁵ Holsclaw et al.,⁴⁶ and Maruotti and Rocci⁴⁷ for other examples. These works statistically address the difficulty raised by the connection between HHMM or NHMM and generalized linear, or linear, regression models.

In this paper, we model T-cell trajectories using HMM’s after a preliminary clustering step, which accounts for the observable features of these trajectories. In particular, we consider one-state, two-state, and three-state homogeneous HMM’s. Biologically, these states pertain to persistent and non-persistent (or “go” and “stop”) states, with a middle persistent behavior in the three-state case.

Another factor of interest in lymphocyte motility is the effect that neigboring proximal T cells have on one another’s motility since there is evidence that T cells interact with each other. Other work has focused on the effect of chemokines on T-cell motility and conceptualizing cell trajectories as features of environmental obstacles, see the work of Beltman et al.⁴⁸ The possibility of calibrating the configuration of an environment by proxy of the resultant cellular motility is intriguing, see the work of Read et al.²⁸ In the literature, Miller et al.¹ considered T cells interacting with one another in terms of abrupt changes in shape, velocity, and direction. Méhes and Vicsek²⁶ also explained the influence of other T cells on the behavior of a specific T cell by defining collective T-cell motion. To the best of our knowledge, our proposed model is unique in explicitly capturing/allowing for this phenomenon.

Here, to measure the effect of T cells in close proximity to each other, we consider NHMM’s, in which the transition matrices depend on the information about the existence of at least one T cell in the neighborhood of a given (index) T cell.

The paper is organized as follows. We first discuss the experiments that produced the T-cell motility data we analyze in the ‘Observing real T-cell motility through live 3D imaging section. In the ‘Methods section, we present an innovative methodology for clustering spatio-temporal 3D T-cell trajectories. We introduce new measures of trajectories, relative persistence (RP) and compactness measure (CM), to differentiate trajectories based on their persistency and underlying distributions of their increments. We then explain the implementation of the HHMM and NHMM for the angular change and step length. The proposed methodology is implemented in the ‘Data analysis section to analyze T-cell trajectories from the experiment described in the ‘Observing real T-cell motility through live 3D imaging section. In the ‘Simulation studies section, we further investigate and evaluate the proposed models with a comprehensive simulation study. Conclusion and future research directions are presented in the ‘Discussion section.

Observing real T-cell motility through live 3D imaging

In this paper, we focus on the following experiment. Primary cytotoxic (CD8 $^{+}$ ) T cells derived from the OT-I transgenic mouse model were embedded within a 3D collagen matrix and then imaged through confocal microscopy every 20 seconds for 1 hour; see the work of Galeano et al.⁶ for methodological details. Spatio-temporal 3D trajectories of T-cell motility within a finite imaging volume of dimension $390 \times 390 \times 65 μ$ m, the observation space, were obtained following image analysis with the Imaris software (Bitplane AG, Zurich, Switzerland). Figure 1 shows 110 real T cells trajectories from a representative observation space.

This methodology allows for the tracking of individual cell trajectories. Trajectories of less than 10 minutes duration were omitted. T cells are free to enter and leave the imaging volume in any direction. This experiment was repeated independently five times, collecting a total of 47,235 3D cell positions distributed across 764 trajectories. Given that the experiments are assumed to be homogenous, we artificially divided the data into a training dataset, experiments 1, 3, and 5, and an out-of-sample dataset, experiments 2 and 4, such that there were nearly equal numbers of T cells in each dataset. We have a random selection for the out-of-sample data and training data while the other choice might be finding all possible selections for the training data, say $(\binom{5}{3})$ , and fitting the models separately. Then, averaging over estimates of model parameters in different set-ups gives the estimate of model parameters that is beyond the scope of this work. In this analysis, we keep all T cells inside an experiment together instead of having them randomly chosen for the training dataset to discover how T cells interact with each other in proximity with neighboring T cells. The training data contained 352 T-cell trajectories with 24,173 total recorded 3D positions while the out-of-sample dataset contained 412 T-cell trajectories with 23,062 total recorded positions.

Cells are not tracked outside of the observation space; a cell’s re-entry will register as a new track. We investigated the potential experimental bias of such re-entries via simulation studies below and, in short, we determine that this issue does not influence our conclusions.

The data is available per request and also through the supplementary materials.

Methods

Suppose $T_{i} = {(x_{i, j}, y_{i, j}, z_{i, j}); 1 \leq j \leq n_{i}}$ constitutes the Cartesian representation of the $i$ ’th T cell’s trajectory where $n_{i}$ is the number of observations (positions) per trajectory, for $i = 1, \dots, N$ , where $N$ is the total number of trajectories observed. Also, let

I_{i, j} = (Δ x_{i, j}, Δ y_{i, j}, Δ z_{i, j}) = (x_{i, j + 1} - x_{i, j}, y_{i, j + 1} - y_{i, j}, z_{i, j + 1} - z_{i, j})

(1)

where

Δ x_{i, j}

Δ y_{i, j}

, and

Δ z_{i, j}

are the jth increment of trajectory

i

, for

j = 1, \dots, n_{i} - 1

, in

X

Y

, and

Z

directions. Quantifying the directionality and speed of motile cells is more easily performed by considering a spherical representation of their increments, which we now present. The idea is that some measures of T-cell movement can be better defined using the spherical coordinates compared to the Cartesian coordinates. The angle between two subsequent increments is defined as

ψ_{i, j} = \arccos (\frac{⟨ I_{i, j}, I_{i, j + 1} ⟩}{‖ I_{i, j} ‖ ‖ I_{i, j + 1} ‖})

(2)

where

⟨ q, s ⟩ = q_{1} s_{1} + q_{2} s_{2} + q_{3} s_{3}

, the inner product in

R^{3}

‖ s ‖ = \sqrt{s_{1}^{2} + s_{2}^{2} + s_{3}^{2}}

“step length”, we use

r_{i, j} = ‖ I_{i, j} ‖

to simplify the notation, and

0 \leq ψ_{i, j} \leq π

, for

j = 1, \dots, n_{i} - 2

q = (q_{1}, q_{2}, q_{3})

, and

s = (s_{1}, s_{2}, s_{3})

. Using (2),

n_{i} - 2

angles are obtained for each trajectory. The

ψ_{i, j}

’s are obtained using the inner product and cosine of the angle between increments, which shows the percentage of

j + 1

th increment,

I_{i, j + 1}

, in the same direction as the

j

th one. This provides an approach to differentiate T cells based on their persistence.

The increment $I_{i, j}$ is fully specified in spherical coordinates by the set $(ψ_{i, j}, ‖ I_{i, j} ‖, θ_{i, j}),$ where $θ_{i, j}$ is the radial angle between $I_{i, j}$ and $I_{i, j - 1}$ , which we now define. Let $u_{i, j - 1}$ , and $v_{i, j - 1}$ be the unit vectors obtained by completing $I_{i, j - 1} / ‖ I_{i, j - 1} ‖$ , $(1, 0, 0)^{⊤}$ and $(0, 1, 0)^{⊤}$ to an orthonormal basis of $R^{3}$ using the Gram–Schmidt procedure. The vectors $u_{i, j - 1}$ and $v_{i, j - 1}$ are well defined as long as $I_{i, j - 1}$ does not lie in the plane spanned by $(1, 0, 0)^{⊤}$ and $(0, 1, 0)^{⊤}$ . In the unlikely case that it is, we replace $(1, 0, 0)^{⊤}$ , or $(0, 1, 0)^{⊤}$ if $I_{i, j - 1} / ‖ I_{i, j - 1} ‖ = (1, 0, 0)^{⊤}$ , with $(0, 0, 1)^{⊤}$ and perform the same procedure. We then define the radial angle as the angle in standard polar coordinates of the point $(⟨ I_{i, j}, u_{i, j - 1} ⟩, ⟨ I_{i, j}, v_{i, j - 1} ⟩)$ in $R^{2}$ , which is given by

θ_{i, j} = \arctan 2 (⟨ I_{i, j}, v_{i, j - 1} ⟩, ⟨ I_{i, j}, u_{i, j - 1} ⟩)

where

\arctan 2

is the angle in radians between the positive

x

-axis of a plane and the point given by

(⟨ I_{i, j}, u_{i, j - 1} ⟩, ⟨ I_{i, j}, v_{i, j - 1} ⟩)

, for

j = 2, \dots, n_{i} - 1

. In other words, the radial angle is the polar angle obtained by projecting

I_{i, j}

onto the plane perpendicular to

I_{i, j - 1}

. In contrast,

ψ_{i, j}

is the angle between the vectors

I_{i, j}

and

I_{i, j - 1}

. We do not use the assumption of isotropy, following the same distribution in

X

Y

, or

Z

direction, in our work.

Clustering spatio-temporal 3D trajectories to capture between cell heterogeneity

T-cell trajectories differ in terms of step length and also the persistency of their movement. In order to take into account these differences and improve our understanding of T cells’ movements, we propose clustering T cells based on measures of RP and track compactness. The methodology we propose for clustering 3D spatio-temporal trajectories is independent of their direction or distance from other moving T cells nearby. We propose below a new physically interpretable measure of persistency in order to categorize trajectories with mostly persistent behavior. We define RP of the $i$ ’th trajectory as

{RP}_{i} = 100 \times [\frac{1}{n_{i} - 2} \sum_{j = 1}^{n_{i} - 2} I_{ψ_{0}} (ψ_{i, j})]

(3)

where

I_{ψ_{0}} (ψ) = 1

ψ \leq ψ_{0}

, and zero otherwise.

{RP}_{i}

then calculates the percentage of increments for which the

i

’th T cell exhibits an angular change less than

ψ_{0}

. Intuitively, a T cell’s trajectory exhibits persistence if it does not frequently turn back on its path. Textor et al.⁷ defined cells with a “short time persistency” as having average turning angles lower than

π / 2

; we adopt the same threshold here; i.e.

ψ_{0} = π / 2

. RP then calculates the percentage of movements a moving object does not turn back on its own path; i.e. RP obtains the rate of persistency for a trajectory inside the viewing region.

In addition to the persistence, we also consider clustering trajectories based on the tail properties of their increment displacement distributions. Considerable attention has been devoted to differentiating between T cells undergoing Brownian motion and/or Lévy flight; see the literary works of Banigan et al.¹³ and Mokhtari et al.¹⁴ Under 3D Brownian motion $X$ , $Y$ and $Z$ increments, $Δ X$ , $Δ Y$ , and $Δ Z$ , are normally distributed, while for Lévy flight movement $X$ , $Y$ , and $Z$ increment displacement distributions are long-tailed. In general, estimating the tail properties of a distribution is difficult, especially without large samples; see the work of Huisman et al.⁴⁹ for a discussion. We propose here to differentiate T cells based on a general measure of the compactness of their trajectories that is sensitive to heavy versus light-tailed increments. In particular, we introduce the CM for the increments in the $X$ direction

{CM}_{i}^{x} = \frac{Q_{3 i x} - Q_{1 i x}}{max (X X_{i}) - min (X X_{i})}

(4)

where

{X X}_{i} = (Δ x_{i, 1}, \dots, Δ x_{i, n_{i} - 2})

Q_{1 i x}

Q_{3 i x}

max (X X_{i})

, and

min (X X_{i})

are the observed increments for the

i

’th trajectory in

X

direction, first quartile, third quartile, maximum and minimum of

{X X}_{i}

, respectively. We use subscript

x

to show

{CM}_{i}^{x}

is obtained on

X

axis.

{CM}_{i}^{y}

and

{CM}_{i}^{z}

are similarly defined for increments in

Y

and

Z

directions. The assumption of isotropy is removed in the sense that there is no need for trajectories to behave the same on different axes. Suppose one generates 100 increments which is a typical number of cell observations in a live-imaging experiment, from a heavy-tailed, for instance, Cauchy distribution. In this case, one would expect only a few increments to register as outliers. This is why we propose to use the range of the increments in the definition of the CM. Having the range in the denominator of the CM helps distinguish light-tailed from heavy-tailed distributions by bringing the CM to the same scale.

We design a classifier based on prediction intervals (PI’s) to compare the underlying distribution of observed increments in $X$ , $Y$ , and $Z$ directions with the normal distribution. To this end, we obtain PI’s of ${CM}_{i}^{x}$ , ${CM}_{i}^{y}$ , and ${CM}_{i}^{z}$ under the assumption of normality using a parametric bootstrap, see the work of Efron and Tibshirani,⁵⁰ as follows.

Based on the data, the sample means and variances are calculated, i.e. ${\bar{X}}_{i}$ , $s_{x_{i}}^{2}$ , ${\bar{Y}}_{i}$ , $s_{y_{i}}^{2}$ , and ${\bar{Z}}_{i}$ , $s_{z_{i}}^{2}$ , sample mean and variance of increments in $X$ , $Y$ , and $Z$ directions, respectively, see Equation (1) for definition. Then, conditionally on the data, assume $Δ X_{i, j}^{b} \sim N ({\bar{X}}_{i}, s_{x_{i}}^{2})$ , $Δ Y_{i, j}^{b} \sim N ({\bar{Y}}_{i}, s_{y_{i}}^{2})$ , and $Δ Z_{i, j}^{b} \sim N ({\bar{Z}}_{i}, s_{z_{i}}^{2})$ are independently identically distributed random variables, for $j = 1, \dots, n_{i} - 1$ , $b = 1, \dots, B$ , and $B = 10, 000$ . Also, under the Brownian motion assumption, suppose increments in $X$ , $Y$ , and $Z$ directions are independent of each other. Define $A_{i}$ as the classifier for the $i$ ’th trajectory. Let $A_{i} = 0$ , for the short-increment group, as would be generated under Brownian motion. Initially, $n_{i} - 1$ sample units are generated in $X$ , $Y$ , and $Z$ directions. The CM for $B$ sets of increments in different directions are calculated and the 95 $%$ PI’s are constructed separately in each direction. Small values of ${CM}_{i}^{x}, {CM}_{i}^{y}$ , or ${CM}_{i}^{z}$ indicate heavier tails and departure from normal increments. If any of the observed values of ${CM}_{i}^{x}$ , ${CM}_{i}^{y}$ , or ${CM}_{i}^{z}$ in a trajectory is less than the lower bound of their corresponding intervals, the trajectory is classified as having long increments indicating the trajectory exhibits Lévy flight type motion, $A_{i} = 1$ . It is worth mentioning that the values of ${CM}_{i}^{x}$ , ${CM}_{i}^{y}$ , or ${CM}_{i}^{z}$ and ${RP}_{i}$ ’s might be obtained from T cells with different level of accuracy due to their inherent difference in the length of stay in the viewing region, in practice, our data is unbalanced.

Spatio-temporal 3D trajectories are clustered considering ${RP}_{i}$ ’s and ${CM}_{i}$ ’s. Table 1 shows our general rule of the proposed clustering scheme where $q_{l}$ and $q_{b}$ are coming from k-means clustering; see the supplemental information for the example. Given the data, there is a need for introducing more clusters. In our application below, we take $q_{l} = q_{b} = 87 %$ , which is the median of ${RP}_{i}$ from the training data. The choice of $q_{l}$ and $q_{b}$ is due to our limitation with the number of trajectories. While the introduction of a clustering scheme increases the accuracy of model fit because of the introduction of mixture models into the data analysis, having more detailed clusters is the scope of our future research; see the literary works of Maruotti and Rocci,⁴⁷ and Hastie et al.⁵¹

Table 1.

Our general rule of the proposed clustering scheme for trajectories. HP, NP, LI, and SI refer to highly persistent, non-persistent, long increments, and short increments, respectively.

Cluster	LI	SI
HP	${RP}_{i} \geq q_{l}$ and $A_{i} = 1$	${RP}_{i} \geq q_{b}$ and $A_{i} = 0$
NP	${RP}_{i} < q_{l}$ and $A_{i} = 1$	${RP}_{i} < q_{b}$ and $A_{i} = 0$

HP: highly persistent; NP: non-persistent; RP: relative persistence; LI: long increment; SI: short increment

Modeling angular change, step lengths, and radial angle

Given the initial clustering performed as described above, the goal is now to model the spherical coordinates, step length, angular change, and also the radial angle, of trajectories belonging to each cluster. If $C$ clusters are produced, let the cluster membership for $i$ th trajectory, $i = 1, \dots, N$ , be defined by $c_{i} \in {1, \dots, C}$ . The first step length in each trajectory is removed from the analysis due to the unavailability of its corresponding angles.

We model the behavior of trajectories in each cluster by NHMM’s and HHMM’s. In the first step, we define statistical models on spherical coordinates to find the corresponding likelihood function. To this end, $M$ hidden states for $ψ_{i, j}$ , angular change, and, accordingly, $\log (r_{i, j})$ , the logarithm of step length, are considered. The hidden state is defined based on the behavior of angular change. In the case of two-state models, the hidden states can be generally interpreted as persistent and non-persistent motion, while the states in the three-state models can be interpreted as persistent, moderately persistent, and non-persistent motion. Let $f_{B_{m}}^{c_{i}}$ denote the density of the Beta $(α_{m}^{c_{i}}, β_{m}^{c_{i}})$ distribution where $α_{m}^{c_{i}}$ and $β_{m}^{c_{i}}$ are unknown, for $m = 1, \dots, M$ and $c_{i} \in {1, \dots, C}$ . The hidden state at time $j = 2, \dots, n_{i} - 1$ for $i$ ’th trajectory is denoted by $H_{i, j}^{c_{i}}$ . Given $H_{i, j}^{c_{i}} = m$ , we assume

η_{i, j} = \frac{ψ_{i, j}}{π} \sim f_{B_{m}}^{c_{i}}

(5)

for

j = 1, \dots, n_{i} - 1

and

i = 1, \dots, N_{c}

The beta distribution is used due to its flexibility for modeling observations taking values in the $(0, 1)$ interval. The initial probabilities, ${π π}^{c_{i}} = [π_{1}^{c_{i}}, \dots, π_{M}^{c_{i}}]$ , are assumed to be the same for trajectories in the same cluster where $π_{m}^{c_{i}} = P (H_{i, 1}^{c_{i}} = m)$ , for $m = 1, \dots, M$ . Following Patterson et al.³³ and Schienbein and Gruler,⁸ we assume $ψ_{i, j}$ , $θ_{i, j}$ , and $r_{i, j}$ are independent of each other given the $H_{i, j}^{c_{i}}$ , however, $η_{i, j}$ depends on the previous movement, $η_{i, j - 1}$ , by the Markovian property, for $j = 2, \dots, n_{i} - 1$ and $i = 1, \dots, N^{c}$ . We leave the more advanced case of correlated spherical coordinates for our future work. Figure 2 shows the distribution of the radial angle obtained from 352 trajectories in the training dataset, from which it is clear that $θ θ$ is roughly uniformly distributed, which suggests T cells scan their environment isotropically. Also, the clustering scheme introduced in the ‘Clustering spatio-temporal 3D trajectories to capture between cell heterogeneity section is independent of the direction which supports taking $θ_{i j} \sim Uniform (0, 2 π)$ in each cluster. We give more evidence for this in the supplemental information, Torkashvand,⁵² and in the ‘Data analysis section this assumption is investigated further.

Let ${G G}_{c, i}$ be the matrix containing densities of the mixture components of the HHMM/NHMM

{G G}_{c, i} (η_{i, j}; {α α}^{c_{i}}, {β β}^{c_{i}}) = d i a g [f_{B_{1}}^{c_{i}} (η_{i, j}), \dots, f_{B_{M}}^{c_{i}} (η_{i, j})]

(6)

i.e.

{G G}_{c, i} (.; {α α}^{c_{i}}, {β β}^{c_{i}})

is a diagonal matrix which defines the distribution of the scaled angular change,

η_{i, j}

’s, in different states,

m = 1, \dots, M

, for the

i

th trajectory of the

c

th cluster at the

j

th time point,

j = 2, \dots, n_{i} - 1

. The unknown parameters are

{α α}^{c_{i}} = (α_{1}^{c_{i}}, \dots, α_{M}^{c_{i}})

and

{β β}^{c_{i}} = (β_{1}^{c_{i}}, \dots, β_{M}^{c_{i}})

. Let

{a a}_{i, 1} = {π π}^{c_{i}} {P P}_{c_{i}}^{(u_{i 1}^{c_{i}})}

and

{a a}_{i, j} = {a a}_{i, j - 1} {P P}_{c_{i}}^{(u_{i, j}^{c_{i}})} {G G}_{c, i} (η_{i, j}; {α α}^{c_{i}}, {β β}^{c_{i}})

be a vector of probabilities in different states starting from the first location of the trajectory under the Markovian assumption, for

j = 2, \dots, n_{i} - 1

and

i = 1, \dots, N_{c}

. The likelihood function is then given by

L (Γ^{c}) = \prod_{i = 1}^{N_{c}} [{(\frac{1}{2 π})}^{n_{i} - 2} {a a}_{i, n_{i} - 1} 1 1]

(7)

where

Γ^{c} = (α_{1}^{c}, β_{1}^{c}, \dots, α_{M}^{c}, β_{M}^{c}, π_{1}^{c}, \dots, π_{M}^{c}, {P P}_{c})

and

{P P}_{c}

contains the possible transition matrices,

P_{c}^{(1)}, P_{c}^{(2)}, \dots, P_{c}^{(K)}

, considered in the analysis potentially based on some known categorical covariate

u_{i, j}^{c_{i}} = 1, \dots, K

, for

j = 1, \dots, n_{i} - 1

i = 1, \dots, N^{c}

c = 1, \dots, C

, and

1 1 = [1, \dots, 1]_{1 \times M}^{'}

. The value of

{a a}_{i, n_{i} - 1}

gives the hidden Markov probability at different states obtained from the

i

’th trajectory in the last step. The covariate

u_{i, j}^{c_{i}}

is a coding covariate defined based on the context. It is worth mentioning that

u_{i, j}^{c_{i}}

is utilized when we deal with the non-homogenous HMM, it then used to select the transition matrix,

{P P}^{(u_{i j}^{c_{i}})}

For the T-cell trajectories, the interest is to see the effect of other T cells on an index T cell in close proximity. To this end, we consider two separate scenarios. In the HHMM, we assume that the T cells’ behavior are independent of each other where only a single transition matrix is considered so that $K = 1$ . For the NHMM, we study the possibility of interaction between T cells; see the work of Patterson et al.³¹ If there exist no T cells in a neighborhood of radius $ρ$ , we assume T cells change their type of movement following the transition matrix of $P_{c_{i}}^{(1)}$ , i.e. $u_{i, j}^{c_{i}} = 1$ . Otherwise, the transition matrix will be $P_{c_{i}}^{(2)}$ , i.e. $u_{i, j}^{c_{i}} = 2$ , so that $K = 2$ . The behavior of T cells might get affected because of having a different number of T cells which is beyond the scope of this work and we leave it for future research. We performed some preliminary analysis to choose the neighborhood radius $ρ$ empirically and took $ρ = 22 μ m$ to balance between too many or too few T cells in the neighborhood. We note that two T cells within a neighborhood of radius $ρ$ are close to contacting each other. The value of $ρ$ is consistent with observations of T-cell interactions in biological experiments.

Due to the complicated form of the likelihood in (7), we apply a Bayesian approach using Hamiltonian Monte Carlo, Neal et al.,⁵³ to estimate the model parameters. Following Berger et al.,⁵⁴ Torkashvand et al.,⁵⁵ and Schad et al.,⁵⁶ we select informative prior distributions. It is worth mentioning that initially, we performed the analysis implementing non-informative prior distributions, but the parameter estimates did not converge.

Informative uniform prior distributions are selected on $α_{m}^{c}$ and $β_{m}^{c}$ , for $m = 1, \dots, M$ and $c = 1, \dots, C$ ; e.g. $α_{1}^{c} \sim Uniform (0, 10)$ and $β_{1}^{c} \sim Uniform (0, 4)$ . Dirichlet prior distributions are defined with the data-driven parameters for the initial probabilities and also rows of the transition matrices. In the ‘Data analysis section, we describe the prior distributions imposed on the model parameters of our example of T-cell trajectories in detail.

We now consider modeling the step lengths corresponding to each trajectory. Firstly, we apply a log transformation to the step length and consider models for $R_{i, j} = \log (r_{i, j})$ , since these are approximately normally distributed; see Figure 6. We also considered modeling the step length using a gamma distribution; see the supplemental information. These variables are then modeled in groups based on local decoding of the NHMM and HHMM models for $η_{i, j}$ ,⁵⁷ p. 87 for $j = 2, \dots, n_{i} - 1$ and $i = 1, \dots, N$ , as follows. Let

{b b}_{i j} (Γ^{c_{i}}) = [b_{i j, 1}, \dots, b_{i j, M}] = {π π}^{c_{i}} [\prod_{k = 1}^{j} {P P}_{c_{i}}^{(u_{i k})}] {G G}_{c_{i}} (η_{i, j})

where

{b b}_{i j} (Γ^{c_{i}})

is a vector of mixture components under the Markovian assumption, for

j = 2, \dots, n_{i} - 1

and

c_{i} \in {1, \dots, C}

. Consequently, the probability states given the observed angular change are

P (H_{i, j}^{c_{i}} = m | η_{i, j}) = \frac{b_{i j, m}}{{b b}_{i j} 1 1}

for

m = 1, \dots, M

. A T cell is considered in the

m^{*}

state if

P (H_{i, j}^{c_{i}} = m^{*} | η_{i, j}) = max_{1 \leq m \leq M} P (H_{i, j}^{c_{i}} = m | η_{i, j})

(8)

The logarithm of step length, $R_{i, j}$ , is assigned to the particular state $m$ using (8), for $j = 2, \dots, n_{i} - 1$ , $i = 1, \dots, N_{c}$ , and $c_{i} \in {1, \dots, C}$ . In this paper, we model ${R R}_{m}^{c} = {R_{i, j} | H_{i, j}^{c_{i}} = m, j = 2, \dots, n_{i} - 1, i = 1, \dots, N_{c}}$ , the logarithm of step length at state $m$ in cluster $c$ , with the normal distribution, $N (μ_{m}^{c}, σ_{m}^{c})$ , for $m = 1, \dots, M$ and $c = 1, \dots, C$ . The fitting of normal distributions to ${R R}_{m}^{c}$ facilitates comparisons and interpretation across clusters and hidden states, for $m = 1, \dots, M$ and $c = 1, \dots, C$ . For purpose of comparison, we use the gamma distribution with similar flexibilities to model $r_{i, j}$ ’s, for $j = 2, \dots, n_{i} - 1, i = 1, \dots, N_{c}$ ; see the supplemental information.

For purpose of Bayesian modeling using rstan package,⁵⁸ the informative prior distributions are imposed on $μ_{m}^{c}$ , e.g. $Unif (- 5, 5)$ , and $σ_{m}^{c}$ , e.g. $Unif (0, 5)$ , for $m = 1, \dots, M$ and $c = 1, \dots, C$ . In the ‘Modeling the logarithm of step length section, we explain the choice of prior distributions for our example of T-cell movements further.

Data analysis

Initially, we apply the proposed methodology to the training dataset. The R code is available through the supplemental materials. T-cell trajectories obtained from three, of five observed, experiments. More advanced cases for the training There are $N = 352$ (i.e. 117 + 132 + 103) trajectories available in the training dataset. Each T cell stays in the observation space for a minimum of 28 steps and a maximum of 179 steps, where a step is considered one change in a T cell’s location. The experiments are sampled for the training data and the out-of-sample data randomly while the more advanced form of this procedure is left for the future research direction.

Clustering 3D spatio-temporal T-cell trajectories

Generally speaking, human brain tends to choose more than required number of clusters; see the literary works of Kim etal.,⁵⁹ and Tseng and Wong.⁶⁰ In real applications, it is required to consider the biologist advice on the possible number of groups. Our interest, for now, is to study the Lévy type of movement versus the Brownian motion while persistence versus non-persistence is of importance. Following (3) and (4), RP and CM are obtained for each trajectory. The ${RP}_{i}$ is calculated for each trajectory by the consideration of $ψ_{0} = \frac{π}{2}$ . This threshold is chosen because of our notion of persistency of the movement of a T cell on its own path; see the work of Textor et al.⁷ The CM is used to cluster trajectories initially to having long increments versus short increments groups; see the supplemental information. The within-cluster sum of squared error for Lévy type of movement and Brownian motion suggests that having two subcategory for each cluster makes sense; see the supplemental information. Figure 3 shows the histogram of RP obtained from the training dataset. The general scheme proposed in Table 1, obtained by applying $k$ -means algorithm to Lévy type of movement and Brownian motion separately, results in clusters where only two-state NHMM converged. Applying the gamma distribution for the step length, the Bayesian model converges in all clusters while the lognormal density diverged even after manipulating the prior distribution on hyper parameters; see the ‘Modeling the angular change section for the detail on modeling and prior distributions. This leads to choose another method for clustering trajectories. In order to get more accurate insight of data, because of highly skewed distribution of RP, and also the limited number of trajectories with $RP \leq 80 %$ , the training data is then partitioned into highly persistent group ( $RP \geq 87 %$ ) and non-persistent group ( $RP < 87 %$ ) where $87 %$ is the median of RP observed in the training dataset.

Consistent with our discussion in the ‘Clustering spatio-temporal 3D trajectories to capture between cell heterogeneity section, clusters are defined as being highly persistent/long increment (HP/LI), highly persistent/short increment (HP/SI), non-persistent/long increment (NP/LI), and non-persistent/short increment (NP/SI) using RP and CM; for mathematical convenience, clusters are numbered $c = 1, \dots, 4$ , though below we refer to them by their abbreviation (e.g. HP/LI) in discussion to aid readability.

There are 93, 84, 137, and 38 trajectories in each cluster, respectively. Figure 4 shows a T-cell trajectory with HP/LI behavior. Figures corresponding to trajectories from clusters of HP/SI, NP/LI, and NP/SI are presented in the supplemental information. Let ${η η}^{c} = {η_{i, j} | j = 2, \dots, n_{i} - 1, i = 1, \dots, N_{c}}$ and ${R R}^{c} = {R_{i, j} | j = 2, \dots, n_{i} - 1, i = 1, \dots, N_{c}}$ , for $c = 1, \dots, 4$ . Clusters HP/LI and HP/SI both capture highly persistent cells, and hence we expect them to have similar ${η η}^{c}$ distributions. Likewise for NP/LI and NP/SI, both capturing non-persistent cells. For ${R R}^{c}$ , reflecting step length distributions, we expect HP/LI and NP/LI to share similar distributions, and likewise for HP/SI and NP/SI where for HP/LI and NP/LI we expect higher variation compared to HP/SI and NP/SI based on the CM presented in (4).

Modeling the 3D spatio-temporal T-cell trajectories

Having observed the distribution of radial angles exhibited across the entire training dataset to be approximately uniform, we modeled these angles using a uniform distribution in (7). We inferred that, having clustered trajectories based on RP and CM, these stratified groups of trajectories each still exhibit this uniformity; see Figure 2 and the supplemental information.

Modeling the angular change

As such, we explicitly model cluster-specific values for ${η η}^{c}$ jointly with ${R R}^{c}$ , but not radial angles. Due to the label switching problem in mixture models, see the work of Gelman et al.,⁶¹ HHMM and NHMM are estimated only for ${η η}^{c}$ , for $c = 1, \dots, C$ , we leave the more advanced case for our future research direction. The label-switching problem in our context means the probability of observing a random variable under each component of either the HHMM or NHMM is very close such that the STAN model cannot converge. We will, however, show that it is necessary to use more complex models in our data analysis.

Among the possible settings for the number of states of the HMM models and clusters, Table 2 gives different models applied to the training data when we consider (a) two clusters, HP and NP; (b) four clusters, HP/LI, HP/SI, NP/LI, and NP/SI; and (c) without clustering. In addition, we consider a different number of states for the HMM’s, i.e. two states: persistent and non-persistent; three states: persistent, moderately persistent, and non-persistent; and one state: no within cell switching behavior.

Table 2.

Different models fitted to trajectories considering two clusters (HP and NP), four clusters (HP/LI, HP/SI, NP/LI, and NP/SI), and without applying the clustering scheme.

Fitted models
Clustering scheme	One-state	Two-state	Three-state
Two clusters	No HMM	HHMM and NHMM	HHMM and NHMM
Four clusters	No HMM	HHMM and NHMM	HHMM and NHMM
Without clustering	No HMM	HHMM and NHMM	HHMM and NHMM

HP: highly persistent; NP: non-persistent; LI: long increment; SI: short increment; HMM: hidden Markov model; NHMM: non-homogenous hidden Markov model; HHMM: homogenous hidden Markov model

Informative uniform prior distributions are defined on $α_{m}^{c}$ and $β_{m}^{c}$ , for $m = 1, \dots, M$ and $c = 1, \dots, C$ . For two-state models, it is expected that the beta distribution has a probability mass below 0.5 for the persistent movement while for the non-persistent one, the behavior is reversed. Accordingly, it is assumed $α_{1}^{c} \sim Uniform (0, 10)$ , $β_{1}^{c} \sim Uniform (0, 4)$ , $α_{2}^{c} \sim Uniform (0, 4)$ , and $β_{2}^{c} \sim Uniform (0, 10)$ where $(α_{m}^{c}, β_{m}^{c})$ are non-persistent and highly persistent parameters, respectively, for $m = 1, 2$ and $c = 1, \dots, C$ with $C = 2, 4$ . A similar assumption is made when three-state models are examined in the data analysis. Uniform prior distributions are defined on $α_{m}^{c}$ and $β_{m}^{c}$ , for $m = 1, 2, 3$ and $c = 1, \dots, 4$ . To this end, it is assumed $α_{1}^{c} \sim Uniform (0, 10)$ , $β_{1}^{c} \sim Uniform (0, 4)$ , $α_{2}^{c} \sim Uniform (0, 10)$ , $β_{2}^{c} \sim Uniform (0, 10)$ , $α_{3}^{c} \sim Uniform (0, 4)$ , and finally $β_{3}^{c} \sim Uniform (0, 10)$ where $(α_{m}^{c}, β_{m}^{c})$ are non-persistent, middle, and persistent parameters, respectively, for $m = 1, 2, 3$ and $c = 1, \dots, 4$ . $α_{2}^{c}$ is set to be equal to $β_{2}^{c}$ . This is done in order to produce a hidden state modeling angular change in between persistence and non-persistence, a middle persistent state. The informative uniform prior distributions $α^{c} \sim Uniform (0, 10)$ and $β^{c} \sim Uniform (0, 10)$ are considered for the simple one-state model.

Dirichlet prior distributions are considered with the data-driven parameters for the initial probabilities and also rows of the transition matrices. These parameters were taken to be the percentage of data with the same feature, for example.^54–56 For $M = 2$ , a dirichlet distribution, $D (d_{1}^{c}, d_{2}^{c})$ , is considered for $π^{c} = [π_{1}^{c}, π_{2}^{c}]$ where $d_{l}^{c}$ is the percentage of ${ψ ψ}^{c}$ in $[\frac{2 - l}{2} π, \frac{2 - l + 1}{2} π)$ with $π$ included in the first state, for $l = 1, 2$ and $c = 1, \dots, C$ with $C = 2, 4$ . Dirichlet prior distributions are considered for each row of the transition matrices where the parameters come from the percentage of the training data with the same feature. For $M = 3$ , the initial probability $π^{c} = [π_{1}^{c}, π_{2}^{c}, π_{3}^{c}]$ , a dirichlet distribution, $D (d_{1}^{c}, d_{2}^{c}, d_{3}^{c})$ , is considered where $d_{l}^{c}$ is the percentage of ${ψ ψ}^{c}$ in $[\frac{3 - l}{3} π, \frac{3 - l + 1}{3} π)$ , for $l = 1, 2, 3$ and $c = 1, \dots, 4$ , with $π$ included in the first state. Similar to two-state HHMM, we define dirichlet prior distributions on each row of transition matrices with the parameters coming from the percentage of the training data with the same feature. For the priors corresponding to two-state or three-state NHMM, the information about T cells nearby is taken into account to set parameters of the dirichlet prior distributions for the transition matrices.

Each model was estimated using the rstan package⁵⁸ in R version 3.1.1. The results obtained from different clusters and also without clustering produce an Rhat, a statistic defined in the STAN model, equal to one, which indicates the convergence of the Bayesian model. The trace plot of parameters of the training data with two-state HHMM is given in the supplemental information. The estimates of initial probabilities and also transition matrices of the fitted two-state HHMM and NHMM of the scaled angular change under different schemes are given in corresponding Tables in the supplemental information.

The influence of the existence of a T cell in the neighborhood is measured by considering two different transition matrices for each of the four clusters on the type of movement, i.e. persistent and non-persistent; see Tables in the supplemental information. The observed changes in transition matrices show the same pattern for the long increment clusters. In HP/LI and NP/LI, the existence of T cells in a neighborhood increases the probability of staying in non-persistent type of movement while for HP/SI, the result is reversed. It is worth mentioning that the unusual value of $P_{2}^{2}$ is estimated using a small sample size. We do not observe this relationship between the estimated matrices when only considering two clusters. Results given in the supplemental information show that there is not a big difference between $P_{1}^{c}$ and $P_{2}^{c}$ , for $c = 1, 2$ . Similar to two-state NHMM, the existence of a T cell in the neighborhood has the same pattern in long clusters (not shown here).

For the purpose of model selection and comparison, the Watanabe-Akaike information criterion (WAIC); see the work of Gelman et al.,⁶² of the fitted models of $η η$ are displayed in Table 3. The three-state NHMM four-cluster model is preferred over other competing models in terms of WAIC, followed by three-state HHMM four-cluster models. The WAIC of two-state HHMM and NHMM four cluster models are approximately equal to those from the three-state HHMM and NHMM four cluster models. The fact that WAIC generally improves with more clustering and hidden states suggests that we may need more complex models to take into account the inherent differences between and within T cells motion. The WAIC computed from STAN returned NaN for the models without clustering. We attribute this to the fact that the posterior likelihoods in the unclustered case are so small that the log posterior is not well defined.

Table 3.

WAIC of different models fitted to $η η$ of the training dataset.

Models	WAIC
Two-state HHMM, 2 clusters	$-$ 137.62
Two-state NHMM, 2 clusters	$-$ 138.43
Two-state HHMM, 4 clusters	$-$ 461.56
Two-state NHMM, 4 clusters	$-$ 463.71
Three-state HHMM, 2 clusters	$-$ 149.24
Three-state NHMM, 2 clusters	$-$ 150.05
Three-state HHMM, 4 clusters	$-$ 477.04
Three-state NHMM, 4 clusters	$- 479.61$

Bold value indicates the smallest WAIC.WAIC: Watanabe-Akaike information criterion; NHMM: non-homogenous hidden Markov models; HHMM: homogenous hidden Markov models.

Modeling the logarithm of step length

Next, we model ${R R}^{c}$ in each cluster, for $c = 1, \dots, C$ with $C = 2, 4$ , and also $R R$ of the training data. To this end, ${η η}^{c}$ is decoded in each cluster using (8), for $c = 1, \dots, C$ with $C = 2, 4$ , and also $η η$ of the training data. Following (5), Bayesian analysis is implemented and normal distributions are fitted to ${R R}^{c}$ , for $c = 1, \dots, C$ with $C = 2, 4$ , and $R R$ of the training data in different states. The initial observation of data shows the log of step length has its mean in $(- 1, 1)$ and variance in $(0, 1.5)$ . Accordingly, the informative prior distributions are applied to $μ_{m}^{c}$ , $Unif (- 5, 5)$ , and $σ_{m}^{c}$ , $Unif (0, 5)$ , for $m = 1, \dots, M$ and $c = 1, \dots, C$ . The Bayesian model converged with Rhat equal to one. The trace plots of parameters for two-state NHMM of the training data showing the chains are well mixed are displayed in the supplemental information. The weight of each hidden state, $w_{m}^{c}$ , obtained by decoding the angular change in each cluster, the estimates of model parameters for the distributions of angular change and the logarithm of step length, $α_{m}^{c}$ , $β_{m}^{c}$ , $μ_{m}^{c}$ , and $σ_{m}^{c}$ , are obtained; see the supplemental information. We present results here for four clusters and without clustering for the two-state NHMM and HHMM models; persistent and non-persistent type of movements. In addition, we present similar information for two clusters in the supplemental information. The normal distribution for the logarithm of step length in non-persistent mode is shifted to the left with higher variance. The difference between the weights and model parameters for two-state HHMM and NHMM was negligible aside from the $(α_{1}^{c}, β_{1}^{c})$ column, specifically the $α_{1}^{c}$ entry. In the ‘Simulation studies section, it is shown the NHMM version gives slightly better prediction of motility metrics of interest.

The proposed models were also evaluated by comparison to the angles and step lengths obtained from out-of-sample data. Two experiments with $N_{o}$ = 412 (i.e. 182 + 230) trajectories are considered as out-of-sample data. Clusters of trajectories of out-of-sample data are obtained following the method of the ‘Clustering spatio-temporal 3D trajectories to capture between cell heterogeneity section. Figure 5 shows the fitted two-state NHMM on angular changes obtained from out-of-sample data and corresponding P-P plots for four clusters. Here, the P-P plots show the empirical cumulative distribution function (ECDF) obtained from the out-of-sample data plotted versus the cumulative distribution function obtained from the fitted models. They are a better fit to the scaled angular change of these clusters compared to the one obtained from the training data or the two clusters in the supplemental information. Based on P-P plots, two-state and three-state HHMM and NHMM have almost the same performance in capturing the true underlying behavior of the angular change within each cluster. The P-P plots in the supplemental information and information criteria (not shown here) indicate that the one-state model is inadequate.

Figure 6 shows the fitted two-state NHMM on the values of $R R$ obtained from the corresponding out-of-sample data. The P-P plots show the proposed model can capture the behavior of step length for the training data and also in four clusters using a two-state HHMM.

The fitted one-state models for $R R$ were inadequate for out-of-sample data when fitted models are obtained with or without clustering scheme imposed to the training dataset in the supplemental information. Also, the three-state HHMM and NHMM do not perform well for ${R R}^{c}$ in clusters of out-of-sample data in the supplemental information, for $c = 1, \dots, 4$ . After decoding, the weights of the non-persistent hidden state for HP/LI and HP/SI are close to zero which shows the three-state HHMM and NHMM suffer from overfitting. Due to this, we focus on two-state HHMM and NHMM in the data analysis and simulation studies below.

In this paper, we have a random selection for the out-of-sample data and training data. The other choice might be finding all possible selections for the training data, say $(\binom{5}{3})$ , and fitting the models separately, averaging over estimates of model parameters, a more accurate model might be obtained. This is not the purpose of the current work.

Simulation studies

The simulation studies that we present in this section were implemented to evaluate the predictive properties of the proposed models. We aimed to compare the 3D motility metrics obtained from the out-of-sample and the training T-cell trajectories with the metrics obtained from simulated T-cell trajectories obtained from HHMM and NHMM for purpose of model comparison. It is worth mentioning that we have been able to build the 3D trajectories except for the first increment due to not having the corresponding angles; however, we work with the origin as the initial locational observation for each T cell to solve the issue.

To this end, $R = 500$ sets of T-cell trajectories are simulated under different scenarios. In particular, we used the fitted two-state HHMM and NHMM models with (a) four clusters, (b) two clusters, and (c) without clustering. In each case, the radial angle was simulated from a uniform distribution on $[0, 2 π]$ . We empirically evaluate the legitimacy of modeling radial angles as being uniformly distributed, by contrasting simulated trajectories drawing radial angles from a uniform distribution with those radial angles from real T-cell trajectories. The estimated model parameters given in corresponding Tables in the supplemental information were used to simulate T-cell trajectories.

We simulated T-cell trajectories so that their time of entry into, and duration within the observation space match the trajectories in the training data. The covariate information about the proximity of other T cells was taken from the experimental data. Initially, the cluster membership of the simulated trajectories are randomly allocated based on the observed percentage of T cells in each cluster of the training data. Then, trajectories are simulated using estimates of model parameters in the corresponding cluster in the case of two and four clusters.

Similarly, trajectories are generated using models fitted to the training data without clustering. An initial step is generated from the origin, after which the next $X$ , $Y$ , and $Z$ increments are built subsequently from the model. Figure 7a shows 140 collapsed time simulated trajectories obtained from the two-state NHMM four-cluster model. Compared with the real trajectories given in Figure 7b, the simulated trajectories appear qualitatively similar.

As we discussed in the ‘Observing real T-cell motility through live 3D imaging section, T cells cannot be tracked when they leave the observation space. We are interested to see if we remove the limitation of dimension how our analysis is affected. To this end, T-cell trajectories are simulated repeatedly until they remain within the observation region for the same number of steps as the corresponding trajectory in the training sample. Figures 8(a and b) show the scaled angular change and logarithm of step length, respectively, obtained from the out-of-sample data fitted by the distributions of simulated values from two-state HHMM when simulated T cells are recorded only if they stay inside the observation space (yellow line), and without this restriction (red line). Our observation shows there is almost no difference between the distributions of step lengths and angular changes in these two settings, indicating that the effect of the finite viewing volume does not noticeably affect the step length and angular changes observed in our simulation studies. Figure 8b also suggests that the proposed model under-predict the value of $R$ . The effect of this under-prediction can be tracked later in motility metrics. It appears then that our sample is sufficiently large to capture the cells’ overall behavior.

Recalling that T-cell trajectories are often summarized using motility metrics, we computed root mean squared displacement (RMSD) over time, outreach ratio, displacement ratio, Hurst index, asphericity, and straightness of simulated trajectories under different models and also the out-of-sample and the training data using MotilityLab package⁶³ in R version 3.1.1. Figures 9 to 11 and Figures 26 to 40 in the supplemental information show the observed values of these measures obtained by averaging over all trajectories in the training dataset and out-of-sample dataset compared to the ones obtained from simulated trajectories using the proposed models. We also compared fitting lognormal and gamma distributions to the step length.⁶⁴

Based on these metrics, the proposed clustering scheme improves the overall performance of the two-state HHMM and NHMM. Moreover, the two-state NHMM slightly improves over the two-state HHMM in predicting the motility metrics, which validates the importance of the information provided by heterogenous transition matrices. Hence, T cells in close proximity to one another influence their motility.

We find little difference between the accuracy of the predicted metrics when the observed radial angle is used to generate T-cell trajectories in the simulation studies. This justifies the choice of uniform distribution for the radial angle; see, for example, Figures 9 to 11 and Figures 26 to 34 in the supplemental information. In order to keep the diagrams clean and informative, we only show the error bars of training data, the behavior of out-of-sample data is similar in terms of variation; see Figures 35 to 40 in the supplemental information.

The clusters that we consider in our analysis are not clusters proposed by the $k$ -means clustering, we have then compared the performance of $k$ -means derived clusters with our proposed modification. Due to having only the gamma distribution converged for the step length in these clusters, we lose accuracy of our model for some metrics, not shown here. This is consistent with our observations of using the gamma distribution for modeling the step length.

Discussion

We propose a novel methodology that explicitly accounts for heterogeneous motility characteristics of real T cells moving in physiologically realistic 3D compartments. We believe that these heterogeneities could, at a population level, profoundly impact T-cell search performance. Our method is based on an initial clustering of T-cell trajectories that aims to separate persistent cells exhibiting Lévy flight type behavior from cells with less persistent, more Brownian, type motion. In an ongoing project, we plan to compare the performance of this scheme with the model based clustering scheme available in the literature. After this clustering step, HMMs are applied to a spherical representation of the trajectories, where the model estimation is conducted with a Bayesian modeling approach. We found in both in-sample evaluations using information criteria and out-of-sample comparisons utilizing motility metrics that generally the best performing models among those considered were NHMMs that (a) account for potential contact between neighboring cells, (b) use two states, which can be thought of as persistent and non-persistent type of movements, and (3) use four preliminary clusters to better account for between-cell heterogeneity in trajectories. In summary, these results suggest that the T cells studied indeed exhibit heterogeneity of motion in two ways: Both between cells, i.e. that there are inherent differences in motion across cells, which clustering captures, as well as within cell tracks, i.e. that cells seem to exhibit a latent switching behavior, which is captured by HMMs. However, the motility metrics suggest that the trajectories may be more inhomogeneous than the model allows. The multivariate normal distribution might be a solution to this issue which we are working on as the continuation of this project. Our understanding of T-cell trajectories is that they display different types of movement under observation, specifically in terms of persistent versus non-persistent, as well as in their step lengths. It seems in terms of stratifying according to short versus long increments, we obtain small though meaningful differences.

Both these aspects of heterogeneity are important to overall population motility dynamics. This is clear from the analysis because the models that performed the best in terms of predicting out-of-sample motility metrics, and in in-sample model evaluations (information criteria) were those using multiple states and clusters.

Our models are able to capture the first few steps of T-cell motion quite well, and, as is apparent from analyzing Figures 9 to 11, this accuracy persists first several minutes for measures like straightness, asphericity, Hurst index, displacement ratio, and outreach ratio. However, accuracy degrades after several minutes of observation of the cells as measured by motility metrics. The accuracy of RMSD is not as good as the other motility metrics, still improved over other models, though this measure has been criticized in the literature; see the work of Textor et al.⁷ We have considered several reasons for this which cannot be easily accounted for within this framework, and which we consider as potential directions for future work. One is that the collagen environment in which the cells are imaged is not homogeneous: it contains some dense “fibers” as well as empty pockets. Invariably when cells are observed over an extended period, they encounter such inhomogeneities, which affect their long term motion and can explain the loss in predictive power of our models over longer time horizons; see the work of Jones et al.⁶⁵

Additionally, our analysis used data only describing the center, or sometimes called centroid, of each cell, which we considered as the cell’s position. T cells are far from being perfect spheres, and from one image to the next their morphology can change dramatically. This affects efforts to try and measure their center. In terms of the centroid data, this can mean that cells that are essentially motionless can appear to have erratic angular changes. This, to an extent, is captured by our model, since motionless cells are then more easily categorized by the HMM as exhibiting non-persistent motion, but a more rigorous method including measurement error to try and filter out this effect is worthy of further work. Also, although we used a decoding approach to estimate the hidden state of the step length following modeling of the angular change, in future work we aim to model the spherical coordinates jointly conditional on the assumed latent/hidden variable that accounts for T-cell switching behavior. Our preliminary results show that the correlation among the step length and the angular change is different in distinct states and clusters. So far, it does not increase the accuracy of the prediction of motility metrics substantially. Our results of this ongoing research, not shown here, indicate $ψ_{i, j}$ ’s is the leading term in T-cell movement.

In spite of these areas of future work, we consider our proposed method a valuable contribution to the T-cell motility modeling literature, especially as it considers between-cell and within-cell heterogeneity of motion. Also, the methodology proposed in this paper can be implemented in oncology studies as well as any type of 3D movements.

Supplemental Material

sj-csv-1-smm-10.1177_09622802231172041 - Supplemental material for Modeling three-dimensional T-cell motility using clustering and hidden Markov models

Supplemental material, sj-csv-1-smm-10.1177_09622802231172041 for Modeling three-dimensional T-cell motility using clustering and hidden Markov models by Elaheh Torkashvand in Statistical Methods in Medical Research

Supplemental Material

sj-csv-2-smm-10.1177_09622802231172041 - Supplemental material for Modeling three-dimensional T-cell motility using clustering and hidden Markov models

Supplemental material, sj-csv-2-smm-10.1177_09622802231172041 for Modeling three-dimensional T-cell motility using clustering and hidden Markov models by Elaheh Torkashvand in Statistical Methods in Medical Research

Supplemental Material

sj-csv-3-smm-10.1177_09622802231172041 - Supplemental material for Modeling three-dimensional T-cell motility using clustering and hidden Markov models

Supplemental material, sj-csv-3-smm-10.1177_09622802231172041 for Modeling three-dimensional T-cell motility using clustering and hidden Markov models by Elaheh Torkashvand in Statistical Methods in Medical Research

Supplemental Material

sj-csv-4-smm-10.1177_09622802231172041 - Supplemental material for Modeling three-dimensional T-cell motility using clustering and hidden Markov models

Supplemental material, sj-csv-4-smm-10.1177_09622802231172041 for Modeling three-dimensional T-cell motility using clustering and hidden Markov models by Elaheh Torkashvand in Statistical Methods in Medical Research

Supplemental Material

sj-csv-5-smm-10.1177_09622802231172041 - Supplemental material for Modeling three-dimensional T-cell motility using clustering and hidden Markov models

Supplemental material, sj-csv-5-smm-10.1177_09622802231172041 for Modeling three-dimensional T-cell motility using clustering and hidden Markov models by Elaheh Torkashvand in Statistical Methods in Medical Research

Supplemental Material

sj-pdf-6-smm-10.1177_09622802231172041 - Supplemental material for Modeling three-dimensional T-cell motility using clustering and hidden Markov models

Supplemental material, sj-pdf-6-smm-10.1177_09622802231172041 for Modeling three-dimensional T-cell motility using clustering and hidden Markov models by Elaheh Torkashvand in Statistical Methods in Medical Research

Footnotes

Acknowledgments

Elaheh Fatima Torkashvand acknowledges the significant support of Dr. Joel A. Dubin and Dr. Gregory Rice for this work and also the continuation of this research. We would like to thank STAN group, especially Bob Carpenter, and WestGrid support staff, Ali Kerrache, for their support during this project. Mark N. Read is supported by David and Judith Coffey LifeLab. Elaheh Fatima Torkashvand also thank Jorge Luis Galeano Niño and Maté Biro from EMBL Australia for providing the data and biological comments of this work.

Author's note

Elaheh Torkashvand is also affiliated at College of Public Health, Ohio State University, Columbus, OH, USA.

Data availability statement

The data that supports the findings of this study are available from the Journal of Statistical Methods in Medical Research through

Declaration of conflicting interests

The author declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author received the following financial support for the research, authorship and/or publication of this article: This work was supported by Australian Research Council project grant DP180102458. Gregory Rice and Joel A. Dubin acknowledge the support of the Discovery Grant from Natural Sciences and Engineering Research Council (NSERC) of Canada, Grant RGPIN-2016-10477, RGPIN-2016-10476, and RGPIN-2014-05911.

ORCID iD

Elaheh Torkashvand

Supplemental information

The Supplemental material includes figures and tables supporting the validity of our choice of HHMM, NHMM, and clustering. A few additional results are also included. The supplemental information is available online.

References

Miller

Wei

Parker

, et al. Two-photon imaging of lymphocyte motility and antigen response in intact lymph node. Science 2002; 296: 1869–1873.

Masopust

Schenkel

. The integration of t cell migration, differentiation and function. Nat Rev Immunol 2013; 13: 309.

Krummel

Bartumeus

Gérard

. T cell migration, search strategies and mechanisms. Nat Rev Immunol 2016; 16: 193.

Kienle

Lämmermann

. Neutrophil swarming: an essential process of the neutrophil tissue response. Immunol Rev 2016; 273: 76–93.

Weninger

Biro

Jain

. Leukocyte migration in the interstitial space of non-lymphoid organs. Nat Rev Immunol 2014; 14: 232.

Galeano

Kwan

RYQ

Weninger

, et al. Antigen-specific T cells fully conserve antitumour function following cryopreservation. Immunol Cell Biology 2016; 94: 411–418.

Textor

Peixoto

Henrickson

, et al. Defining the quantitative limits of intravital two-photon lymphocyte tracking. Proc Natl Acad Sci USA 2011; 108: 12401–12406.

Schienbein

Gruler

. Langevin equation, Fokker-Planck equation and cell migration. Bull Math Biol 1993; 55: 585–608.

Niculescu

Textor

De Boer

. Crawling and gliding: a computational model for shape-driven cell migration. PLoS Comput Biol 2015; 11: e1004280.

10.

Manolopoulou

Matheu

Cahalan

, et al. Bayesian spatio-dynamic modeling in cell motility studies: Learning nonlinear taxic fields guiding the immune response. J Am Stat Assoc 2012; 107: 855–865.

11.

Textor

Sinn

De Boer

. Analytical results on the beauchemin model of lymphocyte migration. BMC Bioinform 2013; 14(6): S10.

12.

Beauchemin

Dixit

Perelson

. Characterizing T cell movement within lymph nodes in the absence of antigen. J Immunol 2007; 178: 5505–5512.

13.

Banigan

Harris

Christian

, et al. Heterogeneous cd8+ T cell migration in the lymph node in the absence of inflammation revealed by quantitative migration analysis. PLoS Comput Biology 2015; 11: e1004058.

14.

Mokhtari

Mech

Zitzmann

, et al. Automated characterization and parameter-free classification of cell tracks based on local migration behavior. PLoS One 2013; 8: e80808.

15.

Ester

Kriegel

Sander

, et al. A density-based algorithm for discovering clusters in large spatial databases with noise. In: Kdd. vol. 96; 1996. p. 226–231.

16.

Kalnis

Mamoulis

Bakiras

. On discovering moving clusters in spatio-temporal data. In: International Symposium on Spatial and Temporal Databases. Springer; 2005. pp. 364–381.

17.

Lee

Han

Whang

. Trajectory clustering: a partition-and-group framework. In: Proceedings of the 2007 ACM SIGMOD international conference on Management of data. ACM; 2007. pp. 593–604.

18.

Jebara

Song

Thadani

. Spectral clustering and embedding with hidden markov models. In: European Conference on Machine Learning. Springer; 2007. pp. 164–175.

19.

Zhang

Deng

Van de Weghe

. Clustering spatio-temporal trajectories based on kernel density estimation. In: International Conference on Computational Science and Its Applications. Springer; 2014. pp. 298–311.

20.

Ester

Kriegel

Sander

, et al. Incremental clustering for mining in a data warehousing environment. In: VLDB. vol. 98. Citeseer; 1998. p. 323–333. Sydney.

21.

Wang

. Mining spatial-temporal clusters from geo-databases. In International Conference on Advanced Data Mining and Applications 2006. pp. 263–270. Xi'an, China: Springer.

22.

Yeh

Chen

. Profiling moving objects by dividing and clustering trajectories spatiotemporally. IEEE Trans Knowl Data Eng 2013; 25: 2615–2628.

23.

Pan

Wang

, et al. Mining regular behaviors based on multidimensional trajectories. Expert Syst Appl 2016; 66: 106–113.

24.

Plank

Auger-Méthé

Codling

. Lévy or not? analysing positional data from animal movement paths. In: Dispersal, individual movement and spatial ecology. Springer; 2013. p. 33–52.

25.

Dustin

. Stop and go traffic to tune T cell responses. Immunity 2004; 21: 305–314.

26.

Méhes

Vicsek

. Collective motion of cells: from experiments to models. Integr Biol (Camb) 2014; 6: 831–854.

27.

Metzner

Mark

Steinwachs

, et al. Superstatistical analysis and modelling of heterogeneous random walks. Nat Commun 2015; 6: 7516.

28.

Read

Bailey

Timmis

, et al. Leukocyte motility models assessed through simulation and multi-objective optimization-based model selection. PLoS Comput Biology 2016; 12: e1005082.

29.

Bartolucci

Farcomeni

. A discrete time event-history approach to informative drop-out in mixed latent markov models with covariates. Biometrics 2015; 71: 80–89.

30.

Raffa

Dubin

. Multivariate longitudinal data analysis with mixed effects hidden markov models. Biometrics 2015; 71: 821–831.

31.

Patterson

Basson

Bravington

, et al. Classifying movement behaviour in relation to environmental conditions using hidden markov models. J Anim Ecol 2009; 78: 1113–1123.

32.

Langrock

King

Matthiopoulos

, et al. Flexible and practical modeling of animal telemetry data: hidden markov models and extensions. Ecology 2012; 93: 2336–2342.

33.

Patterson

Parton

Langrock

, et al. Statistical modelling of individual animal movement: an overview of key methods and a discussion of practical challenges. ArXiv e-prints 2016;1603.07511.

34.

Pedersen

Righton

Thygesen

, et al. Geolocation of north sea cod (gadus morhua) using hidden markov models and behavioural switching. Can J Fish Aquat Sci 2008; 65: 2367–2377.

35.

Thygesen

Pedersen

Madsen

. Geolocating fish using hidden markov models and data storage tags. In: Tagging and tracking of marine animals with electronic devices. Springer; 2009. p. 277–293.

36.

Vermard

Rivot

Mahévas

, et al. Identifying fishing trip behaviour and estimating fishing effort from vms data using bayesian hidden markov models. Ecol Modell 2010; 221: 1757–1769.

37.

Mitani

Sato

Ito

, et al. A method for reconstructing three-dimensional dive profiles of marine mammals using geomagnetic intensity data: results from two lactating weddell seals. Polar Biol 2003; 26: 311–317.

38.

Liu

Zidek

Trites

, et al. Bayesian data fusion approaches to predicting spatial tracks: Application to marine mammals. Ann Appl Stat 2016; 10: 1517–1546.

39.

Laplanche

Marques

Thomas

. Tracking marine mammals in 3d using electronic tag data. Methods in Ecology and Evolution 2015; 6: 987–996.

40.

Orellana

Matrai

Leck

, et al. Marine microgels as a source of cloud condensation nuclei in the high arctic. Proc Natl Acad Sci USA 2011; 108: 13612–13617.

41.

Sims

Southall

Humphries

, et al. Scaling laws of marine predator search behaviour. Nature 2008; 451: 1098.

42.

Humphries

Queiroz

Dyer

JRM

, et al. Environmental context explains lévy and brownian movement patterns of marine predators. Nature 2010; 465: 1066.

43.

Spezia

. Bayesian variable selection in non-homogeneous hidden markov models through an evolutionary monte carlo method. Comput Stat Data Anal 2020; 143: 106840.

44.

Koki

Meligkotsidou

Vrontos

. Forecasting under model uncertainty: Non-homogeneous hidden markov models with pòlya-gamma data augmentation. J Forecast 2020; 39: 580–598.

45.

Meligkotsidou

Dellaportas

. Forecasting with non-homogeneous hidden markov models. Stat Comput 2011; 21: 439–449.

46.

Holsclaw

Greene

Robertson

, et al. Bayesian nonhomogeneous markov models via pólya-gamma data augmentation with applications to rainfall modeling. Ann Appl Stat 2017; 11: 393–426.

47.

Maruotti

Rocci

. A mixed non-homogeneous hidden markov model for categorical data, with application to alcohol consumption. Stat Med 2012; 31: 871–886.

48.

Beltman

Marée

AFM

Lynch

, et al. Lymph node topology dictates T cell migration behavior. J Exp Med 2007; 204: 771–780.

49.

Huisman

Koedijk

Kool

CJM

, et al. Tail-index estimates in small samples. J Bus Econ Stat 2001; 19: 208–216.

50.

Efron

Tibshirani

. An introduction to the bootstrap. CRC press, 1993.

51.

Hastie

Tibshirani

Friedman

, et al. The elements of statistical learning: data mining, inference, and prediction. vol. 2. Springer, 2009: 485–579.

52.

Torkashvand

. Supplement to “modeling three dimensional T cell motility using clustering and hidden markov models”. Stat Methods Med Res 2023.

53.

Neal

, et al. Mcmc using hamiltonian dynamics. Handbook of Markov Chain Monte Carlo 2011; 2.

54.

Berger

Moreno

Pericchi

, et al. An overview of robust Bayesian analysis. Test 1994; 3: 5–124.

55.

Torkashvand

Jozani

Torabi

. Constrained Bayes estimation in small area models with functional measurement error. Test 2016; 25: 710–730.

56.

Schad

Betancourt

Vasishth

. Toward a principled bayesian workflow in cognitive science. arXiv preprint arXiv:190412765 2019.

57.

Zucchini

MacDonald

Langrock

. Hidden Markov models for time series: an introduction using R. 150. CRC press, 2016.

58.

Stan Development Team. RStan: the R interface to Stan 2018;R package version 2.17.3.

59.

Kim

Ogden

Kim

. A method to identify differential expression profiles of time-course gene data with fourier transformation. BMC bioinformatics 2013; 14: 310.

60.

Tseng

Wong

. Tight clustering: a resampling-based approach for identifying stable and tight patterns in data. Biometrics 2005; 61: 10–16.

61.

Gelman

Carlin

Stern

, et al. Bayesian data analysis. 2. FL: CRC press Boca Raton, 2014.

62.

Gelman

Hwang

Vehtari

. Understanding predictive information criteria for Bayesian models. Stat Comput 2014; 24: 997–1016.

63.

Textor

Dannenberg

Berry

, et al. Motilitylab: Quantitative analysis of motion 2016;R package version 0.2-5.

64.

Cho

Bowman

North

. A comparison of gamma and lognormal distributions for characterizing satellite rain rates from the tropical rainfall measuring mission. J Appl Meteorol 2004; 43: 1586–1597.

65.

Jones

CAR

Cibula

Feng

, et al. Micromechanics of cellularized biopolymer networks. Proc Natl Acad Sci USA 2015; 112: E5117–E5122.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.43 MB

0.56 MB

0.54 MB

0.38 MB

0.74 MB

2.94 MB