Migration of individual microvessel endothelial cells: stochastic model and parameter measurement

Cell migration, including directed taxis, is a phenomenon exhibited by many cell types in numerous physiological processes. The need for rigorous, quantitative evaluation of motility is clear. Only quantitative assessment will provide the basis for objective comparisons between motile cell types and migration under various environmental conditions, as well as provide a rational framework for elucidating the underlying biochemical mechanisms of cell locomotion. Dunn and Brown (1987) have previously discussed at great length this need for quantitative evaluation of motility, as well as the desired features of measured quantities. The most appropriate measurements for the evaluation of motility and/or taxis in particular situations may vary, given the wide spectrum of cell types and physiological contexts. In this work we address those systems in which individual cell trajectories appear to be the descriptive level of greatest relevance. One such instance is the encounter and destruction of foreign targets by white blood cells in the immune and inflammatory host defense responses. Another is the development of new microvascular networks by microvessel endothelial cells (MEG) during angiogenesis.

In this paper we propose a stochastic mathematical model for random motility and chemotaxis of single cells. This model allows simulation of individual cell migration paths, and we have used it elsewhere to analyze the role of migration of MEC in angiogenesis (Stokes and Lauffen-burger, 1991). Our approach to modeling and simulation of cell trajectories is through phenomenological analogy to Brownian motion of molecules and inert particles. Clearly, the fundamental mechanisms by which cells move over two-dimensional surfaces or through three-dimensional matrices are radically different from the thermal motion of molecules. Nevertheless, observation of the motility of individual cells reveals comparable random walk-like behavior, indicating a similar stochastic nature and suggesting that a related mathematical description might suffice. Correspondingly, Dunn and Brown (1987) have proposed using the Omstein-Uhlenbeck (O-U) process, traditionally used to describe Brownian motion (Uhlen-beck and Omstein, 1930; Doob, 1942), to quantify individual cell random motility. The O-U process describes the velocity of a particle in terms of two processes: stochastic fluctuations in velocity (speed or direction, or both), and deterministic resistance to the current velocity. The stochastic term encompasses all of the probabilistic processes that might affect cell velocity (e.g. random fluctuations in motile sensing and response mechanisms, thermal noise). The deterministic term represents the decay of the current velocity, providing susceptibility to the random fluctuations and affecting persistence of the motion. Dunn and Brown utilized a discrete description of the O-U process to analyze the motility of chick heart fibroblasts. Using statistical analysis of the autocorrelation functions of cell displacements during discrete time intervals, they demonstrated that this process successfully characterizes fibroblast random motility. Their results support the use of the O-U process as a model of cell motility. In the present work, we extend the approach of Drum and Brown: first by using a continuous version of the O-U process; and second, by adding a directional bias term for cell movement in a chemoattractant gradient. We apply our extended model to analyze experimental data of MEC migrating in response to uniform concentrations and concentration gradients of the angiogenic stimulus acidic fibroblast growth factor (aFGF).

We begin by describing the continuous version of the 0-U process and our extension that provides a description of chemotaxis. Analysis of the statistical properties of this OU process reveals that the two parameters characterizing random fluctuations and velocity decay (α and β, respectively) are explicitly related to two commonly used and easily interpretable average cell movement quantities - speed, S, and persistence time in velocity, P_v. Speed and persistence time have also been defined in earlier probabilistic models for cell motility (Alt, 1980; Dunn, 1983; Othmer et al. 1988; Rivero et al. 1989). We validate our model’s applicability to MEC random migration, and measure α and β (and consequently obtain S and P_v) for migration in several chemical environments, including aFGF. An additional parameter, the chemotactic responsiveness, x, is introduced to account for directed motion of the cell in response to a chemical concentration gradient, x could not be measured directly for individual MEC because detectable aFGF concentration gradients could not be sustained over the length of time (several days) necessary to obtain statistically meaningful cell path lengths for these slowly moving cells. However, we have previously demonstrated using the under-agarose assay, that MEC populations exhibit a significant chemotactic response to aFGF and we have measured the population chemotaxis coefficient, χ (Stokes et al. 1990). We show here that we can estimate a value of the chemotactic responsiveness for individual cells, k, from the population chemotaxis coefficient, χ using model-generated cell paths. What our model implies about the nature of cell motility is addressed in the Discussion.

Human subcutaneous adipose tissue was obtained from human donors according to our Institutional Review Board protocol. The cells were isolated using the procedure of Wagner and Matthews (1975) as modified by Williams (1987). The cultures were incubated at 37 °C in a humidified, 5 % CO2 atmosphere. Cultures were fed three times per week with medium consisting of medium 199 with Earle’s salts (M199; Gibco, Grand Island, New York), 10 % heat-inactivated fetal bovine serum (FBS, Hazelton, Denver, Pennsylvania), 1.5 mM L-glutamine, 20 μgml^-1 of unpurified aFGF, and 90 μg ml^-1 heparin (porcine, Sigma; Thornton et al. 1983) (culture medium), and split at confluence using 0.25 % (w/v) trypsin with 0.09% ethylenediamine tetraacetic acid in normal saline. Cells used during the assays did not exceed two passages.

Cell cultures were analyzed for the presence and distribution of factor VUI-related antigen using the materials and methods of the Dako Peroxidase-Antiperoxidase Kit (Santa Barbara, California). Cultures that expressed factor VUI-related antigen were used in the studies.

Tissue culture polystyrene Petri dishes (35 mm; Coming Glass Works, Coming, New York) were pretreated with a 1% gelatin solution in 0.9 % NaCl overnight, and washed with Medium 199. A confluent culture of MEC was removed from a T25 culture flask with trypsin as above, and resuspended in a final volume of 2 ml of M199 containing 10% FBS. One to two drops (∼0.2ml) of the suspension were added to 4 ml assay medium and plated in a pretreated Petri dish. Assay medium was made by diluting 10 x concentrated M199 with Earle’s salts (does not contain NaHCO_a) to lx with distilled water and adding 0.1% BSA and 10mM Hepes. pH was corrected to 7.4 with NaOH, and the osmolarity was corrected by replacing the missing NaHCO₃ with NaCl. The cell density was approximately 400 cells cm^-2. The cells were allowed to settle and attach to the dish surface for at least two hours prior to videotaping.

Pre-formed agarose overlays were used in some experiments to create the same environment as we have previously used in the under-agarose assay for population migration measurements (Rupnick et al. 1988; Stokes et al. 1990). Overlays were formed by pouring a warm (40°C) solution of 0.5 % agarose (Sigma) in assay medium into an untreated Petri dish and allowing it to gel. After the MEC in a second dish were attached, the medium in the dish was aspirated until only a very thin film of fluid remained on the cells. The gelled agarose block was carefully removed from the first Petri dish by lifting it with a Pasteur pipet tip and laid onto the cells in the second dish. Several drops of medium were placed onto the agarose to prevent drying and to fill in around the edges. The use of agarose overlays was discontinued because the top of the agarose became pitted with time, marring the video image and making image analysis difficult.

The medium or agarose in a Petri dish to be used was covered with 2 ml of light mineral oil (Sigma) to protect the assay from environmental contaminants, while allowing for the flow of air across the dish and not impeding the optics of the microscope. The Petri dish was placed in a warmed Leiden chamber (Ince et al. 1983), a type of microscope stage incubator with a temperature controller, and the temperature was maintained at 37 °C in the dish. Filtered compressed air was blown across the surface of the liquid through the ports in the chamber to help maintain a constant temperature across the radius by convection.

The cells were observed with phase-contrast optics on a Nikon Diaphot inverted microscope (Nikon Inc., Garden City, New York). A video camera (DAGE-MTI Inc., Michigan City, Indiana) was attached to the microscope, leading to a time-lapse video recorder (JVC model BR-9000U). A single field was taped for up to 72 h using time-lapse at 1/240 real time; 15-20 cells per field were present.

We analyze the cell motility in terms of a stochastic mathematical model. As noted in the Introduction, the model is the Omstein-Uhlenbeck process for random motility, modified to account for the directional bias caused by sensitivity to chemoattractant gradients. The O-U process is the simplest type of continuous, autocorrelated, random motion process. It describes a Markov process, meaning that the motion at a given time depends only on the current state of the cell and not on earlier states. That is, there is no memory to the motion, and later steps are not correlated to earlier steps. While this seems very simple, Dunn and Brown (1987) showed that, at least for fibroblasts, this simple Markov process is sufficient to describe cell random motility. We chose the O-U process because it is the simplest stochastic model that can describe the random motility process as a persistent random walk, and there was no reason a priori to choose a more complicated stochastic equation. The implications of this model for the nature of cell motility are addressed further in the Discussion. We expand our random motility model, the O-U process, to account also for the directional bias in velocity caused by a cell reacting to the presence of a chemoattractant gradient. The random walk now becomes a biased random walk. In this model of directional bias, we assume that the cell effectively responds to a spatial gradient of bound receptors across its surface, although we make no statement about the cell’s actual underlying mechanism of perception (see, for example, Lauffenburger et al. 1988). We also assume that the effect of a chemical gradient on cell velocity change increases as the cell is directed further down the gradient. This form is based on previous findings on chemotaxis in leukocytes, though little data are available to indicate whether MEC behave in the same way.

We now describe the details of the mathematical model. One equation (equation (2)) describes the velocity (rate of movement and direction) of a single model cell as it evolves with time. A second equation (equation (3)) is then used to determine the cell’s position from its velocity, resulting in a calculated trajectory for the model cell. From these equations, the mean square displacement of a cell with the model characteristics can be calculated (equations (5) and (6)). We compare the predicted mean square displacement from the model with that measured from MEC cell paths as explained above to: (1) examine the validity of the model, and (2) determine values of the model parameters for the MEC migration.

The geometry of our system is illustrated in Fig. 1. The rate of change of cell velocity in our model is the sum of three components: the first two components represent the random motility behavior, and the third represents the chemotactic, or directed, behavior. For random motility, the first term describes the deterministic resistance to the cell’s motion, and the second represents the random accelerations or fluctuations of the movement. The third component that describes the chemotactic behavior provides a directional bias in the presence of an attractant gradient. Mathematically, these three components are summed to give the following stochastic differential equation for the rate of change of the cell velocity, ⁠:

Carrying out this integration will yield an actual cell path for the cell migrating with the velocity given by equation (2) (e.g. see Fig. 4, below).

We can see from equation (2) that the velocity process for random motility (when ⁠) is characterized by two parameters: a represents the magnitude of the random fluctuations, and β represents the rate constant for decay of the current velocity. We can view β as the magnitude of the resistance to the motion, or the magnitude of susceptibility of the motion to the random fluctuations, leading to a decrease of persistence in velocity as β increases in value. We emphasize that for cell movement these fluctuations are not the same as the thermal fluctuations leading to molecular Brownian motion, but rather are due to intrinsic cellular mechanisms such as receptor-binding fluctuations (Tranquillo and Lauffenburger, 1987).

We define K as the chemotactic responsiveness. It is related, though probably not in a trivial way, to the cell population chemotaxis coefficient, χ (Tranquillo et al. 1988; Farrell et al. 1990; Stokes eí al. 1990).

Equations (2) and (4) provide a description of the velocity process of a single cell moving with the average properties represented by the random motility parameters α and β and the chemotaxis parameter K. Each cell described by equation (2) with a given set of parameters, α, βand K, has a different path or trajectory (owing to different realizations of the random noise process, W) but the same ‘degree’ of randomness, velocity decay rate constant, and sensitivity to attractant gradients. Thus, our modified O-U process provides a quantitative framework in which to interpret cell motility. For this framework to be meaningful, we must first demonstrate that it is a valid description of actual cell trajectories. If shown, we can utilize this framework to obtain values for the parameters α, β and K, which then provide a quantitative description of the characteristics of the actual cell’s intrinsic motility properties.

This states that a migrating cell that can be described by the parameter values α and β, and moving for a duration T, is expected to be displaced by the distance D. By the Pythagorean theorem, the right-hand side of equation (5) is multiplied by 2 or 3 for migration in two or three dimensions, respectively. We compare the model to actual cell migration by making experimental measurements of ⟨D²⟩ versus T as explained above, and comparing the experimental curve to that given by equation (5) with various values of a and fl. In practice, nonlinear regression using the IMSL subroutine ZXSSQ is used to match equation (5) to the data and thereby obtain the values of the parameters that best describe the data.

To relate these fundamental parameters (α and β) to previously defined and perhaps more intuitive quantities, we define two new parameters, the speed, S, and persistence time, P_v. These are not new, independent parameters, but as a pair contain the same information as the pair α and β. The speed of the cell is a measure of the magnitude of the velocity, or its rate of movement. S contains no directional information. The persistence time is a measure of persistence m velocity, or the average time a cell maintains a given velocity. P_v contains information on both direction and rate of movement. Doob (1942) calculated that the expected square speed for random motility (no chemotaxis, ⁠) for the velocity process in equation (2): is⟨s²⟩ =nα/2β, where n is the dimension of the motion.

The most obvious definition of a characteristic speed, S, is thus the root mean square speed (in two dimensions, n=2). The appropriate definition of persistence time is not as obvious, but can be deduced by examining the meaning of the parameter fl. To reiterate, β is the decay rate constant of the current velocity. Small β means little resistance to the motion or little susceptibility to random fluctuations, which will result in greater time spent at the same velocity. Conversely, large β means great resistance to the motion and hence less persistence. β has dimensions of time^-1, and so its inverse, 1/β, should be a measure of the time of persistence in velocity, P_v. This definition of P_v as persistence in motion (both direction and magnitude of the velocity) is similar to the persistence time defined by Rivero et al. (1989) in their cell population model derivation. Othmer et al. (1988), Alt (1980) and Dunn (1983) all defined persistence time as persistence in direction only, with definitions of cell speed similar to ours.

Though S and P_v are derived quantities based on the parameters α and β we will report results primarily in terms of S and P_v because of their more appealing definitions. We emphasize, however, that the fundamental parameters defined in the model are those giving the magnitude of fluctuations and velocity decay rate or resistance to motion; that is, α and β.

Values for the chemotactic responsiveness K were found using computer simulation of cell paths according to equations (2)-(4), and comparing measurements from these paths with measurements of MEC chemotaxis in the under-agarose assay reported elsewhere (Stokes et al. 1990). More precisely, we determined values of the parameter grouping ka_o/(α and β), which contains the chemotactic responsiveness K, and a_o is the attractant concentration at the attractant source (see Fig. 1). We call this grouping δ, the dimensionless chemotactic responsiveness. This grouping emerges from dimensional analysis of the equations as shown in the Appendix.

To determine values of the dimensionless chemotactic responsiveness, δ, statistics of the simulations were compared with experimental results from Stokes et al. (1990) for MEC migration in gradients of aFGF. First, theoretical cell paths were simulated for various values of chemotactic responsiveness. For a given path, this is done by numerically calculating a realization of velocity, ⁠, for a given realization of the noise process, ⁠, for some length of time using equations (2) and (4). This process is described in the Appendix. Then the position of the cell, i, is calculated using equation (3). Examples of simulated cell paths are shown in Fig. 4 (below).

Sets of 100 cell paths were calculated for δ=0 (no chemotactic response), 0.1, 0.3, 1.0 and 3.0 (increasing responses) using different realizations of the random noise for each path. The chemotactic index, CI, was calculated for each path. CI is most simply defined as the net distance travelled towards an attractant source or target divided by the total distance travelled (McCutcheon, 1946). It can be written approximately either in terms of cell velocities: (where is the chemotactic velocity, the component in the direction of attractant due to bias); or in terms of distances: CI=(net distance travelled up the gradient)/(total distance travelled). Othmer et al. (1988) have derived an expression for CI that takes into account the effect of finite persistence time compared to the total time of one’s measurement. The correction is small if the total time is much greater than persistence time, as was the case in our simulations, and thus this correction was ignored. Experimentally estimated values of CI can be calculated because can also be defined in terms of a previously measured parameter for chemotaxis, the chemotaxis coefficient χ. This is defined in the cell motility model of Alt (1980), and quantifies the net drift of a cell population due to directionally biased movement in an attractant gradient. Rivero et al. (1989) have derived the relation: ⁠, so CI=χ ▿a/S. This relation has been verified experimentally by Farrell et al. (1990) for alveolar macrophages responding to C5a. In our case, χ has been measured for migration of MEC in aFGF gradients (Stokes et al. 1990). To estimate δ, then, we compared the CI measured from many simulated cell paths to that calculated from the experimental values of χusing CI=χ ▿a/S. The value at which the predicted and calculated values of CI agree is the one for which ô predicts individual cell paths with the same directed component as that measured in the cell population experiments.

Several typical cell paths for cells migrating in assay medium are shown in Fig. 2. The paths reveal that MEC generally travel with smooth trajectories, and that they persist in a given direction for a time on the order of hours. The duration of all paths obtained ranged from 18 to 51 h. While a field was typically taped for 72 h, no cell paths could be traced for that long because the cells either moved off screen or collided with other cells. The occasional cell that did not move was ignored in this study.

These paths represent the movement of the centroid of the cell area. For these slowly moving cells in which direction changes are infrequent and smooth, the variable selected to represent cell position (area centroid, nuclear area centroid, dry mass centroid etc.) is not expected to be of much importance. For faster cells that display pronounced lamellipodia, however, this choice may be of more significance. For some cells such as leukocytes it is common to see a new protrusion of the cytoplasm in a new direction while the nucleus continues temporarily in the previous direction of cell movement. One’s choice of definition of position must therefore be carefully examined. For example, Brown and Dunn (1989) have compared the changes in spread area and total dry mass of migrating chick heart fibroblasts.

The predictions of the model were compared with experimental measurements using plots of mean square displacement versus time increment. In Fig. 3, the mean square displacement, ⟨D² ⟩, is plotted against time increment T for both the experimental data (continuous curve) and a curve calculated from equation (6) (broken curve) using the best S and P_v values obtained from nonlinear regression of the data. The two examples in the figure illustrate that MEC migration paths follow a function of the form of equation (6) very closely for small to intermediate values of T (compared to the full duration of the cell path). As seen there, the data can become unpredictable for large time increments for which ⟨D² ⟩ is the average of only a few points (7¹ approaching the full duration of the cell path). This is to be expected, because at longer time increments there are fewer and fewer displacement increments to average. For instance, at T equal to the full duration of the cell path, there is only one point, whereas for 7’=lh there are many increments. Because of this phenomenon, values of S and P_v reported here are calculated from the experimental data ⟨D² ⟩ vs T for T up to half the total duration of the path. All paths observed showed agreement between experimental data and the model similar to that in the examples in Fig. 3 for T up to or greater than this value.

a and f, and consequently S and P_v, were determined for MEC migrating in assay medium (no stimulants or inhibitors) both with and without agarose overlays. The results are given in Table 1 for both sets of parameters. The data consist of evaluations of 10 cell paths with overlays and 8 without. With the agarose overlay, the average cell speed is S=16.4±1.3 μmh^-1 (mean±standard error). Without the overlay, S=29.6±5.5 μm^-1, significantly different from the former value (P<0.02). In the latter group, one cell in particular (number 14) lies far outside the range of two standard deviations of the mean, with a speed of 63 μmh^-1. If cell 14 is removed from the average, S=24.9±2.8μmh^-1. This remains significantly different from S with the agarose overlay (P<0.01). The average persistence time for the migration with agarose overlays is P_V =6.0±1.0h. Without overlays, P_V =5.4±1.0h, not significantly different from the former group. With cell 14 removed, P_V =5.5±1.0h, again not significantly different. These data indicate that the agarose overlay hindered the speed of the cells, though it had no effect on the persistence time.

S and P_V were also evaluated for seven cells migrating in culture medium containing unpurified aFGF, heparin and FBS. These data are given in Table 2. For these cells, S=41.6±4.8/anh^-1 and P_v=2.9±l.lh. These values are significantly different than those for motility, both with overlays (for S, P<0.1; for P_V, P<0.1) and without (for S, P<0.02; for P_v, P<0.1). This result indicates that the rate of movement of the cells is increased while the persistence time is decreased, compared to control. If persistence time is loosely interpreted in terms of persistence in direction only (constant speed), then an interesting result of increasing one parameter (here S) while decreasing the other (P_V) is to maintain the mean free path of a cell (the length of a ‘straight’ segment) approximately constant. The mean free path equals the rate of movement times the time the cell moves in one direction, or SP_V.

We have provided the data for the individual cells in Tables 1 and 2 as well as the means in order to demonstrate the variability exhibited by MEC in both speed and persistence time when migrating in a given environment. The relative variation in persistence time within a group was somewhat greater than that for speed. Several factors contribute to the variability in addition to natural variation among cells. First, the length of time a cell path is observed limits the accuracy to which its speed and persistence time can be determined. The more multiples of persistence time a cell is observed, the more accurate will be the evaluation of S and P_v. The length of time a given cell can be observed is constrained by experimental factors including cell-cell collisions and cells moving out of the field of view. In addition, the greater variability in P_v at least partly results from the functional dependence on P_v in equation (6). The curves predicted from equation (6) are less sensitive to P_v than to S (that is, P_v can vary over a larger range than S with less change in the ⟨D² ⟩ versus T curve (Fig. 3)) so the values of P_V reported are more sensitive to tracking error than is S.

The values of the fundamental motility parameters α and βare also given in Tables 1 and 2 for each cell. The magnitude of random fluctuations, a, which is related to cell speed and not to persistence time (recall S= √ (α and β) and P_V=1/β), changed with the addition of the agarose overlay. The velocity decay rate constant, β, is apparently not affected (since P_v was constant). The addition of aFGF affects both the randomness and velocity decay rate, with both α and βincreasing.

More cell paths were tracked than appear in Tables 1 and 2, but the tracks were of too short a duration to analyze in terms of this mathematical model. The analysis to obtain S and P_v is only appropriate when the cell paths are measured for at least several persistence times. At the extreme, for instance, if the cell is moving at constant speed and is followed for less than a full persistence time then the track is essentially straight and the analysis would predict that P_v=∞. Clearly this is incorrect and is only obtained because no turns were made in that section of path. While not quite so drastic, if the path is followed for only a couple of persistence times, P_v will still be incorrectly large because the true tortuosity is not captured in the data. Therefore, in the data reported in Tables 1 and 2 only paths that are at least three times the duration of the predicted P_v (duration/P_v>3) were kept.

The value of the random motility coefficient, p, was calculated from equation (7) for each cell track as given in Tables 1 and 2. To reiterate, p is analogous to the diffusion coefficient, and quantifies the rate at which a cell population will disperse or spread. For migration in control medium, the mean value of μ is 2.3±0.6×l0^-9 cm²s^-1 with agarose overlays and 6.9±2.6×l0^-9cm²s^-1 without overlays (mean calculated from values of q, not from means of S and P_v). Again, if cell 14 is eliminated (it has μ=2.4×l0∼⁸cm²s^-1, about five times the mean), the mean without overlays is 4.6±0.9×l0^-9cm²s^-1. These values are similar to n values previously determined directly from population migration experiments in assay medium in the linear under-agarose assay. In that assay, the random motility coefficient was about 5×l0∼⁹cm²s^- for several MEC lines isolated from human subcutaneous fat (Stokes et al. 1990).

As shown in Table 2, while S increased and P_v decreased in culture medium compared to assay medium, the random motility coefficient increased: μ=7.1±2.7×l0^-9 cm²s^-1. While we have not measured μ directly for population migration in culture medium, we have done so for migration in aFGF with and without heparin in assay medium (Stokes et al. 1990). The calculated value of μ from S and P_v is very similar to that in aFGF at 10⁻¹⁰ to 10⁻⁸M with 1μgml^- heparin, where μ =7×l0^-9 to 9×l0^-9 cm² s^-1. These results indicate that the measurements of S and P_v for individual cells are in good agreement with earlier population measurements, even based upon the relatively small number of cell trajectories analyzed here. Further, the correspondence of the results from the two models reinforces the theoretical basis of the correspondence between the models for individual cell movement and movement of cell populations.

Estimation of the dimensionless chemotactic responsiveness, δ, was accomplished using computer simulation of cell paths according to equations (A3) and (A4) (see Appendix). Fig. 4 illustrates several cell paths for the same white noise realization, showing the changes that occur as the chemotactic responsiveness δ is increased exponentially from 0 to 3. The same noise realization was used for all paths, so the calculations were identical except for the value of δ. In order to obtain an average chemotactic index, CI, for a given value of <5, simulations were run for 100 cells for each δ, and the CI values calculated for the 100 cell paths were averaged. The initial direction of each cell was random. This was repeated three times, giving three sets of 100 cells each for each <5. These results are presented as a plot of CI versus log δ in Fig. 5. Superimposed on the plot are the values of CI generated analytically from the relationship CI=χ ▿a/S from the population migration experiments of Stokes et al. (1990). In those experiments, the maximum chemotactic response was measured when a_o=10⁻¹°M, with ▿a =3.5×l0^-16M μm^-1, which gave χ=2600 cm² s^-1. Using S’=40μmh^-1, which might be expected in the presence of an angiogenic stimulus, CI=χ ▿a/S#x223C;0.08 in those experiments in that gradient. The large open diamond in Fig. 5 represents that experimental datum. It is superimposed on the simulations curve at the place where CI=0.08. The gradient is substantially larger in our simulations, ▿a∼2×10⁻¹⁴M μm^-1. For this Va value and again assuming S=40 1mh^-1, the same chemotactic ability (same /) in this larger gradient should give a CI∼0.5 in the simulations. This is represented by the large open square in Fig. 5, superimposed on the simulations curve at the point where CI=0.5. Reading off the corresponding value for log δ, the simulations predict that values of δ∼1.5 could be expected for comparable chemotactic ability as in the aFGF experiments. Lower values of δ are expected for less than optimal attractant concentrations, and higher values are not expected, on the basis of the experimental measurements of χ. With S=40 μmh^-1, δ =1.5 corresponds to values of κ a₀ of about 2400μm²h^-2. There is no particular reason to specify a₀ in the simulations to obtain a value for κ itself, so we only report the product Ka₀.

A quantitative framework for analyzing the migration of individual cells has been proposed, based on the stochastic O-U process introduced for this problem by Dunn and Brown (1987). In addition to using the continuous version of the O-U process instead of the discrete version used by Dunn and Brown, we have extended this framework to include chemotaxis. In this model, random motility is described as a combination of random fluctuations (quantitatively characterized by their variance, α) and a deterministic decay rate of current velocity (quantitatively characterized by the parameter β). Definitions of speed, S, and persistence time in velocity, P_v, were based on the fundamental parameters α and β. We must emphasize that S and P_v contain the same motility information as the pair α and β, simply written in a different way. Thus, random motility is fully described by only two parameters, the fundamental parameters α and β, or, equivalently, the defined parameters S and P_v. A third term was added to the O-U process to account for the effects of a chemical gradient on the velocity of the cell. This term is postulated to be proportional to the deviation of current direction of cell movement from the steepest gradient and to the attractant gradient magnitude (see Fig. 1), and is characterized quantitatively by the chemotactic responsiveness parameter, κ.

As shown by the match of the O-U process to data for (D²) versus T (Fig. 3) for MEC random migration, the O-U process appears to be an appropriate description of this motility. We have pointed out that the agreement of the MEC motility measurements and the model predictions in terms of average properties is only the first level of comparison on which to verify the applicability of a model. The more detailed analysis of autocorrelations that Dunn and Brown (1987) conducted on their data lends support to the use of the O-U process to describe cell motility, however, and could be applied to our data as well. That work is beyond the scope of the present paper. The reason such further analysis is desirable is because other forms of stochastic equations can be formulated to describe cell motility and chemotaxis (e.g. see Tranquillo and Lauffenburger, 1987) and some method is required to understand which model is more representative of a given cell’s locomotion properties. The first level of distinguishing among models is the average properties such as (D²) as we have demonstrated here, and the next level is the statistical properties that can demonstrate whether sequential movement steps are correlated in time and how strongly. Use of autocorrelation structures also may be required to relate observable quantities (e.g. S and P_v) to motility mechanisms.

While quantities defined as speed and persistence time have been measured previously (Glasgow et al. 1989; Farrell et al. 1990; Vicker et al. 1986; Wilkinson et al. 1984; Zigmond et al. 1985), the under lying stochastic framework provided by the O-U process provides two advantages. First, S and P_v are defined in terms of quantities (cr and ft) that are meant to describe the mechanisms of the motility. While the actual mechanical or biochemical processes that a and fl should represent are not presently defined, these parameters provide a direct means of relating understandable and observable quantities (S and P_v) to basic motility mechanisms. Possible mechanistic interpretations of a and fl are discussed below. Second, having the underlying framework that describes the velocity process rather than deriving S and P_v from probabilities of moving and turning (e.g. see Rivero et al. 1989) allows us to simulate cell paths with the model that have the same average properties as those measured experimentally. How to do so is demonstrated in the Appendix in our estimation of the chemotactic responsiveness κ (Fig. 4). This is very useful if we desire to investigate the roles of motility in various physiological processes. In fact, we have used this migration model with the values of α, β and κ reported herein in a simulation model of network formation in angiogenesis to test our hypotheses about the roles of MEC motility and chemotaxis in that process (Stokes and Lauffenburger, 1991). It has also been used to analyze the effects of macrophage chemotaxis on cell/target encounter rates (Charnick et al. 1991). Measurements of parameters defined in previous cell population and individual cell models do not provide this ability.

Our results for MEC random migration reveal the effects of several different environmental variables. Without added chemical stimuli, MEC exhibited different rates of movement when migrating with and without an agarose overlay. Persistence time, in contrast, was the same in both situations. In terms of the stochastic model parameters, the magnitude of accelerations (α) was greatly reduced by the overlay while the velocity decay rate (β) was unaffected. These results suggest that the agarose overlay created a more difficult environment in which to generate effective locomotory accelerations.

Results for MEC motility in culture medium were somewhat surprising. Compared to the control, speed increased by about 50 % while persistence time decreased by 50 %. This means that the MEC moved more quickly but also changed direction more often. This is the first time that a simultaneous increase in speed has been accompanied by a decrease in persistence time in the reports of several cell types in which these parameters have been measured (Farrell et al. 1990; Glasgow et al. 1989; Vicker et al. 1986; Wilkinson et al. 1984; Zigmond et al. 1985). One consequence of this result may be interpreted using the population random motility coefficient, u. We showed that the value of u increased in culture medium compared to control. Since u is a measure of dispersion or ‘area coverage’ of a migrating cell population, the result indicates increased area coverage or spreading of a population in the culture medium. This predicts the faster growth of new vascular networks in such conditions compared to that in control conditions. At the same time, however, the opposing variations of α and β should yield distinct network morphologies for stimulated versus control conditions (see Stokes and Lauffenburger, 1991), demonstrating the importance of the stochastic cell path parameters.

For chemotaxis, the comparison of simulations to experiment predict that the dimensionless chemotactic responsiveness δ (=κ A₀/S²) is about 1.5; ô can be thought of as a ratio of the chemotactic component of a cell’s migration to its random motility component. With the addition of a chemotactic factor, any increase in speed would require a proportionally squared increase in xu,, to provide the same relative contribution to the resulting movement from chemotaxis. This predicts that increases in cell speed without accompanying increases in chemotactic responsiveness are detrimental to the accomplishment of covering distance up a chemical gradient. We must note that we have not, to this point, demonstrated that the functional form of the chemotactic component in equation (4) is a valid description of MEC chemotactic directional bias in gradients, but use it as a first approximation.

The relationships between (α/β) and (S, P_V) infer what the velocity equation might represent in terms of the mechanism of cell migration. Again, S and P_v are simply combinations of α and βcontain the same motility information in a different form. For random motility (when κ=0), equation (2) indicates that certain experimental treatments might affect the parameters a (degree of randomness) and β (velocity decay rate constant) independently of the other. In terms of observed speed (S=√α/β) and persistence time (P_v=l/β), however, some treatments that affect αshould only affect cell speed, but treatments that affect (3 should affect both speed and persistence time directly. All of data reported here are consistent with this prediction. Our data for MEC migrating in control (no agarose) and culture media indicate that the addition of MEC growth factors increased both α and β. This means that while the decay rate of current velocity increased (β increased), the random fluctuations also got much larger (a increased). The latter seems intuitively reasonable, since several chemical factors (αFGF, serum, heparin) were added that could stimulate the cells through receptor binding, adding mechanistic steps at which random fluctuations might occur. The decreased persistence time might then represent the faster rate of turnover of motile machinery or transfer of signals throughout the cell caused by its ‘stimulated state.’ Comparing control experiments both with and without agarose, only a changed, decreasing in the presence of agarose overlay. This decrease indicates less randomness at some step in the stimulus sensing or response machinery. In other studies, Wilkinson et al. (1984) have characterized a treatment of neutrophils that affects cell speed but not the persistence time. Farrell et al. (1990) have shown that certain concentrations of the complement factor C5a increase both the speed and persistence time of alveolar macrophages. All of these measurements are consistent with the model, where both S and P_v are affected (α and βchanging) or only S is affected (only α changing). We are not aware of any treatments of any cell type in which the persistence time alone was affected by a chemical treatment, which would contradict our model.

We can consider the molecular processes that might be represented by α and β. The small and rapidly changing fluctuations or random kicks represented by √ αdW(t) accumulate to determine the current motile state, or speed and direction, of a cell described by the O-U process, α determines the strength of these fluctuations, while β measures how quickly their effect on the motility fades away. Dunn and Brown (1987) suggested that the accumulation of the random fluctuations represents the changes in the assembly of some component(s) of the locomotory machinery. The parameter a would then represent the strength of fluctuation while βmight represent the rate of turnover of the assembled structures. Alternatively or additionally, α could represent any or all random fluctuations in the multi-step molecular mechanism determining cell motility. This could include fluctuations in rates of receptor-attractant binding on different areas of the cell surface (stimulus sensing), fluctuations in local concentration of attractant (stimulus), or fluctuations in transmission of signal between cell-surface binding events and cytoskeletal motility machinery events (stimulus transduction). For example, Tranquillo and Lauffenburger (1987) have analyzed the effects of random fluctuations in local attractant concentration on cell motility behavior, β would then be the integrative process by which such fluctuations are damped out. Our results suggest that migration in a confined space (under agarose) decreases the effects of random fluctuations, while addition of growth factors increases randomness and the rate at which such randomness is dissipated within the cell.

The chemotactic responsiveness κ can be interpreted similarly, κ represents the ability of the cell to respond to a perceived concentration gradient in its environment. The assumption inherent in the model is that cells can effectively sense spatial gradients (i.e. temporally changing chemical concentrations are not necessary, though the cell’s sensing mechanism may involve temporal aspects). Thus, κ represents the ability of the cell to respond to a differential signal across its length caused by a gradient of bound cell-surface receptors. In a more sophisticated treatment, κ would not be a constant but would be some function of the gradient of bound receptors and the attractant concentration, to account for down-regulation and saturation of cell surface receptors at high concentrations (similar to the treatment by Tranquillo et al. (1988) of the chemotaxis coefficient, χ). This was not attempted in this work because such functions are poorly quantified in the endothelial cell in relation to cell motility at this time.

To understand better the functional significance of cell motility in tissue function, we can compare the values of S and P_v for MEC with those for other cell types. Glasgow et al. (1988) and Farrell et al. (1990) have shown that typical speeds and persistence times for guinea pig and rat alveolar macrophages in response to chemotactic peptides and C5a are 120 to 150 μmh^-1 and 0.3 to 0.5 h, respectively. These compare to our measurements of S=20 to 40 μmh^-1 and P_v∼3 to 6h for MEC. In population migration experiments (and as calculated from speed and persistence time), the random motility coefficients (/i) measured for alveolar macrophages are around 10⁻⁸ cm² s^-1 (Glasgow et al. 1988; Farrell et al. 1990), very similar to those for MEC (Rupnick et al. 1988; Stokes et al. 1990). Therefore, though the rate of dispersion of the macrophages as characterized by u. is approximately the same as that of MEC, the macrophage moves much quicker than the MEC but also makes many more turns. These differences appear to have a functional basis: the function of the alveolar macrophage is to patrol the surface of the lungs, scavenging for foreign particles, while the moltility of endothelial cells is involved in the development of new capillaries and in the re-endothelialization of denuded areas of larger vessels. For the macrophage, the ability to cover a lot of area quite densely is probably of more use in locating foreign particles, whereas the endothelial cell involved in forming new blood vessels should not turn too often or reasonable networks with useful lengths of vessels might not be formed. For another white blood cell, the neutrophil, speed and persistence time are around 1200 μm h^-1 and 0.25 h, respectively (Zigmond et al. 1985). Their random motility coefficient is in the range of 10⁻⁷ to 2×l0⁶cm²s^-1 as measured in the under-agarose assay in response to chemotactic peptides (Tranqüillo et al. 1988). This is one or two orders of magnitude higher than that of MEC and alveolar macrophages, indicating that the rate of dispersion of neutrophils is significantly faster than that of MEC and macrophages. This difference is likely to be important in the inflammatory response, reflecting the time in which neutrophils must respond to infectious agents in order to contain them. The ability to make these comparisons and infer cell-function relationships is possible because of the quantitative nature of the parameters reported. In particular, these parameters represent the intrinsic motility behavior of the various cell types, and have the advantage of not depending on extraneous aspects of the assay system used to measure them.

In summary, we have proposed a mathematical model of single cell random motility and chemotaxis that provides: (1) a quantitative framework in which to assess these cell functions, and (2) the ability mathematically to simulate cell paths for use in mathematical models of physiological processes in which cell migration is involved. We have demonstrated that the continuous O-U process is an appropriate model for MEC random motility at the level of average properties, and have demonstrated how the chemotactic sensitivity can be estimated from measurements of population migration. The use of quantitative models to evaluate cell migration should provide an objective basis upon which to evaluate mechanistic relationships among underlying physiochemical properties, emergent motility properties, and the functional significance of the motility in tissue organization and function.

We thank Paul DiMilla for development work on the image analysis system, Ralph Nossal for helpful discussions, and gratefully acknowledge funding from NIH grant GM41476.

In order to simulate cell paths under chemotaxis conditions, an attractant gradient, ▿a, must be specified. We used one plane of a three-dimensional spherical steadystate gradient created by a sphere (disk in two dimensions) at concentration a₀ with radius R (see Fig. 1, main text). The attractant concentration, a(r,t), and the attractant concentration gradient, ▿a(r, t), for this situation are Va(r, t)=a_oR/r², and Va(r, t)=−a_oR/r², for r>R, with the distance r measured from the source center. We use this gradient for two reasons. First, it is plausibly relevant to in vivo tissue angiogenesis situations (such as the rabbit cornea assay (Auerbach et al. 1974) and the chick chorioallantoic membrane assay (Gimbrone et al. 1974)). Second, it yields a very simple and easily interpretable form of the dimensionless chemotactic responsiveness, Ô, defined below. Any other gradient could have been chosen.

W indicates the Weiner process in terms of the dimensionless time (i.e. mean=0, variance =| τ|). Comparison of equation (Al) with equation (2) reveals that a three-parameter model (α, β and κ) has been reduced to a one-parameter model (δ). τis the only dimensionless parameter that emerged, and will be defined as the dimensionless chemotactic responsiveness. <5 can be thought of as a comparison of the chemotactic bias or drift tendency, KU_O, with the random or dispersive tendency of the migration, S².

Ñ (0,/i) indicates a Gaussian distributed random variable with mean of zero and variance h, where h is equal to the dimensionless time step size.

The method of solution is illustrated in Fig. Al. The necessary time step size for convergence of the solution was determined following the procedure of Wright (1974). At the ith time step, the deterministic terms of equation (A4) (the first and third terms) are calculated. Then the value of the random noise, 7(0,A)” is added to the deterministic part. The new position of the cell is subsequently calculated from the velocity according to equation (A4). This process is successively repeated to obtain a realization of and for a given realization of the noise process Ñ (0,h), as illustrated in Fig. Al.