1. Introduction
Coordinated cell migration in the human body is essential throughout the whole life. Cells need to be in the correct positions during the development of an embryo or for the formation of organs, as well as for keeping them functioning by means of regeneration and repair [
1,
2]. In addition, migration plays a key part in the immune response of the organism [
3,
4]. However, the negative aspects of cell migration can be observed in pathological processes, for example, during cancer invasion [
2]. Endothelial cells form in many vessels a barrier between the bloodstream and the surrounding tissues. Here, migration helps in stabilizing this barrier as well as repairing it if needed. Furthermore, a coordinated movement of endothelial cells is fundamental during the formation of new vessels [
1]. On the pathophysiological side, there are many diseases like atherosclerosis, edema formation or respiratory distress syndrome, which arise from the dysfunction of the aforementioned barrier [
1].
When modeling cell migration, several biological properties and environmental conditions have to be accounted for. For instance, cells are sensitive to chemotactic stimuli [
4] but also to mechanical cues like the stiffness of the substrate, which influences mean cell speeds [
5,
6]. In collective settings, the regulation of adhesion molecules mediates the connection to neighboring cells as well as the adhesion to the extracellular matrix [
1]. While endothelial cells can migrate individually under certain conditions, they also move in groups or sheets [
1]. In a 2D tissue culture, human umbilical vein endothelial cells (HUVECs) also show collective movement patterns. Thereby, cells not only move around while staying connected but also influence each other in their behavior [
7]. This generates the complex dynamics of active individual and collective cell migration.
Here, we analyze the dynamics of confluent endothelial cell migration and quantify changes in the mean squared velocity over time. We also determine the effect of the average cell density and its changes on mean cell velocities within the monolayer by coupling them in a combined mathematical model.
2. Materials and Methods
HUVECs were isolated from umbilical cords according to the previously described protocol [
8]. The cells were then collected from T25 flasks, counted and filled into a microscope slide, which was placed in the incubator. After 20 h, and 4 h before the start of the experiment was scheduled, the cell nuclei were stained with Hoechst 33342 (final concentration 50 ng/mL).
The live-cell imaging microscope Zeiss Axio Observer.Z1 / 7 enabled us to observe the cells at 37 °C and 5% carbon dioxide for 48 h. By taking several micrographs of a large area shortly after one another and stitching them together, we gained pictures of a connected cell layer, comprising an area of more than 42 mm2. We acquired phase contrast and fluorescence images of the cells every 10 min to capture the signals of the stained cell nuclei. For this study, we performed ten experiments with the cells from five different umbilical cords.
To analyze the micrographs with up to 50.000 cell signals per frame, we developed an automated image segmentation software. Using the programming language Python, we were able to extract the positions of the labeled cell nuclei, which were then connected for each cell to a trajectory as a function of time using the Crocker-Gier algorithm [
9]. Final trajectories (
Figure 1) of at least nine time steps were taken into account while the distance between consecutive positions must not exceed 10-12 µm, depending on the cell density. Cell velocities were calculated as the difference quotient of two consecutive cell positions.
The parameters and evidence of the models were estimated with the help of Bayesian data analysis (for details see [
4]). A numerical analysis was performed with the MultiNest-sampling algorithm [
10,
11] in its Python implementation [
12]. Uniform priors were applied for parameters
N0,
β,
t0,
A,
A1 and
A2. For all remaining parameters (
α0,
α1,
σ,
τ,
τ1,
τ2,
A0,
c and two parameters adjusting the experimental error), Jeffreys’ priors were used. The agreement of data and model was assessed using a Gaussian likelihood function. The model means and scatterings were obtained by averaging over the resulting parameter samplings.
3. Results
To obtain a first overview of the temporal development of cell movements, we calculated the mean squared velocity for each time step as the average of all cells that could be tracked for the entire observation time (
Figure 2).
The cells decrease their mean squared velocity over time (
Figure 2). However, this decline in velocities varies strongly between the ten experiments. We also noted different cell numbers at the beginning of the imaging as well as differing growth in the observed cell layers (
Figure 2, insets). With the cells moving around and proliferating vividly, the cell arrangement usually gets more homogenous and macroscopically ordered over time. The resulting overall higher cell densities confine cell movements and therefore decrease mean squared velocities. To highlight this inverse relation, we chose the same axis limits for the velocities and densities for all experiments in
Figure 2. If cell densities are high, then the corresponding velocities are, in comparison with the other experiments, quite low, and vice versa (e.g.,
Figure 2 a as opposed to c).
To develop a model describing the decline in mean cell velocities, the effect of several simultaneously ongoing biological processes should be quantified separately. First, we found that the cells rapidly slowed down at the start of the experiment. This is followed by a second, slower decline before the mean squared velocity seems to converge to a certain limit greater than zero for long times depending on the cell density. In some experiments (
Figure 2c,d) the two declines seem to happen sequentially, whereas there are other examples (
Figure 2h) that show a much smoother transition. We speculate that these two processes actually occur simultaneously with differing parameters since this seems to be more likely from a biological view. Hence, the little plateau, as observed especially in experiments 3 and 4 (
Figure 2c,d) after 10 h, could be caused by a third phenomenon. Based on this consideration, we described the general decline using two exponential functions, with parameters
A1,2 and
τ1,2 and a basal squared velocity
A0.
As shown in
Figure 2, this fit, obtained using Bayesian analysis, well captures the effects of the two processes that are slowing down the cells with the time scales of
τ1 and
τ2, ranging between ~16 and ~142 min and between ~655 and ~1780 min, respectively. These declines could be caused by the continuous cell aging and strengthening of cell-cell contacts. However, there are still data points that deviate from the model of Equation (1) during time ranges mainly within the first 20 h. Up to now, we did not include cell division activities in our model, which should be related to the velocities of cells: On the one hand, cells can move around more freely in areas with low density. They are less prone to collide with each other or have more space to change their directions. On the other hand, temporal phases of very active cell proliferation seem to be marked by a transiently constant or only very slowly decreasing cell velocity (see e.g.,
Figure 2c,d,i). Although this effect is short-lived, it shows that the separation of dividing cells can lead to accelerations in the system. By contrast, these cell divisions can also cause an overall deceleration due to an increasing confinement in densely packed cell clusters. That is why the mean squared velocities in experiments 3 and 4 show a transient plateau (
Figure 2c,d) in their decline or why
τ2, as a means to describe the second, slower decline, might differ in these experiments.
In order to quantify the connection between mean cell velocity and cell density, we calculated the relative changes in cell numbers for each time step (
Figure 3).
The cells usually go through one phase of more intense cell division, which builds up at first and reaches a maximum before it subsides again. However, the onset, duration and strength of this phase differ strongly between the experiments. Occasionally, the maximum of cell division lies probably even before the start of the experimental observation. We can describe the relative increase
α(
t) in the cell numbers
N per time
t with
Here, α0 describes the basal cell division, α1 the strength, σ the duration and t0 the time point of the maximum of the cell division phase. The exponent β changes the shape of the curve with β = 2 being the special case of a Gaussian function.
Since changes in cell density can modify the velocities of cells, we couple the cell division activity of Equation (2) linearly with factor
c to the mean squared velocity of Equation (1). This temporarily raises the mean squared velocity during strong increases in cell numbers for
c > 0.
Depending on the time point of the maximum cell division rate, the cell density can be calculated by integrating Equation (2), leading to the following result:
Here, the lower incomplete gamma function is used according to Equation (5)
and
N0 denotes the initial cell number at
t = 0. The obtained parameters of the simultaneous fit, calculated individually for each experiment, indicate that the trend of the velocities consists of a faster and a slower decline (
Figure 4). After about one to three hours (elapsed time equals two
τ1), the cells have slowed down considerably (
Figure 4a). After that time, the mean squared velocity decreases only slowly, with
τ2 usually being greater than 10 h (
Figure 4b). The exception from these findings is the ninth experiment, where the two declines are merged and therefore difficult to separate. The first two experiments show a particularly slow decline since the cell density had always been high and the cells from this umbilical cord showed less activity. At large, the obtained values for the parameters are roughly of the same magnitude compared to the counterparts of the isolated fit (Equation (1) and
Figure 2). This indicates that there are at least two time-dependent phases, with the cell division causing the time-independent alterations of the overall dynamics. This third influence is especially strong in experiments 3, 4 and 9, as indicated by parameter
c (
Figure 4c), where the mean squared velocities were also the highest. The first experiment shows slightly more cell division than the second while the mean squared velocities are almost the same (
Figure 2a,b). Interestingly, the optimal parameters for the combined fit include a higher
c and subsequently higher
τ1 and
τ2 (
Figure 4) instead of these three parameters just being smaller. This might emphasize the importance of considering cell division in the model.
Based on the estimated parameters, we plotted the fits for the mean squared velocity and density (
Figure 5).
A comparison by eye shows that the model is suitable to describe the data of the cell densities as long as the cell numbers do not decrease (
Figure 5a,b). Furthermore, the fit of the mean squared velocities was improved by the coupling to cell division activities, especially in experiments 3, 4 and 9. To verify that indeed two time scales are necessary to describe the overall decline, we compared the resulting Bayesian evidences with two alternative models. Model 1 (12 parameters) uses one exponential, Model 2 (13 parameters) uses one stretched exponential function in addition to the cell division term (
Table 1). So far, Model 3 (14 parameters) was applied (
Figure 5), which uses, in comparison with Model 1, an additional exponential function (2 in total).
The comparison of the evidences of the applied models shows a clear preference for Model 3 for nine out of ten experiments (
Figure 6). This endorses the theory that there is usually a decline consisting of two parts with different temporal dynamics. Only the ninth experiment could also be described with just one exponential function. These cells showed the strongest and longest-lasting cell division phase out of all experiments, which led to a more than two-fold increase in cell density. This influence was incorporated into the fit for the velocity via a high coupling parameter
c (
Figure 4c). Thus, the effect of the proliferation-independent cell movements is less recognizable, especially because the processes occur largely simultaneously. This might have also caused the small difference between
τ1 and
τ2 and its higher error (
Figure 4a,b).