Organs on chip approach: a tool to evaluate cancer -immune cells interactions

In this paper we discuss the applicability of numerical descriptors and statistical physics concepts to characterize complex biological systems observed at microscopic level through organ on chip approach. To this end, we employ data collected on a microfluidic platform in which leukocytes can move through suitably built channels toward their target. Leukocyte behavior is recorded by standard time lapse imaging. In particular, we analyze three groups of human peripheral blood mononuclear cells (PBMC): heterozygous mutants (in which only one copy of the FPR1 gene is normal), homozygous mutants (in which both alleles encoding FPR1 are loss-of-function variants) and cells from ‘wild type’ donors (with normal expression of FPR1). We characterize the migration of these cells providing a quantitative confirmation of the essential role of FPR1 in cancer chemotherapy response. Indeed wild type PBMC perform biased random walks toward chemotherapy-treated cancer cells establishing persistent interactions with them. Conversely, heterozygous mutants present a weaker bias in their motion and homozygous mutants perform rather uncorrelated random walks, both failing to engage with their targets. We next focus on wild type cells and study the interactions of leukocytes with cancerous cells developing a novel heuristic procedure, inspired by Lyapunov stability in dynamical systems.

n being the number of steps performed up to the time considered and t = τn.
The distribution P(r − r ′ ) qualitatively controls the resulting random walk, possibly giving rise to deterministic walks (e.g. P(r − r ′ ) = δ r−r ′ ,k , k = 0, corresponding to a ballistic motion) or to completely stochastic walks (e.g. P(r − r ′ ) = δ |r−r ′ |,r , corresponding to an isotropic motion where steps have fixed lengthr), etc 1 .
As suggested by Eq. S3, the knowledge of the specific distribution P(∆x, ∆y) possibly allows to get an explicit expression for p((x, y), n).
For instance, one can show that, when diffusion is isotropic, i.e. it is equal in the x and y directions, any distribution P(∆x, ∆y) fulfilling the central limit theorem asymptotically (resuming the continuous time description) leads to the wellknown diffusive limit characterized by the normal distribution  Figure S1.
Step length distributions of CA cells. We repeat the analysis of Fig.1 focusing on CA cells. Experimental data with standard errors, fitting functions and the best-fit coefficients are reported in each panel. The distribution is exponential along negative x and Gaussian along positive x, as in the CC case. On the other hand, for the motion along the y direction, step lengths are distributed exponentially regardless of their sign. This suggests that the drift is weaker than the CC case and it acts mainly along the x direction. Step length distributions of AA cells. We repeat the analysis of Fig.1 focusing on AA cells. Experimental data with standard errors, fitting functions and the best-fit coefficients are reported in each panel. Along both x and y directions the distributions are exponential indicating that the motion displays no significant drift and the cells are expected to perform a random motion with no bias.

3/6 2 Application of Lyapunov stability criteria
We present a detailed mathematical explanation of our application of Lyapunov stability criteria used to analyze PBMCs behavior.
As anticipated in the main text (see Sec. 2.1), we associate to each PBMC an equation of motion in the formẋ = F(x(t)), where x(t) is the vectorial position of the cell as a function of time and F is the driving force responsible for the motion. Then, we consider each TC position as a fixed point x 0 of the equation of motion (i.e., the points where F(x 0 ) = 0, namely where the velocity is zero), appointing them as candidate attractors. A key remark to support this point is that, as expected from stability theory, the velocities of CC cells sensibly reduce when they approach the targeted TC (see Fig. 5). However, we need to relax the deterministic requirement F(x 0 ) = 0 toward F(x 0 ) ≈ 0 due to the presence of noise. This means that, for each PBMC trajectory, we effectively consider as a fixed point the one that is closest to the TC position. Then, in order to characterize this as a fixed point, we need to estimate the four entries of the related Jacobian matrix A. To this aim we can not rely on analytical methods because, as underlined in the main text, we do not have the explicit expression for F(x). On the other hand, we have experimental data forẋ versus x and we can directly fit them, as long as the linear condition is verified at least in the spatial region close to x 0 .
Following this route we study in details each trajectory recorded for CC cells approaching a TC, namely the evolution of the position x traced by each single CC cell in the neighborhood of x 0 . This analysis is accomplished in the phase space (ẋ, x), whereẋ is the velocity of the cell, calculated in either Cartesian (data not shown) or polar coordinates, while x is the position of the cell, calculated in either Cartesian (data not shown) or polar coordinates. Assuming a smooth motion close to TC, we can effectively quantify CC's behavior close to the TC (i.e. for t = 1, 2, ...) by fitting the data points ofẋ versus x seeking for a linear relation that mathematically can be represented as v y =ẏ = A y,x x + A y,y y, for Cartesian coordinates, or as v r =ṙ = A r,r r + A r,θ θ , in polar coordinates. Focusing on the latter case, we can recast the related equations compactly exploiting a matricial notation aṡ

4/6 3 Calculus of the Jacobian matrix elements
In the heuristic approach proposed and applied in the main text (see Sec. 2.1 and 2.2) we identify a TC with a candidate attractor, that is, we identify the TC position with x 0 , we estimate the four elements of the related Jacobian matrix A fitting the experimental data ofẋ versus x in the linear regime close to x 0 for each PBMC approaching TC and we compute the matrix eigenvalues in order to infer the stability of the convergence of the PBMC to the TC. In polar coordinates (r, θ ) this means that we fit the experimental data of (v r , v θ ) versus (r, θ ) in the linear regime close to (r 0 , θ 0 ) for each PBMC approaching TC.
Here we present the fitting analysis of the velocity data for the stable and metastable case shown in Sec. 2.2 from which we can extract the elements of the Jacobian matrix and the related Lyapunov coefficients.
In each case we plot the polar components of the velocity v r and v θ of the immune cell versus the polar coordinates r and θ within the linearization range.
In the stable case, the velocities show a linear behavior up to a distance of ∼ 20µm from the TC (see Fig. S3). As pointed out in the main text (see Sec. 2.2), the radial velocity decreases linearly with r and θ , while the angular velocity increases linearly for decreasing r → 0 and θ → θ 0 . The fitting functions in this linearization range are y = p 1 x + p 2 , where y = (v r , v θ ), x = (r, θ ) and p 1 are the elements of the Jacobian matrix. From the diagonalization of this matrix we obtain two negative eigenvalues (λ 1 , λ 2 ) = (−0.02 ± 0.01, −4.45 ± 0.01) that confirm the stable convergence of the trajectory to the fixed point (r 0 , θ 0 ) (see blue dot in Fig. 8, left panel).
In the metastable case, the velocities show a linear behavior up to a distance of ∼ 49µm from the TC (see Fig. S4). As shown in the main text (see Sec. 2.2), radial and angular velocities decrease and increase linearly with r respectively, while they present a parabolic behavior respect to θ in proximity of the fixed point chosen (i.e. near the TC), due to the fact that the PBMC first approaches and, then, leaves the tumor cell. In the first two cases (Fig. S4, upper panels), the fitting functions are y = p 1 x + p 2 , where y = (v r , v θ ), x = (r, θ ) and p 1 are two elements of the Jacobian matrix. In the case of the parabolic behavior, we linearize the fitting function y = p 1 x 2 + p 2 x + p 3 around the fixed point x 0 where F(x 0 ) ≈ 0, obtaining the other two elements of the Jacobian matrix, whose diagonalization gives the positive and negative eigenvalues (λ 1 , λ 2 ) = (−0.02 ± 0.01, 1.04 ± 0.01) that confirm the metastable convergence of the trajectory to the fixed point (r 0 , θ 0 ) (see yellow dot in Fig. 8, left panel).   Figure S4. Analysis of the metastable case. The polar components of the velocity v r and v θ of the immune cell showing metastable convergence (see Fig. 7, panels B1 − B3) are plotted versus the polar coordinates r and θ within the linearization range. The fitting functions are y = p 1 x + p 2 (linear fit) and y = p 1 x 2 + p 2 x + p 3 (quadratic fit) and the best-fit coefficients are reported for each panel, where p 1 of the linear fit and 2p 1 x 0 + p 2 (with F(x 0 ) ≈ 0) obtained from the linearized quadratic fit are the elements of the Jacobian matrix. The resulting R 2 are, respectively, R 2 v r,r = 0.996, R 2 v θ ,r = 0.995, R 2 v r,θ = 0.972, R 2 v θ ,θ = 0.969 and confirm the locally linear behavior of the PBMC trajectory in proximity of the TC (attractor point). Experimental data with standard errors are always grouped in bins with the same number of points in each bin.