Physical Biomarkers of Disease Progression: On-Chip Monitoring of Changes in Mechanobiology of Colorectal Cancer Cells

Disease can induce changes to subcellular components, altering cell phenotype and leading to measurable bulk-material mechanical properties. The mechanical phenotyping of single cells therefore offers many potential diagnostic applications. Cells are viscoelastic and their response to an applied stress is highly dependent on the magnitude and timescale of the actuation. Microfluidics can be used to measure cell deformability over a wide range of flow conditions, operating two distinct flow regimes (shear and inertial) which can expose subtle mechanical properties arising from subcellular components. Here, we investigate the deformability of three colorectal cancer (CRC) cell lines using a range of flow conditions. These cell lines offer a model for CRC metastatic progression; SW480 derived from primary adenocarcinoma, HT29 from a more advanced primary tumor and SW620 from lymph-node metastasis. HL60 (leukemia cells) were also studied as a model circulatory cell, offering a non-epithelial comparison. We demonstrate that microfluidic induced flow deformation can be used to robustly detect mechanical changes associated with CRC progression. We also show that single-cell multivariate analysis, utilising deformation and relaxation dynamics, offers potential to distinguish these different cell types. These results point to the benefit of multiparameter determination for improving detection and accuracy of disease stage diagnosis.

shows an example deformation event with overlayed images of a cell at various positions in the extensional flow junction, where maximum deformation occurs at the SP. The three CRC cells were deformed in a shear and inertia dominant flow regimes. Images of deformation at the SP are provided in Fig. 1e, highlighting the differences in initial size and deformation between cell lines and regimes. By probing both regimes we show that these cell lines are best differentiated from each other using shear-dominant flow conditions. By tracking the deformation and relaxation of the cells, multiple characteristic parameters were extracted. The results were also compared to HL60 (human leukemia) cells, used as a representative for a non-epithelial cell type and showed a softer phenotype compared to the CRC cell lines (Fig. 1e). The CRC cell lines showed increased softness with metastatic progression and the elastic modulus determined for each cell type gave values comparable to previous works using AFM 37,38 . These results consolidate the liklelyhood that cytoskeletal changes during development of metastatic phenotypes leads to mechanical changes of the cell. Further, single cell analysis was performed on four datasets associated to SW480, HT29, SW620 and HL60 cells to confirm the statistical significance of the different parameters for distinguishing cell types, indicating that multiple parameters are needed to accurately separate each of these cell types. The single cell analysis points the way for future determinations of inherent cell heterogeneity and its role in response to therapy.

Results
High speed imaging of the three colorectal cancer cell lines was used to capture the maximum deformation index DI of the cells at the SP for a range of flow rates, Q, in the inertia and shear dominant deformation regimes. In the shear dominant regime the Reynolds number, Re, was less than 6 for the entire range of flow rates. In the inertia dominant regime, Re > 40 for the entire range of flow rates. The deformation index was normalized by the initial diameter of the cells, A, to account for variation in cell size 26 . Figure 2a shows DI/A as a function of Q, in a shear dominant regime where the suspension buffer viscosity was 33 cP. Here, DI/A increased asymptotically toward a maximum deformation value DI A ( / ) max , which was determined by fitting an exponential function. There is a systematic increase in DI/A for SW620 compared to SW480 indicating that under shear dominant conditions the SW620 are significantly more deformable than the SW480 cells. Interestingly, the HT29 display properties similar Size normalized deformation index (DI/A) for three colorectal cancer cell lines over a range of flow rates (µl/min). DI A / was averaged from datasets from 3 separate experiments combined. (a) The flow regime was shear dominant, the viscosity of the cell suspension buffer was 33 cP. The total number of events measured was: 93 < n < 931 for SW480, 160 < n < 596 for HT29 and 280 < n < 734 for SW620. (b) The flow regime was inertia dominant, the viscosity of the cell suspension buffer was 1 cP. The total number of events measured was: 30 < n < 603 for SW480, 47 < n < 619 for HT29 and 30 < n < 450 for SW620.
www.nature.com/scientificreports www.nature.com/scientificreports/ to those of SW620 at low flow rates, for μ < Q l min 30 / , but undergo less deformation with increasing flow rate and end up displaying properties close to those of the SW480 for higher flow rates, ≥ µ Q l min 40 / . Figure 2b shows DI/A as a function of Q, in a inertia-dominant regime where the suspension buffer viscosity was 1 cP. An abrupt change in behavior can be seen at μ = Q l min 400 / , which was identified as the yield stress of the cells and is associated to cytoskeletal breakdown 32 . For each cell line there is a linear relationship between DI/A and Q for μ < Q l min 300 / . For flow rates μ ≥ Q l min 400 / there is also a linear trend, however the gradient increases for both SW480 and SW620. For μ < Q l min 400 / , SW620 are the most deformable and SW480 are the least deformable. However, for μ ≥ Q l min 400 / the cell lines are less distinguishable from each other. Figure S1 shows the unnormalized data, DI as a function of Q, for the three cell lines in both shear and inertia-dominant regimes.
Averaged deformation traces were measured for each of the 3 cell lines, supplementary videos 1-3 show example deformation events. Figure 3 shows the average strain ε as a function of time for SW480, HT29 and SW620 cells as they deform and recover in the extensional flow device. Cells were deformed in a shear-dominant and low velocity regime (μ μ = . ± . = cP Q l 33 4 0 3 and 5 /min). Multiple parameters were determined from the averaged deformation traces of N = 56 SW480, N = 49 HT29 and N = 50 SW620, which are labelled in the traces in Fig. 3 and values are summarised in Table 1. A is the initial cell diameter, ε max in the maximum strain, ε 0 is the initial strain before entering the SP and τ d is the deformation time as the cell deforms. After cells reached their maximum strain ε max near to the SP, the strain began to decrease exponentially as the cells traversed the outlet channel. A cell relaxation time, τ r , was extracted from an exponential fit, as well as a final strain ε ∞ from extrapolation of the fit. The strain values are vectors due to the inlet channels being perpendicular to the outlet channels, hence some ε 0 values are slightly below 0 due to small deformations induced by shear channel confinement. Thus, if cells recover their original shape after deformation the magnitude of initial strain should be the same as the magnitude of final strain ( ε ε | | = | | ∞ 0 ). The areas shaded in pink in Fig. 3 are a width of ε | | 2 0 , and denote whether cells recover their initial strain.
The elastic modulus E was extracted using the Kelvin-Voigt model (Eq. 4). Figure S2 shows the velocity profiles calculated along the central axis of the extensional flow junction, that were fitted with a sine function and   Table 1 www.nature.com/scientificreports www.nature.com/scientificreports/ used to fit the Kelvin-Voigt model to the deformation trace (shown in red in Figure S2). As a comparison, control data from the averaged deformation trace of N = 50 HL60 cells is also included in Table 1 32 .
Single cell analysis (SCA) was undertaken on the datasets extracted from tracking ε as a function of time for each of the cell lines. Two sample t-tests were performed on 5 of the parameters to quantify the statistical significance between the different cell lines when using these parameters 40 . Figure 4 shows a bar graph of the average A, ε max , ε ∞ , τ r of the 4 cell lines and indicates the level of significance: where p > 0.05 is not significant (ns), 0.01 < p < 0.05 is significant (*), 0.001 < p < 0.01 is very significant (**), 0.0001 < p < 0.001 is extremely significant (***) and p < 0.0001 is extremely significant (****), raw values can be found in Table S1. The data for ε 0 is shown in Figure S6.
Linear Discriminant analysis (LDA) was also used on the four datasets to test the abilities of these parameters for classifying the different cell types (Fig. 5). LDA is a supervised multivariate method that obtains a linear combination of the parameters that best separates the different cell lines. When trained on a 4-class dataset (4 cell types), where each linear discriminant (LD) maximizes the separation of a pair of classes, using all the LDs scores for the final classification. A k-fold validation test was applied to the data, where a random fraction of the data was used to train the LDA model. This model was then applied to the remaining data to test the models ability to correctly classify the cell types based on the given parameters. Here, a 5-fold validation test was applied, where the loadings and scores of the LDs are shown in Fig. 5 Table 2. LD1, LD3 and LD5 showed the best separation between HL60 and the three CRC cell lines, with individual scores shown in Figure 5bi, iii and v. These discriminants were all marked by a higher maximum strain ε max of the non-adherent cells compared to the adherent cells, giving good classification of HL60 compared to the CRC cell lines. The final strain ε ∞ discriminated HT29 and SW620 from HL60, and initial diameter showed good separation of HL60 from HT29 and SW480.

and the confusion matrix is shown in
Between the CRC cell lines, LD2 and LD6 showed the best separation between SW480 cells with HT29 and SW620 cells respectively. LD2 showed higher relaxation time τ r of SW620 and HT29 cells when compared to SW480 cells, accompanied by higher final strain ε ∞ and maximum strain ε max but lower initial strain ε 0 and diameter A. LD6 and LD4, comparing SW620 cells with either SW480 or HT29 cells, mainly classified according to diameter A, showing that SW620 cells were smaller. When using 5-fold LDA classification, 82% of the HL60, 71% www.nature.com/scientificreports www.nature.com/scientificreports/ of SW480 and 85% of SW620 cells were correctly classified by the model. HT29 cells were more difficult to accurately classify, 36% correctly classified with 39% incorrectly classified at SW480 and 20% as SW620. Of the 4-class dataset, the average classification rate was ~69%.

Discussion
Three colorectal cancer cell lines representing different stages of disease progression from Dukes primary stage B (SW480) through to metastasis to the local lymphatic system (SW620) Dukes stage C were studied in the inertial and shear regimes. HT29s represent an intermediate state of Dukes stage C, primary tumour. In the inertial regime, (Fig. 2b) for Q < 300 μl/min the size normalized deformation, DI/A, was greatest for the SW620 and least for the SW480. For Q > 300 μl/min, the cell lines become indistinguishable by DI/A suggesting that Q ≅ 400 μl/min represents the cell yield stress. In Armistead et al. (2019) we showed that cell softening due to disruption of the cytoskeleton using LatA was only measurable at flow rates below the yield stress of HL60 cells. This corroborates that cytoskeletal changes are associated with CRC progression, and cytoskeletal breakdown above the yield stress results in no measurable changes to deformability for the higher flow rates.
In the shear regime the SW620 cells were the most deformable for all values of Q and the SW480's the least (Fig. 2a). The HT29 cells tended to be similar to those of the SW620s at low flow rates but their DI/A values approached those of the less deformable SW480's for higher flow rates, ≥ 40 μl/min. The secondary tumor cells being softer than the primary cells is a result seen in previous works 37,38 . Several papers report that metastatic SW620 cells up-regulate genes associated with cytoskeletal changes, particularly related to the actin, which accompany increased motility, enhanced invasive potential, higher proliferation and reduced adhesion compared to the non-metastatic SW480 [41][42][43][44][45] .
We note, that had we characterized the deformation without adjusting for cell size, these cells would not be distinguished from each other based on their DI index alone in either the inertial or shear regimes ( Figure S1). In the inertial regime, the SW480 showed increased DI compared to HT29 and SW620 for Q > 400 µl/min. However, the SW480 are the largest cell (Table 1) and this would make them appear more deformable above the yield stress associated with subcellular disruption of actin filamnets (<400 µl/min). Figure S6 shows the average DI of these cell lines when deformed at 5 µl/min in the shear regime. Here, SW480 had the smallest DI and HT29 the largest
We have previously shown that a low-strain and shear dominant regime is most sensitive to cytoskeletal changes 32 . Thus, to study mechanical changes in our CRC model system we performed deformation tracking and multiparameter analysis of SW480, HT29 and SW620 deforming at Q = 5 µl/min, in the shear regime. The general trait of SW620 being softer and more deformable than SW480 cells was also evident in the parameters determined from the deformation traces (Table 1) with the SW620 having a lower mean elastic modulus than SW480's, and a higher maximum strain ε max . Palmieri et al. 37 noted that SW480 cells have two appearances in culture, an epithelial-type morphology (E-type) and a rounded morphology 37 . Using AFM they found the elastic modulus of SW480 E-type to be 1060 Pa and SW480 R-type to be 580 Pa. The elastic modulus determined by AFM for SW480 R-type cells was within error the same as that determined here using microfluidics. Further, both AFM and microfluidics show the metastatic cells to be softer than the non-metastatic (SW480). Additionally, the elastic modulus determined for the SW620 cells (Table 1) is within error of the value reported using AFM 37 . Tsikritis et al. 2015 also measured the elastic modulus of SW480 and SW620 using AFM, finding SW620 to be ~3 fold softer 38 .
SW480 and HT29 are closest in terms of initial diameter compared to SW620, however ε max and the elastic modulus of HT29 and SW620 are comparable whereas SW480 is significantly higher. These results suggest that as colorectal cancer progresses from Duke stage B to C (SW480 to HT29), requiring the cells to migrate toward the outer lining of the bowel, the cells undergo structural changes and become softer. The similarity between HT29 and SW620 is a possible indicator that as the cells then move from the outer lining to a secondary tumour site in the lymph nodes, changes to the cell structure are less essential. Further, the SW480's also recovered their original strain values after deformation, whereas SW620 and HT29 both recovered to a final strain of ε > . ∞ 0 04 which is significantly higher than their initial strain ~0. This could indicate an additional slower relaxation process occurring over a longer timescale than our experiments measure, or that the cells have an associated "plastic" deformation. Cell deformation is commonly thought to be a viscoelastic response, however permanent plastic deformations have been seen due to breakdown of the cytoskeletal scaffold 46 . The short relaxation time for the SW620 and HT29 might reflect a more active cytoskeleton. For comparison we note that the HL60 cells, also fully recover to their strain free values.
Single cell analysis was performed on the datasets acquired for obtaining the deformation traces of N = 50 HL60 cells, N = 56 SW480 cells, N = 49 HT29 cells and N = 50 SW620 cells. The significance of the parameters A, ε max , ε ∞ , τ r and ε 0 were tested using two sample t-tests to obtain p-values, shown in Table S2 and summarised by Fig. 4 and Figure S7. We first note that the initial strain ε 0 of the four cell lines are within error and show no significant difference ( Figure S7), additionally the deformation time τ d showed no change between cell lines and thus was not studied using SCA. Statistical significance is shown in separating all the cell lines from each other using their initial diameter A (Fig. 4a) and their maximum strain (Fig. 4b) These results show that even though SW480 and HT29 cannot be separated by their initial size, they can be identified by their maximum strain. It also corroborates previous discussions that SW620 and HT29 have similar deformability. Fitting the cell recovery with an exponential function was used to extrapolate values of final strain ε ∞ (Fig. 4c) and their relaxation time (Fig. 4d). The final strain ε ∞ between SW480 and the other CRC cell lines is extremely significant (p < 0.0001) whereas HT29 and SW620 have no significant difference in ε ∞ . These results continue to show the trend that SW480 are mechanically different to the later stage CRC cell lines HT29 and SW620, which show indistinguishable deformation and relaxation parameters.
Results show that τ r cannot be used to distinguish between any of the CRC cell lines (all have p > 0.05), but all three CRC cell lines compared to HL60 are extremely significant. HL60 are the most deformable and have a relaxation time ~3 fold larger than the CRC cell lines. The deformation and relaxation parameters (ε ε ∞ , max ) were able to distinguish SW480 from SW620 and HT20, and the data suggests SW480 is stiffer than the later stage CRC cell lines. Mechanical changes in CRC cell lines has been attributed to changes in the actin. This suggests that τ r may be more sensitive to the mechanics of the nucleus as opposed to the cytoskeleton. Table S3 shows the nuclear diameter and nuclear ratio (A nucleus /A cell ) of the four cell lines. The nucleus is known to be mechanically stiffer compared to the rest of the cell 47,48 , and HL60 have a smaller nucleus and nuclear diameter compared to the CRC cell lines. However, the responses are likely more complex due to coupling between τ r and ε ∞ . As the extrapolated value ε ∞ did not recover to the initial value (ε ε ≠ | | ∞ 0 ), this suggests longer relaxation processes are at play which could offer more insight into the mechanical response of cells to an applied stress. Overrall, the statistical testing shows that whilst no single parameter can significantly distinguish all cell lines from each other, multiparameter analysis can lead to more classification.
Finally, a 5-fold validation test using LDA showed that HL60, SW480 and SW620 had moderately good classification rates (<71%) However, only 36% HT29 were classified correctly with the majority incorrectly classified as SW480 (39%) or SW620 (20%). HT29 are the intermediate step of the CRC model and were generally harder to distinguish, this may be indicative of these cells having intermediate properties between the two. For instance, results showed that HT29 have comparable initial size to SW480 but deformation and relaxation parameters comparable to SW620 (Fig. 4). The relatively low sample size is also likely to hinder classification rates (n < 56). Overrall, the average correct classification rate of the four cell types was 69% which highlights the need for larger datasets combined with multiparameter analysis for accurate classification of different cell types.

conclusion
Three CRC cell lines, representing different states of disease progression, were studied in the shear and inertia-dominant flow regimes. In the shear dominant regime, we measured deformation traces and determined multiple characteristic parameters, including: maximum strain ε max , elastic modulus E and relaxation time τ r . Interestingly, the elastic modulus values of each cell line were of the same order of magnitude as previous AFM Scientific RepoRtS | (2020) 10:3254 | https://doi.org/10.1038/s41598-020-59952-x www.nature.com/scientificreports www.nature.com/scientificreports/ measurements, in spite of the different modes and timescales of operation. These results showed a general increase in deformability with CRC progression using ε max and elastic modulus. The cell types could only be differentiated under certain flow conditions, indicating that in the inertial regime only low-strains (below the yield stress) are sensitive to cytoskeletal changes whereas high-strains (above the yield stress) are not. Additionally, at a critical flow rate in the inertia-dominant regime we saw an increase in gradient of DI as a function of flow rate, this occurs due to the breakdown of the cells internal structure, and as such the mechanical properties probed beyond this point would depend only on the viscous properties of the cytoplasm. This corrobated previous works that show that changes to actin structure associated with metastatic progression result in the cells becoming softer.
Our results are the first example of using microfluidic deformation to distinguish between non-metastatic and metastatic CRCs and support the expectation that metastatic cells are more deformable than non-metastatic cells due to cytoskeletal changes. Further, we found that multiple parameters were needed to be able to distinguish the four cell types from each other. We observed that SW620 and HT29 are more deformable and softer than SW480 and do not recover their original strain suggesting that they undergo an additional slower relaxation process, occurring over a time period too long to be captured in our experiments. Single cell and multiple parameter analysis showed changes in the mechanical properties of CRC cells attributed to specific sub-structural changes, and that these can also be used to classify different cell types. Results show that a single-cell high-throughput technique combined with multiparameter analysis will allow us to advance our understanding of cancer progression, and to accurately classify heterogeneous samples of disease states.

Materials and Methods
Microfluidic devices. Microdevices were fabricated using standard rapid prototyping and soft lithography techniques, using a Silcon master as a mold for fabricating microfluidic devices in polydimethylsiloxane (PDMS). Piranha wet etch (H 2 SO 4 & H 2 O 2 ) was used to clean a 3-inch Silicon wafer, which was then rinsed with deionized water. A 25 µm layer of photoresist SU-8 2025 (Microchem, Warickshire, UK) was coated onto the wafer.
Channels were etched into the SU-8 using direct-write laser lithography (MicroWriter ML ™ , Durham Magneto Optics), using a 375 nm wavelength laser. PDMS base and cross-linking agent (SYLGARD 184) were mixed at a 1:10 ratio and poured onto the Silicon master, resulting in a negative replica of the SU-8 structures. This was cured for 1 hr at 75 o C, forming a hydrophobic elastomer which was peeled away from the master. Device inlet and outlet holes were punched using a Biopsy Puncher, then the PDMS was sealed to a glass slide by treatment with Oxygen plasma. Microdevices had a channel height of 25 µm and inlet and outlet channels at the extensional flow junction (Fig. 1a) had a width of 35 µm.
Characterising flow regime. Cells were deformed in a microfluidic device with a cross-slot geometry, cells passed through the stagnation point SP of an extensional-flow junction where their subsequent maximum deformation was measured using the deformation index DI = H/W (Fig. 1b). A compressive force F C and a shear force F S act on the cell, where the total force on the cell at the SP can be defined as = + F F F T S C . Equation (1) was used to determine F C , where ρ is the density of the suspension media, U is the fluid velocity, A p is the cross sectional area of the cell and C D is the drag coefficient 49,50 . Equation (2) was used to determine F S , where where μ is the viscosity of the suspension media, r is the cell radius and γ is the strain rate 18,27 . By tailoring the suspension medium viscosity and the flow rate, the flow regime can be either shear or inertia dominant. This can be described using the Reynolds number Re, where Re > 40 represents the change from a shear to an inertial regime 22,42,43 . A more detailed description and calculations of the flow regimes can be found in ref. 32 .
S 2 experimental procedure. To capture cell deformation on-chip, microfluidic devices were positioned on an inverted bright-field microscope (Eclipse Ti-U, Nikon, Japan), using a 10x objective (with an additional 1.5x magnification for flow rates < µ Q l min 100 / ). High speed microscopy (Photron, Tokyo, Japan) was used at frame rates of 7500-260,000 fps and exposure times of 0.37-6.67 µs. To capture images at high frame rates and low exposure times, an additional light source was mounted above the set-up. Off-line automated image analysis was performed using ImageJ and Matlab, which tracked the shape and position of each cell as a function of time. Utilising a mathematical image processing algorithm adapted from flagellar image tracking 51 . This allowed parameters such as initial diameter, velocity of deformation and other cell shape parameters (Fig. 1b) to be extracted, including the time evolution of the shape dynamics.
Cell events were selected which underwent the same applied stress whilst passing through the extensional flow junction, allowing accurate comparison of deformability between samples. This was achieved by looking at the velocity profile of each cell traversing the extensional flow junction. Exclusion of events was based on the change in velocity, ∆v, (Eq. 3) where v inlet is cell velocity in the inlet channel, and v min is the minimum cell velocity in the cross flow junction. If a cell perfectly deforms whilst trapped at the SP inlet m in inlet calculation of cell elastic modulus. The Kelvin-Voigt model, comprising of a linear spring and a viscous dashpot in parallel, was fitted to the strain as a function of time to determine the elastic modulus of the four cell types. The velocity profile of the cell as it passes the extensional flow junction varies approximately as a sine function, for a period T, where ω = π T 2 is the angular frequency. In the shear-dominant regime, the deformation velocity is proportional to applied stress σ t ( ) (Eq. 1) and thus the stress also varies as a sine function. This leads to the adapted trigonometric analytical solution of the Kelvin Voigt model, Eq. 4. Where ε t ( ) is the strain, E is the elastic modulus associated with the linear spring, η is the viscosity associated with the dashpot, σ 0 is the maximum stress on the cell at the SP 32 .
Et 0 2 2 2 2 2 2 2 2 2 2 cell culture. The SW480, SW620 and HT29 cell lines were provided by St James's University Hospital, and had been STR (short tandem repeat) profiled to authenticate cell identity. Cells were cultured in Dulbecco's Modified Eagle Medium (DMEM F-12, Glibco) supplemented with 10% Fetal Bovine Serum (Sigma), 2 mM Glutamax (Thermo Fisher Scientific) and Penicillin 100 units/ml Streptomycin 100 µg/ml (Sigma). Passage numbers were below 50 for all experiments. All three are adherent cell lines and were detached by incubating in TrypLE (Thermo Fisher Scientific) for 5 mins. They were then centrifuged at 100 g for 4 mins and then resuspended in DMEM or PBS with 0.5% (w/v) methyl cellulose (Sigma, UK). The HL60 cell line was purchased as a frozen stock (ECACC European Collection of Authenticated Cell Cultures, 98070106) and cultured in RPMI (Roswell Park Memorial Institute) growth media supplemented with 10% Fetal Bovine Serum (Sigma), 2 mM Glutamax (Thermo Fisher Scientific) and Penicillin 100 units/mL Streptomycin 100 µg/mL (Sigma). HL60 are a non-adherent cell line, centrifuging at 100 g for 4 mins was sufficient to visibly pellet the cells which were then gently resuspended in the desired suspension medium. Cells were either suspended in RPMI or re-suspended in PBS with or 0.50% (w/v) methyl cellulose (Sigma, UK) to increase viscosity. The viscosity of the cell suspension buffer (PBS with 0.5% methyl cellulose) was measured using a Rheometrics SR-500 Dynamic Stress Rheometer in the parallel plate configuration (diameter of 25 mm).