A parametric study of mechanoporation through microfluidic design to modulate shear, compressive, and adhesion forces and loading rates

Abstract

Efficient and reproducible intracellular delivery is critical for manufacturing next generation cell therapies. Mechanoporation employs mechanical forces, including shear loading, adhesion, and compressive strain, to transiently permeabilize cell membranes and enable cargo transport. However, the influence of microsecond-scale unsteady forces and the origins of variability in delivery and viability remain insufficiently characterized. Here, we performed a parametric investigation of microfluidic mechanoporation using parallelized channel designs of varied widths to systematically modulate pre-compression shear loading and strain rates under constant volumetric flow. Narrow channels were found to promote a more uniform pre-constriction loading and compressive dynamics, leading to improved reproducibility of delivery outcomes. High-speed video analysis revealed greater cell focusing and computational fluid dynamics (CFD) confirmed higher pre-constriction shear loading rates and higher asymmetric biaxial forces prior to ridges, yielding a substantial improvement in delivery efficiency in both adherent B16F10 melanoma cells and suspended T-cells. Modulating cell–surface adhesion by adjusting surface chemistry showed that adhesive coatings slightly increase delivery efficiency at the expense of viability. Changing cell stiffness with pharmacological softening caused a decline in delivery efficiency. These trends indicate that mechanoporation outcomes are governed more strongly by the kinetics of loading dictated by fluid-driven acceleration and strain rate rather than by absolute strain or adhesion magnitude. Principal component and multivariate analyses identified two significant predictors of delivery and viability: strain rate and Basset–Boussinesq history (BBH) forces. Both predictors were consistently elevated in narrow multichannel architectures that showed higher delivery and lower viability. Together, these findings demonstrate that narrow channel designs establish a geometry-driven acceleration regime characterized by elevated strain rates and BBH forces that enhances delivery efficiency while imposing a viability tradeoff.

---

Introduction

Mechanoporation has emerged as a promising technique for intracellular delivery of diverse cargo such as messenger RNA (mRNA),^1,2 DNA,³ CRISPR/Cas9,^3,4 nanoparticles,⁵ and radiotracers.⁶ Rapid progress in mechanoporation has led to the development of multiple microfluidic strategies that induce transient membrane permeabilization through mechanical deformation. However, a significant knowledge gap remains in understanding how microfluidic channel geometry and ultra-fast mechanical loading govern cargo delivery and cell viability. This challenge is compounded by heterogeneity in cell size, mechanical properties, and trajectories, resulting in broad distributions of loading conditions even within nominally identical experiments. Microfluidic mechanoporation platforms may also experience channel clogging due to cell aggregation, debris accumulation, or adhesion to constrictions, particularly under highly confined or high-throughput conditions. Such clogging disrupts continuous operation and introduces uncontrolled, time-dependent changes in local flow fields and loading conditions, further complicating interpretation and reproducibility of mechanoporation outcomes. Consequently, mechanoporation performance remains difficult to predict or standardize across devices and operating regimes. Few systematic studies have examined how microfluidic geometry modulates both the magnitude and kinetics of mechanical loading such as strain, strain rate, and hydrodynamic forces and how these parameters relate to successful delivery versus irreversible cell damage. Improving how these loading metrics are defined and quantified is therefore critical for advancing mechanoporation platforms.

Existing mechanoporation techniques include microinjection,^7,8 shear stretching,^9,10 compression of cells through micro-constrictions,^5,11,12 and viscoelastic stretching.^13,14 Each method demonstrates specific advantages to improve delivery efficacy, reduce damage to cells and improve the overall production of modified cells for applications such as cell-based therapies. These strategies can be broadly categorized by their reliance on different force paradigms. Methods based on cell stretching, vortex shedding¹⁵ and hydrodynamic extensional forces¹⁶ leverage shear to facilitate diffusion-dominated delivery. These forces are applied over 100 s of microseconds to 10s of milliseconds generating shear rates exceeding 10 000 s⁻¹. Typical shear stresses¹⁷ range from 0.1 to 10 kPa. In contrast, compression devices^4–6,12 achieve high delivery efficiency through rapid compressive deformation ranging from microseconds to milliseconds that can cause convective transport of cargo. While these approaches are often described as compression-dominated, cells in such geometries are also subject to rapid acceleration, deceleration, and extensional components of deformation. Limited literature has studied extensional and compressive deformation modes independently.^18–25 Their coupled, transient manifestation during mechanoporation remains poorly resolved. Molecular dynamics simulations indicate that increasing membrane stretching rate promotes multi-pore formation²¹ and that symmetric biaxial stretching reduces yield strain thereby leading to early onset of rupture²⁶ compared to asymmetric loading. Viscoelastic additives^2,14,27 can increase membrane shear loading rates and cargo exchange,¹⁴ however, direct correlations between applied forces, deformation kinetics, and mechanoporation outcomes in confined microchannels remain limited.

Despite the importance of loading kinetics, most mechanoporation studies estimate hydrodynamic forces using steady-state metrics, often neglecting unsteady contributions to reduce computational complexity. Experimentally, resolving such unsteady effects is further constrained by the spatial and temporal resolution of high-speed imaging. Basset–Boussinesq history (BBH) forces capture the time-integrated viscous resistance associated with rapid changes in slip velocity, providing a physically grounded metric for unsteady hydrodynamic loading not captured by steady shear stress. While BBH forces are well-established in particle-laden flow theory, their relevance in biological microfluidic systems remains largely unexplored.^28–30 Because BBH depends on the time history of slip acceleration, confinement-induced distortion of local velocity fields due to the passing cells or blocked cells, is expected to amplify BBH contributions in smaller channels, even when bulk flow remains laminar. In this study, we examine whether BBH-related loading metrics correlate with mechanoporation performance under conditions of rapid acceleration and confinement.

Previously, we have reported a unique platform for promoting convective cargo exchange by employing high strain rates through geometric confinement.^4,12,31 Here, we build on this work by systematically varying channel width to modulate flow acceleration, strain rate, and unsteady hydrodynamic loading while directly recording cell trajectories and delivery outcomes (Fig. 1(a)). By reducing width (denoted as W in Fig. 1(b)), we increase the particle-to-channel size ratio, thereby enhancing geometric confinement. In confined laminar flows with finite-sized particles, higher blockage ratios distort local velocity profiles and sharpen ridge-entry velocity gradients.³² We therefore hypothesized that narrower channels would elevate strain rates as well as history-dependent hydrodynamic loading (BBH), leading to improved delivery efficiency but potentially increased viability loss.

Fig. 1 Overview of design changes and force framework for μfluidic channels with varying widths. (a) Experimental setup for mechanoporation showing syringe pump-driven flow. (b) Channel architectures: single-channel (C01, 500 μm width (W), 25 μm height (H_c)), five-channel array (C05, W = 100 μm, H_c = 25 μm), and ten-channel array (C10, W = 50 μm width, H_c = 25 μm). Chevron-shaped ridges (45°) were patterned within each channel. Designs were chosen to preserve comparable theoretical flow resistance across devices. (c) Center schematic shows the breakdown of forces acting on cells as they transit through the channel. Left inset (pre-ridge) shows a zoomed version of membrane stretching and loading changes at the onset of compression due to rapid acceleration of fluid result in boundary layer interactions and right inset (under-ridge) shows combined shear, adhesion and compression force behaviour when the call is undergoing compression.

To interpret the physical processes underlying cell deformation in these devices, we developed a force-based framework that segments cell-geometry interactions into three operational stages (Fig. 1(c)). (i) As cells approach ridge constrictions, rapid changes in local fluid velocity generate transient hydrodynamic loading associated with acceleration and shear gradients. (ii) Upon entering the ridge region, cells experience rapid compressive deformation, accompanied by adhesion and inertial effects. (iii) Finally, cells exit the ridge region while recovering from these rapid forces and compressions exchanging cargo using forced convection.^4,12

Among several trials, a multivariate analysis of flow dynamics and mechanoporation was conducted using viability, delivery, recovery as well as staining index as metrics. High-speed imaging and computational fluid dynamics (CFD) simulations were employed to track cell trajectories and correlate these with local mechanical environments. Together, these design variants and mechanistic considerations form the foundation for analysing how channel width influences mechanoporation performance while highlighting the role of loading kinetics and unsteady hydrodynamic effects.

Materials and methods

Microfluidic channel design

Microfluidic devices with different number of channels were designed using LayoutEditor®. The device architectures comprised of a single wide channel (C01), a five-channel array (C05), and a ten-channel array (C10) (Fig. 1(b)). Channel widths were adjusted such that the total cross-sectional area available for flow was comparable across designs, enabling direct translation of flow conditions and performance metrics between workflows. Inlet and outlet branching networks were analytically modelled to obtain uniform flow partitioning and matched mean fluid velocity across parallel channels, following established design principles.³³ Similar hydrodynamic flow resistance across all designs was validated with incompressible laminar flow simulations at multiple flow rates using Simscale CFD®. Reynnold's numbers for all channels and flow rates were computed to ensure laminar flow behaviour in all channels (Table S1). To ensure consistent mechanical loading, gutters were removed from ridge regions so that all cells experienced compressive deformation at the low gap heights.⁴ This design choice also negated the need for trajectory-focusing serpentine inlet geometries. The number of ridges (6 per channel) was also fixed across all designs based on prior optimization of the cell-VECT platform.⁴ All ridges were chevron-shaped and oriented at 45° relative to the flow direction to promote lateral focusing toward the channel center while permitting displacement of cell aggregates toward the channel edges.³¹

Microfluidic device fabrication

Silicon wafers of diameter 4″ were used as permanent molds for replica molding poly-dimethyl siloxane (PDMS) (Sylgard™ 184 Elastomer, Dow Corporation, Corning, NY). The mold was fabricated using a 2-step photolithography process using SU8. Final mold dimensions were characterized using profilometry. To fabricate the devices, uncured PDMS was mixed with curing agent in the ratio 10 : 1 and cured on the wafer for 90 minutes at 85 °C. Cured PDMS was peeled from the silicon wafers and 1 mm holes were punched at the inlets and outlets. Finally, the devices were cut and bonded to a microscopic glass slide using plasma treatment (Harrick Plasma, Ithaca, NY) for 50 seconds and baking the device for 60 minutes at 85 °C.

For varying surface adhesion interactions, the fabricated devices were plasma treated for 50 seconds and then coated with a solution of 2% v/v APTES in ethanol or 1% v/v Pluronic in PBS to simulate adhesive and non-adhesive conditions, respectively. The coated devices were filled with PBS and stored at 4 °C for use within 12 hours. For most optimization studies, the surface was not passivated with any coating.

Computational fluid dynamics simulations for force calculations

Three-dimensional fluid flow simulations were performed in COMSOL Multiphysics using the laminar flow interface. The steady-state incompressible Navier–Stokes equations were solved for a Newtonian fluid, assuming constant density (ρ = 1000 kg m⁻³) and viscosity (μ = 0.001 Pa s).

Isolated channels were modelled with a stationery particle of diameter 10 μm located approximately 1 μm upstream from the ridge. The inlet velocity was prescribed as a constant value calculated from the total volumetric flow rate of the device, and a zero-pressure condition was applied at the outlet.

The resulting stationary problem was solved using an algebraic multigrid (AMG) iterative solver using the GMRES method with a restart after 50 iterations and left preconditioning. Krylov subspace recycling was enabled through the GCRO-DP method with 25 eigenvectors. The relative convergence tolerance for the nonlinear stationary solver was set to 10⁻³. The resultant hydrodynamic force was calculated as F_i = ∫∫_Spσ_ij·n_jdS_p where F is the drag force, S_p, denotes the particle surface, n is the unit normal vector to the particle surface (pointing into the fluid), and dS_p is the differential surface element. Spatial discretization was performed using a mixed finite-element formulation with quadratic velocity and linear pressure shape functions (P2–P1). A boundary-layer mesh was applied along solid boundaries to improve near-wall resolution and to facilitate solver convergence. To minimize the grid and discretization errors in particle force calculation, the “fine” mesh scheme was chosen (Table 1).

Table 1 Mesh independence analysis. The mesh factor indicates the fraction of particle radius used to discretize the particle surface

Scheme	Max (μm)	Min (μm)	Mesh factor	Mesh count	F x (nN)
Coarse	10	3	0.3	204 055	943.6
Coarse	10	3	0.2	214 794	968.3
Coarse	10	3	0.1	268 815	989.69
Coarse	10	3	0.05	456 939	997.41
Normal	6.8	2	0.05	686 898	990.03
Fine	5.4	1	0.05	1 272 964	996.87
Coarse	10	3	0.02	1 618 987	998.29

---

PDMS device deformation characterization

PDMS : curing agent ratio of 10 : 1 and curing temperature of 85 °C was used to create multiple devices. PBS was flown through fabricated devices using an Elveflow pressure controller combined with a flow sensor. Multiple fabrication replicates of devices were tested. Next, experimental flow rate (Q) and pressure drop (ΔP) measurements were imported into R to isolate average flow rate at a chosen pressure. Experimental ΔP–Q relationships were used to fit Gervais model³⁴ over data identifying two different parameters α_c and α_r. α_c was used to determine deformation in regions of channel with height H_c and α_r in regions under constrictions with height H. Goodness of fit was analyzed using root mean square error (RMSE) and Bayesian information criterion (BIC).

Cell culture

B16F10 cells were purchased from ATCC and cultured in DMEM (Gibco, Thermo Fisher Scientific, Waltham, MA) with 10% fetal bovine serum (FBS), 1% penicillin–streptomycin and 1% L-glutamine supplement as GlutaMAX (Gibco, Thermo Fisher Scientific, Waltham, MA). For each sub-culture, trypsin-0.25% EDTA was added to the culture flask and cells were allowed to lift off for 15 minutes in an incubator with 5% CO₂ at 37 °C. The cells were then diluted in culture media to obtain a final dilution ratio of 1 : 5. For primary T-cells, PBMCs from healthy consenting donors were purchased from AllCells (Alameda, CA) and cryopreserved. Upon thawing, T-cells were isolated with EasySep Human T-cell isolation kit (STEM CELL technologies, Vancouver, Canada) and stimulated with Dynabead Human T-Activator CD3/CD28 (Thermo Fisher Scientific, Waltham, MA) in a 1 : 1 ratio, according to the manufacturer protocol. T cells were expanded in X-VIVO 15 media (Lonza, Greenwood, SC) with 10% FBS (Corning, MA), 1% penicillin–streptomycin (Lonza, Greenwood, SC), 1% GlutaMAX (Gibco, Thermo Fisher Scientific, Waltham, MA) and 100 ng mL⁻¹ IL-2 (PeptroTech, Rocky Hill, NJ). Cells were maintained at a concentration of 1.0 × 10⁶ cells per mL and incubated at 37 °C with 5% CO₂. 24-Hours after stimulation, Dynabeads® were removed by pipetting and brief incubation on DynaMag® magnet. T-cells were reactivated every 6 days from the last activation with a 24-hour Dynabead stimulation and delivery experiments were conducted on day-10 after isolation.

Delivery of FITC-dextran

Cells were counted from a culture sample using NucleoCounter NC200 while recording cell concentrations and viabilities. Cells were then resuspended in serum-free media Opti-MEM™ (Gibco, ThermoFisher Scientific, Waltham, MA) at 2.0 × 10⁶ cells per mL for FITC dextran delivery, alongside 0.4 mg mL⁻¹ 250 kDa FITC-conjugated Dextran (Sigma-Aldrich). The final suspension was passed through a 20 μm cell strainer (PluriSelect, El Cajon, CA) cell strainer for B16F10 and a 15 μm cell strainer for T-cells to ensure a homogenous suspension. The suspension was stored on ice before running through the microfluidic channels. Prior to processing the cells, the channels were cleaned with PBS to ensure the absence of debris and air bubbles in the channel.

We used a syringe pump (Harvard Systems, MA) to infuse the cell suspension into the microfluidic device at the desired flow rate. After collecting the cells from the outlets, they were resuspended in 10× volume of DPBS(+/+) and washed twice by centrifuging at 1000 RPM for 10 minutes. Once washed, we resuspended the cells in DPBS(+/+) to quantify their fluorescence. The collected samples were split into two parts and viability dye 7-AAD was added to one of the samples.

Flow cytometry

Fluorescence was measured using a FACS Melody analyser (BD Biosciences, CA, USA) and a Cytoflex S Cytometer (Beckman Coulter, CA, USA). Before samples were run through the cytometer, a QC cycle was run each time with beads and a compensation cycle was run including unstained cells, cells with FITC molecules in suspension and, thermal shock induced dead cells with 7-AAD. Cytometry data was analysed using FlowJo to identify positively and negatively stained populations. Viability was calculated as the percentage of cells with a negative stain when incubated with a membrane-impermeable label 7-AAD and PI. Delivery efficiency was calculated as the percentage of viable cells that shifted to a higher FITC intensity region compared to an unprocessed control. Staining index was calculated as where Median⁺ represents the median FITC intensity for stained population and Median⁻ represents the median FITC intensity of unstained population. SD⁻ represents the standard deviation in FITC intensity distribution of control population that was not processed through a microfluidic device, denoted as ND (no device) in the rest of this article. Significance testing was performed using R and JASP. Firstly, the homogeneity of variance was tested using Levene's test followed by ANOVA to compare delivery and viability across channel designs. Finally, post hoc Bonferroni and Games–Howell tests were performed for pairwise comparisons for equal and unequal number of samples respectively.

High speed video tracking and analysis

Cells suspended in opti-MEM without any cargo molecules were flown through the channels. High-speed camera (Phantom Vision Research, NJ, USA) imaging was performed using an inverted bright-field microscope (Eclipse Ti, Nikon, Japan). Videos of cells were taken with a 500-pixel by 500-pixel region at ultra-high speed (26 000 fps). Specifically, we recorded videos at the 1st, 3rd and 6th constrictions in each design variant for up to 4 seconds per video. For trajectory analysis, frames were analysed using TrackMate plugin (https://rsb.info.nih.gov/ij/) in FIJI after subtracting a median background. The identified cell trajectories were exported to CSVs. A custom script was developed using R and MATLAB to quantify the cell kinematics from image analysis.

Derivation of mechanistic descriptors and multivariate analysis setup

The modelled channel heights were subsequently used to calculate the average strain imposed under the ridge. For each trajectory, transition times between pre-ridge and under-ridge were used to calculate instantaneous strain rates ([small epsi, Greek, dot above]) and force descriptors. To compute BBH forces, a pre-ridge region was defined for each trajectory as a window spanning one cell diameter upstream of ridge entry. This region was defined to capture the dominant acceleration phase experienced by the cell as it approached geometric confinement. Using the known channel and ridge heights, the local fluid velocity at the start and end of this pre-ridge region was determined from CFD, and a linear acceleration profile was interpolated between these points. Tracked cell velocities and the interpolated fluid velocity profile were then used to compute slip velocities and particle Reynolds numbers along each trajectory. The time integral of the slip velocity derivative was subsequently evaluated to compute the BBH force according to the classical formulation^28,35
Various subsets of resulting variable set were analysed by a custom R script to identify PCA models with adequate sampling and cluster separation. Evaluation criteria included Kaiser–Meyer–Olkin (KMO)³⁶ measure of sampling adequacy (>0.6) and Davies–Bouldin index (DB index).³⁷ The final model was selected as the subset that also showed significant separation by MANOVA (Pillai's trace, p < 0.05).

Results

CFD simulations and flow experiments enable estimation of cell strain in PDMS devices

To determine the extent of PDMS deformation and resulting strain experienced by cells in the microchannels, we first established baseline flow resistance from rigid channel models (Fig. 2(a)). Experimental pressure values (ΔP) were converted from mbar to Pascals and plotted against flow rate (Q) to account for deformation of PDMS devices across flow rates up to 1000 μL min⁻¹. The experimental ΔP–Q relationships deviated from rigid channel predictions and were fit using the Gervais model.³⁴ From these fits, deformed channel heights (H_c,deformed) and deformed constriction heights (H_deformed) were extracted at each flow rate. As predicted by Gervais model, our fits showed a dependence not only the flow rate (Fig. 2(b) and S1) but also along the length of the channel (Fig. 2(c)). For any chosen architecture, deformation increased with increase in flow rates (Table 2). In all channels, the first ridge experienced the highest deformation due to the highest-pressure difference. In C01, constriction height was determined to be 12.77 μm at the 1st ridge under 250 μL min⁻¹, while in C05 the corresponding deformed height was 10.26 μm, and 7.69 μm in C10 (Fig. 2(c)). Notably, the constriction height profile along the channel for C05 at 100 μL min⁻¹ was broadly similar to C10 at 250 μL min⁻¹ (Table 2, Fig. S1) across all ridges. These estimated channel heights were subsequently used to calculate the average cell strain imposed under the ridge regions. This was done by finding the average constriction height across all the ridges in a microchannel and dividing by the average cell diameter. For the 10-channel design (C10), the fitted models predicted an average constriction height of 6.6 μm for B16F10 cells and 6.2 μm for T-cells, based on averages over all 6 constrictions. Corresponding average compressive strains in C10 were computed to be 0.58 for B16F10 and 0.44 for T-cells.

Fig. 2 Constriction height estimation in PDMS devices (a) shows the simulated and measured flow rates in the channel under constant pressures (simulations assuming rigid PDMS device of gap height 4.5 μm). The channels C01, C05 and C10 were designed such that the theoretical flow resistance remains same amongst the three. Measured P–Q relationship was fitted with Gervais model to identify deformed channel heights as shown in (b), (c). (b) Shows the derived channel heights for C05 for a fabricated constriction height of 4.5 μm. Most optimization experiments were performed at flow rates ≤500 μL min⁻¹. (c) Shows heights for different channel designs with fabricated constriction heights of 4.5 μm.

Table 2 Summary of deformed channel dimensions as computed using the modified Gervais model for a channel with rigid constriction height 4.5 μm

	Ridge number	Flow rate
		100 μL min^−1	250 μL min^−1	350 μL min^−1
C01	1	9.55	12.77	14.38
	2	9.03	11.96	13.41
	3	8.46	11.07	12.34
	4	7.80	10.03	11.11
	5	7.02	8.80	9.66
	6	5.98	7.15	7.70
C05	1	7.77	10.26	11.58
	2	7.42	9.65	10.82
	3	7.01	8.97	9.97
	4	6.55	8.19	9.03
	5	6.01	7.30	7.95
	6	5.35	6.17	6.58
C10	1	6.21	7.69	8.44
	2	5.99	7.31	7.99
	3	5.74	6.89	7.49
	4	5.47	6.42	6.92
	5	5.19	5.93	6.32
	6	4.84	5.25	5.48

---

Narrow parallelized channels increase delivery efficiency in both adherent and suspension cell lines

We tested the three channel width variants by evaluating uptake of 250 kDa FITC-dextran using flow cytometry (Fig. 3(a)). The delivery efficiency was defined as the percentage of FITC-positive cells, and viability was quantified by 7-AAD staining immediately after microfluidic processing. All experiments discussed in this study were conducted at a fixed total flow rate of 250 μL min⁻¹ unless stated otherwise.

Fig. 3 Improved delivery efficacy across suspension and adherent cell lines (a) representative FITC distributions of B16F10 cells processed through the three design variants. Narrower channels yielded higher number of cells with FITC+ signal (b) and (c) show delivery efficiencies in green (% cells staining positive for FITC) and viabilities in red (% cells staining negatively for 7-AAD) of recovered cell populations after microfluidic processing for T-cells and B16F10 cells respectively. Bar plots represent means across replicates and error bars represent standard error margin (SEM) across replicates. Bonferroni adjusted p-values are denoted on the charts for statistical significance (ns for p ≥ 0.05, * for p < 0.05, ** for p < 0.01, *** for p < 0.001).

Narrower channel widths resulted in a two-fold increase in delivery efficiency for both adherent (B16F10) and suspension (T-cells) cell populations. Specifically, delivery efficiency increased from 23.6 ± 8.9% to 69.8 ± 3.6% in T-cells (Fig. 3(b)) and from 20.3 ± 4.3% to 55.3 ± 3.0% in B16F10 cells (Fig. 3(c)). Compared to C01, both C05 and C10 devices significantly increased delivery efficiency in both cell types (p < 0.01). While delivery efficiency between C05 and C10 was not significantly different for B16F10 cells (p = 9.70 × 10⁻¹), a significant difference was observed between C05 and C10 in T-cells (p = 3.6× 10⁻³). These delivery gains were accompanied by a decrease in viability. Mean viabilities dropped from 89.9 ± 4.2% for C01 to 58.8 ± 15.4% for C10 in case of T-cells, and from 79.7 ± 6.8% for C01 to 66.0 ± 20.0.% for C10 in case of B16F10 (Fig. 3(b) and (c)). For B16F10 cells, C05 and C10 design did not produce statistically significant difference in viability however, for T-cells, Bonferroni-adjusted p-value revealed a significant decrease (p = 1.37 × 10⁻²). Importantly, despite the higher fraction of FITC+ cells, the staining index remained comparable across all designs (Fig. S2), suggesting that the amount of cargo per cell did not differ substantially.

Narrow multi-channel designs promote uniform compression dynamics and reduced occlusion

Increasing the number of parallel channels while decreasing channel width also improved reproducibility of delivery efficiency across multiple device replicates, as reflected by a lower coefficient of variation (CV) in delivery efficiency as compared to C01 (Fig. 4(a)). Comparing delivery efficiency, for T-cells (N = 3), CV decreased from 0.65 in C01 to 0.11 in C10. For B16F10 cells (N = 4), CV decreased from 0.37 in C01 to 0.26 in C10. While the variability in delivery efficiency decreased with narrower channels, the variability in viability increased. We hypothesized that this suggested higher populations of cells experiencing rapid compressions and variability in viability stemming from heterogeneity in cell mechanics.

Fig. 4 Uniform compression dynamics of cells in narrow-width channels promote scalability for mechanoporation (a) the inter-trial variability goes down with multi-channel designs for both T-cells (N = 3) as well as B16F10 (N = 4) as shown by coefficient of variation plot. (b) Shows the variation of delivery in B16F10 cells across 5 outlets of C05 with a unique single inlet 5-outlet design where H = 4.5 μm. (c) Images of cell trajectories in channels of varying widths and respective FITC distributions of B16F10 cells (d) computed probability density functions of positions at which the cells hit the ridges in 3 different channel designs across the length of the channel. The density function becomes narrower moving from C01 to C10. Moreover, we also see a shift in trajectories towards the center of the channel as we move towards the 6th constriction. (e) Shows a schematic describing the image analysis method to identify ridge areas and occlusion areas used for quantification of occlusion burden. (f) Shows normalized occlusion areas with respect to ridge areas under two different flow rates. Specially for the lower flow rate of 100 μL min⁻¹, C05 has significantly less occlusion than C01. For both conditions C10 has significantly occlusion than C01.

To evaluate scalability, delivery efficiency and viability of B16F10 cells were measured across the five outlets of C05 with a fabricated constriction height of 4.5 μm (Fig. 4(b)). Both delivery and viability remained consistent across outlets, with no statistically significant outlet-to-outlet variation (ANOVA p = 0.89 for delivery, p = 0.19 for viability). To further validate scalability, a 20-channel design (W = 50 μm) was fabricated and tested alongside the 10-channel design. Delivery efficiency and viability were comparable to those obtained in C10 operated at same flow rate per channel (Fig. S3), confirming that performance is preserved with increasing channel number.

To assess the uniformity of compression, we captured high speed videos of cells traversing the channel with a high-speed camera. Imaging revealed distinct trajectory patterns across channel widths (Fig. 4(c)). We first used the trajectories to find out the intersection point of ridge and cell trajectory. The length of the ridge (R_L) was used to normalize the distance of intersection points from channel centre (defined as x = 0). Finally, we visualized computed probability density functions (PDFs) (Fig. 4(d)) of the normalized positions at which cells met the ridges. Compression pathways in C01 were distributed broadly across the channel width, whereas C05 and C10 exhibited narrowly distributed trajectories. Moreover, as cells progressed towards the sixth constriction, the PDFs shifted toward the channel centreline, with the effect being most pronounced in C10. The distribution of trajectories confirmed our hypothesis that reduced variability in delivery efficiency stemmed from more uniform compression dynamics of the entire cell population processed using the microfluidic device.

We also computed the burden of channel occlusion by normalizing visible 2D occluded areas post-processing. Images were taken for all ridges and then normalized with the 2D ridge area for each channel. We found that occlusions were reduced in parallelized designs (Fig. 4(e)). At any given flow rate, normalized clog fraction was significantly greater in C01 than in C10 (p <0.01). For all channels, higher flow rate of 350 μL min⁻¹ further decreased the occlusion area fraction in comparison with 100 μL min⁻¹.

Narrow channels consistently outperform wide channel designs across conditions with persistent delivery-viability trade-offs

To further test whether the trends in delivery and viability metrics of the three designs are consistent across fluid dynamics and biomechanical parameters, we tested 3 different flow rates, 3 different surface chemistries, and pharmacologically softened cells (Fig. 5). A summary of deformed ridge heights as computed from the fitted model discussed previously is provided in Table 2.

Fig. 5 Observing delivery efficiency and viability trends with multiple flow rates, surface chemistries, and cell stiffness: (a) and (b) comparing delivery efficacies at different flow rates, we found that delivery efficacies improved with higher flow rates for all 3 channel widths. This confirmed that mechanoporation efficacies are governed more by fluid forces that increased by varying flow rates from 100 μL min⁻¹ to 350 μL min⁻¹. (c) and (d) Highlight the effect of adhesion forces. Higher adhesion forces lead to higher delivery efficiency in melanoma cells, but the increase is minimal compared to channel width modulation. (e) Shows the changes in cell diameter after treating B16F10 melanoma cells with 10 μM cytoD for 60 minutes. Black line represents the mean of multiple trials. (f) Shows non-significant changes in cell viability before cells were processed through the microfluidic device. Black bar represents the mean of multiple trials. (g) No significant changes in delivery performance were observed with C01 but a significant drop in delivery efficiency was observed with C05. (h) No statistically significant change in viability was observed with cytoD induced cell softening while device geometry had significant effects. All bars represent means with the SD plotted as error bars. * p < 0.05, ** p <0.01, *** p < 0.001.

Although increasing flow rate enlarged the estimated deformed gap heights, delivery increased and viability decreased across designs. This counter-intuitive observation indicates that effective confinement alone is insufficient to explain the observed trends (Fig. 5(a) and (b); N = 3). For C01 there was a non-significant difference between delivery efficiencies increasing from 9.3 ± 3.5% at 100 μL min⁻¹ to 20.7 ± 3.2% at 250 μL min⁻¹ (p = 0.27) and plateaued at 19.3 ± 3.6% for 350 μL min⁻¹ (p = 0.43). In case of C05, delivery efficiencies increased non significantly from 29.2 ± 4.5% at 100 μL min⁻¹ to 42.0 ± 8.9% at 250 μL min⁻¹ (p = 0.16) however, significantly increased to 55.3 ± 3.0% for 350 μL min⁻¹ (p = 4.21 × 10⁻⁴). For C10, across all flow rates, delivery efficiencies remained significantly higher than C01 (p < 0.01). It also increased in a statistically significant manner from 29.6 ± 4.1% at 100 μL min⁻¹ to 55.6 ± 7.3% at 350 μL min⁻¹. Viabilities also decreased monotonously with higher flows, reaching as low as 30.4 ± 6.7% in B16F10 cells at flow rate of 350 μL min⁻¹. The delivery performance retained the relative order of designs (C10 ≥ C05 > C01) which was the inverse of viability trends (C01 > C05 ≥ C10). Since C10 and C05 did not produce significant differences in terms of delivery efficiency at 250 μL min⁻¹, C10 was not tested for the following studies to compare similar pre-ridge forces.

To study whether surface interactions affected delivery efficiency and viability in B16F10 cells, we treated devices reducing adhesion with Pluronic passivation and increasing adhesion with APTES molecules (Fig. 5(c) and (d)). Decreasing adhesion decreased delivery efficiency from 19.6% to 5.4% in C01 (p < 0.01) and from 46.6% to 37.8% in C05 (p < 0.05). Viability increased under Pluronic treatments for both C01 and C05. Increasing adhesion with APTES passivation did not significantly alter delivery efficiency or viability compared to standard control surface treatment (p > 0.05). For both surface chemistries, delivery remained significantly higher in C05 compared to C01 (p < 0.01).

In order to investigate the effect of cell properties on the delivery-viability trade-off, we treated cells with Cytochalasin D (cytoD), known to reduce cell stiffness by disrupting actin polymerization.²² No significant change in cell size or viability was observed before mechanoporation (Fig. 5(e) and (f)). Delivery efficiency reduced in B16F10 cells across all designs with 10 μM cytoD treatment (Fig. 5(g)). In C01, delivery decreased from 27.7 ± 5.9% in untreated cells to 17.8 ± 5.9% with cytoD. In C05, delivery decreased from 57.2% to 35.0% (p = 1.11 × 10⁻³). We also observed that softer cytoD treated cells had more viability when processed through C05 but no significant changes when processed through C01 (Fig. 5(h)). When a smaller dose of 2 μM was tested for C05 design, no significant change in delivery performance as observed (data not shown) as expected from the robustness of B16F10 cells to cytochalasin-D treatments.³⁸ Importantly, the width-dependent trend of higher delivery in C05 relative to C01 was preserved for both untreated and treated cells amongst all cell lines tested.

High-speed tracking analysis reveals sharp velocity and acceleration transitions near constrictions

To quantify cell kinematics across constriction regions, high-speed tracking was performed, and trajectories were analyzed to obtain cell positions (x,y coordinates) and frame numbers. Instantaneous velocity and acceleration along the channel were computed with trajectories having more than three points in the pre-ridge, under-ridge, or post-ridge regions (Fig. 6(a)). Across all designs, cells exhibited distinct velocity and acceleration profiles relative to ridge location. In the pre-compression region (I) and post-compression region (V), normalized velocities remained similar for all three designs while cell acceleration remained negligible (Fig. S5(c)). As cells approached the constriction, velocities peaked sharply, with C10 showing the highest average increase compared to C01 and C05 (Fig. 6(b)). Acceleration magnitude profiles (Fig. 6(c)) also revealed transient spikes near constrictions with C10 accelerating as well as decelerating faster than the other two channel designs.

Fig. 6 Computing cell velocities and acceleration in microfluidic channel (a) a schematic highlighting different regions of the channel for high-speed tracking analysis and relative position of the particle for which simulations were used to calculate forces. (b) and (c) Shows the tracked relative velocities of cells in the channel and acceleration of the particle respectively. The velocities in the shown regions were computed relative to the initial cell velocity in channel far away from the ridge. The graphs portray higher velocities and accelerations experienced by cells in multi-channel devices.

Principal component analysis reveals that dynamic loading kinetics, rather than static strain magnitude, distinguish high-delivery narrow channel designs

To gain insight into the mechanistic differences in cellular loading across microfluidic channel designs, multivariate analysis was performed using experimentally derived mechanoporation metrics and hydrodynamic force parameters. Because both strain rate and BBH forces increase with smaller constriction heights, these kinetic metrics are expected to co-vary across designs and flow rates. We therefore interpret BBH not as an independent geometric variable, but as a history-dependent descriptor of unsteady hydrodynamic loading that complements strain rate in capturing the ridge-entry loading regime.

Together with mechanoporation metrics delivery, viability and staining index, the computed hydrodynamic parameters formed the input dataset for multivariate analysis. A summarized list of all variables and computation is provided in Table S2. The sample dataset comprised of B16F10 mechanoporation across flow rates and device designs (N = 68). The first two principal components explained 45.4% and 23.8% of total variance, respectively, capturing 69.2% cumulative variance (Fig. 7(a)). A summary of variable contributions to the principal components is provided in Fig. S6. PC1 was dominated by dynamic loading-related variables including strain rate, BBH forces, and Stokes drag, along with delivery efficiency and viability drop. In contrast, PC2 was more strongly associated with deformation magnitude including strain fraction and deformed channel height. This separation indicates that dynamic loading descriptors and static deformation metrics capture distinct sources of variation within the dataset. Projection of samples into PC1–PC2 space revealed that the device classes (C01, C05, and C10) were separated primarily along PC1, while showing broader spread along PC2 (Fig. 7(a)). Cluster quality analysis yielded a silhouette score of 0.30, indicating moderate but biologically realistic separation across designs. Visualizing strain rate against delivery efficiency and viability revealed a clear correlation of strain rate and BBH forces on determining delivery efficiency (Fig. 7(b)).

Fig. 7 Multivariate analysis of fluid dynamics and forces in channels (a) shows PCA analysis of ∼68 samples with varying flow rates and device designs. PC1 on x-axis explained 45.4% of variance in data and PC2 explained 23.8% variance. The 3 different channel designs – C01, C05 and C10 formed distinct clusters. (b) Shows the plots of delivery efficiency and viability plotted against the strain rate. For delivery efficiency (FITC+ cells %), we observed distinct trends based on BBH forces as well as strain rates. For viability, there wasn't any clear pattern with the two biomechanical parameters. (c) Shows the computed BBH forces to capture unsteady viscous effects from high-speed videos of all 3 channel designs at a flow rate of 350 μL min⁻¹. The forces were calculated using estimated slip velocity and a discrete quadrature scheme with damping for high particle Reynold's numbers (d) shows the simulated force component values in pre-ridge region (II) at positions highlighted previously in Fig. 5. Simulations revealed an interesting difference in shear force components experienced by a rigid cell in the middle position of channel. Not only does the total magnitude of these forces increase from C01 to C10, but the parallel component of shear force along the ridge increases more than 2.5× resulting in a more biaxial shear stretch on the cells in narrow width devices. (e) The schematics represent the particle, fluid streamlines and the directions of force components in different channel designs. Force components are not drawn to scale and are representative only.

Correlation analysis links strain rates and pre-ridge BBH forces to delivery and viability outcomes

Correlation analysis (Fig. S6(c)) confirmed strong positive associations between strain rate and delivery efficiency (r > 0.7) as well as BBH force and delivery efficiency (r > 0.6). Viability correlated negatively with delivery (|r| > 0.6) while remaining weakly related to other mechanistic variables (|r| < 0.5). The staining index (ΔSI), reflecting the relative amount of cargo delivered per cell, showed weak or no correlation with strain, strain rate, or force-based parameters (|r| < 0.3). Similarly, the number of cargo molecules per cell in buffer (defined as the concentration of cargo particles divided by concentration of cells) was not significantly associated with either delivery efficiency or viability, consistent with our hypothesis that efficiency increases reflected a greater fraction of permeabilized cells rather than higher uptake per cell. Because BBH forces and strain rate arise from the same channel deformation-driven confinement regime, we assessed their joint predictive value using regression models on a proportional (log–log) scale. Modelling delivery as a function of log([small epsi, Greek, dot above]) alone explained a substantial fraction of variance in delivery efficiency (R² = 0.44), and inclusion of log(F_BBH) significantly improved model performance (R² = 0.49; BIC decreased from 645 to 642.1). In contrast, adding deformed constriction height (H_deformed) to the model did not further improve fit (R² = 0.49), indicating that effective gap size is largely accounted for by mechanistic loading descriptors derived from trajectory-resolved kinematics. Together, these results indicate that deformation-driven changes in channel gap heights are effectively encoded by strain rate and BBH forces, and that unsteady loading captured by BBH contributes beyond proportional strain rate effects.

CFD-derived stress components indicate asymmetric biaxial pre-ridge loading in narrow-width devices

Visualization of computed BBH forces along normalized cell trajectories showed a sharp increase in pre-ridge regions for narrow channels, consistent with higher acceleration fields (Fig. 7(d)). Since strain rate and BBH forces emerged as the strongest predictors of delivery and viability in the PCA and multivariate analyses, we next examined whether channel width alone can also reshape the pre-ridge hydrodynamic force landscape independent of PDMS deformation. To quantify the spatial distribution of forces experienced by cells, we performed rigid-channel COMSOL simulations by placing rigid spherical probes at center, middle, and side positions within the channel cross-section (Fig. 6(a)). Laminar flow analysis revealed that geometric confinement in narrow channels produced substantially larger force components especially in the parallel direction to the ridge. Decomposition of these vectors showed that the parallel shear component increased by more than 2.5-fold in C10 relative to C01, indicating a distinct shear-biased and biaxial loading environment. (Fig. 7(d) and (e)). Although the total peak shear magnitude increased modestly (≈10%) from C01 to C10 (365.3 nN to 395.9 nN; Δ ≈ 30 nN), this absolute change is within a biologically meaningful range for membrane-level mechanotransduction and may act as a confounding contributor alongside strain-rate and BBH-dominated loading kinetics. Importantly, this force shift suggests a mechanistic basis for how channel width can alter pre-ridge hydrodynamic loading even when constriction height is held fixed in rigid models.

Discussion

Across experiments, narrow, parallelized channels (C05, C10) produced higher delivery with reduced occlusion and lower inter-trial variability compared to wide single-channel devices (C01). Because channel width, PDMS deformation, effective gap height, and kinematic loading descriptors co-varied in this study, the multichannel advantage should not be attributed to any single factor in isolation. Conditions with broadly similar estimated deformed gap-height profiles, such as C05 at 100 μL min⁻¹ and C10 at 250 μL min⁻¹, still exhibited different delivery outcomes, supporting the claim that consistent static confinement alone does not explain the performance difference. The PCA structure was particularly informative because device classes separated primarily along PC1, the axis most strongly associated with strain rate and BBH-related loading, while remaining broadly distributed along PC2. This pattern suggests that the main architecture-dependent shift was not simply toward greater static constriction, but toward a distinct dynamic loading regime. While BBH forces and strain rate correlated strongly with delivery and viability across all geometries, absolute strain alone could not account for the observed trends. Staining index showed only weak correlation indicating that performance improvements arose from an increased fraction of permeabilized cells rather than enhanced uptake per cell.

High-speed imaging and CFD revealed that loading kinetics, specifically, pre-ridge acceleration fields and cell velocities in the channel varied significantly between the designs. Narrow channels and the associated reduction in constriction heights of channels elevated strain rates and generated larger BBH forces as fluid rapidly accelerated upstream of the ridge. Delivery increased with strain rate up to ∼0.5 × 10⁴–1 × 10⁴ s⁻¹, after which viability declined without further gains. This defines a favorable operating regime within the conditions tested here, in which rapid loading promotes permeabilization while excessive rates lead to irreversible damage. Importantly, this range should not be interpreted as a universal optimum. Rather, within the architectures and operating conditions tested here, these quantities serve as transferable mechanistic descriptors because they are computed from effective cell deformation and trajectory-resolved kinematics rather than nominal device geometry alone. The current study links these descriptors to bulk delivery and viability outcomes after averaging across tracked cell populations; future single-cell-resolved measurements will therefore be needed to establish sharper quantitative targets. It is also worth noting that the optimum values may shift due to charged interactions of cargo molecules like nanoparticles with cell membranes resulting in changes to yield strain.

Perturbations to adhesion and cytoskeletal stiffness altered absolute delivery and viability but did not change the parametric dependencies across geometries. None of these shifts displaced narrow channels from their performance advantage. The limited effect of surface-adhesion modulation on delivery, despite measurable effects on viability, suggests that wall adhesion/friction plays a secondary role relative to bulk cell stretching and compression in governing membrane permeabilization under the conditions tested here. Pharmacological softening slightly decreased delivery while increasing occlusion in narrow channels reinforcing that device-imposed loading kinetics dominate over cell-intrinsic wall-friction effects and mechanical perturbations in determining outcome.

Although the multichannel designs maintain proportional scaling between per-channel flow rate and channel width, the particle-to-channel size ratio increases substantially as width decreases. In confined laminar flows, finite-sized particles perturb the surrounding velocity field more strongly particularly when blockage is high, producing sharper spatial velocity gradients and enhanced slip-acceleration transients at ridge entry. These confinement-driven distortions provide a mechanistic basis for the observed elevation in strain rate and BBH-related loading in narrow channel designs. Importantly, all channel Reynolds numbers remained below 40 in our study indicating that these effects arise from finite-size confinement rather than transition to turbulence. Regression analyses revealed that gap-mediated effects are largely encoded by mechanistic loading descriptors derived from trajectory-resolved kinematics, and that BBH captures aspects of the loading history not fully described by proportional strain-rate scaling. The prominence of BBH forces reflects operation in a regime characterized by rapid, sub-millisecond accelerations at ridge entry. Prior studies typically lacked either the temporal resolution or spatial coverage to capture transient spikes in particle velocity to compute BBH forces. Even here, resolving BBH required ridge-by-ridge imaging and segmentation of separate regions of the device. The inability to image the entire device simultaneously highlights the need for future high-speed, full-device imaging to map BBH contributions with higher fidelity.

Rigid-channel CFD clarifies the role of channel width beyond its effect on deformation. Simulations revealed that narrowing channel width reshapes the pre-ridge hydrodynamic force landscape, biasing forces toward the direction parallel to the ridge and producing a more asymmetric biaxial loading environment. The total increase in peak shear magnitude of ∼30 nN falls within a biologically relevant range and may act as a confounding contributor alongside strain-rate and BBH-dominated kinetics. Importantly, these results also provide a physical basis for why BBH serves as a meaningful surrogate for pre-ridge hydrodynamic changes that accompany channel narrowing.

Together, the study suggests that under the ultra-fast deformation conditions studied here (strain rates up to 10⁴ s⁻¹ and Reynolds numbers below 40), non-steady state metrics including pre-ridge BBH forces and under-ridge strain rate were the most predictive parameters for both delivery and viability. These results motivate the field to move beyond quasi-static descriptors and to more systematically incorporate transient, history-dependent loading metrics when comparing mechanoporation platforms. More broadly, this framework provides practical design guidance for microfluidic platforms that rely on repeated extensional or converging–diverging geometries, emphasizing deliberate control of ridge-entry acceleration fields and force partitioning rather than strain magnitude alone to systematically tune delivery efficiency and cell survival across operating regimes and cell types.

Conclusion

These findings provide actionable design guidance by establishing that narrow channels generate stable trajectories and loading kinetics without relying on serpentine flow-focusing. Narrow geometries amplified both total BBH force magnitude and shear-biased components parallel to ridge axes, forming an asymmetric biaxial loading environment. Prior membrane-failure studies²⁶ report that biaxial loading accelerates rupture compared to uniaxial stretching, explaining the viability penalty that accompanies acceleration-dominated loading. These force components represent newly identifiable tuning knobs for future optimization. Early observations suggest that ridge-angle modification may rebalance shear and compression components, reducing biaxial stress without sacrificing acceleration. Systematic exploration of ridge shape, angle, and spacing offers a clear path for next-generation mechanoporation platforms.

Together, these findings show that narrow, parallelized microchannels (Re < 40) enhance mechanoporation by generating high strain-rate and BBH-force environments while simultaneously inducing biaxial stresses that limit viability. The proposed framework provides a mechanistically grounded basis for refining future device designs while motivating validation across additional cell types and cargo classes as a step toward predictive capability.

Author contributions

A. G. conceptualized the study, performed experiments, collected and analysed data, and wrote the manuscript. C. B. and J. Z. performed experiments and analysed data. H. A. and P. S. contributed to study design and data analysis. G. P., A. A., and T. S. acquired funding, supervised the research, and edited the manuscript. All authors reviewed and approved the final manuscript.

Conflicts of interest

TS and AA receive income, hold intellectual property, and are equity holders of CellFE, Inc., a company that is commercializing microfluidic cell engineering technology. The remaining authors declare no competing financial interests.

Data availability

The data supporting this article are provided within the main manuscript and the supplementary information (SI), including processed datasets and parameter tables. All parameters required to reproduce the CFD simulations are described in the methods section, and representative simulation outputs are included in the figures and SI. MATLAB scripts used for kinematic analyses are available in a GitHub repository – https://github.com/ag-git-0/track-to-kinematics.git/. Raw high-speed imaging datasets and full-resolution flow cytometry files are available from the corresponding author upon reasonable request. Requests for materials should be directed to Todd Sulchek (E-mail: todd.sulchek@me.gatech.edu).

Supplementary information is available. See DOI: https://doi.org/10.1039/d6lc00143b.

Acknowledgements

This work was supported by the National Institutes of Health under award numbers R01GM144727 and R01EB030367. We acknowledge the Petit Institute of Bioengineering Cytometry Core for technical support. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

References

1. J. Loo, I. Sicher, A. Goff, O. Kim, N. Clary, A. Alexeev, T. Sulchek, A. Zamarayeva, S. Han and M. Calero-Garcia, Sci. Rep., 2021, 11, 21407.
2. C. Kwon and A. J. Chung, Lab Chip, 2023, 23, 1758–1767.
3. L. R. Bohrer, N. E. Stone, A. T. Wright, S. Han, I. Sicher, T. A. Sulchek, R. F. Mullins and B. A. Tucker, Stem Cells, 2023, 41, 1037–1046.
4. T. Yu, N. Jhita, P. G. Shankles, A. Fedano, N. B. Kramer, S. S. Raikar and T. Sulchek, Lab Chip, 2023, 23(22), 4804–4820.
5. H. Nejadnik, K. O. Jung, A. J. Theruvath, L. Kiru, A. Liu, W. Wu, T. Sulchek, G. Pratx and H. E. Daldrup-Link, Theranostics, 2020, 10, 6024–6034.
6. K. O. Jung, A. J. Theruvath, H. Nejadnik, A. Liu, L. Xing, T. Sulchek, H. E. Daldrup-Link and G. Pratx, Sci. Rep., 2022, 12(1), 2955.
7. A. Adamo and K. F. Jensen, Lab Chip, 2008, 8, 1258–1261.
8. H. G. Dixit, R. Starr, M. L. Dundon, P. I. Pairs, X. Yang, Y. Zhang, D. Nampe, C. B. Ballas, H. Tsutsui, S. J. Forman, C. E. Brown and M. P. Rao, Nano Lett., 2020, 20, 860–867.
9. D. M. Hallow, R. A. Seeger, P. P. Kamaev, G. R. Prado, M. C. LaPlaca and M. R. Prausnitz, Biotechnol. Bioeng., 2008, 99, 846–854.
10. M. E. Kizer, Y. Deng, G. Kang, P. E. Mikael, X. Wang and A. J. Chung, Lab Chip, 2019, 19, 1747–1754.
11. K. O. Jung, T. J. Kim, J. H. Yu, S. Rhee, W. Zhao, B. Ha, K. Red-Horse, S. S. Gambhir and G. Pratx, Nat. Biomed. Eng., 2020, 4, 835–844.
12. A. Liu, M. Islam, N. Stone, V. Varadarajan, J. Jeong, S. Bowie, P. Qiu, E. K. Waller, A. Alexeev and T. Sulchek, Mater. Today, 2018, 21, 703–712.
13. Y. Deng, S. P. Davis, F. Yang, K. S. Paulsen, M. Kumar, R. Sinnott DeVaux, X. Wang, D. S. Conklin, A. Oberai, J. I. Herschkowitz and A. J. Chung, Small, 2017, 13, 1700705.
14. D. Sevenler and M. Toner, Nat. Commun., 2024, 15, 115.
15. J. A. Jarrell, A. A. Twite, K. H. W. J. Lau, M. N. Kashani, A. A. Lievano, J. Acevedo, C. Priest, J. Nieva, D. Gottlieb and R. S. Pawell, Sci. Rep., 2019, 9, 3214–3214.
16. G. Kang, D. W. Carlson, T. H. Kang, S. K. Lee, S. J. Haward, I. Choi, A. Q. Shen and A. J. Chung, ACS Nano, 2020, 14, 3048–3058.
17. A. Sharei, J. Zoldan, A. Adamo, W. Y. Sim, N. Cho, E. Jackson, S. Mao, S. Schneider, M.-J. Han, A. Lytton-Jean, P. A. Basto, S. Jhunjhunwala, J. Lee, D. A. Heller, J. W. Kang, G. C. Hartoularos, K.-S. Kim, D. G. Anderson, R. Langer and K. F. Jensen, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 2082–2087.
18. E. Bar-Kochba, M. T. Scimone, J. B. Estrada and C. Franck, Sci. Rep., 2016, 6, 30550.
19. M. Urbanska, H. E. Muñoz, J. Shaw Bagnall, O. Otto, S. R. Manalis, D. Di Carlo and J. Guck, Nat. Methods, 2020, 17, 587–593.
20. Y. Qiang, J. Liu, M. Dao, S. Suresh and E. Du, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 19828–19834.
21. T. Shigematsu, K. Koshiyama and S. Wada, Sci. Rep., 2015, 5, 15369–15369.
22. J. R. Lange, J. Steinwachs, T. Kolb, L. Lautscham, I. Harder, G. Whyte and B. Fabry, Biophys. J., 2015, 109, 26–34.
23. M. Razizadeh, M. Nikfar, R. Paul and Y. Liu, Biophys. J., 2020, 119, 471–482.
24. J. Rahimzadeh, F. Meng, F. Sachs, J. Wang, D. Verma and S. Z. Hua, Am. J. Physiol., 2011, 301, C646–C652.
25. A. M. Esfahani, J. Rosenbohm, B. T. Safa, N. V. Lavrik, G. Minnick, Q. Zhou, F. Kong, X. Jin, E. Kim, Y. Liu, Y. Lu, J. Y. Lim, J. K. Wahl, M. Dao, C. Huang and R. Yang, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2019347118.
26. M. A. Murphy, S. Mun, M. F. Horstemeyer, M. I. Baskes, A. Bakhtiary, M. C. LaPlaca, S. R. Gwaltney, L. N. Williams and R. Prabhu, J. Biomol. Struct. Dyn., 2019, 37, 1346–1359.
27. M. I. Maremonti, V. Panzetta, P. A. Netti and F. Causa, J. Nanobiotechnol., 2024, 22, 441.
28. D. Jaganathan, S. G. Prasath, R. Govindarajan and V. Vasan, Front. Phys., 2023, 11, 1167338.
29. D. E. Peterson and B. R. Adams, Chem. Eng. Res. Des., 2025, 221, 11–21.
30. A. J. Dorgan and E. Loth, Int. J. Multiphase Flow, 2007, 33, 833–848.
31. A. Liu, T. Yu, K. Young, N. Stone, S. Hanasoge, T. J. Kirby, V. Varadarajan, N. Colonna, J. Liu, A. Raj, J. Lammerding, A. Alexeev and T. Sulchek, Small, 2020, 16, 1903857.
32. E. S. Asmolov, A. L. Dubov, T. V. Nizkaya, J. Harting and O. I. Vinogradova, J. Fluid Mech., 2018, 840, 613–630.
33. K. M. Young, P. G. Shankles, T. Chen, K. Ahkee, S. Bules and T. Sulchek, Biomicrofluidics, 2022, 16, 034104.
34. T. Gervais, J. El-Ali, A. Günther and K. F. Jensen, Lab Chip, 2006, 6, 500.
35. A. B. Basset, A treatise on hydrodynamics: with numerous examples, Cambridge, 1888.
36. H. F. Kaiser and J. Rice, Educ. Psychol. Meas., 1974, 34, 111–117.
37. D. L. Davies and D. W. Bouldin, IEEE Trans. Pattern Anal. Mach. Intell., 1979, 1, 224–227.
38. Y. Ujihara, D. Ono, K. Nishitsuji, M. Ito, S. Sugita and M. Nakamura, Cell. Mol. Bioeng., 2021, 14, 309–320.

---