Skip to main content

ORIGINAL RESEARCH article

Front. Bioeng. Biotechnol., 23 November 2021
Sec. Biosensors and Biomolecular Electronics
This article is part of the Research Topic Advances and Innovations in Microfluidics View all 6 articles

Computational Modelling and Big Data Analysis of Flow and Drug Transport in Microfluidic Systems: A Spheroid-on-a-Chip Study

Sina KheiriSina Kheiri1Eugenia Kumacheva,Eugenia Kumacheva2,3Edmond W.K. Young,
Edmond W.K. Young1,3*
  • 1Department of Mechanical and Industrial Engineering, University of Toronto, Toronto, ON, Canada
  • 2Department of Chemistry, University of Toronto, Toronto, ON, Canada
  • 3Institute of Biomedical Engineering, University of Toronto, Toronto, ON, Canada

Microfluidic tumour spheroid-on-a-chip platforms enable control of spheroid size and their microenvironment and offer the capability of high-throughput drug screening, but drug supply to spheroids is a complex process that depends on a combination of mechanical, biochemical, and biophysical factors. To account for these coupled effects, many microfluidic device designs and operating conditions must be considered and optimized in a time- and labour-intensive trial-and-error process. Computational modelling facilitates a systematic exploration of a large design parameter space via in silico simulations, but the majority of in silico models apply only a small set of conditions or parametric levels. Novel approaches to computational modelling are needed to explore large parameter spaces and accelerate the optimization of spheroid-on-a-chip and other organ-on-a-chip designs. Here, we report an efficient computational approach for simulating fluid flow and transport of drugs in a high-throughput arrayed cancer spheroid-on-a-chip platform. Our strategy combines four key factors: i) governing physical equations; ii) parametric sweeping; iii) parallel computing; and iv) extensive dataset analysis, thereby enabling a complete “full-factorial” exploration of the design parameter space in combinatorial fashion. The simulations were conducted in a time-efficient manner without requiring massive computational time. As a case study, we simulated >15,000 microfluidic device designs and flow conditions for a representative multicellular spheroids-on-a-chip arrayed device, thus acquiring a single dataset consisting of ∼10 billion datapoints in ∼95 GBs. To validate our computational model, we performed physical experiments in a representative spheroid-on-a-chip device that showed excellent agreement between experimental and simulated data. This study offers a computational strategy to accelerate the optimization of microfluidic device designs and provide insight on the flow and drug transport in spheroid-on-a-chip and other biomicrofluidic platforms.

Introduction

Over the past few decades, microfluidics (MFs) has emerged as a powerful platform for fundamental and applied research in cell biology, soft robotics, medicine, drug screening, materials science, and analytical chemistry (Whitesides, 2006; Tian and Finehout, 2009; Young et al., 2012; Pak et al., 2015; Rajendra et al., 2015; Moore et al., 2017; Humayun et al., 2018; Khuu et al., 2019; Gevorkian et al., 2021). For all applications, the design of an MF device with optimized geometry is a vital first step of the engineering process, but the optimization process is usually conducted in a trial-and-error manner and is, therefore, time- and labour-intensive (Chakraborty, 2010; Huang, 2013). Many factors have to be considered and tested, such as the geometry and dimensions of microchannels or microwells, operating conditions such as fluid flow rates, and fluid properties. It is highly desirable to utilize a rational approach to the design process that would reduce the number of experimental trials and make it more time-efficient.

Computational modelling enables systematic exploration of a large parameter space via in silico simulations aimed at the optimization of the design and operation of an MF device prior to experiments (Enderling and Rejniak, 2013). In silico simulations have been extensively used for lab-on-a-chip platforms to simulate fluid flow and transport phenomena, including droplet formation (Sontti and Atta, 2017, 2020; Mohamed et al., 2019), micro-mixing (Zhang et al., 2008; Suh and Kang, 2010), bacteria growth and culture (Hohne et al., 2009; Westerwalbesloh et al., 2015; Kim et al., 2019; Kheiri et al., 2020), particle sorting and separation (Han et al., 2015; Lu and Xuan, 2015; Amin Arefi et al., 2020; Fallahi et al., 2020), biomechanical forces (Kim et al., 2015; Rousset et al., 2017), and drug delivery (Hossain et al., 2012; Soltani and Chen, 2012; Soltani et al., 2016; d’Esposito et al., 2018). In the specific case of MF devices for cell biology studies, computational fluid dynamics (CFD) was used to investigate molecular transport in a bilayer membrane-based MF device and examine the effect of flow-induced shear stress on endothelial cells (Wong et al., 2017). The CFD model provided guidelines for the design of MF devices operating within a range of physiological shear stresses while ensuring efficient transport of biomolecules through the membrane. In another computational study, a mathematical model coupled with numerical simulations led to an optimal design of MF devices for studies of endothelial cell migration and angiogenesis (Kuzmic et al., 2019). A computational framework was established to elucidate the effect of MF device geometry on cell migration and angiogenesis in a double-channel design with interconnected migration ports. Several studies have focused on the simulations of drug delivery to multicellular tumour spheroids (MCTSs) in a spheroid-on-a-chip platform (Kim et al., 2008; Moshksayan et al., 2018). More specifically, these simulations explored the effect of MF device geometry including microwell and microchannel dimensions on drug supply and uptake by MCTSs. Recently, a numerical study of the design of a tumour-on-a-chip platform was performed to determine its optimal performance in screening multidrug combinations (Hajari et al., 2021).

In the vast majority of numerical studies, the in-silico models are generally applied to a small set of conditions or parametric levels and employed only for “one-factor-at-a-time” (OFAT) experimental designs. For example, in modelling and simulating bacteria growth and glucose transport to bacterial colonies on-chip, the effect of glucose concentration was simulated while all other parameters were maintained invariant (Westerwalbesloh et al., 2015). Similarly, the effect of fluid flow and viscosity was modelled to study cancer cell loading in microwells of the MF device (Hajari et al., 2021). Inlet fluid flow rate and viscosity were varied systematically across three different values and thus nine parameter combinations were used overall, but the geometric characteristics of the MF device and MCTS properties were left constant to reduce the parameter space. Generally, computational models were treated as in silico analogs of experimental setups to simulate an experiment that could have been conducted using standard OFAT or factorial experimental designs. While such numerical studies hold great promise by i) simulating the output result without running an experiment and ii) providing a greater analytical power than is possible with experiments, these models would offer an even greater benefit if a larger parameter space could be explored. By definition, full factorial design (FFD) investigates the effects of all the selected factors and their interactions on the outcome(s) of the experiment (Das and Dewanjee, 2018). The term “full-factorial” refers to the fact that once the various factors of interest and the different levels within each factor have been selected, every conceivable combination of levels in the parameter space will be simulated and analyzed; this is distinct from “fractional factorial design,” where only a fraction or subset of the parameter space is tested and analyzed. Such simulations may reveal the most effective design and operating parameters that have not been considered or experimentally explored. For spheroid-on-a-chip MF platforms, the need for optimization is particularly important, as modifications of the geometry of the MF devices can significantly affect cell culture microenvironment and, subsequently, impact cell aggregation and MCTS structure and fate (Young and Beebe, 2010; Ong et al., 2017; Zuchowska et al., 2017). An approach in which MCTS culture and drug delivery conditions are “screened” in a time- and labour-efficient manner would be highly advantageous, as numerous operating conditions and geometrical constraints create a large parameter space in MF spheroid-on-a-chip platforms. Here, we report an efficient computational strategy for simulating fluid flow and transport of low-molecular-weight solutes (drugs) in an MF spheroid-on-a-chip platform (Wang et al., 2016). The platform is comprised of a large array of MCTSs compartmentalized in cylinder-shaped microwells connected with a common channel used for the continuous supply of cell culture medium and/or drug solution to the MCTSs (Ota et al., 2011; Schneider et al., 2013; Wang et al., 2016). We show that based on (>15,000) MF device designs and flow conditions acquired in a single dataset across a large parameter space (operating conditions, geometrical constraints, and different material porosity), these simulations can be effectively utilized to identify the most effective MF device design. The computational strategy combined and utilized four key factors: i) the governing physical equations; ii) parametric sweeping; iii) parallel computing; and iv) extensive dataset analysis, which all together enable a complete “full-factorial” exploration of the design parameter space in a combinatorial fashion. Importantly, the simulations were conducted in a time-efficient manner without requiring massive computational time. The computational model integrated multiple physical interfaces to examine advective and diffusive modes of drug delivery to MCTSs in MF devices with various geometries and examined the impact of MF device design on drug transport. This study demonstrates the power of the computational strategy and offers insights on drug transport phenomena in MF devices under a broad range of experimental conditions and MF device designs.

Materials and Methods

The overall workflow of the computational framework (Figure 1) consists of five major stages: i) definition of the MF device geometry or design selection, ii) construction of the in-silico model, which involves mesh generation, specification of the governing equations, and selection of the values or levels of all factors, iii) parametric sweeping across the entire parameter space, iv) parallel computing, and v) post-processing and analysis. The computational model was built using COMSOL Multiphysics 5.5 software (COMSOL, Stockholm, Sweden), a finite-element-based commercial CFD package. To capture the set of governing equations needed to model flow and molecular transport, two distinct COMSOL physics interfaces were coupled together, namely 1) the Laminar Flow (spf) interface and 2) the Transport of Diluted Species (tds) interface. An automatic solution-adaptive mesh refinement tool was used to regenerate the mesh and maintain mesh quality within an acceptable range for all simulations (Supplementary Material). The criterion for convergence in the simulations was selected to be a drop of four orders of magnitude (104).

FIGURE 1
www.frontiersin.org

FIGURE 1. Computational modelling workflow for microfluidic design optimization, involving: (1) design selection, (2) model construction including mesh generation, specification of governing equations, and selection of parametric levels, (3) parametric sweeping, (4) parallel computing, and (5) post-processing and analysis.

MF Device Geometry: Our Spheroid-on-a-Chip Model

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

FIGURE 2
www.frontiersin.org

FIGURE 2. Spheroid-on-a-chip device design and geometry. (A) 3D illustration of the array-based spheroid-on-chip microfluidic device used in this study. (B) Illustration of single supply channel (SSC) and double supply channel (DSC) designs of the spheroid-on-a-chip device. (C,D) Geometric parameters considered in the design optimization of the SSC and DSC spheroid-on-a-chip designs.

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

TABLE 1
www.frontiersin.org

TABLE 1. Geometrical, material, and operating conditions used in current spheroid-on-a-chip design study.

Governing Physical Equations

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:

ct=S+D2cuc(1)

where  c is the species concentration (mol/m3),  S is the supply influx (mol/m3⋅s),  u is the fluid velocity (m/s), and D is the diffusion coefficient of the molecular species (m2/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)
ρ[ut+uu]=p+μ2u+ρf(3)

where  ρ is the fluid density (kg/m3); 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/m3). 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)=Qm(4)

where εp is the porosity (or void fraction) of the matrix domain,  v is the velocity inside the matrix domain, and Qm is a mass source within the matrix domain (kg/m3·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 (m2) 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=Dc(6)

where  J is the diffusive flux (mol/m2⋅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:

ct+J+vc=R(7)

where R is the drug reaction rate (mol/m3⋅s). This reaction rate was set to zero in all our simulations.

Parametric Sweeping

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

Post-Processing and Analysis

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

Results and Discussion

Velocity and Pressure Distribution in a Single Row of Microwells

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.

FIGURE 3
www.frontiersin.org

FIGURE 3. Representative results of flow and transport simulations in COMSOL. (A) A single unit of the double-supply channel (DSC) design of the spheroid-on-a-chip device, consisting of the microwell for spheroid formation and culture, two supply channels, two connecting channels, and the spheroid itself (shown in yellow). (B,C) Velocity streamlines in the DSC device. (D) Velocity vectors showing fluid flow circulation overlaying streamlines of flow in the microwell. (E) Velocity magnitude of fluid flow on the surface of the spheroid, shown as arrow surface and surface contours. (F) Velocity vectors within the spheroid (treated as porous medium).

Experimental Validation of the Computational Model

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.

FIGURE 4
www.frontiersin.org

FIGURE 4. Doxorubicin (DOX) uptake by spheroid under two different flow rates (n = 5). (A) A fragment of the 100 μm-diameter arrays of MCF-7 breast cancer spheroids after loading. Scale bar is 1 mm. (B) Brightfield and fluorescence (NucGreen®/NucBlue) images of formed MCF-7 breast cancer spheroids after 72 h cell culture. Scale bar is 50 µm. (C) Immunofluorescence staining of tumour spheroids after 72 h of cell culture by E-cadherin antibody (green), DAPI (blue) and F-actin (red). Scale bar is 25 µm. (D) Fluorescence images of spheroids after 2 h perfusion of 10 µM DOX under different flow conditions. Scale bar is 100 µm. (E) The defined geometry and porous medium domain in COMSOL Multiphysics. (F) Results of drug penetration within the spheroid are visualized as arrows based on the total flux of molecules. (G) Results of total flux of DOX evaluated in the midplane of spheroid when Q = 0.01 ml/h and 0.02 ml/h.

Using the Computational Model as a Design Tool

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×105) 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.

Drug Permeation Depends on Spheroid Porosity

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 Q1 and Q2 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.

FIGURE 5
www.frontiersin.org

FIGURE 5. Total molecular flux and velocity magnitudes for different spheroid porosities. (A) A single unit of the double-supply channel (DSC) design of the spheroid-on-a-chip device, where the defined cross-section represents the sagittal yz-plane. (B) 3D contours of velocity magnitude and 3D vectors of drug flux around the spheroid when φ = 0.9. (C) Velocity magnitudes (2D surface plot) and vectors of drug flux within the sagittal yz-plane for φ = 0.2, 0.5, 0.7, and 0.9.

Drug Permeation Depends on Mode of Transport

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 (WR = HR = WL = HL = 25 μm, Q1 = Q2 = 0.01 ml/h, ε=0.5, Cdox = 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.

FIGURE 6
www.frontiersin.org

FIGURE 6. 3D volume contours of diffusion and advection fluxes for different microwell sizes of R = 80, 120, and 160 μm (WR = HR = WL = HL = 25 µm).

Drug Delivery Efficiency and Biomechanical Forces

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 (Q1 = Q2) (Figure 7A), and ii) unequal (Q1 = 2 × Q2) (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 HR = HL = WR = WL = 60 μm (maximum connecting channel size), iii) varying HR and HL independently (WR = WL = 60 μm), iv) varying WR and WL independently (HR = HL = 60 μm). The “reference” DSC design was chosen arbitrarily from the set of simulated designs to allow convenient comparison between different designs.

FIGURE 7
www.frontiersin.org

FIGURE 7. Concentration vs time profiles simulated at the spheroid core up to t = 120 min under different geometric conditions and spheroid porosities (φ) for (A) Q1 = Q2 and (B) Q1 = 2 × Q2. Legend: “SSC Design”—entire range of SSC designs; “Reference DSC Design”—HR = HL = WR = WL = 60 μm; “Variations in H”—WR = WL = 60 μm are constant and H varies independently; “Variations in W”—HR = HL = 60 μm are constant and W varies independently.

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 (HR and HL) showed a larger effect on drug delivery efficiency compared to changes in the width of the connecter (WR and WL). 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 (Q1 = 2 × Q2). 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 (Q1 = 2 × Q2) than for equal flow rates (Q1 = Q2). For example, when φ=0.6 and R = 120 μm, drug concentration at the spheroid core reached a plateau at 1×103mol/m3 at t = 60 min when Q1 = 2 × Q2, but drug concentration did not reach this plateau for Q1 = Q2, 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).

Post-Processing and Analysis: Using Design Heatmaps and GUIs as Visualization Tools

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 (Q1 = 2 × Q2) (Figure 8A) and equal supply flow rates (Q1 = Q2) (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.

FIGURE 8
www.frontiersin.org

FIGURE 8. Hierarchical clustering analysis of MF design space. Hierarchically clustered “design heatmaps” for (A) Q1 = 2×Q2 and (B) Q1 = Q2, show normalized values of concentration, shear rate, velocity, and advective and diffusive fluxes as five key metrics (columns). Independent device designs (rows) are designated by a naming convention, where P = Porosity, R = Chamber radius, Wr = Connector width—right, Hr = Connector height—right, Wl = Connector width–left, Hl = Connector height—left (e.g., P3. Wr6. R8.Hr6. Wl6. Hl6 denotes p = 0.3, Wr = 60 μm, R = 80 μm, Hr = 60 μm, WL = 60 μm, HL = 60 µm). Hierarchical clusters were computed with Euclidean distance and farthest point algorithm in Python 3.0 (Supplementary Material). (C,D) Illustrative examples of device design selection. Green boxes–specific design that achieves both low and high shear by changing flow rate ratio only, with all other metrics constant; Purple boxes–specific design that achieves both advection-dominant and diffusion-dominant drug penetration by changing flow rate ratio only, with all other metrics constant.

Hierarchical cluster analysis revealed that for both Q1 = 2 × Q2 and Q1 = Q2, 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.Wr6.Hr6.Wl6.Hl2) satisfied all these criteria (Figure 8C), since high shear stress can be achieved with Q1 = 2Q2 (Figure 8A) while low shear stress can be achieved with Q1 = Q2, 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.Wr.Hr6.Wl3.Hl2) could achieve advective-dominant transport when Q1 = Q2, 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 Q1 = 2 × Q2 (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 Mr > 1,000 is advection-dominant, as opposed to the transport of low molecular mass drugs with Mr < 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.

Conclusion

We have described a novel computational framework involving governing physical equations of fluid flow and transport, parametric sweeping, parallel computing, and post-processing and analysis that enables simultaneous simulation of >15,000 flow scenarios and designs of a spheroid-on-a-chip MF platform. While the approach is general and can be applied to any MF design, our study focused on one particular spheroid-on-a-chip configuration, which we used as a case study to illustrate the power and efficiency of our approach. The proposed simulation model has the potential to be extended to studies of other types of MF organ-on-a-chip devices. Moreover, the proposed framework will be helpful as a pipeline for more complex CFD modelling and optimization near future.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

SK: Conceptualization, data curation, formal analysis, investigation, methodology, writing–original draft, writing–review and editing. EK: Conceptualization, funding acquisition, supervision, resources, writing–review and editing. EY: Conceptualization, formal analysis, funding acquisition, methodology, project admin, resources, supervision, writing–review and editing.

Funding

We acknowledge the financial support from the Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grant (RGPIN-2019-5885) and a Dean’s Catalyst Professorship to EY, and financial support from the NSERC CREATE Training Program in Organ-on-a-Chip Engineering and Entrepreneurship (TOeP, 482073-2016) and an Ontario Graduate Scholarship to SK. We also acknowledge the support from CMC Microsystems for access to their CAD Compute Cluster. EK is grateful to Canada Research Chair (Tier 1) NSERC program.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

Some figures were created with BioRender.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2021.781566/full#supplementary-material

References

Afzal, A., Ansari, Z., Faizabadi, A. R., and Ramis, M. K. (2017). Parallelization Strategies for Computational Fluid Dynamics Software: State of the Art Review. Arch. Computat Methods Eng. 24, 337–363. doi:10.1007/s11831-016-9165-4

CrossRef Full Text | Google Scholar

Amin Arefi, S. M., Tony Yang, C. W., Sin, D. D., and Feng, J. J. (2020). Simulation of Nanoparticle Transport and Adsorption in a Microfluidic Lung-On-A-Chip Device. Biomicrofluidics. 14, 044117. doi:10.1063/5.0011353

PubMed Abstract | CrossRef Full Text | Google Scholar

Chakraborty, S. (2010). Microfluidics and Microfabrication. Editor S. Chakraborty (Boston, MA: Springer US). doi:10.1007/978-1-4419-1543-6

CrossRef Full Text | Google Scholar

Chary, S. R., and Jain, R. K. (1989). Direct Measurement of Interstitial Convection and Diffusion of Albumin in Normal and Neoplastic Tissues by Fluorescence Photobleaching. Proc. Natl. Acad. Sci. 86, 5385–5389. doi:10.1073/PNAS.86.14.5385

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Z., Kheiri, S., Gevorkian, A., Young, E. W. K., Andre, V., Deisenroth, T., et al. (2021). Microfluidic Arrays of Dermal Spheroids: a Screening Platform for Active Ingredients of Skincare Products. Lab. Chip. 21, 3952–3962. doi:10.1039/d1lc00619c

PubMed Abstract | CrossRef Full Text | Google Scholar

Cohen, D. E., Schneider, T., Wang, M., and Chiu, D. T. (2010). Self-digitization of Sample Volumes. Anal. Chem. 82, 5707–5717. doi:10.1021/ac100713u

PubMed Abstract | CrossRef Full Text | Google Scholar

Das, A. K., and Dewanjee, S. (2018). “Chapter 3 - Optimization of Extraction Using Mathematical Models and Computation,”, eds. S. D. Sarker, and L. B. T.-C. P. Nahar (Elsevier), 75–106. doi:10.1016/B978-0-12-812364-5.00003-1

CrossRef Full Text | Google Scholar

d’Esposito, A., Sweeney, P. W., Ali, M., Saleh, M., Ramasawmy, R., Roberts, T. A., et al. (2018). Computational Fluid Dynamics With Imaging of Cleared Tissue and of In Vivo Perfusion Predicts Drug Uptake and Treatment Responses in Tumours. Nat. Biomed. Eng. 2, 773–787. doi:10.1038/s41551-018-0306-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Dewhirst, M. W., and Secomb, T. W. (2017). Transport of Drugs From Blood Vessels to Tumour Tissue. Nat. Rev. Cancer. 17, 738–750. doi:10.1038/nrc.2017.93

PubMed Abstract | CrossRef Full Text | Google Scholar

Enderling, H., and Rejniak, K. A. (2013). Simulating Cancer: Computational Models in Oncology. Front. Oncol. 3, 13–14. doi:10.3389/fonc.2013.00233

PubMed Abstract | CrossRef Full Text | Google Scholar

Fallahi, H., Zhang, J., Nicholls, J., Phan, H.-P., and Nguyen, N.-T. (2020). Stretchable Inertial Microfluidic Device for Tunable Particle Separation. Anal. Chem. 92, 12473–12480. doi:10.1021/acs.analchem.0c02294

PubMed Abstract | CrossRef Full Text | Google Scholar

Gevorkian, A., Morozova, S. M., Kheiri, S., Khuu, N., Chen, H., Young, E., et al. (2021). Actuation of Three-Dimensional-Printed Nanocolloidal Hydrogel With Structural Anisotropy. Adv. Funct. Mater. 31, 2010743. doi:10.1002/adfm.202010743

CrossRef Full Text | Google Scholar

Hajari, M. A., Baheri Islami, S., and Chen, X. (2021). A Numerical Study on Tumor-On-Chip Performance and its Optimization for Nanodrug-Based Combination Therapy. Biomech. Model. Mechanobiol. 20, 983–1002. doi:10.1007/s10237-021-01426-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, X., Feng, Y., Cao, Q., and Li, L. (2015). Three-Dimensional Analysis and Enhancement of Continuous Magnetic Separation of Particles in Microfluidics. Microfluid. Nanofluid. 18, 1209–1220. doi:10.1007/s10404-014-1516-6

CrossRef Full Text | Google Scholar

Hohne, D. N., Younger, J. G., and Solomon, M. J. (2009). Flexible Microfluidic Device for Mechanical Property Characterization of Soft Viscoelastic Solids Such as Bacterial Biofilms. Langmuir. 25, 7743–7751. doi:10.1021/la803413x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hossain, S. S., Hossainy, S. F. A., Bazilevs, Y., Calo, V. M., and Hughes, T. J. R. (2012). Mathematical Modeling of Coupled Drug and Drug-Encapsulated Nanoparticle Transport in Patient-Specific Coronary Artery walls. Comput. Mech. 49, 213–242. doi:10.1007/s00466-011-0633-2

CrossRef Full Text | Google Scholar

Huang, H., Yu, Y., Hu, Y., He, X., Berk Usta, O., and Yarmush, M. L. (2017). Generation and Manipulation of Hydrogel Microcapsules by Droplet-Based Microfluidics for Mammalian Cell Culture. Lab. Chip. 17, 1913–1932. doi:10.1039/C7LC00262A

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, T. J. (2013). Introduction to BioMEMS, by Albert Folch . CRC Press, Boca Raton, FL, 2012, 528 Pages. ISBN 978-1-43-981839-8 (Hardcover). Microsc. Microanal. 19, 506. doi:10.1017/S1431927613000081

CrossRef Full Text | Google Scholar

Humayun, M., Chow, C.-W., and Young, E. W. K. (2018). Microfluidic Lung Airway-On-A-Chip With Arrayable Suspended Gels for Studying Epithelial and Smooth Muscle Cell Interactions. Lab. Chip. 18, 1298–1309. doi:10.1039/C7LC01357D

PubMed Abstract | CrossRef Full Text | Google Scholar

Hwang, C. M., Sant, S., Masaeli, M., Kachouie, N. N., Zamanian, B., Lee, S.-H., et al. (2010). Fabrication of Three-Dimensional Porous Cell-Laden Hydrogel for Tissue Engineering. Biofabrication 2, 035003. doi:10.1088/1758-5082/2/3/035003

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaminski, T. S., and Garstecki, P. (2017). Controlled Droplet Microfluidic Systems for Multistep Chemical and Biological Assays. Chem. Soc. Rev. 46, 6210–6226. doi:10.1039/c5cs00717h

PubMed Abstract | CrossRef Full Text | Google Scholar

Kheiri, S., Mohamed, M. G. A., Amereh, M., Roberts, D., and Kim, K. (2020). Antibacterial Efficiency Assessment of Polymer-Nanoparticle Composites Using a High-Throughput Microfluidic Platform. Mater. Sci. Eng. C. 111, 110754. doi:10.1016/j.msec.2020.110754

PubMed Abstract | CrossRef Full Text | Google Scholar

Khuu, N., Alizadehgiashi, M., Gevorkian, A., Galati, E., Yan, N., and Kumacheva, E. (2019). Temperature-Mediated Microfluidic Extrusion of Structurally Anisotropic Hydrogels. Adv. Mater. Technol. 4, 1800627. doi:10.1002/admt.201800627

CrossRef Full Text | Google Scholar

Kim, J.-Y., Fluri, D. A., Marchan, R., Boonen, K., Mohanty, S., Singh, P., et al. (2015). 3D Spherical Microtissues and Microfluidic Technology for Multi-Tissue Experiments and Analysis. J. Biotechnol. 205, 24–35. doi:10.1016/j.jbiotec.2015.01.003

CrossRef Full Text | Google Scholar

Kim, M.-C., Wang, Z., Lam, R. H. W., and Thorsen, T. (2008). Building a Better Cell Trap: Applying Lagrangian Modeling to the Design of Microfluidic Devices for Cell Biology. J. Appl. Phys. 103, 044701. doi:10.1063/1.2840059

CrossRef Full Text | Google Scholar

Kim, S., Lee, S., Kim, J.-K., Chung, H. J., and Jeon, J. S. (2019). Microfluidic-Based Observation of Local Bacterial Density Under Antimicrobial Concentration Gradient for Rapid Antibiotic Susceptibility Testing. Biomicrofluidics 13, 014108. doi:10.1063/1.5066558

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuzmic, N., Moore, T., Devadas, D., and Young, E. W. K. (2019). Modelling of Endothelial Cell Migration and Angiogenesis in Microfluidic Cell Culture Systems. Biomech. Model. Mechanobiol. 18, 717–731. doi:10.1007/s10237-018-01111-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, K., Kim, C., Young Yang, J., Lee, H., Ahn, B., Xu, L., et al. (2012). Gravity-Oriented Microfluidic Device for Uniform and Massive Cell Spheroid Formation. Biomicrofluidics 6, 014114. doi:10.1063/1.3687409

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Khuu, N., Prince, E., Tao, H., Zhang, N., Chen, Z., et al. (2021). Matrix Stiffness-Regulated Growth of Breast Tumor Spheroids and Their Response to Chemotherapy. Biomacromolecules 22, 419–429. doi:10.1021/acs.biomac.0c01287

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, W., Wang, J.-C., and Wang, J. (2015). Controllable Organization and High Throughput Production of Recoverable 3D Tumors Using Pneumatic Microfluidics. Lab. Chip. 15, 1195–1204. doi:10.1039/c4lc01242a

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, X., and Xuan, X. (2015). Continuous Microfluidic Particle Separation via Elasto-Inertial Pinched Flow Fractionation. Anal. Chem. 87, 6389–6396. doi:10.1021/acs.analchem.5b01432

PubMed Abstract | CrossRef Full Text | Google Scholar

McMillan, K. S., Boyd, M., and Zagnoni, M. (2016). Transitioning From Multi-phase to Single-phase Microfluidics for Long-Term Culture and Treatment of Multicellular Spheroids. Lab. Chip. 16, 3548–3557. doi:10.1039/c6lc00884d

PubMed Abstract | CrossRef Full Text | Google Scholar

Mehta, G., Hsiao, A. Y., Ingram, M., Luker, G. D., and Takayama, S. (2012). Opportunities and Challenges for Use of Tumor Spheroids as Models to Test Drug Delivery and Efficacy. J. Controlled Release 164, 192–204. doi:10.1016/j.jconrel.2012.04.045

CrossRef Full Text | Google Scholar

Minchinton, A. I., and Tannock, I. F. (2006). Drug Penetration in Solid Tumours. Nat. Rev. Cancer 6, 583–592. doi:10.1038/nrc1893

PubMed Abstract | CrossRef Full Text | Google Scholar

Mohamed, M. G. A., Kheiri, S., Islam, S., Kumar, H., Yang, A., and Kim, K. (2019). An Integrated Microfluidic Flow-Focusing Platform for On-Chip Fabrication and Filtration of Cell-Laden Microgels. Lab. Chip. 19, 1621–1632. doi:10.1039/C9LC00073A

PubMed Abstract | CrossRef Full Text | Google Scholar

Moore, T. A., Brodersen, P., and Young, E. W. K. (2017). Multiple Myeloma Cell Drug Responses Differ in Thermoplastic vs PDMS Microfluidic Devices. Anal. Chem. 89, 11391–11398. doi:10.1021/acs.analchem.7b02351

PubMed Abstract | CrossRef Full Text | Google Scholar

Moshksayan, K., Kashaninejad, N., Warkiani, M. E., Lock, J. G., Moghadas, H., Firoozabadi, B., et al. (2018). Spheroids-on-a-Chip: Recent Advances and Design Considerations in Microfluidic Platforms for Spheroid Formation and Culture. Sensors Actuators B: Chem. 263, 151–176. doi:10.1016/j.snb.2018.01.223

CrossRef Full Text | Google Scholar

Ong, L. J. Y., Islam, A., DasGupta, R., Iyer, N. G., Leo, H. L., and Toh, Y.-C. (2017). A 3D Printed Microfluidic Perfusion Device for Multicellular Spheroid Cultures. Biofabrication 9, 045005. doi:10.1088/1758-5090/aa8858

PubMed Abstract | CrossRef Full Text | Google Scholar

Ota, H., Kodama, T., and Miki, N. (2011). Rapid Formation of Size-Controlled Three Dimensional Hetero-Cell Aggregates Using Micro-Rotation Flow for Spheroid Study. Biomicrofluidics 5, 034105. doi:10.1063/1.3609969

PubMed Abstract | CrossRef Full Text | Google Scholar

Ota, H., Yamamoto, R., Deguchi, K., Tanaka, Y., Kazoe, Y., Sato, Y., et al. (2010). Three-Dimensional Spheroid-Forming Lab-On-A-Chip Using Micro-Rotational Flow. Sensors Actuators B: Chem. 147, 359–365. doi:10.1016/j.snb.2009.11.061

CrossRef Full Text | Google Scholar

Pak, C., Callander, N. S., Young, E. W. K., Titz, B., Kim, K., Saha, S., et al. (2015). MicroC3: an Ex Vivo Microfluidic Cis-Coculture Assay to Test Chemosensitivity and Resistance of Patient Multiple Myeloma Cells. Integr. Biol. 7, 643–654. doi:10.1039/C5IB00071H

CrossRef Full Text | Google Scholar

Patra, B., Peng, C.-C., Liao, W.-H., Lee, C.-H., and Tung, Y.-C. (2016). Drug Testing and Flow Cytometry Analysis on a Large Number of Uniform Sized Tumor Spheroids Using a Microfluidic Device. Sci. Rep. 6, 1–12. doi:10.1038/srep21061

PubMed Abstract | CrossRef Full Text | Google Scholar

Prince, E., Kheiri, S., Wang, Y., Xu, F., Cruickshank, J., Topolskaia, V., et al. (2021). Microfluidic Arrays of Breast Tumor Spheroids for Drug Screening and Personalized Cancer Therapies. Adv. Healthc. Mater., 2101085. doi:10.1002/adhm.202101085

CrossRef Full Text | Google Scholar

Rajendra, V., Therien-Aubin, H., Abolhasani, M., Villalabos, M., and Kumacheva, E. (2015). An Exploratory Microfluidic Approach to Photopolymerized Polymer-Inorganic Nanocomposite Films. Macromol. Mater. Eng. 300, 1071–1078. doi:10.1002/mame.201500221

CrossRef Full Text | Google Scholar

Rousset, N., Monet, F., and Gervais, T. (2017). Simulation-assisted Design of Microfluidic Sample Traps for Optimal Trapping and Culture of Non-Adherent Single Cells, Tissues, and Spheroids. Sci. Rep. 7, 245. doi:10.1038/s41598-017-00229-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Sabhachandani, P., Motwani, V., Cohen, N., Sarkar, S., Torchilin, V., and Konry, T. (2016). Generation and Functional Assessment of 3D Multicellular Spheroids in Droplet Based Microfluidics Platform. Lab. Chip. 16, 497–505. doi:10.1039/C5LC01139F

PubMed Abstract | CrossRef Full Text | Google Scholar

Schneider, T., Yen, G. S., Thompson, A. M., Burnham, D. R., and Chiu, D. T. (2013). Self-Digitization of Samples Into a High-Density Microfluidic Bottom-Well Array. Anal. Chem. 85, 10417–10423. doi:10.1021/ac402383n

PubMed Abstract | CrossRef Full Text | Google Scholar

Soltani, M., and Chen, P. (2012). Effect of Tumor Shape and Size on Drug Delivery to Solid Tumors. J. Biol. Eng. 6, 4. doi:10.1186/1754-1611-6-4

CrossRef Full Text | Google Scholar

Soltani, M., Sefidgar, M., Bazmara, H., Marcus, C., Subramaniam, R. M., and Rahmim, A. (2016). Effect of Tumor Size on Drug Delivery to Lung Tumors. IEEE Nucl. Sci. Symp. Med. Imaging Conf. Nss/mic. 2015, 1–15. doi:10.1109/NSSMIC.2015.7582238

CrossRef Full Text | Google Scholar

Sontti, S. G., and Atta, A. (2017). CFD Analysis of Microfluidic Droplet Formation in Non-Newtonian Liquid. Chem. Eng. J. 330, 245–261. doi:10.1016/j.cej.2017.07.097

CrossRef Full Text | Google Scholar

Sontti, S. G., and Atta, A. (2020). Numerical Insights on Controlled Droplet Formation in a Microfluidic Flow-Focusing Device. Ind. Eng. Chem. Res. 59, 3702–3716. doi:10.1021/acs.iecr.9b02137

CrossRef Full Text | Google Scholar

Steuperaert, M., Falvo D’Urso Labate, G., Debbaut, C., De Wever, O., Vanhove, C., Ceelen, W., et al. (2017). Mathematical Modeling of Intraperitoneal Drug Delivery: Simulation of Drug Distribution in a Single Tumor Nodule. Drug Deliv. 24, 491–501. doi:10.1080/10717544.2016.1269848

PubMed Abstract | CrossRef Full Text | Google Scholar

Suh, Y. K., and Kang, S. (2010). A Review on Mixing in Microfluidics. Micromachines. 1, 82–111. doi:10.3390/mi1030082

CrossRef Full Text | Google Scholar

Tian, W. C., and Finehout, E. (2009). Microfluidics for Biological Applications. Boston, MA: Springer US. doi:10.1007/978-0-387-09480-9

CrossRef Full Text | Google Scholar

Wang, Y., Li, Y., Thérien-Aubin, H., Ma, J., Zandstra, P. W., and Kumacheva, E. (2016). Two-dimensional Arrays of Cell-Laden Polymer Hydrogel Modules. Biomicrofluidics 10, 014110. doi:10.1063/1.4940430

PubMed Abstract | CrossRef Full Text | Google Scholar

Westerwalbesloh, C., Grünberger, A., Stute, B., Weber, S., Wiechert, W., Kohlheyer, D., et al. (2015). Modeling and CFD Simulation of Nutrient Distribution in Picoliter Bioreactors for Bacterial Growth Studies on Single-Cell Level. Lab. Chip. 15, 4177–4186. doi:10.1039/C5LC00646E

PubMed Abstract | CrossRef Full Text | Google Scholar

Whitesides, G. M. (2006). The Origins and the Future of Microfluidics. Nature 442, 368–373. doi:10.1038/nature05058

PubMed Abstract | CrossRef Full Text | Google Scholar

Wong, J. F., Young, E. W. K., and Simmons, C. A. (2017). Computational Analysis of Integrated Biosensing and Shear Flow in a Microfluidic Vascular Model. AIP Adv. 7, 115116. doi:10.1063/1.5006655

CrossRef Full Text | Google Scholar

Young, E. W. K., and Beebe, D. J. (2010). Fundamentals of Microfluidic Cell Culture in Controlled Microenvironments. Chem. Soc. Rev. 39, 1036. doi:10.1039/b909900j

PubMed Abstract | CrossRef Full Text | Google Scholar

Young, E. W. K., Pak, C., Kahl, B. S., Yang, D. T., Callander, N. S., Miyamoto, S., et al. (2012). Microscale Functional Cytomics for Studying Hematologic Cancers. Blood 119, e76–e85. doi:10.1182/blood-2011-10-384347

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Z., Yim, C., Lin, M., and Cao, X. (2008). Quantitative Characterization of Micromixing Simulation. Biomicrofluidics 2, 034104. doi:10.1063/1.2966454

PubMed Abstract | CrossRef Full Text | Google Scholar

Zuchowska, A., Jastrzebska, E., Zukowski, K., Chudy, M., Dybko, A., and Brzozka, Z. (2017). A549 and MRC-5 Cell Aggregation in a MicrofluidicLab-On-A-Chipsystem. Biomicrofluidics 11, 024110. doi:10.1063/1.4979104

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: organ-on-a-chip, cancer spheroids, drug delivery, full-factorial experiments, microfluidic design space exploration, hierarchical clustering

Citation: Kheiri S, Kumacheva E and Young EWK (2021) Computational Modelling and Big Data Analysis of Flow and Drug Transport in Microfluidic Systems: A Spheroid-on-a-Chip Study. Front. Bioeng. Biotechnol. 9:781566. doi: 10.3389/fbioe.2021.781566

Received: 22 September 2021; Accepted: 27 October 2021;
Published: 23 November 2021.

Edited by:

Z. Hugh Fan, University of Florida, United States

Reviewed by:

Andrzej Przekwas, CFD Research Corporation, United States
Ting-Yuan Tu, National Cheng Kung University, Taiwan

Copyright © 2021 Kheiri, Kumacheva and Young. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Edmond W.K. Young, eyoung@mie.utoronto.ca

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.