The capability to parameterize shapes is of essential importance in biomechanics to identify cells, to track their motion, and to quantify deformation. While various shape descriptors have already been investigated to study the morphology and migration of adherent cells, little is known of how the mathematical definition of a contour impacts the outcome of rheological experiments on cells in suspension. In microfluidic systems, hydrodynamic stress distributions induce time-dependent cell deformation that needs to be quantified to determine viscoelastic properties. Here, we compared nine different shape descriptors to characterize the deformation of suspended cells in an extensional as well as shear flow using dynamic real-time deformability cytometry. While stress relaxation depends on the amplitude and duration of stress, our results demonstrate that steady-state deformation can be predicted from single cell traces even for translocation times shorter than their characteristic time. Implementing an analytical simulation, performing experiments, and testing various data analysis strategies, we compared single cell and ensemble studies to address the question of computational costs vs experimental accuracy. Results indicate that high-throughput viscoelastic measurements of cells in suspension can be performed on an ensemble scale as long as the characteristic time matches the dimensions of the microfluidic system. Finally, we introduced a score to evaluate the shape descriptor-dependent effect size for cell deformation after cytoskeletal modifications. We provide evidence that single cell analysis in an extensional flow provides the highest sensitivity independent of shape parametrization, while inverse Haralick's circularity is mostly applicable to study cells in shear flow.
Mechanical properties of biological systems, e.g., cells and tissues, provide a direct link to their function and can be seen as a biomarker1–5 in cell development, cell state determination, and disease screening.2,4–7 As rheological tests are based on the application of external stresses while monitoring strain, a correct definition of cell deformation is of key importance.5,8,9 Various shape descriptors enable parametric analyses based on geometrical features like size, perimeter, circularity, and membrane curvature. Their potential has already been demonstrated in establishing a link between the cell function and shape for the migration of adherent cells and to investigate the potential of morphological features for a quantitative analysis of cell microscopy images in general.8,9 In contrast, a systematic analysis of shape descriptors for cells in suspension is currently elusive. The lack of standardization makes it difficult to compare viscoelastic properties quantified by different microfluidic methods.
Earlier work investigated red blood cell (RBC) shapes and dynamics in microchannels of different size and at different flow rates.10,11 Experiments have been compared to simulations assuming the RBC membrane as a triangulated network of springs where fluid flow has been modeled by dissipative particle dynamics10 or an immersed boundary algorithm was used, solving the evolution of the moving deformable membrane of RBCs overlayed on a static Eulerian domain.11 While simulations predict well-defined RBC states at given hydrodynamic conditions, experimental data show a broad distribution of shapes.10
RBCs deformed into the shape of parachutes have been analyzed for their shear modulus, determined from a confinement parameter, the Taylor deformation, of a cell-enclosing ellipse.11 The shear modulus of RBCs was shown to increase with aging. The role of deformation and shape of RBCs has also been investigated for patients suffering from diabetes with normal and high levels of cholesterol where Babu reports that deformation is reduced in individuals with high cholesterol.12 Membrane bending of RBCs is also highly relevant for the investigation of malaria disease pathogenesis. Koch et al. demonstrated that the malaria parasite locally reduces the bending modulus to facilitate invasion.13 RBCs under pathological conditions have been compared within different microchannel geometries based on the shape parameters like Taylor deformation and axes ratio of a bounding box, where hyperbolic-shaped constrictions yielded the highest sensitivity to small changes in deformability due to the strong extensional flow.14
High-throughput assays to quantify the mechanical properties of suspended cells have long been limited to RBCs and other soft particles like hydrogel beads since hydrodynamic shear and normal stress were insufficient to induce a finite deformation in nucleated cells. The advent of technologies like deformability cytometry3 and real-time deformability cytometry15–17 expanded the applicability of microfluidic systems to stiffer objects like white blood cells by exploiting the possibilities of inertia-driven flows or by increasing the viscosity of the measurement buffer. These approaches enabled mechanical phenotyping of nearly all cell types and has been used in basic and translational research, e.g., to identify and characterize hematopoietic and mesenchymal stem cells,15,18–20 to monitor the onset and progression of diseases21 and as a quality control in translational medicine.22 While these studies and others contributed a lot to our understanding of life sciences from a biophysical perspective, a systematic description of how shape quantification impacts on deformation analysis is still lacking.
In our work, we compared nine different shape descriptors to characterize the deformation of suspended cells in an extensional as well as shear flow, which originates from microfluidic geometries found in almost any lab-on-chip design. A detailed understanding of how hydrodynamics and shape parametrization impacts on stress distribution and strain quantification, respectively, is of utmost importance for the interpretation of data from any rheological assay that investigates cells in motion. We previously had shed light on the role of hydrodynamic stress for the shape evolution of cells17 and introduced a framework for statistical data evaluation.23 Now we expand our analyses toward the comparison of nine parameters describing cell shapes that are being used in stress relaxation or creep-compliance experiments.3,5,15,24–27 Those descriptors include cell area, front radius, axes ratio, Taylor deformation of a bounding box, principal axes ratio, rescaled circularity, inverse Haralick's circularity, circular variance, and inertia ratio.
Within the present study, we first monitored the evolution of cell shape in response to an extensional as well as shear flow. While stress relaxation dynamics differ for all shape descriptors, our results demonstrate that steady-state deformation can be predicted from single cell traces even for translocation times shorter than its characteristic time. This is an essential requirement to extract material properties from stress–strain experiments. Second, we compared single cell vs ensemble studies to address the question of computational effort vs experimental accuracy. Our data indicate the necessity of analyzing single cell traces for short or unknown relaxation times, i.e., for cells with unknown material properties. Finally, we investigated how cytoskeletal modifications impact on cell deformation represented by different shape descriptors. As a main contributor to mechanical stability, actin polymerization was inhibited by cytochalasin D (CytoD), and the size of this effect has been calculated as a function of flow geometry, i.e., extensional vs shear flow, of analysis strategy, i.e., single cell vs ensemble and of the shape descriptors investigated within this work. Single cell analysis in response to an extensional flow provides highest sensitivity independent of how deformation is quantified while inverse Haralick's circularity provides highest score in effect size for cells in shear flow. This experiment is of fundamental importance in mechanobiology to link alterations on a molecular scale to a cellular phenotype.
MATERIALS AND METHODS
A microfluidic system described in Refs. 15 and 17 is used. Briefly, a narrow microfluidic channel of 30 × 30 μm2 cross section and 300 μm length fabricated facilitating soft-lithography is assembled on a standard inverted microscope as part of the AcCellerator system (Zellmechanik Dresden). While cells in suspension pass the constriction, they are illuminated by a high-power pulsed LED (L1, Zellmechanik Dresden) and imaged using a high-speed CMOS camera (MC1362, Mikrotron) as well as a 40x objective (A-Plan, NA 0.65, Zeiss) leading to a resolution of 0.34 μm per pixel. Sample suspension as well as a sheath flow for hydrodynamic focusing into the channel are connected via tubing to and controlled by a two-channel syringe pump (Nemesys, Cetoni). An outlet is connected to a waste reservoir.
All experiments have been performed using the model system HL60, a myeloid precursor suspension cell line (courtesy of Dan and Ada Olins). Cells are cultured in an RPMI-1640 medium (BioWest) with 10% FCS (Gibco), 1% penicillin/streptomycin (BioWest), and 2 mM L-Glutamin (BioWest) in a standard incubator at 37 °C, 5% CO2, and 95% air. Cells are split approximately every 48 h by centrifugation for 5 min at 200 rcf (Allegra X-15R, Beckman Coulter). The supernatant is discarded and the cells are re-suspended in medium. Concentration is adjusted to approximately 1.5 × 105 cells/ml. Measurements were performed during the log phase, approximately 36 h after splitting. Viability of cells has been assessed to ∼95% using Trypan blue. Cells responded negatively to a test for Mycoplasma infection (MycoSPY®, Biontex).
Dynamic real-time deformability cytometry (dRT-DC) assays are carried out as reported earlier.17 Briefly, cells are centrifuged for 5 min at 200 rcf (Allegra X-15R, Beckman Coulter) and re-suspended in phosphate-buffered saline (PBS−/−, without Ca2+/Mg2+, BioWest) complemented with 1% (w/v) methylcellulose (Sigma-Aldrich) as the measurement buffer. First, the microfluidic chip is filled with the buffer in sheath flow. Second, 100 μl of the sample solution are drawn into the tubing and connected to the sample inlet of the chip. The channel is flushed with the two flows at a flow rate of 100 nl s−1 each. After dosing 10 μl with the two 1 ml syringes (1 ml Luer-Lock syringe, BD), the target flow rate of 8 nl s−1 (2 nl s−1 for sample and 6 nl s−1 for sheath flow) is set. Dynamic RT-DC acquisition is started after 12 min of equilibration time. Proper alignment of cells in the center of the channel is guaranteed by hydrodynamic focusing.
Inhibition of actin polymerization
Cytochalasin D (Sigma-Aldrich) dissolved in dimethylsulfoxide (DMSO) is added at a concentration of 1 μM to HL60 cells (approximately 1 × 106 cells per milliliter).28 Cells are incubated for 10 min at 37 °C, centrifuged for 5 min at 200 rcf (Allegra X-15R, Beckman Coulter), and re-suspended in PBS−/− complemented with 1% (w/v) methylcellulose (Sigma-Aldrich) as the measurement buffer. Corresponding vehicle is treated with 0.25% (v/v) DMSO while the control condition is untreated.
Cell shape analysis
Cell deformation analysis is carried out in a custom software for dynamic cell tracking17 adopted from Ref. 15. The camera acquires images at 3000 frames per second within a field of view of 1280 × 100 pixels, equivalent to 435 × 34 μm2, covering the complete channel length including inlet and outlet regions [Fig. 1(a)]. To reduce the processing bandwidth, the tracking algorithm analyzes a smaller region-of-interest (ROI) within the full field of view by cropping of 250 × 100 pixels around the moving cell. This moving ROI encapsulating the tracked cell starts at the left boundary and is shifted stepwise by 120 pixels to the right (with increasing x) if the cell's center of mass (COM) passes a threshold in the cropped image at approximately 70% of ROI width. Image acquisition and processing is implemented in LabView (LabView 2018, National Instruments) utilizing functions of the computer vision library OpenCV (https://opencv.org/) compiled into dynamic link libraries. After background subtraction, a binary image is generated by comparing the 8-bit image to a threshold of gray scale value 6. The binarized image enables first a border-following algorithm29 and second a convex hull algorithm30 to find contour coordinates, which are saved and used to calculate a rescaled circularity as a deformation measure,
where c, A, and P represent the circularity, the projected cell area, and the perimeter, respectively. Image acquisition, processing, and cell shape analysis is done in real-time during the experiment.
Before entering the channel, cells move slowly and, consequently, shear forces are small resulting into undeformed cell shapes close to a circle (), whereas cells are stretched into an ellipsoidal shape at the inlet due to an acceleration in the direction of flow [Fig. 1(b)]. Inside the channel, cells follow the Poiseuille flow profile and mimic a bullet-shaped steady-state at the channel end.
For each cell and at each position within the field of view, the bright-field images, contour coordinates, as well as area and circularity are stored during the experiment for later analysis. In a post-processing step, cell contour coordinates are used to calculate the time traces of the remaining shape descriptors.
Cell shape descriptors
Deformation of a cell is quantified by multiple shape descriptors. In the following, we present the definition of the cell area, the normalized front radius, the axes ratio, the Taylor deformation, the principal axes ratio, the rescaled circularity, inverse Haralick's circularity, circular variance, and the inertia ratio.9,31,32
All calculations are based on a convex hull cell contour described by a closed polygon with N points where and . The cell is aligned in the direction of flow, which we assign the x direction [Fig. 1(a)].
During image acquisition, cell contour is initially stored in Cartesian coordinates. For the shape descriptors’ principal axes ratio, inverse Haralick's circularity, and circular variance, the contour is interpreted as an angle-dependent radius function and needs to be transformed into polar coordinates. To ensure equidistant sampling in space, the contour is piecewise linearly interpolated to 500 points in total over in Cartesian coordinates . Interpolation is necessary to ensure that calculation of shape descriptors is insensitive against the number of contour points. Here, a 500 point interpolation yields an error below 0.1% (Fig. S1 in the supplementary material). Interpolated points are transformed into polar coordinates ,
where represents the coordinates of center of mass with
Area: Cell area A is defined as the cell area enclosed by its convex hull contour [Eq. (3)].
Normalized front radius: The normalized front radius represents the curvature of the cell's front part in the direction of flow [Fig. 1(a)]. For calculation, contour points under an angle of 90° of the leading part relative to the center of mass are considered [, Fig. 2(a)]. A least squares approximation enables calculating the center of the corresponding circle, from which the front radius can be derived,
If the minimization algorithm fails to find the global minimum, which is examined from a comparison of the returned with an upper threshold, set to , the front radius at the current cell position needs to be excluded.
For a size-independent representation, front radius is normalized to the equivalent radius of a cell with circular shape and of same size (). The front radius has a local minimum at the inlet in contrast to all other shape descriptors but Haralick's circularity. For comparison, we analyze the inverse normalized front radius here, referred to as the normalized front radius ,
Axes ratio: The axes ratio is defined by the side lengths of the bounding box fully enclosing the cell [Fig. 2(a)] and is calculated from length along the axis in the direction of flow divided by length along the perpendicular axis,
Taylor deformation: Taylor deformation is defined as the difference of the two axis lengths of the bounding box divided by their sum [Fig. 2(a)],
Principal axes ratio: The principal axes ratio is defined by the two principal axes lengths of a cell contour and is calculated from a covariance matrix C based on interpolated cell contour points relative to its center of mass [Fig. S1 in the supplementary material, cf. Eq. (2)],
is defined as the ratio of both eigenvalues of C [Fig. 2(a), square roots are depicted since eigenvalues have the dimension length to the power of two]. To ensure correct assignment, corresponding eigenvectors are filtered for their orientation.
Here, the principal axes ratio is calculated as the eigenvalue corresponding to the eigenvector dominated by the axis in direction of flow () divided by its perpendicular counterpart (),
Rescaled circularity: The rescaled circularity is calculated from the cell circularity c according to Eq. (1). For a completely circular cell, the rescaled circularity equals zero.
Inverse Haralick's circularity: Haralick's circularity is based on the radial distribution of all cell contour points and is calculated from the interpolated radii mean divided by their standard deviation33 [cf. Eq. (2)]. In contrast to most shape descriptors, possesses a local minimum at the channel inlet and we, therefore, use the inverse Haralick's circularity in this work, featuring improved numerical stability since the standard deviation is shifted from the denominator to the nominator [Fig. 2(b)],
Circular variance: The circular variance relates the radial variance to the squared radii mean [Eq. (2)],
Circular variance is the square of inverse Haralick's circularity.
Inertia ratio: The inertia ratio relates the area moment of inertia along the two axes of minimal and maximal moment of area [Fig. 2(a), depicted are the fourth roots since the moment of inertia has the dimension length to the power of four]. First, the area moments of inertia , , and , related to the x axis, y axis, and their diagonal, respectively, are calculated for contour points relative to their center of mass ,
The eigenvalues and eigenvectors of the area moment of inertia tensor T are determined (describing minimal and maximal moment of area and corresponding axes).
For correct assignment of extremal axes of inertia ratio, corresponding eigenvectors are filtered for orientation with respect to the direction of flow. is then calculated as the inertia dominated by axis in the direction of flow () divided by its perpendicular counterpart (),
Data filtering strategy
The tracking algorithm derives a convex hull contour of the cell from every raw image. Goodness of the contour fit is determined by comparing the convex hull area to the raw cell area [Fig. 2(a), first panel]. A cell trace is considered to be valid if at least 80% of all images inside the channel possess a convex hull area that is not more than 5% bigger than the raw cell area. If a cell trace does not cover the entire length of the channel, it is excluded from data analysis.
Cells are filtered for projected cell areas within [75, 300] μm2 in order to exclude debris as well as cell clusters. Since the stress on cell surface depends on the cell velocity, only cells with a mean velocity are considered for analysis, with standard deviation .
With all our experiments being carried out at a flow rate of 8 nl s−1, we find a range of 17.6 ± 1.8 mm s−1 for a typical experiment with 1806 cells acquired. After filtering for cell size and velocity, 1766 cells remain in the sample. This corresponds to 2.2% of the cells being excluded [Figs. 3(a) and 3(b)]
The shape dynamics of a cell translocating a microfluidic channel shows a superimposed response to an inlet stress due to its initial acceleration and to a constant stress inside the constriction. Shape mode decomposition utilizing Fourier transformation published earlier enables us to disentangle cellular response to both contributions.17 Briefly, cell contours are transformed from Cartesian into polar coordinates with respect to their center of mass.
From the Fourier transformation of this angle-dependent radius function the first ten Fourier coefficients (FCs) are considered. By symmetry arguments only, it was shown that cell response to inlet stress is reflected by the even Fourier components, whereas the odd Fourier components represent cellular response to the stress inside the channel. To discriminate between the two responses, the shapes are reconstructed from even and odd FC subsets (including FC zero) separately. Consequently, every shape of a cell at any position can be decomposed into two contours representing the response to inlet and channel stresses.
Time-dependent shape descriptor analysis
Deformation traces , with d being a time-dependent scalar representing all shape descriptors considered in this work, are the basis for characteristic time analysis, which is needed, e.g., to calculate viscoelastic material properties.
Here, we extract the characteristic time by fitting an exponential function to using a custom script implemented in Matlab R2018b (Mathworks),
with time t and free fitting parameters , , and representing the deformation amplitude of any shape descriptor discussed in this work, the characteristic time, and the steady-state offset deformation, respectively. This approach is used for fitting of traces independently of the respective analysis strategy, i.e., direct ensemble, ensemble analysis of inlet or channel effects, and single cell analysis of both effects.
Inside the channel, the exponential behavior of can be understood from a Kelvin–Voigt model with characteristic time,
because of the presence of a step stress. Here, corresponds to the viscosity and E to Young's modulus of the cell.
For direct ensemble analysis, we overlay all single traces into a master curve by calculating the median at each point in time inside the channel, using interpolation to synchronize between uneven sampling times of the traces. Fitting the master curve is performed for points in time within a range between 0% and 95% inside the channel to ensure that characteristic times are not affected by the outlet [cf. Fig. 3(b)]. Of note, since the position is defined by the center of mass, a cell can be considered inside the channel while the leading edge already crosses the outlet.
For ensemble analysis of inlet effects, we overlay all single traces into a master curve after Fourier decomposition and reconstruction from even FC. Subsequently, fitting according to Eq. (14) is performed as described above. Channel traces reconstructed from odd FC require the introduction of a lag time corresponding to 10% of the translocation time for fitting. Here, shape descriptors stay constant for a couple of frames before a response to channel stress is observed [Fig. 3(b)].
Single cell analysis is performed for traces reconstructed from even and odd FC analog to the ensemble traces described above. However, the beginning and end point for fitting need to be found for every single trace individually because the responses vary in their behavior over time. For inlet traces, the position of the minimum (if applicable) is determined from a central moving median filter of seven points length, i.e., three points before and after the current time point.
Traces are fitted from channel inlet until minimum position or channel end if no minimum exists. For channel stress traces, the initial lag phase is searched dynamically and traces are fitted from the end of the lag phase until channel outlet.
Quality of fit is assessed from the coefficient of determination , where we consider all traces with for data analysis.17 The coefficient of determination relates the residual sum of squares () to the total sum of squares (),
where is the th data point, is the corresponding fit function value, and is the data mean.
For statistical analysis, a linear mixed-effects model (LMM) is fitted to the dataset consisting of three biological replicates with paired control (control or DMSO vehicle) and treatment (CytoD), as introduced earlier.23 First, a linear mixed-effects model is used considering the treatment as the fixed effect and differences between replicates as the random effect. A second model including random effects only represents the null hypothesis. Utilizing a model comparison by maximizing the likelihood function, a -value can be extracted. This strategy is used for statistical analysis of single cell measurements where datasets consist of thousands of cells. The error of effect sizes is estimated from errors for intercept and treatment from LMM. Propagation of uncertainty is implemented utilizing a first order Taylor expansion.
Statistical analysis for ensemble measurements is performed by a -test as the number of observations is not sufficient for fitting a linear mixed-effects model. Effect size is represented by normalizing the treatment parameter (CytoD) by the vehicle parameter (control or DMSO). Effect size diminished by unity is expected to be zero for no effect. A one-parametric -test on this difference is used to test against zero and provides the corresponding -value.
Cellular shape descriptors reveal characteristic dynamics
For investigating how different shape descriptors affect the calculation of mechanical properties, suspended cells are monitored during the translocation of a microfluidic channel [Fig. 1(a)]. We use dynamic real-time deformability cytometry (dRT-DC) to analyze cell shape within a moving region-of-interest (ROI) at ten different positions between channel inlet and outlet17 (see the Materials and Methods section). At a typical flow rate of 8 nl s−1, the high temporal resolution of our experimental system enables us to observe the dynamic cell response to the hydrodynamic stress distribution for 98 frames [Fig. 1(b)]. The number of acquired frames differs for individual cells and depends strongly on the flow rate.
Derivation of cell material properties from strain requires knowledge of the stress distribution inside our microfluidic system. While stress on a cell can be obtained from an analytical as well as numerical model,6,34 multiple shape descriptors are available to quantify strain. Within this work, we investigated nine different shape descriptors as a measure for cell deformation: projected cell area, normalized front radius, Taylor deformation, axes ratio, principal axes ratio, rescaled circularity, inverse Haralick's circularity, circular variance, and inertia ratio [Fig. 2(a), see the Materials and Methods section]. All of these are geometric parameters for characterizing two-dimensional shapes by a scalar, represent different shape features, and possess different properties. Cell rotation inside the channel, e.g., impacts on normalized front radius, Taylor deformation, and axes ratio, whereas the other shape descriptors are rotational invariant, either by definition like area, rescaled circularity, inverse Haralick's circularity, and circular variance, or due to reference to principal axes of cell orientation like principal axes ratio and inertia ratio.
When monitoring the translocation of a single HL60 cell through a microfluidic constriction, the nine shape descriptors reveal characteristic dynamics [Fig. 2(b)]. Whereas area stays constant over time as expected, shape descriptors based on cell axes like axes ratio, Taylor deformation, principal axes ratio, and inertia ratio show a maximum at inlet position and are monotonous decreasing toward the channel outlet. We also observe a decrease in normalized front radius, indicating that curvature of the leading part of a cell decreases, and, hence, radius increases inside the channel. In contrast, rescaled circularity, inverse Haralick's circularity, and circular variance reveal a minimum amplitude between the maximum at channel inlet and steady-state at channel outlet. Here, we define the steady-state as a condition where cell deformation is time-independent.
Characteristic times depend on cellular shape descriptors
A cell responds to two different stress distributions inside a microfluidic channel: a peak stress at the inlet due to a reduction in cross section and subsequent acceleration, as well as a constant stress inside the channel originating from the Stokes flow.17 This superposition results in complex cellular dynamics reflected in the behavior of the geometrical descriptors.
In general, shape descriptors can be discriminated in quantities revealing a minimum in their deformation trace and those with a monotonic behavior [Fig. 2(b)]. Focusing on a deformation measure derived from circularity, we have recently shown that Fourier decomposition of shape modes enables us to disentangle the cell response to the inlet stress from the effect of the channel stress.17 Here, we expanded that analytical framework to all nine shape descriptors introduced above [Fig. 2(a)]. Performing dRT-DC, traces from more than 4800 HL60 cells in three independent experimental replicates were acquired, shapes decomposed and reconstructed from even and odd Fourier components individually. These represent cell response to inlet and channel stress, respectively (see the Materials and Methods section).
Shape descriptors reconstructed from only even Fourier coefficients reveal a peak at the channel entrance and a constant amplitude toward the outlet [Fig. 3(a) and Fig. S2(a) in the supplementary material]. In contrast, most shape descriptors reconstructed from only odd Fourier coefficients possess a constant magnitude at channel inlet, which increases until the cell leaves the constriction [Fig. 3(b) and Fig. S2(b) in the supplementary material]. Remarkably, axes ratio and Taylor deformation are both insensitive to channel stress and stay constant around one and zero, respectively. As expected, projected cell area does not change as a function of time and is, therefore, excluded from further analysis [Figs. S2(a) and S2(b) in the supplementary material]. For simplicity, data for axes ratio and circular variance are displayed in the supplementary material only due to their mathematical similarity to Taylor deformation and inverse Haralick's circularity, respectively [Figs. S2(a) and S2(b) in the supplementary material].
We analyzed the traces of all shape descriptors, first, by obtaining an ensemble average and, second, on a single cell level. The ensemble average [Fig. 3(a), bold green line, and Fig. 3(b), bold blue line] has been calculated from the median over all single cell traces [Fig. 3(a), green traces, and Fig. 3(b), blue traces] for each of the 85 interpolated positions within the microfluidic system and separated for inlet and channel effects. Fitting an exponential function [Eq. (14)] to the ensemble trace yields two characteristic times, one for the cell response to inlet stress [, Fig. 3(a), yellow lines] and one for the cell response to the constant channel stress [, Fig. 3(b), red lines]. For the latter, can be interpreted as the response of a cell [Eq. (15)],
to a creeping-flow experiment exploiting a Kelvin–Voigt model.27,35–37 With E corresponding to the cellular Young's modulus and representing its viscosity, the characteristic channel time contains the full viscoelastic information of the cell in the presence of a step stress. The transition into the channel, i.e., the deformation trace in response to the peak stress at the inlet, can also be described by an exponential function. While the time-dependent and non-uniform stress distribution renders the derivation of Young's modulus and viscosity from analytical models challenging, the corresponding characteristic time can still be understood as a quantity that contains the full viscoelastic information of the cell in response to the peak stress at the channel inlet. Both and represent two material parameters and are investigated in this work.
Next, we analyzed our data on a single cell level [Figs. 3(a) and 3(b), green and blue traces] by fitting an exponential function [Eq. (14)] to every inlet and channel trace [Fig. 3(c), inset top left panel]. Goodness of the exponential fit is assessed by the coefficient of determination (, see the Materials and Methods section). Only fits with are included in our analysis, as we already demonstrated that this threshold provides a good data representation.17 Comparison between different conditions further requires that experimental timescales given by the channel translocation time match the characteristic timescales and of individual cells.
This is specifically of importance when analyses require cells being in steady-state, which corresponds to five times the characteristic time (equivalent to more than 99% of asymptotic value of an exponential function for ). This implies that a relaxation time of, e.g., 5 ms requires a translocation time of 25 ms.
For investigating the steady-state condition, we introduced an amplitude fraction as the ratio between the shape descriptor amplitude at the channel outlet and the steady-state amplitude , both corrected for the amplitude at the channel inlet , all obtained from the exponential fit,
If the amplitude fraction of an experiment exceeds (Fig. S3 in the supplementary material), the characteristic time from the fit is considered to be a valid representation of the cellular relaxation time. Utilizing the exponential dependency between the amplitude and time domain, Eq. (17) can be rewritten by means of the relevant time scales [Eqs. (1)–(5) in the supplementary material],
Equation (18) enables the definition of an upper threshold of valid characteristic times for single cells by measuring a median translocation time and imposing a minimal amplitude fraction .
In a typical experimental setting with a 300 μm long channel and applying a flow rate of 8 nl s−1, we find a mean cell velocity of 17.6 mm s−1 for HL60 cells and a channel translocation time of 17.0 ms (Fig. S4 in the supplementary material, top left panel). For , i.e., at least 50% of a given steady-state amplitude is found inside the channel [Fig. 3(c), inset top left panel]; the threshold in characteristic time is given by , which corresponds to the condition [Eq. (18)].
In our work, we analyzed single cell traces from three biological replicates by Fourier decomposition and subsequent shape reconstruction from odd as well as even coefficients. Applying a minimal amplitude fraction as a selection criterion for valid exponential fits, the resulting characteristic times are described by a lognormal distribution [Fig. 3(c) and Fig. S4 in the supplementary material, one replicate is presented]. In general, we find , which can be explained by the short peak stress at the inlet and the longer constant stress in the channel. Interestingly, traces described by the normalized front radius demonstrate an opposite behavior, i.e., , a result that potentially originates from its sensitivity toward cell contour curvature in the direction of flow. In contrast, traces quantified by axis ratio and Taylor deformation are nearly insensitive toward the constant channel stress and the amplitude of the shape descriptors remains unaltered [Fig. 3(b) and Fig. S2(b) in the supplementary material].
Finally, we compare the single cell analyses to the results from the ensemble fit and find that the ensemble characteristic time is equal or smaller than the mode of the single cell distribution in case of the inlet response [Fig. 3(c) and Fig. S4 in the supplementary material, yellow lines] while the channel traces reveal an opposite trend [Fig. 3(c) and Fig. S4 in the supplementary material, red lines]. These qualitative differences can be understood from the different filtering strategies. While single cell analysis enables evaluation of every single trace, the individual steady-state information is hidden in the intrinsic averaging for the ensemble mean. However, on an ensemble level, Eq. (18) must be fulfilled as well, which is true for all presented conditions assuming [Figs. 3(a) and 3(b)]. Interestingly, the relationship between single cell analyses, on the one hand, and inlet vs channel response, on the other, seems to be valid also for adherent cells being brought into suspension [Fig. S5 in the supplementary material].
Measurement bias introduced by single cell analysis
We analyzed the impact of different filtering strategies on characteristic times from single cell traces captured inside the channel. As stated above, traces with are excluded from the analysis. This leads to a shift in the median characteristic time toward a lower value and introduces a systematic error into single cell analysis. Depending on , the bias can be smaller [Fig. 4(a), top] or bigger [Fig. 4(a), bottom].
We investigated this effect by simulating a characteristic time distribution defined by a lognormal function of median time and a standard deviation of , with [Fig. S6(a) in the supplementary material]. Varying the translocation time and normalization to enables calculation of the filtered fraction of the sample distribution, which corresponds to 50% of all events of the full distribution for and to 99% of all events for [Fig. 4(b), left panel]. Here, we assumed a minimal amplitude fraction of . A comparison between and , where is experimentally not accessible, reveals an apparent deviation for [Fig. 4(b), right panel]. Filtering for traces in single cell analysis with translocation times equal to at least the median relaxation time can lead to errors of up to 10%, while we find an error of only 1.5% for [Fig. 4(b), right panel, top inset].
From these calculations, the measurement bias in single cell analysis can be estimated. An evaluation of the experimentally observed median characteristic times from the eight shape descriptors indicates a range of for inlet data and for channel response [Fig. 3(c) and Fig. S4 in the supplementary material]. Normalizing a typical translocation time to each median characteristic time yields a ratio from 3.5 to 9.2 for inlet and 2.4 to 4.9 for channel response. This results in an error in single cell analysis between −1.5% and 0.9% [Fig. 4(b), right panel, right inset].
Interestingly, for or , we obtain an error of approximately 20%, which has the same magnitude as the estimated experimental error, e.g., due to pixilation effects of our image acquisition.38 Practically, these considerations imply first that the systematic bias due to the inaccessibility of the true median value of the characteristic time is independent of the shape descriptor and can be neglected. Second, our simulations verify that the steady-state of a cell can be extrapolated from translocation times as short (and even shorter) as the characteristic time and that Young's modulus as well as the viscosity can be calculated from such traces.
Measurement bias introduced by ensemble analysis
Next, we investigated the conditions under which the ensemble average is a good representation of the underlying single cell distribution. Although single cell traces are characterized by an exponential function, the median or mean of many such exponential functions does not necessarily possess exponential behavior. In fact, summation of non-linear functions does not result in a function of the same type without loss of generality.
In a simulation, we addressed this issue and generated characteristic time distributions as described above. Corresponding amplitude values of the traces were taken from a lognormal distribution with center at unity and standard deviation of [Fig. S6(b) in the supplementary material]. From each pair of characteristic times and deformation amplitudes, an exponential curve is created, either with a negative [Fig. 4(c), left panel, inset, yellow line] or a positive slope [Fig. 4(c), left panel, inset, red line], representing inlet and channel response, respectively. An ensemble response is calculated from the median over all traces. When calculating the amplitude fraction of the ensemble response while varying the translocation time relative to the median characteristic time of the single trace distribution, we observe for an increasing translocation time or channel length an increase in [Fig. 4(c), left panel]. For , we obtain , or if we impose , meaning that at least 50% of deformation amplitude are supported by acquired data points within the channel, the translocation time of cells should be chosen to be at least approximately 70% of the expected median relaxation time.
In line with experimental assays, we extract the ensemble relaxation time by fitting an exponential function to the median response. Comparing the ensemble relaxation time to the median of the real single trace characteristic time distribution , both with respect to , we find in almost all conditions a small but obvious underestimation of [Fig. 4(c), center panel]. In general, this effect is more pronounced for the inlet response, in particular, for short translocation times. An absolute error below 11% or 9% is observed for the response to inlet or channel stress, respectively, if is kept between and [Fig. 4(c), center panel, top inset]. The fact of having a finite error even for long translocations times and the apparent differences between inlet and channel computations originates from the non-trivial problem of summation over non-linear functions.
Looking at the experimental data, we find ensemble characteristic times for all eight shape descriptors of for the inlet and for the channel stress response [Figs. 3(a) and 3(b) and Fig. S2 in the supplementary material]. Normalization of the translocation time to each ensemble characteristic time yields a ratio from 3.9 to 11.3 for inlet and from 1.9 to 4.6 for channel response. This results in an absolute error below 7% for inlet and from −7% to 4% for channel, respectively [Fig. 4(c), center panel, right inset]. A comparison of these values to the systematic error introduced by excluding traces from single cell analysis shows that ensemble averaging results in larger bias but still well below other uncertainties, e.g., pixilation effects.
Since experiments do not permit access to the full characteristic time distribution in single cell analysis nor obtaining the unbiased ensemble characteristic time, we finally compare to the normalized translocation time [Fig. 4(c), right panel]. For , i.e., for short translocation times, we find an overestimation of compared to single cell analysis as many traces are included in the ensemble characteristic time that do not reach the minimal amplitude fraction. For long translocation times, we find little deviation between the ensemble average and the filtered median of the single cell analysis, as expected. Interestingly, for intermediate translocation times, , is slightly smaller than for both inlet and channel measurements.
Impact of shape descriptors on effect size after cytoskeletal modifications
Utilizing the results on comparing single cell and ensemble analysis for all eight shape descriptors, we perform measurements on HL60 cells treated with cytochalasin D (CytoD) inhibiting actin polymerization. Cells were incubated with 1 μM CytoD and the corresponding vehicle control [0.25% (v/v) DMSO], respectively (see the Materials and Methods section). Data were acquired by dRT-DC and evaluated as described above (Figs. S7–S10 in the supplementary material). For comparison, also a direct analysis of all time traces, i.e., without Fourier decomposition and reconstruction, was performed [Fig. 5(a) and Fig. S10 in the supplementary material, gray data points]. Here, the cell traces were averaged yielding one ensemble response of deformation expressed by a shape descriptor (Fig. S11 in the supplementary material). The latter approach is computationally less expensive than all other analyses and yields results rapidly but does not allow for a discrimination of inlet and channel effects. Since rescaled circularity, inverse Haralick's circularity, and circular variance reveal a minimum in amplitude inside the channel, data cannot be processed for these shape descriptors without Fourier decomposition.
We define the effect size of cytoskeletal modifications as the ratio between the characteristic times after CytoD treatment and the vehicle control, where indicates no effect. With the exception of the single cell channel data represented by the principal axes ratio [Fig. 5(a), lower panel, blue data points], all shape descriptors yield an effect size > 1 of varying amplitude and statistical significance. For the direct ensemble analysis, principal axes ratio reveals the largest effect size, while Fourier-based ensemble analysis of the inlet and channel data is best described by rescaled circularity and inertia ratio, respectively. In contrast, single cell data after Fourier reconstruction from the inlet and channel have the largest effect size when being evaluated by circular variance and inverse Haralick's circularity, respectively. In general, we find that shape descriptors calculated from ensemble data reveal larger effect sizes than single cell data for the inlet [Fig. 5(a), upper panel, yellow and green data points] as well as for channel measurements [Fig. 5(a), lower panel, red and blue data points]. It has to be emphasized that this result can be generalized for varying experimental conditions and is, e.g., robust against alterations in the hydrodynamic stress level [Fig. S12 in the supplementary material].
Finally, we combine the results of our effect size analysis and define a score for each shape descriptor,
which takes into account the magnitude of the effect size and the significance level given by the -value. The score provides a measure, how well a combination of shape descriptor and analysis strategy is suited for probing characteristic times for compound screening. Comparing all approaches presented in this work yields the highest score for single cell analysis at the inlet, a result that could originate from the peak stress at the channel entrance [Fig. 5(b)]. Interestingly, tracking and evaluating individual traces inside the channel also reveal higher scores compared to ensemble averages and emphasizes the importance of single cell analysis to extract quantitative information from compound-response assays. Being computationally less expensive, we also compare inlet and channel characteristic times calculated from ensemble averages. Here, higher scores for traces representing inlet effects were observed.
A detailed look at the different shape descriptors indicates that ensemble analysis at the inlet is best performed using circular variance (score 3.5) and the corresponding channel evaluation should be done by inertia ratio (score 3.5). In contrast, all shape descriptors with the exception of inertia ratio yield high scores for single cell analysis at the inlet (scores between 4.1 and 6.1) while channel traces should be evaluated using normalized front radius (score 4.1) and inverse Haralick's circularity (score 4.9). Ensemble analysis without shape mode decomposition yields the highest scores for principal axes ratio (score 2.5) and inertia ratio (score 2.4) with amplitudes that are equal to or lower than scores using Fourier decomposition and symmetry-based shape reconstruction.
Assessing how well a shape descriptor is suited for the interpretation of cell mechanical measurements independently of the type of analysis, a mean score of analysis strategies was calculated for inlet and channel effects individually and also combined [Fig. 5(b)]. Circular variance (score 4.8) and inverse Haralick's circularity (score 4.3) yield the largest mean score for inlet stress, whereas impact of channel stress is best described by inverse Haralick's circularity (score 3.6) and normalized front radius (score 3.5). The highest overall mean score is reached by inverse Haralick's circularity (score 3.9) and circular variance (score 3.6).
Utilizing dynamic real-time deformability cytometry (dRT-DC), we captured the dynamics of cells passing a microfluidic channel and their deformation in response to two overlapping stress distributions.17 The resulting image sequences allow us to extract various time-dependent parameters from dRT-DC data,23 which we limited within this work to scalar quantities.9 Typical examples of these shape descriptors used in microfluidic methods are axes ratio3,25,26 and circularity15,27 that represent cellular strain and are thus suitable to determine mechanical properties if the corresponding hydrodynamic stress distribution can be estimated.
While a detailed comparison of shape descriptors for suspended cells in microfluidic systems is elusive, adherent cells have already been characterized in great detail.8,9 Here, a link between morphology and cell function was established by shape and motion analysis of motile cells in combination with a standardization of shape parameters. For example, for keratocytes, a close correlation between cell state, shape, and the speed of motion has been found that depends on the self-organization of the membrane and the cytoskeleton. In consequence, cell shape can be predicted by a model based only on cell area, membrane tension, and actin content.8 Shape descriptors further allow for quantitative analysis of microscopy images to study cell morphological features in general.9
Our study aimed to expand shape analysis to suspended cells inside an extensional as well as in a shear flow. In total, nine different shape descriptors have been analyzed and compared: cell area, front radius, axes ratio, Taylor deformation, principal axes ratio, rescaled circularity, inverse Haralick's circularity, circular variance, and inertia ratio. We utilized Fourier decomposition of cellular shape modes and reconstruction from even and odd coefficients separately to disentangle the superposition of an extensional (inlet) and Poiseuille flow (constriction) inside a microfluidic channel of finite length. In both hydrodynamic geometries, cell response follows an exponential function with characteristic times and amplitudes depending on the respective shape descriptor.
Within this analytical framework, we addressed three questions: First, we investigated the conditions required for cells to reach a steady-state deformation while translocating the microfluidic channel, which is a precondition to determine cell elasticity.6,34 The steady-state essentially depends on the viscoelastic property of a cell, i.e., the characteristic time that is a priori unknown. Since constriction length is usually fixed in microfluidic systems, the only possibility (at least for volumetric flows of constant viscosity) to control the position of steady-state deformation is by translocation time via changing the flow rate. However, flow rates can only be controlled in a narrow range to ensure sufficiently high shear stress to induce cell deformation and sufficiently small cell velocities to allow dynamic tracking and reduce image blurring. Here, we introduced an amplitude fraction as a measure of how much of the steady-state amplitude is found inside the constriction. By performing dRT-DC and analyzing the resulting dynamics for different shape descriptors in combination with analytical simulations, we find for the error in being below 20%. This implies that steady-state deformation of cells can be predicted from dRT-DC traces for translocation times even shorter than the characteristic time. For , the error is only 11% and smaller.
Second, we compared single cell vs ensemble studies, which is relevant to estimate the experimental and computational effort vs data accuracy. For short translocation times, i.e., , we find an overestimation of the ensemble characteristic time as traces are included into that do not reach a steady-state. That emphasizes the need for single cell analyses if of the sample is unknown. In contrast, for long translocation times, little differences are found between single cell and ensemble studies. When comparing typical translocation times to experimental characteristic times for HL60 cells in this study, we find a ratio of approximately 5 (depending on the shape descriptor), which results in an error of 5% for ensemble and 1.5% for single cell data, which is far below other systematic errors, e.g., due to finite pixel size of our camera system.
Finally, we performed a drug-response assay using cytochalasin D since testing the sensitivity of a method toward cytoskeletal alterations is of essential importance for its application in mechanobiology. Here, we focused on filamentous actin and how the inhibition of its polymerization is reflected by alterations in different shape descriptors. Calculating a score from effect size and statistical significance, we demonstrated that single cell analysis at the inlet provides highest sensitivity almost independent of how deformation is quantified. In contrast, when measurements are limited to the constant flow profile inside the channel to compute elasticity and viscosity, inverse Haralick's circularity provides the highest score. This capability to provide a quantitative decision criterion for the selection of a shape descriptor in dependence of the hydrodynamic geometry might be of essential importance in any project where mechanical properties are being used to study the response of suspended cells to unknown compounds.
In summary, we have demonstrated the potential of different shape descriptors for quantifying the dynamics of a cell passing a microfluidic constriction. Analysis has been performed under different conditions in response to a drug inhibiting actin polymerization. We performed the data analysis on three different levels of depth, spanning a range from directly calculating the shape descriptors to performing Fourier analysis to disentangle inlet and channel contributions before investigating ensemble and single cell response, separately. Importantly, our studies demonstrate that a steady-state deformation of cells is not required to calculate viscoelastic properties and to compare the stress response of cells under different conditions, e.g., drug treatment, as along as a sufficient part of the shape dynamics inside a microfluidic channel is known and can be quantified.
Our investigations include not only an in-depth comparison of a variety of shape descriptors for viscoelastic cell characterization in microfluidic systems but also provide experimentalists with a framework to estimate systematic errors introduced by data analysis strategy and translocation time, i.e., acquisition time per cell.
See the supplementary material for supporting information and equations.
We gratefully acknowledge support from the Bundesministerium für Bildung und Forschung, Germany (ZIK grant to O.O. under Grant Agreement No. 03Z22CN11), the Deutsches Zentrum für Herz-Kreislauf-Forschung, Germany (Postdoc start-up grant to O.O. under Grant Agreement No. 81X3400107), and the Deutsche Forschungsgemeinschaft (Project No. 374031971-TRR240).
Conflict of Interest
O.O. is co-founder of Zellmechanik Dresden distributing the technology for real-time deformability cytometry. Other authors have no conflicts to disclose.
O.O. and B.F. designed the experimental assay. D.B. performed sample preparation. B.F. established and performed the measurements as well as data analysis. B.F. performed simulations. O.O. supervised the project. O.O. and B.F. wrote the manuscript. All authors reviewed the manuscript.
The data that support the findings of this study are available from the corresponding author upon reasonable request.