1 Introduction

The field of cell manipulation in microfluidic systems has grown tremendously in recent years producing multiple technical solutions to perform unit operations such as concentration, trapping, washing, alignment, stretching, sorting and separation (Lee et al. 2017; Nilsson et al. 2009). To perform these operations, various physical principles, such as dielectrophoresis, magnetophoresis, acoustophoresis, deterministic lateral displacement, optical and inertial focusing among others, have been utilized (Pethig et al. 2010; Pethig 2017; Alnaimat et al. 2018; Zhao et al. 2016; Antfolk and Laurell 2017; McGrath et al. 2014; Kampmann et al. 2018; Martel and Toner 2014; Lenshof et al. 2017). These techniques commonly rely on the physical properties of the manipulated cells such as size, density, shape, deformability, compressibility, polarizability, magnetic susceptibility and refractive index. However, not all of these physical properties are commonly known nor found in the literature. Although methods to measure some of these parameters at single cell level have been presented, they do not always provide the population-based information needed to define optimal separation conditions when isolating a specific phenotype from a complex background. High-speed measurements of some biophysical properties at single cell level have been reported for deformability (Otto et al. 2015; Hur et al. 2011) and electrical impedance (Spencer et al. 2014), where data collection in ensemble provides population-based information. In the field of acoustophoresis, efforts to measure the acoustophoretic contrast factor (Bruus 2012; Settnes and Bruus 2012; Karlsen and Bruus 2015) of single cells by tracking cell trajectories in known acoustic fields have been reported (Augustsson et al. 2010; Hartono et al. 2011; Wang et al. 2019). However, such measurements are suffering from large error estimates since single- or few-cell measurements poorly represent the statistics of a population.

Acoustophoretic manipulation of microbeads and cells is a well-established microfluidic technique (Evander and Nilsson 2012; Augustsson and Laurell 2012; Nordin and Laurell 2012; Manneberg et al. 2008; Augustsson et al. 2012; Ding et al. 2014; Petersson et al. 2007) that exploits (i) the size, (ii) the density and the compressibility of a cell relative to the carrier fluid. Although the size distribution of the sample can easily be obtained by microscopy or using Coulter Counter instruments and the density is known to some extent from the literature for the most common cell types, the compressibility of cells is rarely found. More seldom, the correlation between size and compressibility/density has been investigated.

Manipulating and separating samples in biological matrices with high precision and accuracy is fundamental in clinical studies, e.g. aiming to investigate cell differentiation or performing phenotype-specific separation (Titushkin and Cho 2006; Cross et al. 2007; Dao et al. 2003; Remmerbach et al. 2009; Constantino 2011). It is, therefore, important to know the biophysical properties of cells in given biological states. In acoustophoresis, this calls for new methods for measuring the compressibility of specific cell phenotypes.

Microparticle manipulation by acoustophoresis utilizes the acoustic radiation force (Bruus 2012; Settnes and Bruus 2012; Karlsen and Bruus 2015)

$$\begin{aligned} {\varvec{F}}^\mathrm {rad}_{}({\varvec{x}})=-{\varvec{\nabla }}U^\mathrm {rad}_{}({\varvec{x}})\,, \end{aligned}$$
(1)

where \(U^\mathrm {rad}_{}({\varvec{x}})\) denotes the acoustic radiation potential at a position \({\varvec{x}}\) (\({\varvec{\nabla }}\) is the spatial gradient). The acoustic radiation potential is

$$\begin{aligned} U^\mathrm {rad}_{}({\varvec{x}})=\frac{4}{3}\pi \,r^3\left[ \frac{1}{2}f_0^{}({\tilde{\kappa }})\,\kappa _\mathrm {f}^{}\langle p^2_\mathrm {in}\rangle -\frac{3}{4}f_1^{}({\tilde{\rho }})\,\rho _\mathrm {f}^{}\langle v^2_\mathrm {in}\rangle \right] \,, \end{aligned}$$
(2)

in which r is the microparticle radius, \(\kappa _\mathrm {f}^{}\) and \(\rho _\mathrm {f}^{}\) are the compressibility and the density of the surrounding fluid, \(\langle p^2_\mathrm {in}\rangle\) and \(\langle v^2_\mathrm {in}\rangle\) are the time-averages of the incoming pressure and velocity waves, \(f_0^{}\) and \(f_1^{}\) are the monopole and the dipole scattering coefficients, respectively,

$$\begin{aligned} f_0^{}({\tilde{\kappa }})&=\,1-{\tilde{\kappa }}\,, \end{aligned}$$
(3a)
$$\begin{aligned} f_1^{}({\tilde{\rho }})&=\,\frac{2({\tilde{\rho }}-1)}{2{\tilde{\rho }}+1}\,, \end{aligned}$$
(3b)

in which \({\tilde{\kappa }}=\kappa /\kappa _\mathrm {f}^{}\) and \({\tilde{\rho }}=\rho /\rho _\mathrm {f}^{}\) are the particle/fluid compressibility and density ratios.

In the case of a one-dimensional standing pressure wave in a microchannel of width W, see Fig. 1a, the acoustic radiation force in the direction of the wave propagation can be simplified to

$$\begin{aligned} F^\mathrm {rad}_y(y)=\frac{4}{3}\pi \,r^3\frac{\pi }{W}\varPhi ({\tilde{\kappa }},{\tilde{\rho }}) E_\mathrm {ac}^{}\sin \left( \frac{2\pi \,y}{W}\right) \,, \end{aligned}$$
(4)

where \(\varPhi ({\tilde{\kappa }},{\tilde{\rho }})=f_0^{}({\tilde{\kappa }})+\tfrac{3}{2}f_1^{}({\tilde{\rho }})\) is the acoustophoretic contrast factor, \(E_\mathrm {ac}^{}\) is the acoustic energy density and y is the direction of the pressure-wave propagation, i.e. in the direction of the channel width.

Several strategies have been proposed to measure the acoustophysical properties of cells or microbeads by exploiting acoustophoretic principles. The most common approaches utilize the measurement of the acoustophoretic terminal velocity (ATV) of single cells in a planar standing wave. The theoretical expression for the ATV is obtained by dividing the acoustic radiation force by the drag coefficient

$$\begin{aligned} v^\mathrm {rad}_y(y)=\frac{1}{6\,\pi \,\eta (y,r)\,r}F^\mathrm {rad}_y(y)\,, \end{aligned}$$
(5)

where \(\eta (y,r)\) is the apparent dynamic viscosity for a particle of radius r in a suspending fluid in confined space. The quantity \(\eta (y,r)\) accounts for the fluid bulk viscosity and the particle–wall interaction (See Supplemental Material Particle-Wall section, Figs. S1–S8). Therefore, the ATV for the microparticle in a planar standing wave is

$$\begin{aligned} v^\mathrm {rad}_y(y,{\varvec{p}})=\frac{2\,\varPhi ({\tilde{\kappa }},{\tilde{\rho }})\,r^2}{9\,\eta (y,r)}\frac{\pi \,E_\mathrm {ac}^{}}{W}\sin \left( \frac{2\,\pi \,y}{W}\right) \,, \end{aligned}$$
(6)

where \({\varvec{p}}=[r,{\tilde{\kappa }},{\tilde{\rho }},E_\mathrm {ac}^{}]\) is a set of parameters measured or assumed in the experiments. For example, by determining the trajectory, e.g. by time-lapse imaging, of a microparticle migrating in the acoustic standing wave, the ATV can be measured and provided that the microparticle radius is known, the acoustic contrast factor can be estimated (Hartono et al. 2011; Augustsson et al. 2010). Furthermore, the compressibility of the specific target cell can be either estimated from the contrast factor if the cell density is known or by assuming its density from literature. However, calculating compressibility from single-cell ATV measurements using average cell-density data introduces uncertainties due to the inherent distribution of the cell population. Alternatively, by repeating the ATV measurement on the same cell in a suspending medium of different density and/or compressibility, the cell compressibility and density could be derived. However, the complexity in undertaking such measurements does not offer a truly feasible method since the measurement on the single cell must be performed in two different buffers and then repeated for a large number of cells to obtain statistically meaningful data. Most critical to ATV-based methods is the fact that the cell-size estimation under such conditions introduce major errors.

An alternative method, iso-acoustic focusing, to measure acoustophysical properties of cells independently from the radius was proposed by Augustsson et al. (2016). Here the acoustophoretic migration of cells occurs in a density gradient which is generated by the lateral diffusion of two different media co-flowing in the microchannel. Under these conditions, each cell reaches the so-called isoacoustic point, where the acoustic impedance of a cell equals the acoustic impedance of the surrounding medium and, therefore, no acoustic radiation force acts on the cells. The cells thus maintain a fixed lateral position in the density gradient. By (i) measuring the iso-acoustic position of cells in the continuously flowing density gradient and (ii) having phenotype-specific fluorescence labelling of the different cell types, the acoustic impedance of cells carrying a specific label could be estimated in large numbers, providing population-based data independent of cell radius. However, it should be noted that using this method, the acoustic impedance of the cells was derived, but not the compressibility.

More recently a method to estimate the average compressibility of cell populations was proposed by Cushing et al.(2017). This method is based on (i) the measurement of the cell density by titrating the conditions for neutral buoyancy followed (ii) by measurement of the speed-of-sound of the cell suspension at different dilution rates. Utilizing the fact that (i) the Wood equation (Wood and Lindsay 1956), for neutrally buoyant suspensions, yields a linear relationship between the speed of sound and the cell volume fraction and (ii) given a known suspending buffer, a linear regression of the speed of sound measurement as function of the cell volume fraction provides an average measurement of the cell compressibility. A complicating fact is, however, that very large number of cells are needed for performing the speed-of-sound measurements and also the fact that the entire measurement procedure runs over a full day makes the method less convenient.

To circumvent the above shortcomings, this paper describes an alternative and novel method that enables rapid measurement of cell compressibility based on size distributions derived from Coulter Counter measurements. The method requires modest sample volumes in terms of cell number (\(\sim 10^4_{}\) cells) and is based on the statistical analysis of cell-number count at the outlets of an acoustophoretic device as function of the acoustic energy density.

The details of the statistical physical modeling and analysis is outlined in “Materials and methods”. It utilizes the measurement of the Side-Stream Recovery \(\mathrm {SSR}\), i.e. the ratio between the microparticles collected at the side-stream outlet and the total microparticles collected at the outlets (center + side). The theoretical model of \(\mathrm {SSR}\) is subsequently used to fit the experimental separation data and (i) to calibrate the device, i.e. to estimate the device effectivity factor \(\alpha\) for which the acoustic energy density is proportional to the piezo-actuator voltage squared, i.e. \(E_\mathrm {ac}^{}=\alpha \,V^2_{}\) and (ii) to measure the acoustic contrast factor and derive the cell compressibility of specific cell phenotypes.

Using this method, compressibility data for a range of clinically relevant cell types were derived, i.e. lymphocytes, granulocytes, monocytes, as well as mesenchymal stem cells, hematopoietic stem cells, and several common tumor cell lines.

2 Materials and methods

2.1 Theory

Modelling Equation (6) is the basis for the development of the separation model and the data analysis addressed in this paper. Additionally, the contribution of the acoustic streaming in the Rayleigh limit (Muller et al. 2012; Barnkob et al. 2012) at the microchannel half-height \(z=H/2\), i.e. where the particles enter the acoustophoretic separation zone, is considered. The contribution from acoustic-induced streaming is taken into account since the calibration microbeads are of two different sizes. With these assumptions, the overall microparticle velocity (radiation and streaming-driven velocity) is

$$\begin{aligned} v_y^{}(y,{\varvec{p}})=\left[ \frac{\pi }{W}\,\mu (r,{\tilde{\kappa }},{\tilde{\rho }},\eta _{}^{})+\frac{3}{4\,\rho _\mathrm {f}\,c_\mathrm {f}}\right] E_\mathrm {ac}^{}\sin \left( \frac{2\,\pi \,y}{W}\right) \,, \end{aligned}$$
(7)

where \(\mu (r,{\tilde{\kappa }},{\tilde{\rho }},\eta _{}^{})=2\varPhi ({\tilde{\kappa }},{\tilde{\rho }})r^2/(9\eta )\) is the acoustophoretic mobility for the particle. When the modeling and the analysis of particle populations, i.e. distributions, is considered, a method accounting for the spatial and parameter distributions of the particles must be considered. Here we briefly review the principle of the analysis method, while a more in-depth theoretical description of the method can be found in Garofalo (2018). This method is based on the moment analysis of the Fokker–Planck operator acting on a probability density conditioned to an initial distribution. In this regard, the equations for the modeling of acoustophoresis have to be expressed as function of spatial position \({\varvec{x}}\) and given a set of experimental parameters \({\varvec{p}}\).

Fig. 1
figure 1

(Color Online) Comparison between the results of the Monte Carlo simulations (symbols (c), bins (d)) and the analysis method (lines (c, d)) for a 1D-resonance model (a). Particle distribution in the microfluidic channel for the selected value of the acoustic energy density (b). Side-Stream recovery as function of the acoustic energy density for the two populations (c). Y-marginal distributions for the two particle populations at the outlet section for different values of the acoustic energy density (d): (left half) bins and mixture-prediction, (right half) single Gaussians predictions, vertical dashed lines are the average entering position for the microparticles

Full size image

First, we assume that the particle population can be described using a single Gaussian distribution in the coordinate+parameter space \({\varvec{q}}=[{\varvec{x}},{\varvec{p}}]^T_{}\). This hypothesis will be amended at the end of this section by introducing Gaussian-mixture models able to describe more complex microparticle distributions. The coordinates \({\varvec{x}}\) describe the position in the channel, while the parameters \({\varvec{p}}\) account for the microparticle properties, the fluid properties, and all the parameters associated with hydrodynamics and acoustics. By assuming a Gaussian distribution, the mean and the covariance are necessary to describe the distribution as it evolves during the separation. The initial distribution (particle number density) can be written as

$$\begin{aligned} p({\varvec{q}}\,|\,{\varvec{q}}_0^{},t_0^{})=\mathcal {N}({\varvec{q}}\,;\,{\varvec{m}}_0^{},{\varvec{s}}_0^{})\,, \end{aligned}$$
(8)

where

$$\begin{aligned} \mathcal {N}({\varvec{q}}\,;\,{\varvec{m}},{\varvec{s}})=\frac{1}{\sqrt{2\,\pi }^d|{\varvec{s}}|^{\tfrac{1}{2}}_{}}\exp \left[ \frac{1}{2}({\varvec{q}}-{\varvec{m}})^T_{}{\varvec{s}}^{-1}({\varvec{q}}-{\varvec{m}})\right] \, \end{aligned}$$
(9)

is a multivariate normal distribution with mean \({\varvec{m}}\) and covariance \({\varvec{s}}\), d is the number of dimensions of the coordinate+parameter space. Specifically, we assume that the microparticles enter the separation channel at half-height and a lateral position that depends on the inlet flow rate ratio (see Garofalo (2018) for a full specification of these conditions). The initial condition in the coordinate+parameter space translates into \({\varvec{q}}_0^{}=[{\varvec{x}}_0^{},{\varvec{p}}_{}^{}]^T\) and accounts for (i) the initial spatial position \({\varvec{x}}_0\) at the inlet section of the separation channel, and (ii) the values for the parameters \({\varvec{p}}\) that are constants during the migration. The acoustophoresis separation modifies the average position and the dispersion (or covariance) of the Gaussian distribution in the coordinate+parameter space according to the system (Garofalo 2018)

$$\begin{aligned} {\varvec{\dot{m}}}&={\varvec{f}}({\varvec{m}})\,,\quad {\varvec{m}}(0)={\varvec{m}}_0^{} \end{aligned}$$
(10a)
$$\begin{aligned} {\varvec{\dot{s}}}&={\varvec{J}}({\varvec{m}})\cdot {\varvec{s}}+{\varvec{s}}\cdot {\varvec{J}}^T_{}({\varvec{m}})\,, \quad {\varvec{s}}(0)={\varvec{s}}_0^{}, \end{aligned}$$
(10b)

where \({\varvec{f}}=[{\varvec{f^x}}_{},{\varvec{0}}]^T\) and \({\varvec{J}}=\partial _{{\varvec{q}}}^{}{\varvec{f}}\) are given by the specific model chosen to describe the separation protocol, i.e. Eqs. (47) in Garofalo (2018). As explained therein, the drift \({\varvec{f}}\) includes the contribution of the hydrodynamic, the acoustic, and the (constant) parameter distributions to separation and dispersion. Specifically in this paper the velocity induced by the radiation force and the streaming, i.e. (7) is used. We assume that the particle distribution during acoustophoresis can be approximated by

$$\begin{aligned} p({\varvec{q}},t\,|\,{\varvec{q}}_0^{},t_0^{})=\mathcal {N}[{\varvec{q}}\,|\,{\varvec{m}}(t\,|\,{\varvec{m}}_0^{}),{\varvec{s}}(t\,|\,{\varvec{m}}_0^{},{\varvec{s}}_0^{})]\,. \end{aligned}$$
(11)

Since we are interested in quantifying the splitting in the y-direction of the particles arriving at the outlet section, it is necessary to consider the y-marginal of the particle distribution

$$\begin{aligned} p^y_{}(y,t\,|\,m^y_{},s^{yy}_{})=\mathcal {N}[y\,|\,m^y_{}(t\,|\,{\varvec{m}}_0^{}),s^{yy}_{}(t\,|\,{\varvec{m}}_0^{},{\varvec{s}}_0^{})]\,. \end{aligned}$$
(12)

Note that, despite only the mean \(m^y\) and variance \(s^{yy}\) of the y-marginal are needed, these quantities must be computed by the complete evolution of the mean \({\varvec{m}}\) and covariance \({\varvec{s}}\) of the multivariate distribution. A theoretical approach based on an effective velocity has been illustrated in Garofalo (2018) (see Eq. (50) in the reference). It aimed to reduce the complexity of the problem by lowering the number of spatial coordinates. Therein the comparison between the Monte Carlo simulations and the evolution of an one-dimensional model based on Eq. (10) shows that these can be used to build a simplified and yet accurate model of the separation data.

As we are aiming to count the fraction of particles arriving to the side outlet, we introduce the Side-Stream Recovery (SSR) as the cumulative of the y-marginal with respect to the flow-splitting abscissa

$$\begin{aligned} \mathrm {SSR}({\varvec{m}}_0^{},{\varvec{s}}_0^{})=\frac{1}{2}\,\mathrm {erfc}\left[ -\frac{y_\mathrm {fs}^{}-m^y_{}(t^*_{}\,|\,{\varvec{m}}_0^{})}{2\sqrt{s^{yy}_{}(t^*_{}\,|\,{\varvec{m}}_0^{},{\varvec{s}}_0^{})}}\right] \,, \end{aligned}$$
(13)

where \(y_\mathrm {fs}^{}\) is the flow-splitting abscissa, \(t^*_{}\) is the separation time, i.e. for which \(m^x_{}(t^*_{}\,|\,{\varvec{m}}_0^{})=L\) and L is the length of the device, and whose calculation is addressed in Garofalo (2018).

In the case, when it is not possible to describe the particle size distribution by a single Gaussian, a mixture of K Gaussians

$$\begin{aligned} p({\varvec{q}},t\,|\,{\varvec{q}}_0^{},t_0^{})=\sum _{k=1}^K w_k^{}\,\mathcal {N}[{\varvec{q}}\,|\,{\varvec{m}}_k^{}(t\,|\,{\varvec{m}}_{k,0}^{}),{\varvec{s}}_k^{}(t\,|\,{\varvec{m}}_{k,0}^{},{\varvec{s}}_{k,0}^{})]\, \end{aligned}$$
(14)

can be used. Therefore, the theory so far illustrated is still valid provided that the definitions

$$\begin{aligned} p^y_{}(y\,|\,{\varvec{m}}_0^{},{\varvec{s}}_0^{})=\sum _{k=1}^K w_k^{} \,\mathcal {N}[y\,|\,m_k^{yy}(t\,|\,{\varvec{m}}_{k,0}^{}),s_k^{yy}(t\,|\,{\varvec{m}}_{k,0}^{},{\varvec{s}}_{k,0}^{})] \end{aligned}$$
(15)

for the y-marginal and

$$\begin{aligned} \mathrm {SSR}({\varvec{m}}_0^{},{\varvec{s}}_0^{})=\sum _{k=1}^K \frac{w_k^{}}{2}\,\mathrm {erfc}\left[ -\frac{y_\mathrm {fs}^{}-m^y_k(t^*_k\,|\,{\varvec{m}}_{k,0}^{})}{2\sqrt{s^{yy}_k(t^*_k\,|\,{\varvec{m}}_{k,0}^{},{\varvec{s}}_{k,0}^{})}}\right] \, \end{aligned}$$
(16)

for the SSR are used. Here the index k indicates the kth Gaussian in the mixture, and \(w_k^{}\) are weights derived from spatial and parameter distributions. Specifically in this paper, these weight factors are derived by the microparticle-size histogram.

Figure 1 illustrates the results obtained by applying the analysis above to a model of acoustophoresis where the superposition of three Gaussians per particle population is considered. Panel (d) shows the two spatial distributions, i.e. the y-marginal (15), of microparticles (blue and orange) at the outlet section for different acoustic energy density \(E_\mathrm {ac}^{}\). The vertical lines define the flow-splitting abscissa and divide the portions of fluid towards the side (green shaded area) and the center (red shaded area) outlet. The bins represent the Monte Carlo simulations, while the lines are the results of the above-illustrated analysis when the distribution is described by the superposition of three Gaussians for each population. The evolution of these three Gaussians are depicted on the right. As the acoustic energy density increases the microparticle distributions move towards the microchannel center. The orange-colored particles have higher mobility than the blue and they cross the flow-splitting abscissa for a lower value of the acoustic energy density. Panel (c) shows the results of this process in terms of \(\mathrm {SSR}\) as function of \(E_\mathrm {ac}^{}\) for the two particle populations: symbols are the Monte Carlo simulations, while lines are the results obtained using (16). Finally, panel (b) depicts the spatial positions of the particles within the microchannel for the selected value of acoustic energy density \(E_\mathrm {ac}^{}=3.5\,\hbox {J m}^{-3}\), for which it has the maximum possible separation of the two particle populations.

Data analysis The theoretical SSR can be used to fit experimental SSR data and thus determine unknown experimental quantities related to either (i) the device operational conditions, specifically the relationship between the acoustic energy density and the applied voltage \(E_{\mathrm {ac}}=\alpha \,V^2\), or (ii) the microparticle properties, i.e. the compressibility. The determination of \(\alpha\) is performed by fitting the experimental SSR for the calibration microbeads using the cost function

$$\begin{aligned} f_\mathrm {cal}(\alpha )&=\sum _{i=1}^N\varDelta \mathrm {SSR}^2_{}(\alpha \,|\,V_i,...)\nonumber \\&=\sum _{i=1}^N \left[ \mathrm {SSR}^\mathrm {th}(\alpha \,|\,V_i^{},...)-\mathrm {SSR}^\mathrm {cal}_i \right] ^2\,, \end{aligned}$$
(17)

where the ellipsis indicate a set of fixed parameters which includes the physical properties of the calibration microbeads. The device effectivity \(\alpha\) is sought by minimizing the cost function

$$\begin{aligned} \alpha ^*_{}=\arg \underset{\alpha }{\min }\,f_\mathrm {cal}^{}(\alpha )\,. \end{aligned}$$
(18)

When \(\alpha\) is known from the calibration procedure, unknown microparticle properties can be used as fitting parameters and thus estimated. Specifically for determining the microparticle compressibility, the measurement is performed using the cost function

$$\begin{aligned} f_\mathrm {meas}(\kappa )=\sum _{i=1}^N \left[ \mathrm {SSR}^\mathrm {th}(\kappa \,|\,\alpha ^*_{},V_i^{},...)-\mathrm {SSR}^\mathrm {meas}_i \right] ^2\,, \end{aligned}$$
(19)

that subjected to a minimization condition yields the measured compressibility

$$\begin{aligned} \kappa ^*_{}=\arg \underset{\kappa }{\min }\,f_\mathrm {meas}^{}(\kappa )\,. \end{aligned}$$
(20)

As the experimental procedure could be performed on devices with different performance, we rescale the experimental data by computing the (square root of) acoustic energy density corresponding to the different \(\alpha ^*_{}\).

2.2 Materials and experiments

Polymer microbeads Polystyrene beads (Sigma-Aldrich, St. Louis, MO, USA) in sizes \(2r\simeq 5\,~\upmu\hbox {m}\) (\(\mathrm {PS5}\upmu\hbox{m}\)) and \(2r\simeq 7~\upmu\hbox{m}\) (\(\mathrm {PS7}\upmu\hbox{m}\)) were used as reference beads. The density and the compressibility of these microbeads are given in literature Cushing et al. (2017).

Ethical statement Sampling of blood, peripheral blood progenitor cell (PBPC) products and bone marrow was approved by the Regional Ethical Review Board at Lund University, and samples were obtained from donors after informed consent at Lund University Hospital.

Blood samples Blood samples were collected from healthy volunteers in Vacutainer tubes (BD Bioscience, San Jose, CA, USA), containing ethylenediaminetetraacetic acid (EDTA) as anticoagulant. Samples were stored at RT and used within the day of collection. Blood was selectively lysed with BD Pharm Lyse lysing solution (BD Bioscience, San Jose, CA, USA) for viable, non-fixed white blood cell (WBC) samples.

CD34+ Hematopoietic stem cells. Mononuclear cells (MNC) were isolated from PBPC products by Ficoll density gradient centrifugation using Ficoll-Paque PREMIUM (1.078 g/mL, GE Healthcare Bio-Sciences AB, Uppsala, Sweden). PBPC products were collected from healthy donors undergoing standard mobilization treatment, as reported previously (Urbansky et al. 2016). Further isolation of CD34+ hematopoietic stem cells (HSC) from MNC was done with the MACS human CD34 MicroBead Kit (Miltenyi Biotec, Bergisch Gladbach, Germany) according to the manufacturer’s instructions.

Mesenchymal stromal cells from human bone marrow MNCs were isolated from fresh bone marrow aspirates using Ficoll density gradient centrifugation with Ficoll-Paque PREMIUM (GE Healthcare) according to the manufacturer’s instructions. Mesenchymal stromal cells (MSC) in passage 1 were obtained by cultivating bone marrow MNC in xenofree medium as described by Schallmoser (Schallmoser et al. 2007).

Cultured cancer cell lines The human cancer cell lines DU145 (prostate cancer) and MCF7 (breast cancer) were obtained from ATCC (American Type Culture Collection, Manassas, Virginia, USA), while SH-SY5Y (neuroblastoma) was purchased from DSMZ (Leibniz Institute, German Collection of Microorganisms and Cell Cultures, Brunswick, Germany). Cells were cultured according to the supplier’s recommendation. The head and neck squamous carcinoma cell line LU-HNSCC-25 and LU-HNSCC-26 were established in Lars Ekblad’s laboratory (Lund University, Lund, Sweden). Dulbecco’s modified Eagle’s medium supplemented with 10% fetal bovine serum, 100 units/mL penicillin, and 100 units/mL streptomycin sulfate was used and the cells were grown at \({37^{\circ }}\hbox {C}\) under humidified atmosphere with 5% \({\mathrm {CO}_{2}}\).

Acoustophoretic separation The acoustophoretic setup including the microfluidic chip design (Augustsson et al. 2012), the acoustic actuation and the fluidic setup has been described previously (Urbansky et al. 2017). Shortly, the chip comprised a prefocusing channel (\(300~\upmu\hbox{m}~\times ~150~\upmu\hbox {m}\) cross section) where particles are ensured to be in the same flow vector using 2D acoustic focusing with a 5 MHz transducer. This step effectively allows particles to be separated only on their acoustophysical properties in the separation step and not by their position in the flow. The separation step comprises a channel (\(375~\upmu\hbox{m}~\times ~150~\upmu\hbox {m}\) cross section) where the particles are subsequently being focused from the side outlets into the central outlet with increasing voltage applied to the 2 MHz transducer. To determine the acoustophoretic properties of unknown beads or cells, an equal mixture of the \(\mathrm {PS5}\upmu\hbox{m}\) and \(\mathrm {PS7}\upmu\hbox {m}\) reference beads as well as microbeads/cells of interest with a final concentration of \(150~\times ~ 10^{3}\,\hbox {mL}^{-1}\) was prepared. The sample was injected at \(Q=100\,\upmu\hbox {L}~\hbox {min}^{-1}\) into the acoustophoresis chip, the central-inlet buffer was infused at \(Q=300\,\upmu\hbox {L}~\hbox {min}^{-1}\), and the amplitude of the \(2\,\hbox {MHz}\) transducer for the separation channel was varied to generate the experimental SSR data points.

Flow cytometry Particle numbers for the side and center fractions were recorded on a FACS Canto II flow cytometer (BD Bioscience) and analyzed using FlowJo software (Tree Star Inc., Ashland, OR, USA). Polystyrene beads were identified by their characteristic forward scatter (FSC) and side scatter (SSC) signal. FITC-marked melamine resin microparticles were identified based on their FSC/SSC signal and their positive FITC signal. In all the microbead experiments, cytometry data were gated for singlets. In cell experiments, propidium iodide (PI, BD) was added before flow cytometer analysis for dead cell exclusion. Cancer cells and white blood cells were identified as PI-negative cells and based on the FSC/SSC signal. HSC cells were further labelled with CD45-FITC (clone HI30, BD) and CD34-PE (clone 581, BD). The antibody staining was added before acoustic separation and incubated for 15 min in the dark at room temperature. For information on gating strategies and how FACS plots were interpreted, see supplemental information Fig. S9.

Size determination The size distribution of each bead type and cell type used for the acoustophoretic experiments was measured with a Multisizer 3 Coulter Counter (Beckman Coulter, Brea, CA, USA).

3 Results

Calibration The calibration procedure, i.e. the determination of the device effectivity \(\alpha\), was performed using the SSR separation data for polystyrene microbeads with nominal diameters \(2r=5\,\upmu\hbox{m}\) and \(2r=7~\upmu\hbox {m}\). The physical parameters for the calibration microbeads have been summarized in Table 1. The parameters \(w^r_k\), \(m^r_k\) and \(\sigma ^{r}_k\) are the weight, the mean radius and the standard deviation for the kth Gaussian used for the reconstruction of the size distribution, i.e. \(p(r)=\sum _k w^r_k\mathcal {N}(r;m^r_k,\sigma ^{r}_k)\), \(\kappa\) is the compressibility and \(\rho\) is the density reported (Cushing et al. 2017).

Table 1 Compressibility \(\kappa\), density \(\rho\) (Cushing et al. 2017) and Gaussian mixture weights \(w_k^r\), mean \(m_k^r\), and standard deviation \(\sigma ^{rr}_k\) for polystyrene microbeads
Full size table
Fig. 2
figure 2

(Color Online) Size histograms for \(\mathrm {PS}5\upmu \hbox {m}\) (blue) and \(\mathrm {PS}7\upmu \hbox {m}\) (orange) microbeads. Symbols are the Coulter Counter data while the lines are the interpolation with the Gaussian mixtures (parameters in Table 1). Each Gaussian is plotted in dashed line

Full size image
Fig. 3
figure 3

(Color Online) a Calibration functions \(f_\mathrm {cal}^{}\) as function of \(\alpha\) (symbols) for three different calibration experiments. Continuous lines are interpolating parabolas. Vertical dashed lines (abscissae) are the minima \(\alpha ^*_h\) each corresponding to the hth experiment. b Side-Stream Recovery \(\mathrm {SSR}\) for \(\mathrm {PS5}\upmu\hbox {m}\) (blue) and \(\mathrm {PS7}\upmu\hbox {m}\) (orange) microbeads as function of the square-root of the acoustic energy density \(E_\mathrm {ac}^{1/2}\). Symbols are the experimental data and lines are the fitting curves. The measured device effectivity was found to be \(\alpha ^{}_{}=2.475(3)\,\hbox {J m}^{-3}\,\hbox {V}^{-2}\). Dashed lines refer to the worst-case error \(\kappa \mp \sigma ^\kappa _{},\rho \pm \sigma ^\rho _{},\alpha \pm \sigma ^\alpha _{}\)

Full size image

Figure 2 shows the Coulter Counter data (symbols) for the two particles, i.e. \(\mathrm {PS}5\upmu \hbox {m}\) (blue) and \(\mathrm {PS}7\upmu \hbox {m}\) (orange), with the corresponding interpolations (solid lines) derived by mixing the two Gaussians (dashed lines). The fitting procedure provides the Gaussian weights \(w_k^r\) to be used in (16) in the form

$$\begin{aligned} \mathrm {SSR}(V)=\sum _{k=1}^K \frac{w_k^{r}}{2}\,\mathrm {erfc}\left[ -\frac{y_\mathrm {fs}^{}-m^y_k(t^*_k\,|\,m_{k}^{r})}{2\sqrt{s^{yy}_k(t^*_k\,|\,m_{k}^{r},s_{k}^{rr})}}\right] \,. \end{aligned}$$
(21)

The theoretical SSR is then fitted to the experimental SSR (separation data) with \(\alpha\) as the fitting parameter by minimizing the calibration function for both microbead sizes

$$\begin{aligned} f_\mathrm {cal}^{}(\alpha )=\sum _{i=1}^N \varDelta \mathrm {SSR}^2(\alpha \,|\,V_i,\mathrm {PS5}\upmu \hbox {m})+ \varDelta \mathrm {SSR}^2(\alpha \,|\,V_i,\mathrm {PS7}\upmu \hbox {m})\,. \end{aligned}$$
(22)

Panel (a) in Fig. 3 shows some examples of calibration functions (symbols) as a function of the device effectivity \(\alpha\). Solid lines are interpolating parabolas. The calibration function is computed by accounting for the error contributions of both \(\mathrm {PS}5\upmu\hbox {m}\) and \(\mathrm {PS}7\upmu\hbox {m}\) microbeads. The vertical dashed lines are the abscissae corresponding to the minima for the different calibration functions.

A second-order statistical analysis on the different values \(\alpha ^*_{}\) gives the mean \(\alpha\) and the standard deviation \(\sigma ^\alpha _{}\). These are used to compute an averaged \(\mathrm {SSR}\) and to check the quality of this interpolation against the worst-case errors for the three experiments. An example of this procedure is shown in Fig. 3b, where the data for the three experiments have been rescaled with the corresponding acoustic energy density \(E_{\mathrm {ac},h}^{}=\alpha ^*_h\,V_{h}^2\) for the hth experiment.

The calibration procedure was also tested on different devices, accounting for the fact that bulk acoustophoresis assemblies commonly vary in performance from device to device, see values reported in Table 2.

Table 2 Operating frequency f, device effectivity \(\alpha\) and average residuals for two different devices
Full size table

Validation To validate the general applicability of the SSR method, we measured the compressibility of melamine microbeads with nominal diameter \(2r=4\,\upmu \hbox {m}\) (\(\mathrm {MA4}\upmu \hbox {m}\)) using the two devices in Table 2.

Fig. 4
figure 4

(Color Online) a Size histograms for \(\mathrm {MA}4\upmu\hbox {m}\) microbeads: symbols are the Coulter Counter data while the lines are the interpolation with Gaussian mixture with the parameters in Table 3. Each Gaussian is plotted in dashed line. b Side-Stream Recovery \(\mathrm {SSR}\) as function of the square-root of the acoustic energy density \(E_\mathrm {ac}^{1/2}\): (symbols) experimental data (\(N=3\)) and (lines) theory. The average device effectivity is \(\alpha =2.487(2)\,\hbox {J m}^{-3}\,V^2\), while the derived compressibility for MA is \(\kappa =108(6)\,\hbox {TPa}^{-1}\). Dashed lines refer to the worst-case error \(\kappa {\mp }\sigma ^\kappa ,\rho \pm \sigma ^\rho _{},\alpha \pm \sigma ^\alpha\)

Full size image
Table 3 Density \(\rho\) and Gaussian mixture parameters weights \(w_k^r\), mean \(m_k^r\), and standard deviation \(\sigma ^{rr}_k\) for melamine microbeads
Full size table

Panel (a) in Fig. 4 shows the size-histogram for the \(\mathrm {MA4}\upmu\hbox {m}\) microbeads. When the size-histogram is interpolated with a Gaussian mixture (parameters in Table 3), we observed the occurrence of a bimodal distribution that features a main peak at \(m_1^r=1.775\,\upmu\hbox {m}\) and a secondary peak at \(m_3^r=2.226\,\upmu\hbox {m}\). This secondary peak is due to microbead pairing in doublets and it can be checked by computing the equivalent radius of doublets \(2^{1/3}m_1^r\simeq 1.256\times 1.775\,\upmu \hbox {m}\simeq 2.229\,\upmu \hbox {m}\) which is indeed very close to \(m_3^r\). For this reason, we excluded the third Gaussian from the Gaussian mixture and we used the density provided by the manufacturer, \(\rho =1510(15)\,\hbox {kg m}^{-3}\). Based on these input data, device 1 yielded a compressibility of \(\kappa _{}^{}=108(6)\,\hbox {TPa}^{-1}\) while for device 2, \(\kappa _{}^{}=111(7)\,\hbox {TPa}^{-1}\). These are intermediate values compared to the value reported by Cushing et al. (2017) and that used by Settnes et al. (2012).

Measurements of cell compressibility When the device effectivity \(\alpha ^*_{}\) has been estimated, the compressibility of cells can be measured by analyzing the separation data for a mixed suspension of reference microbeads and cells. The spiking of reference beads into the sample ensures accurate estimation of the sample properties by simultaneous calibration of the device effectivity \(\alpha ^*_{h}\). The compressibility of several different cell types was measured using the developed method. In the following, we illustrate the measurement of the compressibility \(\kappa ^\mathrm {DU145}\) for the prostate cancer cell line DU–145.

Fig. 5
figure 5

(Color Online) Size histogram for DU–145: (symbols & shaded area) average and standard deviation for Coulter Counter data (three measurements), and (lines) the interpolation with Gaussian mixture with parameters reported in Table 4

Full size image
Table 4 Prostate cancer cell line DU–145 parameters for the approximation of the size histogram distribution by a Gaussian mixture
Full size table
Fig. 6
figure 6

(Color Online) Side-Stream Recovery \(\mathrm {SSR}\) as function of the square-root of the acoustic energy density \(E_\mathrm {ac}^{1/2}\): (symbol) experimental data, and (lines) theory. The average device effectivity is \(\alpha ^{}_{}=2.475(3)\,\hbox {J m}^{-3}\,V^{-2}\), while the derived compressibility for DU–145 is \(\kappa ^\mathrm {DU145}_{}=423(4)\,\hbox {TPa}^{-1}\). Dashed lines refer to the worst-case error \(\kappa {\mp }\sigma ^\kappa ,\rho \pm \sigma ^\rho ,\alpha \pm \sigma ^{\alpha }\)

Full size image

Figure 5 shows the average Coulter Counter data for DU–145 (symbols and shaded areas) based on three repeats interpolated with the Gaussian mixture (continuous line) whose parameters are reported in Table 4. The symbols represent a downsample of the average measurements while the shaded area indicates the real standard deviation around the average. The dashed lines refer to the three Gaussians in the mixture.

The results of the calibration procedure together with the simultaneous measurements of the DU-145 cells are illustrated in Fig. 6. The experimental data have been linearly interpolated to account for the different device effectivities in the three repeats, and subsequently average values (symbols) with the corresponding standard deviations (error bars) have been derived by statistical analysis of the interpolated values. The value of the compressibility was derived by assuming our previously measured density for the DU–145 cell line, \(\rho ^\mathrm {DU145}=1062(10)\,\hbox {kg m}^{-3}\), for which we assume a larger uncertainty (\(\pm 1\%\)) to account for the possible difference in cell lines between the present work and Cushing et al. (2017). The resulting compressibility for the DU–145 cell line was found to be \(\kappa ^\mathrm {DU145}_{}=423(4)\,\hbox {TPa}^{-1}\), where the uncertainty was estimated by adapting the results for the worst-case SSR curves (dashed lines).

The SSR method was subsequently used to estimate the compressibility for a wide range of cell types commonly used in cell biology or in clinical work, Table 5. Note that the compressibility for white blood cells (WBC) as a single population was measured and also using a Gaussian mixture the compressibility for WBC subpopulations: lymphocytes, granulocytes and monocytes was measured.

Table 5 List of the values for the device effectivity \(\alpha\), assumed density \(\rho\) and estimated compressibility \(\kappa\) of different cell populations
Full size table

4 Discussion

The SSR method has shown a good ability to predict the compressibility of polymer microbeads and cell populations. Also, the compressibility of WBCs and the subpopulations has been computed with a very low error when compared to literature values. The measured compressibility value for DU–145 and MCF–7 cancer cells was found to be approximately 10% larger than those measured by the neutrally buoyant speed-of-sound method (Cushing et al. 2017), but still within ranges reported by other authors (Hartono et al. 2011). The reason of this discrepancy is most likely ascribed to different handling of the cells during the experiments. In the SSR method, the environmental conditions (pH, isotonicity, etc.) can be chosen to preserve physiological conditions. Furthermore, the SSR method only measures live cells by discriminating dead or dying cells in the cytometry gating. In the macroscopic method of Cushing (Cushing et al. 2017), the measurement process runs over a full day, involving steps such as density titration and centrifugation, that could potentially alter the cell viability. The outcome of the measurement includes both live and dead cells. However, the neutrally buoyant speed-of-sound method by Cushing et al. (2017) has the advantage of independently measuring the compressibility and the density.

When compared to macroscopic or single cell measurements, the SSR method offers the following benefits. It requires a much smaller sample (ca 10–30 k cells for each SSR curve) than the macroscopic method which requires \(100\times 10^6\) cells, and at the same time provides more detailed information on the distribution of the physical parameters across the cell populations and possibly their correlations. Single cell measurements inherently requires smaller sample volume but, therefore, do not provide statistically meaningful information. The measurement is also biased by a large uncertainty when the radius is determined by optical setups for cell-tracking experiments.

5 Conclusions

In this paper, we present a new method for the measurement of cell compressibility based on statistical modeling and data analysis of free-flow acoustophoretic separation experiments. The experimental data are obtained by measuring the cell/microbead counts at the outlet sections of an acoustophoresis separation while sequentially increasing the acoustic energy density. The magnitude of the migration speed depends on the particle properties; hence, a connection between the separation data, the acoustic energy density and the physical properties of the particles can be established in terms of the statistics for the cell/bead count at the different outlets.

The SSR method represents a compromise between single-cell analysis and the neutrally buoyant speed-of-sound method. Indeed, it requires minimal cell number and obtains statistically meaningful measurements based on population data. Most importantly the measurements are derived for live cells.