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





[corresponding author:]fabio.garofalo@bme.lth.se
Department of Biomedical Engineering, Lund University,
Ole Römers Väg 3 S-22363, Lund, Sweden









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.


Quantifying Acoustophoretic Separation of Microparticle Populations by Mean-and-Covariance Dynamics for Gaussians in Mixture Models
    Fabio Garofalo
    December 30, 2023
===================================================================================================================================





§ INTRODUCTION

Acoustofluidics is a microfluidic technology that using acoustic waves is able to perform separation of microbeads and viable manipulation of cells <cit.>.
Indeed, by exploiting the interaction between acoustic pressure waves and a carrier-fluid suspension of microbeads/cells at microscale <cit.>, acoustophoresis <cit.> is able to trap <cit.>, wash <cit.>, concentrate <cit.>, align <cit.> and separate the suspended microparticles <cit.>.
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 <cit.>.

The microparticle physical parameters that appear in the acoustophoretic force expression <cit.> 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 <cit.>.
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 <cit.>.
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 <cit.>, 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 <cit.>, and recently for the quantification of dispersion in acoustophoretic models <cit.>.
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 <cit.>, 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 <cit.>.
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 <ref> (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 <ref> (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 <ref> 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 <cit.>, 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.



§ THEORY



 §.§ Mean-and-Covariance Dynamics

Let us consider the nonlinear stochastic differential equation in the Itô sense <cit.>


    dQ(t)   = f(Q,t) dt+√(2) σ(Q,t)·dW(t) ,
    Q(t_0^)   = Q_0 ,


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


    E[W(τ)]   = W̅= 0 ,
    cov[W(τ)]   = τ I ,


into the displacement dQ for the states Q.
In equation (<ref>), I 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 (<ref>) represents the initial condition for the realizations of Eq. (<ref>) in terms of the realizations Q_0 that is distributed according to a probability density function .

Equation (<ref>) corresponds to the (forward) Fokker-Planck equation for the probability density ρ(,t |0 t_0) <cit.>

    ∂^_t ρ(,t | 0,t^_0)=ℒ^_FP(,t) ρ(,t | 0,t_0) ,

conditioned for  by

    ρ(,t_0 | 0,t_0)=∫_ℚδ(-0) ρ(0,t_0) d0 .

In equation (<ref>) the Fokker-Planck forward operator (assuming Einstein notation)

    ℒ^_FP(,t)g=
    -∂^_h[f^h_(,t) g]+∂^2_hk[ε^hk_(,t) g] ,

includes the drift f(,t), and the diffusion matrix

    ε(,t)=σ^T_(,t) σ^_(,t) .

This relation can be used to derive the diffusion contribution to the Fokker-Planck operator ℒ_FP when the Itô process Eq. (<ref>) is known, as well as to construct the Itô process when the Fokker-Planck operator is given <cit.>.
The second derivation is performed by computing σ as the Choleski decomposition of the diffusion matrix ε, which is indeed defined by Eq. (<ref>).

Alongside the Fokker-Planck forward operator defined in Eq. (<ref>) is possible to introduce the backward operator <cit.>

    ℒ_FP^bw(,t)=f^h_(,t)∂^_h+ε^hk_(,t)∂^2_hk .

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

The dynamics of the first-order moment m is derived from Eq. (<ref>) multiplying by q^l, integrating over the  ℚ and using the definition of the backward operator ( and ), obtaining

    ṁ^l_(t | 0,t_0)=⟨ℒ_FP^bw(,t)q^l_ | 0,t_0| ⟩,

where  is meant the expectation of g at time t for a distribution that at time t_0 “occupied” the states 0, or conditioned to Eq. (<ref>). Noting that  and ,  equation (<ref>) can be rewritten in terms of the drift

    ṁ^h_(t | 0,t_0)=⟨f^h_(,t) | 0,t_0| ⟩,

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

    ṡ^hk_(t | 0,t_0)   = ⟨f^h_(,t)(q^k-m^k) | 0,t_0|+⟩
       + ⟨f^k_(,t)(q^h-m^h) | 0,t_0|+⟩
       + 2⟨ε^hk_(,t) | 0,t_0| ⟩.

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


    f^h(,t)   = f^h(m,t)+∂_l f^h(q^l_-m^l_) ,
    ε^hk(,t)   = ε^hk(m,t)+∂_l ε^hk(q^l_-m^l_) ,


and substituting these expansions into Eqs. (<ref>) and (<ref>) it has


    ṁ^h_(t | 0,t_0)   = f^h_[m(t | 0,t_0),t] ,
    ṡ^hk_(t | 0,t_0)   = ∂^_l f^h_[m(t | 0,t_0),t]) s^lk_(t | 0,t_0)+
       + ∂^_l f^k_[m(t | 0,t_0),t]) s^hl_(t | 0,t_0)+
       + 2 ε^hk_[m(t | 0,t_0),t] ,


that using vector analysis notation, read as


    ṁ(t | 0,t_0)   = f[m(t | 0,t_0^),t] ,
    ṡ   = J[m(t | 0,t_0^),t]·s(t | 0,t_0)+
       + s(t | 0,t_0)·J^T_[m(t | 0,t_0^),t]+
       +  2 ε[m(t | 0,t_0^),t] ,


where J(m,t)=∂_qf|_m,t. The initial condition for this set of ODE is


    m(t_0^ | 0,t_0^)   = E[Q_0^] ,
    s(t_0^ | 0,t_0^)   = cov[Q_0^,Q_0^] .


Equations (<ref>), or equivalently Eqs. (<ref>), represent the set of differential equations here briefly addressed as , while Eqs. (<ref>) 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

    ∂_q^f|_m,t^=∇_x,gf^x_·([    I    0;    0 ∂_pg ]) ,

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.



 §.§ 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


    ρ^x_ (x,t | 0,t_0)   = ∫_ℙρ(,t | 0,t_0) dp ,
    ρ^ p_(p,t | 0,t_0)   = ∫_𝕏ρ(,t | 0,t_0) dx ,


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

    f(,t)=f(x,p,t)=[f^x_(x,p,t),f^p_(x,p,t)]_^T .

Similarly the diffusion matrix can be split as

    ε(,t)=ε(x,p,t)=([ ε^xx_(x,p,t)            0;            0 ε^pp_(x,p,t) ])

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 ε^xx=0.
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. (<ref>) becomes


    Ẋ(t)   = f^x_(X,P,t) ,
    Ṗ(t)   = 0 ,


and the initial condition is


    X(t_0)   = X_0 ,
    P(t_0)   = P_0 ,


with the initial sample such that


    X_0^   ∼ 𝒩( · | m^x_0,s^xx_0) ,
    P_0^   ∼ 𝒩( · | m^p_0,s^pp_0) ,


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

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


    ṁ^x_   = f^x_(m^x_,m^p_0,t) ,
    ṡ^xp_   = ∂_xf^x_·s^xp_+∂_pf^x_·s^pp_0 ,
    ṡ^xx_   = 2 ∂_xf^x_·s^xx_+2 ∂_pf^x_·s^xp_ ,


and the initial condition is given by


    m^x_(t_0 | 0,t_0^)   = m^x_0 ,
    s^xp_(t_0 | 0,t_0^)   = 0 ,
    s^xx_(t_0 | 0,t_0^)   = s^xx_0 .




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,


    ρ_^(,t | 0,t_0^)=∑_k∈𝒦_^w_k^ 𝒩[ | m^_k(t | 0,t_0^),s^_k(t | 0,t_0^)] ,


namely a gaussian mixture model, where 𝒦 is a set of gaussians that span the state space, and w_k^ 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 x and p.
Extending the representation Eq. (<ref>) in the case of multiple populations and considering the initial parameter configuration independent on the initial spatial positions, the parameter marginal for the h-th population can be approximated as

    ρ^p_h(p)=∑_k∈𝒦_hw_hk^ 𝒩(p | m^p_hk,s^pp_hk) ,

where w_hk, m^p_hk and s^pp_hk are the weight, the means and the covariance associated with the k-th gaussian in the h-th population (in the following the subscript “h” is meant to address the population, while the subscript “k” refers to the gaussian).
The solution of Eqs. (<ref>) allows then to approximate the spatial marginal for the h-th population

    ρ^x_h(x,t | 0,t_0^)=∑_k∈𝒦_hw_hk 𝒩[x | m^x_hk(t | 0,t_0^),s^xx_hk(t | 0,t_0^)] .




Equations (<ref>), the approximation (<ref>), and the reconstruction (<ref>) 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.
























































§ EXAMPLES




 §.§ 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  ()

    Ẋ(t)=R_^2sin[2π X(t)] ,

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


    X_k^(0)   ∼ 𝒩( · | m^x_0,k,s^xx_0,k) ,
    
    R_k^   ∼ 𝒩( · | m^r_k,s^rr_k) ,


with spatial averages m_0,k^x and variances s_0,k^xx, and with radius averages m_k^r and variances s_k^rr.
Note that once the initial distribution is generated, Eq. (<ref>) 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 k-th gaussian associated with Eq. (<ref>) read as


    ṁ^x_k   = (m_k^r)^2 sin(2π m^x_k) ,
    ṡ^xr_k   = 2π (m_k^r)^2 cos(2π m^x_k) s^xr_k+
       + 2 m_k^r sin(2π m^x_k) s^rr_k ,
    ṡ^xx_k   = 4π (m_k^r)^2 cos(2π m^x_k) s^xx_k+
       +4 m_k^r sin(2π m^x_k) s^xr_k .


The numerical solutions of this equations are thus compared with ensemble simulations with 10^5 total particles generated by Eq. (<ref>) 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. (<ref>) and Eq. (<ref>) were performed by the Matlab routine  with suitable options as to ensure convergence and accuracy.

Figure <ref> 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  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.



 §.§ 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 <cit.>.
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. <ref>.
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 <cit.>.



The 1D model equation with radius distribution read as 

    Ẏ_h^(t)=ϵ μ(R_h^,κ̃_h^,ρ̃_h^,η_f^) sin[2π Y_h^(t)/W] ,

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

    μ(r,κ̃,ρ̃,η_f^)=
    2 Φ(κ̃,ρ̃) r^2_/9 η_f^ ,

where r is the particle radius, η_f^ is the fluid viscosity and the contrast factor

    Φ(κ̃,ρ̃)=(1-κ̃)+3/22 (ρ̃-1)/2 ρ̃+1 ,

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 <cit.>.
The radius distributions for the RBCs and the WBCs are given in terms of gaussian mixtures

    R_hk^∼ 𝒩( · | m^r_hk,s^rr_hk) ,   h=RBC ,WBC ,

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

    Y_hk^(0)∼ 𝒩( · | m^y_0,hk,s^yy_0,hk) ,   h=RBC ,WBC ,

where m^y_0,hk=0.05 W and s^xx_0,hk=(5× 10^-3 W)^2_, meaning that they have the same starting position and the initial spread.





The MCDGM equations corresponding to the ensemble Eq. (<ref>) can be derived by Eqs. (<ref>) applying the transformation Eq. (<ref>). The equations obtained are formally identical to Eqs. (<ref>)


    ṁ^y_hk   = ϵ μ_hk^ sin(2π m^y_hk/W) ,
    ṡ^yr_hk   = ϵ μ_hk^ 2π/W cos(2π m^y_hk/W) s^yr_hk+
       +ϵ μ'_hk sin(2π m^y_hk/W) s^rr_hk ,
    ṡ^yy_hk   = 2 ϵ μ_hk^2π/Wcos(2π m^y_hk/W) s^yy_hk+
       +2 ϵ μ'_hk sin(2π m^y_hk/W) s^yr_hk .


where  is the mobility and  is the derivative of the mobility with respect to the radius calculated, both calculated for the average radius m_hk^r.
The integration of Eq. (<ref>) and Eqs. (<ref>) were performed by using the Matlab routine   with suitable parameters to ensure convergence and accuracy.

Figure <ref> shows the comparison of spatial marginal ρ^y_ 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 y/W≃ 0.4 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. <ref>), 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.



 §.§ 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

    q_in^=Q_s,in^/Q_c,in^ ,

where Q_s,in^ is the inlet flowrate at sides and Q_c,in^ 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

    q_in^(y_fj^)=2∫_0^y_fj^∫_0^H u_ax(y,z) dydz/Q-2∫_0^y_fj^∫_0^Hu_ax(y,z) dydz ,

where H is the height of the microchannel and u_ax(y,z) is the axial velocity field considered constant for the entire channel length L.

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

    q_out^=Q_s,out^/Q_c,out^ ,

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

    q_out^(y_fs^)=2∫_0^y_fs^∫_0^H u_ax(y,z) dydz/Q-2∫_0^y_fs^∫_0^Hu_ax(y,z) dydz .

resulting thus

    y_sep^=y_fs^.

The axial flow u_ax^(y,z) can be computed by considering the Poisson problem for the microchannel cross-section

    ∇_⊥^2 u_ax^=const ,     (y,z)∈[0,W]×[0,H]

where , and the constant such that the normalization condition

    ∫_0^H∫_0^W u_ax^(y,z) dydz=Q ,

is verified.



Figure <ref>(a) shows the axial velocity profile obtained by solving Eq. (<ref>) with the constraint Eq. (<ref>) for the case , , and . The vertical dotted lines as well as the inlet position are computed by setting  in Eq. (<ref>). These are typical dimensions, flowrate and flowrate ratios used in acoustophoretic separation experiments <cit.>.



  
*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 h-th population that flow into the side-stream can be computed by considering the spatial marginal along the width of the channel, i.e. y-direction,

    ρ^y_h(y,t | 0,t_0^)=∑_k∈𝒦_h^w_hk 𝒩[y | m^y_hk(t | 0,t_0^),s^yy_hk(t | 0,t_0^)] ,

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

    SSR_h=∑_k∈𝒦_h^ w_hk SSR_hk^ ,

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

    SSR_hk^=1/2 erfc[
    m^y_hk(t_sep^ | 0,t_0^)-y^_sep/√(2 s^yy_hk(t | 0,t_0^))] ,

where t_sep^ is the separation time, and y_sep^=y_fs^ 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 L, that for the simulations is L=4.3 cm.

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

    CSR=1-SSR ,

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


    Ẋ_h^(t)   = u_ax^[Y_h^(t),Z_h^(t)] ,
    Ẏ_h^(t)   = ϵ μ(R_h^,κ̃_h^,ρ̃_h^,η_f^) sin[2π Y_h^(t)/W] ,
    Ż_h^(t)   = -μ_g(R_h^,ρ_h^,ρ_f^,η_f^) ,


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

    μ_g^(r,ρ_^,ρ_f,η_f^)=2 g r^2/9 η_f^(ρ_^-ρ_f^) .

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 

    E_ac^=α V_pp^2 ,

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 α=8.364.

Here we consider up to four different types of microparticles, h=PS5, PS7, WBC, RBC.
The histograms for  and  are those shown in Sec. <ref>, while  and  are polystyrene particle with diameters  () and  (), respectively, with a standard deviation (assuming a single kernel) σ^r_h=7% m^r_h, the compressibility  and the density  are given by <cit.>.



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

    V=1/y_fs^-y_in^∫_y_in^^y_fs^u(y',z_in^) dy ,

is the average particle velocity between the inlet position y_in at height z_in^ and the abscissa where the side and center outlet split,
and (ii) a 2D model corresponding to disregard the equation for the z-component in Eq. (<ref>).
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. (<ref>) is a good approximation.
For the case (i) the MCDGM equations reduces to Eqs. (<ref>), while for the case (ii) there are additional equations to Eqs. (<ref>)


    ṁ^x_hk   = u_ax^(m^y_hk,z^0_hk) 
    ṡ^xx_hk   = 2 ∂_y^u_ax^(m^y_hk,z^0_hk) s^xy_hk ,
    ṡ^xy_hk   = ϵ μ_hk^2π/Wcos(2π m^y_hk/W)+
       + ϵ μ'_hk sin(2π m^y_hk/W) s^xr_hk+
       + ∂_y^u_ax^(m^y_hk,z^0_hk) s^yy_hk ,
    ṡ^xr_hk   = ∂_y^u_ax^(m^y_hk,z^0_hk) s^yr_hk ,


where z^0_hk is the vertical position of the k-th kernel in the h-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

    m^y_0,hk=1/2 y_fj^ ,     m^z_0,hk=1/2H .

The initial position in the y-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


    X_hk^(0)   ∼ δ(·) ,
    
    Y_hk^(0)   ∼ 𝒩( · ;m^y_0,hk,s^yy_0,hk) ,
    
    Z_hk^(0)   ∼ 𝒩( · ;m^z_0,hk,s^zz_0,hk) ,


where δ is the Dirac-delta distribution,  and  are the average positions given above, and  with  the variances.
The initial distribution is shown in Fig. <ref>(b).

Figures <ref> show the side-stream recovery as function of the applied voltage when , two SIP concentrations and three cases of  for four different types of microparticles, two polymer microbeads PS5 and PS7 and two cells WBC and RBC.
The particle ensemble simulations are indicated with the symbols, the MCDGM method applied to the 1D model is indicated by the solid lines, while it is indicated with dashed lines for the 2D model.
In all the cases both the 1D and 2D models can approximate the numerical data quite well, showing the approximation properties of the MCDGM method, its robustness in terms of the parameter variations, and the validity of the effective velocity assumption Eq. (<ref>).



  
*No-Prefocusing Strategy.
For non-prefocused particle streams the particle ensemble simulations are initializated by the conditions


    X_hk^(0)   ∼ δ( · ) ,
    
    Y_hk^(0)   ∼ 𝒰( · ;m^r_hk,y_fs^-m^r_hk) ,
    
    Z_hk^(0)   ∼ 𝒰( · ;m^r_hk,H-m^r_hk) ,


where  is a uniform distribution between a and b. Note that a small portion of the cross-section has been excluded from the particle distribution, and this corresponds to the fact that the particles cannot have a distant from the walls less than the average radius.
For that regards the starting position of the gaussians chosen a resolution n_y^ in the y-direction, the resolution in the z-direction is , it has for the average positions


    m^y_hk   = h Δ y ,     h=1 ,... ,n_y-1^
    
    m^z_hk   = k Δ z ,     k=1 ,... ,n_z-1^


where  and , while for the variances

    σ^y_hk=Δ y log 2 ,    σ^z_hk=Δ z log 2 ,

where  is a factor chosen as to accommodate for the smootheness of the spatial distribution.
The initial distribution is shown in Fig. <ref>(c).



Figures <ref> shows the side-stream recovery as function of the applied voltage when 
for non-prefocused particle streams for two different polymer microbeads PS5 and PS7.
The particle ensemble simulations are indicated with the symbols, the MCDGM method applied to the 1D model is indicated by the solid lines.
Also for the case of non-prefocused particle streams the MCDGM method in both the 1D and 2D versions can approximate the numerical data quite well.
This confirms again the relialbility of the MCDGM method when applied to free-flow acoustophoretic separations.





 §.§ Inferring Size Histograms

Let us suppose one wants to determine the weights for the gaussians that span the parameter space by performing separation experiments at different voltages.
The theoretical value for the side-stream recovery parametrized with respect to the unknown weights w^p_ can be written as

    SSR(V; w^p_)=∑_k∈𝒦w^ p_k SSR_k^(V) ,

where V is the voltage. When a set of measurements SSR^exp_l as function of the voltage applied on the transducer V_l is available, the set of unknown weights can be determine by requesting that the distance between the measurements and the values provided by MCDGM is minimum

    ∑_l[∑_k∈𝒦F_kl^ w^ p_k-SSR^exp_l]^2=min ,

where the matrix F_kl^ is

    F_kl^=SSR_k^(V_l^) ,

and it can be viewed as the transfer function for the k-th gaussian when a voltage V_l^ is applied.
The problem Eq. (<ref>) can be reformulated as rectangular linear equation

    ∑_k∈𝒦F_kl^ w^ p_k≈SSR^exp_l ,

and the solution is given by

    w^p_=F^+_SSR^exp_ ,

where  is the Moore-Penrose pseudodinverse. Alternatively, the Matlab routine  or the Mathematica routine  perform the same calculations by providing the matrix F and the measurements SSR^exp_.

Figure <ref> shows the results obtained for the size histogram estimation of WBCs from synthetic SSR experiments generated by particle ensemble simulations at different SIP% with α dependent on the SIP%: , , and . The different concentrations were chosen to enlarge the dynamic voltage range, both the full range and the dynamic range were used to determine the histograms.
The gaussians where chosen so that 𝒦=25 in a range r=2-8  m, namely


    m^r_k   = r_min^+k Δ r ,     k=0...𝒦 ,
    σ^r_k   = Δ r log 2 .


with  and Δ r=0.25  m.
As it can see from the figure, with the aid of MCDGM method is possible to estimate the radius distribution from the SSR measurements as function of the voltage. The discrepancies between the histogram used to generate the SSR measurements and that computed from the MCDGM method applied to the dynamic range are small. Using the full range seems to provide subsampled histograms. An exception seems to occur for  where the dynamic range gives an additional “bump” for smaller radii.



§ DISCUSSION

The MCDGM method have been applied to a variety of situations. In all of these case the method has proven its reliability and robustsness in terms of varying the simulation parameters, providing a good approximation of the spatial marginals, the prediction of the SSR, and in the estimation of the histograms.

With respect to the particle ensemble simulations the MCDGM method can generate results much more faster.
The relative computational costs have been estimated by assuming as reference the 1D-MCDGM equations and resulted to be
7-10× for the 2D-MCDGM equations, and 10k-15k× for the particle ensemble simulations.
The exact cost depend on the number of gaussians used in the MCDGM method or the number of particles in the ensemble simulations, the dimensionality of the problem, i.e. 1D, 2D or 3D. It is however possible to claim that the speedup of the MCDGM method over the particle ensemble simulations is about three order of magnitude.

The advantage of this speedup can seem not beneficial for the cases presented in Sec. <ref>(A)-(B) for which the analytical solution of the lateral movement can be applied to particle ensembles <cit.>.
It's however remarkable that when the MCDGM method is applied to the simple models in Sec. <ref>(A)-(B), analytical solutions similar to that presented in <cit.> for the dispersion problem are available.
Therefore, the analysis of the computational speedup should been performed based on the analytical solutions, but the results of the comparison are trivial since the MCDGM method can capture the behavior of the PDF with just a few gaussians.
The advantage of using the MCDGM method is however obvious when applied to the models investigated in Sec. <ref>(C)-(D),
especially in the case of histogram estimations where a significant number of gaussians must be used to have a fine sampling of the parameter space.





§ CONCLUSIONS

A method for quantifying the acoustophoretic separation of microparticle populations with continously-distributed parameters has been introduced.
The method has been applied to an one-dimensional abstract model of acoustophoretic separation, where the particles were considerd to have a radius distribution.
The approximation property of the method has been illustrated by comparing the statistics for the particle ensemble simulations and those computed by using the MCDGM method.
The application of the method to a model related with previously published experiments of WBC/RBC separation has shown its robustness with respect to distributions and changes in the fluid properties occurring in real-world applications.
Furthermore, the method has been employed to quantify free-flow acoustophoretic separation performance with and without prefocusing of the particle streams, and for the estimation of size histograms from separation performance data.
In all of the cases here investigated the application of the MCDGM method to the model equations has shown very good results in terms of approximation of the numerical data from particle-ensemble simulations and in the estimation of the size histograms.

For that regards future comparisons with experimental data, the MCDGM method promises undisputed advantages for the experimentalists in terms of analysis of the experimental outcomes.
Firstly, it is model-free, meaning that it is possible to increase the complexity of the physical model to obtain a more refined representation and a better consistency with the experimental data.
The complexity can reach the level of numerically synthesized velocity fields, acoustic fields, and precomputed scattering laws on particles of arbitrary shape, while the applicability of the MCDGM method is still guaranteed.
Secondly, it can be adapted to perform both hydrodynamics and acoustics calibration of acoustophoretic devices, so that the model inputs can be actually measured instead of being derived from approximate theoretical estimates and guesses such as those used in this manuscript.

The introduction of the MCDGM method and its application to free-flow acoustophoresis represents a breakthrough for the assessment of the separation performance in acoustophoretic device.
It possesses unprecedented features such as incorporating and estimating parameter- and spatial-distributions, very low computational cost compared to particle ensemble simulations, and effective/practical dimensional reduction.
This means that the numerical implementations of this method are suitable to be executed on single-board computers, enabling thus for ultra-compact applications which embed control, calibration, and analysis algorithms on the same processing unit.















This work was supported by the Knut and Alice Wallenberg
Foundation (Grant No. KAW 2012.0023).


























































