Inter-Particle E ﬀ ects with a Large Population in Acoustoﬂuidics

: The ultrasonic manipulation of cells and bioparticles in a large population is a maturing technology. There is an unmet demand for improved theoretical understanding of the particle–particle interactions at a high concentration. In this study, a semi-analytical method combining the Jacobi–Anger expansion and two-dimensional finite element solution of the scattering problem is proposed to calculate the acoustic radiation forces acting on massive compressible particles. Acoustic interactions on arrangements of up to several tens of particles are investigated. The particle radius ranges from the Rayleigh scattering limit ( ka « 1) to the Mie scattering region ( ka ≈ 1). The results show that the oscillatory spatial distribution of the secondary radiation force is related to the relative size of co-existing particles, not the absolute value (for particles with the same radius). In addition, the acoustic interaction is non-transmissible for a group of identical particles. For a large number of equidistant particles arranged along a line, the critical separation distance for the attraction force decreases as the number of particles increases, but eventually plateaus (for 16 particles). The range of attraction for the formed cluster is stabilized when the number of aggregated particles reaches a certain value.


Introduction
Biological or biophysical detection of cells and macromolecules is essential for clinical diagnostics, while physical or chemical processing of particles in fluids plays another important role [1,2]. Portable devices that keep the targets viable in a few cubic millimeters of nutrient medium offer advantages such as high sensitivity, short processing time and reduced sample volume. As ultrasound is inherently compatible to fluid media, the marriage of acoustics and microfluidic dynamics promises the development of acoustofluidic biosensors and processing technologies [3][4][5][6][7]. The acoustic wave is no longer examined as the agent that correlates the measured quantity (e.g., biosensors using surface acoustic wave) but works as an enhancer of functional performance by spatially manipulating cells or biomolecules. Bioparticles with the size from sub-micrometer to a few millimeters are successfully trapped, separated and transported in water-based solutions or body fluids [8][9][10][11][12]. Miniature bulk wave transducers and piezoelectric substrates with deposited interdigital transducers are widely used to couple the solid displacements with microfluidic medium [13][14][15]. The direct effect of incident wave is the acoustic radiation force (ARF) field, which applies a time-averaged force on each particle.
For a standing wave, the generated ARF field is a potential field. Cells and latex with larger density and lower compressibility relative to the acoustic medium are pushed into the pressure nodes (potential minima), whereas ultrasound contrast agents having negative acoustic contrast factor move to the pressure anti-nodes. For example, in the cytometer employing a one-dimensional standing wave, the ARF focuses cells into a narrow stream along a slim pressure node, which usually locates at the centre of a micro-channel [16]. By building up more delicate acoustic potential field in two-dimensional (2D) or three-dimensional (3D) space, the ARF will be able to tune the arrangement of particles [17][18][19]. The underlying physics of the ARF in an acoustic standing, traveling or hybrid wave is readily understood. The linear momentum transfer of ultrasound to a single particle, namely the primary ARF, has been extensively investigated. The influences of shape, size and compressibility of the micro-object, together with the incident wave type and medium viscosity have particularly been considered in previous theoretical analyses [20][21][22][23][24]. Diverse calculation methods, including closed-form analytical solutions, infinite series solutions, finite element method (FEM) and boundary element method, have been proposed [22,[24][25][26][27][28][29]. These seminal studies acknowledge the framework of primary ARF calculation: the surface integral of radiation stress tensor over a certain obstacle, which is only determined by the primary incident wave. Therefore, the major task is to solve the scattering field precisely and efficiently, and it still remains an active research subject with increasingly complex mathematics addressing complex geometry, incident waves, and surrounding boundaries [30][31][32][33][34].
Micro-particles trapped initially at the pressure nodal planes of a standing wave eventually aggregate [17]. Ultrasound aggregation of bioparticle suspensions demonstrates potential for immunoassays of molecular or particulate antigens [2,7]. It is also of help to the perfusion process of cell cultivation, in which the medium is exchanged continuously, while the cell retention in the suspension culture is achieved by using a standing wave [35]. The group behaviours of micro-particles suspended in an acoustic field are more complicated than the individual behaviour, i.e., a single target being driven to the pressure node. The aforementioned primary ARF model captures only the time-averaged mechanical effects originating from the incident wave on a single particle. However, the presence of other particles will produce a series of scattering events, which will generate an additional component of the radiation stress tensor on the current target. The time-average force resulting from the acoustic interactions of multiple particles is referred as the secondary ARF, the acoustic interaction force or the Bjerknes force. Micro-particles clumping together, especially at a pressure node, is attributed to the attraction effect of the secondary ARF.
Unlike the well-established primary ARF theory of solid particle, enrichment of the deep physics for the secondary ARF, either in theoretical models or in fundamental experiments, is still on the way [36,37]. Most of the related studies focus on the acoustic interaction effects of pulsating bubbles, in which the force is treated as a body force, far from the traditional theoretical framework of ARF [38][39][40]. The derived conclusions all stem from the assumption of an incompressible fluid medium. Considering the physical aspect of the incident wave, the employed simplification limits the distance between the particles to be much smaller than the wavelength. However, the operating frequency of the current acoustofluidic devices is usually set at several mega-Hertz, which makes these analyses inapplicable. Employing the widely adopted concept of acoustic radiation stress tensor, Embleton proposed a preliminary analysis of the acoustic inter-particle force between two spheres with the axial line parallel to the wave propagation direction [41]. Zhuk provided a rigorous solution for a pair of rigid spheres suspended in an ideal fluid with plane waves traveling along or perpendicular to the line joining the particle centres [42,43]. This method is an extension of the mathematical treatment for a single target used by King [20]. Doinikov paved the way for implementing acoustic stress tensor in the calculation of acoustic interaction force acting on bubbles [44][45][46]. The bubble dynamics for different pressure amplitudes, particle sizes, and separation distances were analysed, which gave plausible explanations for the formation of bubble grapes and acoustic streamers observed in the corresponding experiments. Later, Doinikov extended the interaction force theory to a general form that allows for particles with different nature (bubble-drop, bubble-solid, and drop-drop) [47][48][49]. Departing from here, the studies of the secondary ARF theory returned to solving the acoustic scattering problem, but this time with multiple obstacles.
More tedious calculations and mathematical difficulties are involved while solving the secondary ARF. Multipole expansions must be performed not only for the primary scattering field of co-existing particles but also for the multiple re-scattering events. For simplicity, most of the early studies took into account only the primary scattering wave from multiple particles, while Zheng and Apfel performed a calculation of the re-scattering field that was accurate up to the second order [50]. Within the Rayleigh scattering limit, the re-scattering contribution to the probe point is dominated by the primary scattering wave (having undergone only a single prior scattering event at the source point) and the re-scattered field can be approximated with a combination of the monopole and dipole scattering. Silva and Bruus proposed a closed-form mean-field theory for the acoustic interaction force of the Rayleigh scatterers in inviscid fluid [51]. For the general case (free of any restriction on the particle size, separation distance or wavelength), the multipole series expansion is essential to solve the multi-scattering problem, and analytically deriving the scattering coefficients from the boundary conditions is a significant challenge for the state-of-the-art mathematics. Accordingly, numerical schemes based on the weighted-residue method have been proposed to determine the coefficients of expansions [23,[52][53][54]. Recently, more accurate algorithms for the secondary ARF on elastic spheres have been proposed [27,28,55,56], in which the fully coupled scattering fields are simultaneously solved by the FEM, spontaneously including the high order scattering and multipole terms. With this numerical method, we successfully explained why particle agglomeration only happens at the pressure node and predicted the dominant range of attractive inter-particle force [55], which agreed well with the measurement.
Even though all the aforementioned calculation methods claim the ability to address the complex case of multi-particle structure, the implementation in specific analyses has only been achieved for the pairwise interaction. The complexity of partial-wave expansion based method increases enormously for more scattering events (particle number and pre-scattering process), while the 3D fully coupled FEM analysis of multiple targets requires considerable computational resources. A sophisticated algorithm for analysis of the secondary ARF of massive co-existing particles is still further off. Limited by the available computational tools, the group behaviours of bioparticles on the pressure node plane could not be fully investigated. For acoustofluidic devices, moving from the proof-of-concept demonstrations to application-driven product designs requires a powerful but convenient tool to understand the particle-to-particle interactions at high particle concentrations.
Here we introduce an efficient hybrid algorithm to model the ARF for a number of co-existing compressible particles. The multi-scattering and re-scattering problems are solved within a simplified FEM framework. By adopting the Jacobi-Anger expansion, the 3D scattering problem can be reduced to a 2D one, which greatly saves time and computation cost. Then the relationship between the surface integral of the acoustic radiation stress tensor and the line integral is analytically derived to figure out the ARF. With a commercial FEM package and the derived analytical equation, we investigate the effects of massive particle interactions between compressible spheres on their group behaviours. Considering that the secondary ARF overweighs the primary one only at the pressure node, we focus on the targets located there.

ARF Acting on Massive Cylindrically Symmetrical Targets
Once the momentum carried by the incident sound wave is changed by the scattering of an obstacle, the ARF arises according to the conservation law. When the obstacle size is much larger than the viscous penetration depth, viscous dissipation of the incident acoustic wave can be neglected. Accordingly, the integral formula for the ARF acting on an arbitrary object in an ideal compressible fluid is applicable to most of the acoustofluidic operations (below the GHz range in water) [23]: where Ω is the unperturbed surface of the object, ρ 0 is the medium density, κ 0 is the compressibility of fluid medium, P 1 and V 1 represent the first-order pressure and velocity of the total sound field, respectively, → n is the outward unit vector normal to surface Ω, and < > denotes the average operator over a time period. Equation (1) is derived from the time-average perturbation theory, which makes it possible to calculate the ARF accurately up to the second-order terms in the Mach number by using the velocity field potential satisfying the linear wave equation. With the first-order quantities derived from a finite element solution of the scattering problems, the time-averaged radiation stress tensor in Equation (1) is obtained. However, applying the surface integral to massive particles not only requires solving the multi-scattering problem efficiently, but also relies on the precise calculation across the fluid-particle interface. When using 3D FEM, the mesh elements on the particle surface should be fine enough to ensure the integral accuracy, which makes the computational cost unaffordable.
The depth of micro-channel employed in various acoustofluidic devices is usually less than 100 µm, and the particle agglomeration only form a single-layer structure. Correspondingly, the particle trajectories in the x-y plane are of interest. In Figure 1a, let n cylindrically symmetric compressible bodies be freely suspended at the same height in a fluid irradiated by an acoustic plane wave. The cylindrically symmetric bodies may have different geometrical configurations. Assume that the plane wave has been excited for a sufficiently long time so that transient processes have died out and the particles oscillate with the acoustic frequency. The shear wave inside the solid body is ignored. Thus, the particles and outside acoustic medium can be treated as ideal compressible fluids. The first-order quantities of the scattering and total fields are directly obtained by solving the discretized Helmholtz equation together with a suitable boundary condition. In the following, we will show that the x-and y-direction ARF on a cylindrically symmetric body can be derived by a summation of 2D FEM calculations. Consider a target locating at the origin of the global Cartesian coordinate system and a plane wave coming from the x-axis (Figure 1b). According to the Jacobi-Anger expansion, the time-independent component of the incident plane wave can be expanded into a series of cylindrical sound waves as follows: (2) where i is the imaginary unit, ε m = 2 − δ 0m , δ 0m is the Kronecker delta, J m (kr) is the m-th Bessel function, k is the wave number, and r and ϕ are cylindrical coordinate variables. Similarly, the time-independent part of the scattered wave field can be expanded as a sum of cylindrical Hankel functions in the following form: where the dimensionless scattering coefficients a m for each partial wave is position dependent, different from the traditional expression of a cylinder. It should be notified that the scattering field of other co-existing particles are presented in their respective coordinates according to Equation (2). By using the translational addition theorem [57], acoustic parameters of the scattering field can be unified to the global coordinate system, where the target particle locates at the origin. With Equations (2) and (3), the scattering of massive cylindrically symmetric targets can be solved in a lower-dimensional space (y-r space) by FEM. The total sound pressure p t of any point on the target surface S can be written as where p t_m is the m-th order component of the total pressure field, P m (y, r) is the m-th order component of the total pressure field derived by 2D FEM, and l is the target profile on the r-y (or x-y) plane. The complex velocity potential is related to the complex pressure. According to Equation (4), we can obtain the total velocity along the radial, tangential and y-axis directions at any point on the target surface, respectively, where v y_m (v r_m or v ϕ_m ) is the m-th order component of the y-axis (radial or tangential) velocity field, V y_m (V r_m or V ϕ_m ) is the m-th order component of the y-axis (radial or tangential) velocity field derived by the 2D scattering simplification, and Substituting Equations (4)-(6) into Equation (2), the time-averaged ARF can be calculated. Theoretically, for the incident plane wave with the wave vector in the x-y plane, the z-axis ARF is always zero. In the next part, we will further simplify the surface integral on a cylindrically symmetric target to a line integral over the target profile on the x-y plane.

Force in the y-Direction
According to Equations (2) and (6), the y-axis component of the ARF has the following form: n r v r v y + n y v 2 y ds.
For convenience, the above equation is divided into three parts to perform the integral calculation. The first part is related to P 2 1 : Due to the orthogonality of trigonometric functions, the corresponding components vanish when m n. Thus, we have For the second part related to V 2 1 , Continuing in the same manner, we consider the line integral of the third part: Combine Equations (9)-(11), the y-axis component of the ARF can be expressed as a series expansion of line integrals: where F rad y,m = − πκ 0 2 l n y r(l)P 2 m (l)dl is the m-th order component of the y-axis ARF.

Force in the x-Direction
Following the procedures of Section 2.2, the x-axis component of the ARF has the following form (Equation (14)). Due to the second-order cross terms in velocity, calculation of the x-axis ARF component is more complicated than the y-axis one, however, the process of transforming the surface integral into the line integral is the same. Here, we skip the details and only give the final expression Equation (15).
2 v 2 r n r cos ϕ + v r v ϕ n r sin ϕ + v r v y n y cos ϕ + v y v ϕ n y sin ϕ ds, F rad where

Method Validation
The commercially available FEM package COMSOL Multiphysics is used to perform the scattering event calculations. 2D model is built in the Pressure acoustics module, so the adopted acoustic parameters in the calculation are density and sound speed. A standing wave field is formed by the superposition of two traveling waves traveling in opposite directions. A sphere with radius a is located at a distance λ/8 to a pressure node. A spherical radiation boundary applied on the peripheries of the acoustic medium (radius of λ) delimits the computational domain, hence, no waves are reflected by the model boundary. The pressure and normal velocity are continuous across the interface between the particle and fluid domain. Free triangular elements are used to mesh the calculation domain. The incident wave is expanded into a series of cylindrical waves by selecting the plane wave expansion node. Then, F rad y,m and F rad x,m are derived by a simple summation. A mesh refinement study is performed, which shows that the integral is accurate enough when the maximum element size of the sphere is smaller than a/120 and at least six elements are maintained in a wavelength of the acoustic medium. The acoustic parameters of silica bead used in the calculation are listed in Table 1. A finite number of terms in the cylindrical wave expansion are used in Equations (13) and (16) to calculate the ARF and the truncation order of the series expansion is m. The relative change in the force caused by increasing the series length is where m n=0 F rad n and m+1 n=0 F rad n are the calculated ARFs with the truncation orders m and m + 1, respectively. Figure 2a,b shows the calculation convergence for a single sphere with various dimensionless radii ka and different angles Φ of incident wave (relative to the y-axis), respectively. It can be seen that the calculation converges rapidly for all cases and five terms are suffice to achieve a desired accuracy of 0.01%. The smaller the value of ka is, the faster the calculation converges. If the wave is incident from the y-direction, using only one term is acceptable.  Table 2 compares the results of our semi-analytical method to both the 3D FEM approach and the analytical theory (Yosioka and Kawasimato [21]). The ARF calculated by the proposed semi-analytical (2D FEM) approach shows good agreement with the 3D FEM method throughout the ka range considered (difference smaller than 1%). Within the Rayleigh scattering limit, the 2D FEM provides quality results with a similar agreement with the analytical one as the 3D FEM calculation. Since the assumption of small scatterer is adopted in Yosioka and Kawasimato's analytical model, the analytical prediction overestimates the force for larger particles. In the 2D FEM, the computational time varies from 20 to 30 s with the particle size, while it exceeds 20 min in the 3D FEM(All solved on a laptop with a 2.4GHz Intel Core i5-9300H CPU with 8G RAM). The significant advantage of our approach is the convenience of modeling the force acting on a target with an arbitrary geometry. We therefore calculated the ARF of a revolving body, which consists of two hemispheres with a radius of 10 µm and a cylinder with a height of 10 µm. In this case, the direction of the incident wave will affect the magnitude of ARF. When the wave vector does not coincide with the symmetric axis of the target, the traditional FEM should be implemented in a 3D domain, taking. Table 3 illustrates the predicted force for various angles between the wave vector and the y-axis. Our 2D FEM approach maintained the same level of accuracy (177,742 degrees of freedom are solved in the 2D FEM, and 3,769,394 degrees of freedom in the 3D FEM), while the computational time was reduced to 20 s (more than one hour in the 3D FEM), and hence provided a powerful mathematical tool for analyzing the acoustic interactions of massive targets. For two spheres separated by a different distance along z-axis and placed at the center of a pressure node and antinode, the secondary ARF acting on each particle has the same magnitude but opposite direction due to the symmetric configuration. The normalized second-order ARF acting on one target is shown in Figure 3. The difference between the 3D and 2D FEM calculations is negligible except for the critical separation distance (near l = 0.5λ), in which case the repulsive force will change to attraction with a small magnitude.

Acoustic Interactions Between Massive Particles
Using the proposed hybrid computational scheme, a series of interaction cases between compressible spheres are hereby investigated, which are closely related to the well-known experimental observations. As shown in Figure 4a, suppose a pair of particles is located in a standing wave field at a distance d to the pressure node and separated from each other by a distance l. Our previous study found that the secondary ARF acting on a pair of particles located at the pressure node (d = 0) becomes attractive when the relative distance l reaches a critical value [53]. Here, we show the phase diagram of the secondary ARF (along center lines of the targets) as a function of the spatial parameters l and d (Figure 4b). The calculation shows that ka does not significantly affect the value of the critical distance (where the ARF changes direction). In the range spanning from the Rayleigh scatter limit (ka << 1) to the Mie scattering region (ka ≈ 1), the critical distance changes only by 0.01λ. The critical distance at the antinode is the smallest, about 0.44λ, while the maximum value, about 0.62λ, occurs at the pressure node. When a pair of particles is approaching the pressure node, the effective range of the attractive secondary ARF expands, which provides a reasonable explanation for the experimental observations that particle agglomeration occurs near the pressure node [53]. In the following analysis, we focus on the acoustic interactions at the pressure node. The acoustic interaction of the two spheres is always symmetric, inducing an equal in value but opposite in direction secondary ARF. The asymmetric acoustic interactions are then analyzed to find the general characteristics of the secondary ARF. For three identical particles with different relative distances l 1 and l 2 (Figure 5a), Figure 5b,c shows the normalized y-axis secondary ARF acting on the upper and middle particle. The relative position of the lower particle does not significantly change the force on the upper one, but the middle particle is significantly affected. Both re-scattered waves from the upper and lower particles contribute to the secondary ARF of the middle one. It can be concluded that the effective secondary ARF acting on a target depends mainly on its direct neighbors, i.e., the acoustic interaction is non-transmissible and the indirect re-scattered waves from distant particles do not contribute much to the secondary ARF and can be neglected. When the particles are equidistant (l 1 = l 2 ), the middle particle is balanced, and the critical distance for particle aggregation around the center one is about 0.66λ, and a slight change to this value is made by varying ka (Figure 5d). Since the relative position of the lower particle does not significantly change the force on the upper one (Figure 5b), the relative movement between the upper and middle particles resembles to the case of two particles where one is fixed. According to Silva and Bruce's analytical model [49], the interaction force of Rayleigh scatters (ka << 1) only depends on distance between the source and the probe. Our numerical results show that the constraint of ka << 1 can be extended to ka < 2 for Rayleigh limit, which applies to the majority targets in acoustofuidics. Figure 5e shows the phase diagram of the secondary ARF acting on the upper sphere for various l 1 and l 2 . The critical distance for the upper sphere fluctuates slightly with l 2 . As for the interaction between two particles, ka has little effect on the critical distance for the interactions of three identical particles. Figure 5f shows the change of the normalized secondary ARF acting on the upper sphere when the size of lower particle varies. The particles are equidistant, and the dimensionless radii of the upper and middle targets are kept at ka = 0.17. Increasing ka of the lower sphere from 0.085 to 0.340, changes of the secondary ARF on the upper sphere increases from one to three times within a wavelength separation distance. Within the first attraction range, the secondary ARF acting on the upper sphere strengthens for the lower particle with larger ka, which is different from the effects of ka on particles with the same size. We may conclude that the absolute size of co-existing particles (same radius) does not affect the law of acoustic interactions, but the relative size of co-existing particles dominates the group behavior. We further increase the number of co-existing particles to four and investigate the acoustic interactions for the case of identical particle size and relative distance. The only variable is the distance l between adjacent particles. Due to the symmetric configuration, the secondary ARFs acting on the two internal targets have the same amplitude but opposite directions, so do the ARFs of the outer two particles. Figure 6a shows that the secondary ARF acting on the internal particles changes direction three times (l = 0.31λ, 0.59λ and 0.86λ, respectively) within a wavelength separation distance. However, the force of external particles changes direction only at l = 0.63λ (Figure 6b). In addition, the secondary forces on the external particles are almost an order of magnitude larger than those on the internal particles. According to the distributions of secondary ARF, the group behaviors of four particles are more complicated. When the relative distance is less than 0.31λ, the particles tend to gather at the center, whereas for 0.31λ < l ≤ 0.59 the upper two particles and the lower two particles will agglomerate with each other, forming two clusters. When 0.59λ < l ≤ 0.63λ, the particles would theoretically gather at the center, however, the force is then relatively small, which may not induce obvious particle movement. For l > 0.63, the external particles will be repelled, and no cluster will form. In real-life ultrasonic trapping and aggregation applications, the number of co-existing targets exceeds significantly those analyzed above (i.e., three or four particles). For the case where many identical particles are uniformly distributed, it is found that the closer the particle is to the center, the smaller the distance it needs to be attracted. Figure 7 shows the phase diagram of the innermost particles. Below the boundary line, all the particles will aggregate at the center. For larger number of particles, the critical distance for forming a single cluster decreases monotonically and reaches a steady value of 0.07λ, when the number exceeds 16 at ka = 0.170. Increasing ka from 0.085 to 0.340 has a negligible effect on the critical distance. This finding explains the experiment well, in which clusters formed with a smaller relative distance than for a pair of particles [53].
Dispersed particles do not aggregate at the same time. In experiments, particles with smaller distance aggregate first, and then attract other particles to form a larger cluster. Therefore, we investigate the acoustic interactions between a cluster and a single particle, as illustrated in Figure 8. For a certain value of ka, the attractive range of the formed cluster decreases with the increase in cluster size. The critical distance exhibits slight periodic fluctuation after the particle number reaches a certain value (for ka = 0.170, 0.127 and 0.085, n = 15, 20 and 27, respectively). Thus, the range of attraction stabilizes when the cluster has grown to a certain size.

Conclusions
We presented a hybrid scheme for calculating the ARF acting on massive co-existing particles in a standing wave field. By employing the partial-wave expansions of spherical functions, the integral of the acoustic stress tensor acting on the target surface can be analytically derived by a summation of line integrals over the particle profile, and the related scattered field is numerically calculated by 2D FEM. The proposed method is efficient and convenient for investigating the group behavior of targets with a large population. With the proposed method, we analyzed a series of interaction cases between compressible spheres, which are closely related to the observed trapping or agglomeration phenomena in a non-viscous standing wave. The size of studied targets ranged from the Rayleigh scattering limit (ka << 1) to the Mie scattering region (ka ≈ 1). From the calculations, several significant principles of acoustic interactions have been revealed. (1) Acoustic interaction varies with the particle population size but stabilizes after a certain population is reached (smaller than 20 for the line alignment). (2) For particles with the same size, the effective secondary ARF acting on a target is mainly related to its direct neighbors and the acoustic interaction is non-transmissible and the indirect re-scattered waves from distant particles do not contribute significantly to the secondary ARF. (3) For particles with the same size, the value of dimensionless radius ka has little influence on the normalized secondary ARF. (4) The relative radius of co-existing particles greatly affects the acoustic interactions. This study is the first detailed analysis of the massive acoustic interactions (number of particles up to several tens), and the results will contribute to the understanding of complex group behaviors of massive particles in standing-wave interactions. The proposed hybrid calculation method will provide guidance for designing particle separation or agglomeration apparatuses using secondary ARF.
Author Contributions: Conceptualization, K.J. and J.C.; methodology, Y.W.; software, Y.W.; validation, Y.W. and L.L.; writing-original draft preparation, K.J.; writing-review and editing, J.C. and K.Y. All authors have read and agreed to the published version of the manuscript.