Array-type MF devices have been extensively used for spheroid formation and drug screening (Liu et al., 2015; Sabhachandani et al., 2016; Moshksayan et al., 2018). They can accommodate hundreds of cell-laden droplets in the microwells (Wang et al., 2016; Huang et al., 2017). As a case study, our computational study focused on a particular design shown schematically in Figure 2. The MF device consists of parallel rows of cylindrical microwells positioned on a side of a microchannel “supply line.” The parallel rows of channels are connected to a single inlet and a single outlet via serial divisions. Cell-laden droplets are compartmentalized in the microwells, and subsequently, transform into MCTSs embedded in a hydrogel scaffold. The geometry of the microwells determines droplet size, which in turn, controls MCTS dimensions (Kaminski and Garstecki, 2017). This MF platform has the following advantages: i) the capability to form uniformly sized MCTSs, ii) the ability to monitor individual MCTSs in real-time, and iii) long-term MCTS culture, enabled by the use of a suitable hydrogel scaffold. Because droplet generation is controlled by the relationship between the channel geometry and the microwell aspect ratio (Schneider et al., 2013), the design of such MF devices to study, for example, the biomechanical forces exerted on the MCTSs or drug delivery to the MCTSs is a challenging task. While the effect of MF device geometry on droplet formation has been studied (Cohen et al., 2010), the effect of these constraints on fluid flow and transport of biological molecules such as drugs or nutrients has yet to be examined (Lee et al., 2012; Schneider et al., 2013; Wang et al., 2016).

Here, our goal was to simulate the effects of geometrical and non-geometrical constraints on the transport of drugs and nutrients to the arrays of MCTSs in the spheroid-on-a-chip MF platform illustrated in Figure 2A. In all simulations, MCTS diameter was 100 µm. To demonstrate the versatility of our computational strategy, we studied both the single-supply-channel (SSC) design (shown in Figure 2A) and a modified double-supply-channel (DSC) design that has not been reported (Figure 2B). Five different geometric parameters of the MF device were examined, each with a specific range of values (Table 1 ).

Molecular transport was modeled as time-dependent and three-dimensional, and included advection and mass transfer in the porous domain. Molecular transport was thus governed by the advection-diffusion equation: ∂ c ∂ t = S + D ∇ 2 c − u ⋅ ∇ c (1) where   c is the species concentration (mol/m³),   S is the supply influx (mol/m³⋅s),   u is the fluid velocity (m/s), and D is the diffusion coefficient of the molecular species (m²/s). Fluid flow was considered time-dependent, three-dimensional, laminar, and incompressible, and was thus governed by the continuity and Navier-Stokes equations: ∂ ρ ∂ t + ∇ . ( ρ u ) = 0 (2) ρ [ ∂ u ∂ t + u ⋅ ∇ u ] = − ∇ p + μ ∇ 2 u + ρ f (3) where   ρ is the fluid density (kg/m³); u is the fluid velocity (m/s),   p is the pressure (Pa), μ is the dynamic viscosity of the fluid (Pa⋅s), and f is the body force per unit volume (N/m³). Eqs 2 and 3 represent conservation of mass (continuity equation) and momentum (Navier-Stokes equations), respectively. Here we assumed the body force term was zero as we ignored gravity and assumed there were no other external field-induced forces present. To study drug transport and penetration into spheroids, the spheroids in our simulations were treated as porous media (defined as a “matrix domain” in COMSOL) surrounded by the primary fluid domain where MF flow was occurring. Porous media flow within the spheroid was governed by: ∂ ∂ t ( ρ ε p ) + ∇ . ( ρ v ) = Q m (4) where ε p is the porosity (or void fraction) of the matrix domain,  v is the velocity inside the matrix domain, and Q m is a mass source within the matrix domain (kg/m³·s). The velocity field within the porous domain can be calculated using Darcy’s law: v = − κ μ ∇ p (5) where κ and μ are the permeability of the porous material (m²) and dynamic viscosity of the fluid material (Pa⋅s), respectively. Solving the velocity field within the matrix domain determined the advection mode of drug transport. The diffusion mode of transport, in contrast, was calculated using Fick’s law: J = − D ∇ c (6) where   J is the diffusive flux (mol/m²⋅s), which points in the negative direction of the concentration gradient. Finally, drug transport was solved through mass conservation in which both advection and diffusion modes of transport were included using the following: ∂ c ∂ t + ∇ ⋅ J + v ⋅ ∇ c = R (7) where R is the drug reaction rate (mol/m³⋅s). This reaction rate was set to zero in all our simulations.

To tackle this large combinatorial problem of >15,000 modelling scenarios, we employed parallel computing using the CAD Compute Cluster from CMC Microsystems and divided the various scenarios into smaller subsets to be run simultaneously on different parallel processors (Afzal et al., 2017). Using this large computational resource (eight nodes, Dual 14-core, 2.4–3.3 GHz CPU), the total computational time was reduced ∼10-fold compared to running the same CFD simulation on a single desktop workstation (7-core, 2.1–3.2 GHz CPU).

For all simulations, pressure, velocity, and concentrations were analyzed and evaluated in 2D and 3D configurations. The 2D and 3D plots were generated using COMSOL Multiphysics plotting functions. For 1D plotting and investigation, the data were extracted and exported from the CFD model using volume averaged evaluation functions. Five key metrics of flow and transport in the spheroid domain were extracted: concentration, velocity, shear rate, advective flux, and diffusive flux. These metrics were normalized between 0 and 1 before being compiled into heatmaps. Clustered heatmaps were generated based on the exported data and post-processing by the Seaborn library and clustering functions in Python (Supplementary Material). For convenience of describing specific device designs, we established a naming convention based on the concatenation of the six design parameters selected for the given design. An executable standalone application was also designed and created using the application builder module in COMSOL Compiler (Supplementary Material).

We conducted numerical simulations with our model to examine the effects of different parameters on biomechanical forces, fluid transport, and drug uptake in and around tumour spheroids cultured on-chip. First, a single row of the device, which includes 50 microwells, was modeled to determine the steady state pressure gradient and velocity field in the microchannel (Supplementary Figure S4). Using a prescribed inlet flow rate of 0.01 ml/h, we found the maximum pressure in the microchannel to be 17.5 Pa and confirmed that the pressure drop between any two adjacent microwells was consistent along the entire length of the microchannel (Supplementary Figure S4A). Fluid velocity inside the microwells was ∼10× lower compared to the fluid velocity in the main microchannel (Supplementary Figure S4B ). Notably, the velocity magnitudes inside the chambers were found to be ∼0.02 mm/s (Supplementary Figure S4C), consistent with typical physiological conditions in tumour tissue (Chary and Jain, 1989).

To explore a larger design space, we examined the DSC design consisting of two parallel supply microchannels on either side of the spheroid-containing microwells (Figure 3A). Note the DSC design consisted of the same set of geometrical parameters and range of values as the SSC design. Flow rates for both channels were first prescribed at 0.01 ml/h. Under this condition, velocity magnitudes inside the spheroid chambers were ∼0.02 mm/s, similar to the SSC design (Figure 3B). We examined fluid flow circulation inside the microwell (Figures 3C,D ) and found that the streamlines formed two flow vortices that surrounded the spheroid. These vortices in the DSC device were not observed in the SSC design ( Supplementary Figure S2 ). Furthermore, adding the second supply microchannel not only improved circulation of the fluid inside the microwell, but also improved the uniform distribution of surface shear forces on the spheroid (Figure 3E) compared to the SSC design (Supplementary Figure S2). We next simulated fluid flow inside the spheroids (treated as a porous medium) to better understand the effects of adding the second supply microchannel to flow through the internal spaces of the spheroid (Figure 3F). We found that fluid motion had two high velocity regions corresponding to the locations closest to the supply microchannels on the two sides, while the rest of the spheroid observed velocity magnitudes that were an order of magnitude lower. Taken together, these initial simulations revealed several major differences between SSC and DSC designs that motivated further exploration into ways to optimize the design to achieve maximum drug transport efficiency for the purposes of developing a rapid high-throughput spheroid drug screening tool. Before conducting such a massive design space exploration, however, we experimentally validated our computational model by using the actual spheroid-on-a-chip platform and the select experimental flow conditions.

We conducted a series of experiments to validate our computational model, using an SSC array-type MF device to form MCF-7 breast cancer MCTSs and examine the effect of different flow rates on drug uptake. The experimental MF device contained 200 MCTSs in the microwells that were organized in four parallel rows (Figure 4A). The microwells were connected to a common channel used to supply cell culture medium and drug treatment [for more details on fabrication, MCTS formation, and cell staining protocols see Supplementary Material and (Chen et al., 2021; Prince et al., 2021)]. The cell suspension occupied 100% of the microwells in the form of cell-laden droplets; MCTSs formed within 72 h of cell culture. The viability of breast cancer cells in the MCTSs was evaluated using NucBlue (Hoechst 33342, nuclei) and NucGreen (dead cells), showing ∼89% viability throughout the entire spheroid array (Figure 4B ). We also confirmed the expression of E-cadherin and F-actin via immunostaining to examine proper localization of cell-cell junctions and the presence of actin cytoskeleton, respectively (Figure 4C). Next, MCTSs were treated with 10-µM Doxorubicin (Dox) under two different flow rates, Q, of 0.01 and 0.02 ml/h. The results showed that after 2 h of drug treatment by perfusion, drug uptake at Q = 0.02 ml/h was significantly higher than the uptake at Q = 0.01 ml/h (Figure 4D). To validate our computational model, we simulated the experimental geometry and flow conditions in silico using the same corresponding geometry and conditions ( Figure 4E). As predicted, the fluxes of drug molecules into the MCTSs at Q = 0.02 ml/h were significantly higher than at Q = 0.01 ml/h (Figure 4F ). The increase in flow rate resulted in approximately five-fold increase in the total flux of the drug. Moreover, we compared and analyzed the penetration of drug molecules into the MCTSs using the CFD model (Supplementary Material). The results revealed that the increase in the flow rate of drug treatment significantly enhanced the penetration of the MCTSs with a drug (Figure 4G ). Thus, the experiments validated our model, thereby enabling us to study changes in drug transport under different conditions and uncover optimal MF device design parameters for a particular application.

To better understand and compare the effects of device design on drug delivery and biomechanical forces on the spheroids, we used a systematic computational approach that allowed the exploration of >15,000 device designs and flow scenarios in a single simulation dataset. Our computational model covered a large parameter space representing a full-factorial experimental design, the results of which can be used during the design stage of the MF device to study the impact of design changes on drug transport. In total our complete dataset includes 10 billion datapoints (∼15,000 scenarios × ( ∼ 7 × 10 5 ) finite elements per fluid flow and porous medium domain) resulting in ∼95 GBs of data. In the following sections, we describe results from the massive dataset, where we simulated the drug transport of Dox, a widely used drug for tumor therapy (Mehta et al., 2012), as our model drug.

Our parametric “design screen” of >15,000 simulations enabled us to systematically investigate the effects of individual parameters on drug uptake a posteriori. One physical property of spheroids that plays an important role in drug delivery is the spheroid porosity (Soltani and Chen, 2012; Soltani et al., 2016). Here, we focused on the effects of spheroid porosity on the velocity field and drug delivery inside the spheroid when the inlet flow rates Q₁ and Q₂ were both 0.01 ml/h, and considered advective and diffusive modes of transport separately. We studied the effects of porosity by extracting the simulated results in DSC designs and four different porosities, ε = 0.2 , 0.5, 0.7, and 0.9, while keeping the geometry, fluid flow, and drug properties the same. We found that changes in porosity increased the total flux of the drug by ∼17× when porosity increased from ε = 0.2 , to ε = 0.9 (Figure 5C). However, velocity magnitudes (as viewed on the sagittal yz plane) were similar across all porosities. Thus, spheroid porosity clearly influences the drug/nutrient delivery, but has minimal effect on the velocity magnitudes within the spheroid.

Diffusive and advective fluxes are concurrent processes in biological tissues where fluid flow exists, such as the interstitial flow within the stroma. Studying both diffusion and advection transport modes in the spheroid-on-a-chip provides insight into the dynamics of drug uptake while also opening the door to applying our computational models for novel drug design based on transport efficiency (Minchinton and Tannock, 2006; Dewhirst and Secomb, 2017). The delivery of adequate concentrations of anticancer drugs to cancer site strongly depends on the structure and drug transport mechanism in tissue. While the effects of structure on tumor response has been studied (Steuperaert et al., 2017; Li et al., 2021), the effect of transport-related factors that play an important role in the delivery of anticancer drugs are much less studied (Dewhirst and Secomb, 2017). Our model considered both diffusion and advection modes of transport separately, facilitating the evaluation of drug transport mechanisms in our different device designs.

Using our model, we assessed the drug transport modes for Dox under various conditions (see summary of conditions in Supplementary Table S1). First, we arbitrarily selected one device design (W_R = H_R = W_L = H_L = 25 μm, Q₁ = Q₂ = 0.01 ml/h, ε = 0.5 , C_dox = 10 µM) and varied only the chamber radius from 80, 120, and 160 µm (the complete dataset is summarized in the sections below). When we compared the advection and diffusion transport modes by visualizing 3D volume contours of flux, we found significant differences between the rates of change in diffusion and advection fluxes (Figure 6 ). Closer inspection of the graphs revealed that the transport of drug molecules was diffusion-dominant for smaller microwells (80 μm) but advection-dominant for larger microwells (160 μm). Note that we refer to a transport mode as “dominant” if the volume-averaged flux of the dominant mode is (or exceeds) 10x greater than the volume-averaged flux of the other mode. If the flux ratio between modes is less than 10, we consider it a mixed-mode scenario. Our simulations showed that even a single change in one geometric feature (e.g., chamber radius) can substantially affect the drug uptake mechanisms. These results can be further used as a convenient tool to assist pharmacodynamics studies and help identify limiting factors for drug delivery to spheroids cultured on-chip.

As stated above, adding a second supply microchannel to the device creates vortices inside the microwell and significantly affects the fluid flow around the spheroid. To examine the impact of device geometry and flow patterns on drug transport, we studied drug penetration and uptake by measuring the drug concentration at the center of the spheroid over time using a 3D point probe in COMSOL.

Drug concentration time profiles at the center of the spheroid were plotted for a myriad of different conditions (Figure 7). We applied two different flow conditions of interest, corresponding to when supply channel flow rates were i) equal (Q₁ = Q₂) (Figure 7A), and ii) unequal (Q₁ = 2 × Q₂) (Figure 7B). For each of these flow conditions, the dataset was further divided by varying spheroid porosity ε (rows) and chamber radius R (columns). Finally, for each specific ε and R, we plotted four different groups of curves that allowed us to compare effects of geometry: i) all SSC designs, ii) our “reference” DSC design where H_R = H_L = W_R = W_L = 60 μm (maximum connecting channel size), iii) varying H_R and H_L independently (W_R = W_L = 60 μm), iv) varying W_R and W_L independently (H_R = H_L = 60 μm). The “reference” DSC design was chosen arbitrarily from the set of simulated designs to allow convenient comparison between different designs.

For equal supply channel flow rates (Figure 7A ), drug delivery was clearly more efficient in all DSC designs compared to the SSC designs, as the drug concentration at the spheroid core was higher for all DSC designs compared to SSC designs in the same amount of time. The reference DSC design was the most efficient, as the spheroid core reached maximum concentration in the shortest amount of time due to the maximum cross-sectional area of the connecting microchannels. Additionally, changing the height of the connector (H_R and H_L) showed a larger effect on drug delivery efficiency compared to changes in the width of the connecter (W_R and W_L). This effect was more apparent in the simulation with a smaller microwell diameter and less porous spheroids. This can be attributed to the fact that by increasing the size of the microwell, the dominant drug transport mode changes from diffusion to advection (as mentioned above), while by increasing the connector height, more drugs will enter the spheroid at a higher location that is closer to the spheroid core.

To explore the effect of flow rate on drug delivery, the same device designs were simulated under the same conditions, but with one supply channel flow rate at twice the flow rate of the other (Q₁ = 2 × Q₂). The same trends of drug delivery efficiency were observed (Figure 7B ). Interestingly, we observed that unequal flow rates in the supply microchannels further improved drug delivery efficiency compared to equal flow rates, as the amount of drug reaching the spheroid core was higher in the same amount of time for unequal flow rates (Q₁ = 2 × Q₂) than for equal flow rates (Q₁ = Q₂). For example, when φ = 0.6 and R = 120 μm, drug concentration at the spheroid core reached a plateau at 1 × 10 − 3 mol / m 3 at t = 60 min when Q₁ = 2 × Q₂, but drug concentration did not reach this plateau for Q₁ = Q₂, even after t = 120 min. This difference was likely caused by a pressure gradient between two connectors that was generated by the unequal flow rates, ultimately leading to more drugs being actively driven inside the microwell (Supplementary Material).

Because of the massive dataset generated by our computational design screen, design space exploration became complicated, cumbersome, and time-consuming. To reduce the complexity of the design space and summarize the CFD simulation results in a more convenient fashion, we employed novel classification and clustering methods to sort the various conditions and flow scenarios into groups and reveal underlying patterns in the design parameters ( Figure 8). To evaluate a given device design, we examined the volume-averaged value of five key metrics of flow and transport in the spheroid domain: concentration, velocity, shear rate, advective flux, and diffusive flux. Hierarchically clustered heatmaps were generated for the cases of unequal supply flow rates (Q₁ = 2 × Q₂) (Figure 8A ) and equal supply flow rates (Q₁ = Q₂) (Figure 8B ). These heatmaps summarized the five-key metrics for all 15,000 simulations, enabling us to rapidly and succinctly compare different design effects on the transport of drugs in our principal spheroid-on-a-chip layout.

Hierarchical cluster analysis revealed that for both Q₁ = 2 × Q₂ and Q₁ = Q₂, shear rate, velocity, and advective flux were grouped together as flow and transport metrics, while both concentration and diffusive flux were differentiated from this primary cluster with only moderate similarities between them. As expected, this clustering predicted that modifications of the device design will affect fluid velocity, shear stress, and advective flux with the same tendency, and that these metrics cannot be varied or controlled separately. On the other hand, the heatmap appeared to show that specific designs can have high concentration and low diffusive flux or vice versa, which illustrates the potential to tune the design based on these two metrics. In terms of clustering of the designs, the heatmap showed that generally, in SSC devices, the drug uptake is mainly diffusion-dominant, and the exerted shear stress and velocity on the spheroids are lower compared to the DSC device design. However, the total amount of drug transport in DSC devices appeared to be higher compared to the SSC device, which is in line with the results in Figure 7.

Our hierarchically clustered heatmaps improve our understanding of the effects of different parameters on the metrics of flow and transport in the spheroid-on-a-chip layout and facilitate the selection of design parameters prior to any fabrication or experimentation. When designing MF devices, choosing the right range of parameters requires substantial effort and time. Hierarchical clustering helps narrow the design parameter space and empowers the researcher to quickly find the optimal design solution for a given application. To illustrate, we considered a hypothetical experiment that involved studying the effects of shear stress on drug delivery, and asked whether we could select one specific device design from the heatmap that could a) achieve rapid drug penetration (i.e., high drug concentration in the spheroid core in the shortest amount of time), b) permit both high and low applied shear stress on spheroids by changing operating conditions only, and c) keep all other key metrics relatively invariant to create a controlled experiment. Indeed, when we examined the heatmaps in Figures 8A,B closely, we discovered that the design (P9.R8.W_r6.H_r6.W_l6.H_l2) satisfied all these criteria (Figure 8C), since high shear stress can be achieved with Q₁ = 2Q₂ (Figure 8A) while low shear stress can be achieved with Q₁ = Q₂, all in the same geometry. Similarly, as another example, we considered a hypothetical experiment where we examined the effects of advective-dominant vs mixed advective-diffusion transport modes on drug penetration into spheroids. We were able to identify that the design (P9.R8.W_r.H_r6.W_l3.H_l2) could achieve advective-dominant transport when Q₁ = Q_2, where advection is ∼23× greater compared to diffusive flux. To switch to the mixed mode of transport where advective and diffusion fluxes contribute similarly (i.e., diffusive flux is only 2.3× greater than the advective flux), the same design can be used with Q₁ = 2 × Q₂ (Figure 8D). Hence, generating hierarchically clustered heatmaps enabled us to summarize and predict drug delivery distribution and efficiency and identify a range of parameters for a particular set of experiments or applications. Thus, our big data “design heatmap” can ultimately offer a rational approach to the design process that helps us select appropriate designs for proposed experiments in a more time-efficient manner.

Since the primary objective of our computational modelling framework (including our post-processing methods) is to provide the MF engineer with a useful design tool to reduce the design iteration time, guide experiments, and offer operating conditions for end-users of the spheroid-on-a-chip devices, we ventured to create a standalone user-friendly application that can be installed and executed on different devices and operating systems, including Windows, macOS, and Linux (Supplementary Material). This standalone app, which we call “SINAI” (Spheroid-chip Investigation by Numerical modelling-Application Interface), enables engineers, biologists, and other MF researchers to modify and run the simulations directly—even without the need for parallel computing resources—and with full control over geometrical, biomaterial, and fluid parameters. The app allows the sharing and dissemination of computational modelling knowledge widely and openly, including to those who do not have a COMSOL Multiphysics license. We envision that this full-factorial computational framework will have the potential to accelerate MF device design for future spheroid-on-a-chip systems, and more broadly for other organ-on-a-chip systems.

We note that every computational model is by its nature a simplification of the entire system; however, in comparison to experiments, our model covered a comprehensive range of flow rates of liquids, microwell and supplying channel dimensions, and MCTS porosity to provide a broad spectrum of possible combinations that may result in desired or undesired effects for MCTS formation and drug delivery. Typically, the process of parametrical sweeping and its range is limited to the computational techniques and capacity. In our in-silico model, we utilized parallel computing on a Cloud cluster to improve the performance of this sizeable computational task in a reasonable time. Sweeping and post-processing over the parameter space must be performed for each independent design, which will lead to thousands of design modifications that represent refinements of the principal design. It is conceivable to repeat the process for as many different conceptual designs as needed for the given design problem. Indeed, our approach cannot generate completely novel design concepts automatically: the microfluidics design engineer must continue to play the creative role of generating innovative design concepts that can then be refined extensively using our modelling workflow. In the case described here, we conceptualized SSC and DSC configurations first, and then utilized our modelling and post-processing approach to determine optimized design and operating parameters.

Our work focused on the interplay between both diffusive and advective modes of drug transport and the design of MF devices, taking into account the geometry of the MF device and porosity of MCTSs. The advective component of drug transport plays an essential role in drug delivery in vivo, as the molecular transport of drugs with a molecular mass M_r > 1,000 is advection-dominant, as opposed to the transport of low molecular mass drugs with M_r < 1,000, which are diffusive-dominant (Dewhirst and Secomb, 2017). To examine this effect, we evaluated the dominant drug transport mechanisms in MF devices with various designs. The CFD model facilitated the assessment of the dominant transport mechanisms of drugs in MF devices with different designs.

Our study shows that it is more favourable to use MF devices with a DSC design that enable more rapid delivery of the model drug (DOX) to the MCTSs, compared to the MF device with a single supplying channel. The leading cause of this effect is attributed to the fact that flow in a two-channel MF device creates a pressure gradient in the microwell that results in a more efficient delivery of drug molecules to the MCTSs than in a single-channel MF device. Given the complexity of finding the relationship between various coupled (interdependent) geometrical parameters and non-geometrical conditions, as well as the significance of parametrizing the model by integrating the data into experimental data, we developed hierarchically clustered heatmaps to provide guidance for interpreting data and designing and designing new experiments for spheroid-on-a-chip platforms.