Quantifying Acoustophoretic Separation of Microparticle Populations by Mean-and-Covariance Dynamics for Gaussians in Mixture Models

Quantifying Acoustophoretic Separation of Microparticle Populations by Mean-and-Covariance Dynamics for Gaussians in Mixture Models


A method for the quantification of acoustophoretic separation and dispersion for microparticle populations featuring continuously distributed physical parameters is presented. The derivation of the method starts by (i) considering the equation of motion for a particle ensemble in the coordinate+parameter space, (ii) performing moment analysis on the transport equation for the probability density function (PDF), and (iii) expanding up to the first-order the drift (and the diffusion coefficient) around the mean of the PDF. Following these steps, a system of ordinary differential equations for the evolution of the mean and the covariance in the coordinate+parameter space is derived. These differential equations enable for the approximation of the acoustophoretic separation dynamics of particle ensembles by using a gaussian mixture for which the mean and the covariance of each gaussian evolve according to the mean-and-covariance dynamics. The approximation property of this method is shown by comparison with direct numerical simulations of particle ensembles in the cases of prototypical models of acoustophoretic and free-flow acoustophoretic separations for which the particle populations are distributed according to the radius. Furthermore, the indicators for quantifying free-flow acoustophoretic separation performance are introduced, and a method for the inference of particle-histogram parameters is illustrated.

I Introduction

Acoustofluidics is a microfluidic technology that using acoustic waves is able to perform separation of microbeads and viable manipulation of cells Burguillos et al. (2013); Wiklund (2012). Indeed, by exploiting the interaction between acoustic pressure waves and a carrier-fluid suspension of microbeads/cells at microscale Bruus (2012); Settnes and Bruus (2012); Karlsen and Bruus (2015), acoustophoresis Bruus et al. (2011) is able to trap Evander and Nilsson (2012), wash Augustsson and Laurell (2012), concentrate Nordin and Laurell (2012), align Manneberg et al. (2008) and separate the suspended microparticles Augustsson et al. (2012); Ding et al. (2014); Petersson et al. (2007). The ability to separate microparticles is based on the different particle properties, such as compressibility and density. Specifically for cells, the different physical properties are associated with biological differentiation, type-uniformity and pathological conditions Titushkin and Cho (2006); Cross et al. (2007); Dao et al. (2003); Remmerbach et al. (2009); Constantino (2011).

The microparticle physical parameters that appear in the acoustophoretic force expression Bruus (2012); Settnes and Bruus (2012); Karlsen and Bruus (2015) are not well-represented by unique values, e.g. single values for the radius, the compressibility and the density, but they occur as distributions for the microparticle populations. Therefore, a model for the quantification of acoustophoresis must incorporate a mechanism that, taking into account for the statistics of the samples, allows to predict a continuous differentiation in the microparticle population trajectories and thus in the separation performance. However, the present models of acoustophoretic trajectories rely on statisticsless descriptions that do not quantify the impact of the continuously distributed particle parameters on the separation performance Muller et al. (2012, 2013); Garofalo (2014). Furthermore, since the acoustophoresis outcomes are directly related to the distribution of the physical properties, it is of interest to establish if (i) assuming the knowledge of the device features by performing hydrodynamic and acoustic calibration, and (ii) measuring the separation performance is possible to determine the distribution of the physical parameters for the particle population.

A possible and straightforward solution to overcome the drawbacks of the present models is to evolve particle ensembles that are normally distributed in both parameter and space Simon et al. (2017). However, the limitation of this kind of techniques becomes apparent when the parameter and/or the spatial distributions are not gaussians, and even more in parameter-estimation procedures which, being based on multiple calculations, must be extremely cheap in terms of the computational cost associated with a single calculation, i.e. the evolution of a single gaussian.

A more convenient method that (i) evolves the mean and the covariance of a normally distributed ensemble Garofalo (2014), and (ii) approximates the particle distribution by using a mixture model with gaussian kernels is proposed in this paper. For that, the proposed method can be addressed as “mean-and-covariance dynamics for gaussians in mixture models”, or briefly MCDGM.

A method for the evolution of the mean and the covariance of particle ensembles can be traced back to the stochastic linearization methods, that are widely used in mechanics Editor (2012); Socha (2007), and recently for the quantification of dispersion in acoustophoretic models Garofalo (2014). Stochastic linearization methods can include higher-order moments, but then closure assumptions are needed and the reconstructed PDF can violate the positivity assumption. The only difference is that in modeling acoustophoresis of microparticle populations, the thermal fluctuation, i.e. Brownian motion, can be neglected as this becomes relevant only for nanoparticles. For completeness, in this paper the derivation of the mean-and-covariance dynamics retains the diffusion term, that is dropped when the method is applied to acoustophoresis of microparticle population. The mean-and-covariance dynamics can be also framed within the moment analysis techniques Brenner (1993), that, together with the PDF reconstruction, have been used in the context of quantifying the dispersion in microfluidic devices, such as in Deterministic Lateral Displacement separators Cerbelli et al. (2013, 2015). Since the MCDGM method can be used to approximate the actual particle distribution at the outlet section of the device, it also provides the indicators necessary to quantify the acoustophoretic separation performance.

In order to illustrate the application of the MCDGM method for the quantification of acoustophoretic separation of microparticle populations, this manuscript is organized as follow. Section II (A) reviews the derivation of the mean-and-covariance dynamics by (i) introducing the equation of motion for a particle ensemble in the state space, (ii) introducing the associated transport equation, and (iii) performing moment analysis with linearization of the drift and the diffusion around the mean of the PDF. Section II (B) (i) translates the mean-and-covariance dynamics from the state space to the spatial+parameter space by providing the explicit expressions for the evolution of the spatial average and the spatial/mixed-covariance of a single gaussian and (ii) introduces the gaussian mixture approximation for the parameter marginal and for the reconstruction of the spatial marginal. Section III specializes the MCDGM method to the study of acoustophoretic separation by showing the comparisons with particle ensemble simulations for (A) a prototypical model of acoustophoretic separation, (B) the buffer-dependent separation of RBC and WBC similar to that presented in Urbansky et al. (2017), and (C) the 3D simulations for free-flow acoustophoresis in a rectangular microchannel. Finally (D) the application of the method in the estimation of particle size histogram is illustrated.

Ii Theory

ii.1 Mean-and-Covariance Dynamics

Let us consider the nonlinear stochastic differential equation in the Itô sense Risken (1996); Frank (2010)


where with , are the realizations of the random process in the -dimensional state-space, is the drift, and is the standard deviation matrix. The latter trasforms the differential of the multivariate Wiener process defined by


into the displacement for the states . In equation (2), is the unit tensor, is a time-translation, and are the expected value and the cross-covariance, respectively. The covariance can be written as . Equation (1b) represents the initial condition for the realizations of Eq. (1a) in terms of the realizations that is distributed according to a probability density function .

Equation (1) corresponds to the (forward) Fokker-Planck equation for the probability density  Risken (1996)


conditioned for by


In equation (3) the Fokker-Planck forward operator (assuming Einstein notation)


includes the drift , and the diffusion matrix


This relation can be used to derive the diffusion contribution to the Fokker-Planck operator when the Itô process Eq. (1a) is known, as well as to construct the Itô process when the Fokker-Planck operator is given Risken (1996); Frank (2010). The second derivation is performed by computing as the Choleski decomposition of the diffusion matrix , which is indeed defined by Eq. (6).

Alongside the Fokker-Planck forward operator defined in Eq. (5) is possible to introduce the backward operator Risken (1996)


as the state-space adjoint of the Fokker-Planck forward operator , defined by , where and/or satisfy certain regularity conditions for .

The dynamics of the first-order moment is derived from Eq. (3) multiplying by , integrating over the state-space and using the definition of the backward operator ( and ), obtaining


where is meant the expectation of at time for a distribution that at time “occupied” the states , or conditioned to Eq. (4). Noting that and , equation (8) can be rewritten in terms of the drift


and this equation shows that the dynamics of the first-order moment is independent on the diffusion matrix. An analogous derivation can be performed for computing the dynamics of the covariance, that results


Expanding in Taylor series up to the first-order the drift and the diffusion matrix


and substituting these expansions into Eqs. (9) and (II.1) it has


that using vector analysis notation, read as


where . The initial condition for this set of ODE is


Equations (12), or equivalently Eqs. (13), represent the set of differential equations here briefly addressed as mean-and-covariance dynamics, while Eqs. (14) are the corresponding initial conditions.

It must be noted that in the case when the drift has an implicit dependence on the parameters, i.e. , the Jacobian of the drift transforms according to


this representation is useful when the parameters appear in functions of the device or suspension features, e.g. particle compressibility and density in the acoustophoretic contrast factor.

ii.2 Dynamics of Microparticle Populations

As stated in the introduction, in order to quantify the separation and dispersion of microparticle populations it is necessary to devise a method that accounts for the statistics of the sample. The minimum requirement for this model is that it should be able to deal with normally distributed statistics. This restriction is discussed and amended at the end of this section by using a gaussian mixture to approximate arbitrary PDFs.

As first step, the state-space is split into a coordinate subspace and a parameter subspace, namely with , , and such that . As a consequence of this splitting, the coordinate marginal and the parameter marginal are


respectively, where the coordinate marginal is the actual distribution of the particle as it is seen in space, while the the parameter marginal is the distribution over the parameters that is constant in time, i.e. stationary. Therefore, the drift can be separated into spatial and parameter components


Similarly the diffusion matrix can be split as


where no cross-correlation for the diffusion of particles is allowed between the coordinate subspace and the parameter subspace, i.e. . If the particles do not undergo the action of the Brownian motion, that is the case of microparticles, also . A convenient choice that constrains the parameter marginal to be constant in time is and . This choice is not unique, for example and gives the same results in terms of mean-and-covariance dynamics, but reformulates the particle ensemble dynamics in term of SDE instead of ODE. Here, we opt for the ODE form of the particle ensemble dynamics.

With the assumptions so far introduced and dropping SDE notation in favor of ODE notation, Eq. (1a) becomes


and the initial condition is


with the initial sample such that


where is a multivariate normal distribution with mean and covariance . Since, the parameter marginal is time-independent, it is immaterial to write in place of and the same holds for the variance .

Figure 1: (Color Online) Simulations results for the prototypical model Eq. (27). Probability density function (a) from ensemble simulations (grayscale) and approximated by mean-and-covariance dynamics (black line ) for four gaussians (colored ellipsis). Spatial marginal (ax) and radius marginal (ar) from ensemble simulations (gray bins) and approximated by mean-and-covariance dynamics (black lines) for four gaussians. Dynamics of the first-order moments (b) and the spatial dispersions (c) for the four kernels (same colors as in panel (a)).

Note that the deterministic process Eq. (19) retains the statistics information about the parameter distribution by including a sample that is distributed according to Eqs. (21). The mean-and-covariance dynamics associated with Eq. (19) can be computed by applying Eq. (13) and considering that and because of the stationariety of the parameter marginal (omitting conditionals)


and the initial condition is given by


The assumption of normally distributed states can be limiting for describing the PDF, therefore it is proposed to approximate the distribution of the generic population in the state space as a superposition of gaussians,


namely a gaussian mixture model, where is a set of gaussians that span the state space, and are weights such that . Note that adopting this representation allows for the introduction of correlations even in the case when , inasmuch when the gaussians span the state space different weights can be assigned to different locations and thus the PDF exhibits a combined dependence on both and . Extending the representation Eq. (24) in the case of multiple populations and considering the initial parameter configuration independent on the initial spatial positions, the parameter marginal for the -th population can be approximated as


where , and are the weight, the means and the covariance associated with the -th gaussian in the -th population (in the following the subscript “” is meant to address the population, while the subscript “” refers to the gaussian). The solution of Eqs. (22) allows then to approximate the spatial marginal for the -th population


Equations (22), the approximation (25), and the reconstruction (26) forms the MCDGM method that applied to acoustphoresis models is used to approximate the dynamics of microparticle populations during acoustophoretic separation and to derive the separation indicators.

Iii Examples

iii.1 Minimal Working Model

In order to illustrate the basic features of the MCDGM method when applied to acoustophoresis, we consider an one-dimensional prototypical model for a single particle population with a radius distribution. Therefore, we can assume that the ensemble dynamics is given by ()


where are the particle positions, and are corresponding the particle radii. The initial conditions and the radius distribution are generated by the weighted superposition of four normally distributed random processes


with spatial averages and variances , and with radius averages and variances . Note that once the initial distribution is generated, Eq. (27) does not retain any information about the four processes that have generated the ensemble. However in the particle ensemble simulations, a tag corresponding to the generating process is applied to the particle to reconstruct the dynamics of the statistics for the single gaussians.

The mean-and-covariance dynamics for the -th gaussian associated with Eq. (27) read as


The numerical solutions of this equations are thus compared with ensemble simulations with total particles generated by Eq. (28) with un-normalized weights , all starting at position with spatial standard deviation . Four average radius were considered with standard deviation , and cross-covariance . Integration of Eq. (27) and Eq. (29) were performed by the Matlab routine ode45 with suitable options as to ensure convergence and accuracy.

Figure 1 reports the comparisons of the direct numerical simulations for the particle ensemble and the mean-and-covariance dynamics for the four gaussians. Panel (a) compares the probability density function in the state space at for the ensemble simulations (grayscale) and for the gaussian mixture: the black line is the isolevel at , while the colored ellipsis are the confidence ellipsis for the four gaussians. Panel (ax) plots the coordinate marginal for the particle-ensemble simulations (gray bins) and that reconstructed by the gaussian mixture (black line). The colors for the gaussians correspond to the colors for the confidence ellipsis in panel (a). Finally, panel (ar) shows the radius-marginal that is stationary in time and corresponds to the radius distribution. As it can see, the MCDGM method provides a good approximation for the four gaussians as well as the PDF. This can be appreciated by the quantitative comparisons in Fig. 1(a)-(b) where the average population positions, namely the first-order moments (panel a), and the spatial dispersion (panel b) for the four gaussians are plotted as function of the time. The symbols correspond to the data extracted from the particle-ensemble simulations by using the particle tags, while the lines are computed by the MCDGM method. In all of the four cases, MCDGM provides a good approximation of the statistics for the marginals of the four gaussians, and as a consequence of the spatial marginal.

iii.2 Buffer-Dependent Separation Performance in Acoustophoresis of Blood Components

Adjustments in the carrier-fluid properties to enhance the separation performance has been successfully employed in acoustophoretic separation involving diluted blood samples Urbansky et al. (2017). The authors, in place of using pure Phosphate Buffer Saline (PBS) in which RBCs/WBCs separation was highly unefficient, employed PBS and Stock Isotonic Percoll (SIP) at different dilution rates () to change the fluid properties and consequently the acoustic contrast factor for both the RBCs and the WBCs. Because of the specific cell properties, as the concentration of SIP is increased the mobility of WBCs decreases as much as that of RBCs, see Fig. 2. About the mobility of WBC and RBC are almost equivalent, and in the correspondence of the acoustic mobility for the WBC population approaches zero and thus the two populations can be successfully separated being the WBCs segregated in the correspondence of the inlet position. In the following, these experiments are simulated by considering a 1D model of separation, radius distribution derived from Coulter Counter measurements, and density and compressibility measurements adapted from Cushing et al. (2017).

Figure 2: (Color Online) Mobility of the RBCs (red) and the WBCs (black) as function of the SIP concentration . Continuous line is the mobility computed for the average radius of the particle populations by Eq. (31), while the areas indicate the ranges . ()

The 1D model equation with radius distribution read as


with . In this equation, where is the acoustic energy density, and is the channel width. The acoustophoretic mobility is given by


where is the particle radius, is the fluid viscosity and the contrast factor


is a function of the particle/fluid compressibility ratio , and the particle/fluid density ratio . The fluid compressibility and density are considered as function of the SIP concentration using the polynomial interpolations described in Urbansky et al. (2017). The radius distributions for the RBCs and the WBCs are given in terms of gaussian mixtures


which are shown in Fig. 3 and for which the caption reports the gaussian mixture parameters and the physical parameters. The initial spatial distributions for the two cell types are


where and , meaning that they have the same starting position and the initial spread.

Figure 3: (Color Online) Radius distribution from Coulter Counter (symbols) and approximated by gaussian-mixtures (black) for RBC and WBC populations. The gaussian-mixture parameters and particle physical parameters are:
Figure 4: (Color Online) Simulations results for the WBC/RBC separation model Eq. (30) with population distribution given in Fig. 3. WBC (left) and RBC (right) -marginal from ensemble simulations (bins) and MCDGM (lines) at different times as the concentration varies: (L1-R1), (L2-R2), (L3-R3), (L4-R4).

The MCDGM equations corresponding to the ensemble Eq. (30) can be derived by Eqs. (22) applying the transformation Eq. (15). The equations obtained are formally identical to Eqs. (29)


where is the mobility and is the derivative of the mobility with respect to the radius calculated, both calculated for the average radius . The integration of Eq. (30) and Eqs. (35) were performed by using the Matlab routine ode45 with suitable parameters to ensure convergence and accuracy.

Figure 4 shows the comparison of spatial marginal resulting from particle-ensemble simulations (bins) and the MCDGM method (lines) for WBC (black) and RBC (red) for different times and at different SIP concentrations. For pure PBS buffer, i.e. , the mobility of WBCs is higher than that of RBCs and the two gaussians corresponding to the largest radii of the WBCs already moved at the channel center-line for . At the same time-instant the RBCs, to which the gaussian for the smallest radius contributes for the largest part, are still located at or better in the range . For , there is not an appreciable difference between the mobility of the WBCs and that of the RBCs (see Fig. 2), and for the gaussians for the two largest WBC radii are located at , while the largest part of RBCs occupy the region . For , there is a dramatic change in mobility for WBCs that now for are located in the range , while the RBCs are in the region . For , the WBCs reaches the isoacoustic concentration so that they remain close to the initial point for times . The RBCs still have an appreciable mobility for this SIP concentration and for they occupy the region . Finally, also in this case one can appreciate the good approximation properties of the MCDGM method when compared with the ensemble simulations for microparticle distributions and mobility values occurring in real-world applications.

iii.3 Free-Flow Acoustophoretic Separation

So far, although the approximation properties of the MCDGM method have been illustrated, the method has not been applied to any model corresponding to a real-world case of acoustophoretic separation, namely free-flow acoustophoretic separations. In order to do this, (i) one needs to introduce a model for the axial flow that takes into account for the hydrodynamics parameters, such as the overall flowrate and the side/center flowrate ratios, and (ii) it is necessary to develop further the MCDGM method to introduce the separation indicators. In this section we investigate the approximation properties of the MCDGM method and study the reliability of the method when some of the separation parameters vary while adopting different prefocusing strategies.

Axial Flow Model.

Let us consider the inlet flowrate ratio


where is the inlet flowrate at sides and is the inlet flowrate at the center, is then the total flowrate. Given the flowrate ratio , it is possible to estimate the position of the streamline separating the side and the center inlet streams (“fj” stands for flow-joining) by assuming


where is the height of the microchannel and is the axial velocity field considered constant for the entire channel length .

Similarly, one can estimate where the particle are separated into the outlet and center streams, by considering the outlet flowrate ratio


or in terms of the position of the streamline separating the side and the center outlet streams (“fs” stands for flow-splitting)


resulting thus


The axial flow can be computed by considering the Poisson problem for the microchannel cross-section


where , and the constant such that the normalization condition


is verified.

Figure 5: (Color Online) Axial velocity profile for a rectangular cross section with , , and as function of the cross-section dimensionless position (a). Gaussian starting positions and corresponding distribution for prefocused particle streams (b) and non-prefocused particle streams (c). Vertical dotted lines are in the correspondence of .

Figure 5(a) shows the axial velocity profile obtained by solving Eq. (41) with the constraint Eq. (42) for the case , , and . The vertical dotted lines as well as the inlet position are computed by setting in Eq. (37). These are typical dimensions, flowrate and flowrate ratios used in acoustophoretic separation experiments Urbansky et al. (2017).

Separation Indicators.

The usual procedure to quantify the performance during acoustophoretic separation experiments is to measure the fraction of particles in the side and the center outlets downstream the separation channel by varying the voltage applied to the transducer. This can be done by collecting the samples and counting the particles with either a Coulter Counter or a FACS machines, here the same quantification method is adopted by developing further the analysis of the MCDGM method.

The fraction of particles belonging to the -th population that flow into the side-stream can be computed by considering the spatial marginal along the width of the channel, i.e. -direction,


and defining the side-stream recovery, henceforth , for the -th population as the associated cumulative (omitting conditionals)


in which the side-stream recovery for the -th gaussian of the -th population is given by


where is the separation time, and is the separation abscissa. These two parameters are constants that we assume depending solely on the flow conditions. Additionally, the separation time depends on the channel length , that for the simulations is .

When the fraction of particle in the center stream is measured, in place of using the SSR the center stream recovery


can be used. For the present example we use exclusively the side-stream recovery.

Model Equations.

For the ensemble simulations of free-flow acoustophoretic separation, we consider the three-dimensional model


that takes into account for the axial flow, acoustophoresis and gravity. The gravitational mobility is


Since the aim is to show how the separation performance depend on the (measured) voltage on the transducer, we adopt a model that is linear with the square-voltage for the energy density in the factor


where the factor should depend on the experimental conditions such as fluid properties, temperature, and generally on the system features. In the present paper the value is fixed .

Here we consider up to four different types of microparticles, . The histograms for and are those shown in Sec. III.2, while and are polystyrene particle with diameters () and (), respectively, with a standard deviation (assuming a single kernel) , the compressibility and the density are given by Cushing et al. (2017).

Figure 6: (Color Online) Side-stream recovery for prefocused streams as function of the voltage for , (left) and (right), different values of , and for different microparticles: beads (blue), beads (orange), (black), and (red). Symbols are computed from ensemble simulations, lines from the MCDGM method applied to Eqs. (35)-(51) (dashed), or Eqs. (35) with correction Eq. (50) (continuous).

The MCDGM method applied to Eq. (47) yields cumbersome equations, so here we restrict the MCDGM analysis to two cases: (i) a corrected plug-flow model for which where


is the average particle velocity between the inlet position at height and the abscissa where the side and center outlet split, and (ii) a 2D model corresponding to disregard the equation for the -component in Eq. (47). We expect that for either moderate/weak buoyant forces (as the overwhelming majority of the cases for polymer microbeads and cells) or fast passages in the separation channel, neglecting the vertical component in Eqs. (47) is a good approximation. For the case (i) the MCDGM equations reduces to Eqs. (35), while for the case (ii) there are additional equations to Eqs. (35)


where is the vertical position of the -th kernel in the -th population.

Prefocusing Strategy.

For prefocused streams at the inlet section of the separation channel since it is expected that the particles focus at one quarter of the channel width and at half-height, it can assume


The initial position in the -th direction is computed by the calculation illustrated in the previous paragraph for , and it is . The initial conditions for the particle ensemble simulations are