Unraveling the vascular fate of deformable circulating tumor cells via a hierarchical computational model

TL;DR

This study presents a hierarchical computational model combining fluid dynamics and cell mechanics to predict how circulating tumor cells of varying deformability interact with blood vessels, influencing metastasis.

Contribution

It introduces a novel coupled Lattice-Boltzmann and Immersed Boundary Method model to simulate CTC vascular transport and adhesion based on cell deformability.

Findings

Soft CTCs behave like RBCs in large vessels and tend to adhere firmly.

Rigid CTCs are pushed away from vessel walls in large vessels.

In small vessels, soft CTCs can rapidly establish firm adhesion.

Abstract

Distant spreading of primary lesions is modulated by the vascular dynamics of circulating tumor cells (CTCs) and their ability to establish metastatic niches. While the mechanisms regulating CTC homing in specific tissues are yet to be elucidated, it is well documented that CTCs possess different size, biological properties and deformability. A computational model is presented to predict the vascular transport and adhesion of CTCs in whole blood. A Lattice-Boltzmann method, which is employed to solve the Navier-Stokes equation for the plasma flow, is coupled with an Immersed Boundary Method. The vascular dynamics of a CTC is assessed in large and small microcapillaries. The CTC shear modulus k ctc is varied returning CTCs that are stiffer, softer and equally deformable as compared to RBCs. In large microcapillaries, soft CTCs behave similarly to RBCs and move away from the vessel walls; whereas rigid CTCs are pushed laterally by the fast moving RBCs and interact with the vessel walls. Three adhesion behaviors are observed, firm adhesion, rolling and crawling over the vessel walls, depending on the CTC stiffness. On the contrary, in small microcapillaries, rigid CTCs are pushed downstream by a compact train of RBCs and cannot establish any firm interaction with the vessel walls; whereas soft CTCs are squeezed between the vessel wall and the RBC train and rapidly establish firm adhesion. These findings document the relevance of cell deformability in CTC vascular adhesion and provide insights on the mechanisms regulating metastasis formation in different vascular districts.

Tables

Table 1. Table 1: Evolution of the Taylor parameter D D \rm D as a function of the capillary number Ca Ca \rm Ca .

Ca	$D$ Lac (2004)	$D$ $\left(r=10,L=160\right)$	$\%$ Diff.	$D$ $\left(r=5,L=80\right)$	$\%$ Diff.
0.0375	0.0798	0.0794	0.04	0.0790	0.080
0.075	0.1594	0.1487	1.07	0.1461	1.33
0.15	0.2718	0.2594	1.24	0.2491	2.49
0.3	0.4053	0.3904	1.49	0.3749	3.04

Table 2. Table 2: Evolution of the inclination angle ϕ / π italic-ϕ 𝜋 \phi/\pi as a function of the capillary number Ca Ca \rm Ca .

Ca	$\varphi /\pi$ Lac (2004)	$\varphi /\pi$ $\left(r=5,L=80\right)$	$\%$ Diff.
0.0375	0.21739	0.21929	0.19
0.075	0.18391	0.19694	1.3030
0.15	0.15826	0.16401	0.5750

Table 3. Table 3: Measured axial and transversal RBC diameters under a stretching test.

$F$(pN)	$D_a$($\mu$m) Mills (2004)	$D_a$($\mu$m)	$D_t$($\mu$m) Mills (2004)	$D_t$($\mu$m)
20	9.65672 $\pm$ 1.08955	9.19057	6.61194 $\pm$ 0.73134	6.31599
40	11.41791 $\pm$ 1.50746	9.95102	5.58209 $\pm$ 1.22388	6.08232
60	12.23881 $\pm$ 1.61194	10.98340	5.09851 $\pm$ 1.56716	5.81418
100	13.64179 $\pm$ 2.01309	12.19840	4.44179 $\pm$ 1.26866	5.09710
120	14.08955 $\pm$ 2.14925	12.69840	4.40299 $\pm$ 1.20896	4.70220
140	14.58209 $\pm$ 2.31343	13.77490	4.50746 $\pm$ 1.22388	4.52800
160	14.83582 $\pm$ 2.32836	14.58190	4.46269 $\pm$ 1.01256	4.30200
180	15.34328 $\pm$ 2.41791	15.2200	4.25373 $\pm$ 0.94030	4.14313

Table 4. Table 4: Table listing the parameters used in the simulations.

Parameters	Symbol	Value
Tube height, size	$H=W,D$	$25,15\mu$m
Tube segment lenght	$L$	$60\mu$m
RBC radius	$r_{\mathrm{rbc}}$	$4\mu$m
CTC radius	$r_{\mathrm{ctc}}$	$12\mu$m
RBC count	$N_{\mathrm{rbc}}$	Variable
Hematocrit	$H_t$	$N_{\mathrm{rbc}}{V}_{\mathrm{rbc}}/V_{\mathrm{cap}}$
Plasma kinetic viscosity	$\nu$	$1.2\times {10}^{-6}$ m2/s
RBC stiffness modulus	$k_{\mathrm{rbc}}$	$10\mu$ N/m
RBC bending modulus	$k_{\mathrm{ctc}}^b$	$2.4\times {10}^{-19}$ J
RBC volume constraint	$k_{\mathrm{ctc}}^v$	$10$
CTC stiffness modulus	$k_{\mathrm{ctc}}$	$5-200\mu$ N/m
Reynolds number	$Re$	$2.5\times {10}^{-}2$
Center velocity (no cells)	$u_0$	$\nu Re/H$
Shear rate	$\dot{\gamma }$	$4u_0/H$
Pressure gradient	$f=\Delta p/L$	$16\nu \rho u_0/H^2$
Lattice resolution	$\Delta x$	$0.5714\mu$m
Time step	$\Delta t$	$\Delta x^2/6\nu$
Capillary number	$\mathrm{Ca}$	$\nu \dot{\gamma }{r}_{\mathrm{rbc}}/k_{\mathrm{rbc}}$
Adhesive number	$Ad$	$\sigma /\nu \dot{\gamma }{r}_{\mathrm{rbc}}$
Receptor-ligand resting lenght	$\lambda$	$50$ nm

Equations

$$f_i(X+c_i\Delta t,t$$

$$-\frac{\Delta t}{\tau }[f_i(X,t)-f_i^{\mathrm{eq}}(X,t)]+\Delta t{F}_i,$$

$$f_i^{\mathrm{eq}}(X,t)={\omega }_i\rho \left[1+\frac{1}{c_s^2}\left(c_i\cdot u\right)+\frac{1}{2c_s^4}(c_i\cdot u{)}^2-\frac{1}{2c_s^2}{u}^2\right],$$

$$F_i={\omega }_i(1-\frac{1}{2\tau })\left(\frac{c_i-u}{c_s^2}+\frac{c_i\cdot u}{c_s^4}{c}_i\right)\cdot f_{\mathrm{ib}}.$$

$$c_i=\{\begin{array}{ll}(0,0,0),& i=0,\\ (\pm 1,0,0),(0,\pm 1,0),(0,0,\pm 1),& i=1-6,\\ (\pm 1,\pm 1,0),(\pm 1,0,\pm 1),(0,\pm 1,\pm 1),& i=7-12,\\ (\pm 1,\mp 1,0),(\pm 1,0,\mp 1),(0,\pm 1\mp 1),& i=13-18\end{array}$$

$$f_{\mathrm{ib}}\left(X,t\right)=\sum _i{F}_i(t)D(x_i(t)-X)A_i,$$

$${\dot{x}}_i(t)=\sum _{X}u(X,t)D(x_i(t)-X),$$

$${\phi }_2(x)=\{\begin{array}{ll}1-|x|,& \text{if}0\le |x|\le 1,\\ 0,& \text{if}|x|>1\end{array}$$

$${\phi }_4(x)=\{\begin{array}{ll}\frac{1}{4}\left(1+\cos (\frac{\pi x}{2})\right),& \text{if}|x|\le 2,\\ 0,& \text{if}|x|>2\end{array}$$

$$F_{n,m}^{\mathrm{s}}=-\mathrm{k}\sum _{n^{\prime },m^{\prime }}\frac{d_{n^{\prime },m^{\prime }}-d_{n,m}^0}{d_{n,m}^0}\frac{(x_{n^{\prime },m^{\prime }}-x_{n,m})}{d_{n,m}},$$

$$F_{n,m}^{\mathrm{b}}=-{\mathrm{k}}^{\mathrm{b}}\sum _j\left(\frac{{\theta }_j-{\theta }_j^0}{{\theta }_j^0}\right){\hat{n}}_j,$$

$$F_{n,m}^{\mathrm{v}}=-{\mathrm{k}}^{\mathrm{v}}\left(1-\frac{V}{V^0}\right)A_{n,m}{n}_{n,m},$$

$$F_{n,m}^{\mathrm{a}}=-{\mathrm{k}}^{\mathrm{a}}\left(1-\frac{A}{A_0}\right){\zeta }_{n,m},$$

$$U_M(r)=D_e(e^{2\beta (r_0-r)}-2e^{\beta (r_0-r)}),$$

$$F_{n,m}^{\mathrm{inte}}=\{\begin{array}{ll}-2\beta D_e{U}_M\frac{x_{n,m}-x_{n^{\prime },m^{\prime }}}{|x_{n,m}-x_{n^{\prime },m^{\prime }}|},& \text{if}{d}_{n,m}<d_{\mathrm{cut}}\\ 0,& \text{otherwise}.\end{array}$$

$$F_{n,m}^{\mathrm{adh}}=\{\begin{array}{ll}-\sigma (l_{n,m}-\lambda )\frac{x_{n,m}-x_{\mathrm{wall}}}{l_{n,m}},& \text{if}{d}_{n,m}^{\mathrm{w}}<d_{\mathrm{wall}}\\ 0,& \text{otherwise}.\end{array}$$

$$\{\begin{array}{l}x=a\alpha \sin (\theta )\cos (\varphi )\\ y=\frac{a}{2}(0.207+2.003{\sin }^2(\theta )-1.123{\sin }^4(\theta ))\cos (\theta )\\ z=a\alpha \sin (\theta )\sin (\varphi )\end{array}$$

Full text

*Methods – * A computational model is presented to predict the vascular transport and adhesion of CTCs in whole blood. A Lattice-Boltzmann method, which is employed to solve the Navier-Stokes equation for the plasma flow, is coupled with an Immersed Boundary Method.

*Results – * The vascular dynamics of a CTC is assessed in large and small microcapillaries. The CTC shear modulus ${\rm k}_{\rm ctc}$ is varied returning CTCs that are stiffer, softer and equally deformable as compared to RBCs. In large microcapillaries, soft CTCs behave similarly to RBCs and move away from the vessel walls; whereas rigid CTCs are pushed laterally by the fast moving RBCs and interact with the vessel walls. Three adhesion behaviors are observed – firm adhesion, rolling and crawling over the vessel walls – depending on the CTC stiffness. On the contrary, in small microcapillaries, rigid CTCs are pushed downstream by a compact train of RBCs and cannot establish any firm interaction with the vessel walls; whereas soft CTCs are squeezed between the vessel wall and the RBC train and rapidly establish firm adhesion.

*Concluisons – * These findings document the relevance of cell deformability in CTC vascular adhesion and provide insights on the mechanisms regulating metastasis formation in different vascular districts.

1 Introduction

The shedding into the vascular network of so-called ‘circulating tumor cells’ (CTCs) is the main mechanism by which malignant masses colonize distant organs and tissues [1, 2, 3]. After leaving the primary cancerous lesion, following a complex set of biological adaptations, CTCs face the blood stream and, as any other blood cell, are transported away along the vascular network. Although the mechanism by which CTCs select their final homing tissue is not yet fully understood, experimental evidence supports the notion that CTCs tend to more efficiently interact with the vessel walls and eventually extravasate in microcapillaries, with diameters in the range of a few tens of micrometers [4]. The vascular transport, adhesion and subsequent extravasation of CTCs are regulated by local hemodynamics and biological conditions and affected by the cell size and deformability. Individual CTCs present an average radius ranging from $5$ to over $15\mu$m [5]. More interestingly, multiple studies have documented significant differences in deformability among cancer cells [6, 7, 8]. In general, atomic force microscopy, optical and magnetic tweezer-based assays, and micropipette aspiration studies have demonstrated that malignant cells are more deformable than their healthy counterparts. For instance, Hrynkiewicz and his group tested multiple cells lines using scanning force microscopy and documented a difference of about $1$ order of magnitude between healthy and cancer cells[9]. Using primary cells from humans, oral carcinoma cells were found to be about $3$ times more compliant than healthy cells [10]. Similar observations were provided by the group of Gimzewski [11]. Recently, the group of Manalis has elegantly compared the deformability of cancer cells (lung, breast and prostate cancer) directly to that of blood cells (erythrocytes, leukocytes, and peripheral monocytes), using a suspended microchannel resonator and measuring cell passage times through a constriction [12]. The work concludes that CTC deformability can be larger, comparable or lower than that documented for blood cells.

Blood flow dynamics in microvessels is primarily ruled by a non-Newtonian effect called the Fåhræus-Lindqvist effect [13], which is characterized by the RBC migration away from the walls and progressive accumulation in the vessel core. At the continuum scale, the complex rheological behaviour of blood can be mathematically treated by adopting proper constitutive models, such as three-dimensional neo-Hookean or viscoelastic laws [14, 15]. However, within a microvascular network, cell-cell collisions and cell-cell adhesive interactions need to be explicitly accounted for. As such, mesoscale models are required in order to analyze the spatial and temporal evolution of blood flow and its cellular component [16]. Starting with the pioneering work of Pozrikidis and colleagues [17], who described blood cells using immersed boundary methods (IBM), different computational approaches have been presented in the open literature to capture the dynamics of a multitude of deformable cells. The group of Karniadakis has adopted a dissipative particle dynamics method (DPD) in which the RBC membrane is represented by a coarse-grained spring network [18, 19]. More recently, Gompper and colleagues introduced a particle-based mesoscopic simulation technique, called the smoothed dissipative particle dynamics (SDPD) method, which combines smoothed particle hydrodynamics and dissipative particle dynamics [20, 21].

IBM coupled on one side with a fluid solver and on the other side with a structural solver represents a potent tool for predicting the transport of cells and particles in microcapillaries [22]. The structural solver serves to capture the deformation of the immersed object over time under different hydrodynamic stresses. This could be a network of viscoelastic springs or a membrane with a complex rheological behavior solved by a Finite Element procedure [23, 24, 25]. The latter requires heavy parallelization as compared to the spring network limiting the maximum number of immersed objects.

In recent years the Lattice Boltzmann method (LBM) has been widely used as a fluid solver in IBM schemes [26, 27, 28]. In 2D problems, LBM-IBM approach with elastic spring networks has been extensively documented in simulation of transport of RBCs, cells and particles in microcapillaries [29, 30, 31]. In 3D simulations, community codes such as LAMMPAS (large-scale atomic/molecular massively parallel simulator) and ESPResSo (Extensible Simulation Package for RESearch on SOft matter) have been recentely proposed [32, 33].

In the literature, only a few works are focused on predicting the microvascular transport of circulating tumor cells. For instance, Rejniak used a $2D$ Immersed Boundary method to study the interaction of a single tumor cell with the vascular endothelium in pure plasma [34]. Yan and colleagues studied the vascular adhesion of a rigid, spherical cell in either a curved or a straight capillary, using a LBM method for describing blood flow, without accounting for the presence of RBCs[35]. More recently, the same group adopted a DPD computational scheme for predicting the transport of an individual tumor cell, initially released at the blood vessel wall, in the presence of RBCs [36]. This study demonstrates that RBCs enhance CTC adhesion in small capillaries whereas, in large vessels, CTC can be more easily detached away from the wall, especially at higher hematocrits. In the same work, a preliminary analysis of CTC deformability was also provided, considering a cell shear modulus about $25,2.5$ and $0.25$ time larger than that of RBCs. It was concluded that softer cells can engage a larger number of ligand-receptor bonds upon adhesion with the vascular walls.

In the present work, the effect of CTC deformability on metastasis formation is analyzed using a hierarchical computational model, where Lattice Boltzmann and Immersed Boundary methods are combined together. The Lattice-Boltzmann method is employed to solve the Navier-Stokes equation governing the pure plasma flow [37, 27, 28, 26, 38]; whereas, the Immersed Boundary Method [22] is adopted to describe the deformation and transport of RBCs and CTCs. The cell membranes are discretized as an ensemble of linear elastic springs, connecting neighboring membrane points. Following the seminal work of Hammer and colleagues [39, 40], an adhesive potential is also included to describe vascular adhesion, as mediated by the formation of individual receptor-ligand bonds treated as linear elastic springs [39, 41, 30]. The aim of this work is to elucidate the role of CTC deformability on vascular margination and subsequent adhesion in the presence of whole blood flow.

2 Computational method

The presented hierarchical computational model relies on the combination of a fluid solver for the incompressible Navier-Stokes equation, based on the three-dimensional D3Q19 Lattice-Boltzmann Method (LBM) [37, 27, 29], and a structural solver for the dynamics of the deformable membrane, based on an Immersed Boundary Method (IBM) [22]. Details of the physical model, governing equations and numerical implementation are given in the sequel.

3 Lattice Boltzmann method

The LBM introduces a number of $N$ populations $\{f_{i}\}$,$(i=0,\dots,N-1)$ streaming along a regular lattice in discrete time steps. These populations can be regarded as mesoscopic particles propagating and colliding. The evolution of the $N$ populations is given by the Lattice-Boltzmann equations [37], which takes the form:

[TABLE]

where $X$ is the spatial coordinate on a Cartesian regular lattice and $t$ is the time coordinate; $\{c_{i}\},(i=0,\dots,N-1)$ is the set of discrete velocities; $\Delta t$ is the time step; and $\tau$ is the relaxation time. The kinematic viscosity $\nu=c^{2}_{s}(\tau-\frac{1}{2})\Delta t$ of the flow is related to the relaxation time $\tau$, being $c_{s}=\frac{1}{\sqrt{3}}\frac{\Delta x}{\Delta t}$ the reticular speed of sound. The local equilibrium density functions $\{f^{\rm eq}_{i}\}$ are expressed by the Maxwell-Boltzmann distribution:

[TABLE]

where $\{\omega_{i}\},(i=0,\dots,N-1)$ are the lattice weights, depending on the underlying lattice structure; $\rho$ is the density and $u$ is the velocity field. A forcing term $f_{\rm ib}$, having the dimension of a body force density, can be incorporated via $F_{i}$ as:

[TABLE]

It is worth noting that the fluid interacts with an immersed object only via the forcing $f_{\rm ib}$. In this regard, the body force density $f_{\rm ib}$ is the term linking togheter the LB and IB modules. From equations (1-3), macroscopic quantities can be recovered respectively as the fluid density $\rho=\sum_{i}f_{i}$ and velocity $\rho u=\sum_{i}c_{i}f_{i}+\Delta tF_{i}/2$.

On the three-dimensional square lattice with $N=19$ speeds (D3Q19) [42], the set of discrete velocities is given by:

[TABLE]

with the weight, $\omega_{i}=1/18$ for $i=1,\dots,6$, $\omega_{i}=1/36$ for $i=7,\dots,18$, and $\omega_{0}=1/3$.

3.1 Immersed Boundary method

In the IB module, two independent meshes are considered to approximate respectively the immersed membranes and the fluid domain. The structure is discretized by a moving Lagrangian mesh in which the position of each node is $x_{i}(t)$, while the fluid is discretized by a fixed Eulerian mesh. The different steps in the computational algorithm are the following [22]. First, compute the total particle forces $F_{i}(t)$ acting on the Lagrangian point $x_{i}(t)$ of the immersed object (see Fig.1). These forces account for the internal elastic forces, interaction forces, adhesive forces and define the biological and mechanical behavior of the immersed object. Second, spread forces from the Lagrangian to the Eulerian mesh via the interpolation:

[TABLE]

where the index $i$ ranges over all Lagrangian points $x_{i}$ inside the interpolation stencil surrounding the Eulerian point $X$. $A_{i}$ is the area element associated with the Lagrangian node $x_{i}$. The operator $D$ is the discretized delta function. Third, solve the Lattice-Boltzmann equation for the fluid and find the velocity vector $u(X,t)$. Fourth, interpolate the fluid velocity to derive the velocity at each boundary node:

[TABLE]

where $X$ ranges over all Eulerian points inside the interpolation stencil surrounding $x_{i}$. Fifth, update the particle position as $x_{i}(t+1)=x_{i}(t)+\dot{x}_{i}(t)\Delta t$.

The choice of the discretized delta function is not unique and characterize the size of the interpolation stencil. Let $r=x_{i}-X$, the Dirac delta function $D(r)$ is factorized as $D(r)=\varphi_{d}(r_{x})\varphi_{d}(r_{y})\varphi_{d}(r_{z})$, where $\varphi_{d}$ define an interpolation kernel and $d$ is the support of the stencil. Two kernel functions have been used in the simulations:

[TABLE]

and

[TABLE]

The choice of $\varphi_{2}$ leads to an interpolation stencil consisting of $2^{3}$ Eulerian points for the interpolations of the forces and velocities [15, 43], while $\varphi_{4}$ leads to an interpolation stencil containing $4^{3}$ points.

3.2 Constitutive model for the membrane of the immersed object

A regular mesh is used to approximate the surface of the immersed object. This could be a spherical or biconcave capsule if it represents a CTC or a RBC, respectively. Meridional and azimuthal angles are used to describe the surface of the immersed object. The advantage of using a parametric description of the mesh is that the local connectivity of each node of the network, i.e. its neighborhood structure, is simply induced by the equations of the geometric surface. This is particularly useful when the elastic forces are computed. Each Lagrangian node on the surface of the immersed object $x_{n,m}$ is identyfied by two parameters $0\leq n,m\leq N$ , such that the total number of points on the surface is $N^{2}$. The stretching force $F^{\rm s}$ acting on each Lagrangian node $x_{n,m}$ is defined by summing up the forces on the springs connected to its neighboring nodes $x_{n^{\prime},m^{\prime}}\in\{x_{n\pm 1,m},x_{n,m\pm 1}\}$ given by [44, 45]:

[TABLE]

where $d_{n,m}=|x_{n^{\prime},m^{\prime}}-x_{n,m}|$ is the distance between node $x_{n,m}$ and its neighbor $x_{n^{\prime},m^{\prime}}$ and ${\rm k}$ is the shear elastic modulus. A simple bending term can be considered as:

[TABLE]

where index $j$ runs over adjacent edges $\{(n,m),(n\pm 1,m\pm 1)\}$, $\theta_{j}$ are the angles between adjacent triangles sourranding point $x_{n,m}$ and $\hat{n}_{j}$ is the outward unit vector related to edge $j$. To enforce membrane incompressibility, a constraint on the total volume is needed:

[TABLE]

where ${\rm k}^{\rm v}$ is the volume constraint factor, $V^{0}$ and $V$ are the volumes of the immersed object in the reference and current configuration respectively, $n_{n,m}$ is the outward unit normal vector associated to the Lagrangian point $x_{n,m}$ and $A_{n,m}$ is the area of the triangular element defined by the three mesh points $\{x_{n,m},x_{n+1,m},x_{n,m+1}\}$. To evaluate the normal unit vector and the area $A_{n,m}$ consider $\xi_{1}=x_{n+1,m}-x_{n,m}$ and $\xi_{2}=x_{n,m+1}-x_{n,m}$, then the unit outward normal vector is $n_{n,m}=\xi_{1}\wedge\xi_{2}/|\xi_{1}\wedge\xi_{2}|$ and the area of the triangular element is $A_{n,m}=|\xi_{1}\wedge\xi_{2}|/2$. The volume is computed at each time iteration via the discrete Green theorem as $V=\sum_{T}n_{T}\cdot c_{T}A_{T}/3$, where the summation index runs over all triangles $T=\{x_{n,m},x_{n+1,m},x_{n,m+1}\}$, $n_{T}$ the unit normal vector to the surface associated to triangle $T$, $A_{T}$ is the area of the triangular element, and $c_{T}=(x_{n,m}+x_{n+1,m}+x_{n,m+1})/3$ is the baricenter of the triangle $T$. Another constraint is imposed on the conservation of the total surface of the cell:

[TABLE]

where ${\rm k}^{\rm area}$ is the area constraint factor and $\zeta_{n,m}$ is the unit vector pointing from the centroid of triangle $T=\{x_{n,m},x_{n+1,m},x_{n,m+1}\}$ to the vertex $x_{n,m}$ . This constraint is applied in the case of RBCs and CTCs to account for their cytoskeleton [46, 47].

3.3 Particle-particle interactions

The cell-cell interactions can be approximated with the Morse potential given by [13]:

[TABLE]

where $r$ is the distance between two cells and $D_{e}$ is the energy well depth. Note that the calibration of the Morse potential can be done following the same procedure in [41] returning $D_{e}=0.04k_{B}T$, where $k_{B}T=4.14\times 10^{-21}N\cdot m$. Let $d_{n,m}=|x_{n,m}-x_{n^{\prime},m^{\prime}}|$ be the distance between two nodes on different structures. The interaction force between node $x_{n,m}$ on the surface of the current cell and node $x_{n^{\prime},m^{\prime}}$ on another cell is given by:

[TABLE]

The Morse interaction potential is implemented between two nodes of separate cells if they are within a cutoff distance $d_{\rm cut}$. This type of interaction consists of a high short-range repulsive force when $r<r_{0}$ and a low long-range attractive force for $r>r_{0}$. Parameters used are: $\beta=3.84\ \mu$m*-1*, $r_{0}=0.5\ \mu$m, and $d_{\rm cut}=1.5\ \mu$m [13, 48].

3.4 Wall-particle interactions

Ligand and receptor molecules are distributed over the cell and blood vessel surfaces, respectively. Ligand molecules are modeled as linear springs which tend to establish bonds with receptors on the vascular wall. Let $x_{n,m}$ be a point on the particle surface, $x_{\rm wall}$ be the normal projection of $x_{n,m}$ on the wall of the channel, and $l_{n,m}=|x_{n,m}-x_{\rm wall}|$ be the distance between a point on the cell surface and the corresponding point on the wall, when the distance between the surface of the capsule and the wall is less than a critical distance $d_{\rm wall}$, the adhesive force acting on $x_{n,m}$ is defined as [34, 36]:

[TABLE]

where $d^{\rm w}_{n,m}=d(x_{n,m},\rm wall)$ is the distance of $x_{n,m}$ from the wall, $\sigma$ is the adhesion constant, $\lambda$ is the equilibrium separation distance for the spring and $d_{\rm wall}=2\lambda$ is the critical distance. It is here important to recall that the mathematical model for cell-wall adhesion mediated by receptor-ligand interactions was originally presented by Hammer and colleagues and applied to study the rolling dynamics of leukocytes on endothelial surfaces [39]. Since then, this Adhesive-Dynamics model (AD) has been extended to study many biologically relevant problems, such as the hydrodynamic recruitment of rolling leukocytes [40], platelet-surface and platlet-platlet interactions [49, 50]. Note also that the current model does not account for the stochastic formation and rupture of ligand receptor bonds, which can be readily included following previous works by the authors and others [39, 40, 30].

4 Results and discussion

The proposed hierarchical computational model combines a Lattice-Boltzmann (LBM) algorithm, for solving the fluid flow, with an Immersed Boundary method (IBM), for determining particle-fluid and particle-wall interactions. As such, the computational model can efficiently deal with multiple scales and different biophysical problems, spanning from the ligand-receptor adhesive interactions (molecular scale) to the deformation of cell membranes (mesoscopic scale) and the transport of multiple red blood cells (RBCs) in a capillary flow (macroscopic scale). This is schematically presented in Fig.2. The computational model is first validated against known test cases: a deformable spherical capsule in a linear shear flow; the stretching of a red blood cell under uniaxial loading. Finally, the model is applied to document the vascular transport and adhesion dynamics of a single circulating tumor cells (CTCs), in whole blood flow, in microvessels of different calibers.

It is worth noting that in the present paper, only one CTC is considered. Given the low abundance of CTCs in blood, this condition is indeed physiologically sound. Other blood cells, such as leukocytes, platelets, monocytes and so on, are not explicitly modeled in this problem in that they are far less abundant than RBCs. Indeed, RBCs represent up to $95\%$ of all cellular components of blood and are responsible of the peculiar blood rheology. In small vessels, ranging from $10$ to $40\;\mu\rm m$, the volume fraction of RBCs varies from $15$ to about $35\%$ [18]. In the present work, the hematocrit has been fixed to the average value of $20\%$. Indeed, higher values for the hematocrit are associated with higher computational burden. The fluid solver used in this model relies on the Lattice-Boltzmann method which is very accurate and effective in terms of single cost per time iteration. However, a Lattice-Boltzmann solver is explicit in time and, as such, requires small time steps (in the order of $10^{-7}$ sec) to ensure stability. The effect of the mesh resolution on the accuracy of the solution is provided for the first test case.

4.1 A deformable spherical capsule in a linear shear flow: Test case 1

A spherical capsule of radius $r$ is placed in the center of a cubic box of size $L$ and subjected to a linear shear flow, realized by moving the top and bottom walls with velocities $u_{0}$ and $-u_{0}$, respectively (Fig.3a). The size of the fluid mesh is $L=16r$. $6400$ Lagrangian points are used to discretize the capsule membrane. The fluid is considered as quasi-stationary in the so-called Stokes regime with a Reynolds number $Re={10}^{-2}$ and a shear rate $\dot{\gamma}=\nu Re/r^{2}$, where $\nu=1/6$ is the lattice viscosity. Thus, the magnitude of the imposed velocity $u_{0}$ is given as $u_{0}=\dot{\gamma}L/2$. The capsule has a shear elastic modulus $k$. The capillary number ${\rm Ca}$ $({\rm Ca}=\dot{\gamma}\nu r/k)$, which represents the relative effect of the viscous drag over surface tension, provides a measure of the capsule deformability: ${\rm Ca}=0$ means a rigid capsule ($k\rightarrow+\infty$); ${\rm Ca}>0$ means a deformable capsule.

In a linear shear flow, the capsule rotates and deforms assuming eventually the shape of an ellipsoid (Fig.3b,c), with $r_{a}$ and $r_{b}$ being the two principal axes of the ellipsoid. Fig.3b gives the steady state, normalized velocity distribution $u/u_{0}$ within the box, for ${\rm Ca}=0.075$. Four recirculation areas are clearly developing around the elongated capsule. The deformation of the capsule is quantified via the Taylor parameter $D$ ($D=(r_{a}-r_{b})/(r_{a}+r_{b})$). Deformed configurations of the capsule at steady state are shown in Fig.3c for three representative values of the capillary number, namely ${\rm Ca}=0,0.075$ and $0.15$. Indeed, the larger is the capillary number, the larger is the capsule deformation $D$. This is also presented quantitatively in Fig.3d, where the Taylors number $D$ is plotted versus the non-dimensional time $\dot{\gamma}t$, for different values of the capillary number ${\rm Ca}(=0.0375,0.075,0.15$ and $0.3)$. $D$ increases with time until a steady state configuration is reached for $\dot{\gamma}t\geq 3$. For longer times, the capsule membrane rotates along its own shape with a inclination angle $\phi$, as predicted by the known tank-treading condition [51]. In Fig.3d, the present solution is compared with data from a neo-hookean membrane model, solved using the Boundary Integral method (BIM), by Lac and colleagues [51] for two values of the mesh resolution, namely low ($r=5,L=80$, dotted lines) and high ($r=10,L=160$, solid lines) resolution. The presented results are in good agreement with the BIM data, for different ${\rm Ca}$ values [51]. Notice that the agreement between the two solutions improves as the mesh resolution increases. This appears to be particularly relevant at large ${\rm Ca}$, which corresponds to more deformable capsule. The difference between the two numerical solutions is plotted as a function of ${\rm Ca}$ in Fig.3e, at steady state. This difference grows with the capillary number ${\rm Ca}$, as previously pointed out [45], but it stays well below $4\%$ for all considered cases. Table.1 summarizes the $D$ values for different ${\rm Ca}$ numbers obtained by BIM simulations and the current method, demonstrating that the difference ranges from $0.047\%$ for ${\rm Ca}=0.0375$ to $3.04\%$ for ${\rm Ca}=0.3$. In Table.2, are reported the values of the capsule angle $\phi/\pi$ for ${\rm Ca}=0.0375,0.075$ and $0.15$. The angle $\phi$ is computed as

$\phi=0.5\times\arctan\left(2r_{a}r_{b}/(r^{2}_{b}-r^{2}_{a})\right)$.

4.2 The stretching of a red blood cell under uniaxial loading: Test case 2

A single red blood cell (RBC) is stretched longitudinally by applying a force $\rm F$ at two opposite sites of the cell membrane (Fig.4a). After a transient phase, the elastic reaction force ${\rm F}^{\rm el}$ arising at the cell membrane balances out the external applied force $\rm F$ so that a steady deformation is achieved. This case serves to predict the stretching of a RBC in a pulling test realized using an optical tweezer [46, 47]. Briefly, two silica microbeads are attached at opposite sites of the cell membrane: one bead is anchored to the surface of a glass slide, while the other one is trapped by a laser beam. By moving the bead attached to the glass slide, a well defined strain is applied to the cell. At equilibrium, the diameter ${\rm D}_{a}$ in the pulling direction (axial) and the diameter ${\rm D}_{t}$ orthogonal to the pulling direction (transverse) are measured for each value of the applied force. This is documented by solid dots in Fig.4b. Experimentally the traction force $\rm F$ ranges from [math] to $180$ pN.

The $3D$ shape of the RBC is given by the following parametric equations, for $0\leq\theta\leq 2\pi,0\leq\phi\leq\pi$:

[TABLE]

where $a=0.74{\rm r}_{\rm rbc}$, with ${\rm r}_{\rm rbc}=3.91\ \mu$m being the equivalent RBC radius and $\alpha=1.39$. Note that this parametric equation is equivalent to the well known Evans-Fung formula in cartesian coordinates [46, 52, 53]. The RBC shear modulus is taken as ${\rm k}_{\rm rbc}=8.3\ \mu$N$/$m [47]. The solid lines in Fig.4b report the values of ${\rm D}_{a}$ and ${\rm D}_{t}$ computed at equilibrium for different values of $\rm F$ via the present hierarchical model. As expected, the diameter ${\rm D}_{a}$ increases while ${\rm D}_{t}$ decreases with the applied force $\rm F$. The computed axial and transverse diameters ${\rm D}_{a}$ and ${\rm D}_{t}$ are in good agreement with the experimental data [47]. Fig.4c shows representative equilibrium configurations of the RBC under different applied forces. The proposed hierarchical model captures correctly the mechanical deformation of RBCs with an overall error lower than $10\%$. Values of the diameters $D_{a}$ and $D_{t}$ for different applied forces $F$ are reported in Table 3 and compared with experimental data by [47]. The RBC bending modulus ${\rm k}^{b}_{\rm rbc}$ and volume constraint ${\rm k}^{v}_{\rm rbc}$ can be estimated as previously documented in [18, 47], returning ${\rm k}^{b}_{\rm rbc}=2.4\times 10^{-19}$ J ($=0.0016$ in LB units) and ${\rm k}^{v}_{\rm rbc}=10$ (LB units).

4.3 Margination dynamics of circulating tumor cells with different deformability

Cancer spreading to distant tissues (metastasis) involves the shedding in the circulation of malignant cells from a primary lesion; their vascular transport, adhesion, extravasation and proliferation [54, 2, 3, 4]. During their journey, CTCs reach peripheral vascular beds, with blood vessels characterized by a diameter ranging between $10$ and $40\ \mu$m. In these small vessels, the presence of RBCs would favor the margination of CTCs towards the vascular endothelium, just like for leukocytes in an inflamed vessel. While the shear modulus of a RBC is in the range of ${\rm k}_{\rm rbc}=5-10\ \mu$N $/$m [46], CTC deformability can vary from fractions of $1$ kPa (soft) up to $100$ kPa, which is $20$ times higher than a RBC (rigid), under physiological conditions [55, 12]. In this section, the effect of cell deformability on CTC vascular dynamics is analyzed.

The proposed hierarchical computational model is applied to study the transport of an ensemble of RBCs and a single CTC within a capillary for an hematocrit of $20\%$, which is a physiological value in the microcirculation. It is well known that, due to their deformability, RBCs move away from the walls and migrate to the center of the capillary [13, 19, 56]. This results in the formation of a region depleted of cells next to the wall, which is called the cell-free layer (CFL). This non-Newtonian effect known as the Fåhræus-Lindqvist effect [19] is responsible for the modulation of the blood viscosity.

Following [41, 32], a capillary with a square cross section of $H=W$ and lenght $L$ is considered. The equivalent radius of the RBC is ${\rm r}_{rbc}=3.9\ \mu$m, so that $V_{\rm rbc}=94.1\mu$m3. Following the same equation described in the previous section, RBCs are modelled as biconcave membranes.

A number $N_{\rm rbc}$ of vectors, representing the RBC and CTC centers of mass, are uniformely generated in space and positioned in the volume $V_{\rm cap}$. Each RBC has initially a random orientation with respect to the $x$- and $z$- axes. For each RBC, the number of Lagrangian points is $4900$. Periodic boundary conditions are imposed in the flow direction at the inlet and outlet sections of the capillary, while no-slip velocity boundary conditions are prescribed on the remaining walls. Bounce-back boundary conditions are employed to treat the no-slip velocity conditions at the walls. The Reynolds number is fixed to be $Re=2.5\times{10}^{-2}$. The blood flow is driven by a constant body force density $\rho f$, which is equivalent to prescribe a pressure difference over the lenght of the capillary given by $\Delta P/L=16u_{0}\nu\rho/D^{2}$, where $u_{0}=\nu Re/H$ is the peak velocity in the flow direction, $D=4HW/2(H+W)=H$ is the hydraulic diameter of the channel. The non-dimensional time $\dot{\gamma}t$ is considered, where $\dot{\gamma}=4u_{0}/H$ is the shear rate. The capillary number is defined as ${\rm Ca}=\nu\dot{\gamma}{\rm r}_{\rm rbc}/{\rm k}_{\rm rbc}$ and fixed to $0.026$ for the RBCs. The Lattice resolution is $\Delta x=0.5714\ \mu$m and the Lattice-Boltzmann viscosity is given by $\nu=\frac{1}{6}\Delta x^{2}{\Delta t}^{-1}$. Note that, in dimensional units, the viscosity is equal to $\nu=1.2\times{10}^{-6}$ m2s*-1* and the plasma density is $\rho=1000$ Kg$/$m3. A capillary with a square cross section of $H=W=25\ \mu$m and lenght $L=60\ \mu$m is considered. The CTC is modelled as a spherical capsule having a radius ${\rm r}_{\rm ctc}=12\ \mu$m and discretized with $6400$ Lagrangian points. A Poiseuille flow with Reynolds number $Re=0.25\times{10}^{-2}$ is assumed. Periodic boundary conditions are prescribed at the front and back walls of the capillary, while no-slip velocity boundary conditions are imposed on the remaining walls. Table4 collects the values of all physical parameters used in this simulation. The initial position of the CTC center of mass is $(x_{0},y_{0},z_{0})=(W/2,L-2{\rm r}_{\rm ctc},H/2-{\rm r}_{\rm ctc})$.

Three different CTC deformability values are considered: a CTC softer than RBCs (${\rm k}_{\rm ctc}=0.5{\rm k}_{\rm rbc}$: SOFT); a CTC stiffer than RBCs (${\rm k}_{\rm ctc}=10{\rm k}_{\rm rbc}$: RIGID); a CTC as deformable as a RBC (${\rm k}_{\rm ctc}={\rm k}_{\rm rbc}$: EQUAL). Fig.5a,c show the RBC distribution and the CTC location within a longitudinal section of the capillary, at time $\dot{\gamma}t=100$, for the ‘soft’, ‘rigid’ and ‘equal’ cases. On the right, fluid velocity profiles are shown. Specifically, the different velocity profiles are referred to the classical Poiseuille case in the absence of RBCs (black line); the time-averaged velocity profile for ${\rm H_{t}}=20\%$ with RBCs and no CTC (red line); and the time-averaged velocity profile for ${\rm H_{t}}=20\%$ with RBCs and the CTC (blue line). The presence of the CTC does not change significantly the time-averaged velocity profile, whereas the addition of RBCs flattens the classical Poiseuille parabolic profile, as well documented in the literature. In the ‘rigid’ case (${\rm k}_{\rm ctc}=10{\rm k}_{\rm rbc}$), malignant cells are rapidly pushed out from the center of the capillary and confined to move within the CFL (Fig.5b). This indeed increases the likelihood of building adhesive interactions with the wall. For the other two cases (${\rm k}_{\rm ctc}=0.5{\rm k}_{\rm rbc}$ and ${\rm k}_{\rm ctc}={\rm k}_{\rm rbc}$), malignant cells are not observed to marginate within the considered simulation time. This is due to the fact that RBCs and CTCs would move similarly in the channel, deforming under flow and moving away from the walls (Fig.5b,c).

In Fig.5d, the CTC trajectories are presented for the three considered cases. In the ‘soft’ and ‘equal’ cases (red and blue lines), malignant cells migrate towards the centerline $z=H/2$ and stay within the capillary core without interacting with the vessel walls throughout the simulation time. Conversely, in the ‘rigid’ case, malignant cells deviate from the streamlines and, eventually, reach the capillary wall (margination). For these simulations, the cell-wall adhesion potential is turned off and a moderate repulsive force is included only to prevent body compenetration.

4.4 Adhesion dynamics of circulating tumor cells in a large microcapillary: the $25\mu\rm m$ case

After marginating towards the vessel walls, CTCs could firmly adhere to the endothelial cells, if proper conditions are met. The margination and vascular adhesion are fundamental steps in the cascade of events regulating the extravasation of both leukocytes, in inflammation, and CTCs, in metastasis. Adhesive interactions are governed by receptor molecules, expressed on the vascular endothelial cells and ligand molecules, distributed over the CTC membrane (Fig.2). These molecular interactions operate only, and only if, CTCs are sufficiently close to the wall, namely closer than a critical distance set to $d_{\rm wall}=100\ \rm nm$. The molecular bonds are computationally treated as linear springs, whose strength is dictated by the adhesive number $\rm Ad=\sigma/\nu\dot{\gamma}{\rm r}_{\rm ctc}$. Simulations are performed by fixing a value of the adhesive strenght $\sigma$ and varying the CTC stiffness ${\rm k}_{\rm ctc}$. In this section, the proposed hierarchical computational model is applied to predict the vascular adhesion dynamics of CTCs as a function of their deformability.

Based on the results of the previous paragraph, the CTC stiffness is here assumed to be high enough to allow rapid margination, namely ${\rm k}_{\rm ctc}=10,15,20$ and $25{\rm k}_{\rm rbc}$. For soft CTCs, margination would not occur in capillaries with a tube size of $D=25\ \mu$m, within the time of the simulations. Indeed, this implies that soft CTCs could interact with the vascular walls only in capillaries slightly larger, or even smaller, than the CTC tube size, as shown in the sequel. The relative position of CTCs, the size of the adhesion area and the adhesion forces are monitored over time. The CTC dynamics is presented in Fig.6a,b in terms of the vertical coordinate $z_{\rm ctc}(y)$ and corresponding velocity component $v^{z}_{\rm ctc}(y)$ as a function of the position along the $y$-axis, and of the horizontal coordinate $y_{\rm ctc}(t)$ and corresponding velocity component $v^{y}_{\rm ctc}(t)$ as a function of time (Fig.6c,d). The $z$-direction is normal to the flow, whereas the $y$-direction is aligned with the flow.

In Fig.6a the black line corresponds to ‘soft’ CTCs (${\rm k}_{\rm ctc}={\rm k}_{\rm rbc}$) for which margination does not occur within the considered simulation time. In this case, the vertical position $z_{\rm ctc}$ of the CTC oscillates sligthly around the centerline of the capillary ($z=H/2$) while the cell is transported downstream. This is also confirmed by the variation over time of $y_{\rm ctc}$, as shown in Fig.6c. The black line grows steadily over time implying that the CTC is steadily moving along the flow direction. As such, the velocity components $v^{z}_{\rm ctc}$ along $z$, in Fig.6b and $v^{y}_{\rm ctc}$ along $y$ in Fig.6d are, respectively, nearly zero and quasi-constant but larger than zero. Correspondingly, the size of the adhesion area $A_{\rm ctc}$ and adhesion force $F^{\rm adh}_{z}$ are both null (Fig.6e and f).

A totally different behavior is documented for ‘rigid’ CTCs (${\rm k}_{\rm ctc}=25{\rm k}_{\rm rbc}$), as per the red lines in Figs.6. The cell is pushed downstream by the flow but also laterally towards the vessel walls until a stable adhesion is established (point A). This is documented in Figs.6a and c by the reduction of $z_{\rm ctc}$ till zero (interaction with the lower vessel wall) and final constant value of $y_{\rm ctc}$. Also the velocity goes to zero, after a significant spike in $v^{z}_{\rm ctc}$ due to the margination process (point A in Fig.6b). The area of adhesion $A_{\rm ctc}$ grows upon interaction with the vessel wall and stays quasi-constant over time (Fig.6e). Similar observations apply for the adhesion force $F^{\rm adh}$ (Fig.6f). This is consistent with the stable adhesion of a relatively rigid cell that does not squeeze down onto the wall. The case ${\rm k}_{\rm ctc}=10{\rm k}_{\rm rbc}$ is depicted in Figs.6 by blue lines. The behavior is quite similar to that of ${\rm k}_{\rm ctc}=25{\rm k}_{\rm rbc}$, whereby the cell moves downstream and laterally towards the vessel wall and starts interacting with its surface. However, the margination process occurs over a longer time (see Fig.6a and c), with a smoother variation in the velocities (see Figs.6b and d). Interestingly, and differently from the more rigid case, the CTC preserves a non-zero $v^{y}_{\rm ctc}$ velocity, which implies that firm adhesion is not established but the cell is rather rolling steadly over the vessel wall (point B). Similarly, the size of the adhesion area and the value of the adhesion force do not change significantly over time after an initial increase (see Figs.6e and f). The stable rolling is supported by the continuous rupture and formation of ligand-receptor bonds, respectively, of the tail and leading edge of the adhesive area. Finally, the case ${\rm k}_{\rm ctc}=15{\rm k}_{\rm rbc}$ is depicted by the green line, in Figs.6. Under this condition, the CTC exhibits an even slower approach to the vessel wall. Also adhesion does not appear to be firm and complete suggesting a ‘crawling’ behavior over the vessel wall (point C). The velocity $v^{y}_{\rm ctc}$ is close to zero but not null (Fig.6d), the adhesion area grows steadily with time just like for the adhesion force documenting a progressive CTC flattening over the wall (see Figs.6e and f). Indeed, the higher CTC deformability favors its continuous deformation and conformation to the vessel walls under hemodynamic forces.

Figs.7 shows the RBC distribution and velocity field around the CTC under firm adhesion (Fig.7a) and steady rolling (Fig.7b) at different time points. Fig.7a shows images of the CTC undergoing firm adhesion on the vessel wall, after margination at different time points (${\rm k}_{\rm ctc}=25{\rm k}_{\rm rbc}$). When the CTC reaches the substrate, the adhesion force dominates over the lift hydrodynamic force. In this case the cell firmly adheres to the wall and experiences small variations in configurations due to the complex whole blood flow dynamics. The streamlines of the velocity $u/u_{0}$ show the initiation of a local recirculation area around the adhering CTC. The adhered CTC becomes an obstacle for the red blood cells, which are continuously hitting over the trailing edge of the CTC, and might detach the cell if adhesion is insufficiently strong. In Fig.7b, a CTC undergoing stable rolling dynamics on the wall is shown (${\rm k}_{\rm ctc}=10{\rm k}_{\rm rbc}$). The CTC moves along the vessel wall with a constant size of the adhesion area at the interface between the particle and the substrate. The rolling behaviour is confirmed by following the red spot on the cell membrane over time. It can be appreciated the increase of the contact area over time which is associated with the progressive flattening of the CTC over the vessel wall. Although a simplified adhesive model is used based on elastic springs, it is capable to predict different adhesive regimes. In particular, it is observed that, if sufficiently soft, CTCs can detach from the wall of a capillary.

4.5 Adhesion dynamics of circulating tumor cells in a small microcapillary: the $15\mu\rm m$ case

The analysis of the margination and adhesion process is here conducted in a capillary having a tube size comparable with the size of the CTC. The relative position of CTCs, the size of the adhesion area and the adhesion forces are monitored over time. The CTC dynamics is presented in Fig.8a,b in terms of the vertical coordinate $z_{\rm ctc}(y)$ and corresponding velocity component $v^{z}_{\rm ctc}(y)$ as a function of the position along the $y$-axis, and of the horizontal coordinate $y_{\rm ctc}(t)$ and corresponding velocity component $v^{y}_{\rm ctc}(t)$ as a function of time. The red and blue lines in Fig.8a,b correspond to ‘rigid’ CTCs, with ${\rm k}_{\rm ctc}=25{\rm k}_{\rm rbc}$ and ${\rm k}_{\rm ctc}=20{\rm k}_{\rm rbc}$. For these values, the vertical position $z_{\rm ctc}$ remains close to the centerline of the capillary ($z=H/2$). Correspondingly, $y_{\rm ctc}$ (red and blue lines in Fig.8c) is linearly increasing in time. As such, the velocity components $v^{z}_{\rm ctc}$ along $z$, in Fig.8b and $v^{y}_{\rm ctc}$ along y in Fig.8d are, respectively, nearly zero and quasi-constant but larger than zero. Also, for ‘rigid’ CTCs (${\rm k}_{\rm ctc}=25{\rm k}_{\rm rbc},{\rm k}_{\rm ctc}=20{\rm k}_{\rm rbc}$), the size of the adhesion area $A_{\rm ctc}$ and adhesion force $F^{\rm adh}$ are both null (Fig.8e and f). In the ‘rigid’ case, CTCs are transported along the flow direction by the fast moving RBCs and do not appear to be interacting with the vessel walls. Very different is the behavior observed for a ‘soft’ CTCs. This is shown in Fig.8a (green line: ${\rm k}_{\rm ctc}=10{\rm k}_{\rm rbc}$ ; violet line: ${\rm k}_{\rm ctc}=5{\rm k}_{\rm rbc}$ ; and black line: ${\rm k}_{\rm ctc}={\rm k}_{\rm rbc}$). $z_{\rm ctc}$ increases till the top capillary wall is reached. Spikes in the velocity component $v^{z}_{\rm ctc}$ (green, violet and black line in Fig.8b) demonstrate the interaction with the top wall and the subsequent establishment of adhesive interactions. In Fig.8c the components in the flow direction $y_{\rm ctc}$ for ‘soft’ CTCs, slightly deviate from the straight line and the velocity component $v^{y}_{\rm ctc}$ decreases, as in Fig.8d (green, violet and black line). This implies that the cell is no longer in the fluid phase and has established an adhesive interaction with the vessel walls. In Fig.8e,f, the variations of the contact area and adhesion force are shown for ‘soft’ CTCs (green, violet and black line). In all the cases, the contact area and the adhesion force increases with time. This is due to the interaction with RBCs and hydrodynamic lift force, which affects the formation of the contact area.

Fig9.a shows images of a CTC undergoing ‘train’ dynamics in a small capillary ($D=15\ \mu$m) , at different time points for ${\rm k}_{\rm ctc}=25{\rm k}_{\rm rbc}$. The ‘rigid’ CTC is transported downstream without marginating. This results in a stable dynamics in which the CTC becomes a moving obstacle inside the capillary. The RBCs pile up behind the cell and form a dense aggregate [57]. Fig9.b shows representative images of the dynamics of a ‘soft’ CTC (${\rm k}_{\rm ctc}={\rm k}_{\rm rbc}$) in a small capillary. Due to its deformability the cell interacts with the vessel walls establishing adhesion. This results in a partial occlusion of the capillary. RBCs tend to pile up behind the cell, continuously pushing the cell against the wall, until they are free to pass again.

It can be concluded that in small capillaries of size comparable with that of the CTC, the stiffness of the cell is responsible for a transition from a ‘train’ dynamics, in which the CTC moves along the capillary without interacting with the walls (‘rigid’ case), to margination and subsequent adhesion dynamics (‘soft’ case).

5 Conclusions

Using a combined Lattice Boltzmann-Immersed Boundary method, the transport of deformable CTCs in a whole blood capillary flow has been analyzed in terms of cell displacements and velocities, and cell interactions with the vessel walls. The evolution over time of the area of adhesion and adhesion forces exchanged at the cell-wall interface has been documented for CTC stiffer, softer and equally deformable as compared to RBCs. It has been demonstrated that the interaction between deformable CTCs and RBCs is crucial in shaping the metastatic process.

Rigid CTCs have been observed to marginate rapidly within $25\ \mu$m microcapillaries and efficiently interact with the vessel walls as they are pushed laterally by the RBCs. On the contrary, in smaller $15\ \mu$m microcapillaries, rigid CTCs cannot establish firm interactions with the vessel walls. This should be ascribed to the significant flow obstruction induced by a rigid CTC adhering in such a small vessel. In other words, the fast moving RBCs in $15\ \mu$m microcapillaries form a compact train that constantly pushes and dislodge downstream any obstacle, such as a rigid CTC. Different adhesive regimes have been predicted for the rigid CTCs depending on their relative stiffness to RBCs. Very rigid CTC would firmly adhere, if proper local hemodynamic and biophysical conditions ar met. Intermediate rigid CTCs would roll over the vascular walls, whereas CTCs that are slightly stiffer than RBCs could crawl over the surface as a combination of rolling and progressive squeezing against the wall. Soft CTCs have been observed to deform and navigate together with the RBCs in the core of the blood vessel in $25\ \mu$m microcapillaries. Very differently, in smaller $15\ \mu$m microcapillaries, soft CTCs can deform and squeeze progressively within the train of fast moving RBCs and the vessel walls. This indeed increases the surface area of the CTC exposed to the vessel wall and inevitably favors firm adhesion. As a consequence, soft CTCs are expected to have a higher longevity in blood and, possibly, the ability to evade more efficiently than rigid CTC the recognition by cells of the immune system. These findings highlight the role of CTC deformability in defining the metastatic potential of cancer cells.

6 Acknowledgements

This project was partially supported by the European Research Council, under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement no. 616695, by the Italian Association for Cancer Research (AIRC) under the individual investigator grant no. 17664, and by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 754490.

7 Conflicts of interest

Dr. Lenarda, Dr. Coclite and Dr. Decuzzi declare that they have no conflicts of interest.

8 Ethical standards

No animal studies or experiments with human samples were carried out by the authors for this article.