An acoustic spheroid-on-chip platform for long-term culturing of 3D cell aggregates

: Microfluidic lab-on-chip devices are widely being developed for chemical and biological studies. One of the most commonly used types of these chips is perfusion microwells for culturing multicellular spheroids. The main challenge in such systems is the formation of substantial necrotic and hypoxic zones within the cultured spheroids. Herein, we propose a novel acoustofluidic integrated platform to tackle this bottleneck problem. We show that such an approach enhances cell viability and shrinks necrotic and hypoxic zones in these spheroid-on-a-chip platforms without the need to increase the flow rate, leading to a significant reduction in costly reagents' consumption. Proof-of-concept, designing procedures, and finite element numerical simulation are discussed in details. Also, the effects of acoustic and hydrodynamic parameters on the cultured cells are investigated. The results show that by increasing acoustic boundary displacement amplitude ( 𝑑 (cid:2868) ), the spheroid’s proliferating zone enlarges greatly. Moreover, it is shown that by implementing 𝑑 (cid:2868) =0.5 nm, the required flow rate to maintain the necrotic zone below 13% will be decreased 12 times compared to non-acoustic chips.


Introduction
Cancer is the second leading cause of death globally, accounting for more than 8 million deaths per year [1]. One of the first crucial steps in the battle against this disease is understanding the underlying processes happening during tumorigenesis and finding the different affecting factors [2]. A popular method to mimic the tumor's microenvironment and study its behaviors is to employ tumor spheroids [3]. Spheroids are three-dimensional (3D) cellular masses formed from the aggregation of cells under special culturing conditions [4]. The morphology, growth, and cell interactions in spheroids are similar to those in actual tumors, making spheroid culture both an appropriate in vitro model of tumor and a valuable tool for testing drugs and treatment efficiency [5,6].
For many years, methods such as hanging drops and 96-well plates have been used to form and study tumor spheroids [7]. However, the static condition in these methods leads to accelerated depletion of nutrients and accumulation of waste, adversely affecting spheroid growth and leading to false results on drug candidates [8]. and culture these spheroids on microfluidic chips is the integration of microchannels and microwells [13,17].
The cell suspension is first introduced at the inlet for a spheroid generation, reaching the microwells through interconnected microchannels [18]. Then, the trapped cells inside the microwells start to aggregate and subsequently, a multicellular spheroid structure is formed by self-assembling of cells [19]. The continuous, long-term flow of culture media through the channels provides the cells with their required nutrition and removes cells' waste from the system [20].
Depending on the availability of oxygen and nutrients, there are two critical regions in a tumor spheroid, namely necrotic and quiescent regions [21]. The huge lack of oxygen and nutrients in the necrotic region leads to cell death [22]. The quiescent region is where the lack of oxygen and nutrients causes the cells to secrete growth factors that may lead to angiogenesis [20,23,24]. Barisam et al. numerically studied necrotic core and quiescent zone in the spheroids cultured in U-shaped barrier microfluidic chips and investigated the effect of several hydrodynamic parameters on multicellular aggregates in such platforms [23]. In another work, toroidal and spherical 3D cellular aggregates were simulated numerically in both microwells and U-shaped barrier chips [20]. It was shown that although the U-barrier design provides spheroids with a better concentration of oxygen, it also exposes cells to higher values of fluid shear stress, which is a negative point in culturing spheroid on microchips. Grimes et al. investigated the effect of spheroid diameter on necrotic and hypoxic development both numerically and experimentally and concluded that the increase in diameter would decrease the percentage of the proliferating zone [25]. Recently, Im et al. studied the hypoxic zone of stem cell spheroids and their gene expression experimentally [26].
Consequently, based on the literature, changing the microchip's geometry, increasing the media's flow rate, and decreasing spheroids' diameters are among the conventionally proposed solutions to battle the growth of necrotic and quiescent zones [13]. Changing the geometry is not always a good solution because of the difficulties in the fabrication of other designs and challenges corresponding to their operation and function. These challenges can induce high shear stresses on cells and make initial cell seeding even more arduous [13,20,23]. Moreover, culturing small spheroids is not always favorable since they may not resemble the in-vivo geometry of cancerous tumors.
Increasing the flow rate is another possible solution to prevent necrotic cores of in vitro spheroid culture platforms. While this approach improves the nutrition distribution inside the spheroid, it suffers from serious downsides: 1) Higher flow rates introduce more culture media and drugs in long-term culturing and drug screening. This increment in reagent consumption might be too expensive in some cases, and it is in contrast with the goal of using microfluidic technology, which is less consumption of materials and reagents [27]. 2) A higher value of fluid shear stress will be imposed on the spheroid, which can damage the plasma membranes of the cells forming the spheroid and may lead to cellular damage and death [28]. 3) By increasing the flow rate, the lift force on the spheroid rises, and the spheroid may move to the outlet and get wasted [29].
On the other hand, acoustofluidics is an emerging field that can address the potential problems associated with conventional microfluidic platforms. A comprehensive review of the application of bulk and surface acoustic waves in microfluidics has been conducted in the acoustofluidics tutorial series [30] and other reviews with a principal focus on surface acoustic waves [31,32]. Many studies have targeted microfluidic acoustic waves for biological purposes, especially for cell/particle manipulations. For example, Li et al. have utilised surface acoustic waves to separate circulating tumor cells from blood samples [33]. Using concentrated surface acoustic waves, Shilton et al. managed to aggregate and separate particles and mix fluids within a microfluidic droplet [34]. Incorporating standing surface acoustic waves into their microfluidic chip, Li et al. co-cultured cancer cells and endothelial cells in their novel study [35]. In a relevant work, Greco et al. managed to highly increase the cell-proliferation rate in a microfluidic chip with surface acoustic waves compared to its static condition [36]. Regarding the modeling of acoustofluidics, Muller et al. adopted perturbation theory to model the acoustic streaming phenomenon inside a microfluidic chip [37]. Raghavan et al. developed a numerical model to simulate the acoustic streaming effect inside a droplet. Although their model predicted the experimental velocity trends inside the droplet, it was not accurate as the velocity across the height of the droplet obtained from simulations was an order of magnitude lower than the one obtained from experiments. [38].
To address the aforementioned challenges in spheroid-on-chip microsystems, in this study, we took a step forward towards integrating acoustic microfluidics with conventional spheroid-on-chip platforms as a novel technique to decrease consumption of the reagents and decrease the associated shear stress on the cultured cells. The proposed model can also shrink the necrotic and quiescent zones in spheroids while preventing spheroid eluding the well and avoiding intense fluid shear stresses in the system. The present concept is feasible for fabrication and operation and can also be the starting point of efficient, low-cost, and long-term acoustic cell culture platforms.
The acoustic field inside a microchannel directly affects the flow pattern, thereby indirectly influencing the mass transport processes that are vital to living spheroids on the chip [39]. So, here, the aim is to utilise acoustic waves to improve culturing conditions of spheroids on microchips and reduce the consumption of culture media in these devices. To do so, the numerical simulation of the system is discussed in detail, and the effects of acoustic and hydrodynamics parameters on cultured cells are investigated. Moreover, a comparison between on-chip spheroid culturing with and without acoustofluidic integration is carried out to better illustrate the benefits of integrating acoustic field into spheroidon-chip systems.

Geometry and model description
The geometry of the system is illustrated in Figure 1. The spheroid is cultured in a microwell located beneath the perfusion microchannel. The dimensions of the microfluidic platform and the 3D cell aggregate are listed in Table 1. Culture media passes through the perfusion channel and leaves the system via the outlet. At the bottom surface of the microwell, an ultrasonic piezoelectric transducer is located, which vibrates at its resonance frequency and propagates acoustic waves to the flow when attached to electrical voltage. The boundary vibration disturbs the flow and brings about a more homogenous culture media through the acoustic streaming phenomenon.

Governing equations
Based on the physics of the problem, a 2D computational domain has been considered for the simulation. All the parameters and properties of the governing equations are listed in Table 2. It is assumed that the culture medium's properties are similar to those of water at 37 •C [23]. Moreover, in this study, dissolved oxygen and glucose in the culture media are the main parameters of interest that will be discussed in details.

Microfluidic flow
We consider an incompressible and steady laminar flow of culture medium inside the channel and around the spheroid. So, the continuity and momentum equations can be written as follows [23]: where ⃗ is velocity vector, p is the pressure of the fluid and and are density and viscosity, respectively. ⃗ represents the net value of the volume forces, which is due to the forces that the acoustic field exerts on the flow. The flow velocity inside the spheroid is set to be zero, and its outer surface works as an impermeable wall for the medium flow.

Transport of dilute species
In order to model the spheroid, a continuum modeling approach was implemented. The other model, namely the discrete approach, is suitable for the simulation of cellular aggregates' growth since it has control over cell interactions. The continuum approach is suitable for macroscopic analysis with a focus on epigenetic variations such as environment-related variables. As such, the general form of the diluted species mass transfer equation applies [23]: where represents the concertation of oxygen/glucose, t is time and . is the respective diffusion coefficient of either oxygen or glucose. models oxygen/glucose consumption by living cells in the computational domain and is defined by Michaelis-Menten reaction equation as follows [40]: where is the maximum reaction rate and is Michaelis-Menten constant. To model the culture media, the steady-state form mass transfer without the effect of cells can be used. As such, the transport phenomenon of culture media flow is purely governed by diffusion and convection, and Eq. 3 can be simplified as follows: Inside the spheroid, there are no flows. Only oxygen/glucose consumption by cells and diffusion of species within the spheroid's tissue exist. Consequently, the convective term of Eq. 3 should be omitted, and finally, Eq. 6 controls the mass transport inside the spheroid's inner domain:

Acoustic
Considering a compressible fluid, the governing equations are as follows: Eq. (7A) is the well-known continuity equation for a compressible fluid, and Eq. (7B) is the general Navier-Stokes equation by considering the transient (the term added in the left-hand side of the equation) and compressibility (added to the right-hand side of the equation, in which stands for viscosity ratio) effects. As this nonlinear equation does not have an analytical solution, perturbation theory is a conventional method to obtain a reasonable approximate solution. In this method, a quiescent fluid with constant density , constant pressure and zero velocity is considered before the wave incidence. As the wave influences the fluid, it disturbs the flow. Considering very small changes in the flow properties, the following first-order approximation can be made: Considering the process to be isentropic, pressure can be written as a function of density: In which, is the isentropic speed of sound defined as = . Substituting equations (8) and (9) into equation (7) and solving for the first-order terms results in: In obtaining Eq. 10 the product of 1 st order terms have been neglected as their value is small compared to the other ones.
Navier-Stokes equations are nonlinear in general, and approximating them by a firstorder term may introduce some errors into the final results. As a result, it is common practice to continue Eq. 8 to second-order terms: (For the sake of simplicity, second-order energy equation has been neglected [37].) Using the definition of the sound velocity, it is possible to write the pressure as a function of density as follows: Inserting equations 11 and 12 into the Navier-Stokes equations (Eq. 7) and solving for the second-order terms yields: In the derivation of Eq. 13, the product of higher-order terms has also been neglected. Generally, second-order terms are negligible compared to the first-order ones, except for the harmonic-time-dependent cases in which the time average of the terms in Eq. 10 disappears; however, some of the terms in Eq. 13 have a non-zero time average [41]. Defining < > as the time average function for the variable ( ): and considering the following harmonic time-dependence for the acoustic variables: Upon inserting these terms into Eq. (13), this set of equations converts to the following ones: in which < > is the streaming velocity.

Microfluidic flow
Fully-developed flow and zero pressure were imposed at the inlet and the outlet, respectively. No-slip boundary condition was considered on all walls and also at the interface of the spheroid and the culture medium.

Transport of dilute species
At the inlet, constant concentrations of oxygen ( ) and glucose ( ) inflow was applied. At the outlet, the gradient of concentration for both glucose and oxygen was set to be zero. All walls are impermeable, and as a result, a no-flux boundary condition was used for them. At the medium culture's interface with the cell aggregates, concentration jump boundary condition is applied by utilising a solubility coefficient (S) as in equations 17 and 18 [23]: Conservation of fluxes was also fixed at the interface of the culture media and the spheroid [23]: . . (20) Also, concentration conditions that describe the necrotic and the quiescent zones are as follows [42,43]: For the necrotic zone: For the quiescent zone:

Acoustic
The first-order velocity boundary condition is as follows: In Eq. (20), is the amplitude of the wall normal displacement and is the stress tensor which is based on the first-order velocity. For the first-order temperature, all the boundaries are considered to be isothermal, and for the second-order velocity equation, the no-slip boundary condition was applied on the walls, and a fully developed boundary condition was considered for the inflow and outflow, as described in Section 2.3.1 All constants and parameters of the governing equation and their boundary conditions are presented in Table 2.

Numerical method
The abovementioned equations were discretised and solved utilising an FEM code. All the domains in all physics were meshed with triangular elements. A boundary layer mesh was generated for all the wall boundaries to completely capture the acoustic streaming effects taking place within the viscous penetration depth. For the convergence criteria, a residual value of 10 was set for Navier-Stokes and continuity equations, 10 was set for the transport equation, and 10 was used for thermoviscous acoustic equations.
The numerical procedure to solve the present multiphysics problem is as follows: Frst, thermoviscous equations were solved in all the domains to find the first-order pressure and velocity. Having obtained the mentioned fields, they are inserted in Eq. (15) as source terms to obtain the second-order terms, which are the fluid dynamic velocity and pressure of the domain. Finally, the obtained velocity and pressures are used to solve the diluted species transport equations and find the concentrations of oxygen and glucose within the domain.

Mesh-independent study
Specific care was given for the mesh generation process pertinent to the thermoviscous and laminar flow solutions [37]. The viscous penetration depth, which is defined by the relation = √(2 / ) = 0.52 , was divided into 6 layers near the wall to fully capture the acoustic streaming phenomenon. In the well and in the free stream domain, maximum element size was set to 10 and 20 times the viscous penetration depth, respectively.
According to Muller et al. [37], streaming velocity was the last quantity to converge. As a result, in this study, a grid independence study is conducted on the shear stress, which is derivative of velocity and more sensitive to velocity variations and a better indicator of mesh independence. Accordingly, a grid study on shear stress aims not only thermoviscous acoustic solution but also laminar flow. To do so, the maximum value of fluid shear stress (FSS) on the spheroid was considered against mesh refinement. On the upper half of the spheroid, the flow of the fluid in the perfusion channel is mainly responsible for the FSS, and for the lower half, the acoustic field is. Another point of this consideration is to find out that if the magnitude of FSS after acoustic application does not surpass the allowable range mentioned in the literature [44]. In Figure 2, details of this investigation are presented. Figure 2. Mesh study graph for laminar flow and mass transport equations. As the grid size reduced in a step-by-step manner, both maximum shear stress and average oxygen concentration tend to be constant and independent of the number of the grids.
Mesh study's error remains below 3 percent in this section. To prove that the mass transport solution is not dependent on generated meshes, the average magnitude of oxygen's concentration in the whole spheroid was measured while the number of mesh elements increased. Simulation outcomes proved the point that the solution remains constant, independent of the number of meshes. Figure 2 shows the convergence of grid study. With 40368 elements, the error will be less than 1 percent compared to the finer sizes of meshes.

Validation of the study
We have previously validated the accuracy of our proposed numerical model for laminar flow and mass transport simulations [20,23]. For the validation of the acoustic model, the model of Muller et al. [37] was regenerated and simulated, whose geometry is shown in Figure 3. The width ( ) and height (ℎ) of the channel are 380 µm and 160 µm, respectively. For the first-order temperature boundary condition, all boundaries are considered to be isothermal. Therefore, the first-order velocity conforms to the following equations: Side Walls: . (± /2, ) = ± (24) Bottom and top walls: ( , ±ℎ/2) = 0 The comparison of the results for the first-order y-component velocity of this study and the one conducted Muller et al. [37] is shown in Figure 4.  Figure 3 between this study and the one obtained by Muller et al. [37].
This line graph is plotted on the dashed line of Figure 3, located /4 to the right of the domain center. Since the sharp gradients take place in the boundary's vicinity, only the near-boundary region is shown in Figure 4. The two graphs almost coincide, which validates the numerical method and boundary conditions used in this study.

Conventional spheroid-on-chip platform (No Acoustic)
In this section, the model was simulated without acoustic integration to obtain a broad view of the physics of the problem. Figure 5 inside the microwell and around the spheroid is negligible. In Figures 5(b) and 5(c), oxygen and glucose distributions and their fluxes are observable, respectively. Diffusion is mainly responsible for transporting nutritious species to the cell aggregate, while the convection term in the well and adjacent to the spheroid is not strong enough.

Acoustic spheroid-on-chip platform
After introducing acoustic ( = 0.5 and = 1 ) to the system, for the same flow rate of 1 µl/min, flow's streamlines start to form some vortexes near the microwell, which enhances convention inside the microwell (Figure 6(a)). Oxygen and glucose concentration distributions and their fluxes also are presented in Figure 6 In the following subsections, the effect of two principal factors, namely boundary displacement amplitude and the inlet flow rate, on the flow pattern inside the chip and the concentration level of glucose and oxygen in the spheroid vicinity will be studied.

Boundary displacement amplitude
The effect of boundary displacement amplitude ( ), which is controlled by the properties of the piezoelectric substrate and the actuation power, on the spheroid culturing is studied in this section.
represents the displacement amplitude of the actuated boundary and in practical setups is usually measured to be in the range of 0.1-0.5 nm [45]. In this section, frequency and flow rate are both kept constant at 1 MHz and 1 µl/min, respectively. Figure 7(a) shows the effect of on glucose and oxygen concentrations within the spheroid. As the amplitude of boundary vibration increases, the acoustic streaming effects escalate, which influences the fluid flow pattern more intensely and consequently, the convection in the microwell around the spheroid is improved. As a result of this disturbance in the flow pattern, the nutrient gradient in the vicinity of the spheroid will be decreased, which brings about a more homogenous culture medium in the spheroid proximity. Therefore, the mass transfer in the vicinity of the spheroid will improve, and oxygen and glucose can more easily and effectively diffuse to the core of the cell aggregate. Interestingly, for oxygen, increasing from 0.2 nm to 0.5 nm leads to a 228% rise in the average concentration of oxygen in the spheroid. For glucose, although the amplification is not so huge, still employing acoustics improves the distribution of glucose in the spheroid. Improving oxygen/glucose diffusion into the core of the spheroid ends in the shrinkage of both necrotic and quiescent zones and consequently growth of the proliferation zone.
(a) (b) (c) Figure 7. (a) Effect of on glucose and oxygen concentrations inside the spheroid. The higher the , the more oxygen and glucose will be found inside the cell aggregate, and the better cells in the inner side and the core will be nourished. (b) While increasing is beneficial due to enhancement of oxygen and glucose concentrations, it is not also causing any side effects such as high magnitudes of fluid shear stresses or lift forces. Here, represents the maximum value of shear stress which is exerted to the peripheral boundary of the spheroid. (c) Graphical illustration of oxygen's proliferation zone after applying the acoustic field. Before applying acoustic to the microchip, the proliferation zone had only a 32.7% share of the spheroid. With the acoustic field, this share is increased to 66.8%.
While introducing acoustics to the chip can vividly enhance the nutrient-enrichment of the spheroid, it is important to make sure that the integration of the acoustics is not damaging cells by imposing a high magnitude of shear stress on them or inhibiting the spheroid entrapment in the microwell due to exerting lift force on it. From the hydrodynamic point of view [46], the maximum allowable fluid shear stress on the cells is 0.5 / . Also, the extra produced lift force due to the acoustic field should not exceed the downward force which is the net of weight and the spheroid's buoyance force. Figure 7 proves that from the hydrodynamic point of view is within the safe range, and acoustic neither hamper the spheroid culturing process nor adversely affect cell viability due to shear stress. Figure 7(c) illustrates how increasing eliminates unwanted quiescent/necrotic zones and causes the proliferating zone for the oxygen to grow inside the spheroid.

Flow rate
Flow rate is the main controllable parameter that can directly govern the hydrodynamic field inside microchannels. Higher flow rates correspond with more mass of the species in the computational domain and consequently higher average concertation of oxygen/glucose inside the spheroid. Figure 8(a) shows how flow rate affects concentration and compares the average glucose/oxygen concentrations in acoustic and non-acoustic systems as a function of flow rate. As expected, there is an upward trend between the flow rate and the average concentration. For oxygen, for example, a 300% increase in the flow rate from 1 µl/min to 4 µl/min leads to the rise of average concentration from 0.013 mM to 0.028 mM. But, keeping the flow rate at 1 µl/min and just introducing the acoustic field to the domain, oxygen's average concertation will be pumped up to 0.046 mM. Study the effect of flow rate and acoustic integration to the system on maximum fluid shear stress and the lift force. Acoustic wave is amplifying these two parameters but they are still in the safe range for this application. (c) A comparison in necrotic/quiescent shrinkage with flow rate between acoustic and no-acoustic modes. Increasing the flow rate helps the growth of the proliferation zone. Acoustic field solely can do better without the need for any increases in the flow rate. Approximately, acoustic improves the proliferating zone by 100% percent at each flow rate.
Since acoustic wave inside the microwell is enhancing the hydrodynamic flow around the sphere, it might be interesting to also compare lift force and maximum shear stress in both acoustic ( =0.5 nm)/non-acoustic ( =0 nm) modes and also investigate how flow rate is affecting it. As shown in figure 8(b), it is true that the addition of acoustic to the system is exerting more shear stress and lift force on the cell aggregate, but still, it is way much below the thresholds. Figure 8(c) demonstrates that how necrotic and quiescent zones shrinkage as the flow rate is increased. When the acoustic field is introduced to the system, proliferating zone's share from the spheroid raised interestingly, and also, it can be seen that with acoustic ( =0.5 nm), the downward trend in the shrinkage of the necrotic zone is also sharper against the increase in the flow rate. In Figure 7(c) it is good to compare the first (with acoustic, =0.5 nm, Q=1 µl/min) with the last (no acoustic, q=12 µl/min) graphical illustrations. Based on simulation results, a low flow rate of 1 µl/min with the acoustic boundary displacement of 0.5 nm, results in 66.8% of proliferation zone, 19.9% of the quiescent zone, and 13.3% of necrotic zone. To achieve a similar percentage of zones without acoustic waves, the flow rate should be increased 12 times. This proves the point that acoustic improves cell viability astonishingly on the chip without the need for consuming more reagents. About 100% increase was observed in the area of proliferating zone after applying the acoustic field at each flow rate.

Conclusions
In this research, the concept of acoustic spheroid-on-chip microfluidic platforms was introduced and studied. The goal was to engineering a microsystem to evade the pitfalls of conventional spheroid-on-chip platforms. By implementing the proposed approach, a better cell viability rate is achievable without the need for utilising complex geometries, increasing the flow rate, or even culturing spheroids with smaller radii.
First, the proposed geometry and the concept of the model were described. Next, the numerical method and governing equations were presented in detail. And finally, simulation results were reported and discussed. It was observed that by implementing acoustic waves, the pattern of the flow, the distribution of concertation and mass fluxes throughout the microfluidic chip were changed. As a result of the acoustic streaming phenomenon, fluidic convection was enhanced in the microwell and around the spheroid. Consequently, more oxygen/glucose was transported to the vicinity of the spheroid, which itself ended up in higher concentrations of nutritious species inside the cell aggregate.
The effect of acoustic power on spheroid culturing was studied via the wall displacement amplitude ( ). It was shown that increasing this parameter results in the extension of the proliferating zone. Compared to the non-acoustic mode, with a constant, small flow rate, the share of proliferating zone in the spheroid can be doubled with the aid of acoustics. The effect of flow rate was also investigated. It was concluded that although higher flow rates cause shrinkage in both necrotic and quiescent zones, the introduction of acoustic field to the system enhances this effect without the need for consuming more culture media or treating reagents.
The possible drawback of acoustic integration is the production of high magnitudes of fluid shear stresses and lift forces. We numerically calculated both shear stress and lift forces and proved that they were in the safe range with no adverse effect on the cell aggregate. In summary, the proposed acoustic microfluidic cell