# Deciphering circulating tumor cells binding in a microfluidic system thanks to a parameterized mathematical model

Giorgia Ciavolella, Julien Granet, Jacky Goetz, Nael Osmani, Christèle Etchegaray, Annabelle Collin

<figure><img src="image_1.png" /></figure>

## To cite this version:

Giorgia Ciavolella, Julien Granet, Jacky Goetz, Nael Osmani, Christèle Etchegaray, et al.. Deciphering circulating tumor cells binding in a microfluidic system thanks to a parameterized mathematical model. Journal of Theoretical Biology, 2025, 600, pp.112029. 10.1016/j.jtbi.2024.112029. hal-04253803v3

HAL Id: hal-04253803

https://inria.hal.science/hal-04253803v3

Submitted on 8 Jan 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Deciphering circulating tumor cells binding in a
microfluidic system thanks to a parameterized
mathematical model

Giorgia Ciavolella<fnref n="1" /> Julien Granet<fnref n="1" /> Jacky G. Goetz<fnref n="2" /> Naël Osmani<fnref n="2" />
Christèle Etchegaray<fnref n="1" /> Annabelle Collin<fnref n="3" />

October 23, 2024

Abstract

The spread of metastases is a crucial process in which some questions remain unanswered. In this work, we focus on tumor cells circulating in the bloodstream, the so-called Circulating Tumor Cells (CTCs). Our aim is to characterize their trajectories under the influence of hemodynamic and adhesion forces. We focus on already available *in vitro* measurements performed with a microfluidic device corresponding to the trajectories of CTCs – without or with different protein depletions – interacting with an endothelial layer. A key difficulty is the weak knowledge of the fluid velocity that has to be reconstructed. Our strategy combines a differential equation model – a Poiseuille model for the fluid velocity and an ODE system for the cell adhesion model – and a robust and well-designed calibration procedure. The parameterized model quantifies the strong influence of fluid velocity on adhesion and confirms the expected role of several proteins in the deceleration of CTCs. Finally, it enables the generation of synthetic cells, even for unobserved experimental conditions, opening the way to a digital twin for flowing cells with adhesion.

2010 Mathematics Subject Classification. 62-07; 65L09; 76Z99; 97M60;
Keywords and phrases. Differential equations; Parameter estimation; Circulating tumor cells;
Biological data

# 1 Introduction

One of the most important and deadly features of solid tumors is the increased ability of cancer cells to migrate and invade other organs, which is called metastatic spread. In the last 70 years the number of cancer deaths registered with metastasis has tripled [1]. Although the incidence of individual types of cancer varies greatly, metastasis remains the main cause of cancer-related deaths [2].

<footnote n="1">¹Inria, Univ. Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251, F-33400 Talence, France</footnote>

<footnote n="2">²INSERM UMR_S 1109, Univ. Strasbourg, FMTS, Équipe labellisée Ligue Contre le Cancer, F-67000 Strasbourg, France.</footnote>

<footnote n="3">³Laboratoire de Mathématiques Jean Leray, Nantes Université, F-44100 Nantes, France</footnote>
The blood and lymphatic circulations are used as a means of transport to reach distant organs. Tumor cells that have previously detached from a primary tumor can invade the surrounding extracellular matrix. Successful intravasation into the vessels means that cancer cells can now leave the original site. Inside the blood vessels, hostile conditions prevail. Circulating Tumor Cells (CTCs) are subjected to physical stresses that include hydrodynamic flow and loss of attachment to a substrate, as well as other obstacles involving the human immune system (and platelets) [3, 4]. These factors lead to a significant decrease in the number of CTCs and also to their eventual clustering. The remaining single cells or small cell clusters eventually extravasate, reaching a secondary site where they either stay dormant or form a new tumor [5].

CTCs receive much research interest due to their therapeutic potential in liquid biopsy [6, 7]. Indeed, they could allow to monitor tumor heterogeneity or response to a treatment, but also to detect the minimal residual disease, and serve as a prognosis biomarker or as a target for personalized therapies [8–10]. However, the detection, identification and characterization of CTCs present important challenges due to their heterogeneity and low abundance [11]. From the biological standpoint, understanding the key steps involved in CTCs arrest on the endothelial wall is crucial to explain secondary tumour locations. Indeed, the possibility of extravasating is permitted by CTCs arrest and firm adhesion to the vascular endothelium, phenomena that need further insights [12–16].

Studies previously pursued by biologists Follain et al. [13] and Osmani et al. [14] have deepened into the mechanical cues that promote CTCs successful arrest and extravasation using statistical approaches. In [13], they have shown that an optimal flow is required for CTCs to arrest on the endothelium of the vascular wall. Furthermore, in [14], they have identified the adhesion receptors at play. Early adhesion is mediated by the glycoprotein CD44, involved in a weak form of bonds, while integrin ITGB1 favours stabilization of the adhesions. These findings originate from both *in vitro* and *in vivo* experiments on the metastatic D2A1 cell line, and were validated in a fully syngeneic experimental mouse metastasis model. The *in vitro* experiments consist in using a microfluidic channel with controlled fluid velocity (simulating a blood vessel) into which tumor cells are injected. They allowed to establish that both CD44 and ITGB1 favour cell arrest, while only ITGB1 is necessary for stable cell adhesion. The resulting data were validated *in vivo* using intravital imaging in zebrafish embryos where CTCs were pumped by the heart along the vascular architecture. In that setting, it was further proved that CD44 mediates early cell arrest (5 minutes after cell injection), while ITGB1 is mainly required for stable cell adhesion (3 hours after injection). These two molecules were identified *in vitro* among several others such as ITGA3 and ITGA4 that were not found significant, and ITGA5 that was identified as the partner subunit of ITGB1. Furthermore, ITGB3 was found to behave as CD44 in mediating cell arrest. In summary, these recent works point out the interplay of the flow and the CD44/ITGB3 and ITGB1/ITGA5 molecules on CTCs transient or stable adhesion to the endothelial wall. However, they do not inform on how their whole dynamics may be affected even before arrest. In the present work, we aim to exploit the *in vitro* data from [14] to derive cell trajectories and to interpret them in the light of **deterministic mathematical modeling**. This will allow to decipher the respective **roles of the shear flow** and of the identified molecules of interest in the early **interaction between CTCs and the endothelial wall**.

Various theoretical models of cell adhesion have been developed over time. First studies have focused on the binding dynamics of either single bonds or clusters sharing a constant or varying load [17–21]. In the case of inflow cell dynamics, several biological questions can be addressed, such as the emergence of several cell displacement regimes (freely-flowing, rolling, slipping, stationary arrest) with possible bistability or shear-threshold effect between them. A related issue concerns the bonds response to hydrodynamic forces, with *catch bonds* whose lifetime increases with load, *slip bonds* for which it decreases exponentially with load (so-called Bell's
law), or a combination of both depending on the shear rate.

First computational approaches allowed to describe a hard sphere submitted to hydrodynamic forces and stochastic binding interaction with the wall [22–28]. In [29, 30], the adhesion of a rolling sphere is described following the membrane approaching the wall at the front, and detaching at the rear, for catch-slip bonds. In [31], both translational and rotational motions of a spherical cell are affected by elastic bonds. This allows to explain the interplay between rolling and slipping, and provides a numerical state diagram of leukocyte motion. Theoretical models can also take ligands positions into account, thus enabling bonds tilting and subsequent cell sliding [32–36]. In the spirit of [36] but in the absence of space structure, stochastic and deterministic models are developed for a particle cell in [37]. Although minimal in the hydrodynamic description, these models have less parameters and are therefore more suited to calibration with experimental data.

Theoretical frameworks have been confronted to microfluidics experiments on CTCs in [38–40]. In particular, in [40], the authors perform microfluidics experiments to study the effect of the shear rate on the dynamics of breast cancer cells interacting with an EpCAM-coated wall. Three regimes were observed (freely-flowing, firmly adhering, and rolling/slipping). Experimental data consisted in trajectories and in stopping times and lengths that were used to empirically calibrate a model based on [31]. More precisely, the cell-wall gap, the typical adhesion force and the spring constant were sequentially identified by numerical investigations. Then, the cell velocity during capture was well fitted by a decreasing exponential function, yielding a typical decreasing time characteristic of the cell-wall interaction.

In this work, we aim to capture the role of hydrodynamic forces and of different adhesion proteins during the first phase of CTCs interaction with the endothelial wall using a mathematical model. First, we exploit the movies from the *in vitro* experiments from [14] and we derive cell trajectories by tracking. Next, we perform a statistical analysis of the resulting mean velocities to highlight significantly different mean dynamics. Then, we build a mathematical model to explain these differences, relying on the full reconstruction of the fluid velocity. We consider a Poiseuille model for the fluid velocity, and weakly couple it to a modification of the adhesion model proposed in [36, 37]. The cell velocity depends on both the fluid velocity and the bonds density, while the binding dynamics takes into account bonds formation, adhesion growth, and unbinding. A well-designed and robust calibration based on nonlinear mixed-effects models [41] is used to fit the model to the *in vitro* experiments carried out by Osmani and collaborators, see [14, 42].

The work is arranged as follows. Section 2 contains the main information about the biological data: see Subsection 2.1 for protocol details and Subsection 2.2 for data presentation. Trajectories and velocities of 278 cells with or without proteins depletions were extracted from these data. Section 3 is devoted to methods. After a brief statistical analysis of the data in Subsection 3.1, which shows the statistically significant slowing behavior of CTC velocities over time in most cases, Subsections 3.2 and 3.3 present the mathematical modeling. In Subsection 3.4, the parameter estimation strategy is presented. The strategy to obtain a digital twin after calibration is presented in Subsection 3.5. The results showing the good agreement of the model with the data and the obtained synthetic cells are presented in Section 4. A discussion is presented in Section 5, and finally, conclusions are given in Section 6.

## 2 Data

In this section, we present the experimental data, beginning with their acquisition and ending with their extraction. First, in Subsection 2.1 we present the experimental protocol. Then, in Subsection 2.2, we show what kind of data are obtained, and briefly present the tracking
<figure><img src="image_3.png" /><figcaption>Figure 1: Experimental setup. A rectangular microfluidic channel (black box) contains an endothelial layer (turquoise cells) and a fluid containing **DA21** tumor cells (grey cells). The fluid is pumped by a peristaltic pump (left in the figure) with a pressure gradient $G(1 + \xi_f \cos(\omega_f t + \varphi))$ at the entrance. A camera is placed to record the motion of cells located in a focal plane at a distance $h_f^m$ from the endothelial layer (grey cells surrounded by a red contour), see [42].</figcaption></figure>

techniques used to extract the trajectory and velocity of 278 cells. We then present the resulting cell velocities.

## 2.1 Protocol

This study relies on the experimental *in vitro* data obtained in [14], to which we refer to for further details. Human Umbilical Vein Endothelial cells (HUVEC, Promocell) were seeded at 30 000 cells per channel in a rectangular microfluidic channel (IBIDI) of length $L = 1.7 \times 10^{-2}$ m, of width $l = 3.8 \times 10^{-3}$ m and height $h = 4.0 \times 10^{-4}$ m. Medium was changed twice a day until they reach maximal confluency (3 to 4 days). **DA21** mouse breast carcinoma cells were with siRNA using Lipofectamine RNAiMAX (Thermo Fisher) following the manufacturer's instructions. Experiments were performed between 72 hours and 96 hours post-transfection. 3 days after siRNA transfection, **DA21** cells were resuspended at a concentration of $10^6$ cells/ml in a Hépès-buffered cell culture medium and perfused into the channel using a REGLO Digital MS-2/12 peristaltic pump (Ismatec), Tygon LMT-55 3-stop tubing (IDEX), 0.5 and 1.6 mm silicon tubing and elbow Luer connectors (IBIDI).

In the setup of the pump, the mean value of the entry pressure gradient is fixed and denoted by $G$. The fluid velocity generated – which contains oscillations due to the pump – depends on the position in the channel and is not measured. A cMOS camera (IDS) is placed to record the motion of cells located in a focal plane at a distance $h_f^m$ from the endothelial layer. The experimental data consist of timelapse movie acquired at a rate of 24 frames per seconds for 2 minutes, on a rectangle of width $\ell_{cam} = 5.63 \times 10^{-4}$ m and of height $h_{cam} = 2.99 \times 10^{-4}$ m. The setup is shown in Figure 1.
## 2.2 Data availability

An example of a video image is shown in Figure 2. The cells forming the endothelial layer are seen in the background, whereas moving CTCs that are not in the focal plane of the camera are seen in the foreground. These CTCs have a well-defined shape, so that their trajectories can easily be followed while they appear in the video. Most cells are smoothly transported through the fluid. Sometimes, a cell stops on the endothelial layer. This arrest can be stable, meaning that the cell remains attached to the endothelial layer, or unstable, i.e. the cell can detach due to a collision with other cells or under the effect of the hydrodynamic forces. However, in this work we focus only on non-arrested cells.

Three experiments are considered, each of them marked by a control case representing the unmodified biological conditions (D2A1 cells treated with a CTL siRNA):

(1) Experiment 1: Case ITGB1, CD44

* siCTL (ITGB1, CD44): control group;

* siITGB1: depletion of ITGB1 (D2A1 cells treated with a siRNA targeting ITGB1);

* siCD44: depletion of CD44 (D2A1 cells treated with a siRNA targeting CD44);

(2) Experiment 2: Case ITGA5

* siCTL (ITGA5): control group;

* siITGA5: depletion of ITGA5 (D2A1 cells treated with a siRNA targeting ITGA5);

(3) Experiment 3: Case ITGB3

* siCTL (ITGB3): control group;

* siITGB3: depletion of ITGB3 (D2A1 cells treated with a siRNA targeting ITGB3)

The different experimental subgroups are summarized in Table 1. Experiments have been realized keeping the fluid at a controlled pressure gradient by the peristaltic pump. Three possible values of pressure gradient cohorts have been considered: $G^{(1)} = 50.33 \text{ Pa.m}^{-1}$, $G^{(2)} = 100.66 \text{ Pa.m}^{-1}$ and $G^{(3)} = 201.32 \text{ Pa.m}^{-1}$. For each of these cohorts, control and depletion of adhesion proteins gives rise to seven sub-cohorts, see Western-blot results in Supplementary materials in [14].

We collect cell trajectories from the 17 different videos (Experiment 1: 3 different pressure gradients cohorts × 3 different proteins expressions sub-cohorts; Experiment 2: 2 cohorts × 2 sub-cohorts; Experiment 3: 2 cohorts × 2 sub-cohorts). A preliminary analysis of data in Experiment 1 and the discovery of not reliable and significative results for the higher pressure gradient cohort are the reason why in the other two experiments we have not considered this value. For the two first values of the pressure gradient, we use a semi-automatic tracker called *Channel and Spatial Reliability Tracker (CSRT)* that consists in first manually designing a box surrounding the cell of interest and second automatically recording the box evolution at each frame of the video [43], see Figure 2. In the videos realised with the highest pressure gradient, the mean fluid velocity is too high for the tracker to automatically track the CTCs. For this reason, they were tracked manually. By deriving the trajectories using a first-order scheme, we can directly determine the cell velocities. Note that the tracking procedure captures only the translational motion. Therefore, the data used in this work do not allow discussing possible CTCs rolling.

Table 1 summarizes the number of cells considered in each of the 17 videos. For each velocity cohort, at least 20 cells were considered. Some cells were significantly different from the others since they had very high initial velocities. We assumed that they either collided with another
(a) Experiment 1

<table><thead><tr><th>Protein modification</th><th>siCTL</th><th>siITGB1</th><th>siCD44</th><th>Total</th></tr></thead><tbody><tr><th>G</th><td></td><td></td><td></td><td></td></tr><tr><th>50.33 Pa.m<sup>-1</sup></th><td>15 (15)</td><td>15 (16)</td><td>12 (14)</td><td>42 (45)</td></tr><tr><th>100.66 Pa.m<sup>-1</sup></th><td>14</td><td>20</td><td>9 (11)</td><td>43 (45)</td></tr><tr><th>201.32 Pa.m<sup>-1</sup></th><td>24 (29)</td><td>14</td><td>14 (16)</td><td>52 (59)</td></tr><tr><th>Total</th><td>53 (58)</td><td>49 (50)</td><td>35 (41)</td><td>137 (149)</td></tr></tbody></table>

(b) Experiment 2

<table><thead><tr><th>Protein modification</th><th>siCTL</th><th>siITGA5</th><th>Total</th></tr></thead><tbody><tr><th>G</th><td></td><td></td><td></td></tr><tr><th>50.33 Pa.m<sup>-1</sup></th><td>11 (17)</td><td>12 (14)</td><td>23 (31)</td></tr><tr><th>100.66 Pa.m<sup>-1</sup></th><td>7 (11)</td><td>18 (19)</td><td>25 (30)</td></tr><tr><th>Total</th><td>18 (28)</td><td>30 (33)</td><td>48 (61)</td></tr></tbody></table>

(c) Experiment 3

<table><caption>Table 1: Total number of cells considered in each experiment, cohort and sub-cohort after removing the outliers. The values in parenthesis correspond to the number of tracked cells.</caption><thead><tr><th>Protein modification</th><th>siCTL</th><th>siITGB3</th><th>Total</th></tr></thead><tbody><tr><th>G</th><td></td><td></td><td></td></tr><tr><th>50.33 Pa.m<sup>-1</sup></th><td>11 (14)</td><td>13 (19)</td><td>24 (33)</td></tr><tr><th>100.66 Pa.m<sup>-1</sup></th><td>6 (12)</td><td>21 (23)</td><td>27 (35)</td></tr><tr><th>Total</th><td>17 (26)</td><td>34 (42)</td><td>51 (68)</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Figure 2: Example of the use of the CSRT tracker on three different cells of the video corresponding to the control case (siCTL) of Experiment 1, at a fluid velocity corresponding to the smallest pressure gradient $G^{(1)} = 50.33 Pa.m^{-1}$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 3: All extracted cell velocities over time in Experiment 1 (137 cells, without outliers). Time is given in seconds, and velocities in micrometers per seconds. Cell velocities are given in transparent color and weighted mean across cells is shown in normal color. The choice of weighted means is motivated by a better visualisation of the results, since the noise is reduced. The three figures correspond to the three different pressure gradient cohorts: top left for $G^{(1)}$, top right for $G^{(2)}$ and bottom for $G^{(3)}$. The red curves stand for the siCTL cases, green for siITGB1 cases, and blue for siCD44 cases. For each figure, the straight lines correspond to the linear regressions of the weighted means of cell velocities. The initial time of each cell velocity is adjusted in order to synchronize the oscillations. The data are normalized respectively by 100, 200 and 400 µm.s⁻¹.</figcaption></figure>
cell before entering the video time-lapse, or that they were located in a different focal plane, and we excluded them. These outliers make up only a small portion of the data, since they counted for at most 9 for each velocity cohort (42 for 278 cells ~15%), see Table 1 for details.

The extracted velocities over time cells are given in Figure 3 for Experiment 1, as an example, thus for 137 cells out of 236. Similar pictures can be obtained for the other two experiments. Individual cell velocities are in transparent color and weighted mean across cells is shown in normal color. The choice of weighted means is motivated by a better visualisation of the results, since the noise is reduced. Indeed, the raw data are preprocessed as follows. First, the velocities are normalized with respect to the fluid pressure cohort, to remove the linear dependence on the fluid velocity. Each velocity in transparent color is corrected for phase shift to synchronize the oscillations. To do so, we use the estimated cell velocities aligning their maximum points over two periods of oscillation. Finally, spurious data are filtered, to deal with the additional noise brought by the first-order derivation of velocities from positions. Indeed, some velocity values can be artificially large or low and perturb the analysis. To deal with this difficulty, we compute weighted means assigning a null weight on values above 0.85 µm.s⁻¹ and below 0.35 µm.s⁻¹. The red curves correspond to siCTL cells, green to siITGB1 cells, and blue to siCD44 cells. For each figure, the straight lines correspond to the linear regressions of the weighted means of velocities after adjustment of the initial time in order to synchronize the oscillations. The data are normalized by 100, 200 and 400 µm.s⁻¹. We normalize the data in order to facilitate the comparison between the different cohorts. The values considered will be explained later in Section 4.2. However, it should be noted that the ratio between the selected velocities is equal to the ratio of the pressure gradients.

# 3 Methods

This study aims at deciphering the influence of hydrodynamic and adhesion forces on the dynamics of CTCs moving in interaction with the wall of a microfluidic device. The methods regarding the experiments (cell culture, protein depletion, microfluidic design, etc.) can be found in [14]. First, we perform a brief statistical examination of the tracked cell velocities in Subsection 3.1. Second, a fluid velocity model is presented in Subsection 3.2. Third, a model for CTCs velocity under both hydrodynamic transport and adhesion to the wall is derived in Subsection 3.3. In Subsection 3.4, we fit the model to the experimental data using a well-designed technique for parameter estimation. Finally, in Subsection 3.5, the calibrated model for the first experiment allows to build a digital twin of flowing cells in adhesion.

All statistical analyses were performed with R. For the t-tests, we use the t_test function of the library rstatix.

## 3.1 Quick statistical analysis of the data

We perform a quick statistical analysis of these velocities between the different cohorts and subcohorts in each experiment. First, we run t-tests of the mean velocity values to determine if the differences between cohorts and subcohorts are significant. Second, we run linear regressions on the velocity values and use them to determine for which cohort and subcohort the observed decreases are significant.

## 3.2 Fluid velocity modeling

In this subsection, we derive a model for the fluid dynamics in the microfluidic device. When the viscous effects of the fluid prevail over convection, the Navier-Stokes equations can be reduced
to a Poiseuille equation. In that case, the flow shows a parabolic profile at each time, with a maximal velocity in the center of the channel decreasing to zero at the walls. In case of a time-independent pressure gradient, the Poiseuille regime is valid when the fluid verifies the following properties :

(1) it is incompressible and Newtonian ;

(2) the gravitational effect on the fluid is negligible;

(3) its flow is laminar ;

(4) and its velocity profile does not evolve over the pipe's length denoted by L in what follows.

Conditions 1 and 2 are allowed when working with a microfluidic device where the fluid is mainly comparable to water. Condition 3, can be checked by calculating the Reynolds number given by

$$Re = \frac{\rho Q D_h}{\mu S},$$

where $\rho$ is the fluid density, $Q$ the volumetric flow rate, $\mu$ the dynamic viscosity, $D_h = \frac{2(l \times h)}{l+h}$ the hydraulic diameter of a fully submerged rectangular channel and $S = l \times h$ the cross-section surface. The density can be taken as $\rho = 1.00 \times 10^3$ kg.m$^{-3}$. Based on the experiments and procedure of Osmani and coworkers in [14, 42], we have $D_h = 7.24 \times 10^{-4}$ m, $S = 1.52 \times 10^{-6}$ m$^2$ and $Q \le 5.67 \times 10^{-9}$ m$^3$.s$^{-1}$. For the dynamic viscosity, one may refer to [44] (Table 2) to obtain a close approximation of the value (different medium, but similar composition when no FBS is added). The corresponding value is $\mu = 7.31 \times 10^{-4}$ Pa.s.

It follows that $Re \le 3.75$. This value is much smaller than the critical Reynolds number for the transition from a laminar to a turbulent state, that is equal to 2600 in the case of a rectangular tube with a width eight times larger than the height, see [45]. Finally, in order to verify Condition 4, we must determine the hydrodynamic entrance length of our microfluidic device (for more details, see [46], Chapter 8, Section 8.1). For rectangular channels at laminar flow, a formula has been derived in [47]. The hydrodynamic entrance $\ell$ (in meters) is then given as a non-linear function of the aspect ratio $AR = \frac{h}{l}$ and the Reynolds number $Re$. Using the upper bound of $Re$ determined previously, a quick computation gives $\ell \le 1.53 \times 10^{-3}$ m. We can thus consider the Poiseuille flow to be fully developed (e.g. independent of the length $L$) if we observe at least $1.53 \times 10^{-3}$ m away from the pipe inlet. According to the protocol given in [42] the data was collected as close as possible to the center of the device lengthwise, which is at about $8.50 \times 10^{-3}$ m meters away from the inlet.

In our case, however, the hypothesis of a time-independent pressure gradient is not valid. In fact, the fluid dynamics is affected by the angular velocity of the pump rotor, leading to a time-dependent oscillatory perturbation term to the pressure gradient term, which we model as follows

$$G(1 + \xi_f \cos(\omega_f t + \varphi)),$$

where $\xi_f$ is the multiplicative correction amplitude, $\omega_f$ the angular velocity, and $\varphi$ the cell-dependent phase shift.

When working with an oscillating pressure gradient, the condition for the establishment of a parabolic velocity profile is strongly tied to the frequency of the oscillation relative to the viscosity of the fluid. Such relation is given through a dimensionless coefficient $Wo = \frac{h}{2} \sqrt{\frac{\omega_f \rho}{\mu}}$ introduced by Womersley in [48], which has to be inferior to 1 when the mean value of the pressure gradient is zero. A non-zero mean value will however relax such constraint, and using
the results from [49] along with supplementary observations, one can readily show a parabolic profile is obtained in our case, see Remark 3.1 and Supplementary Material A.3 for more details. The fluid velocity can therefore be written as

$$ \frac{G}{2\nu\rho}h_f^m(h-h_f^m) + \xi_f \frac{G}{\rho\omega_f}i \left( \frac{\sinh(Zh_f^m) + \sinh(Z(h-h_f^m))}{\sinh(Zh)} - 1 \right) e^{i(\omega_f t + \varphi)}, \quad (1) $$

where $i$ is the imaginary unit and $Z = (1+i)\sqrt{\frac{\omega_f}{2\nu}}$ with $\nu = \frac{\mu}{\rho}$ the kinematic viscosity. We recall that $h_f^m$ is the distance between the wall and cells in the focal plane and $h = 4.00 \times 10^{-4}$ m is the channel height. Taking the real part of this solution and performing computations (based on the linearity of the system and the principle of superposition), the fluid velocity in our context writes

$$ u_f(t) = \bar{u}_f(h_f^m) + i_f(h_f^m, \xi_f, \omega_f) \cos(\omega_f t + \varphi) + r_f(h_f^m, \xi_f, \omega_f) \sin(\omega_f t + \varphi), \quad (2) $$

where

$$ \begin{aligned} \bar{u}_f(h_f^m) &= \frac{G}{2\nu\rho}h_f^m(h-h_f^m), \\ r_f(h_f^m, \xi_f, \omega_f) &= \frac{\xi_f G}{\omega_f \rho} \mathrm{Re} \left( 1 - \frac{\sinh(Zh_f^m) + \sinh(Z(h-h_f^m))}{\sinh(Zh)} \right), \\ i_f(h_f^m, \xi_f, \omega_f) &= \frac{\xi_f G}{\omega_f \rho} \mathrm{Im} \left( 1 - \frac{\sinh(Zh_f^m) + \sinh(Z(h-h_f^m))}{\sinh(Zh)} \right). \end{aligned} $$

Note that $\bar{u}_f$ is the mean fluid velocity, while the unknown parameters are $h_f^m$, $\xi_f$, $\omega_f$ and $\varphi$.

**Remark 3.1.** Numerical approximations for the full device were also performed to test the hypothesis and investigate its limitations. These results confirming our fluid modelling can be found in the Supplementary Materials A. More specifically, we begin by verifying the value of the hydrodynamic entrance length, see Subsection A.1, and we follow with validation of the fluid expression in Subsection A.2. A final subsection A.3 focuses on how the non-zero mean value of the pressure gradient allows us to work with a Poiseuille velocity profile.

## 3.3 Cells velocity modeling

In this subsection, we define a deterministic model for cell motion based on a coupling between fluid velocity and adhesion dynamics, following previous studies [36, 37]. The interest is in capturing the different behaviours induced by varying fluid velocity and the number of expressed proteins. Both changes have an impact on cells velocity.

By denoting $N$ the bonds density and $u_c$ the cell velocity, the model writes $\forall t > 0$,

$$ \left\{ \begin{aligned} N'(t) &= c + (r - d)N, \\ u_c(t) &= u_f(t) - B(u_f(t), u_c(t)) N(t), \end{aligned} \right. \quad (3) $$

together with the initial condition $N(0) = 0$. In System (3), $c, r, d$ are given in $s^{-1}$ and stand respectively for the global binding rate, the growth rate and the unbinding rate. The function $B$ accounts for the velocity decrease arising from a unit adhesion density. All parameters are nonnegative.
**Assumption 3.1.** We assume that the adhesion parameters are time-independent and depend only on the mean fluid velocity:

$$c = c(\bar{u}_f), r = r(\bar{u}_f), \text{ and } d = d(\bar{u}_f).$$

This amounts to neglecting the effects of fluid velocity oscillations on the binding dynamics.

The function B can depend either on fluid velocity or on the cell one. Three models are considered:

1. Constant force model: $B(u_f(t), u_c(t)) = b$, where b (in µm.s⁻¹) quantifies the absolute velocity decrease induced by each unit of bonds density.

2. Fluid-dependent force model: $B(u_f(t), u_c(t)) = b u_f(t)$, with b is the dimensionless proportion of velocity decrease induced by each unit of bonds density.

3. Cell-dependent force model: $B(u_f(t), u_c(t)) = b u_c(t)$, where b is the dimensionless parameter for the friction ratio between bonds stiffness and fluid viscosity.

**Assumption 3.2.** We consider the cell-dependent model for the cell velocity equation:

$$B(u_f(t), u_c(t)) = b u_c(t).$$

Under Assumptions (3.1)-(3.2), System (3) has an explicit solution. If $d - r \neq 0$, we obtain for $t > 0$

$$N(t) = \frac{c}{d-r}(1 - e^{-(d-r)t}), \qquad (4)$$

and then

$$u_c(t) = \frac{u_f(t)}{1 - \frac{bc}{d-r}(e^{-(d-r)t} - 1)}. \qquad (5)$$

Note that at $u_c(0) = u_f(0)$, so that the cell has not formed adhesion bonds at initial time. On the other hand, experimental observations may occur only after the cell has initiated an adhesive interaction with the wall. This is why we introduce an additional parameter $\tau \ge 0$ that stands for the observation time lag. A cell with a small value of $\tau$ would be observed with an initial velocity approaching the fluid velocity, while a cell with a large value of $\tau$ would be entering the observation zone with a lower velocity.

Finally, the percentage of decrease between the cell and the fluid velocities at time $t \ge 0$ is given by the quantity $1 - \frac{u_c}{u_f}(t)$. Its limit as $t \to \infty$ then quantifies the asymptotic cell regime, and is given by

$$d_{\%} := \frac{bc}{d - r + bc}.$$

To conclude, our coupled model of parameters $h_f^m, \xi_f, \omega_f, \varphi, \tau, b, c, r$ and d reads

$$\begin{aligned} u_f(t) &= \bar{u}_f(h_f^m) + i_f(h_f^m, \xi_f, \omega_f) \cos(\omega_f t + \varphi) + r_f(h_f^m, \xi_f, \omega_f) \sin(\omega_f t + \varphi), \\ u_c(t) &= u_f(t)/(1 - \frac{bc}{d-r}(e^{-(d-r)(t+\tau)} - 1)). \end{aligned} \qquad (6)$$
## 3.4 Parameters estimation

We calibrate the model using a well-adapted estimation procedure. The main difficulties in fitting our model to the data are (a) the little information on the fluid velocity, (b) the data noise (see Figure 3), and (c) the fact that the adhesion parameters are strongly correlated with the fluid parameters. Therefore, we choose a mixed-effects parameter estimation procedure. The nonlinear mixed-effects model consists of pooling all subjects in a population and estimating a global distribution of uncertainties in the population to compensate for identifiability problems [41]. For example, the parameters of each cell $i$ could be divided into two types of uncertainties: a first part that is the same for all cells (denoted by $\theta_{pop}$ for the parameter $\theta^i$) and corresponds to the fixed effect, and a second part that represents individual variability (denoted by $\theta_{ind}^i$) and corresponds to random effects, i.e. $\theta^i = \theta_{pop} + \theta_{ind}^i$. Different covariates can also be added, e.g., for different cohorts of the population. A nonlinear mixed-effect estimation algorithm – called the stochastic approximation expectation maximization (SAEM) algorithm [50] – is implemented in the software Monolix [51]. Thanks to the R package lixoftConnectors, we could easily run Monolix using R. The code and extracted cell velocities are available here (a recent version of Monolix is required).

Since the mean fluid velocity is too high for the tracker to automatically track the CTC as the highest pressure gradient, we do not track the cells in all the 3 experiments but only for Experiment 1. However, it is useful to have a case with 3 velocities to make the estimation procedure more robust. This motivates us to propose an estimation strategy that is decomposed into two parts:

(1) Cell modeling choice and fluid reconstruction with estimation of the fluid parameters ($h_f^m$, $\omega_f$, $\xi_f$), the adhesion parameters ($bc$, $d-r$) and the time parameters ($\phi$ and $\tau$) using Experiment 1, as this is the configuration with the most available data: 3 velocities were tracked and 3 different cells (siCTL, siITGB1, siCD44);

(2) Estimation of only one fluid parameter ($h_f^m$ corresponding to the distance between the wall and the cells in the focal plane) and the adhesion parameters ($bc$, $d-r$) as well as the time parameters ($\phi$ and $\tau$) using all experiments.

In this second step, the other fluid parameters are set to the values obtained in the first step and the parameter priors were set to the estimated values in the first step. See the following two subsections for more details on each of these steps.

### 3.4.1 Fluid parameters estimations thanks to Experiment 1

In this first step, we recall that the data of Experiment 1 are used (137 trajectories).

**Fluid parameters** The only known fluid parameter is the mean pressure gradient given by $G^{(c)} = 2^{c-1}G^{(1)}$, where $\cdot^{(c)}$ denotes here and in the following the velocity cohort for $c \in \{1, 2, 3\}$, and $G^{(1)} = 50.33 \text{ Pa.m}^{-1}$. The fluid setting is going to be the same for all cells meaning that $h_f^m$, $\omega_f$ and $\xi_f$ are estimated without considering individual variability. In addition always to avoid identification problems, we strongly incorporate the information we have between the three velocity cohorts by considering

$$h_f^m^{(c)} = h_f^m^{(1)} \text{ m}, \quad \omega_f^{(c)} = 2^{c-1} \omega_f^{(1)} \text{ rad.s}^{-1}, \quad \xi_f^{(c)} = 2^{1-c} \xi_f^{(1)},$$

for $c \in \{1, 2, 3\}$. The hypothesis on $\omega_f$ is obvious and the hypothesis on $\xi_f$ implies that the oscillation amplitude of the pressure gradient is constant as $G^{(c)}\xi_f^{(c)} = G^{(1)}\xi_f^{(1)}$. Only the phase shift $\phi$ depends on the considered cell.
This implies that 142 parameters only for fluid velocity must be estimated: 3 fixed effects for $h_f^{m(1)}$, $\omega_f^{(1)}$, $\xi_f^{(1)}$, 2 fixed effects for $\varphi$ (mean and standard deviation) and 137 random effects for $\varphi$.

**Adhesion parameters** Since $b$, $c$, $r$, and $d$ are strongly paired, they can not be identified independently from the observations. We will then estimate only $bc$ and $d-r$.

The fact that we use a mixed effect approach on Experiment 1 has the great advantage that we can constrain the values of the fluid parameters, but it introduces a bias in the estimation of the individual variability of the parameters, since the fixed effects are the same for the 3 velocity cohorts and for the 3 types of cells. It is then possible to consider cohort covariates and here we decide to consider

$$ \log(bc^i) = \log(bc_{pop}) + \beta_{bc}^{(2)}[\text{if } c = 2] + \beta_{bc}^{(3)}[\text{if } c = 3] + bc_{ind}^i. $$

For the adhesion parameters $bc$, $d-r$ and $\tau$, the behaviour of the individual cell is integrated meaning that the number of parameters to be estimated is 420: 3 × 2 fixed effects (mean and std), 3 covariates for $bc$ and 3 × 137 random effects.

**Distribution laws** For parameter distributions, we consider logitnormal distributions for $h_f^m$, $\xi_f^m$, $\varphi$ and $\tau$ to keep them respectively in 0–1.5 × 10⁻⁵ m, 0–1, 0–2π, and 0–10 s. We consider lognormal distributions for the other parameters. Furthermore, whenever mixed effects are estimated, prior values and prior standard deviations are given to the SAEM algorithm. We will denote by the superscript $(\cdot)^*$ the priors and by $(\cdot)^s$ the prior standard deviations. We consider

$$ (h_f^{m(1)})^* = 7.5 \times 10^{-6} \text{ m}, (\omega_f^{(1)})^* = 12 \text{ rad.s}^{-1}, (\xi_f^{(1)})^* = 0.3, (\varphi)^* = \pi, (\varphi)^s = 1, $$

$$ (bc)^* = 0.5 \text{ s}^{-1}, (bc)^s = 1, (d-r)^* = 0.5 \text{ s}^{-1}, (d-r)^s = 1, (\tau)^* = 5 \text{ s}, (\tau)^s = 1. $$

**Model error** We consider a constant error model and we estimate the error model standard deviation starting from the prior value 25. Adding this value to the estimation of the fluid and adhesion parameters, 563 parameters are estimated during the procedure.

### 3.4.2 Estimation of parameters coupling all experiments

To facilitate comparison with the first estimation, we only indicate the differences:

* all the trajectory data are used,

* $\omega_f^{(1)}$ and $\xi_f^{(1)}$ are fixed with the estimated parameters of the first step,

* a covariate corresponding to the experiment is added for the estimation of $h_f^{m(1)}$ such that the value for Experiment 1 is fixed to the already estimated value,

* the priors values of $h_f^{m(1)}$, $bc$, $d-r$, $\varphi$ and $\tau$ are set to the mean estimates of the first step,

* and there is no random effect on $d-r$.

This second step leads to the estimation of 720 parameters: 2 covariates for $h_f^{m(1)}$ and 3 for $bc$, 1 fixed effect for $d-r$, 3 × 2 fixed effects (mean and standard deviation) and 3 × 236 random effects for $bc$, $\varphi$ and $\tau$.
## 3.5 Digital twin

The fact that our approach is based on a mathematical model allows the reproduction of a digital twin of the experimental setup – once the parameters are estimated – even for fluid configurations that were not considered during the experiments. The procedure has been tested for Experiment 1. For that purpose, we use the reconstructed fluid model at a chosen pressure gradient, set $\tau$, $\phi = 0$, and consider the estimated values of the adhesion parameters to simulate synthetic cell velocities in each subcohort for 2000 cells. Parameters with individual variability are assumed to follow a Gaussian distribution with estimated mean and standard deviation. This allows to mathematically reconstruct the velocity behavior of the cells over time from their entry into the microfluidic channel. The limitation of the experiments to videos from a specific area can thus be overcome.

We can go further by building a phenomenological law for parameters with individual variability as a function of the fluid velocity. By considering a linear interpolation and using the estimated standard deviations of these parameters we can generate synthetic cells in unobserved velocity configurations and obtain statistical quantities about the digital cohort.

# 4 Results

## 4.1 Quick statistical analysis of the data

The *p-values* resulting from the t-tests are given in Table 2. The diagonal blocks correspond to comparison between different protein modifications at the same fluid pressure, while the upper diagonal blocks account for comparison of subcohorts having the same protein modification and different fluid pressures. Significant differences can be seen between different protein modifications at the same fluid pressure and between the same protein modifications at different fluid pressure. Indeed, *p-values* are smaller than $4 \times 10^{-3}$, except between between siITGB1(3) and siCD44(3) ($p = 0.6$) and siITGB3(1) and siITGB3(2) ($p = 0.2$). Thus, at high fluid pressure gradients, the differences depleting the first or second adhesion protein are not relevant, but at smaller pressure gradients they are informative.

The outputs of the linear regressions can be found in Table 3. The intercepts – corresponding to the value of cell velocity at time $t = 0$ – increase with fluid pressure for the same protein modification, and for a given fluid pressure it increases with respect to protein modifications or it remains stable. Anyways, intercepts are always smaller than 1, value corresponding to the normalisation by 100, 200 and 400 µm.s$^{-1}$. The slope estimate and its *p-value* – related to the adhesions effects during the observation duration – show a significant decrease from the intercept value for almost all cases except for siCTL(1) and siITGA5(2) in Experiment 2 and siCTL(1) and siCTL(2) in Experiment 3.

## 4.2 Parameters estimation

**Fluid parameters** In what follows, when the parameters depend only on the velocity cohort, the 3 estimated values are given in a vector where the $c^{th}$ value corresponds to the $c^{th}$ cohort for $c \in \{1, 2, 3\}$. The results of the first estimation (using the 3 velocities but only Experiment 1, see Subsection 3.4.1 for details) allows to estimate the following values for the fluid parameters

$$ \xi_f^{(c)} = (0.27, 0.13, 0.07), \quad \omega_f^{(c)} = (12.2, 24.4, 48.8) \text{ rad.s}^{-1}, $$

and $h_f^{m(1)} = 7.2 \times 10^{-6}$ m corresponding them to Experiment 1. Using the second estimation, we obtain for $h_f^{m(1)} = 6.8 \times 10^{-6}$ m for Experiment 2 and $h_f^{m(1)} = 6.0 \times 10^{-6}$ m for Experiment 3.
(a) Experiment 1

<table><thead><tr><th></th><th>siITGB1<sup>(1)</sup></th><th>siCD44<sup>(1)</sup></th><th>siCTL<sup>(2)</sup></th><th>siITGB1<sup>(2)</sup></th><th>siCD44<sup>(2)</sup></th><th>siCTL<sup>(3)</sup></th><th>siITGB1<sup>(3)</sup></th><th>siCD44<sup>(3)</sup></th></tr></thead><tbody><tr><th>siCTL<sup>(1)</sup></th><td>&lt; 10<sup>-7</sup></td><td>&lt; 10<sup>-7</sup></td><td>&lt; 10<sup>-7</sup></td><td>-</td><td>-</td><td>&lt; 10<sup>-7</sup></td><td>-</td><td>-</td></tr><tr><th>siITGB1<sup>(1)</sup></th><td>-</td><td>2.5 &times; 10<sup>-3</sup></td><td>-</td><td>1.3 &times; 10<sup>-5</sup></td><td>-</td><td>-</td><td>&lt; 10<sup>-7</sup></td><td>-</td></tr><tr><th>siCD44<sup>(1)</sup></th><td>-</td><td>-</td><td>-</td><td>-</td><td>3.2 &times; 10<sup>-3</sup></td><td>-</td><td>-</td><td>&lt; 10<sup>-7</sup></td></tr><tr><th>siCTL<sup>(2)</sup></th><td>-</td><td>-</td><td>-</td><td>&lt; 10<sup>-7</sup></td><td>&lt; 10<sup>-7</sup></td><td>&lt; 10<sup>-7</sup></td><td>-</td><td>-</td></tr><tr><th>siITGB1<sup>(2)</sup></th><td>-</td><td>-</td><td>-</td><td>-</td><td>5.3 &times; 10<sup>-4</sup></td><td>-</td><td>&lt; 10<sup>-7</sup></td><td>-</td></tr><tr><th>siCD44<sup>(2)</sup></th><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>&lt; 10<sup>-7</sup></td></tr><tr><th>siCTL<sup>(3)</sup></th><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>4.1 &times; 10<sup>-6</sup></td><td>&lt; 10<sup>-7</sup></td></tr><tr><th>siITGB1<sup>(3)</sup></th><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>0.6</td></tr></tbody></table>

(b) Experiment 2

<table><thead><tr><th></th><th>siITGA5<sup>(1)</sup></th><th>siCTL<sup>(2)</sup></th><th>siITGA5<sup>(2)</sup></th></tr></thead><tbody><tr><th>siCTL<sup>(1)</sup></th><td>&lt; 10<sup>-7</sup></td><td>3.3 &times; 10<sup>-3</sup></td><td>-</td></tr><tr><th>siITGA5<sup>(1)</sup></th><td>-</td><td>-</td><td>1.2 &times; 10<sup>-4</sup></td></tr><tr><th>siCTL<sup>(2)</sup></th><td>-</td><td>-</td><td>&lt; 10<sup>-7</sup></td></tr></tbody></table>

(c) Experiment 3

<table><caption>**Table 2:** Resulting *p*-values from the *t*-tests performed on the mean velocity values between different protein modifications at the same fluid pressure and between the same protein modification but at different fluid pressures. Superscripts refer to the experiment number with the fluid pressure values as in Subsection 2.2. The bold values are the *p*-values inferior to $4 \times 10^{-3}$.</caption><thead><tr><th></th><th>siITGB3<sup>(1)</sup></th><th>siCTL<sup>(2)</sup></th><th>siITGB3<sup>(2)</sup></th></tr></thead><tbody><tr><th>siCTL<sup>(1)</sup></th><td>6.9 &times; 10<sup>-5</sup></td><td>6.4 &times; 10<sup>-5</sup></td><td>-</td></tr><tr><th>siITGB3<sup>(1)</sup></th><td>-</td><td>-</td><td>0.2</td></tr><tr><th>siCTL<sup>(2)</sup></th><td>-</td><td>-</td><td>9.8 &times; 10<sup>-6</sup></td></tr></tbody></table>
(a) Experiment 1

<table><thead><tr><th>Experiment</th><th>Intercept estimate</th><th>Slope estimate</th><th>Slope p-value</th></tr></thead><tbody><tr><td>siCTL<sup>(1)</sup></td><td>0.49</td><td>-0.006</td><td>&lt; 10<sup>-5</sup></td></tr><tr><td>siITGB1<sup>(1)</sup></td><td>0.57</td><td>-0.007</td><td>&lt; 10<sup>-5</sup></td></tr><tr><td>siCD44<sup>(1)</sup></td><td>0.59</td><td>-0.0035</td><td>1.5 &times; 10<sup>-2</sup></td></tr><tr><td>siCTL<sup>(2)</sup></td><td>0.54</td><td>-0.01</td><td>1.6 &times; 10<sup>-5</sup></td></tr><tr><td>siITGB1<sup>(2)</sup></td><td>0.62</td><td>-0.02</td><td>&lt; 10<sup>-5</sup></td></tr><tr><td>siCD44<sup>(2)</sup></td><td>0.61</td><td>-0.009</td><td>4.9 &times; 10<sup>-2</sup></td></tr><tr><td>siCTL<sup>(3)</sup></td><td>0.65</td><td>-0.027</td><td>1.1 &times; 10<sup>-3</sup></td></tr><tr><td>siITGB1<sup>(3)</sup></td><td>0.71</td><td>-0.054</td><td>&lt; 10<sup>-5</sup></td></tr><tr><td>siCD44<sup>(3)</sup></td><td>0.70</td><td>-0.053</td><td>&lt; 10<sup>-5</sup></td></tr></tbody></table>

(b) Experiment 2

<table><thead><tr><th>Experiment</th><th>Intercept estimate</th><th>Slope estimate</th><th>Slope p-value</th></tr></thead><tbody><tr><td>siCTL<sup>(1)</sup></td><td>0.43</td><td>-0.002</td><td>2.4 &times; 10<sup>-1</sup></td></tr><tr><td>siITGA5<sup>(1)</sup></td><td>0.52</td><td>-0.006</td><td>2.9 &times; 10<sup>-3</sup></td></tr><tr><td>siCTL<sup>(2)</sup></td><td>0.46</td><td>-0.011</td><td>1.8 &times; 10<sup>-2</sup></td></tr><tr><td>siITGA5<sup>(2)</sup></td><td>0.57</td><td>-0.0016</td><td>7 &times; 10<sup>-1</sup></td></tr></tbody></table>

(c) Experiment 3

<table><caption>**Table 3:** Linear regressions of the velocity values for all cohorts and subcohorts. Column 2: Intercept estimates (cell velocity value in the regression when $t = 0$). Column 3: Slope estimates. Column 4: *p-value* of the slope estimates. Superscripts refer to the experiment number with the fluid pressure values as in Subsection 2.2. The bold values are the *p-values* inferior to $5 \times 10^{-2}$.</caption><thead><tr><th>Experiment</th><th>Intercept estimate</th><th>Slope estimate</th><th>Slope p-value</th></tr></thead><tbody><tr><td>siCTL<sup>(1)</sup></td><td>0.42</td><td>-0.0005</td><td>7.2 &times; 10<sup>-1</sup></td></tr><tr><td>siITGB3<sup>(1)</sup></td><td>0.54</td><td>-0.014</td><td>&lt; 10<sup>-5</sup></td></tr><tr><td>siCTL<sup>(2)</sup></td><td>0.45</td><td>0.008</td><td>1.0 &times; 10<sup>-1</sup></td></tr><tr><td>siITGB3<sup>(2)</sup></td><td>0.52</td><td>-0.014</td><td>&lt; 10<sup>-5</sup></td></tr></tbody></table>
Using these estimated values, the mean velocity values can be computed using the first term of Equation (1) and we obtain

$$ \bar{u}_f^{(c)} = (99.25, 198.5, 397.0) \, \mu\text{m.s}^{-1} \\ \text{(respectively } \bar{u}_f^{(c)} = (92.9, 185.7, 371.4) \, \mu\text{m.s}^{-1} \text{ and } \bar{u}_f^{(c)} = (83.2, 166.4, 332.8) \, \mu\text{m.s}^{-1}) $$

for Experiment 1 (respectively for Experiment 2 and Experiment 3).

**Remark 4.1.** These estimated mean velocity values justify the normalization considered in Figure 3.

For the adhesion parameters, we give the results obtained at the second estimation. The parameter $d - r$ is estimated at 0.62 and Figure 4 summarizes the individual estimated values. More precisely, it corresponds to the box plots of the percentage of decrease in cell velocity $d_\%$ resulting estimating $bc$ and $d - r$. We also add p-values significance ranges of the t-tests between the estimated values of parameters for different protein modifications at the same fluid pressure. To facilitate the comparison between the same protein modifications at different fluid pressure, Figure 5 shows the same estimated values but sorted by protein modifications instead of fluid pressure gradient. The p-values ranges of the t-tests between the estimated values of parameters at different fluid pressure gradients for the same proteins are shown.

Table 4 summarises the mean and standard deviations of $d_\%$ by cohorts and subcohorts for each experiment. To facilitate the reading of this table, the mean values and the associated standard deviations are plotted in Figure 6.

Finally, Figure 7 shows numerical fits for Experiment 1 compared to the experimental data. Independently of the fluid cohort, two typical behaviours are observed: the CTCs velocity either remains stationary or decreases. Therefore, we show examples of these behaviours for each velocity cohort in the siCTL case only. Velocity values are normalised by $2^{c-1} \times 100 \, \mu\text{m.s}^{-1}$ for each $c \in \{1, 2, 3\}$.

## 4.3 Digital twin

We use the mathematical model calibrated on Experiment 1 to generate synthetic data of cell velocities. Figure 8 shows as an example the time evolution of 2000 cell velocities when $G = G^{(1)} = 50.33 \, \text{Pa.m}^{-1}$ for the 3 protein cases. Using the estimated values of the parameters summarized in Table 4, we fix $d - r$ to the mean 0.62 and randomly select $bc$ with given mean and standard deviation for each cohort. The opaque black lines correspond to the velocities of the 2000 cells and the red line to their mean. Figure 9 depicts the positions where the synthetic cell velocities fall below $35 \, \mu\text{m.s}^{-1}$ for several values of the fluid pressure gradient $G$. The color quantifies the percentage of cells in this state.

# 5 Discussion

## 5.1 Data extraction and statistical analysis

One difficulty in detecting cells was to select the correct ones. Indeed, among the cells with an apparently free trajectory, we had to pick out the CTCs without collisions, without arrests, or without problems in tracking, e.g., due to proximity to other cells. This selection obviously has an impact on the results. This selection could be automated by a more efficient tracker, especially a fully automatic tracker even for large velocities. In addition, both the tracking method and
<figure><img src="image_6.png" /><figcaption>Figure 4: Box plots of the percentage of decrease in cell velocity *d*% for the 3 experiments (red: siCTL for the three experiments, green: siITGB1 and siITGA5, blue: siCD44 and siITGB3) and for the 3 pressure gradient cohorts (top: pressure gradient fixed at *G*<sup>(1)</sup> corresponding to ~100 µm.s<sup>-1</sup>, middle: *G*<sup>(2)</sup> corresponding to ~200 µm.s<sup>-1</sup>, bottom: *G*<sup>(3)</sup> corresponding to ~400 µm.s<sup>-1</sup>). We substitute p-values with symbols such that: ns correspond to p> 10%, * to p&le; 10%, ** to p&le; 5%, *** to p&le; 0.5%, and **** to p&le; 0.05%.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 5: Box plots of the percentage of decrease in cell velocity *d*% for the 3 experiments and for the 3 pressure gradient cohorts (top: siCTL, middle: siCD44 and siITGB3, bottom: siITGB1 and siITGA5). We substitute p-values with symbols such that: ns correspond to p> 10%, * to p≤ 10%, ** to p≤ 5%, *** to p≤ 0.5%, and **** to p≤ 0.05%.</figcaption></figure>
(a) Experiment 1

<table><thead><tr><th>Experiment</th><th>Values</th><th>d%</th></tr></thead><tbody><tr><td>siCTL<sup>(1)</sup></td><td></td><td>0.54 (0.05)</td></tr><tr><td>siITGB1<sup>(1)</sup></td><td></td><td>0.47 (0.06)</td></tr><tr><td>siCD44<sup>(1)</sup></td><td></td><td>0.43 (0.05)</td></tr><tr><td>siCTL<sup>(2)</sup></td><td></td><td>0.49 (0.04)</td></tr><tr><td>siITGB1<sup>(2)</sup></td><td></td><td>0.42 (0.07)</td></tr><tr><td>siCD44<sup>(2)</sup></td><td></td><td>0.41 (0.05)</td></tr><tr><td>siCTL<sup>(3)</sup></td><td></td><td>0.38 (0.09)</td></tr><tr><td>siITGB1<sup>(3)</sup></td><td></td><td>0.38 (0.07)</td></tr><tr><td>siCD44<sup>(3)</sup></td><td></td><td>0.39 (0.07)</td></tr></tbody></table>

(b) Experiment 2

<table><thead><tr><th>Experiment</th><th>Values</th><th>d%</th></tr></thead><tbody><tr><td>siCTL<sup>(1)</sup></td><td></td><td>0.55 (0.06)</td></tr><tr><td>siITGA5<sup>(1)</sup></td><td></td><td>0.46 (0.06)</td></tr><tr><td>siCTL<sup>(2)</sup></td><td></td><td>0.53 (0.04)</td></tr><tr><td>siITGA5<sup>(2)</sup></td><td></td><td>0.39 (0.06)</td></tr></tbody></table>

(c) Experiment 2

<table><caption>**Table 4:** Mean and standard deviation values of *d%* for the all cohorts and subcohorts. The first value corresponds to the mean and the second value in parenthesis to the standard deviation.</caption><thead><tr><th>Experiment</th><th>Values</th><th>d%</th></tr></thead><tbody><tr><td>siCTL<sup>(1)</sup></td><td></td><td>0.50 (0.05)</td></tr><tr><td>siITGB3<sup>(1)</sup></td><td></td><td>0.46 (0.05)</td></tr><tr><td>siCTL<sup>(2)</sup></td><td></td><td>0.44 (0.03)</td></tr><tr><td>siITGB3<sup>(2)</sup></td><td></td><td>0.44 (0.06)</td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>Figure 6: Representation of the mean values of *d*% for the different mean fluid velocities (x-axis, given in micrometers per seconds) and for the different protein modifications in Experiment 1 (solid lines), Experiment 2 (dashed line), Experiment 3 (dotted lines). The error bars correspond to the standard deviations.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 7: Numerical fits for Experiment 1. Time is given in seconds, and velocities in micrometers per seconds. The dot-dash line represents the experimental cells velocity and in solid, the fit that we obtain with Equation (6). Top: Pressure gradient fixed at $G^{(1)}$ corresponding to ~100 µm.s⁻¹, middle: $G^{(2)}$ corresponding to ~200 µm.s⁻¹, bottom: $G^{(3)}$ corresponding to ~400 µm.s⁻¹. Left: CTCs with a stationary regime, right: CTCs with a decreasing velocity. The time interval decreases since CTCs observation duration reduces with higher velocity.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 8: Generation of 2000 synthetic cells when $G = G^{(1)} = 50.33 \text{ Pa.m}^{-1}$ for the 3 protein cases of Experiment 1 with $d-r = 0.62$, $bc \sim \mathcal{N}(m_{bc}, \sigma_{bc})$, where $m_{bc}$ and $\sigma_{bc}$ are the estimated values summarized in Table 4. In opaque black their velocities $u_c$ over 10s, calculated from Model (6) with $\tau = 0$, in red (siCTL), green (siITGB1) and blue (siCD44) the evaluated mean value over time and in dashed black the fluid velocity $u_f$ in Equation (2) with $\phi = 0$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 9: Generation of $10^5$ synthetic cells considering $G = 40, 45, 50, 55$ Pa.m$^{-1}$ and the 3 protein cases of Experiment 1 with $d - r = 0.62$ and $bc \sim N(\mu_{bc}, \sigma_{bc})$, with $\mu_{bc}$ and $\sigma_{bc}$ are taken from Table 4. For each protein modification and fluid pressure gradient $G$, we represent the position at which cells velocity goes below 35 $\mu$m/s. The percentage of cell is given in color.</figcaption></figure>
the experimental setup only resulted in the measurement of translational cell velocities, which
did not provide any information about possible rolling of the cells.

Concerning the statistical analysis of velocity data, significant differences between mean ve-
locities in almost all of the cohorts and subcohorts (*p*-values smaller than $4 \times 10^{-3}$, see Table 2)
illustrate the importance of the fluid velocity and of the CD44/ITGB3 and ITGB1/ITGA5 pro-
teins in the adhesion phenomenon. The mean velocities were not significantly different between
siITGB1(3) and siCD44(3) ($p$ = 0.6), probably due to the higher fluid velocity impeding adhesion,
and siITGB3(1) and siITGB3(2) ($p$ = 0.2). This is also supported by the linear regressions on
the velocities, see Table 3. In each cohort and subcohort, the intercept smaller than 1 indicates
the presence of adhesion, higher for lower fluid velocities and for siCTL cases. The slopes show a
decelerating dynamics for 13 out of 17 subcohorts (*p*-values less than $5 \times 10^{-2}$). These values are
mostly negative but close to zero, indicating the presence of stationary velocity profiles, see Fig-
ure 7. Non-significant *p*-values slopes in Experiment 2 and 3 are related to a lower distance $h_f^m$
of CTCs from the wall. Indeed, in this case cells have already experienced adhesion as it comes
out from the lower intercept estimate values. The goal of the mathematical model – presented
in Subsections 3.2 and 3.3 – was to understand and quantify these initial observations.

5.2 Fluid modeling

As for the fluid modeling, the choice of the pressure gradient under the cosinusoidal form has
an important impact, since it can be shown that in our context the imaginary part *i*<sub>*f*</sub> is larger
than the real part *r*<sub>*f*</sub>, which means that most of the oscillations of the cell velocities are under
the cosinusoidal form. To make more complex assumptions, a better knowledge of the pump is
needed.

As for the estimation of fluid parameters done mostly in the first step of the estimation pro-
cess (see Subsection 3.4.1), the strong correlation between fluid and adhesion parameters leads
to identification problems. For example, we observed that multiple values of the parameter pair
(*bc*, *ū*<sub>*f*</sub>(*h*<sub>*f*</sub><sup>*m*</sup>) resulted in similar fits. We then selected very carefully the priors of the fluid param-
eters using the literature. The prior value of *h*<sub>*f*</sub><sup>*m*</sup> (only parameter appearing in the mean velocity
values *ū*<sub>*f*</sub>) has a strong impact on the results. The value has been selected to be coherent with
the velocity measures with single-particles performed in [13] (see Figure 6). Its estimation led to
mean velocity values very close to the estimated ones in [42], which are 100, 200 and 400 µm.s<sup>-1</sup>.
The prior of the angular velocity *ω*<sub>*f*</sub> and its dependence on the pressure gradient can be derived
directly from Figure 3-Top-Left (10 oscillations 5 s gives ~ 10/5 × 2π ~12.6 s<sup>-1</sup>). As for the
correction amplitude *ξ*<sub>*f*</sub>, which appears in the amplitudes *r*<sub>*f*</sub> and *i*<sub>*f*</sub> in Equation (2), we have
halved its value in each successive cohort. This is confirmed by the amplitudes observed in Fig-
ures 3 and 7. This last hypothesis was also confirmed by looking at the AIC for the case where
*ξ*<sub>*f*</sub> was constant (AIC=130019.6) as opposed to *ξ*<sub>*f*</sub> depending on fluid pressure (AIC=129908.7).
Of course, more complex hypotheses could be tested, but this would also require a better under-
standing of how the pump works.

5.3 Adhesion modeling

**Modeling** The adhesion dynamics is described by an ODE on the adhesion density. The modeling choice is built on Assumptions (3.1)-(3.2). Assumption (3.1) states constant binding rates in each velocity cohort. By doing so, the effects of velocity oscillations on the adhesion dynamics are neglected. Several biophysical studies investigated the relation between the load applied on a cell and its binding dynamics, and introduced catch bonds or slip bonds (see e.g [25–27]). At a lower scale, local binding rates can be given as a function of the membrane
distance to the wall [52], which could be averaged in a macroscopic rate decreasing with cell velocity. Catch bonds are well described in particular for integrins-based adhesions with a load-based reinforcement [53, 54]. In our model, this can be described by either a reproduction rate $r$ or a force parameter $b$ increasing with $u_f$. Finally, increased dissociation rates with shear force (slip bonds) are commonly considered and validated in other experimental settings [20, 52, 55]. However, these characterizations arise from lower scale experiments, depicting the effect of continuous variations in shear forces [56–58]. In our context, it relates to the effects of fluid velocity oscillations on the adhesion dynamics. Some studies have highlighted the importance of oscillatory forces on the adhesion dynamics [59–61], so that it could have been relevant in our modeling approach. In the end, our choice of neglecting fluid oscillations in the binding terms followed from the experimental data. As can be seen in Figure 7, the cell velocity oscillations were very regular, even in decelerating cells and among cells. This strongly suggested that the oscillating fluid forces variations at this timescale were not enough to disrupt the macroscopic cell regime. Assumption (3.1) then enables to derive an explicit solution for the bonds density over time given in Equation (4). It was found during the estimation of parameters that $d-r > 0$. In this case, the bonds density increases exponentially and then saturates at $c/(d-r)$, while the asymptotic cell regime $d\%$ can not reach 1, preventing cell arrest. This shows that our modeling is not suitable to account for arrested cells.

In our model, cell velocity is given by the difference between the fluid velocity and an adhesion term. This formulation is classically found in other modeling approaches, see e.g [29, 31, 36, 37]. In a macroscopic setting, the adhesion term is proportional to the closed bonds density, and involves both geometric features and forces and torques exerted by the bonds. Note that the linearity in adhesion density intertwine adhesion dynamics and force parameters. Rather than making this term explicit, we kept a minimal framework and compared three expressions.

The constant force model, given by $u_c(t) = u_f(t) - bN(t)$, accounts for constant binding forces at the cell scale, see [37]. The fluid-dependent force model writes $u_c(t) = u_f(t)(1 - bN(t))$ and can be seen as a simplification of a model for elastic bonds given in [36]. Finally, the cell-dependent force model writes $u_c(t) = u_f(t) - bu_c(t)N(t)$, leading to Equation (5). This framework can be seen as a macroscopic viewpoint for the average force exerted by an elastic bond over its lifetime. Several studies show matching microscopic viewpoints involving structured bonds density capturing bonds elongation. In [29], the membrane at the cell rear moves away from the wall at a normal velocity proportional to $u_c$. In [34–36], elongation is the product of the bond’s age and the cell velocity. Furthermore, in [34, 35], a scaling limit in fast bonds turnover and rigid forces allowed to justify rigorously a macroscopic adhesion term proportional to $u_c(t)$.

Assumption (3.2) consists in the choice of the cell-dependent force model over the others. Indeed, the associated deceleration term rewrites $\frac{bu_f}{1+bN}N$. From the biological viewpoint, it depicts an adhesion force proportional to $u_f$ (as for catch bonds) that is shared among bonds, as in adhesion clusters [52, 62]. In contrast, the two other models, less realistic, were based on additive individual adhesion forces, either constant or proportional to $u_f$. Finally, this assumption was validated by a better adequation with experimental data in the first estimation step. The fluid-dependent force model (AIC=129910.6) slightly differs from the cell-dependent one (AIC=129908.7) in AIC, suggesting a weak importance of load sharing in the observed regime. Using the constant force model strongly degrades the AIC (AIC=130066.8) showing that this assumption should be excluded.

**Parameters estimation** Concerning the adhesion parameters, we could only estimate $bc$ and $d-r$ for the chosen cell model. In contrast with the fluid whose behaviour is the same in all experiments, cell parameters vary from cell to cell. They were taken as the sum of a fixed population effect and of an individual random term. Several attempts have been realised to obtain
the optimal parameters from which we could interpret the biological phenomenon at study. In particular, we have considered a covariate model for *bc* in order to integrate the cohorts effects. One could also add covariates with respect to subcohorts, but this seems only to complexify the definition of *bc* without improving the AIC. The same procedure can be applied to *d* − *r*. Again we have tried this strategy in the first estimation step adding covariates respect to cohorts, subcohorts and adhesion group and both fixing *d*−*r* in the population and considering individual variability. We got either similar or higher AICs (superior to 129920) in estimated parameters. Consequently, we considered only a fixed effect on *d*−*r*. These choices aim to reduce the number of estimated parameters, while making sure to have enough decelerating cells in each estimation group.

Calibrating mathematical models on noisy and partial data is a challenging task. Both the model and the estimation strategy should be carefully adapted to the nature and scale of the data in order to provide a reliable framework for biological interpretation. Several of our choices aimed at ensuring the relevance of the obtained parameters. First, the preliminary estimation of fluid parameters allowed to reconstruct the better understood experimental system before investigating the cell behaviors. Moreover, the mixed-effect modeling, while increasing the number of parameters, constrains them with the knowledge that the observed cells behave following a common framework. This approach improves identifiability and reduces the risk of obtaining spurious parameter sets in the population. Some results were also indicative of coherent estimations: the obtained cell heights were no more than 1.2 µm apart one from another, and they were in accordance with the expected values from the experimental settings. Furthermore, while the three experiments were conducted independently, the estimation results allowed to recover the known molecular associations across them, at least fully for the lower fluid velocity. Altogether, these estimations provide a credible framework with respect to the biological phenomenon under study.

**Interpretation of the estimated adhesion parameters** A better adhesion parameters interpretability is given by the percentage of decrease in cells velocity *d*%, depending both on *bc* and *d*−*r*. Figure 4 shows significant differences for low velocities in all the three experiments. As fluid velocity increases, we lose significant differences since adhesion is less important. This also justifies the interest in reproducing experiments under lower fluid velocities through the digital twin as in Figure 9.

In Figure 5, a comparison is done over different fluid velocities. In particular we compare control cases and similar proteins over the three experiments. Significant differences are observed whenever adhesion is meaningful, i.e. control (red boxes) and stronger adhesion protein depletion (green boxes). This stands for a more important contribution in the adhesion phenomenon of the first – albeit weak– protein CD44 as well as ITGB3.

Moreover, in both figures we can also compare the different control cases effect (which is
similar as expected) as well as proteins one. Proteins CD44 and ITGA5 (siITGB1 and siITGA5
case respectively) have an active role in the first steps of adhesion, whereas ITGB1 and ITGB3
(siCD44 and siITGB3 case respectively) has a minor impact.

As stressed before, this analysis has been pursued for *d* − *r* fixed. In this context, individ-
ual variability does not change the AIC, instead it can cause the loss of significant differences
among proteins adhesion effects. A constant value indicates that the timescale of deceleration
(or equivalently of global adhesion growth) is identical in each case. In conclusion, the choice
on *d* − *r* estimation is based on better decelerating results in the second estimation step. Note
that deciphering further experiments with a wider fluid velocities range or a longer observation
duration could require relaxing this choice. Consequently, the cell variability between cohorts
and subcohorts is carried by the binding term *bc*.
## 5.4 Comparison with the literature

Conclusions made from the estimation results are tightly linked to the experimental setting, considered cells and proteins, and scale of measures (see e.g [63] for variable behaviors among tumor cells in flow). One can nonetheless discuss their coherence with respect to the literature. A deceleration timescale independent of shear rate was also found in [40] for other adhesion bonds. Note that this is not in contradiction with the concepts of catch and slip bonds, since having an increasing dissociation rate *d* and a decreasing reinforcement rate *r* with respect to the fluid velocity may maintain a constant difference *d* − *r* if the magnitudes are similar. However, the scaling of present experiments were not suited for capturing these features.

In our work, we find that the fast and metastable CD44/ITGB3 bonds are more important for CTCs deceleration *in vitro* than slower and more stable ITGB1/ITGA5 ones, which could not be established in the prior analysis. Nevertheless, the *in vivo* experiments on the zebrafish embryo had underlined that only the first type of bonds was able to induce early cell arrest, indicating a higher adhesion efficacy at faster flows [14]. The importance of CD44 in adhesion is already established in more general settings [64–68]. We further find that the binding parameter *bc* is a decreasing function of fluid velocity in the control case, and weaker in the siITGB1 case, but not in the siCD44 case. As a result it can be considered as a feature of CD44-based adhesion, yielding increasing force parameter or binding rate with decreasing velocities. Note that opposite trends have also been observed in other context, such as an increase in CD44-based adhesion of leukemic cells above a shear stress threshold [69]. It is also known that CD44 has a larger binding rate than ITGB1, but forms weaker and shorter-lived bonds [14]. In combination with our findings, a possible framework could be that, at a given fluid velocity level, CD44 is associated with larger dissociation and reinforcement rates than ITGB1, and that the larger binding rate compensate for the weaker forces. Conversely, at the same fluid velocity level, ITGB1/ITGA5 bonds may be associated with larger forces but much lower binding rates, together with lower dissociation and reinforcement rates yielding close values for *d* − *r*. Integrins being tightly linked with cell cytoskeleton, they are particularly responsive to load and involved in many biological processes involving mechanotransduction [70–72]. They are primarily considered as mediating firm adhesion, and behave as catch-slip bonds within the relevant speed range [53, 54, 68, 73]. Altogether, our findings are therefore compatible with low scale studies on cell adhesion.

Inflow cell adhesion has been investigated in various *in vitro* settings for fibroblasts or immune cells [74–77]. Few microfluidics studies have investigated tumor cells adhesion. In the absence of flow, cell types heterogeneity in adhesion strenght has been pointed out in single cell measures in [78]. The role of flow rate and acceleration on resting tumor cells detachment in a microchannel has been highlighted in [38]. The importance of chemokines in extravasation has been studied in [79, 80]. Finally, two microfluidic approaches investigated the adhesion dynamics of MDA-MB-231 breast cancer cells [40, 81]. In [81], a wide range of shear stresses applied on 500 to 1500 resting cells showed that the shear stress threshold for detachment was unchanged among all conditions, whereas detachment kinetics were different. In [40], cells were flowing and interacting with an anti-EpCAM ligands-coated wall. Sequential fitting of the computational model of [31] on mean translational velocities for several shear rates allowed to identify the average cell height, binding force, and bonds spring constant. Normalized cell velocities for all shear rates were successfully fitted by a generic exponentially decreasing curve, suggesting a strong dependence of the cell velocity magnitude on the fluid velocity, whereas it was not the case for the typical decay time. Several differences exist between our frameworks. From the biological viewpoint, the wall in [40] is passive, while ours is a monolayer of endothelial cells, whose flow-driven active behaviour in cell arrest has been observed [14]. Together with lower fluid velocities, it may explain their measures of smooth velocity decays until cell arrests, that were not observed
in our case. Moreover, the EpCAM-based approach induces a bias towards cells with strong
epithelial identity, which are not necessarily the most aggressive (see e.g [7, 82]). Furthermore,
we worked with partial observations since cell velocities were not measured from their entrance
in the experimental setting. On the other hand, we took into account the temporally oscillating
fluid velocity that needed to be reconstructed and developed a mixed-effects calibration strategy
that was able to deal with the individual cell velocities over time and to estimate the parameters
jointly. In this setting, we obtained insights on the role of the fluid velocity that are consistent
with previous observations [13, 14, 40]. Finally, our original and robust approach allowed not
only to investigate the respective roles of several adhesion proteins in the cell dynamics, but also
to build a digital twin of flowing cells with adhesion in other fluid configurations.

Therefore, our methodology represents a valuable contribution not only to validate experi-
mental results, but also to interpret them more finely and use them to build further knowledge.

**6 Conclusion and perspectives**

In this work, we have attempted to characterise CTCs in the flow and their interaction with the
vessel wall, relying on the *in vitro* experiments performed by Osmani and collaborators in [14, 42].
Whereas previous analyses focused on cell arrest, the use of the CSRT tracker allowed us to record
trajectories and velocities of individual cells. We were able to analyse 3 experiments comprising
different cell cohorts with respect to three different values of fluid pressure gradient (below the
threshold for efficient CTC adhesion found by Osmani and collaborators in [13]) and 7 different
protein expressions (siCTL, the control case relative to each of the three experiment; siITGB1,
depletion of ITGB1, integrin that promotes adhesion stabilisation; siCD44, depletion of CD44,
protein involved in early adhesion; siITGA5, depletion of ITGA5, whose effect is similar to
ITGB1; siITGB3, depletion of ITGB3, with a comparable role to CD44). Statistical analysis
of the mean of the extracted cell velocities and linear regression allowed the observation of a
slowing behaviour over time, see Tables 2 and 3. This shows that adhesion is a continuous-time
phenomenon involving CTCs in a fluid with a velocity below the threshold of 400 µm.s⁻¹.

Since the fluid velocity was not measured directly, we only knew the values of the pressure
gradient generated by the peristaltic pump that made up the device. This lack of data was
compounded by our lack of knowledge about the pump. However, we were able to establish a
Poiseuille regime and describe the fluid velocity as a combination of oscillatory functions induced
by the pump and evident in the tracked cell velocity in Figure 3. We then focus on the modeling
of the cell velocity. The oscillating Poiseuille flow was weakly coupled to a simple ODE model
for cell adhesion that describes the cell velocity as the fluid velocity affected by bond formation
and disruption.

Optimal parameters for our model were not easy to find. Indeed, there are practical problems
with identifiability, mainly due to data noise and little information about the fluid parameters.
Our strategy to overcome this problem is based on a mixed-effects model and careful selection
of fluid parameter priors.

The well-designed parameter estimation has led to very attractive results. The first important
result is that the estimated parameter values allow deciphering the CTC binding and deriving
biological outcomes. We briefly summarize them here: (1) the strong influence of fluid velocity on
protein binding: a low fluid velocity favours a decrease in cell velocity and the formation of bonds,
while a high fluid velocity makes it difficult to observe this adhesion phenomenon even when all
adhesion proteins are expressed, (2) the importance of CD44/ITGB3 in early adhesion compared
to the integrin ITGB1/ITGA5. Similar conclusions have already been assessed for point (1) using
statistical strategies. However, point (2) was only assessed using *in vivo* experiments and could
not be extracted from the *in vitro* experiments prior to our work. This underlines the quality of
the strategy we developed, which is based on mathematical modeling and data assimilation.

This methodology leads to a second important result. Indeed, since it is based on a mathematical model – once the model parameters have been estimated – the digital twin of the setup can be reproduced, making it possible to learn more details about the behavior of the cells, even for configurations that were not considered during the experiments.

In conclusion, it can be noted that this experimental setup is quite simple and far from
reality, but from a mathematical point of view, its adoption allows the selection of a model
and its parameterization. It should be noted that even in this case, where a lot of data can be
extracted, some parameters are not identifiable. With less data, such as only the percentage of
locally arrested cells, as is often the case with *in vivo* data, it would be impossible to identify
the estimated parameters. This first step is therefore crucial. Indeed, these first estimates make
it possible to establish laws for the parameters as a function of fluid velocity and to perform
simulations even for configurations not considered in the experiments, thus paving the way for
*in vivo* data.

As for the perspectives, the first one concerns the improvement of the CTC tracker from the
experimental videos, since it is not fully automatic and has a great need of optimization. The
second is the development of a mathematical model adapted to the *in vivo* experiments. These
data should allow us to incorporate cell arrest into our model. The more complex geometry
will require a more complex model of blood circulation. If we consider a different experimental
setup, this model of flow velocity should be adapted based on the very extensive literature on the
subject. In particular, the fact that blood is a non-Newtonian fluid – as its viscosity depends on
the shear rate – must be taken into account. Finally, from a biological standpoint, it is possible
to use the model presented to study and predict additional molecular modes involved in the
arrest of CTCs at the vascular wall.

**Author contributions** AC and CE designed the study based on the biological data generated by NO in the team of JGG. GC and JG analyzed the data. AC and GC implemented the software code. AC, CE and GC interpreted the results. AC, CE and GC wrote the manuscript and JG the supplementary materials.

Acknowledgements

The work of JGG and NO has been funded by Plan Cancer 2014-2019 (OptoMetaTrap), CNRS IMAG'IN (JGG) and by institutional funds from INSERM and the University of Strasbourg.

Declaration of interest

The authors declare no competing interests.

References

[1] H. Dillekås, M. S Rogers, and O. Straume. Are 90% of deaths from cancer caused by metastases? *Cancer Med.*, 8(12):5574-5576, 2019. URL https://doi.org/10.1002/cam4.2474.

[2] P. S. Steeg. Targeting metastasis. *Nat. Rev. Cancer*, 16(4):201–218, 2016. URL https://doi.org/10.1038/nrc.2016.25.
[3] D. Wirtz, K. Konstantopoulos, and P. C. Searson. The physics of cancer: the role of physical interactions and mechanical forces in metastasis. *Nat. Rev. Cancer*, 11(7):512–522, 2011. URL https://doi.org/10.1038/nrc3080.

[4] G. Follain, D. Herrmann, S. Harlepp, V. Hyenne, N. Osmani, S. C. Warren, P. Timpson, and J. G. Goetz. Fluids and their mechanics in tumour transit: shaping metastasis. *Nat. Rev. Cancer*, 20(2):107–124, 2020. URL https://doi.org/10.1038/s41568-019-0221-x.

[5] A. C. Obenauf and J. Massagué. Surviving at a distance: organ-specific metastasis. *Trends Cancer*, 1(1):76–91, 2015. URL https://doi.org/10.1016/j.trecan.2015.07.009.

[6] S. B. Lim, W. Di Lee, J. Vasudevan, W.-T. Lim, and C. T. Lim. Liquid biopsy: one cell at a time. *NPJ Precis. Oncol.*, 3(1):23, 2019. URL https://doi.org/10.1038/s41698-019-0095-0.

[7] M. Peralta, N. Osmani, and J. G. Goetz. Circulating tumor cells: Towards mechanical phenotyping of metastasis. *iScience*, 2022. URL https://doi.org/10.1016/j.isci.2022.103969.

[8] S. A. Joosse, T. M. Gorges, and K. Pantel. Biology, detection, and clinical implications of circulating tumor cells. *EMBO Mol. Med.*, 7(1):1–11, 2015. URL https://doi.org/10.15252/emmm.201303698.

[9] J. Liu, J. Lian, Y. Chen, X. Zhao, C. Du, Y. Xu, H. Hu, H. Rao, and X. Hong. Circulating tumor cells (ctcs): a unique model of cancer metastases and non-invasive biomarkers of therapeutic response. *Front. Genet.*, 12:734595, 2021. URL https://doi.org/10.3389/fgene.2021.734595.

[10] A. Ring, B. D. Nguyen-Straeuli, A. Wicki, and N. Aceto. Biology, vulnerabilities and clinical applications of circulating tumour cells. *Nat. Rev. Cancer*, 23(2):95–111, 2023. URL https://doi.org/10.1038/s41568-022-00536-4.

[11] Q.-Q. Huang, X.-X. Chen, W. Jiang, S.-L. Jin, X.-Y. Wang, W. Liu, S.-S. Guo, J.-C. Guo, and X.-Z. Zhao. Sensitive and specific detection of circulating tumor cells promotes precision medicine for cancer. *J. Cancer Metastasis Treat.*, 5:34, 2019. URL https://doi.org/10.20517/2394-4722.2018.94.

[12] K. A. Rejniak. Circulating Tumor Cells: When a Solid Tumor Meets a Fluid Microenvironment. In K. A. Rejniak, editor, *Systems Biology of Tumor Microenvironment*, volume 936, pages 93–106. Springer International Publishing, 2016. URL https://doi.org/10.1007/978-3-319-42023-3_5.

[13] G. Follain, N. Osmani, A. S. Azevedo, G. Allio, L. Mercier, M. A. Karreman, G. Solecki, M. J. García-León, O. Lefebvre, N. Fekonja, et al. Hemodynamic forces tune the arrest, adhesion, and extravasation of circulating tumor cells. *Dev. Cell*, 45(1):33–52, 2018. URL https://doi.org/10.1016/j.devcel.2018.02.015.

[14] N. Osmani, G. Follain, M. J. García-León, O. Lefebvre, I. Busnelli, A. Larnicol, S. Harlepp, and J. G. Goetz. Metastatic tumor cells exploit their adhesion repertoire to counteract shear forces during intravascular arrest. *Cell Rep.*, 28(10):2491–2500, 2019. URL https://doi.org/10.1016/j.celrep.2019.07.102.

[15] C. D. Paul, K. Bishop, A. Devine, E. L. Paine, J. R. Staunton, S. M. Thomas, J. R.
Thomas, A. D. Doyle, L. M. Miller Jenkins, N. Y. Morgan, et al. Tissue architectural cues drive organ targeting of tumor cells in zebrafish. *Cell Syst.*, 9(2):187-206, 2019. URL https://doi.org/10.1016/j.cels.2019.07.005.

[16] I. Martínez-Pena, P. Hurtado, N. Carmona-Ule, C. Abuín, A. B. Dávila-Ibáñez, L. Sánchez, M. Abal, A. Chaachou, J. Hernández-Losa, S. R. Cajal, et al. Dissecting breast cancer circulating tumor cells competence via modelling metastasis in zebrafish. *Int. J. Mol. Sci.*, 22(17):9279, 2021. URL https://doi.org/10.3390/ijms22179279.

[17] G. I. Bell. Models for the specific adhesion of cells to cells. *Science*, 200(4342):618-627, 1978. URL https://doi.org/10.1126/science.347575.

[18] G. I. Bell, M. Dembo, and P. Bongrand. Cell adhesion. competition between nonspecific repulsion and specific bonding. *Biophys. J.*, 45(6):1051-1064, 1984. URL https://doi.org/10.1016/S0006-3495(84)84252-6.

[19] T. Erdmann and U. S. Schwarz. Stability of adhesion clusters under constant force. *Phys. Rev. Lett.*, 92(10):108102, 2004. URL https://doi.org/10.1103/PhysRevLett.92.108102.

[20] T. Erdmann and U. S. Schwarz. Adhesion clusters under shared linear loading: a stochastic analysis. *EPL*, 66(4):603, 2004. URL https://doi.org/10.1209/epl/i2003-10239-3.

[21] T. Erdmann and U. S. Schwarz. Impact of receptor-ligand distance on adhesion cluster stability. *Eur. Phys. J. E*, 22:123-137, 2007. URL https://doi.org/10.1140/epje/e2007-00019-8.

[22] D. A. Hammer and D. A. Lauffenburger. A dynamical model for receptor-mediated cell adhesion to surfaces. *Biophys. J.*, 52(3):475-487, 1987. URL https://doi.org/10.1016/S0006-3495(87)83236-8.

[23] M. Dembo, D. C. Torney, K. Saxman, and D. Hammer. The reaction-limited kinetics of membrane-to-surface adhesion and detachment. *Proc. Royal Soc. B*, 234(1274):55-83, 1988. URL https://doi.org/10.1098/rspb.1988.0038.

[24] D. A. Hammer and S. M. Apte. Simulation of cell rolling and adhesion on surfaces in shear flow: general results and analysis of selectin-mediated neutrophil adhesion. *Biophys. J.*, 63 (1):35-57, 1992. URL https://doi.org/10.1016/S0006-3495(92)81577-1.

[25] K.-C. Chang, D. Tees, and D. A. Hammer. The state diagram for cell adhesion under flow: leukocyte rolling and firm adhesion. *Proc. Natl. Acad. Sci.*, 97(21):11262-11267, 2000. URL https://doi.org/10.1073/pnas.200240897.

[26] K. E. Caputo, D. Lee, M. R. King, and D. A. Hammer. Adhesive dynamics simulations of the shear threshold effect for leukocytes. *Biophys. J.*, 92(3):787-797, 2007. URL https://doi.org/10.1529/biophysj.106.082321.

[27] M. T. Beste and D. A. Hammer. Selectin catch-slip kinetics encode shear threshold adhesive behavior of rolling leukocytes. *Proc. Natl. Acad. Sci.*, 105(52):20716-20721, 2008. URL https://doi.org/10.1073/pnas.0808213105.

[28] A. Efremov and J. Cao. Bistability of cell adhesion in shear flow. *Biophys. J.*, 101(5): 1032-1040, 2011. URL https://doi.org/10.1016/j.bpj.2011.07.026.
[29] L. Li, H. Tang, J. Wang, J. Lin, and H. Yao. Rolling adhesion of cell in shear flow: A theoretical model. *J. Mech. Phys. Solids*, 119:369-381, 2018. URL https://doi.org/10.1016/j.jmps.2018.07.013.

[30] L. Li, W. Kang, and J. Wang. Mechanical model for catch-bond-mediated cell adhesion in shear flow. *Int. J. Mol. Sci.*, 21(2):584, 2020. URL https://doi.org/10.3390/ijms21020584.

[31] C. B. Korn and U. S. Schwarz. Dynamic states of cells adhering in shear flow: from slipping to rolling. *Phys. Rev. E*, 77(4):041904, 2008. URL https://doi.org/10.1103/PhysRevE.77.041904.

[32] S. Reboux, G. Richardson, and O. E. Jensen. Bond tilting and sliding friction in a model of cell adhesion. *Proc. Royal Soc. A*, 464(2090):447-467, 2008. URL https://doi.org/10.1098/rspa.2007.0210.

[33] L. Preziosi and G. Vitale. A multiphase model of tumor and tissue growth including cell adhesion and plastic reorganization. *Math. Models Methods Appl. Sci.*, 21(09):1901-1932, 2011. URL https://doi.org/10.1142/S0218202511005593.

[34] V. Milišić and D. Oelz. On the asymptotic regime of a model for friction mediated by transient elastic linkages. *J. Math. Pures Appl.*, 96(5):484-501, 2011. URL https://doi.org/10.1016/j.matpur.2011.03.005.

[35] V. Milišić and D. Oelz. On a structured model for load-dependent reaction kinetics of transient elastic linkages mediating nonlinear friction. *SIAM J. Math. Anal.*, 47(3):2104-2121, 2015. URL http://epubs.siam.org/doi/10.1137/130947052.

[36] B. Grec, B. Maury, N. Meunier, and L. Navoret. A 1d model of leukocyte adhesion coupling bond dynamics with blood velocity. *J. Theor. Biol.*, 452:35-46, 2018. URL https://doi.org/10.1016/j.jtbi.2018.02.021.

[37] C. Etchegaray and N. Meunier. A stochastic model for cell adhesion to the vascular wall. *J. Math. Biol.*, 79:1665-1697, 2019. URL https://doi.org/10.1007/s00285-019-01407-7.

[38] L. S. L. Cheung, X. Zheng, A. Stopa, J. C. Baygents, R. Guzman, J. A. Schroeder, R. L. Heimark, and Y. Zohar. Detachment of captured cancer cells under flow acceleration in a bio-functionalized microchannel. *Lab Chip*, 9(12):1721-1731, 2009. URL https://doi.org/10.1039/B822172C.

[39] L. S. L. Cheung, X. Zheng, L. Wang, R. Guzman, J. A. Schroeder, R. L. Heimark, J. C. Baygents, and Y. Zohar. Kinematics of specifically captured circulating tumor cells in bio-functionalized microchannels. *J. Microelectromech. Syst.*, 19(4):752-763, 2010. URL https://doi.org/10.1109/JMEMS.2010.2052021.

[40] L. S. L. Cheung, X. Zheng, L. Wang, J. C. Baygents, R. Guzman, J. A. Schroeder, R. L. Heimark, and Y. Zohar. Adhesion dynamics of circulating tumor cells under shear flow in a bio-functionalized microchannel. *J. Microelectromech. Syst.*, 21(5):054033, 2011. URL https://doi.org/10.1088/0960-1317/21/5/054033.

[41] M. Lavielle. *Mixed effects models for the population approach: models, tasks, methods and tools.* CRC press, 2014.
[42] N. Osmani, G. Follain, V. Gensbittel, M. J. García-León, S. Harlepp, and J. G. Goetz. Probing intravascular adhesion and extravasation of tumor cells with microfluidics. *Metastasis: Methods and Protocols*, pages 111–132, 2021. URL https://doi.org/10.1007/978-1-0716-1350-4_8.

[43] A. Lukežič, T. Vojíř, L. Zajc, J. Matas, and M. Kristan. Discriminative correlation filter tracker with channel and spatial reliability. *Int. J. Comput. Vis.*, 126(7):671–688, 2018. URL https://doi.org/10.1007/s11263-017-1061-3.

[44] C. Poon. Measuring the density and viscosity of culture media for optimized computational fluid dynamics analysis of in vitro devices. *J. Mech. Behav. Biomed. Mater.*, 126:105024, 2022. URL https://doi.org/10.1016/j.jmbbm.2021.105024.

[45] T. W. Kao and C. Park. Experimental investigations of the stability of channel flows. part 1. flow of a single liquid in a rectangular channel. *J. Fluid Mech.*, 43(1):145–164, 1970. URL https://doi.org/10.1017/S0022112070002288.

[46] F. P. Incropera, D. P. DeWitt, T. L. Bergman, A. S. Lavine, et al. *Fundamentals of heat and mass transfer*, volume 6. Wiley New York, 1996.

[47] G. Ferreira, A. Sucena, L. L. Ferrás, F. T. Pinho, and A. M. Afonso. Hydrodynamic entrance length for laminar flow in microchannels with rectangular cross section. *Fluids*, 6(7):240, 2021. URL https://doi.org/10.3390/fluids6070240.

[48] J. R. Womersley. Method for the calculation of velocity, rate of flow and viscous drag in arteries when the pressure gradient is known. *Physiol. J.*, 127(3):553, 1955. URL https://doi.org/10.1113/jphysiol.1955.sp005276.

[49] H.-L. Ma and C.-H. Kuo. Theoretical analysis of an oscillatory plane poiseuille flow—a link to the design of vortex flow meter. *Phys. Fluids*, 29(5):053602, 2017. URL https://doi.org/10.1063/1.4984006.

[50] E. Kuhn and M. Lavielle. Maximum likelihood estimation in nonlinear mixed effects models. *Comput. Stat. Data Anal.*, 49(4):1020–1038, 2005. URL https://doi.org/10.1016/j.csda.2004.07.002.

[51] Monolix 2021R2, Lixoft SAS, a Simulations Plus company. URL https://monolix.lixoft.com.

[52] C. Zhu, G. Bao, and N. Wang. Cell Mechanics: Mechanical Response, Cell Adhesion, and Molecular Deformation. *Annu. Rev. Biomed. Eng.*, 2(1):189–226, 2000. URL https://doi.org/10.1146/annurev.bioeng.2.1.189.

[53] N. Strohmeyer, M. Bharadwaj, M. Costell, R. Fässler, and D. J. Müller. Fibronectin-bound α5β1 integrins sense load and signal to reinforce adhesion in less than a second. *Nat. Mater*, 16(12):1262–1270, 2017. URL https://doi.org/10.1038/nmat5023.

[54] M. Bharadwaj, N. Strohmeyer, G. P. Colo, J. Helenius, N. Beerenwinkel, H. B. Schiller, R. Fässler, and D. J. Müller. αV-class integrins exert dual roles on α5β1 integrins to strengthen adhesion to fibronectin. *Nat. Commun*, 8(1):14348, 2017. URL https://doi.org/10.1038/ncomms14348.

[55] C. Zhu. Kinetics and mechanics of cell adhesion. *J. Biomech.*, 33(1):23–33, 2000. URL https://doi.org/10.1016/S0021-9290(99)00163-3.
[56] P. Bongrand. Ligand-receptor interactions. *Rep. Prog. Phys.*, 62(6):921-968, 1999. URL https://doi.org/10.1088/0034-4885/62/6/202.

[57] J Weisel. Protein-protein unbinding induced by force: single-molecule studies. *Curr. Opin. Struct. Biol.*, 13(2):227-235, 2003. URL https://doi.org/10.1016/S0959-440X(03)00039-3.

[58] A. Paul and J. Alper. Calculating the force-dependent unbinding rate of biological macromolecular bonds from force-ramp optical trapping assays. *Sci. Rep.*, 12(1):82, 2022. URL https://doi.org/10.1038/s41598-021-03690-1.

[59] D. Kong, B. Ji, and L. Dai. Stability of Adhesion Clusters and Cell Reorientation under Lateral Cyclic Tension. *Biophys. J.*, 95(8):4034-4044, 2008. URL https://doi.org/10.1529/biophysj.108.131342.

[60] L. Irons, M. R. Owen, R. D. O'Dea, and B. S. Brook. Effect of Loading History on Airway Smooth Muscle Cell-Matrix Adhesions. *Biophys. J.*, 114(11):2679-2690, 2018. URL https://doi.org/10.1016/j.bpj.2018.04.026.

[61] L. Irons, H. Huang, M. R. Owen, R. D. O'Dea, G. A. Meininger, and B. S. Brook. Switching behaviour in vascular smooth muscle cell-matrix adhesion during oscillatory loading. *J. Theor. Biol.*, 502:110387, 2020. URL https://doi.org/10.1016/j.jtbi.2020.110387.

[62] R. Changede and M. Sheetz. Integrin and cadherin clusters: A robust way to organize adhesions for cell mechanics. *BioEssays*, 39(1):1-12, 2017. URL https://doi.org/10.1002/bies.201600123.

[63] R Giavazzi, M Foppolo, R Dossi, and A Remuzzi. Rolling and adhesion of human tumor cells on vascular endothelium under physiological flow conditions. *J. Clin. Investig.*, 92(6): 3038-3044, 1993. URL http://doi.org/10.1172/JCI116928.

[64] Kenneth V. Honn and Dean G. Tang. Adhesion molecules and tumor cell interaction with endothelium and subendothelial matrix. *Cancer Metastasis Rev.*, 11(3-4):353-375, 1992. URL http://doi.org/10.1007/BF01307187.

[65] H C DeGrendele, P Estess, L J Picker, and M H Siegelman. CD44 and its ligand hyaluronate mediate rolling under physiologic flow: a novel lymphocyte-endothelial cell primary adhesion pathway. *J. Exp. Med.*, 183(3):1119-1130, 1996. URL https://doi.org/10.1084/jem.183.3.1119.

[66] Y. Katayama, A Hidalgo, J Chang, A Peired, and P S. Frenette. CD44 is a physiological E-selectin ligand on neutrophils. *J. Exp. Med.*, 201(8):1183-1189, 2005. URL https://doi.org/10.1084/jem.20042014.

[67] S. N. Thomas, F. Zhu, R. L. Schnaar, C. S. Alves, and K. Konstantopoulos. Carcinoembryonic Antigen and CD44 Variant Isoforms Cooperate to Mediate Colon Carcinoma Cell Adhesion to E- and L-selectin in Shear Flow. *J. Biol. Chem.*, 283(23):15647-15655, 2008. URL https://doi.org/10.1074/jbc.M800543200.

[68] R. P. McEver and C. Zhu. Rolling Cell Adhesion. *Annu. Rev. Cell Dev. Biol.*, 26(1):363-396, 2010. URL https://doi.org/10.1146/annurev.cellbio.042308.113238.

[69] C. Christophis, I. Taubert, G. R. Meseck, M. Schubert, M. Grunze, A. D. Ho, and
A. Rosenhahn. Shear Stress Regulates Adhesion and Rolling of CD44+ Leukemic and Hematopoietic Progenitor Cells on Hyaluronan. *Biophys. J.*, 101(3):585-593, 2011. URL https://doi.org/10.1016/j.bpj.2011.05.045.

[70] E Puklin-Faucher and M P. Sheetz. The mechanical integrin cycle. *J. Cell Sci.*, 122(2): 179-186, 2009. URL https://doi.org/10.1242/jcs.042127.

[71] P Roca-Cusachs, T Iskratsch, and M P. Sheetz. Finding the weakest link – exploring integrin-mediated mechanical molecular pathways. *J. Cell Sci.*, 2012. URL https://doi.org/10.1242/jcs.095794.

[72] J A Huaman. *Molecular Mechanisms Underlying Migratory and Immune-Evasion Properties of Circulating Tumor Cells*. PhD thesis, City University of New York, 2020.

[73] S. Rakshit and S. Sivasankar. Biomechanics of cell adhesion: how force regulates the lifetime of adhesive bonds at the single molecule level. *Phys. Chem. Chem. Phys.*, 16(6):2211, 2014. URL https://doi.org/10.1039/c3cp53963f.

[74] G. A. Truskey and J. S. Pirone. The effect of fluid shear stress upon cell adhesion to fibronectin-treated surfaces. *J. Biomed. Mater. Res.*, 24(10):1333-1353, 1990. URL https://doi.org/10.1002/jbm.820241006.

[75] A. Yamamoto, S. Mishima, N. Maruyama, and M. Sumita. Quantitative evaluation of cell attachment to glass, polystyrene, and fibronectin- or collagen-coated polystyrene by measurement of cell adhesive shear force and cell detachment energy. *J. Biomed. Mater. Res.*, 50(2):114-124, 2000. URL https://doi.org/10.1002/(SICI)1097-4636(200005)50:2%3C114::AID-JBM4%3E3.0.CO;2-6.

[76] C. Dong and X. X Lei. Biomechanics of cell rolling: shear flow, cell-surface adhesion, and cell deformability. *J. Biomech.*, 33(1):35-43, 2000. URL https://doi.org/10.1016/S0021-9290(99)00174-8.

[77] S. I. Simon and H. L. Goldsmith. Leukocyte Adhesion Dynamics in Shear Flow. *Ann. Biomed. Eng.*, 30(3):315-332, 2002. URL http://doi.org/10.1114/1.1467677.

[78] S. Mao, Q. Zhang, H. Li, W. Zhang, Q. Huang, M. Khan, and J.-M. Lin. Adhesion analysis of single circulating tumor cells on a base layer of endothelial cells using open microfluidics. *Chemical Science*, 9(39):7694-7699, 2018. URL https://doi.org/10.1039/C8SC03027H.

[79] J. W. Song, S. P. Cavnar, A. C. Walker, K. E. Luker, M. Gupta, Y.-C. Tung, G. D. Luker, and S. Takayama. Microfluidic Endothelium for Studying the Intravascular Adhesion of Metastatic Breast Cancer Cells. *PLoS ONE*, 4(6):e5756, 2009. URL https://doi.org/10.1371/journal.pone.0005756.

[80] R. Riahi, Y. L. Yang, H. Kim, L. Jiang, P. K. Wong, and Y. Zohar. A microfluidic model for organ-specific extravasation of circulating tumor cells. *Biomicrofluidics*, 8(2):024103, 2014. URL https://doi.org/10.1063/1.4868301.

[81] P. Rupprecht, L. Golé, J.-P. Rieu, C. Vézy, R. Ferrigno, H. C. Mertani, and C. Rivière. A tapered channel microfluidic device for comprehensive cell adhesion analysis, using measurements of detachment kinetics and shear stress-dependent motion. *Biomicrofluidics*, 6 (1):014107, 2012. URL https://doi.org/10.1063/1.3673802.
[82] Vincent Hyenne, Jacky G. Goetz, and Naël Osmani. Liquid Biopsies: Flowing Biomarkers. *Adv. Exp. Med. Biol.*, 1379:341-368, 2022. URL https://doi.org/10.1007/978-3-031-04039-9_14.

[83] F. Hecht. New development in freefem++. *J. Numer. Math.*, 20(3-4):251-266, 2012. URL https://doi.org/10.1515/jnum-2012-0013.
# Supplementary Materials of « Deciphering circulating tumor cells binding in a microfluidic system thanks to a parameterized mathematical model »

Giorgia Ciavolella<fnref n="1" /> Julien Granet<fnref n="1" /> Jacky G. Goetz<fnref n="2" /> Naël Osmani<fnref n="2" /> Christèle Etchegaray<fnref n="1" /> Annabelle Collin<fnref n="3" />

October 23, 2024

## A Numerical modeling of the fluid

### A.1 Validation of the hydrodynamic entrance length $\ell$

Let us define the domain $\Omega := (0, L) \times (0, h) \times (0, l)$, $\partial\Omega$ its boundary and $\partial\Omega_{in} := \{0\} \times (0, h) \times (0, l)$ the boundary corresponding to the channel inlet. We would like to numerically establish the distance $\ell$ (along the x-axis) at which a uniform velocity profile from $\partial\Omega_{in}$ fully develops into a Poiseuille profile. To do so, we compare the solution of the two following problems: (i) a Poiseuille profile is already developed over the entire domain $\Omega$, (ii) boundary layer equations are used to take into account the entrance effect near the inlet. Let $u_P$ be the fluid velocity anywhere in $\Omega$ according to the Poiseuille equations

$$ \begin{cases} \rho \frac{\partial u_P}{\partial t} - \mu \Delta u_P = G(1 + \xi_f \cos(\omega_f t + \varphi)) & t > 0, \Omega, \\ u_P(t, \cdot) = 0 & t > 0, \partial\Omega, \\ u_P(0, \cdot) = 0 & \Omega, \end{cases} \quad (1) $$

and $u_B$ be the fluid velocity anywhere in $\Omega$ according to the boundary layer equations

$$ \begin{cases} \rho \left( \frac{\partial u_B}{\partial t} + u_B \nabla \cdot u_B \right) - \mu \Delta u_B = G(1 + \xi_f \cos(\omega_f t + \varphi)) & t > 0, \Omega, \\ u_B(t, \cdot) = 0 & t > 0, \partial\Omega \setminus \partial\Omega_{in}, \\ u_B(t, \cdot) = \max_{t>0} u_P & t > 0, \partial\Omega_{in}, \\ u_B(0, \cdot) = 0 & \Omega. \end{cases} \quad (2) $$

<footnote n="1">¹Inria, Univ. Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251, F-33400 Talence, France</footnote>

<footnote n="2">²INSERM UMR_S 1109, Univ. Strasbourg, FMTS, Équipe labellisée Ligue Contre le Cancer, F-67000 Strasbourg, France.</footnote>

<footnote n="3">³Laboratoire de Mathématiques Jean Leray, Nantes Université, F-44100 Nantes, France</footnote>

<footnote n="1">¹Inria, Univ. Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251, F-33400 Talence, France</footnote>

<footnote n="2">²INSERM UMR_S 1109, Univ. Strasbourg, FMTS, Équipe labellisée Ligue Contre le Cancer, F-67000 Strasbourg, France.</footnote>

<footnote n="3">³Laboratoire de Mathématiques Jean Leray, Nantes Université, F-44100 Nantes, France</footnote>
<figure><img src="image_12.png" /><figcaption>Figure 1: Top: Experimental setup. Middle: Relative error $\varepsilon_1$ on the cross-section $(0, L) \times (0, h) \times \{\frac{L}{2}\}$. For $x \ge 5.12 \times 10^{-4}$ m we have, $\varepsilon_1(\cdot) \le 1\%$. For reference, the entry length $l = 1.53 \times 10^{-3}$ m given by [47] is also displayed. Bottom-Left: Relative error $\varepsilon_2$ on the cross-section $\{\frac{L}{2}\} \times (0, h) \times (0, l)$. For $z \in [6.73 \times 10^{-4}$ m, $3.13 \times 10^{-3}$ m] we have, $\varepsilon_2(\cdot) \le 1\%$. Bottom-Right: Relative error $\varepsilon_2$ on the cross-section $\{\frac{L}{2}\} \times \{h_f^m\} \times (0, l)$.</figcaption></figure>
The third equation of Problem (2) corresponds to a uniform velocity profile at the inlet of the channel, which is here equal to the maximum velocity of the fluid with a fully developed Poiseuille profile. This choice was made to be consistent with the steady state hypothesis. Problems 1 and 2 are numerically solved using FreeFEM++ [83]. The domain is discretized on a 120 × 10 × 20 mesh, using P1 Lagrange elements. The timestep Δt is taken to be 0.5 s and the total time of the simulation is 15 s.

For the parameters, we work with the worst case scenario, e.g. the highest Reynolds number possible, with $Re = 3.75$. We recall that $L = 1.70 \times 10^{-2}$ m, $h = 4.00 \times 10^{-4}$ m, $l = 3.80 \times 10^{-3}$ m, $\rho = 1.00 \times 10^3$ kg.m$^{-3}$, $\mu = 7.20 \times 10^{-4}$ Pa.s, $G = 201.32$ Pa.m$^{-1}$, $\omega_f = 4.88 \times 10^1$ rad.s$^{-1}$, and $\varphi = 0$ rad, see Subsections 2.2, 3.2 and 4.2 for details. By convention, a flow can be considered fully developed when its velocity profile matches the asymptotic one with an error margin of less than a percent. Assuming $u_B(t, \cdot) \neq 0$, for all $t > 0$, the relative error in percent is computed as follows

$$ \varepsilon_1(\cdot) := \max_{t>0} \left( \frac{|u_B(t, \cdot) - u_P(t, \cdot)|}{|u_B(t, \cdot)|} \right) \times 100, \quad \Omega. $$

Figure 1-Middle shows the relative error $\varepsilon_1$ over the longitudinal cross-section $(0, L) \times (0, h) \times \{\frac{l}{2}\}$. One may remark the hydrodynamic entrance length $l$ – displayed along the cross-section – seemingly over-predicts the distance at which $\varepsilon_1$ drops below one percent. The numerical entrance length at which $\varepsilon_1 \le 1\%$ is given by $l_{num} = 5.12 \times 10^{-4}$ m. Numerical tests including spatial and temporal convergence were performed to validate the value of $l_{num}$.

From Subsection 3.2, we know the data was obtained $\frac{L}{2} = 8.50 \times 10^{-3}$ m away from the inlet, which is a full order of magnitude above the upper bound $l$ of the hydrodynamic entrance length. Under those circumstances, the velocity profile of $u_B$ closely matches the velocity profile of $u_P$ and is considered as a fully developed Poiseuille flow.

**Remark A.1.** The value $l = 1.53 \times 10^{-3}$ m given in Subsection 5.2 and coming from [47] is superior to our numerical evaluation of the hydrodynamic entrance length $l_{num}$. Furthermore, additional hydrodynamic entrance length comparisons were performed with other formulas available in the literature (see [47], Table 1), and all of them were found to overpredict $l_{num}$. This difference is explained by the relation between the velocity of the fluid far from the inlet and its velocity at the inlet boundary. As shown in [47], their study was performed for a ratio whose maximum is between 1.5 and 2, while in our study we set this maximum at 1. Therefore, their entry velocity is much higher than the asymptotic velocity of the fluid, which is reflected in the reported value of $l$.

## A.2 Validation of the fluid expression

In Subsection A.1, we have shown the device was long enough for a Poiseuille flow to become fully developed and therefore independent on the x-axis. While this reduces Problem (1) into a 2D one, we can simplify it even further by showing the channel is wide enough that we can neglect a change in the z-axis as well. We recall that aspect ratio $AR = \frac{h}{l}$ equals to $1.05 \times 10^{-1}$, which tells us the microfluidic device is roughly 10 times wider than it is tall. To show that 1D reduced model (in y-axis direction) is reasonable, we compare the fluid velocity $u_P$ to the expression $u_f$, which must verify

$$ \left\{ \begin{array}{ll} \rho \frac{\partial u_f}{\partial t} - \mu \frac{\partial^2 u_f}{\partial y^2} = G(1 + \xi_f \cos(\omega_f t + \varphi)) & t > 0, (0, h), \\ u_f(t, \cdot) = 0 & t > 0, (0, h), \\ u_f(0, \cdot) = 0 & (0, h). \end{array} \right. \quad (3) $$
<figure><img src="image_13.png" /><figcaption>Figure 2: Quotient $\bar{C}(y)/\tilde{C}(y)$ for $y \in (0, h)$.</figcaption></figure>

Using the same parameter values and mesh as in Subsection A.1, we solve Problem (3) with FreeFEM++. The same error threshold of one percent is taken, with

$$ \varepsilon_2(\cdot) := \max_{t>0} \left( \frac{|u_P(t, \cdot) - u_f(t, \cdot)|}{|u_P(t, \cdot)|} \right) \times 100, \quad (0, h) \times (0, l). $$

Figure 1-Bottom-Left gives $\varepsilon_2$ on the lateral cross-section $\{\frac{L}{2}\} \times (0, h) \times (0, l)$ corresponding to the section from which the data were obtained. For $z \in [6.46 \times 10^{-4} \text{ m}, 3.15 \times 10^{-3} \text{ m}]$ we have, $\varepsilon_2(\cdot) \le 1\%$. Figure 1-Bottom-Right shows the relative error $\varepsilon_2$ at the camera's focal plane, using the value of $h_f^m$ found in Subsection 4.2. As expected from the observation made in Figure 1-Bottom-Left, the distance between the two inner ticks does not increase nor decrease in any significant way. The sub-1% relative error area thus spans most of the microfluidic device, with $\{\frac{L}{2}\} \times (0, h) \times (6.46 \times 10^{-4}, 3.15 \times 10^{-3})$. Upon closer inspection of the data provided in [14, 42], we know the data was obtained at the middle of the device length wise, but also width wise, as the lateral sides of the device are not visible in the videos. Given the width of the sub-1% area, $u_P$ can therefore be represented by $u_f$.

## A.3 Validation of the parabolic shape for the velocity profile

For a purely oscillating pressure gradient, a parabolic velocity profile can be considered for a Wo up to 1, after which the profile rapidly evolves into a plug-like shape [48]. Counting the number of oscillations in Figure 3, it can be seen that $\omega_f$ depends on the cohort. The highest value is reached at the third cohort and is close to $50 \text{ s}^{-1}$ (~12 oscillations on 1.5 s gives ~12/1.5 $\times 2\pi$ ~50 s$^{-1}$), resulting in a Wo around 1.6. However, it is shown in [49] that a parabolic profile can still emerge for moderately low values of Wo ($\approx 2$) when

$$ |\text{Re}(u_f(t, y))| \gg |\text{Im}(u_f(t, y))| \quad t > 0, \quad y \in (0, h). $$

Using the superposition principle, one can rewrite $u_f$ as a sum of its constant and oscillating part, that we respectively denote by $\bar{u}_f(t, y)$ and $\tilde{u}_f(t, y)$. We now have to check that

$$ |\text{Re}(\tilde{u}_f(t, y))| + |\text{Re}(\bar{u}_f(t, y))| \gg |\text{Im}(\tilde{u}_f(t, y))| + |\text{Im}(\bar{u}_f(t, y))| \quad t > 0, \quad y \in (0, h). \quad (4) $$

From Equation (2), one has $|\text{Re}(\bar{u}_f(t, y))| = \bar{C}(y) > 0$ and $|\text{Im}(\bar{u}_f(t, y))| = 0$, so that the condition rewrites

$$ |\text{Re}(\tilde{u}_f(t, y))| + \bar{C}(y) \gg |\text{Im}(\tilde{u}_f(t, y))| \quad t > 0, \quad y \in (0, h), $$
which is true for all $t$ if $\bar{C}(y)$ is taken large enough. In the most pessimistic case, we have $|\text{Re}(\tilde{u}_f(t, y))| = 0$ and $\max_{t>0} |\text{Im}(\tilde{u}_f(t, y))| = \tilde{C}(y) > 0$. Therefore, it is enough to show that

$$ \frac{\bar{C}(y)}{\tilde{C}(y)} \gg 1. $$

Figure 2 shows the value of $\frac{\bar{C}(y)}{\tilde{C}(y)}$ for $y \in (0, h)$. It can be deduced that Condition (4) is always verified. This means that for $u_f = \bar{u}_f + \tilde{u}_f$, the constant component is large enough to ensure that the magnitude of oscillations does not disrupt the parabolic shape of the velocity profile. This therefore validates our hypothesis.