Data-Driven Model Construction for Anisotropic Dynamics of Active Matter

The dynamics of cellular pattern formation is crucial for understanding embryonic development and tissue morphogenesis. Recent studies have shown that human dermal fibroblasts cultured on liquid crystal elastomers can exhibit an increase in orientational alignment over time, accompanied by cell proliferation, under the influence of the weak guidance of a molecularly aligned substrate. However, a comprehensive understanding of how this order arises remains largely unknown. This knowledge gap may be attributed, in part, to a scarcity of mechanistic models that can capture the temporal progression of the complex nonequilibrium dynamics during the cellular alignment process. The orientational alignment occurs primarily when cells reach a high density near confluence. Therefore, for accurate modeling, it is crucial to take into account both the cell-cell interaction term and the influence from the substrate, acting as a one-body external potential term. To fill in this gap, we develop a hybrid procedure that utilizes statistical learning approaches to extend the state-of-the-art physics models for quantifying both effects. We develop a more efficient way to perform feature selection that avoids testing all feature combinations through simulation. The maximum likelihood estimator of the model was derived and implemented in computationally scalable algorithms for model calibration and simulation. By including these features, such as the non-Gaussian, anisotropic fluctuations, and limiting alignment interaction only to neighboring cells with the same velocity direction, this model quantitatively reproduce the key system-level parameters—the temporal progression of the velocity orientational order parameters and the variability of velocity vectors—whereas models missing any of the features fail to capture these temporally dependent parameters. The computational tools we develop for automating model construction and calibration can be applied to other systems of active matter.

Figure 1

(a) Schematics of the two types of substrates (isotropic and nematic) tested in this study. On the nematic substrate, the molecular alignment (denoted by a red double-sided arrow), is parallel to $x$. (b) A stitched view of a snapshot of long-term cell imaging where cell cytoplasm is labeled in red. (c) Correlation plots in exploratory data analysis, where $x, y$ velocity components are plotted against those of the neighboring average in the previous step. The slope of the fit (red solid lines) is smaller than 1 (black dashed lines). (d) Parameter estimation of the slope parameters, $\omega$, the slope value of the red solid line in (c), and their $95\%$ confidence intervals (shaded area). The statistical tests yield four potential features (green blocks), representing new features added to the classical Vicsek model in fitting the current data set. (e) The grid-based system to search for the neighbors. In order to find the neighbor of the target cell $i$ (circled in red), a radius of $r$ will be searched. The computational procedure is simplified by only searching through neighboring squares connected to the square cell $i$ is located in. This workflow generates various physical parameters via simulation to compare to quantities computed from the experiment.

Figure 2

Parameters estimated from the data set. (a) The standard deviation of the cell velocities is shown in yellow circles and green triangles. The light-colored histogram denotes the cell number count. [(b), (c)] The velocity distribution at a representative time of 20 h. The dashed black lines denote the best-fit normal distribution. The solid orange lines denote the best-fit Laplace distribution. (d) The blue squares denote the polar order parameter, and the polar histograms of the velocity angle at two different time points are plotted in (e) and (f).

Figure 3

A permutation F-test shows that the variance of the fluctuation along $x$ and $y$ directions is indeed unequal at time = 20 h.

Figure 4

Shapiro-Wilk test for normality. The figure shows the $p$-value of individual frames, which is a statistical measurement applied to validate whether or not the observed data follow a normal distribution. A small $p$-value indicates statistical significance. The inset shows a representative quantile-quantile (QQ) plot at time = 20 h to compare the sample quantiles of normalized velocity with the theoretical quantiles from the standard normal distribution. The black curve indicates the linear curve with slope 1 (normal).

Figure 5

Effects of accounting for different neighbor ensembles. [(a)–(c)] Parameters and correlation plots using the neighbor ensemble of the classical Vicsek model for time = 20 h. [(d)–(f)] The scenario where only cells traveling in the same direction are considered. [(b), (c), (e), (f)] The correlation between the cell at $t$ versus the mean of the velocity of its neighbors at $t-1$. (g) Correlations between velocity $v_{i,x}\left(t\right), v_{i,y}\left(t\right)$ of any cell $i$ and mean velocity of the neighboring cells from the previous step ${\overline{v}}_{i,x}(t-1), {\overline{v}}_{i,y}(t-1)$, where the $95\%$ confidence intervals, shown in faint shades, are computed based on Fisher transformation of correlation coefficients. Two scenarios are shown, either by including all neighbors (bottom two curves), or only the same direction neighbors (top two curves). (h) The estimated weights, equivalent to the fitted slopes, ${\widehat{\omega }}_x\left(t\right), {\widehat{\omega }}_y\left(t\right)$ from these two scenarios and their confidence intervals are plotted.

Figure 6

Reproducing orientational order parameters and mean absolute deviation of velocity with different simulation models at a representative radius $r$ = 75 $\mu \mathrm{m}$. The left, middle, and right panels compare simulations using the Gaussian, Laplace, and GGD fluctuation with experimental observations. For each type of fluctuation distribution, simulations of different interaction rules are presented. [(a)–(c)] The orientational order parameter and [(d)–(f)] the average absolute deviation of the velocities.

Figure 7

Distributions of random fluctuations between the observations and simulated models at a representative radius $r=75 \mu \mathrm{m}$ and time = 20 h with the distributions of the fluctuation in the simulation following [(a), (d)] Gaussian, [(b), (e)] Laplace, and [(c), (f)] GGD.

Figure 8

Change of order parameters due to the change of simulation parameters with ${\tau }_x=0.02$ and $w_x=0.7$. The color bar corresponds to the simulated order parameter at long times. The dashed lines and circles denote the actual path of the experimental results plotted against the phase diagram. The colors of the markers correspond to the order parameter as denoted by the color bar. Red and blue X's denote the beginning and the end of the trajectory, respectively.

Figure 9

Parameter estimation with interval (shaded) for a generalized Gaussian distribution following Eq. (8) along the $x$ direction (yellow), and along the $y$ direction (green).

Figure 10

Change of order parameters due to the change of simulation parameters with (a) ${\tau }_x$ = 0.01 and (b) ${\tau }_x$ = 0.04, which cover the range of the observed $\tau$'s. $w_x=0.7$ is assumed for both cases. Red and blue X's denote the beginning and the end of the trajectory.

Figure 11

Several typical simulation progressions, where the first 20 steps are regarded as burn-ins and left out in computing the averages in the phase diagram.

Figure 12

Reproducing [(a)–(c)] orientational order parameter and [(d)–(f)] average absolute deviation of velocities from different simulation models at a representative radius $r=75 \mu \mathrm{m}$ for experiments of cells migrating on a nematic substrate when starting at a higher initial density.

Figure 13

Reproducing [(a)–(c)] orientational order parameter and [(d)–(f)] average absolute deviation of velocities with different simulation models at a representative radius $r=75 \mu \mathrm{m}$ for experiments of cells moving on an isotropic substrate, which do not align, and have isotropic velocity.

Figure 14

Reproducing orientational order parameters with different radii. The left, middle, and right panels compare simulations using the Gaussian, Laplace, and GGD fluctuation with experimental observations. Open square symbols show the order parameter calculated from experimental data; the dark blue solid line, yellow dash-dotted line, and red dashed line represent simulation results of different radii 75, 50, and 25 $\mu \mathrm{m}$.