Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Cell Mol Bioeng. Author manuscript; available in PMC 2010 April 27.
Published in final edited form as:
Cell Mol Bioeng. 2010 March 1; 3(1): 3–19.
doi:  10.1007/s12195-010-0114-2
PMCID: PMC2860333

Application of Population Dynamics to Study Heterotypic Cell Aggregations in the Near-Wall Region of a Shear Flow


Our research focused on the polymorphonuclear neutrophils (PMNs) tethering to the vascular endothelial cells (EC) and the subsequent melanoma cell emboli formation in a shear flow, an important process of tumor cell extravasation from the circulation during metastasis. We applied population balance model based on Smoluchowski coagulation equation to study the heterotypic aggregation between PMNs and melanoma cells in the near-wall region of an in vitro parallel-plate flow chamber, which simulates in vivo cell-substrate adhesion from the vasculatures by combining mathematical modeling and numerical simulations with experimental observations. To the best of our knowledge, a multiscale near-wall aggregation model was developed, for the first time, which incorporated the effects of both cell deformation and general ratios of heterotypic cells on the cell aggregation process. Quantitative agreement was found between numerical predictions and in vitro experiments. The effects of factors, including: intrinsic binding molecule properties, near-wall heterotypic cell concentrations, and cell deformations on the coagulation process, are discussed. Several parameter identification approaches are proposed and validated which, in turn, demonstrate the importance of the reaction coefficient and the critical bond number on the aggregation process.

Keywords: Melanoma cells, Leukocytes, Endothelium, Collisions, Aggregates, Probabilities, Cell adhesion, Shear conditions


Malignant melanoma, as the most deadly form of skin cancer, has a highly metastatic nature. No effective treatment has been found so far to prevent the metastasis process. Polymorphonuclear neutrophils (PMNs) were suggested to be a major contributor that facilitated tumor cell extravasations from the circulation during metastasis.39,40,42,46,47 Several recent studies have also explored how the fluid dynamics in shear flows would affect interactions between PMNs and melanoma cells.26,27 Still, some aspects of the process remained unclear, e.g., how heterotypic cell populations and the hydrodynamic environment would affect receptor–ligand mediated adhesion kinetics between PMNs and melanoma cells, although studies on non-statistically populated PMN–melanoma interactions have recently been reported.8,14,2325 At the microscopic level, PMN-facilitated melanoma cell adhesion to the endothelial cell (EC) monolayer under shear has been proposed as a two-step process,2527 which involves initial PMNs tethering on the EC and subsequent melanoma cells being captured by the tethered PMNs. Thus, a special focus of our study is on how shear conditions affect two populated groups, PMNs and melanoma cells, in terms of PMN–melanoma cell collisions and subsequent aggregations in the near-wall region.

The cell population dynamics under various flow conditions has been a subject of research for decades. Experimental studies employed the cone-plate viscometer assay or the parallel-plate flow chamber assay to understand how fluid dynamics affected cell–cell collision and coagulation.1,4,5,21,31,37,38,44,48,49 Computational fluid dynamics (CFD) simulation and statistical modeling also contributed remarkably to the understanding of the aggregation process.3,1416,21,25,29 Most of the researches focused on cell–cell interactions in a linear uniform shear flow. For example, Laurezi and Diamond21 addressed the platelet and PMN heterotypic aggregation in the cone-plate linear-shear assay by Monte Carlo simulation. Their approach was more recently extended to model the dynamics of PMN–melanoma cell aggregations in the parallel-plate flow chamber.45 While the approach was shown to have a great deal of potential, these studies remained limited and a number of issues need to be addressed further. In particular, it is important to consider the effects of cell deformation, known to play significant roles in the collision and adhesion process.8,9,14 More systematic parameter identifications should be performed instead of an ad hoc approach due to the lack of sufficient experimental data.45 Furthermore, it has been well understood that for the more physiologically relevant case, the ratio between the concentration of PMNs (CP) and that of melanoma cells (CT) is much greater than 1, however, most mathematical models, such as in Wang et al.,45 only work when the ratio CP:CT is strictly 1. It is critical to eliminate the ratio restrictions and develop a framework suitable for more general cases. Indeed, the ratio between heterotypical cell concentrations changes significantly in the near-wall region of a parallel-plate flow chamber, and it cannot keep the 1:1 ratio as in the inlet set by the experiment. The near-wall region defined here is about 8 µm away from the substrate, equivalent to the height of an undeformed PMN tethered on the substrate.

Our goal is to propose, with the help of experimental studies and CFD simulations, a multiscale model, which can incorporate both cell population dynamics on the macro scale and bond formation kinetics on the micro scale, to characterize the PMN–melanoma cell aggregations and to address all the issues mentioned above. In correspondence with the two-step aggregation process, the main components of our modeling efforts include the development of collision rate and adhesion efficiency models which are used in population balance models, and the prediction of the aggregation percentage and relative variables.

To be more specific, we first set up a model to describe the dynamic coagulation processes between PMNs and melanoma cells. This model is based on the population balance equations (PBEs), originally proposed by Smoluchowski,41 which were also named as coagulation equations,2 but in a more specialized form. There are many investigations on applying PBEs to model the aggregation problem in a linear flow, e.g., generated by a cone-plate viscometer.21 However, the collision kernels, which describe aggregation properties including the collision rate (the number of collisions per unit time) and the adhesion efficiency (the probability of aggregation after collision) are mostly not applicable to nonlinear flow (e.g., in a parabolic velocity-profile shear flow developed in a parallel-plate flow chamber). Moreover, the effects of cell deformation on cell–cell aggregation are not taken into account. Hence, we develop a new collision rate formulation to incorporate both cell deformation and general flow conditions. And then, for the adhesion efficiency determination, we used the master equations,6,28,30,51 which describe the time evolution of the probabilities associated with the receptor–ligand bond formation kinetics (see Eq. 11 in later discussion), to avoid the complicated force calculations in Monte Carlo simulation.11,19,20,34,43,45

Our newly developed modeling framework provides an effective tool to study the aggregation between different kinds of deformable cells under general flow conditions at the population level, while eliminating the restriction on the ratio between the heterotypic cell concentrations and accounting for the effect of cell deformation to the aggregation process. We define the aggregation percentage as the percentage of PMNs with an attached tumor cell out of the total PMN population and show that it is both shear rate and shear stress dependent. The understanding of how the aggregation percentage changes with respect to the system parameters can be used to provide theoretical guidance to clinical studies of tumor metastasis, for example, it offers a new insight to the possible control of the aggregation process through the change of PMN population. Our numerical simulation results agree well with literature.26 We also demonstrate the existence of the seemingly contradicting threshold phenomena for receptor–ligand bond formations observed in many in vitro assays and provide an explanation for such observations. The flow affects aggregation by altering important parameters including, but not restricted to, the intrinsic binding molecules property, the extent of tethered cell deformation, the velocity profile of flow driving cells, and the heterotypic cell concentrations near the substrate (and therefore the ratio between the cell concentrations). To assess the relative importance of different parameters clearly, sensitivity analyses are conducted, and we come to a conclusion that the reaction coefficient and the critical bond number for adhesion efficiency play the most important roles in controlling the aggregation process.


Parallel-Plate Flow Experiment

The view of the cells in a parallel-plate flow chamber is shown in Fig. 1. Briefly, a syringe pump (Harvard Apparatus, South Natick, MA) was used to generate a steady parabolic laminar flow field in the flow chamber, where a confluent EC monolayer (as a ligand-binding substrate) was present. The flow channel is 800 µm long (direction of the flow) by 600 µm, with a height around 127 µm. The image focal plane was set on the substrate. The flow chamber was perfused with predetermined PMN and WM9 melanoma cell populations (1 × 106 cells mL−1) at 1:1 ratio. PMNs were pre-stimulated with 1 µM fMLP for 1 min or 1 ng mL−1 IL-8 for 1 h before perfusion into the parallel-plate flow chamber. After allowing PMNs and WM9 cells to reach the near-wall region (under a very slow flow rate for 2 min), shear stresses were adjusted to the experimental range of 0.625–2 dyn cm−2 and kept constant for 6–7 min.

Aggregation between a tethered PMN and a tumor cell in the parallel-plate flow chamber at shear rate 200 s−1 and viscosity 1.0 cP (Flow is from right to left). (a) At 0 s, a tumor cell and a PMN reaching the substrate; (b) At 10 s, collision between ...

Population Balance Model

Our discrete population balance model was built upon the widely accepted continuous population balance equations model (PBEs)2,41 with two additional assumptions. First, no aggregation event occurred in the free stream near the wall. Namely, there were only two types of cells in the near-wall region, PMN monomers and WM9 monomers. Hence, the aggregation only took place between the tumor cell monomers in the free stream near wall and the tethered PMNs on the substrate. Second, the flow was in steady state, namely, the concentrations of PMNs and tumor cells in the near-wall region are constants independent of time.

Let N(i,j,t) denoted the population of cells tethered on the substrate composed of i tumor cells and j PMNs adhered to the substrate at time t. N(i,j,t) would increase

  • by aggregation of the tethered particle composed of i − 1 tumor cells, j PMNs, and one tumor monomer in the near-wall region;
  • by aggregation of the tethered particle composed of i tumor cells, j − 1 PMNs, and one PMN monomer in the near-wall region;

and would decrease

  • by aggregation with a PMN monomer in the near-wall region;
  • by aggregation with a tumor cell monomer in the near-wall region.

Let β(i,j;i′,j′) be the aggregation kernel between one particle composed of i tumor cells, j PMNs and another particle composed of i′ tumor cells, j′ PMNs, which described the number of aggregations between PMNs and tumor cells per unit time per unit concentration. Thus, the population balance model for simulating the aggregation in the near-wall region was given by



where CT and CP are the concentrations of tumor cells and PMNs in the near-wall region, respectively.

PMN–tumor doublets and PMN monomers were the only observed cells that tethered on the substrate, while larger particles composed of more than one PMN and one tumor cell, such as triplets, were very rare.26 Hence, only four types of cells were taken into account in our model, including: free PMN monomers, free tumor cell monomers in the near-wall region, PMN monomers adhered on the substrate and tethered PMN–tumor cell doublets. To simplify notation, in Eqs. (1) and (2), the terms N(1,1,t) and N(0,1,t) were denoted by NPT and NP, respectively, and β(0,1;1,0) was denoted as βPT. The new equations then became:



with an analytic solution given by:



where Tf was PMN tethering frequency that described the rate of PMNs attached to the substrate; NP0andNPT0 were the initial populations of tethered PMNs and PMN–tumor doublets on the substrate, respectively.

Among the parameters in Eqs. (5) and (6), the coagulation kernel βPT, generally defined as the product of collision rate and adhesion efficiency, would be discussed in the following two sections. NP0,NPT0 and Tf were experimentally recorded data in Liang et al.26 The near-wall concentrations CP and CT were evaluated analytically based on the results of Munn et al.,32 which showed that the concentration of cells in the near-wall region of a parallel-plate chamber obeyed the following rule:


where Hdep was the depth of field of the microscope objective used in the experimental observations (8.58 µm); Ψ was the concentration of cells within the region of height Hdep; C was the bulk concentration; vavg was the average settling velocity for cells above Hdep and uavg was the average cell convection velocity in the flowing direction (x direction) of height Hdep. The determination of settling and convection velocity would be discussed in detail in the next section Eqs. (9 and 10). The total length of our parallel-plate flow chamber was 2 cm, and we evaluated the concentration at the midpoint (x = 1 cm) to approximate the average near-wall cell concentration. Using the corresponding bulk concentrations of PMNs and tumor cells, one can find CP and CT easily from Eq. (7).

Development of Collision Rate Model

The collision rate, [beta]PT, was defined as the number of collisions between a tethered PMN and flowing tumor cell monomers per unit time per unit concentration. We evaluated it numerically by a collision model. If particles were driven by the flow, collisions took place due to their different velocity distributions. Following Eq. (7) and the steady state assumption, we presumed that each kind of cell was uniformly distributed in the region of the same height, while the distribution varied with respect to height. Any cell of radius r2 would collide with a target cell of radius r1 if its center reached the boundary of the target’s concentric sphere of radius r1 + r2, a so called collision region (Fig. 2).

The collision between two cells with radius r1 and r2, respectively.

Notice that similar geometric descriptions of the collision region may be used to compute the collision rate for more general particle shapes. For instance, with a given velocity field vrel, if a target cell of any shape is going to collide with another cell during the time interval [t, t + Δt], then the collision rate corresponding to the collision region Ω is the flux (number of cells) that reached the boundary of the collision region from out-side. Hence, the collision kernel was given by


where n is the outward unit normal vector of the collision region.

Therefore, the collision rate between deformable cells can be obtained by integral, given the shapes of the deformed cells and the velocity profile. This formulation was not restricted to a linear shear flow field or rigid spherical cells. This approach, to the best of our knowledge, makes the present work the first attempt to reconcile the population balance model parameter determination that accounts for the cell deformation.

To apply the collision model to our parallel-plate flow chamber assay, the velocity profile was determined first. The trajectory of a cell would depend on gravity, which causes the cell to settle, and flow convective force, which drives the cell moving horizontally. Taking rt to be the radius of a cell; μ to be the fluid viscosity; ρt and ρm to be the cell density and the media density, respectively; h to be the distance from the substrate to the center of the cell, which represented the position of the cell, Davis and Giddings7 proposed that the settling velocity of the cell can be approximated by


If a cell stayed very close to the chamber substrate, its convective velocity was reduced by hydrodynamic wall effects. For h < 4rt, the convective velocity was approximated by the far-field asymptotic formula10,48,49


with near-wall shear rate [gamma with dot above]. We used the velocity profile given by the settling velocity vs in Eq. (9) and the convective velocity vc in Eq. (10) to calculate the collision rate in the parallel-plate flow chamber.

Concerning cell deformability, we assumed that every tumor cell was a rigid spherical particle which did not deform in the flow, while the shape of a tethered PMN was subject to change under various flow conditions. The shape of a tethered PMN was taken as an arc with a fixed volume, as shown in Fig. 3. Since H was significantly smaller than the height of the chamber, we assumed that the deformation of PMNs would not affect the flow conditions.

The collision of a tumor monomer and a deformable tethered-PMN with height H.

The height H was determined by the shear rate [gamma with dot above] and medium viscosity μ according to H = f([gamma with dot above], μ), where f satisfied the boundary condition f(0,μ) = 2rp. Namely, a PMN was a sphere with radius rp without the influence of external forces. In general, the shape of a deformed PMN can be observed by experiments for any given shear rate and viscosity. Combining with the velocity profile around a tethered PMN, we could compute the collision rate βPT by estimating the integral in Eq. (8). The parallel-plate flow chamber assay only counted the number of collisions between PMNs and tumor cells (unpublished data in Liang et al.,26 Table 1) instead of reporting their shapes. Hence, the determination of the function f was accomplished through an inverse procedure. For example, we varied H from 0 to 2rp with a small step size, then computed the corresponding collision number and compared with the experimentally recorded total collision numbers during 5 min. The H providing the best match was taken as the corresponding height under that flow condition. This procedure helped us to develop the height dependence on the flow.

Total collision number.

Adhesion Efficiency Model

The adhesion efficiency, εPT, was defined as the probability that two colliding cells would stick to each other by forming receptor–ligand bonds to overcome the hydrodynamic force acting on the doublet to separate them. We evaluated εPT via the master equation presented below, which described the binding kinetics of a small number of receptors and ligands following the physiological mechanism.6,28,29,35,51 Taking Pn(t) as the probability of existing n bonds at time t, the rate of change of Pn(t) with respect to time was given by


where Nr and Nl were the densities of receptors on tumor cells and ligands on PMNs, respectively; Ac was the contact area between a tethered PMN and an attached tumor cell; Kon(n) was the forward reaction coefficient per unit density with n bonds connecting two cells and Koff(n) was the reverse reaction rate per cell with the existence of n bonds. The first term on the right hand side of Eq. (11) described the generation of n bonds by forming one extra bond from the existed n − 1 bonds; the second term stood for the loss by either the formation or break-up of a bond from the existed n bonds; the third term described the generation through breaking one bond from the existed n + 1 bonds. Additional assumptions were made for further simplification. For example, Ac, Nr and Nl were assumed to be constants under a fixed flow environment; Kon(n)andKoff(n) were assumed to be independent of the existed number of bonds n; and the actual number of formed bonds was significantly smaller than the maximum possible number of bonds N = Ac min(Nr, Nl) that could be formed by the available receptors and ligands.

Bond formation simulations in CFD and the parallel-plate flow chamber assay were analyzed on two different time scales. The bond formation molecular-level dynamics was simulated on the order of microseconds (Δt),14,24,25 while the lab experiments was done within 6 min collecting only the cellular-level data.26 The development of cellular-level bond formation kinetics was thus understood as an equilibrium adhesion model. The adhesion distribution reached the steady state probability distribution satisfying the following system:


with N standing for the size of the system and A = AcNrNlKon/Koff. The equilibrium state distribution solution depended on the initial condition of the system, the contact area Ac, the concentrations of receptors Nr and ligands Nl, the forward and backward reaction rates Kon and Koff, and the size of the system N. A careful analysis of Eq. (11) revealed that the existence of a unique equilibrium state as a probability distribution was independent of the initial condition. Moreover, the steady state was not sensitive to the size of the system when N was large enough, and Pn converged to eA An/n! (Appendix). Therefore, the constant A was the only essential parameter for the determination of the equilibrium state solution.

CFD simulations24,25 provided us with the peak bond force per bond and hydrodynamic force under different flow conditions (Table 2). The critical bond number, denoted as n*([gamma with dot above], μ), was the least number of bonds required for the firm adhesion, which was given by the ratio of the hydrodynamic force over the peak bond force for the shear rate [gamma with dot above] and viscosity μ (Table 2). Then, the probability of aggregation corresponded to the probability of having more than n*([gamma with dot above], μ) bonds at the equilibrium state. Notice that n* was not necessarily a whole number.

Hydrodynamic force and peak bond force under different shear conditions.

We proposed two formulations to evaluate the adhesion efficiency, both of which depended only on A and n*:




Here, [n with macron] was the smallest integer bigger than n, and stood for the largest integer smaller than n. The first formulation Eq. (13) was intuitive: if we summed up all the possibilities that the number of bonds exceeded n*, the bond force would be sufficient to overcome the hydrodynamic force and hence form the firm adhesion. The formulation Eq. (14) accounted for the fractional terms in critical bond number so as to take full consideration of all probabilities. The choice made between these two formulations would be explained in the “Discussion” section.

We proposed to find the reaction parameter A via an inverse method, which shared exactly the same idea as the prediction of the height of a deformed PMN. In Eqs. (5) and (6), except for adhesion efficiency, all the other parameters, including the final attached number of cells, were available. A reaction constant A providing the correct adhesion efficiency for the proper aggregation number would then correspond to the desired value of A.

Tethering Frequency

The tethering frequency, Tf, was a standard concept which reflected the rate at which PMNs adhered to EC monolayer under different flow conditions. Experimental records26 have demonstrated its dependence on shear rate. Other than that, the cell deformability would also affect the cell tethering. The capillary number, Ca, is a non-dimensional parameter which characterizes the relative importance of viscous forces and the cell deformability in the flow system. In our case, the deformability of the cell was characterized by the membrane tension of the PMN. The Capillary number equals the Weber number, We, divided by the Reynolds number Re.


where ρm was the fluid density (1000 kg m−3); d was the PMN cell diameter (~8 µm); σ was the membrane tension (3.1 × 10−5 N m−1); vavg was the average velocity of the cell, approximated by 12γ˙d and [gamma with dot above] was the shear rate near the substrate. A correlation was assumed between tethering efficiency and the capillary number, Ca. The least square fitting method was used in order to fit the tethering frequency as a function of the wall shear rate [gamma with dot above] and the capillary number Ca. The parameters and data used to fit the tethering frequency Tf are listed in Table 3.

Experimental recorded tethering frequency data.

Aggregation Percentage

Aggregation percentage, α, was defined as the number of PMN–tumor cell doublets divided by the total number of tethered PMNs, including those PMNs that have one or more tumor cells attached. The expression of aggregation percentage was the following:

α=the aggregation number of PMNtumor cell doubletsthe total adhered PMNs at the substrate×100%

This ratio reflected the ability of PMNs to facilitate the aggregation of tumor cells to the substrate and the extent of tumor cell aggregation to the substrate. Recorded experimental data would be compared with the numerical predictions of our coagulation model given by Eqs. (5) and (6). For some very rare cases, e.g., more than one tumor cell attached to one PMN, they were counted as two or more aggregations instead of one. The aggregation percentage of our model was then expressed as follows:


The second equality remained true under the assumption that initially there were no tumor cells attached to PMNs.

Determination of Model Parameters and Discussion

A flow chart for conducting the numerical simulation of the population based aggregation model is given in Fig. 4. Following the flow chart, the aggregation percentage under any flow condition could be computed with the given parameters, including the number of initially attached PMN–tumor doublets and PMN monomers (NPT0andNP0), the height of the deformed PMNs (H), the critical bond number (n*) and reaction constant (A) for adhesion efficiency, and the tethering frequency (Tf). These parameters, except for the number of initially tethered PMN–tumor cell doublets (NPT0), which was assumed to be zero, were collected either from literature (NP0,n*andTf)14,25,26 or from the solutions of the inverse problems (H and A) as explained for the determination of the height of deformed PMNs tethered on the substrate. When the viscosity was 0.625 cP, we used curve fitting to predict the corresponding parameters. The optimum values were determined using least square methods and then we made predictions on function values for other inputs. The values of R2 were reported and the residuals were tested for normal distribution hypothesis using the Shapiro–Wilk goodness of fit test with a significance level of 0.05, so as to demonstrate the validity of our prediction.

Flow chart of population balance model simulation.

The Effect of Variation in the Ratio of PMNs vs. Tumor Cells

It should be noticed that we did not have restrictions on CP:CT in the construction of our model. However, CP was not a parameter that appeared explicitly in our formulation. Instead, it controlled the tethering frequency (Tf) and the number of initially tethered PMNs (NP0), both of which were key parameters in our model. Without further knowledge concerning how these three variables were connected, we presumed that the change of CP was linearly proportional to the change of the tethering frequency and the number of initially tethered PMNs. By increasing the concentration of PMNs in the inlet, these two parameters increased to Itf times and Iini times of the increment of CP, respectively. Following Eq. (15), it was easy to analytically prove that the aggregation percentage monotonically increased with Iini/Itf, and when Itf = Iini the increase in CP would not affect the aggregation percentage at all. In other words, if the increase in NP0 was more dramatic in comparison with the increase in tethering frequency, it would promote aggregation and vice versa.


Deformed PMN Cell Height Prediction

Table 4 lists the relation between the height of deformed PMN cells and flow conditions solved by inverse problem. We hypothesized that the height function f([gamma with dot above], μ) expressed in the following way,


where a, b and c are fitted through the data in Table 4, which gives a = 6.574 × 10−6, b = 4.663 × 10−6 and c = 3.133 × 10−6. Notice that the height is on the order of 10−6, so the coefficients are picked with the same order. The R2 value is reported in Table 5. The difference between the prediction function and the experimental data was tested to satisfy a normal distribution with a significance level of 0.05. Thus we have the height prediction function as follows:

The height of a deformed PMN vs. varying shear rates and viscosities.
Lists of data that are estimated.

Tethering Frequency Regression

Nonlinear regression was conducted with the help of least square method to predict the relation among the tethering frequency (Tf), shear rate ([gamma with dot above]) and the Capillary number (Ca) based on Table 3. Using standard fluid dynamics dimensional analysis, a correlation was found between the Capillary number and the tethering frequency, which verified our assumption. Figure 5 shows how Tf regresses on γ˙1.6429Ca0.8046, which gives a good fit to the data. The R2 value of the approximation is 0.9459, and we accept the normality hypothesis of the prediction residuals at a significant value of 0.05 based on the Shapiro–Wilk normality test, which means it is a reasonably good prediction.

Fit of tethering frequency Tf as function of γ˙1.6429Ca0.8046 on a log–log scale.

The relation given by Eq. (17) is further explored in order to understand how shear rate and shear stress would affect the PMN tethering process. We first change shear stresses with varying shear rates ranging from 50 to 500 s−1, as shown in Fig. 6a. Then we manipulate the shear stresses variation under a constant shear rate by varying viscosity in the range from 1.0 to 4.0 cP. Results are plotted in Fig. 6b. Both of these figures demonstrate that shear rate and shear stress have positive influences on the tethering frequency of PMNs.

Effect of shear rate and viscosity respectively on tethering frequency Tf. (a) Solid line represents viscosity of 1 cP, dash–dot line represents viscosity of 2 cP; dash line represents viscosity of 3.2 cP; (b) Solid line represents shear rate ...

Time Evolution of Aggregation Percentage

Figures 7a–7c show the aggregation percentage vs. time for shear rates of 62.5, 100, and 200 s−1, respectively. These figures are illustrations of the time evolution of aggregation percentages of the experimental data in Liang et al.,26 where only the final aggregation percentages at the end of the 300-s time interval are given. Each part is plotted with a fixed shear rate but varying viscosity (i.e., changing shear stresses at a constant shear rate). Figure 7d compares the dependence of aggregation percentage at the end of 300 s on shear stress under various flow conditions mentioned in Figs. 7a–7c. From literature,50 if shear stress is large enough, the receptor–ligand bonds are more likely to break, and hence depress the likelihood of the aggregation. It is then natural for us to expect a monotonic decrease in aggregation percentage. A careful observation of Fig. 7d shows that with a shear stress of 100 dyn cm−2, the percentage increases first before decreasing. Under a fixed 200 dyn cm−2 (varying shear rate), the three observations display relatively large variations (standard derivation is more than 20% of the mean of observations). Hence, we hypothesize that the influence on aggregation comes from the influence of shear rate and shear stress, rather than any one of them alone. With a fixed viscosity, the shear rate makes a negative impact on the aggregation.

The time evolution of aggregation percentage responses with respect to various flow conditions: Varying shear stress by varying viscosity at fixed shear rate, including (a) Shear rate = 62.5 s−1; (b) Shear rate = 100 s−1, (c) Shear rate ...

With viscosity fixed at 0.625 cP, we explored necessary parameters (including H, A, Tf, n*) before simulation using the same method for height function prediction. The R2 values for the fitted parameters are listed in Table 5. The predicted aggregation percentages based on estimated parameters for shear rates 62.5, 100, and 200 s−1, respectively, are comparable with the experimental findings26 as showed in Fig. 8. Results from the numerical simulation suggest that the aggregation between PMNs and tumor cells displayed an increasing trend for the first couple of minutes before reaching a steady-state. At higher shear rates, the time required to reach the plateau would be shorter comparing to lower shear rates. Both of these simulation findings agree with experimental records in literature.26

The time course of aggregation percentage at fixed viscosity 0.625 cP.

In Fig. 9, we compare the aggregation percentage at 300 s under small shear rates at 1 cP. Contrary to the intuition and the observation in Fig. 7d, as the shear rate (stress) increases, we observe that the aggregation percentage trend is also increasing during this period. Such a behavior was also reported in the literature,12,13 and we will provide further explanation in the “Discussion” section.

Aggregation percentage at 300 s under viscosity 1 cP.

Parameter Sensitivity Analysis

To discuss the sensitivity of the aggregation percentage curve quantitatively, the evaluation of the influence of different parameters on the short (60 s), medium (180 s) and long term (300 s) aggregations are compared. Different sets of aggregation percentages are recorded with a change of one important parameter at a time, including the initial condition of simulation, the tethering frequency, the near-wall cell concentration, the height of a deformed PMN, the critical bond number for firm adhesion and the reaction constant for receptor–ligand bond formation. Here, we restrict the variation of parameters within 10%, and we compare the increasing or decreasing ratio of the parameters with respect to that of the aggregation percentage to understand the sensitivity between different parameters. No significant differences are observed among short, medium and long term aggregations. The long-term change ratio of aggregation percentage vs. the change ratio of initial condition of simulation, tethering frequency, near-wall cell concentration and the height of deformed PMNs, respectively, are listed in Fig. 10, while the rest are showed in Fig. 11. The initial condition and tethering frequency are found to have the least effect on aggregation percentage. A greater number of initially tethered cells would contribute positively to the increase in the percentage of PMNs that have tumor cells attached to them, but the tethering frequency would impose a negative impact on aggregation percentage. The critical bond number also has a substantially negative correlation to the PMN–tumor cell aggregation percentage, while the concentration, height of deformed cells and reaction constant are only mildly correlated to aggregation.

Sensitivity of aggregation percentage with respect to the initial condition, tethering frequency, tumor cell concentration and height, when shear rate = 100 s−1, viscosity = 0.625 cP.
Sensitivity of aggregation percentage with respect to the reaction constant and critical bond number for adhesion efficiency when shear rate = 100 s−1, viscosity = 0.625 cP.

The Effect of Variation in the PMNs:Tumor Cell Ratio

We increase the concentration of PMNs in the inlets by 4 times (2 × 106 cells mL−1) to get CP:CT = 4 at the inlet and plot the time course evolution of aggregation percentages. Without information on how parameters change with respect to the variation of CT, our investigation on the effect of the ratio change in the model is achieved by considering a number of different cases: first in case I, both the initially attached PMNs (NP0) and tethering frequency (Tf) of PMNs are simply increased to 4 times of the case that the concentration of PMNs is 5 × 105 cells mL−1 (Fig. 12a); then in cases II, only the NP0 is multiplied by 4 with Tf being fixed (Fig. 12b); while in case III, only Tf is multiplied by 4 and NP0 being fixed (Fig. 12c).

Comparison of the influences of PMN:tumor cell ratio on aggregation percentage at viscosity 0.625 cP at various shear rate. (a) Fixed 4:1 for PMN:tumor cell concentration ratio at the inlet, by increasing both the initial tethered PMN number and tethering ...

Compared to the case of a 1:1 ratio, Case I reveals almost no influences on the aggregation (identical to Fig. 8), while case II leads to a promotion in aggregation percentage and case III demonstrates a reduction in the aggregation.


Choice Between Two Adhesion Efficiency Models

A choice was made between the two adhesion efficiency models, Eqs. (13) and (14), based on the numerical simulation results. The results reported here were numerically implemented according to Eq. (14). The rationale behind such a choice was mainly due to the issue that the number of bonds (as the subscript of Pn(t), the probability density) could only be whole numbers. The second formulation (14) was adopted to control the oscillations that would otherwise appear in the predicted results due to the jumps associated with the step functions. Indeed, because the critical bond values n* changed abruptly due to the integer rounding operation, even though the Eq. (13) inherited the negative correlation to shear rate, its predicted adhesion efficiency exhibited a variance of 200% with respect to average value. To make further clarification, we consider, for example, the case where a very small change occurred to the shear rate. With the viscosity fixed, it was reasonable to expect a small change in the adhesion efficiency. However, if n* happened to change from 3.98 to 4.02, the adhesion efficiency provided (εPT) by Eq. (13) would have a huge jump as the εPT changed from the summation of P4 to PN into the summation of P5 to PN, while (14) would give a 98% probability of aggregation to P4 when n* is 4.02, which helped to avoid the unrealistic oscillation. Thus the second formulation is more reasonable and was taken to yield more gradual changes in the estimation.

Explanation of the Inconsistency at High Shear Rate Low Viscosity

Based on Figs. 7a–7d, we see that the aggregation percentage is modulated by both shear rate and shear stress. As plotted in the Fig. 8, the numerical simulation results at 0.625 cP, 200 s−1 are obviously smaller than the experimental record especially from the third to fifth minutes. Here we outline some possible reasons from both the tethering frequency and the cell concentration point of view to explain the apparent inconsistency.

Even though the dependence of adhesion on tethering frequency does not play the most essential role as explored in sensitivity analysis, given a dramatic change, we would still expect to see its effect on aggregation. The tethering frequency data we used here are experimentally determined. It is the number of both tethered and steady rolling PMNs. But our model treat the case that only the firmly adhered PMNs would contribute to the formation of doublets; a larger tethering frequency was thus introduced in the simulation. Intuitively, the percentage of rolling PMNs increases with shear rate, which could contribute to the underestimation shown in Fig. 8.

The parallel-plate flow chamber experiment was conducted with a notable decrease of cell concentrations throughout the 5-min experiment, because the injected cell concentration decreased from the initial time due to sedimentation before being pumped into the flow chamber. A fixed weighted correction was used for the tumor cell concentrations to reconcile this problem for all flow conditions, but not on the PMNs. The tethering frequency would be affected by the number of PMNs in the close-wall region, which served as the source of tethered PMNs. Averaged over 5 min, the tethering frequency is larger than the recorded constant initially but smaller toward the end of the experiment due to the change of near-wall cell concentration. Thus, the simulation took a larger tethering frequency at the latter half of the simulation procedure so that we would expect a smaller decrease with higher shear rate, and a more pronounced inconsistency also.

Aggregation Threshold Phenomena

Several labs, for example, Hentzen et al.12 in their cone-plate experiments and Hinds et al.13 in their vertical parallel-plate flow chamber experiments, reported the existence of a threshold below which the increasing of shear stress would promote the aggregation between PMNs and tumor cells, seemingly contrary to conventional wisdom of the effect of shear stress upon adhesion for receptor–ligand bonds. This kind of aggregation threshold phenomena has been a focus of study in several labs, e.g., Zhu52 at GaTech. Figure 9 displays the same outcome and proves the existence of such threshold for receptor–ligand bonds. We propose that the existence of the threshold comes from the fact that when shear rate is below the threshold, the number of bonds required to hold two cells together is mostly just 1 or 2. Given sufficient contact time (inversely proportional to the small shear rate), the number of formed bonds is always sufficient to tie two cells together. When shear rate keeps increasing, the critical bond number increases dramatically. With less contact time, it becomes harder to form aggregates and the cells tend to stay separated, which would result in a smaller aggregation percentage. Notice that the threshold phenomena have been commonly proposed as the influence of the shear stress, but we conjecture that the influence of the flow condition might be more correlated with shear rate. Follow up research will help verifying our conjecture.

Influence from PMNs and Tumor Cells Ratio on PMN–Tumor Cells Aggregation

The size difference between PMNs and tumor cells would result in a different cell concentration at the near-wall region.32 Intuitively, the smaller PMNs are more easily transported by flow, while the bigger tumor cells would experience more friction drag and settlements, and hence, tend to settle to the substrate. The ratio of concentration in the near-wall region of parallel-plate flow chamber would be dramatically changed comparing to the inlets. Experimentally, the change of near-wall concentration with respect to different ratios could not be manipulated easily, but we could accomplish it through numerical simulations. Due to the lack of further studies on how cell concentrations in the inlet would affect the tethering frequency and the number of initially tethered PMNs, we only provided a framework which is capable of discussing the importance of the ratio parameter. Figure 12 demonstrates our analytically predicted results that, given the linear proportionality assumption, a larger change in the number of initially tethered PMNs would promote aggregation, while a larger change in tethering frequency restrained the aggregation percentage from increasing. We conjecture that an inflammation environment (with a larger number of NP0) provided more anchors initially for tumor cells to attach, which would promote the aggregation process. However, if a dramatic increase in tethering frequency would follow later, the anchors would then be overstock; so there would be a smaller percentage of tethered PMNs having a tumor cell tethered on it. More concrete conclusions could be drawn with further studies on the relation between CP, Tf and NP0.

Identification and Sensitivity Analysis of Different Parameters

Notice that the parameter identification and prediction has played an important role in our simulation. The determination of parameters we mentioned above is certainly model based; simulation results would be affected by the change of parameters due to the change of models. Hence, identifying the most influential parameters in the aggregation process is more important than determining the numerical values of the parameters. Mathematical analysis of the model can help to illustrate the qualitative importance of different parameters. The working formula for aggregation percentage Eq. (15) was an increasing function of the initial condition NP0, the coagulation kernel βPT, the tumor cell concentration CT, and a decreasing function of the tethering frequency Tf, which agreed with our numerical simulations.

Limitations and Future Model Refinements

Our proposed population balance model as well as the collision model provides results comparable to biological in vitro experiments; it reveals the important factors in the aggregation process, and goes beyond experiments by providing new predictions that have not been achieved by in vitro experiments. However, further improvements can be made in several ways, both experimentally and theoretically. The presence of red blood cells (RBCs) has been proved to affect the transport of white blood cells (WBCs).31,33 To be more specific, the RBCs which constitute approximately 40% of the volume of blood would drive more tumor cells and PMNs to the near-wall region, which increases the concentrations and adjusts the convection trajectories of PMNs and tumor cells. Therefore, RBCs would also play an important role in stimulating coagulation process. Putting red blood ghost cells into the experiment would allow us to simulate more realistically the physiological conditions and to improve the experimental procedure that provides better identification of all parameters. This is the subject of our ongoing research.

Meanwhile, several modeling assumptions could also be relaxed. For instance, in collision rate modeling, although we have incorporated cell deformation as a factor based on the previous numerical studies of cell deformation in shear flow,17,18,22 the shape of deformed white blood cells is more likely to be a tear drop8 instead of a spherical cap as in our current simplified assumption. One of the possible improvements is to solve the fluid–structure coupling problem to evaluate quantities, such as fluid velocity profiles and PMNs’ deformation. A possible procedure could be: with boundary condition for a given shape of the PMN, we evaluate the fluid convection to get hydrodynamic force combined with bond kinetics for the bond force, and then simulate the PMN deformation. The newly deformed PMN provides an updated boundary for our fluid model, and therefore, new force balance results. This coupling process evolves and possibly stops when an equilibrium state is reached. Then, based on Eq. (8), we have a complete collision determination mechanism, which leads to more challenging numerical implementation. In addition, all the rolling cells on the substrate are not described in our population balance model, although the steady rolling PMN might also contribute to the PMN–tumor cell aggregation. Some aggregations may also happen in the free stream, before any PMN or tumor cell lands on the substrate. A thorough consideration of those phenomena would definitely lead to a more sophisticated model. The spatial distribution of cells is also expected to affect the cell aggregations, which is yet to be discussed. Sabelfeld and Kolodko36 proposed the spatially inhomogeneous Smoluchowski equation which helped to remove the constant cell concentration assumption. Similar works in this direction can be pursued in the future.


In this paper, a comprehensive modeling of cell coagulation was accomplished in cooperation with experiments as well as CFD simulations. To our best knowledge, we incorporated for the first time the cell deformation factor under given shear flow conditions into the population balance modeling of the aggregation between tumor cells and EC facilitated by tethered PMNs. Our model was able to produce aggregation percentages that match well with the experimental data. More importantly, it removed the restriction on cell concentration ratios. Our finding illustrates that biologically the effectiveness of increasing PMN population on lowering the aggregation percentage can be improved when the PMN tethering frequency increases proportionally even if the number of initially tethered PMN remains unchanged; or when the PMN tethering frequency undergoes a little change and the number of initially tethered PMN decreases significantly. Our model was also able to predict the aggregation threshold phenomena. A step by step discussion on parameter identification was conducted. We theoretically analyzed the influence of different factors, e.g., tethering frequency, collision rate, and adhesion efficiency, on the cell aggregation. Complementing the limited experimental studies available in the literature, this analysis provided an effective approach to gather more quantitative information on the relative importance of all relevant factors in the complex aggregation process.


The authors thank Dr. Meghan Hoskins for providing simulation data on hydrodynamic force. This work was supported by the National Institutes of Health grant CA-125707 and National Science Foundation grant CBET-0729091.


The reaction constant as a combination of Ac, Nr, Nl, Kon and Koff (s−1)
The contact area between a PMN and a tumor cell (µm2)
Ca, We, Re
The Capillary number, the Weber number and the Reynolds number
The collision number between two kinds of cells (events)
The concentration of PMNs, tumor cells and bulk cell concentration, respectively, in the near-wall region (cells mL−1)
The diameter of an undeformed PMN cell (µm)
Acceleration due to gravity (m s−2)
H, Hdep, h
The height of a deformed PMN cell attached to the substrate, the height defined by the depth of field of the microscope objective used in the experimental observations and the distance from the substrate to the center of the cell, which represented the position of the cell (µm)
Iini, Itf
The increasing multiple of initial condition or tethering frequency
The incoming flux to collision region near PMN per concentration (kg−1 m3)
Kon, Kon(n)
The forward reaction rate per unit density for bond formation (s−1 µm−2)
Koff, Koff(n)
The backward reaction rate per cell for bond breakage (s−1)
Outward unit normal (vector) of collision region
The approximation of the smallest number of bonds required for firm adhesion (bonds)
Nr, Nl
The concentration of receptor and ligand on cells (µm−2)
The number of tethered PMN–tumor cell doublets and tethered PMN monomers on the substrate (cells)
Pn(t), Pn
The probability of having n bonds.
rp, rt
The radius of an undeformed PMN cell and a tumor cell (µm)
The tethering frequency of cells: including firmly adhered cell and rolling cells (events s−1 per view)
vavg, vrel
The average settling velocity for cells above height Hdep and the relative velocity (difference between velocities) of two kinds of cells (µm s−1)
vs0, vs, vc
The free settling velocity, settling velocity and convection velocity of cells in the parabolic flow profile (µm s−1)
The aggregation percentage
βPT, β(i,j;i′,j′)
The coagulation kernel of tumor cells in the near-wall region to the tethered PMNs (µm3 s−1)
[beta]PT, [beta]
The collision rate between tethered PMNs and tumor cell monomers (µm3 s−1)
[gamma with dot above]
The shear rate in the close-wall region (s−1)
The adhesion efficiency between tethered PMNs and tumor cell monomers
The viscosity of the fluid flow (Poise)
ρm, ρt
The fluid density, the tumor cell density (kg m−3)
Membrane tension (N m−1)
The concentration of cells within the region of height Hdep
The collision region


Computing the Collision Rate [beta]PT

We use the spherical coordinates to express all the related parameters so as to evaluate the integral in Eq. (8) with a given velocity profile. In our rectangular coordinates, we suppose the chamber cross section (Fig. 3) is in yz-plane and the x-axis points to the reader, so that

v=(vx,vy,vz) and vx=0,vy=vc,vz=vs.

We take the standard spherical coordinate system with the origin being the center of the sphere where the arc shaped PMN lies on and the two bottom points of PMN having spherical coordinates (r,φ0,π2)and(r,φ0,3π2). Assume that the deformable PMN preserves its volume V=43πrp3, where rp is the radius of an undeformed PMN, we thus have that


which gives a closed form for H, r and [var phi]0 via the following explicit formulae:

r=VπH3+H3 and φ0=arccos(1Hr).

Now consider a tumor cell which is colliding with this PMN. Assume that the contact point has coordinate (r + rt, Θ, [var phi]), where rt is the radius of a tumor cell. Thus, h, the distance between the center of tumor cell and the substrate, is


Our task is to compute the coagulation kernel under different flow conditions. For any given parameters [gamma with dot above] and μ, we compute H by the fitted function f in Eq. (16) first, and compute [var phi]0 and r next. We may then plug in all these parameters and the flow velocity profile into the Eq. (8) to compute the kernel by estimating the integral over the collision region Ω. That is, we estimate


where the integrand is given by


This integral is evaluated numerically via quadrature rules. The parameters that we need in the calculation are listed in Table 6.

Values for specific parameters.

Sensitivity to the Maximum Number of Bonds N

The algebraic system


has a general solution Pn=AnP0n!, for n = 1, 2, 3,…, with A and P0 satisfying


Under a given flow condition, the reaction constant A is fixed, by taking N → ∞,


Therefore, the solution Pn converges as N → ∞ to eAAnn!. By numerical simulation, we could actually find that when N ≥ 300, the solution and the limit are almost identical (Fig. 13). So we can truncate the original system to N = 300.

The steady state solution distribution of adhesion efficiency with the system size taken to be 300, 500, and 1000.


1. Abkarian M, Viallat A. Dynamics of vesicles in a wall-bounded shear flow. Biophys. J. 2005;89:1055–1066. [PubMed]
2. Aldous DJ. Deterministic and stochastic models for coalescence (aggregation and coagulation): a review of the Mean-Field Theory for Probabilities. Bernoulli. 1999;5:3–48.
3. Belval T, Hellums J. Analysis of shear-induced platelet aggregation with population balance mathematics. Biophys. J. 1986;50:479–487. [PubMed]
4. Caputo KE, Hammer DA. Effect of microvillus deformability on leukocyte adhesion explored using adhesive dynamics simulations. Biophys. J. 2005;89:187–200. [PubMed]
5. Chang KC, Tees DF, Hammer DA. The state diagram for cell adhesion under flow: leukocyte rolling and firm adhesion. Proc. Natl Acad. Sci. USA. 2000;97:11262–11267. [PubMed]
6. Chesla SE, Selvaraj P, Zhu C. Measuring two-dimensional receptor–ligand binding kinetics with micropipette. Biophys. J. 1998;75:1553–1572. [PubMed]
7. Davis JM, Giddings JC. Influence of wall-retarded transport of retention and plate eight in field-flow fractionation. Sci. Technnol. 1985;20:699–724.
8. Dong C, Cao J, Struble EJ, Lipowsky H. Mechanics of leukocyte deformation and adhesion to endothelium in shear flow. Ann. Biomed. Eng. 1999;27:298–312. [PubMed]
9. Dong C, Lei X. Biomechanics of cell rolling: shear flow, cell-surface adhesion, and cell deformability. J. Biomech. 2000;33:35–43. [PubMed]
10. Goldman AJ, Cox RG, Brenner H. Slow viscous motion of a sphere parallel to a plane wall. II. Couette flow. Chem. Eng. Sci. 1967;20:653–660.
11. Hammer DA, Apte SM. Simulation of cell rolling and adhesion on surface in shear flow: general results and analysis of selectin-mediated neutrophil adhesion. Biophs. J. 1992;63:35–57. [PubMed]
12. Hentzen ER, Neelamegham S, Kansas GS, Benanti JA, Smith LV, McIntire CW, Simon SI. Sequential binding of CD11a/CD18 and CD11b/CD18 defines neutrophil capture and stable adhesion to intercellular adhesion molecule-1. Blood. 2000;95:911–920. [PubMed]
13. Hinds MT, Park YJ, Jones SA, Giddens DP, Alevriadou BR. Local hemodynamics affect monocytic cell adhesion to a three-dimensional flow model coated with E-selectin. J. Biomech. 2001;34:95–103. [PubMed]
14. Hoskins M, Kunz RF, Bistline J, Dong C. Coupled flow–structure–biochemistry simulations of dynamic systems of blood cells using an adaptive surface tracking method. J. Fluids Struct. 2009;25:936–953. [PMC free article] [PubMed]
15. Huang P, Hellums J. Aggregation and disaggregation kinetics of human blood platelets: Part I. Development and validation of a population balance method. Biophys. J. 1993;65:334–343. [PubMed]
16. Huang P, Hellums J. Aggregation and disaggregation kinetics of human blood platelets: Part II. Shear induced platelet aggregation. Biophys. J. 1993;65:344–353. [PubMed]
17. Khismatullin D, Truskey G. A 3D numerical study of the effect of channel height on leukocyte deformation and adhesion in parallel-plate flow chambers. Microvasc. Res. 2004;68:188–202. [PubMed]
18. Khismatullin D, Truskey G. 3D numerical simulation of receptor-mediated leukocyte adhesion to surfaces: effects of cell deformability and viscoelasticity. Phys. Fluids. 2005;17:53–73.
19. Kolodko A, Sabelfeld K, Wagner W. A stochastic method for solving Smoluchowskis coagulation equation. Math. Comput. Simul. 1999;49:57–79.
20. Kwong D, Tees DFJ, Goldsmith HL. Kinetics and locus of failure of receptor ligand-mediated adhesion between latex spheres. II. Protein–protein bond. Biophys. J. 1996;71:1115–1122. [PubMed]
21. Laurenzi IJ, Diamond SL. Monte Carlo simulation of the heterotypic aggregation kinetics of platelets and neutrophils. Biophys. J. 1999;77:1733–1746. [PubMed]
22. Lei X, Lawrence MB, Dong C. Influence of cell deformation on leukocyte rolling adhesion in shear flow. J. Biomech. Eng. 1991;121:636–643. [PubMed]
23. Liang S, Fu C, Wagner D, Guo H, Zhan D, Dong C, Long M. 2D kinetics of β2 integrin-ICAM-1 bindings between neutrophils and melanoma cells. Am. J. Physiol. 2008;294:743–753. [PMC free article] [PubMed]
24. Liang S, Hoskins M, Dong C. Tumor cell extravasation mediated by leukocyte adhesion is shear rate-dependent on IL-8 signaling. Mol. Cell. Biomech. 2009;7:77–91. [PMC free article] [PubMed]
25. Liang S, Hoskins M, Khanna P, Kunz RF, Dong C. Effects of the tumor-leukocyte microenvironment on melanoma–neutrophil adhesion to the endothelium in a shear flow. Cell. Mol. Bioeng. 2008;1:189–200. [PMC free article] [PubMed]
26. Liang S, Slattery M, Dong C. Shear stress and shear rate differentially affect the multi-step process of leukocyte-facilitated melanoma adhesion. Exp. Cell Res. 2005;310:282–292. [PMC free article] [PubMed]
27. Liang S, Slattery M, Wagner D, Simon SI, Dong C. Hydrodynamic shear rate regulates melanoma-leukocyte aggregations, melanoma adhesion to the endothelium and subsequent extravasation. Ann. Biomed. Eng. 2008;36:661–671. [PMC free article] [PubMed]
28. Long M, Goldsmith HL, Tees DF, Zhu C. Probabilistic modeling of shear-induced formation and breakage of doublets cross-linked by receptor–ligand bonds. Biophys. J. 1999;76:1112–1128. [PubMed]
29. Lyczkowski RW, Alevriadou BR, Horner M, Panchal CB, Shroff SG. Application of multiphase computational fluid dynamics to analyze monocyte adhesion. Ann. Biomed. Eng. 2009;37:1516–1533. [PubMed]
30. McQuarrie DA. Kinetics of small systems. I. J. Phys. Chem. 1963;38:433–436.
31. Melder RJ, Munn LL, Yamada S, Ohkubo C, Jain RK. Selectin- and integrin mediated T-lymphocyte rolling and arrest on TNF-activated endothelium: augmentation by erythrocytes. Biophys. J. 1995;69:2131–2138. [PubMed]
32. Munn LL, Melder RJ, Jain RK. Analysis of cell flux in the parallel plate flow chamber: implications for cell capture studies. Biophys. J. 1994;67:889–895. [PubMed]
33. Munn LL, Melder RJ, Jain RK. Role of erythrocytes inleukocyte-endothelial interactions: mathematical model and experimental validation. Biophys. J. 1996;71:466–478. [PubMed]
34. Piper JW, Swerlick RA, Zhu C. Determining force dependence of two-dimensional receptor–ligand binding affinity by centrifugation. Biophys. J. 1998;74:492–513. [PubMed]
35. Rinker KD, Prabhakar V, Truskey GA. Effect of contact time and force on monocyte adhesion to vascular endothelium. Biophys. J. 2001;80:1722–1732. [PubMed]
36. Sabelfeld K, Kolodko A. Stochastic Lagrangian models and algorithms for spatially inhomogeneous Smoluchowski equation. Math. Comput. Simul. 2003;61:115–137.
37. Shankaran H, Neelamegham S. Nonlinear flow affects hydrodynamic forces and neutrophil adhesion rates in cone-plate viscometers. Biophys. J. 2001;80:2631–2648. [PubMed]
38. Shankaran H, Neelamegham S. Effect of secondary flow on biological experiments in the cone-plate viscometer: Methods for estimating collision frequency, wall shear stress and inter-particle interactions in non-linear flow. Biorheology. 2001;38:275–304. [PubMed]
39. Slattery M, Dong C. Neutrophils influence melanoma adhesion and migration under flow conditions. Int. J. Cancer. 2003;106:713–722. [PMC free article] [PubMed]
40. Slattery M, Liang S, Dong C. Distinct role of hydrodynamic shear in PMN-facilitated melanoma cell extravasation. Am. J. Physiol. 2005;288(4):C831–C839. [PMC free article] [PubMed]
41. Smoluchowski M. Mathematical theory of the kinetics of the coagulation of colloidal solutions. Z. Phys. Chem. 1917;92:129.
42. Starkey JR, Liggitt HD, Jones W, Hosick HL. Influence of migratory blood cells on the attachment of tumor cells to vascular endothelium. Int. J. Cancer. 1984;34:535–543. [PubMed]
43. Tees DFJ, Coenen O, Goldsmith HL. Interaction forces between red cells agglutinated by antibody. IV. Time and force dependence of break-up. Biophys. J. 1993;65:1318–1334. [PubMed]
44. Tees DFJ, Goldsmith HL. Kinetics and locus of failure of receptor–ligand mediated adhesion between latex spheres. I. Proteincarbohydrate bond. Biophys. J. 1996;71:1102–1114. [PubMed]
45. Wang J, Slattery M, Hoskins M, Liang S, Dong C, Du Q. Monte Carlo simulation of heterotypic cell aggregation in nonlinear shear flow. Math. Biosci. Eng. 2006;3:683–696. [PubMed]
46. Welch DR, Schissel DJ, Howrey RP, Aeed PA. Tumor-elicited polymorphonuclear cells, in contrast to “normal” circulating polymorphonuclear cells, stimulate invasive and metastatic potentials of rat mammary adnocardinoma cells. Proc. Natl Acad. Sci. USA. 1989;86:5859–5863. [PubMed]
47. Wu QD, Wang JH, Condron C, Bouchier-Hayer D, Redmond HP. Human neutrophils facilitate tumor cell transendothelial migration. Am. J. Phsyiol. Cell Physiol. 2001;280:814–822. [PubMed]
48. Zhang Y, Neelamegham S. Estimating the efficiency of cell capture and arrest in flow chambers: study of neutrophil binding via E-selectin and ICAM-1. Biophys. J. 2002;83:1934–1952. [PubMed]
49. Zhang Y, Neelamegham S. An analysis tool to quantify the efficiency of cell tethering and firm-adhesion in the parallel-plate flow chamber. J. Immunol. Methods. 2003;278:305–317. [PubMed]
50. Zhu C, Bao G, Wang N. Cell mechanics: mechanical response, cell adhesion, and molecular deformation. Annu. Rev. Biomed. Eng. 2000;2:189–226. [PubMed]
51. Zhu C, Chesla SE. Dissociation of individual molecular bonds under force. In: Simon BB, editor. Advances in Bioengineering. Vol. 36. New York: ASME; 1997. pp. 177–178.
52. Zhu C, McEver R. Catch bonds: physical models and biological functions. Mol. Cell. Biomech. 2005;2(3):91–104. [PubMed]