Margination and adhesion dynamics of tumor cells in a real microvascular network

In tumor metastasis, the margination and adhesion of tumor cells are two critical and closely related steps, which may determine the destination where the tumor cells extravasate to. We performed a direct three-dimensional simulation on the behaviors of the tumor cells in a real microvascular network, by a hybrid method of the smoothed dissipative particle dynamics and immersed boundary method (SDPD-IBM). The tumor cells are found to adhere at the microvascular bifurcations more frequently, and there is a positive correlation between the adhesion of the tumor cells and the wall-directed force from the surrounding red blood cells (RBCs). The larger the wall-directed force is, the closer the tumor cells are marginated towards the wall, and the higher the probability of adhesion behavior happen is. A relatively low or high hematocrit can help to prevent the adhesion of tumor cells, and similarly, increasing the shear rate of blood flow can serve the same purpose. These results suggest that the tumor cells may be more likely to extravasate at the microvascular bifurcations if the blood flow is slow and the hematocrit is moderate.

Introduction Tumor metastasis is the major cause of cancer treatment failure, where tumor cells detach from the primary site, and then spread to distant organs or other parts of the body through blood and lymphatic systems, forming the secondary tumor [1][2][3]. The metastasis through the blood circulation is known as hematogenous metastasis, which was reported to cause 90% death of cancer [4], and the tumor cells entering into blood circulation are so-called the circulating tumor cells (CTCs). There are several critical steps for finishing a tumor metastasis, including intravasation, margination towards the vascular wall, adhesion onto the vascular wall, and extravasation from the circulation to a distant host organ [5]. We here pay our attention to the margination and adhesion of the CTCs in a real microvascular network, for providing a deep understanding of the tumor metastasis.
The margination and adhesion are two critical and closely related steps in tumor metastasis, and generally, the former step is often regarded as a prerequisite for the latter one. Both of behaviors were firstly observed on the leukocytes [6], instead of tumor cells, and hence, most of the subsequent studies have still focused on the leukocytes [7][8][9][10][11][12][13][14]. For example, Firrell et al. [7] and Abbitt et al. [9] suggested that the margination and adhesion of a leukocyte mainly occur at low flow rates based on in vivo and in vitro experiments. Jain et al. [11] demonstrated in the microfluidic experiments that the leukocytes have obvious margination within a range of hematocrit of 20-30%, and the lower or higher hematocrit may cause less margination. Fedosov et al. [12][13][14] also gave the similar dependence of leukocyte margination on the flow rate and hematocrit by numerical simulations. Some other types of cells, such as platelets [15][16][17][18][19] and malaria-infected RBCs [20,21], have been found to have similar margination in the vessel. So far, there are few quantitative studies on the margination and adhesion of tumor cells. King et al. [22] pointed out that the soft CTCs are marginated more slowly than the rigid CTCs, because larger deformation of the soft CTCs caused by the collisions with surrounding RBCs can prolong the time that CTCs arrive at the vessel wall. However, some experiments [23][24][25] have shown that the soft CTCs are more likely to extravasate from the circulation, because the soft CTCs have a large contact surface with the vascular wall and thus have a higher possibility of adhesion than the rigid CTCs. Recently, some simulation work [26,27] has also favored the latter conclusion. A stochastic adhesive model was developed by Hammer and Apte [28], which has been widely used to theoretically and numerically investigate the adhesion of various cells, certainly including the CTCs [27]. The model proposed that the cell adhesion is attributed to the dynamic association and dissociation between the receptors on the cell and the ligands on the vessel wall. Obviously, if the receptors and ligands contact with each other frequently, the cell is more likely to adhere to the vessel wall. But, the contact frequency depends on a lot of parameters, such as the numbers of the receptors and ligands, the shear force of fluid flow, and local force from neighboring cells [5,[29][30][31][32][33]. Apart from these, some studies [34][35][36] have shown that the cell adhesion has a considerable dependence on the microvascular geometry. For example, Liu et al. [37] observed in vivo experiments that breast cancer cells preferentially extravasate from the vessel bifurcation. More recently, Hynes et al. [36] and Pepona et al. [38] studied the metastatic behavior within a complex vasculature via both experiments and simulations, and also confirmed that CTCs are prone to attach at the branch points. However, so far there are few studies on the tumor metastasis in a complex microvascular network at a cellular scale, at which the cell deformation, aggregation and adhesion should be taken into account simultaneously [39][40][41]. Experimental observations at this scale are limited by the reliability of measurements and the complexity of blood flow in microvascular networks [42,43]; numerical modeling also poses great challenges due to insects. The microvascular network is filled only with plasma having no cells at the initial state. In order to make cells flow into the network, a cylindrical inlet with 120 μm in length is added at its left side, and it is arranged by a large number of cells, including RBCs and CTCs. An external force is imposed only in the inlet to generate an inflow, which pushes the cells into the network. Meanwhile, a cylindrical outlet with 20 μm in length is added at its right side for handling the cells that move out from the network. Hence, the total length of the simulation domain is actually 400 μm. The RBC is modeled to be a biconcave with a maximum diameter of 7.82 μm [50], while the CTC is spherical with the diameter of 9 μm. Due to the diversity of CTCs, their sizes are very different. In experiments, Gassmann et al. [51] observed that the diameter of a tumor cell in vivo is only about 5 μm, while Guo et al. [34] gave its diameter about 16 μm. In simulations, a moderate value of 8-10 μm is often chosen for saving computational cost [26,27,52]. The CTC's nucleus is also neglected for the same purpose, saving the computational cost, and this is one of limitations of our cell model. Fortunately, its effect on the CTC mechanics can be compensated to some extent by enhancing the CTC stiffness [53,54]. The SDPD-IBM model, a particle-based method, is used as a numerical method in the present work, which discretizes the computational domain into a lot of particles. They are classified into six types in total, ghost particles, repulsive particles, fluid particles, membrane particles, receptor sites and ligand sites, respectively, as shown in Fig 1C. The ghost and repulsive particles are used to model the solid wall of the microvessels; the fluid particles are used to model the plasma and cell cytosol; the membrane particles are used to model the cell membrane; the receptor sites are scattered on the cell membrane, while the ligand sites are on the microvascular wall. The receptor and ligand sites are bonded for describing the CTC adhesion on the vascular wall. In the present work, there are 29704 ligand sites on the wall, 1176 receptor sites on the CTC's membrane, and no receptor sites on the RBC's membrane. It should be noted that these numbers are not same as the experimental values due to the limitation of computational resources. Hence, a coarse-grained technique is used, that is, a numerical receptor or ligand may represent several tens or hundreds of experimental receptors or ligands. The number of 29704 ligands is calculated by the fixed number density of 3 μm −2 , similar with that used in the work of Zhang et al. [55], while the number of 1176 receptors is similar with that used the work of Xiao et al. [27]. The ghost particles are 93116, the repulsive particles are 181810, and the fluid particles are 88771, respectively. Each RBC has 613 membrane particles, and each CTC has 1176 membrane particles equal to receptor sites for simplicity. . 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.

Overview of CTC metastasis in the network
Take one of CTCs for example, the CTC-4 in Fig 2. At � t ¼ 429, it arrives at the apex of a bifurcation, and a couple of bonds are formed on it. At � t ¼ 497:5, it is ready to move into a branch, and more bonds are formed obviously. At � t ¼ 544:5, it enters into the branch completely, and all the formed bonds are dissociated. Until � 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 D � 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 � X c ¼ 19 (near a bifurcation), about 4 bonds are formed in Fig 3A, and it stays for a time interval of D � t ¼ 30, as shown in Fig 3B. Similarly, when the CTC-4 moves to about � 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 D � 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.
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 � 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 � g ¼ 19:15 and 5.73 are applied to generate the fluid flow with the different mean velocities of � v m ¼ 2:94 and 0.88, respectively. Thus, the average shear rates are calculated to be D v ¼ 0:27 and 0.08, corresponding to the physical values of 272 s −1 and 81 s −1 . 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

PLOS COMPUTATIONAL BIOLOGY
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 � _ g ¼ 0:08 for example, as � X c 2 ½0; 70�, the radial aggregation force is larger than zero, and thus the CTC continuously deviates from the tube centerline. As � X c 2 ½70; 120�, the radial aggregation force is smaller than zero, and thus the CTC moves towards 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 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 � _ g ¼ 0:08 is less than that at � _ g ¼ 0:27 on the whole, but the CTC margination at � _ g ¼ 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 _ g ¼ 7000 s À 1 . The fluid density is ρ = 10 3 kg/ m 3 , and the viscosity is η = 10 −3 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 −3 and 1.75 × 10 −3 N/m. The bending modulus is calculated by E B = 0.02R 2 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 The unstressed formation and dissociation rates are determined by two dimensionless parameters, 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  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].
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.
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 � 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 walldirected 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 9 shows the effect of the shear rate of fluid flow, which is adjusted to be � _ g ¼ 0:048, 0.095 and 0.19 by varying the externally-applied acceleration. The fluid flow has the mean velocity of � 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 � X c ¼ 54, there are 13, 6 and 3 bonds formed for � _ g ¼ 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, � _ g ¼ 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, � _ g ¼ 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 � 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.

Discussion
A direct three-dimensional simulation of tumor cells in a complex microvascular network was carried out to understand the metastasis of tumor cells as realistic as possible. The microvascular network was constructed from a real mesenteric vasculature of a rat, and comprised of bifurcating, merging and winding vessels. The cells include the red blood cells (RBCs) and circulating tumor cells (CTCs), without the platelets due to their small portion in blood. Three mechanical behaviors of cells, the deformation, aggregation and adhesion, were taken into account in the present work.
The margination and adhesion of tumor cells were mainly investigated, as well as the effects of the RBC hematocrit and flow shear rate. The simulation results showed that the tumor cells adhere more frequently at the microvascular bifurcations, and there is a positive correlation between the adhesion and wall-directed force from the surrounding RBCs. The larger the walldirected force is, the closer the tumor cells are marginated towards the wall, and the higher the probability of adhesion behavior happen is. At a low hematocrit, e.g., 10%, the RBCs cannot provide the enough wall-directed force, because of the small RBC number, while at a high hematocrit, e.g., 40%, the wall-directed force is offset by the outer RBCs surrounding the CTC, because more RBCs exist in the vessels. Hence, a moderate hematocrit of 30% that is close to the normal level of hematocrit in a real arteriole, is found to have the largest wall-directed force, and have the largest bond number. Moreover, it is also found that the CTC margination and adhesion are enhanced as the shear rate decreases. At a low shear rate, the cells are transported slowly and aggregated easily, and the wall-directed force lasts long on a CTC. Hence, the CTC margination is more obvious at the low shear rate than at the high shear rate, and the bond number is larger, promoting the CTC adhesion.
To sum up, the present work suggests that the tumor cells may be more likely to adhere at the microvascular bifurcations with low shear rate and moderate hematocrit, because of a high wall-directed force from the surrounding RBCs. This may help to predict the location where tumor cells extravasate from the circulation, and further give new insights into the cancer pathophysiology and its diagnosis and therapy.

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] r � v ¼ 0; ð1Þ where ρ, v and t are the density, velocity and time, respectively; P is the pressure field, η is the shear viscosity, g is an acceleration externally applied to drive the fluid flow, and f is a singular force from the cell membrane. The interaction between the cell and fluid is modeled by the immersed boundary method (IBM). The action on the fluid from the cell is expressed as while the action on the cell membrane from the fluid is shown by the membrane evolution, i.e., where f cell is the force acting on the cell membrane due to the cell deformation, aggregation and adhesion, etc. x is the position of the fluid in the Eulerian system, X is the position of the cell membrane in the Lagrangian system, and (r, s) denotes a curvilinear coordinate on the cell membrane to label a Lagrangian point. Γ is the surface domain occupied by the cell membrane, and O is the computational domain. It should be here pointed out that both the cytosol and plasma are assumed to be a single type of fluid with the same physical properties, and they are incompressible and isothermal Newtonian fluid. The whole blood consisting of the plasma and the cells behaves as non-Newtonian fluid, because the cells are incorporated by the singular force f, and thus the blood viscosity is described to depend on the cell hematocrit [58]. 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 where m, v i and x i are the mass, velocity and position of the fluid particle i, respectively. The force F C is the conservative force describing the compressibility of fluid. F D is the dissipative force describing the viscosity of fluid. F R is the random force to keep a constant Boltzmann temperature of fluid. See S1 Appendix for the detailed formulations of these three forces. F G is an externally-applied force to drive the fluid flow. F B is the force from the cell membrane, defined as where β ik is the weighted coefficient determined by the kernel function of SDPD, and F def k , F agg k and F adh k are the membrane forces due to the cell deformation, aggregation and adhesion, respectively. The membrane particles are evolved by where X k is the position of the membrane particle k. The ghost and repulsive particles are stationary during the whole simulation, as well as the ligand sites. However, the positions of the receptor sites are evolved with the membrane particles, because the membrane particles are assumed to be same as the receptor sites.

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] and the deformation force is calculated as where U s , U b , U a and U v denote the in-plane energy, bending energy, area-restraint energy and volume-restraint energy, respectively. See S1 Appendix for more details about these four energy.
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], where k B T and p j are the Boltzmann temperature and persistence length, respectively. l 0 j is the length of the triangular edge j at the stress-free state, s 0 ¼ l 0 j =l d j , l d j is the maximum length of the triangular edge j. K B , K AG and K AL are the bending coefficient, global area restraint constant and local area restraint constant, respectively.

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], and the aggregation force is given by,

PLOS COMPUTATIONAL BIOLOGY
where N t is the number of triangles, and A m is the area of the triangle m. φ(r mm 0 ) is the Morse potential energy between two facing plane elements of the separate cells, defined as [66] φðr mm 0 Þ ¼ E I ½e 2bðr 0 À r mm 0 Þ À 2e bðr 0 À r mm 0 Þ �; ð16Þ where r mm 0 is the local distance between the two elements, E I is the strength of surface energy, r 0 is the zero-force distance, and β is a scaling factor. The term of (n m � k m )(n m 0 � k m 0 ) is added to consider the effect of curved elements instead of plane elements by the DLVO theory [68], where n is the outward unit normal vector of a triangle, and k is the unit vector parallel to the line connecting the centers of two interacting cells. This model behaves as a weak attractive force at a far distance (r mm 0 > r 0 ), whereas at a near distance (r mm 0 < r 0 ) it behaves as a strong repulsive force.

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, and the adhesion force is thus calculated by where E A is the adhesion strength, N b is the number of formed bonds, x m is the bond length, and λ is an equilibrium distance. When a new bond is formed, it contributes to the adhesion potential energy. On the contrary, when an existing bond is dissociated, the corresponding contribution should be eliminated. 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], where Δt is the time step. k f and k d are the formation and dissociation rates of a bond, defined as and where k 0 f , k 0 d , σ f and σ d are four constants, and k B T is Boltzmann temperature. k 0 f and k 0 d are the unstressed formation and dissociation rates of bonds at the equilibrium distance λ; σ f and σ d are the formation and dissociation strengths within a given reactive distance l r . If a bond is formed, the following two conditions must be satisfied, where ξ 1 is a random number with uniform distribution in [0, 1]. If an existing bond is dissociated, one of the following conditions needs to be satisfied, where ξ 2 is another random number similar to ξ 1 .

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 0 , mass m 0 and energy � 0 . They are set as: l 0 = 2 μm, m 0 = 1.0 × 10 −15 kg, and � 0 = 4.142 × 10 −21 J, respectively. The other physical parameters can be scaled by these three parameters; for example, the shear modulus E S is scaled by ε 0 /l 02 , and the bending modulus E B is directly scaled by ε 0 . 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 (� �) 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.
Supporting information S1 Video. Metastasis of tumor cells in the microvascular network.