Overview of CTC metastasis in the network

Fig 2 shows simulation snapshots of the CTCs in the microvascular network. At the initial state, there are 13 RBCs and one CTC placed into the microvascular inlet. They are driven into the microvascular network by an acceleration of $$ \overline{\mathit{g}}=(9.79,0,0)$$, which generates a mean flow velocity of $$ {\overline{v}}_m=0.735$$. Here, $$ \overline{\mathit{g}}$$ is scaled by ε′/m′l′, i.e., $$ \overline{\mathit{g}}=\mathit{g}/({\epsilon }^{\prime }/m^{\prime }{l}^{\prime })$$, $$ {\overline{v}}_m=v_m/\sqrt{{\epsilon }^{\prime }/m^{\prime }}$$, and $$ \overline{t}=t/l^{\prime }\sqrt{m^{\prime }/{\epsilon }^{\prime }}$$. Two main phenomena are observed: i) the CTCs are deviated from the centerline of microvessels towards to the wall, known as the CTC margination; ii) the CTCs are more likely to be adhered at the bifurcation.

Fig 2

Snapshots of the CTCs in the microvascular network at different time instants $$ \overline{t}$$.

The RBCs are in red, while the CTCs are in purple and also labeled by Arabic numerals for identification. The black arrows show the adhesion force of the CTCs. For more snapshots, see S1 Video.

Take one of CTCs for example, the CTC-4 in Fig 2. At $$ \overline{t}=429$$, it arrives at the apex of a bifurcation, and a couple of bonds are formed on it. At $$ \overline{t}=497.5$$, it is ready to move into a branch, and more bonds are formed obviously. At $$ \overline{t}=544.5$$, it enters into the branch completely, and all the formed bonds are dissociated. Until $$ \overline{t}=574$$, it is still in the branch and there are no bonds formed any more. The other CTCs also present the similar adhesion behavior, such as the CTC-1, 2, 3 and 5. This implies that the CTCs are more likely to be adhered at the bifurcation region, rather than in the branch.

To support this conclusion, a quantitative analysis is performed by examining the variations of bond number and retention time $$ \Delta \overline{t}$$ with respect to the x−position of each CTC, as shown in Fig 3. It is found that the maximum number of bonds appears at a bifurcation region more often, and the CTC also stays for longer time at the bifurcation region. After the CTC enters into the branch completely, the bond number decreases to zero quickly, and it moves fast due to high velocity in a narrow branch. For example, when the CTC-6 arrives at about $$ {\overline{X}}_c=19$$ (near a bifurcation), about 4 bonds are formed in Fig 3A, and it stays for a time interval of $$ \Delta \overline{t}=30$$, as shown in Fig 3B. Similarly, when the CTC-4 moves to about $$ {\overline{X}}_c=55$$ (near another bifurcation), the number of bonds achieves a peak value of about 200 in Fig 3A, and it stays for a time interval of $$ \Delta \overline{t}=17$$, as shown in Fig 3B. After the CTC passes through the bifurcation, the bond number quickly decreases to zero while the velocity of CTC becoming faster. It is clear that this CTC does not experience a firm adhesion, and thus its adhesion is temporary. Generally, the formation of adhesive bonds at a bifurcation causes the longer time interval during which the CTCs stay at that position due to the adhesion force. In turn, the longer the CTC stays at a position, the larger the probability is to form a bond. As a result, it can be concluded that the tumor cells are indeed more likely to adhere at vascular bifurcations, and may further extravasate from the bifurcations into the tissue for finishing the tumor metastasis. This is consistent with the previous experiments [34], indicating that hemodynamic factors at vascular bifurcations play an important role in tumor metastasis.

Fig 3

Adhesion behavior of the CTCs in the microvascular network.

Variations of (A) the bond number and (B) the physical time with respect to the x-position $$ {\overline{X}}_c$$.

Except the adhesion at a bifurcation, a CTC may also be adhered in a branch, e.g., the CTC-5, the CTC-6 in Fig 2, although this adhesion is not the mainstream. It is calculated that the distance between the CTC-5 centroid and the vessel wall ranges from 3.85 to 1.71 as the vessel diameter ranges from 7.7 to 5.57. That is, the CTC is deviated from the centerline of the vessel towards the wall, such that it arrives at the vessel wall completely. The similar margination behavior occurs with other CTCs, which provides more opportunities for the adhesion. The marginated CTC-5, CTC-6, CTC-7 are observed in Fig 2. Therefore, the margination is often regarded as a prerequisite for the adhesion.

Margination of a CTC in a straight tube

To fully understand the margination of CTCs, we investigate the behavior of one CTC and 49 RBCs in a cylindrical tube with the diameter of $$ {\overline{D}}_v=11$$. The RBC hematocrit is calculated to be H_ct = V_RBC/V_tube = 30%, corresponding to a normal level of RBCs in the real arteriole, where V_RBC and V_tube are the volumes of all RBCs and the tube. The accelerations of $$ \overline{\mathit{g}}=19.15$$ and 5.73 are applied to generate the fluid flow with the different mean velocities of $$ {\overline{v}}_m=2.94$$ and 0.88, respectively. Thus, the average shear rates are calculated to be $$ \overline{\dot{\gamma }}={\overline{v}}_m/{\overline{D}}_v=0.27$$ and 0.08, corresponding to the physical values of 272 s⁻¹ and 81 s⁻¹. They are regarded as a high shear rate for the real arteriolar flow and a low shear rate for the venular flow in microcirculation, respectively [42, 56]. In addition, we use a periodic boundary condition in consideration of the simple straight cylindrical tube, such that a long simulation is allowed to be run. In simulations, the cells pass through the tube for 12 rounds, and thus the effect of the initial configurations of the cells can be negligible.

At the initial state, the CTC is placed at the tube centerline, while the RBCs are evenly distributed throughout the tube, as shown in Fig 4A. As the time elapses, the RBCs gradually migrate towards to the tube centerline, and the RBC-free layer is presented around the tube wall, as shown in Fig 4B. Meanwhile, the CTC is expelled away from the tube centerline, and an obvious margination behavior is observed in Fig 4C. This margination behavior is mainly attributed to the aggregation force acting on the CTC from the surrounding RBCs. Fig 5 shows the deviation of the CTC centroid from the tube centerline, and the radial aggregation force acting on the CTC from the RBCs. It is easily found that they are directly related. When the radial aggregation force is larger than zero, i.e., it is wall-directed, the deviation of the CTC centroid from the tube centerline continuously increases, that is, the CTC is pushed towards the vessel wall. In contrary, when this force is less than zero, i.e., it is center-directed, the CTC moves towards the tube centerline. Taking $$ \overline{\dot{\gamma }}=0.08$$ for example, as $$ {\overline{X}}_c\in [0,70]$$, the radial aggregation force is larger than zero, and thus the CTC continuously deviates from the tube centerline. As $$ {\overline{X}}_c\in [70,120]$$, the radial aggregation force is smaller than zero, and thus the CTC moves towards to the tube centerline.

Fig 4

Snapshots of the CTC margination in a straight tube.

These snapshots represent the motion states of RBCs and the CTC in the margination process. (A) The initial state, (B) RBCs move toward the center and (C) the CTC moves toward the wall.

Fig 5

Effect of the interaction force between the CTC and RBCs on margination behavior.

(A) The deviation $$ \Delta {\overline{r}}_c$$ of the CTC centroid from the tube centerline at the shear rates of $$ \overline{\dot{\gamma }}=0.08$$ and 0.27. (B) The radial aggregation force $$ {\overline{\mathit{F}}}^{agg}$$ acting on the CTC from the RBCs, where the positive is the force directing to the tube wall, and the negative is the force directing to the tube centerline.

Another phenomenon is also observed from Fig 5 that as the CTC gets closer to the wall, it is pushed back towards the tube centerline to some extent. Hence, the CTC cannot, in fact, approach to the vessel wall enough to exhibit the adhesion behavior, albeit the margination observed. We call this the incomplete margination, for example, the margination of CTC at $$ {\overline{X}}_c=70$$ under $$ \overline{\dot{\gamma }}=0.08$$, and the margination of CTC at $$ {\overline{X}}_c=190$$ under $$ \overline{\dot{\gamma }}=0.27$$. This is attributed to the local interaction between the CTC and the RBCs near the wall. As the CTC is close enough to the vessel wall, it is still surrounded by several RBCs. These RBCs provide a relatively strong center-directed aggregation force, because the distance between the CTC and RBCs is small enough at this moment. This force pushes the CTC back to the tube centerline to some extent.

In addition, it is also found from Fig 5 that the wall-directed force at $$ \overline{\dot{\gamma }}=0.08$$ is less than that at $$ \overline{\dot{\gamma }}=0.27$$ on the whole, but the CTC margination at $$ \overline{\dot{\gamma }}=0.08$$ is more obvious. This is because the wall-directed force acts on the CTC much longer at the low shear rate than at the high shear rate. As a result, the CTC margination often occurs in the venules where the blood flows slowly although the wall-directed force is small, rather than the arterioles with a high shear rate. That is one of the reasons why the CTC adhesion is more likely to occur in the venules and rarely in the arterioles [34]. Certainly, if the vessel is narrow enough, such as the capillary in Fig 2D, the CTC can reach the microvascular wall completely, and its adhesion may happen.

Adhesion of a deformable capsule on a plate

To fully understand the adhesion behavior and establish reliability of adhesion model, we investigate the adhesion behaviors of a deformable capsule, and compare the simulation results with the work of Zhang et al. [55]. The capsule is spherical with a radius of R = 3.75μm, which is placed into a cubic tube with the size 10R × 6R × 6R. A Couette flow is generated by moving the top wall of the tube with a constant shear rate $$ \dot{\gamma }=7000{\mathrm{s}}^{-1}$$. The fluid density is ρ = 10³ kg/m³, and the viscosity is η = 10⁻³ pa ⋅ s. The Reynolds number is Re = 0.1. The capillary number Ca is set to be 0.005 and 0.015, such that the shear modulus of the capsule is E_S = 5.25 × 10⁻³ and 1.75 × 10⁻³ N/m. The bending modulus is calculated by E_B = 0.02R²E_S, and the dilation modulus is calculated by E_D = 100E_S. In the adhesion model, the equilibrium length is λ = 50 nm, and the reactive distance is l_r = 375 nm. The bond strength is determined by a dimensionless parameter $$ K_{SP}=E_A/\eta R\dot{\gamma }=250$$, that is, E_A = 6.56 × 10⁻³ N/m. The formation strength is σ_f = 0.02E_A, and the dissociation strength is σ_d = 0.98E_A. The unstressed formation and dissociation rates are determined by two dimensionless parameters, $$ K_f=k_f^0/\dot{\gamma }=10$$ and $$ K_d=k_d^0/\dot{\gamma }=$$ 1.0 (or 0.01), respectively. All of these parameters are set to be same as those used by Zhang et al. [55]

Fig 6 shows the three distinct adhesion states under Re = 0.1. At a high shear rate (Ca = 0.015) and a high dissociation rate (K_d = 1.0), the capsule detaches completely from the bottom wall and moves along the fluid flow. This state is called the detachment adhesion, as shown in Fig 6A. As the shear rate decreases (Ca = 0.005) but the dissociation rate remains unchanged (K_d = 1.0), the capsule rolls on the bottom wall, and this state is called the rolling adhesion, as shown in Fig 6B. Finally, if the dissociation rate decreases (K_d = 0.01) but the shear rate remains unchanged (Ca = 0.015), the capsule is firmly adhered onto the wall and does not move any more. This state is called the firm adhesion, as shown in Fig 6C. It is clear that our results show the same three states as those of Zhang et al. [55].

Fig 6

Comparisons of the three adhesion states.

Comparisons of the three distinct adhesion states between our results (A-C) and those obtained by Zhang et al. [55] (D-F) under Re = 0.1. The detachment state (A and D) is obtained under Ca = 0.015 and K_d = 1.0; the rolling adhesion state (B and E) is obtained under Ca = 0.005 and K_d = 1.0; the firm adhesion state (C and F) is obtained under Ca = 0.015 and K_d = 0.01.

Which adhesion state is present depends on the number of bonds formed between the receptors on the cell membrane and the ligands on the bottom wall. Fig 7 shows the time evolutions of the bond number and translational velocity in the three adhesion states. In the detachment adhesion, there are a few bonds formed at the beginning, but they are all dissociated quickly. This causes the capsule to exhibit a translational velocity that decreases first and then increases to move with fluid flow. In the firm adhesion, there are about 70 bonds formed at the steady state, which hold the capsule firmly on the bottom wall with the zero translational velocity. In the rolling adhesion, there are around 20 bonds that are kept for the whole simulation, leading to a translational velocity fluctuating around a constant value. To sum up, the rate of forming the new bonds is slower than that of dissociating the existing bonds in the detachment adhesion, but the opposite is for the firm adhesion. Both the rates are in equilibrium in the rolling adhesion.

Fig 7

Comparisons of (A) bond number and (B) translational velocity in the three adhesion states.

In addition, a difference is noted that the bond number in our work is not exactly same as that in the work of Zhang et al. [55], as observed from Fig 7A. This difference is attributed to the difference of the deformation models used. A discrete particle model is used in our work, but Zhang et al. [55] adopted the Mooney-Rivlin model. Both models are approximated as linear at a small deformation, whereas at a large deformation, our cell has a larger deformation force than theirs if they undergoes the same deformation. That is, our model is strain hardening, while the model used by Zhang et al. [55] is strain softening. Therefore, the capsule in our work looks more rigid in each adhesion state, as shown in Fig 6, and thus they undergo the larger translational velocity in our work, as shown in Fig 7B. In the detachment and rolling adhesion, the rigid capsule is more likely to dissociate the bonds at a high dissociation rate K_d = 1.0, such that the bond number in our work is less than that in the work of Zhang et al. [55] In the firm adhesion, however, the rigid capsule is observed to have the larger contact surface with the bottom wall when its rear is firmly adhered on the bottom wall with K_d = 0.01, leading to more bonds formed in our work than that given by Zhang et al. [55] It can be found from this analysis that the cell deformation, cell adhesion and shear fluid are entangled and interact with each other.

Adhesion of a CTC in a bifurcation

As mentioned above, the CTC adhesion happens more often in a bifurcation region, and hence we here focus on its behaviors in a ‘single’ bifurcation for a deeper understanding. The bifurcation is directly extracted from the microvascular network in Fig 1B, and the effects of RBC hematocrit and shear rate of fluid on the CTC adhesion are investigated in this section.

Fig 8 shows the effect of the RBC hematocrit, which is adjusted to be H_ct = 10–40% by varying the number of RBCs in the microvascular inlet. The fluid flow has a mean velocity of $$ {\overline{v}}_m$$ of 0.735 in all these cases. It is found that the bond number at H_ct ≥ 30% is obviously greater than that at H_ct < 30%. Its peak value is larger than 11 for H_ct ≥ 30%, but between 4 to 7 for H_ct < 30%, as shown in Fig 8B. This implies that the CTC has the strongest adhesion at H_ct = 30%, and this is also supported by the radial aggregation force that achieves the maximum value at H_ct = 30%, as shown in Fig 8C. This result seems to contradict the expectation that an increase of H_ct can enhance the CTC adhesion because of the increasing repulsion from the more RBCs. In fact, this expectation is not true for all the hematocrits. From a low hematocrit, such as H_ct = 10–25%, the radial aggregation force indeed increases as H_ct increases to 30%, as shown in Fig 8C, which in turn promotes the CTC adhesion. This was also observed in the work of Xiao et al. [27], where H_ct ranging from 10 to 30% was considered. However, when the hematocrit varies from a normal value to a higher value, for example from H_ct = 30 to 40%, the radial aggregation force does not increase but instead decreases, and correspondingly the CTC adhesion becomes weak. This is because the RBC-free layer is narrowed down subject to the more RBCs at a high hematocrit. Thus, the CTC is more likely to be surrounded by those RBCs, as shown in Fig 8A, and the outer RBCs give the CTC a center-directed aggregation force, counteracting the wall-directed aggregation force acting on the CTC. This was observed for the leukocytes in the work of Fedosov et al. [12] To sum up, the low hematocrit cannot provide the enough wall-directed aggregation force to promote the CTC adhesion, because of the small RBC number. The high hematocrit cannot also promote the CTC adhesion because the wall-directed aggregation force is offset by the more RBCs surrounding the CTC.

Fig 8

Effect of the RBC hematocrit on the CTC adhesion.

(A) The snapshots of the CTC at about $$ {\overline{X}}_c=54$$, the variations of (B) the bond number and (C) the radial aggregation force from the RBCs, at the RBC hematocrits of H_ct = 10–40%.

Fig 9 shows the effect of the shear rate of fluid flow, which is adjusted to be $$ \overline{\dot{\gamma }}=0.048$$, 0.095 and 0.19 by varying the externally-applied acceleration. The fluid flow has the mean velocity of $$ {\overline{v}}_m=0.37$$, 0.735, 1.47, respectively, and the RBC hematocrit is fixed as H_ct = 20%. A conclusion is drawn that the CTC adhesion becomes weak as the shear rate increases, which was also confirmed in the recent simulation work of Dabagh et al. [26]. For example, at about $$ {\overline{X}}_c=54$$, there are 13, 6 and 3 bonds formed for $$ \overline{\dot{\gamma }}=0.048$$, 0.095 and 0.19, respectively. This is directly attributed to the increase of the translational velocity of the CTC with increasing the shear rate, as shown in Fig 9C. At a low shear rate, $$ \overline{\dot{\gamma }}=0.048$$, the cells move slow so that the RBCs is aggregated obviously, as shown in Fig 9A, leading to a strong wall-directed aggregation force to promote the CTC adhesion. At a high shear rate, $$ \overline{\dot{\gamma }}=0.19$$, the RBCs are distributed more dispersedly, and the less wall-directed aggregation force is provided for the CTC adhesion. Moreover, it is worth mentioning that the CTC obtains the maximum bond number at the bifurcation (about $$ {\overline{X}}_c=54$$), regardless of discussing the effects of hematocrit or shear rate. This again confirms that the CTC is more likely to be adhered at the position of a bifurcation. The similar phenomena have been observed in experiments. Liu et al. [37] found that MDA-MB-231 human breast cancer cells adhere at the microvascular bifurcations, and then extravasate from the blood vessel, by an in vivo experiment. Guo et al. [34] also observed that the tumor cells may more often adhere at the bifurcated regions of vessels, as well as small venules with relatively low shear rates.

Fig 9

Effect of the shear rate of fluid on the CTC adhesion.

(A) The snapshots of the CTC at about $$ {\overline{X}}_c=54$$, the variations of (B) the bond number, and (C) the translational velocity, at the shear rates of $$ \overline{\dot{\gamma }}=0.048$$, 0.095 and 0.19.

SDPD-IBM model for fluid flow

The Navier-Stokes (N-S) equations are adopted to govern the motion of the plasma and cytosol, given by [57]

The SDPD-IBM method is used to discretize Eqs 2 and 4, and it has been already developed well in our previous work [46, 59], and here we only outline its framework for completeness. The fluid particles are evolved by

Discrete elastic model for cell deformation

To characterize the cell deformation, the cell membrane is modeled as a triangular network by connecting all membrane particles. The edges of each triangle are modeled as springs to describe the membrane elasticity, and the angles of any two neighboring triangles are used to describe the membrane bending. Moreover, the membrane area is restrained to be varied in 3%, due to its strong membrane elasticity [60]. The cell volume is also restrained to be varied in 3%, because the material exchange is often not considered in the previous work [61–63]. As a result, the total deformation potential energy U^def of the triangular network is given by [64, 65]

The deformation model is associated with there important macro parameters: the shear modulus E_S, the bending modulus E_B, and the dilation modulus E_D [64],

Morse potential model for cell aggregation

To describe the cell-cell interaction, the Morse potential model proposed by Liu et al. [66] is used, in which the total aggregation energy between the cells is approximated as [67],

Stochastic binding model for cell adhesion

To characterize the CTC adhesion, the stochastic binding model proposed by Hammer and Apte [28] is used, in which an adhesion potential energy is given by,

The formation and dissociation of a bond are stochastic processes, controlled by a formation probability p_f and a dissociation probability p_d, respectively [28],

Numerical details

There are two types of boundary conditions for Eqs 5–8, the solid boundary condition, and inflow/outflow boundary condition. The former one is applied on the microvascular wall for two purposes. One is to improve the accuracy near the microvascular wall by introducing the ghost particles. The other is to avoid the fluid and membrane particles to penetrating the microvascular wall by introducing the repulsive particles. More details about this boundary condition can be found in the work of Liu et al. [69] The inflow/outflow boundary condition is applied in the inlet/outlet of microvascular network, also for two purposes. One is to keep cells continuously moving into the microvascular network, and the other is to settle down the fluid and membrane particles that move out from the microvascular network. This boundary condition is one of advantages of our model, which not only ensures that the system including cells and fluid have reached steady states when they flow into the microvascular network, but also guarantees that the mass and momentum of the system are conserved. More details about this boundary condition can be found in our previous work [59].

A non-dimensional procedure is carried out, by choosing the cut-off radius, particle mass and Boltzmann temperature as the characteristic length l′, mass m′ and energy ϵ′. They are set as: l′ = 2 μm, m′ = 1.0 × 10⁻¹⁵ kg, and ϵ′ = 4.142 × 10⁻²¹ J, respectively. The other physical parameters can be scaled by these three parameters; for example, the shear modulus E_S is scaled by ε′/l′², and the bending modulus E_B is directly scaled by ε′. Table 1 lists the physical parameters and their corresponding simulation values used in the present work. In the present work, we add an overline on the physical quantity ($$ \overline{·}$$) to stand for its dimensionless form. After that, the velocity-Verlet algorithm [70] is employed to numerically solve SDPD-IBM model in Eqs 5–8, because of its advantages of high computational efficiency. This algorithm is parallelized by the technique of message passing interface (MPI) to reduce the computational time. The detailed procedure about this algorithm can be found in our previous work [59]. The computational cost is about 51,600 core ⋅ hour for a typical simulation case.

Table 1

Physical quantities and their characteristic quantities.

Physical quantities	Physical values	Characteristic quantities	Simulation values
Fluid density (ρ)	1.0 × 103kg/m3	m′/l′3	8
Shear viscosity (η)	1.0 × 10−4pa ⋅ s	$$	197
Temperature (T)	300K	ϵ′	1
CTC diameter (Dt)	9.0μm [26]	l′	4.5
RBC diameter (Dr)	7.82μm [50]	l′	3.91
Shear modulus of CTC (ES)	1.0 × 10−6N/m [71]	ϵ′/l′2	9.66 × 102
Shear modulus of RBC (ES)	6.0 × 10−6N/m [72]	ϵ′/l′2	5.794 × 103
Bending modulus of CTC (EB)	1.35 × 10−19J [73]	ϵ′	33
Bending modulus of RBC (EB)	2.07 × 10−19J [72]	ϵ′	50
Dilation modulus of CTC (ED)	1.16 × 10−4N/m†	ϵ′/l′2	1.12 × 105
Dilation modulus of RBC (ED)	1.26 × 10−4N/m†	ϵ′/l′2	1.22 × 105
Equilibrium length (λ)	0.2μm [55]	l′	0.1
Reactive distance (lr)	1.0μm [55]	l′	0.5
Adhesion strength (EA)	3.0 × 10−5N/m [55]	ϵ′/l′2	2.9 × 104
Formation strength (σf)	6.0 × 10−7N/m [55]	ϵ′/l′2	5.8 × 102
Dissociation strength (σd)	7.54 × 10−7N/m [55]	ϵ′/l′2	7.28 × 102
Unstressed formation rate ($$)	1.205 × 106s−1 [55]	$$	1.205 × 103
Unstressed dissociation rate ($$)	3.55 × 105s−1 [55]	$$	3.55 × 102

^† The dilation modulus in simulations is, in general, smaller than real values to save computational cost, as long as it can guarantee that the variation of surface area of a RBC or tumor cell is less than 3%.