Microparticles with tunable, cell-like properties for quantitative acoustic mechanophenotyping

Mechanical properties of biological cells have been shown to correlate with their biomolecular state and function, and therefore methods to measure these properties at scale are of interest. Emerging microfluidic technologies can measure the mechanical properties of cells at rates over 20,000 cells/s, which is more than four orders of magnitude faster than conventional instrumentation. However, precise and repeatable means to calibrate and test these new tools remain lacking, since cells themselves are by nature variable. Commonly, microfluidic tools use rigid polymer microspheres for calibration because they are widely available in cell-similar sizes, but conventional microspheres do not fully capture the physiological range of other mechanical properties that are equally important to device function (e.g., elastic modulus and density). Here, we present for the first time development of monodisperse polyacrylamide microparticles with both tunable elasticity and tunable density. Using these size, elasticity, and density tunable particles, we characterized a custom acoustic microfluidic device that makes single-cell measurements of mechanical properties. We then applied the approach to measure the distribution of the acoustic properties within samples of human leukocytes and showed that the system successfully discriminates lymphocytes from other leukocytes. This initial demonstration shows how the tunable microparticles with properties within the physiologically relevant range can be used in conjunction with microfluidic devices for efficient high-throughput measurements of mechanical properties at single-cell resolution.

. Polyacrylamide precursor formulations for the six MP lots tested in this study, with Lot ID corresponding to IDs found on Fig. 5. "Disp." and "Cont." correspond to dispersed and continuous phase, respectively. "Fab." corresponds to fabrication, i.e., droplet generation rate calculated by dividing dispersed phase flow rate by average precursor droplet volume.

S2. Python-based GUI
Once images of fluorescent particle/cell trajectories were acquired, the 16-bit TIF images were read into the custom Python-based software for image processing, trajectory tracing, acoustic contrast or acoustic energy density calculations, and data concatenation. The typical processing pipeline started by loading the image stack into the software from the "File" tab, along with the preestablished 2D velocity matrix spreadsheet (Fig. S1), which was calculated using previously reported equations [2].
Next, in the "Image Transform" tab, image slices containing no complete (i.e., fluorescent trajectory does not traverse full channel length within field of view) particle trajectories were removed from the stack by clicking on "Delete Image Slice" button. Once all vacant micrographs were removed, the gray-scale contrast was enhanced using one of the contrast adjustment methods (Fig. S2). The transform was saved, then the background was removed with optimal "Gaussian Standard Deviation" and "Approx. Feature Size" values were identified through the rapid transform process. The transformed image was then exported in the native file type of the original image (e.g., 16-bit TIFF). The channel edge coordinates were then entered into the "X-Y Channel Edge Positions [px]" and saved by selecting "Use Discrete Channel Edge Locations" (Fig. S1). For calculating the acoustic contrast of particles or calculating the acoustic energy density, the "Acoustic Analysis" tab or "Energy Density Cal." tabs were selected, respectively. The transformed image was loaded into the active window and channel edges were displayed by selecting "Load Transformed Image" and "Display Channel Edge", respectively. Known values were entered into their corresponding sections. Once the correct values were entered, the fluorescent trajectory was traced by selecting "Trace Particle Trajectory." The trajectory was traced from the binarized image generated via Otsu thresholding [3]. Adjoining foreground pixels were classified as a foreground object and the object must traverse the entire field of view to be classified as a valid MP trajectory (i.e., contain x-coordinates spanning width of camera sensor). Once a foreground object was classified as a valid trajectory, a best fit curve of the transverse centerline was calculated. Elapsed time, longitudinal position, and transverse position data were then calculated and used to calculate the acoustic contrast factor ( ) by selecting "Calculate Acoustic Contrast Factor." Similarly, when calculating the acoustic energy density, the "Calculate Acoustic Energy Density" button was selected. The traced trajectory data, which contains elapsed time, longitudinal location, and transverse position data at pixel-level resolution, was exported by selecting "Export Trajectory Data." User specified input variables that dictated the fit of acoustic contrast or energy density were then exported by selecting "Export Acoustic Contrast Fit Data" or "Export Energy Density Fit Data," respectively ( Fig. S3). Lastly, once all images were analyzed and data exported, the "Format Files" tab was selected. Acoustic contrast or energy density data was then concatenated to common spreadsheets that were divided by the user specified "File Division Group" (Fig. S4).  Mouse scroll through image slices (for . tif stacks) 3. Image slice deletion (for . tif stacks) 4. Image alignment (for . tif stacks), normalized correlation coefficient matching based on user specified ROI in (2) 5. Interactive saved transformed image pane 6. Channel edge detection tools, user drawn edges in (5), user specified pixel locations of one channel edge  Key 25. Input for organizing data output based on longitudinal position 26. User specific file division, previously exported .csv's (fit and trajectory outputs) will be concatenated to a common .xlsx f ile specified in input bars in (27) and/or (28), but separate files will be created depending on division group with value of division group a ppended to each common file (e.g., for "Flow Rate" division, exported files: "Combined_EAC_Fit_Parameters_100.0uL -min.xlsx" with the voltages se parated by sheets) 27. Concatenation of acoustic energy density data 28. Concatenation of acoustic contrast data 25 26 27 28

S3. Acoustic energy density along chip
For each experiment, acoustic energy density was calculated during the calibration steps, prior to measuring the acoustic contrast of the polyacrylamide MPs, as described in the main manuscript in the Microparticle acoustophoretic characterization section. We observed a longitudinal dependence on the calculated acoustic energy density, where the energy density monotonically increased along the length of the channel over the region measured, with the exception of the 24 V condition at 6 mm (Fig. S5). This position-dependent energy density was taken into account in the acoustic contrast calculations.

S4. Analysis of operational parameters effect on calculated acoustic contrast
To assess systematic bias across the operating parameters, we performed a correlation analysis between the calculated acoustic contrast and the corresponding operating parameters used for each MP lot (Fig. S6). If the p-value in the Pearson correlation calculation was greater than 0.01, the operating parameter was considered uncorrelated and did not statistically contribute to the calculated acoustic contrast.

S5. Identification of multiple populations
Common statistical tests were used to assess whether the transverse position data were interpreted as a single population or as two populations, where a high-level overview is outlined in Fig. S7. First, the shape of the transverse position distribution was tested using a Tukey Lambda probability plot correlation coefficient (PPCC), where the shape parameter value corresponding to the maximum of the correlation coefficient provides insight into which type of distribution might best match the data. Common values are shown in Fig. S8. While this is not a definitive test, it provides insight into potential distributions that would model the data more appropriately and assumes symmetry. For example, if the transverse position data was best represented by a normal distribution (i.e., shape parameter, , was approximately 0.14), the population was unlikely to contain multiple subpopulations. Alternatively, a best fit shape parameter that is greater than 0.14 suggests a short-tailed distribution, and values approaching 0.5 suggest a U-shaped distribution (Fig. S8A). High-level decision tree for assessing which distribution more appropriately fits transverse position data. SSMD stands for strictly standardized mean difference and variables , , and correspond to the best fit Tukey Lambda PPCC shape parameter, prominence, and peak prominence of strictly standardized mean difference. Empty brackets "[]" when calculating peak prominence suggest lack of optimal cut point to separate the distribution into two sub-populations.
Following the Tukey Lambda PPCC, we then tested whether the transverse position data could be separated into two distributions and at what optimal value. To do this, the transverse position data were sorted and systematically separated into two distributions, with a group that that was equal to or lower than a cut value and a second group that was greater than a cut value. The cut value was incrementally increased starting at the 1 st percentile to the 99 th percentile of the data. The strictly standardized mean difference between the two separated distributions was calculated at each cut value. The strictly standardized mean difference was plotted with respect to the cut value, and the cut value that resulted in the maximum prominence along the curve was taken as the optimal value to separate the two subpopulations.
Results from the PPCC analysis suggested the distributions of the transverse position data for Donors 1 and 2 were unlikely normally distributed and therefore merited analysis as two populations. Additionally, the distribution appeared asymmetric for Donor 2 given the greater percentage of cells exhibiting high migration, which limited to applicability of the PPCC test. Confirming this, distinct maximum prominence points were present in the strictly standardized mean difference curve for Donors 1 and 2. In contrast, Donor 3 exhibited no cut value having significant prominence supporting interpretation as a unimodal population (Fig. S8B).