Cell trajectory analysis

Goal of our study was to investigate cellular motilities at the individual cell level to determine whether stimulant dose or its gradient is the main regulator of the motility. Time lapse images collected at N time points with time interval ∆t provides the cellular positions at snapshot times. We characterized the motilities of the cells by computing (Fig. 1) i) average total instantaneous velocity of the cells v_tot, which is the velocity independent of movement direction, ii) average instantaneous velocity component along the ligand gradient direction v_d, which corresponds to directed (chemotactic) velocity, and iii) average instantaneous velocity along the direction perpendicular (normal) to gradient direction, v_p.

Schematic diagram for cell velocity components: i) average total instantaneous velocity of the cells v_tot, ii) average instantaneous velocity component along the EGF gradient direction v_d, and iii) average instantaneous velocity along the normal to gradient direction, v_p.

Average total instantaneous velocity v_tot of a cell is computed from its measured locations as vtot=1N−1∑k=2N|r→k−r→k−1|/Δt where k is the frame index of the image time-series data. This average instantaneous velocity is equal to the ratio of the total trajectory length of a cell to the total elapsed time. Note that v_tot is a scalar quantity and does not have direction dependence.

Average instantaneous velocity component along the EGF gradient direction v_d was computed similar to v_tot except that v_d is directional and corresponds to a vector quantity. If the gradient is along the y-direction, then vd=1N−1∑k=2N(yk−yk−1)/Δt where y_k is the y-position of the cell in the k^th frame of the time-series image data.

Lastly the average instantaneous velocity component perpendicular (orthogonal) to the gradient direction v_p was also computed. If the gradient is along y-axis, v_p would be the average instantaneous velocity along the x-direction. As the directed motility along this perpendicular direction is expected to be random, the v_p calculations were used as the negative control to check for the biases in the experiments.

As an additional measure of the directed motion, we analyzed the cell trajectories by computing the directed persistence of cells as well⁹. Directed persistence (P_d) of a cell is defined as the ratio of cell’s net movement along the gradient direction to the total trajectory length over the course of time. If the gradient is along the y-direction, then Pd=(yN−y1)/∑k=2N|r→k−r→k−1|. Directed persistence ranges from −1 (cell persistently moves opposite to the gradient direction) to 1 (cell persistently moves along the gradient direction). Small absolute P_d values correspond to mostly random cell movement.

The focus of our analysis was to characterize the velocity of a cell as a function of the local EGF concentration and its gradient at cell’s location. This uniquely integrated approach enabled us to separate the contributions of these interlinked system variables. We achieved this by experimentally measuring the cell velocities and computing the local EGF distribution in microchip devices in fluid dynamics simulations using the COMSOL Multiphysics software version 5.1. This allowed us to relate the velocities of cells to local EGF concentration and its gradient at the individual cell level. Collected individual cell data was then analyzed using multivariate non-linear regression method to establish the dependence of cell velocities to EGF concentration and EGF gradient. Pursued analysis consisted of the following steps:

Identify the cells in time-lapse image frames and obtain their position: We used the commercial Volocity image analysis software for automated cell tracking. An image sequence was created from the time lapse images (czi format). The “Find Object” task was used to select the cells in the image sequence by selecting an appropriate intensity limit. We then used the “Track” task to follow the motion of selected objects. Settings for the tracking were: a) ignore new objects, b) ignore static object, and c) automatically join broken tracks. Maximum cell to cell distance for tracking was 19 µm. This setup was appropriate to track cells correctly throughout the whole imaging duration. Quality of tracking was further validated visually. Using GFP expressing cells made cell tracking process easier. We note that several groups have used manual tracking for analyzing the cell motility, which can be highly susceptible to errors. For example, Huth et al. showed that manual cell tracking can lead to miscalculation of migration rates by up to 410%⁵¹. Using an automated tracking system enabled us to avoid such errors.

Compute the EGF concentration [L] and gradient (∇L) distribution within microdevices and determine the ([L], ∇L) values at individual cell locations: Time-dependent distribution of chemoattractant was computed using the COMSOL program. Because the used growth media was mostly water based, properties of the bulk media in the devices were considered equivalent to water. Diffusion constant of EGF was chosen as 134 µm²/s⁵². Two types of microdevices which were used to collect the motility data differed in terms of the presence of flow. Therefore, their computational analysis was slightly different.

Double chamber device: EGF applied to one of the side chambers is distributed through the narrow channel to the whole device by diffusion (Fig. 2). As there was no fluid flow, there was no EGF transport through convection. Therefore, EGF distribution was computed using the Fick’s equation for diffusion

where c_i is the concentration of dissolved species (EGF) and D_i is its diffusion constant. As in the experiments, 30 ng/ml EGF was added to one of the side chambers and the other chamber had no EGF. Simulation of ligand distribution in the device revealed that EGF concentration varied in the 4–24 ng/ml range across the narrow channel after 30 minutes of initial stabilization period and ligand distribution resulted in an almost linear gradient across the narrow central channel. This linear gradient remained stable and varied only in a very narrow range of ~20–24 ng/ml/mm over the study period (Fig. 2C). From the experimental trajectory files cell positions were identified and corresponding EGF concentration and EGF gradients were interpolated from COMSOL simulation results with the help of a MATLAB program.

(A) Schematic diagram of the double chamber device, top and side views. (B) EGF distribution across the cross section of the device after 300 min as computed in the simulations. The color bar at the right shows the EGF concentration in ng/ml. (C) EGF concentration across the narrow mixing channel at different time points during the course of time-series experiments as computed in the simulations. Results shown in parts (B) and (C) were computed using the COMSOL program with the conditions detailed in the text. Size of the imaged area was approximately 1 mm by 1 mm. Because it is the region with the most linear gradient, imaged area was centered at the middle of the channel (i.e., around the 0.5 mm point in the graph).

Y-channel device: Experimental conditions for the y-channel device (Fig. 3) included flow through the device. Mass transport of dissolved species can occur through two mechanisms in these devices: Diffusional transport due to concentration gradient and convectional transfer due to the bulk fluid flow. Therefore, diffusional transport of diluted species is coupled to laminar flow. This was accounted for by coupling the Fick’s diffusion equation with the convection mass transfer equation to compute the spatio-temporal distribution of EGF in the y-channel device

where u is the fluid velocity and N_i is the local flux. Because of the continuous pumping, laminar flow was at steady state in our system. At this limit, streamlines along velocity field do not cross each other and this prohibits convectional mass transfer to adjacent fluid layers. Therefore, only the diffusional mass transfer contributes to the formation of EGF gradient in the y-channel device.

(A) Schematic diagram of the y-channel device. Middle inlet was sealed and not used and ligand was injected through the lower inlet. (B) Steady-state EGF distribution within the device when the input EGF concentration was 15 ng/ml. Cellular motilities were imaged at the center axis at four different downstream positions along the device after the y-junction. Size of the imaged area was approximately 1 mm by 1 mm. Dashed lines indicate the imaged positions: blue (5 mm), green (10 mm), red (15 mm), and cyan (20 mm) where the numbers in the paranthesis refer to the distance from the y-junction. (C) EGF concentration along the direction perpendicular to flow at the imaged positions which are marked with dashed lines in part (B). Imaged area was centered at the middle of the channel (i.e., around the 1.5 mm point in the graph (C)). As can be seen, formed gradient at the imaged areas are almost linear and vary along the length of the device. Results shown in parts (B) and (C) were computed using the COMSOL program with the conditions detailed in the text.

To minimize the shear stress on the cells due to fluid flow, a low flow rate of 0.5 µl/min through each of two inlets of the y-channel device was utilized (Fig. 3). At this flow rate, computed shear stress on cells was 0.0022 dyne/cm². Note that this value is about 10 times lower than the values typically used in literature³⁸. To cover a wide EGF concentration and gradient ranges, experiments were performed with several different EGF concentrations through one of the inlets (cf., cell motility experiments section above). Overall, pursued y-channel device based experiments covered a range of 0–35 ng/ml for EGF concentration and 0–37 ng/ml/mm for its gradient. It should be noted that these stimulation ranges overlap with the ranges that were used with the double chamber device.

Using the computed EGF distribution in the devices (cf., step ii), we determined the ([L],∇L) values at individual cell locations which were determined in the first step. A MATLAB program developed in house was used for this step. The [L] and ∇L exposure values for each cell were calculated using the mean position of the cell over its trajectory during the imaging period. Since cells typically moved only short distances during monitoring time (typically less than 50 µm over 200 minutes), choice of mean cell position as the coordinate for ([L],∇L) estimation was reasonable. To test whether this choice biased our conclusions, analysis was repeated by estimating the ([L],∇L) values the cell is exposed to by using the starting coordinates of a cell at the beginning of the imaging time (results not shown). The conclusions of analyses with different coordinate choices agreed with each other.

Perform multivariate non-linear regression analysis to determine the dependence of cell velocities on EGF concentration [L] and its gradient ∇L: As described above (Fig. 1), we computed multiple velocity properties (v_tot, v_d, v_p, persistence) to characterize the cell motion. For each of these properties, we fitted the results for all individual cells to regression models. Four models were considered:

These non-linear models varied in complexity and inclusion of [L] and ∇L as independent system variables. These models consist of constant term C (corresponding to [L] or ∇L independent response), EGF concentration dependent term (C_conc[L]^n_conc), EGF gradient dependent term (C_grad∇L^n_grad), and a cross-term that depends on both EGF concentration and gradient (C_cross[L]^n_c,c∇L^n_g,c). These terms have the power form, allowing for accommodating possible non-linear dependencies. Optimal values of the regression model parameters were obtained using MATLAB curve fitting toolbox where root mean square error (RMSE) was the optimization criteria. Fitted regression models were further examined by ensuring that the residual (i.e., model prediction − experimental data) values for individual cells had a random distribution with no obvious systemic trends.

It should be noted that, since the concentration must be finite to be able to create a gradient, a complete coverage of the ligand concentration and gradient space is not feasible technically. However, our experiments were performed for conditions that made it possible to cover a large portion of the mentioned space. As can be seen in Figure 4, coverage of the variable space was substantial, which increased the confidence of our conclusions.

Coverage of the EGF concentration, [EGF], and EGF gradient, ∇EGF, space in our motility experiments. Each point in the figure corresponds to a cell and reports the EGF concentration and gradient that the cell was exposed to during the experiment. Exposure of the cells to EGF concentration and gradient covered a wide range of the parameter space, which made our study unique and enabled us to investigate the effect of EGF concentration and gradient changes independently.