BMDCs were generated as previously described (12, 13). Briefly, tibias and femurs from BALB/c mice (8–12 weeks old, female) were flushed, and red blood cells were lysed with ammonium-chloride-potassium (ACK) lysis buffer (Gibco). Bone marrow cells were plated in 24 well tissue culture well plates (1 × 10⁶ cells/mL) in complete medium containing RPMI 1640, supplemented with 5% fetal bovine serum (FBS), 1% antibiotic-antimycotic solution, 1% HEPES buffer, and 0.1% 2-mercaptoethanol (all reagents were purchased from Gibco) containing 20 ng/mL recombinant mouse granulocyte-macrophage colony-stimulating factor (GM-CSF; Peprotech). The medium was completely replaced with fresh GM-CSF every two days. On day six, non-adherent and loosely adherent cells were collected by gentle pipetting and transferred to Petri dishes. After one day of culture, immature BMDCs, which appeared as floating cells, were collected. Phenol-red-free RPMI 1640 medium (Gibco) was used for fluorescence microscopy experiments, including cell height measurements and cell viability tests. To generate mature dendritic cells (mDCs), immature DCs (imDCs) were stimulated with 100 ng/mL lipopolysaccharide (LPS, LPS-EB Ultrapure; Invivogen) for 30 min. Cells were carefully washed three times and incubated for 6 h in fresh complete medium, as described in previous studies (12, 13). After incubation, floating cells were harvested as mDCs. All animal experiments were performed according to protocols approved by the Institutional Animal Care and Use Committee of the Ulsan National Institute of Science and Technology (UNISTIACUC-19-15).

The specific surface markers of imDCs and mDCs were characterized by flow cytometry (Cytoflex, Beckman Coulter; Supplementary Figure S2 ). The upregulated expression levels of co-stimulatory molecules, CD86, CD80, and CD40, antigen-presenting molecule MHC class II (I-A/I-E), and the chemokine receptor CCR7 were evaluated as DC maturation markers, and CD11c and CD11b were analyzed as dendritic cell markers. The following antibodies were used in flow cytometry experiments: anti-CD86-FITC (GL1, 1:200), anti-CD40-PE (1C10, 1:200), anti-MHC class II (I-A/I-E)-FITC (M5/114.15.2, 1:200), and anti-CCR7-PE (4B12, 1:200) were purchased from Thermo Fisher Scientific, and anti-CD11b-APC (M1/70, 1:200), anti-CD11c-APC (HL3, 1:200), and anti-CD80-PE (16-10A1, 1:200) were purchased from BD Biosciences. The acquired data were analyzed using the FlowJo software (BD).

Under-agarose migration by gel confinement was performed as described previously (13, 27). Briefly, the gel confiner consisted of a custom-designed PDMS structure and low-melting agarose gel. The PDMS structure consisted of a 10:1 ratio for the main body and a 30:1 ratio for the sticky PDMS-coated bottom. Before casting the gel solution, the PDMS structure was placed in a Petri dish. Low-melting agarose (2.4%) was dissolved in phenol-red free HBSS buffer and heated at 80°C for 20 min, and the resulting solution was cooled at room temperature to 40°C. The same volume of 2× conditioned medium (RPMI 1640, 10% FBS, 2% HEPES buffer, 2% antibiotic-antimycotic solution, and 0.2% 2-mercaptoethanol) was mixed to a final concentration of 1.2% and cast to the PDMS structure. The gel was cured for 20 min at room temperature. The cured gel confinement was incubated overnight in a cell culture incubator with complete medium. Subsequently, to prevent non-specific cell-to-cell interactions in the migration assay, 800 cells in a small drop of cell suspension were seeded on 10 mm diameter coverslips with 20 μg/mL bovine fibronectin-coating. The cells were incubated for 30 min in a cell culture incubator to enable them to settle on the substrate. Subsequently, the cells were carefully covered by gel confinement, and motility was imaged after 1 h. In the experiments to measure cell motility, the field of view (FOV) was located within the radial distance of 1.7 ± 0.6 mm from the center, and there was no significant difference in cell behavior observed at the center and edge positions.

Agarose gel blocks were prepared at concentrations of 1.2% (w/v) with the same gel confinement fabrication. The mechanical properties of the gels were determined using a rheometer (MCR502 WESP; Anton Paar). The gel height was approximately 1.0 mm. The shear moduli (G, Pa) were analyzed using the relationship between shear stress and shear strain before gel disruption, and Young’s moduli (E, kPa) were calculated using the following equation: E = 2(1 + v)G, with the Poisson’s ratio (v) of agarose set as 0.5 (28).

To check the reproducibility of the confined environment, the height of the fluorescence-stained DCs under 2D gel confinement was measured using a laser scanning confocal microscope (A1R, Nikon), as shown in Supplementary Figure S3 . The DC suspensions were stained using DiO (Invitrogen) for three minutes and washed thrice with complete medium. After three minutes of recovery, stained DCs were seeded onto the substrate and covered by gel confinement. The 3D confocal image was acquired using a 100x plan apo lens, and the z interval was 0.5 μm. Cell height was measured manually using a NIS Element (Nikon).

Live-cell imaging was performed using an inverted microscope (Eclipse Ti-E, Nikon) configured with a 10x dry objective lens and sCMOS camera (Flash4.0, Hamamatsu). The cells were imaged for 24 h, and bright-field images were obtained every 1 min. The recorded images were processed using the z-max intensity projection for non-labeled automatic tracking (29). The sample focal plane was focused, and two more image sequences were obtained along the z-axis. Subsequently, rolling-ball background subtraction was performed, and the maximum intensity projection overlaid the z-stacks, facilitating contrast enhancement for cell segmentation. During the experiments, an incubator (Chamlide HK; Live Cell Instrument) maintained the microscope at 37°C with 95% humidity and 5% CO₂.

For tracking, cells in the pre-processed image were automatically detected using the IMARIS (Bitplane) ‘Spots’ function. Cells were identified using intensity thresholding and size estimation as 20 μm diameter particles. The position of the center of mass was tracked, and tracking errors such as misconnections between different objects or misrecognized objects were corrected manually.

We obtained the following raw trajectories of DC migration:348 (imDCs, training dataset), 383 (mDCs, training dataset), 449 (imDCs, testing dataset), and 350 (mDCs, testing dataset). Among them, a few raw trajectories were not tracked perfectly, and as a result, a few data points were absent in the trajectories. Therefore, we polished such trajectories. The number of imperfect trajectories (fraction) was identified as 16 (4.5%, imDC training dataset), 20 (5.22%, mDCs training dataset), 26 (5.79%, imDC testing dataset), and 20 (5.71%, mDCs testing dataset). The missing data points were attributed to tracking errors that mainly arise from cell-to-cell interactions during attachment between two distinct cells and partially from the finite field of view, resulting in the incomplete observation of cell movement near the edge of the observation window. Thus, we manually analyzed the trajectories using the following rule: First, we chose an arbitrary time window (≤ 10 min) and applied the spline interpolation provided in the Python library “scipy” package to the raw data (30). If the untracked period was ≥10 min, we split the trajectory into two before and after attachment. Splitting an imperfectly tracked trajectory into two independent short trajectories does not seriously alter the statistics and the detection of migration patterns of single trajectories because this portion is approximately 4%–6% of the trajectory data. We also evaluated the interval of missing points, such as 20 min or 30 min time windows; however, the results shown in analyzed motility data did not differ qualitatively. If we consider three dynamical modes, the statistical properties of the modes are robust against the interval of missing points owing to the relatively few imperfect trajectories obtained.

After data polishing, we removed inappropriate track samples in our ML analysis that failed to satisfy any of the following three criteria as previously referred (12, 14): Mean track speed ( < V D > t ≥ 1.5   μm / min ), track duration ( T o b s ≥ 60   min ), and maximal displacement from the starting position ( max | R → | ≥ 20   μm ).

We constructed a hybrid ML kernel combining two machines each used as unsupervised and supervised learning. Firstly, the K-means clustering was applied to generate pseudo labels for 1 h long segmented trajectories as a training data set. Then, randomly chosen 2000 trajectories (36.6%) are used for training XGBOOST. The trained XGBOOST predicted the label for unknown trajectories in the segment pool with a 99.6% accuracy.

In our novel kernel design, the classical K-means clustering with silhouette score systematically determines the number of clusters, and the decision-tree-based XGBOOST provides the importance of the feature. The benefits of two distinct types of machines are providing a deeper understanding of trajectories. Thus, our designed kernel possibly adapts to similar single-particle-tracking data analysis with user-defined features.

The silhouette score has been widely used to infer the number of clusters in ML algorithms (31). With the given cluster numbers as a hyperparameter, we measured the mean distance of a given data point within the cluster to which it belongs. For data point i, the mean distance within the cluster is estimated as D W ( i ) = 1 | C I | ∑ j ∈ C I d ( i , j ) ,   where d ( i , j ) represents the Euclidian distance between the i and j   data points in the cluster C I , and | C I | is the size of the cluster. We estimated the dissimilarity of the data point i by measuring the minimum distance between the data point and other clusters, which is defined as D D ( i ) = min { C J } 1 | C J | ∑ j ∈ C J d ( i , j ) , where i ∈ C I , j ∈ C J , and { C J } is the set of all neighboring clusters for J ≠ I . The silhouette score for data point i is defined as s ( i ) = D D ( i ) − D w ( i ) max [ D w ( i ) , D D ( i ) ] . Thus, the average silhouette score lies between [ − 1 , 1 ] . The limit of 1 describes a sample (i.e., cluster) that is well separated from neighboring clusters.

To study the role of myosin II and the Arp2/3 complex, the molecular inhibitors blebbistatin (20 μM, Sigma-Aldrich), CK666 (100 μM, Sigma-Aldrich) and SMIFH2 (12.5 μM, Sigma-Aldrich) were used during DC migration. Each inhibitor was mixed in an uncured gel (40°C) before casting. After gel casting and curing at RT, the gel confiner was incubated in complete medium with the same inhibitor concentration overnight. DMSO (0.1% (v/v); Sigma-Aldrich) was used as a negative control. Subsequently, the DCs located on the substrate were carefully covered with an inhibitor-containing gel confiner. Cell motility assays were conducted after 1 h of incubation with covering by the gel confiner, similar to the approach described above.

Ethidium homodimer (EthD-1; Thermo Fisher Scientific) was used to evaluate DC viability. After 24 h of incubation, DCs under gel confinement were stained with EthD-1 in complete medium (500 nM) in an incubator for 20 min. Fluorescent-stained dead cells were counted manually using epifluorescence microscopy with a 20× Plan Apo lens (Nikon TiE). Cell survival rate was calculated as the number of stained cells per the number of whole cells in the field of view. More than 50 cells were analyzed for each experiment, and three independent experiments were performed.

Unless stated otherwise, all data represent the mean ± SEM of three independent experiments for each condition. Normality was determined using the D’Agostino and Pearson tests. The Mann-Whitney and Kruskal-Wallis tests, along with Dunn’s post hoc test, were used to determine statistical significance. To prevent type I errors due to the large number of samples and confirm statistical significance, randomly selected samples were used to obtain the P-values. The criteria for random selection followed the hypothesized sample size, with a statistical power of 0.8. All statistical analyses were performed using the Origin Pro 2020 software (OriginLab).

To investigate the motile behaviors of DCs, we used primary mouse bone marrow-derived DCs (BMDCs) and characterized the DC phenotype (Materials and Method, Supplementary Figure S2 ). We used a gel confiner for long-term (24 h) monitoring of an unbiased population of freely migrating cells ( Figures 1A, B ). The gel confiner allowed for precise control over the degree of confinement in experiments involving multiple chips ( Supplementary Figure S3 ). The detailed fabrication steps are described in our previous reports (13, 27). Cellular movements were imaged every 1 min using bright-field microscopy with a millimeter-scale broad field of view (1.3 × 1.3 mm²; Supplementary Movies 1 , 2 ). Using live-cell tracking data ( Figure 1C ), we found that mDCs showed faster and more persistent motility than imDCs, which is consistent with the results of previous studies ( Figures 1D, E ) (12, 13). However, outliers were consistently present ( Figure 1F ), and these atypical phenomena were often neglected in the pooled population analysis ( Supplementary Figure S1 ).

Our ML kernel was designed to analyze individual, unperturbed, heterogeneous cellular motility. To construct the ML kernel, we prepared two independent datasets: one for training and the other for testing. The input data were pre-processed to remove trivial errors (Materials and Method). After pre-processing, we obtained trajectories of 301 (imDCs) and 348 (mDCs) for training and 460 (imDCs) and 371 (mDCs) for testing.

To analyze dynamic heterogeneity in a single DC migration, we segmented the trajectory data into 1 h long tracks and obtained a total of 8787 tracks from the testing dataset ( Figure 2A ). For the training dataset, we prepared 5741 segmented tracks. From the segmented tracks, we extracted five features quantifying motility (radius of gyration, end-to-end distance, and average kinetic energy) and directionality (asphericity A and turning angle fluctuation Δ θ ) as the input bases for the machine kernel ( Figure 2B ). The mathematical details of these features that enable quantification of the motility propensity and directionality of DC migration are as follows:

With the five extracted features from the segmented trajectories, we constructed the input data of size 5741×5 for training the machine ( Figures 2A, B ). First, we applied K-means unsupervised clustering (32, 33) to the input data to examine the number of distinct motility patterns of DCs. Next, we used the results from K-means clustering as pseudo-labels and cross-validated these using the supervised learning algorithm XGBOOST (34) (Materials and Method, Figure 2C ). We point out that the combination of the unsupervised and supervised learning is advantageous when the cell-to-cell migration variation originates from temporally heterogeneous motion of a cell. After trained with the input data, the XGBOOST machine enables us to analyze unsegmented raw trajectories and decipher the migration pattern in the 1 h interval time windows along the trajectory ( Figure 2D ). We validated the ML method and its output by demonstrating that pattern classification is robust to variations in feature selection and their possible combinations.

In Figure 2C , we examined the average silhouette score to determine the number of distinct clusters in the migration patterns (31). The plot shows that the silhouette score reached the maximum when the number of distinct clusters was three if we used all five features introduced in our ML analysis. To assess the robustness of the results, we performed ML analysis with various combinations of input features. For this test, we used a training dataset of 5741 segmented trajectories. In Supplementary Figure S4 , we evaluated the silhouette scores for five distinct combinations of the input features. For all cases, we confirmed that the DC migration data most likely consisted of three dynamic modes. We temporarily refer to the three groups as I, II, and III based on their population size.

The XGBOOST machine also allowed us to evaluate the statistical properties of the migration trajectories in each group with feature importance and the distribution of the five features ( Supplementary Figure S5 ). It turns out that asphericity, A, is the most important feature among the five when determining the group assignment ( Supplementary Figure S5A ). It plays a critical role in differentiating group I from groups II and III ( Supplementary Figure S5B ). The features, R g ,     R e t e , and E are effective in distinguishing group II and group III ( Supplementary Figures S5C–E ). These features quantify DC migration motility. In terms of R g ,     R e t e , and E , groups I and II exhibited almost indistinguishable “slow” motility patterns, whereas group III, with high motility, differed from the other two groups. The variance of the turning angles showed a minor contribution, with a similar unimodal profile for all three groups ( Supplementary Figure S5F ). In Table 1 , a few important dynamic features of groups I, II, and III are summarized. Based on averaged dynamic properties, we assigned the dynamic modes of groups I, II, and III to slow-diffusive (SD), slow-persist (SP), and fast-persist (FP) phenotypes, respectively ( Table 1 ; Figure 3 and Supplementary Movies 3 - 5 ). Note that the name has been given by the clustering results in five-dimensional feature space reflecting averaged dynamical properties, thus, slow and fast do not assign from their absolute velocity.

DCs in the SD mode perform anti-persistent walks at a slow average speed, whereas FP mode cells show directional random walks at a fast average speed. The SP mode differed from the other two modes. It had a directional walk analogous to that of FP mode cells. However, the feature properties differed significantly from each other ( Supplementary Figure S5 ). In terms of average speed, the SP mode was similar to the SD mode.

We investigated the dynamic properties of the SD, SP, and FP modes based on commonly used statistical observables, including mean-squared displacements (MSDs) and Van-Hove self-correlation functions. We also examined observables such as the turning angle heat map and density map of the phase space, ( V n , Δ θ n ) . Additionally, we investigated the zigzag-like patterns of DC migration in detail.

To investigate how the maturation status of DCs influences the distribution of the three migration modes, we assigned the migration mode 1 h interval traces over time until the cell finally escaped from the field of view ( Figure 4A and Supplementary Movies 6 , 7 ). The speed and fraction of the three migration modes suggest that DC migration occurs mostly in the SD and SP modes, and their speeds are unaffected, regardless of their maturation status ( Figures 4B, C ). The FP mode was more frequently observed in mDCs than in imDCs.

Cell speed and persistence are controlled by the regulation of actin dynamics (41, 42). Because myosin II-mediated contractility and actin polymerization key mechanism of cell motility (43), we sought to test how actin mediators influence migration mode distribution ( Figure 5A and Supplementary Movies 8 - 15 ). During experiments, they survived similarly to the controls over the next 24 h ( Supplementary Figure S12A ).

Although the overall speed of both imDCs and mDCs was reduced ( Supplementary Figure S12B ), which is consistent with previous studies (14, 44–46), the FP mode was absent, and the mean speed of both SD and SP modes was decreased significantly for both cases of myosin II and Arp2/3 inhibition ( Figure 5B ). Remarkably, myosin II inhibition increased the SP mode of migration in the imDC population ( Figure 5C ), whereas Arp2/3 inhibition increased the SD mode of migration of mDCs. Therefore, after the inhibition, mode fractions seem to converge to a similar value irrespective of maturation status.

Next, we determined the mode lifetime, length of mode continuity, and how molecular inhibition influenced the mode longevity. Interestingly, myosin II inhibition reduced imDC SD mode durations while increasing mDC SD and SP mode lifetimes ( Figure 5D ). Arp2/3 complex inhibition enhanced the lifespan of mDCs in the SD mode but not that of imDCs, demonstrating that mode lifetime dynamics and the control of actin dynamics are highly dependent on DC maturation status.

Unlike myosin II and Arp2/3 inhibition cases, the FP mode was maintained in DCs treated with a Formin inhibitor, SMIFH2 (12.5 μM) ( Figure 5A ). While the speed of the SD and SP modes was not changed much, there was a decrease in the speed of the imDC FP mode, but not in mDC ( Figure 5B ). Analysis of the mode fraction showed that the ratio of SD and SP mode was comparable to Arp2/3 inhibition cases. However, FP mode was preserved only in DCs treated with SMIFH2 ( Figure 5C ). It is noteworthy that SMIFH2 treatment also impacts mode lifetime. A comparable pattern to CK666 treatment was observed when DCs were treated with SMIFH2. However, an increase in the FP mode lifetime of imDC was observed, whereas no such change was detected in mDC ( Figure 5D ).

To study the transition dynamics of the DC migration modes, we constructed a transition matrix ( Figure 6A ) that showed several noteworthy features. First, self-transition predominates, signifying that DCs move with strong mode persistence. Compared with imDCs, mDCs have higher self-transition rates for the SP and FP migration modes and lower self-transition rates for the SD migration mode. This tendency explains why mDCs moved more persistently than imDCs; mDCs avoided SD migration in favor of SP or FP migration. Second, SD was the most recurrent migration mode for both imDCs and mDCs. The highest sum of influx rates was observed for the SD mode ( ∑ i ∈ { SD ,  SP ,  FP } P i → S D ). Third, maturation increases the transition to FP mode. Comparing the FP column of the transition matrix of imDCs and mDCs, we found that every component of the transition from SD, SP, and FP to FP was significantly increased.

Irrespective of maturation, DC migration is dominated by two slow modes (SD and SP), suggesting that DCs usually move at a slow speed with or without directional persistence. The high energy-consuming FP migration mode was observed intermittently, and the frequency of occurrence increased after maturation. In Figure 6B , we plotted the population fraction of mode-preserving (single-mode) and mode-changing (multi-mode) cellular migration tracks. Interestingly, mDCs have a larger fraction of single-mode migration tracks than imDCs, indicating that maturation participates in stabilizing the innate migration of DCs.

Next, we examined the cross-transition dynamics among different migration modes ( Figure 6C ). Mode changes can occur in one or two steps. We emphasize that the cross-mode transition diagram should be distinguished from the total transition matrix, as shown in Figure 6A . For both imDCs and mDCs, a single-step change was the dominant pathway for the transition between the SD and SP modes. For the transition from the SD to FP modes, the one-step pathway was the most common for imDCs. For mDCs, the fractions for the one-step (SD→FP) and two-step (SD→SP→FP) transitions were similar. From the SP to FP transition, the two-step transition (SP→SD→FP) was the primary pathway for imDCs, whereas both one- and two-step transitions were similar for mDCs. Taken together, imDCs predominantly followed the unicyclic SD→FP→SP→SD transition ( Supplementary Movie 16 ). In contrast, mDCs showed no transition directionality.