|Home | About | Journals | Submit | Contact Us | Français|
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.
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,23–25 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,25–27 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,14–16,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.
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.
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
and would decrease
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; 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. 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).
The collision rate, 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).
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 . 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 height H was determined by the shear rate and medium viscosity μ according to H = f(, μ), 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.
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; was the forward reaction coefficient per unit density with n bonds connecting two cells and 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; 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 e−A 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*(, μ), 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 and viscosity μ (Table 2). Then, the probability of aggregation corresponded to the probability of having more than n*(, μ) bonds at the equilibrium state. Notice that n* was not necessarily a whole number.
We proposed two formulations to evaluate the adhesion efficiency, both of which depended only on A and n*:
Here, was the smallest integer bigger than n, and n̠ 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.
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 and 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 and the capillary number Ca. The parameters and data used to fit the tethering frequency Tf are listed in Table 3.
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:
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.
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 , 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 , which was assumed to be zero, were collected either from literature 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.
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 , 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 was more dramatic in comparison with the increase in tethering frequency, it would promote aggregation and vice versa.
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(, μ) 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:
Nonlinear regression was conducted with the help of least square method to predict the relation among the tethering frequency (Tf), shear rate () 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 , 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.
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.
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.
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
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.
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.
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 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 is multiplied by 4 with Tf being fixed (Fig. 12b); while in case III, only Tf is multiplied by 4 and being fixed (Fig. 12c).
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.
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.
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.
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.
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 ) 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 .
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 , the coagulation kernel βPT, the tumor cell concentration CT, and a decreasing function of the tethering frequency Tf, which agreed with our numerical simulations.
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.
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
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 . Assume that the deformable PMN preserves its volume , where rp is the radius of an undeformed PMN, we thus have that
which gives a closed form for H, r and 0 via the following explicit formulae:
Now consider a tumor cell which is colliding with this PMN. Assume that the contact point has coordinate (r + rt, Θ, ), 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 and μ, we compute H by the fitted function f in Eq. (16) first, and compute 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.
The algebraic system
has a general solution , 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 . 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.