A computational study of cell membrane damage and intracellular delivery in a cross-slot microchannel

Abstract

We propose a three-dimensional computational framework to simulate the flow-induced cell membrane damage and the resulting enhanced intracellular mass transport in a cross-slot microchannel. We model the cell as a liquid droplet enclosed by a viscoelastic membrane and solve the cell deformation using a well-tested immersed-boundary lattice-Boltzmann method. The cell membrane damage, which is directly related to the membrane permeability, is considered using continuum damage mechanics. The transport of the diffusive solute into the cell is solved by a lattice-Boltzmann model. After validating the computational framework against several benchmark cases, we consider a cell flowing through a cross-slot microchannel, focusing on the effects of the flow strength, channel fluid viscosity and cell membrane viscosity on the membrane damage and enhanced intracellular transport. Interestingly, we find that under a comparable pressure drop across the device, for cells with low membrane viscosity, the inertial flow regime, which can be achieved by driving a low-viscosity liquid at a high speed, often leads to much larger membrane damage, compared with the high-viscosity low-speed viscous flow regime. However, the enhancement can be significantly reduced or even reversed by an increase of the cell membrane viscosity, which limits cell deformation, particularly in the inertial flow regime. Our computational framework and simulation results may guide the design and optimisation of microfluidic devices, which use cross-slot geometry to disrupt cell membranes to enhance intracellular delivery of solutes.

---

1. Introduction

Intracellular delivery, which refers to the process of introducing exogenous cargoes such as small-molecular drugs, nucleic acids or synthetic nanoparticles into living cells, has drawn increasing attention due to its ground-breaking applications ranging from mRNA vaccines, preparation of CAR-T cells for cancer treatment, to genetic engineering of plants.^1–5 To enhance the internalisation of foreign cargoes in cells, various approaches have been proposed, which are commonly classified into carrier-mediated and membrane disruption-based methods.⁴ The former uses carriers, in various forms including viral vectors, liposomes or cell-penetrating peptides to deliver encapsulated cargoes into cells using their pathways such as infection or fusion.^1,6,7 Membrane disruption-based methods often apply external energy, e.g., in an electrical, thermal or mechanical form, to cells to create structural damage including pores and ruptures in the cellular membrane, enabling foreign cargoes dispersed in a solution to enter the cells.^8–10 While carrier-mediated delivery strategies are limited in the types of cargo–carrier combinations that can be used in specific applications, membrane disruption-based methods are nearly universal, allowing for the delivery of almost any cargo that can be dissolved in a solution.¹

Amongst the numerous membrane disruption-mediated intracellular delivery methods, microfluidic-based membrane disruption, such as cell squeezing through a narrow channel^11–13 or cell stretching in a micro cross-slot,^10,14 have been particularly promising, for their simplicity and superior capabilities in precisely controlling the flow condition and fluid stresses, high processing throughput rate, and potential for fully automated systems.^1,4,15 In those systems, cells are flown through and deformed in well-defined flow geometries, where the cell membrane is disrupted largely via two mechanisms: mechanosensitive (MS) channel opening and membrane pore or rupture formation (i.e., membrane damage).

MS channels that form in the lipid bilayer of the cell membrane are mainly gated by membrane tension, introduced by cell deformation, and therefore their opening and closing can be controlled by the ambient flow strength.¹⁶ Exogenous cargoes can be transported through the membrane during the opening state of MS channels.¹¹ When cells undergo significant deformation, the cell membrane can be damaged through the formation of nano/micro pores or ruptures,^17–21 leading to enhanced membrane permeability to exogenous substances. Unlike MS channels, whose opening and closing are almost instantaneously determined by the transient membrane tension, membrane damage could last for seconds or minutes before the membrane reseals.^11,22

Despite recent developments and the success of microfluidic-based techniques for membrane disruption-mediated intracellular delivery, numerical simulation of the problem at the single cell level, taking into account the flow-induced cell deformation, membrane damage and the associated mechanical weakening, as well as the solute transport inside and around the cell with spatial resolution, is still at a very early stage. Amongst the relevant pioneer studies, Liu and co-workers have built a comprehensive three-dimensional (3D) model which considers the cell deformation using an immersed boundary method.^23–27 The cell membrane damage depends on the local strain and is predicted using molecular dynamics. The membrane damage leads to enhanced membrane porosity, and the cross-membrane mass transport is resolved using a 1D mathematical model. Specifically, the comprehensive model takes into account the recovery of the membrane damage, on a time scale of several minutes. With the model, the group considered the intracellular drug delivery by rapid squeezing in a constricted channel,²⁶ and various other problems such as haemoglobin release from red blood cells.^24,26 Luo and Bai²⁸ considered the release of a diffusive solute from an elastic microcapsule flowing through a constricted channel. The mass transport was resolved by the convection–diffusion equation, and the diffusion coefficient at the membrane was assumed to depend on the membrane tension through an exponential equation. Misbah and co-workers pioneered the coupling of cell deformation and membrane MS channel opening for numerical simulation of cross-membrane solute transport.^29–32 In their model, the membrane permeability depends on the local membrane stress and the temporal evolution of the membrane curvature. They employed a lattice Boltzmann model for the mass transport and developed a highly efficient approach to deal with the three different types of boundary conditions at the deforming cell membrane. The method has been used to study ATP release from RBCs flowing in microchannels^29–31 and capillary networks.³²

In the studies mentioned above, the effect of the cell membrane viscosity, which has been shown to play important roles in the transient cell deformation,^33–38 has not been considered. Besides, the models often neglect the membrane mechanical degradation (i.e., reduction of strain energy density) as a result of the membrane damage. In the present study, we develop a comprehensive computational framework to take into account those important factors. We consider the cell using a recent mechanical model that has been verified against flow experiments conducted on human leukaemia cells in a constricted microchannel. The model takes into account both the cell membrane elasticity and viscosity. The fluid flow and cell deformation are simulated using a well-tested 3D immersed-boundary lattice-Boltzmann method (LBM). We employ continuum damage mechanics (CDM) to model the evolution of cell membrane damage, caused by flow-induced deformation, and its effect on membrane mechanical degradation. The CDM approach has been widely used for simulating the damage process of biological tissues and bioartificial microcapsules.^39–42 The cell membrane permeability depends on the membrane damage, and the transport of the diffusive solute across the deforming cell membrane, inside and outside the cell, is solved using a lattice-Boltzmann model.

We consider the flow-induced cell deformation, membrane damage and associated transport of a diffusive solute into a cell flowing through a cross-slot microchannel. The channel geometry has been successfully applied to disrupt the cell membrane to enhance the intracellular delivery of agents including dextran, DNA and mRNA with high throughput (up to 10⁶ cells per min),^10,14 as it can generate a strong extensional flow with a stagnation point at the centre of the cross-slot region. At high flow speeds when the inertial effect becomes important, the extensional flow transits to a spiral flow.^43,44 Despite the great success in practical applications, there are still a number of fundamental open questions. In the context of generating membrane damage to enhance intracellular delivery, experiments can be conducted in either the inertial or viscous flow regimes, with comparable pressure drop across the whole device. The inertial flow can be achieved by driving a low-viscosity suspension through the channel at a high speed, while the viscous flow regime features a high-viscosity medium flowing at a low speed. Which regime is more effective in mechanically disrupting the cell membrane and enhancing cross-membrane mass transport? What are the effects of cell membrane properties, in particular the membrane viscosity, in this process? The present study aims to address those open questions.

This paper is organised as follows: the flow geometry, the governing equations and the main dimensionless parameters are detailed in Section 2; the numerical methods and their validations are presented in Section 3. In Section 4, we present the findings of our study on cell deformation, membrane damage, and intracellular mass transport of a cell flowing through a cross-slot microchannel under various flow regimes. We then conclude the paper in Section 5.

2. Problem statement

We consider the cell deformation, membrane damage, and intracellular mass transport as a cell flows through a cross-slot microchannel, which is illustrated in Fig. 1. The four branches of the inflow and outflow channels all have a constant square cross-section with a side length of 2l. A three-dimensional Cartesian coordinate is defined with the x-axis along the inflow channel branches, the z-axis along the outflow branches, and x = y = z = 0 at the centre of the channel cross-slot. The fluid flow is governed by the Navier–Stokes equations, and the inlets S_i1, S_i2 and outlets S_o1, S_o2 are set with fully developed laminar flow profiles and have the same flow rate Q. A no-slip boundary condition is imposed on the channel walls.

Fig. 1 (a) Geometry of the cross-slot microchannel. The shadow represents the cross-slot region. Top left inset is the three-dimensional view. (b) Illustration of the flow-induced cell membrane damage that leads to increased membrane permeability and enhanced intracellular delivery.

The cell is initially spherical with radius a and is enclosed by a viscoelastic membrane. The fluids inside and outside the cell have identical density ρ, but different viscosities μ_in and μ_out, respectively. The centre of the cell is initially located within the cross-section S_c, which is 2l away from the inlet S_i1. To mimic experiments where cells are often not perfectly aligned with the centreline of the feeding channel, a small off-centre distance of 0.02l is set along the z-direction.

2.1. Cell mechanical model

The cell is modelled as a liquid droplet enclosed by a viscoelastic membrane with a small bending stiffness. The total membrane stress τ is the sum of the elastic and viscous stresses:

τ = τ^e + τ^v. (1)

The membrane elasticity is governed by the two-dimensional Skalak (SK) law,⁴⁵ with a strain energy function

where G_s is the surface shear elasticity modulus, and I₁, I₂ are the strain invariants with I₁ = λ₁² + λ₂² − 2 and I₂ = (λ₁λ₂)² − 1. λ₁ and λ₂ are the two principal extension ratios. The constant term C is related to the membrane area dilatational modulus K_s by K_s = (1 + 2C)G_s. We set C = 1 in the present study.

The elastic stress tensor τ^e can be calculated from

τ^e = τ^e₁e₁ ⊗ e₁ + τ^e₂e₂ ⊗ e₂, (3)

whereare the two principal stresses. e₁, e₂ are their corresponding directions, which can be determined from the unit eigenvectors of the left Cauchy–Green deformation tensor G = F·F^T. F = ∂x/∂X is the deformation gradient of the deformed cell configuration x with respect to the undeformed configuration X.

The membrane viscous stress τ^v has contributions from the shear viscosity μ_s and area dilatational viscosity ³³:

where D is the membrane strain rate tensor, tr(D) is the area dilatation rate, and P = I − nn is the projection tensor of the deformed membrane, with n as the unit normal vector. In the present study, we set , so that the relaxation times corresponding to the shear and dilatational viscosities (μ_s/G_s and ) are equal.

The bending resistance of the membrane is modelled using Helfrich's formulation⁴⁶

where k_c is the bending modulus, A is the surface area, H is the mean curvature and cur₀ is the spontaneous curvature which is set as 0. In the present study, a small resistance to bending k_c = 0.004G_sa² is used to prevent the formation of membrane wrinkles.

2.2. Membrane damage model

Membrane damage in the form of pores or ruptures appears when cells are subjected to intense mechanical stresses from the surrounding flow. With the damage growing, the membrane mechanically degrades and has less strength to resist deformation. In the present study, we model this process by the CDM, that has been widely applied for the damage of biological tissues and biomaterials.^39–42,47 In the CDM, the local average damage state is represented by a continuum damage variable. For a small patch of the cell membrane, the continuum damage variable can be defined as the ratio of the damaged area δS_D to the total area δS of the membrane patch:where S̃ is the undamaged area. The damage variable d can range from 0 (undamaged state) to 1 (completely damaged).

To account for the mechanical degradation due to damage, we adopt a strain energy function under the CDM:

The damage strain energy release rate iswhich is needed to construct the damage threshold function f(Y, d), which is useful to build the damage criterion and derive the damage evolution equation. In the present study, we choose a simple form of f(Y, d), which was developed for the membrane of a bioartificial microcapsule:^47,48

f(Y, d) = Y − κ(d), (10)

where κ(d) is the function of damage, which employs the model of Marigo⁴⁸:

κ(d) = G_s(Y_D + Y_Cd), (11)

where Y_D is the damage threshold and Y_C is the hardening modulus. Note that such simple functions of damage, where κ only depends on the instantaneous d, have been widely used for the damage of viscoelastic materials.^49,50 More realistic models probably should take into account the rate of evolution of the damage variable.⁵¹

We assume that the damage evolution is irreversible and follows the Karush–Kuhn–Tucker (KKT) loading–unloading conditions:^41,47,52

f ≤ 0, ḋ ≥ 0, fḋ = 0. (12)

To illustrate the evolution of the local membrane damage variable d, we can consider a small membrane element. At a specific time, the damage strain energy release rate Y of the element is determined by its instantaneous deformation (via eqn (9)). If Y ≤ κ(d), the element is under elastic or neutral loading, which will not lead to a change of d (i.e., ḋ = 0). If the deformation is large, Y would exceed κ(d), leading to damage loading of the membrane element, and then the damage variable will increase (i.e., ḋ > 0). To ensure fḋ = 0 of KKT conditions, f must be zero. From eqn (10), this leads to κ(d) = Y_max, where Y_max is the maximum damage strain energy release rate during the damage process. The final damage evolution function isIn the present study, since we are mainly interested in the membrane damage when the cell flows through the channel cross-slot, which often takes place in milliseconds in experiments,^10,14 we have neglected the membrane recovery process that usually requires a much longer time period. To consider the membrane recovery, one could use the model of Nikfar et al.,²⁶ in which the membrane pore size reduces with time following an exponential decay function.

2.3. Model for mass transport

When the cell membrane is damaged, solutes can be transported across the membrane through the damaged area into the cell. The cell membrane permeability P measures the passive diffusion rate of the solute molecules across the membrane, and is proportional to the porosity,^53,54 which is equivalent to the damage variable d of the present model:

P = k_pd. (14)

The term k_p is a permeability coefficient which should depend on the diffusion coefficient of the solute in the bulk, membrane thickness and the size ratio between the solute molecule and the membrane pore.^55,56 In the present study, k_p is set to be a constant for simplicity.

We assume that the transport of the solute follows the convection–diffusion equation

where c is the solute concentration, and u is the fluid velocity. The term D is the diffusion coefficient, describing how quickly the solute diffuses through the bulk fluids inside and outside the cell. For small molecules D could be estimated by the Stokes–Einstein equation.⁵⁷ The mass flux across the cell membrane followswhere c^± and ∂c^±/∂n are respectively the solute concentration and the normal concentration gradient at the membrane surface. The symbols + and − denote the external and internal sides of the membrane, respectively, and n is the unit normal vector pointing from the internal side to the external side. The numerical method to deal with the Neumann boundary condition at the deforming cell membrane will be introduced in Section 3.2.

2.4. Dimensionless parameters

The present problem is governed by the following dimensionless parameters:

• the flow Reynolds number Re = 2ρUl/μ_out, where U is the average flow speed in the inflow/outflow channels;

• the capillary number Ca = μ_outU/G_s, which measures the relative importance of the fluid viscous and membrane elastic forces;

• the viscosity ratio between the fluids inside and outside the cell λ_μ = μ_in/μ_out;

• the dimensionless membrane viscosity η = μ_s/μ_ina;

• the cell confinement ratio a/l, which compares the size of the cell relative to the channel; and

• the dimensionless permeability coefficient.

To quantify the cell deformation, we mainly use the Taylor deformation parameter

where a₁ and a₃ are the maximum dimensions of the cell's y-plane projection along the x- and z-axes, respectively. When the cell is flowing in the spiral flow regime, its long axis is often out of the plane of y = 0 and its deformation deviates from a regular ellipsoid. Therefore, a normalised cell height along the z-axis, a₃/2a, is used as the indicator of the cell elongation in this flow regime.

3. Numerical methods and validations

3.1. Method for flow-induced cell deformation

The flow-induced cell deformation is resolved using a 3D immersed-boundary LBM, that has been extensively verified in our previous studies of the dynamics of microcapsules and cells in various flow geometries.^38,58–64 Here we only give a brief introduction. The fluid flow is governed by the Navier–Stokes equations which are solved using a 3D nineteen-velocity LBM with a grid size of Δx = Δy = Δz = 2l/80. At the walls of the channel branches, a second-order bounce back scheme based on interpolation⁶⁵ is used to maintain the no-slip boundary conditions. A second-order non-equilibrium extrapolation method⁶⁶ is employed to impose the velocity boundary condition at the inlets and outlets.

The fluid–cell interaction is addressed using an immersed-boundary method.⁶⁷ The cell membrane is discretised into 8192 flat triangular elements, which are connected by 4098 nodes, following Ramanujan and Pozrikidis.⁶⁸ The immersed-boundary method ensures that the cell membrane moves at the same velocity as the fluid surrounding it, thereby satisfying the no-slip boundary condition. We use the approach of Yazdani and Bagchi³⁴ to calculate the viscoelastic force of the membrane. This involves a slightly modified mechanical system to approximate eqn (1). Further information on the implementation and validation can be found in Yazdani and Bagchi³⁴ and Wang et al.⁶⁴ The bending force density can be derived from the bending energy formulation (eqn (6)), and we follow the approaches of Garimella and Swartz⁶⁹ and Yazdani and Bagchi.⁷⁰

3.2. Numerical method for mass transport

The mass transport is governed by the convection–diffusion equation, which we solve using a 3D seven-velocity D3Q7 LB model, with a hybrid regularisation (HR) collision scheme. Details can be found in the Appendix A. To deal with the Neumann boundary condition for the mass flux at the cell membrane, we extend the 2D LB model of Zhang and Misbah³⁰ designed for ATP release from RBCs, to 3D simulations.

We describe the implementation of the Neumann boundary conditions with the aid of Fig. 2, which illustrates the boundary lattice points in a 2D cross-section of the present problem. For lattice points close to the cell membrane, represented by the square symbols, the streaming process in the LB model is affected by the existence of the membrane. For example, in Fig. 2 at the interior boundary lattice x, its concentration distribution function along the direction of c_i at time t + Δt, g_i(x, t + Δt), cannot be determined by streaming of the post-collision distribution function , due to the presence of the membrane. Therefore g_i(x, t + Δt) needs to be constructed. Inspired by the half-way bounce-back scheme of Huang et al.⁷¹ and Zhang and Misbah³⁰ proposed that the solute concentration gradient along the direction c_ī at the middle point M, α^mid, can be used to determine g_i(x, t + Δt), following

where ī is the opposite direction of i, and ĉ_i = c_i/|c_i| is the unit vector. The terms u′ and D′ are the dimensionless velocity and diffusion coefficient, respectively, defined as

Fig. 2 Schematic diagram of the lattice points close to the cell membrane. The full and empty squares are the interior and exterior boundary lattice points, respectively. The solid circles are the intersection nodes between the cell membrane and the fluid lattice mesh. The full diamonds represent the middle points of the fluid lattice.

The concentration gradient α^mid can be constructed as

where t₁ and t₂ are the tangential directions and are mutually perpendicular. Following Zhang and Misbah,³⁰ the normal and tangential derivatives on the right-hand side of eqn (20) can be obtained from the following equations:where b is the normalised distance between the boundary lattice point at x and the boundary node B, as shown in Fig. 2. α = ∂c⁻/∂n is the Neumann boundary condition for the solute mass flux at the cell membrane, which can be obtained from eqn (16).

With the cell membrane moving, some lattice points of the fluid domain would shift from one side of the membrane to the other side. The concentration distribution functions of these points are constructed from those of the adjacent points on the same side of the membrane, using the linear interpolation scheme of Zhang and Misbah.³⁰

3.3. Validation of the method for mass transport

We first validate our implementation of the numerical model for mass transport by considering the time-dependent diffusion of a solute across the membrane of a non-deformable spherical capsule in the absence of fluid flow, a problem that has been studied by Amiri and Zhang⁷² with a LBM. The computational domain has dimensions of L × W × H = 250 × 80 × 80, and the spherical capsule has a radius R = 10, placed at the centre of the domain. The diffusion coefficient D is set as 1/80 and 1/160, respectively, inside and outside the capsule. The membrane permeability P = 1/2400. The solute concentration at the top and bottom surfaces of the domain is fixed at zero, and periodic boundary conditions are applied on surfaces along the other two directions. At time t = 0, the initial solute concentration inside the capsule is c₀, and is zero elsewhere in the computational domain. The solute will diffuse through the capsule membrane and eventually be absorbed by the top and bottom surfaces. In Fig. 3, we present the time evolution of the volume-averaged solute concentration c̄ inside the capsule, and compare the result with that of Amiri and Zhang.⁷² Very good agreement has been achieved.

Fig. 3 Time evolution of the volume-averaged solute concentration c̄ inside the capsule. Symbols are the results of Amiri and Zhang.⁷²

To further test the numerical method for mass transport for moving boundary problems, we consider the time-dependent diffusion of a solute into a non-deformable spherical capsule with a permeable membrane that is advected by a flow with a constant speed u₀ = (U, 0, 0), as shown in Fig. 4. The solute concentration outside the capsule remains constant at c_∞, and the capsule membrane permeability P = 3 × 10⁻⁴. Inside the capsule, the initial solute concentration is zero, and the solute diffusion coefficient is D = 0.000781. It can be deduced from the Galilean invariance that the solute concentration inside the capsule is independent of the advection velocity of the surrounding fluid.

Fig. 4 Time-dependent diffusion of a solute into a non-deformable spherical capsule with a permeable membrane that is advected by a flow with a constant speed. The solute concentration outside the capsule remains constant at c_∞, while the solute concentration inside the capsule c(r, t) increases with time towards c_∞.

The volume-averaged solute concentration inside the capsule can be obtained analytically:

The solute concentration c(r, t) within the capsule is given in a number of textbooks (e.g. Bergman et al.⁷³)andThe ζ_n is the positive roots of

In Fig. 5, we present the time evolutions of the volume-averaged solute concentration c̄ within a stationary and a moving capsule. The results are visually identical, and both agree very well with the analytical solution.

Fig. 5 Comparison of the time evolutions of the volume-averaged solute concentration in a moving and a stationary capsule obtained from numerical simulations. Symbols are the analytical results.

3.4. Verification of the method for membrane damage

In the present study, we have employed the model developed by Grandmaison et al.⁴⁷ for the flow-induced membrane damage. To verify our numerical implementation of the model, we consider the membrane damage of an initial spherical capsule with a hyperelastic neo-Hookean membrane in simple shear flow in the Stokes regime, which was studied by Grandmaison et al.⁴⁷ using a boundary integral method. The fluids inside and outside the capsule are assumed to have equal density and viscosity. The governing dimensionless parameter for the capsule deformation is the capillary number, defined as Ca = μ_outKa/G_s, where K is the shear rate. In the membrane damage model, the damage threshold Y_D = 0.2 and the hardening modulus Y_C = 2. After being released into the simple shear flow, the capsule deforms into a steady ellipsoidal shape with a certain inclination with respect to the flow direction, and its membrane rotates around the fluid inside (i.e., tank-treading motion). As the capillary number increases, the capsule becomes more deformed. Membrane damage appears and the maximum damage variable d^∞_max increases almost linearly with Ca until Ca ≈ 0.6. A further increase of Ca results in that d^∞_max grows drastically, and at Ca ∼ 0.72, complete local membrane damage is onset, where d^∞_max = 1. In Fig. 6, we present d^∞_max of the capsule membrane as a function of Ca. Our results agree very well with those of Grandmaison et al.⁴⁷

Fig. 6 Maximum damage value at steady state, d^∞_max, as a function of Ca for Y_D = 0.2 and Y_C = 2. Symbols are the result of Grandmaison et al.⁴⁷

4. A cell flows through a cross-slot microchannel

Using the comprehensive computational model proposed in the present study, in this section, we consider the flow-induced deformation, membrane damage and enhanced transport of a diffusive solute into a modelled 3D cell when it is flowing through a cross-slot microchannel. We attempt to address two open questions: firstly, under similar pressure drop across the device, is the viscous or inertial flow regime more effective in generating cell membrane damage and enhancing intracellular mass transport? Secondly, how are the membrane damage affected by the cell membrane viscosity? In experiments, cells have been flown through cross-slot microchannels at flow speeds up to metres per second,^14,74,75 leading to shear rates well above 10³ s⁻¹, where the cell membrane viscosity is expected to play an important role in determining the cell transient deformation and the associated membrane damage.

To elucidate the effect of flow regime, we conduct numerical simulations using three sets of parameters, corresponding to the same cell flowing in three distinctive flow regimes with increasing flow strength and inertial effect, under comparable pressure drop across the whole device:

• viscous flow regime: Re = 0.4, U = U₀, μ_out = 2μ_in;

• moderate-inertia flow regime: Re = 40, U = 10U₀, μ_out = μ_in/5; and

• spiral flow regime: Re = 80, U = 14U₀, μ_out = μ_in/7.

Note that in practical flow experiments for cell stretching, the cross-slot region often only has a length of tens of micrometres, representing a very small faction of the entire flow path of the whole device, compared with the straight feeding and outlet channels that are usually much longer, in tens of millimetres.^10,74 Therefore, the pressure drop across the whole device is primarily determined by the flow characteristics and lengths of the feeding and outlet channels. For the same device operating in the steady laminar flow regime, the total pressure drop is approximately proportional to the product of the fluid viscosity and average flow speed in straight channels ΔP ∝ μ_outU.⁷⁶ For the three sets of flow parameters considered in the present study, when increasing the average flow speed in the feeding channels, we reduce the fluid viscosity to ensure μ_outU = constant.

The streamlines of the background flow in the channel cross-slot region for the three distinctive flow regimes are presented in Fig. 7. From Fig. 7(a and b), one can see that the flow patterns for the viscous and moderate-inertia flow regimes are qualitatively similar, featuring extensional flows that are symmetric about the x = 0 and z = 0 planes, with a stagnation point at the centre of the cross-slot region. Our simulations suggest that when the flow Reynolds number exceeds 43, a spiral flow onsets due to a symmetry-breaking bifurcation.⁴³Fig. 7(c and d) present the streamlines of the spiral flow at Re = 80. The spiral vortex rotates about the axis of the outlet channels. It is steady and symmetric with respect to the z = 0 plane. In the following sessions, we first consider the flow-induced cell membrane damage in the viscous and moderate-inertia flow regimes, focusing on the effects of inertia and cell membrane viscosity in Sections 4.1 and 4.2, respectively. We then consider the spiral flow regime in Section 4.3.

Fig. 7 Front views of streamlines of the background flow starting from the feeding channel in the plane y = 0 at Re = (a) 0.4, (b) 40, and (c) 80. (d) Top view of streamlines at Re = 80. The colour indicates the velocity magnitude.

4.1. A cell with a hyperelastic membrane

We start from considering a simple cell where its membrane is purely hyperelastic that follows the SK law. We do so mainly for two reasons: firstly, the study serves as an important baseline reference for more realistic cells that have a viscoelastic membrane; secondly, the membrane damage and enhanced cross-membrane transport for a simple hyperelastic cell remains largely unknown. We conduct numerical simulations of a cell with a/l = 0.4, , under the viscous and moderate-inertia flow regimes at Re = 0.4 and 40, respectively.

Fig. 8(a) shows the time evolution of the Taylor shape parameter D_XZ of a cell in the viscous and moderate-inertia flow regimes at Ca = 0.2. In this and other figures of the present study we set t = 0 as the time when a cell begins to enter the cross-slot region, i.e., any membrane element reaches the plane at x = −l. The time when the entire cell leaves the cross-slot region is marked with a circle symbol. To aid the understanding of the result, we also present the instantaneous shapes of the cell in Fig. 9 and 10 for the two flow regimes, respectively. In the viscous flow regime at Re = 0.4, the cell reaches a quasi-steady ellipsoidal shape with its long axis aligned with the direction of extension when it is approaching the stagnation point (see Fig. 9). The Taylor shape parameter remains a constant for a considerable time period. With a moderate inertial effect at Re = 40, after entering the cross-slot region, the cell quickly reaches its maximum elongation at t₂ and then starts to retract, even when the cell is still approaching the stagnation point. This was termed an “overshoot–retract” motion and was analysed in detail in an earlier study by the present authors.³⁸ In the present study, the viscosity of the fluid inside the cell is much higher than that of Lu et al.,³⁸ and therefore only a mild overshoot–retract motion has been observed here. In Fig. 8(a), comparing the cell deformation in the two flow regimes, it is clear that the cell experiences much larger deformation under the effect of inertia.

Fig. 8 Time evolutions of the (a) Taylor deformation parameter D_XZ, (b) area-averaged membrane damage variable d̄ and (c) volume-averaged solute concentration within the cell c̄, when the cell is flowing through the channel cross-slot region in the viscous and moderate-inertia flow regimes, at Re = 0.4 and 40, respectively, with a/l = 0.4, , Ca = 0.2, Y_D = 0.02 and Y_C = 2. t₁, t₂, and t₃ are three dimensionless times at 0.30, 0.61 and 0.90, respectively, when the instantaneous cell profiles will be shown in Fig. 9 and 10. At t = 0 the cell starts to enter the channel cross-slot; the circle symbols mark the moments when the cell completely leaves the region. The cell has a hyperelastic membrane. In (a), the temporal evolutions of D_XZ of the same cell without taking into account the membrane damage are also presented as references.

Fig. 9 Instantaneous membrane damage profiles of the cell of Fig. 8 at Re = 0.4 as seen from (a) the y-axis and (b) the x-axis. The solid lines in (a) are the trajectories of the cell's mass centre. The time instances are provided in Fig. 8.

Fig. 10 Instantaneous membrane damage profiles of the cell of Fig. 8 at Re = 40 as seen from (a) the y-axis and (b) the x-axis. The solid lines in (a) are the trajectories of the cell's mass centre. The time instances are provided in Fig. 8.

We also present the distributions of the membrane damage variable of the cell in the viscous and moderate-inertia flow regimes in Fig. 9 and 10, respectively. In the viscous flow regime, the maximum membrane damage occurs at the membrane points that have the maximum y-axis values (see Fig. 9(a) at t₂), due to the highest membrane tension there (not shown). Interestingly, when the inertial effect becomes significant, as can be seen from Fig. 10(b) at t₂, the locations for maximum membrane damage have shifted to the central regions of the membrane faces that are perpendicular to the direction of the inflows (i.e., x-direction). The result is not surprising. In the inertial flow regime, the membrane elastic force is generated mainly to balance to the fluid inertial force, which is proportional to ρU² and therefore has the maximum magnitude in the core of the inflow that has the highest velocity.

Fig. 8(b) shows the time evolutions of the area-averaged membrane damage variable d̄ of the cell in the two flow regimes. After the cell enters the cross-slot region, the membrane damage increases to peak values, when the cell deformation is at the maximum, and then remains unchanged due to the condition of ḋ ≥ 0 in the present membrane damage model. Because of the greater cell deformation in the moderate-inertia flow regime compared with the viscous regime, d̄_max at Re = 40 is almost four times larger than that of Re = 0.4, which is expected to result in much faster intracellular mass transport. We also consider the effect of the cell membrane damage on its transient deformation. In Fig. 8(a), one can compare the time evolutions of the Taylor shape parameter of the same cell with and without taking into account the membrane damage. It can be found that in the viscous regime, due to the relatively small membrane damage, there is visually no difference between the results. However, in the moderate-inertia region, d̄ reaches the order of 0.1, which has led to a considerable increase of D_XZ, compared with the cell whose membrane cannot be damaged by the fluid flow.

Note that the present computational model can predict the instantaneous field of the solute concentration inside the cell. Two examples are presented in Fig. 11 and 12, corresponding to the cell in Fig. 9 and 10, respectively. From the figures, it is seen that within the cell the regions with relatively high solute concentration are always close to membrane areas with a higher level of membrane damage. Fig. 8(c) presents the time evolutions of the volume-averaged solute concentration within the cell in the two flow regimes. Surprisingly, although the cell at Re = 40 has developed much greater membrane damage, leading to faster cross-membrane mass transport, its internal solute concentration at the exit of the cross-slot region is much lower than that of the cell at Re = 0.4. The intriguing result is mainly due to the difference of the residence time of the cell in the channel cross-slot in the two flow regimes. In the dimensional form, the cell residence time at Re = 40 is only about a tenth of that at Re = 0.4.

Fig. 11 Instantaneous solute concentration inside the cell of Fig. 9 in the middle planes along (a) the y-axis and (b) the x-axis.

Fig. 12 Instantaneous solute concentration inside the cell of Fig. 10 in the middle planes along (a) the y-axis and (b) the x-axis.

With all other cell and flow parameters kept unchanged, we vary the cell membrane shear elasticity G_s and consider its effect on the maximum cell deformation and membrane damage in the channel cross-slot, as well as the volume-averaged solute concentration c̄_l in the cell at the exit of the cross-slot region in the viscous and moderate-inertia flow regimes. The results are presented in Fig. 13.

Fig. 13 Effect of the cell membrane shear elasticity on the (a) maximum Taylor deformation parameter D^max_XZ; (b) maximum area-averaged membrane damage variable d̄_max and (c) volume-averaged solute concentration within the cell when it leaves the cross-slot region at Re = 0.4 and 40. Cell membrane shear elasticity is related to the capillary number by G_s = μ_outU/Ca. Other parameters are the same as those of Fig. 8. For the cell at Re = 40, we only show results for Ca ≤ 0.21, since at even high Ca the cell will have complete local membrane damage with d_max = 1. From Fig. 13, we can find that the general conclusions drawn from Fig. 8 are robust with respect to the cell membrane shear elasticity. For all G_s values considered, the cell always has larger deformation and greater membrane damage when flowing through the channel cross-slot region in the inertial flow regime. However, this does not necessarily lead to increased intracellular solute transport due to the much shorter dimensional residence time in the channel cross-slot.

4.2. Effect of cell membrane viscosity

In this section, we consider a cell with a viscoelastic membrane, focusing on the effect of the membrane viscosity on the membrane damage and the associated intracellular solute transport. The membrane elasticity and its damage follow the same models with identical parameters that are used in Fig. 8 of Section 4.1, and we add the membrane viscosity using eqn (1). We conduct numerical simulations with increasing cell membrane viscosity η from 0 to 40. Note that in an earlier study, membrane viscosity on the order of η = 10 was reported for human leukaemia cells⁶⁴ flowing through a constricted microchannel.

Fig. 14(a–c) respectively show the time evolutions of the Taylor deformation parameter D_XZ, the area-averaged membrane damage variable d̄ and the volume-averaged solute concentration inside the cell c̄ for the cell of Fig. 8 with increasing membrane viscosity in both flow regimes.

Fig. 14 Time evolutions of (a) the Taylor deformation parameter D_XZ; (b) area-averaged membrane damage d̄; and (c) volume-averaged concentration c̄ of the cell of Fig. 8 with increasing membrane viscosity η in the viscous (solid lines) and moderate-inertia (dashed lines) flow regimes, at Re = 0.4 and 40, respectively.

In general, the membrane viscosity slows down the cell deformation. This has much more significant effects when the cell is flowing in the moderate-inertia regime, where the flow speed is much higher and the residence time of the cell in the channel cross-slot is much shorter, compared with the viscous regime. The slow deformation and short residence time have led to the a cell with high membrane viscosity reaching a much lower maximum deformation in the channel cross-slot, as can be seen from Fig. 14(a). The results of the area-averaged membrane damage variable d̄ follow a similar trend, showing that the membrane viscosity can significantly reduce the cell membrane damage in the inertial flow regime, leading to much less solute entering the cell in the channel cross-slot (see Fig. 14(c)). Interestingly, we find that with η ≥ 20 the maximum Taylor deformation parameter and area-averaged membrane damage variable of a cell at Re = 40 have both become lower than those of the same cell at Re = 0.4. This is in contrast to the results observed from hyperelastic cells of Section 4.1, for which the cell deformation and membrane damage are always larger in the inertial flow regime.

4.3. Spiral flow regime

We also consider a cell with a viscoelastic membrane flowing through the channel cross-slot region under the spiral flow regime at Re = 80. The undisturbed background flow has been shown in Fig. 7(c and d). Fig. 15 presents the instantaneous shape deformation and membrane damage profiles of the cell of Fig. 14 with η = 10. In the figure, the trajectory of the cell's mass centre is marked by a solid line. Compared with the viscous and moderate-inertia flow regimes (see Fig. 9 and 10), a distinct feature of the cell motion in the spiral flow regime is that its trajectory is helical. Specifically, rather than directly flowing through the channel cross-slot region via the close proximity of the stagnation point at the centre region in the viscous and moderate-inertia flow regimes, the cell rotates around the z-axis for about 540° before leaving the cross-slot in the spiral flow at Re = 80. This helical trajectory, compared with a directly flow-through scenario, helps to increase the cell's residence time in the cross-slot, giving it more time to deform under the flow extension.

Fig. 15 Instantaneous membrane damage profiles of the cell of Fig. 14 with η = 10 in the spiral flow regime at Re = 80, as seen from (a) the y-axis and (b) the x-axis. The solid lines are the trajectories of the cell's mass centre. t₁, t₂, and t₃ are three dimensionless times at 0.30, 0.77 and 1.50, respectively.

Fig. 16(a–c) present the time evolutions of a₃/(2a), the area-averaged membrane damage variable d̄ and the volume-averaged solute concentration inside the cell of Fig. 15 in all three flow regimes. Under the combined effects of the strongest inertial force and the increased residence time due to the helical flow trajectory, the cell flowing at Re = 80 has the maximum deformation and area-averaged membrane damage. These have also led to a slightly higher volume-averaged solute concentration inside the cell at the exit of the channel cross-slot, compared with the cell at Re = 40.

Fig. 16 Time evolutions of the (a) normalised cell height , (b) area-averaged membrane damage variable d̄ and (c) volume-averaged solute concentration c̄ inside the cell of Fig. 15 in the three flow regimes at Re = 0.4, 40 and 80, respectively. The three time instances t₁, t₂, and t₃ are from Fig. 15 where the instantaneous cell profiles are shown.

For the cell of Fig. 16, we increase its membrane viscosity to η = 40, while keeping all other parameters the same, and present the results in Fig. 17. At Re = 80, despite the strongest inertial force and increased residence time, both the maximum cell deformation and area-averaged membrane damage have fallen below those of the cell flowing at Re = 0.4 in the viscous flow regime. The results therefore confirm our early conclusion that the cell membrane viscosity tends to protect the cell from large deformation and membrane damage when it is flowing through the channel cross-slot, in particular in the inertial flow regime.

Fig. 17 Time evolutions of the (a) normalised cell height , (b) area-averaged membrane damage variable d̄ and (c) volume-averaged solute concentration c̄ inside the cell of Fig. 14 with a higher membrane viscosity η = 40 at Re = 0.4, 40 and 80.

In Fig. 18 we summarise the maximum area-averaged membrane damage of the cell considered in the present study with different cell membrane viscosity values in the three flow regimes that are driven by the same pressure drop across the whole device. The results reveal distinct membrane damage features that are strongly affected by the membrane viscosity.

Fig. 18 Maximum area-averaged membrane damage d̄_max as a function of membrane viscosity η in the three flow regimes at Re = 0.4, 40 and 80 respectively.

The results suggest that for cells with low membrane viscosity, e.g., η ≤ 10, driving low-viscosity suspension at a high speed in the inertial flow regime will lead to higher level of cell deformation and membrane damage. However, the enhancements can be significantly reduced or even reversed with an increase of the cell membrane viscosity. For cells with high membrane viscosity, e.g., η = 40, operating experiments with a high-viscosity suspension can be more effective in enhancing cell deformation and membrane damage.

5. Conclusions

A novel three-dimensional computational framework has been proposed in the present study, which can, for the first time, simultaneously take into account the transient flow-induced deformation of a viscoelastic biological cell, the deformation-induced cell membrane damage, as well as the damage-resulted membrane mechanical weakening and enhanced cross-membrane intracellular mass transport. With the computational framework, we have considered an initially spherical cell flowing through a cross-slot microchannel, focusing on the effects of the flow regimes, under comparable pressure drop across the whole device, and the cell membrane viscosity on the cell deformation, membrane damage and enhanced intracellular mass transport. Our numerical simulations have generated novel and interesting results, suggesting distinct cell dynamics and membrane damage features that strongly depend on the cell membrane viscosity. We find that for cells with relatively low membrane viscosity, e.g., η ≤ 10, operating experiments in the inertial flow regime, that can be achieved by driving a low-viscosity suspension at a high speed, often results in larger cell deformation and membrane damage, due to the effect of the flow inertia. Increasing the cell membrane viscosity slows down the cell deformation, and can lead to considerably smaller cell deformation and membrane damage, when a cell flows through the cross-slot region within a short time period and therefore does not have sufficient time to deform. Thus, for cells with high membrane viscosity, e.g., η = 40, in order to enhance the deformation-induced membrane damage, it could be advantageous to operate experiments in the viscous flow regime, where one drives a high-viscosity suspension at a low speed, offering sufficient residence time for the cell to develop its deformation in the channel cross-slot region. The present simulation results may provide useful guidelines to practical experiments, which flow cells through cross-slot microchannels to temporarily damage cell membranes to enhance intracellular drug delivery.

As a general computational framework for modelling the flow-induced cell membrane damage and the resulted enhanced intracellular mass transport, the present model will need to be validated against carefully designed experiments in future. Such validations are presently not possible, due to the lack of quantitative experimental data, that record the time evolutions of the cell membrane damage and solute concentration inside and around a cell. Once validated, the computational framework may be used in a wide range of applications, ranging from designing microfluidic devices to physically disrupt the cell membrane for enhanced drug delivery, to optimising flow conditions to minimise cell damage in cell printers and other inertial microfluidics.

Author contributions

Y. S. and R. X. L. designed the computational model and the research problem. R. X. L. wrote the codes and conducted the numerical simulations. All authors analysed the results, prepared the figures, wrote and reviewed the manuscript.

Conflicts of interest

There are no conflicts to declare.

Appendix

A. LBM model for mass transport

We use the D3Q7 model (three-dimensional seven discrete micro velocities) to formulate a convection–diffusion LBM model for mass transport. The discrete micro velocities are defined asThe terms Δx and Δt represent the spatial and temporal mesh sizes, respectively.

Similar to the LBM used in the fluid flow, fictive particles with velocity c_i can propagate and collide on a discrete lattice mesh and form a concentration distribution function g_i(x, t) at the position x and time t. The temporal evolution of the concentration distribution function can be divided into the collision and streaming processes:

• the collision process

• and the streaming process

Here g^neq_i is the non-equilibrium distribution and τ_g is the dimensionless relaxation time which is related to the diffusion coefficient by D = (τ_g − 1/2)Δtc_s², where c_s = 1/2(Δx/Δt) is the sound speed in the D3Q7 model. The macro mass concentration c can be obtained from the distribution function by

The LBM becomes unstable when the dimensionless relaxation time is close to 1/2 with the BGK collision. To avoid the numerical instability, we reconstruct the non-equilibrium distribution g^neq_i, following the method of Zhang et al.,⁷⁷ by expanding the non-equilibrium distribution with the Hermite polynomial:

where w_i is the weight factor, H⁽ⁿ⁾ is the Hermite polynomial and a⁽ⁿ⁾₁ is the expansion coefficient of the non-equilibrium distribution. In our simulation, we expand g^neq_i to the first-order, i.e., n = 1. The first-order Hermite polynomials are H⁽¹⁾_iα = c_iα, and the first-order expansion coefficient a⁽¹⁾_1,α can be obtained by directly projecting the non-equilibrium distribution function:This first-order expansion coefficient a⁽¹⁾_1,α can also be calculated from the macroscopic quantity through the Chapman–Enskog analysis:⁷⁷In eqn (32), ∇c is computed with the second-order centred finite differences scheme, and the time derivative term ∂_t(cu_α) is computed with an explicit Euler scheme. It was shown in an earlier study that reconstructing the expansion coefficient a⁽¹⁾_1,α using a hybrid regularisation (HR) approach,⁷⁷ in the form ofcan lead to better numerical stability and accuracy, compared with using eqn (31) or (32) alone. The term β ∈ [0, 1] is a free parameter, and is set to 0.95 in the present study.

Acknowledgements

R. X. L. acknowledges the PhD studentship provided by the QMUL and the Chinese Scholarship Council. The simulations were performed using the high-performance computer clusters of QMUL (funded by the UK EPSRC grant EP/K000128/1, EP/P020194/1 and EP/T022213/1). Y. S. acknowledges financial support from the UK Royal Society (IES/R2/212075).

Notes and references

1. M. P. Stewart, R. Langer and K. F. Jensen, Chem. Rev., 2018, 118, 7409–7531.
2. F. J. Cunningham, N. S. Goh, G. S. Demirer, J. L. Matos and M. P. Landry, Trends Biotechnol., 2018, 36, 882–897.
3. A. Wadhwa, A. Aljabbari, A. Lokras, C. Foged and A. Thakur, Pharmaceutics, 2020, 12, 102.
4. J. Hur and A. J. Chung, Adv. Sci., 2021, 8, 2004595.
5. A.-R. Shokouhi, Y. Chen, H. Z. Yoh, J. Brenker, T. Alan, T. Murayama, K. Suu, Y. Morikawa, N. H. Voelcker and R. Elnathan, Adv. Mater., 2023, 2304122.
6. H. Li, T. Y. Tsui and W. Ma, Int. J. Mol. Sci., 2015, 16, 19518–19536.
7. J. Yang, A. Bahreman, G. Daudey, J. Bussmann, R. C. Olsthoorn and A. Kros, ACS Cent. Sci., 2016, 2, 621–630.
8. D. J. Stevenson, F. J. Gunn-Moore, P. Campbell and K. Dholakia, J. R. Soc., Interface, 2010, 7, 863–871.
9. Y. Cao, E. Ma, S. Cestellos-Blanco, B. Zhang, R. Qiu, Y. Su, J. A. Doudna and P. Yang, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 7899–7904.
10. M. E. Kizer, Y. Deng, G. Kang, P. E. Mikael, X. Wang and A. J. Chung, Lab Chip, 2019, 19, 1747–1754.
11. A. Sharei, J. Zoldan, A. Adamo, W. Y. Sim, N. Cho, E. Jackson, S. Mao, S. Schneider, M.-J. Han and A. Lytton-Jean, et al. , Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 2082–2087.
12. A. Uvizl, R. Goswami, S. D. Gandhi, M. Augsburg, F. Buchholz, J. Guck, J. Mansfeld and S. Girardo, Lab Chip, 2021, 21, 2437–2452.
13. C. Kwon and A. J. Chung, Lab Chip, 2023, 23, 1758–1767.
14. G. Kang, D. W. Carlson, T. H. Kang, S. Lee, S. J. Haward, I. Choi, A. Q. Shen and A. J. Chung, ACS Nano, 2020, 14, 3048–3058.
15. 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.
16. O. S. Pak, Y. N. Young, G. R. Marple, S. Veerapaneni and H. A. Stone, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 9822–9827.
17. E. A. Evans, R. Waugh and L. Melnik, Biophys. J., 1976, 16, 585–595.
18. M. Lokhandwalla and B. Sturtevant, Phys. Med. Biol., 2001, 46, 413.
19. S. T. Braakman, R. M. Pedrigi, A. T. Read, J. A. Smith, W. D. Stamer, C. R. Ethier and D. R. Overby, Exp. Eye Res., 2014, 127, 224–235.
20. J. M. Meacham, K. Durvasula, F. L. Degertekin and A. G. Fedorov, Sci. Rep., 2018, 8, 3727.
21. M. Razizadeh, M. Nikfar, R. Paul and Y. Liu, Biophys. J., 2020, 119, 471–482.
22. P. L. McNeil and R. A. Steinhardt, Annu. Rev. Cell Dev. Biol., 2003, 19, 697.
23. S. Sohrabi and Y. Liu, Artif. Organs, 2017, 41, E80–E91.
24. M. Nikfar, M. Razizadeh, R. Paul and Y. Liu, Microfluid. Nanofluidics, 2020, 24, 1–13.
25. M. Nikfar, M. Razizadeh, J. Zhang, R. Paul, Z. J. Wu and Y. Liu, Artif. Organs, 2020, 44, E348–E368.
26. M. Nikfar, M. Razizadeh, R. Paul, Y. Zhou and Y. Liu, Biomicrofluidics, 2021, 15, 044102.
27. K. Islam, M. Razizadeh and Y. Liu, Phys. Chem. Chem. Phys., 2023, 25, 12308–12321.
28. Z. Y. Luo and B. F. Bai, Phys. Fluids, 2019, 31, 121902.
29. H. Zhang, Z. Shen, B. Hogan, A. I. Barakat and C. Misbah, Biophys. J., 2018, 115, 2218–2229.
30. H. Zhang and C. Misbah, Comput. Fluids, 2019, 187, 46–59.
31. Z. Gou, H. Zhang, M. Abbasi and C. Misbah, Biophys. J., 2021, 120, 4819–4831.
32. Z. Gou, H. Zhang and C. Misbah, J. R. Soc., Interface, 2023, 20, 20230186.
33. D. Barthès-Biesel and H. Sgaier, J. Fluid Mech., 1985, 160, 119–135.
34. A. Yazdani and P. Bagchi, J. Fluid Mech., 2013, 718, 569–595.
35. D. Cordasco and P. Bagchi, Phys. Fluids, 2017, 29, 041901.
36. F. Guglietta, M. Behr, L. Biferale, G. Falcucci and M. Sbragaglia, Soft Matter, 2020, 16, 6191–6205.
37. P. Li and J. Zhang, Cardiovasc. Eng. Technol., 2021, 12, 232–249.
38. R. Lu, Z. Guo, P. Yu and Y. Sui, J. Fluid Mech., 2023, 962, A26.
39. J. Hokanson and S. Yazdani, Mech. Res. Commun., 1997, 24, 151–159.
40. A. Natali, P. Pavan, E. L. Carniel, M. Lucisano and G. Taglialavoro, Med. Eng. Phys., 2005, 27, 209–214.
41. X. Deng, A. Korobenko, J. Yan and Y. Bazilevs, Comput. Methods Appl. Mech. Eng., 2015, 284, 349–372.
42. G. A. Holzapfel and B. Fereidoonnezhad, Biomechanics of living organs, Elsevier, 2017, pp. 101–123.
43. S. J. Haward, R. J. Poole, M. A. Alves, P. J. Oliveira, N. Goldenfeld and A. Q. Shen, Phys. Rev. E, 2016, 93, 031101.
44. K. Kechagidis, B. Owen, L. Guillou, H. Tse, D. Di Carlo and T. Krüger, Microsyst. Nanoeng., 2023, 9, 100.
45. R. Skalak, A. Tozeren, R. Zarda and S. Chien, Biophys. J., 1973, 13, 245–264.
46. O.-Y. Zhong-Can and W. Helfrich, Phys. Rev. A: At., Mol., Opt. Phys., 1989, 39, 5280.
47. N. Grandmaison, D. Brancherie and A.-V. Salsac, J. Fluid Mech., 2021, 914, A25.
48. J.-J. Marigo, C. R. Acad. Sci. Paris, 1981, 292, 1309–1312.
49. J. C. Simo, Comput. Methods Appl. Mech. Eng., 1987, 60, 153–173.
50. E. Peña, B. Calvo, M. Martnez and M. Doblaré, Int. J. Numer. Methods Eng., 2008, 74, 1198–1218.
51. Y. Hamiel, Y. Liu, V. Lyakhovsky, Y. Ben-Zion and D. Lockner, Geophys. J. Int., 2004, 159, 1155–1165.
52. S. Hosseini, J. J. Remmers and R. De Borst, Int. J. Numer. Methods Eng., 2014, 98, 391–398.
53. K. M. Hämäläinen, K. Kontturi, S. Auriola, L. Murtomäki and A. Urtti, J. Controlled Release, 1997, 49, 97–104.
54. O. Hosoya, S. Chono, Y. Saso, K. Juni, K. Morimoto and T. Seki, J. Pharm. Pharmacol., 2004, 56, 1501–1507.
55. E. M. Renkin, J. Gen. Physiol., 1954, 38, 225.
56. X. Lin, Q. Yang, L. Ding and B. Su, ACS Nano, 2015, 9, 11266–11277.
57. J. T. Edward, J. Chem. Educ., 1970, 47, 261.
58. Y. Sui, Y. T. Chew, P. Roy, Y. P. Cheng and H. T. Low, Phys. Fluids, 2008, 20, 112106.
59. Y. Sui, Y.-T. Chew, P. Roy and H.-T. Low, J. Comput. Phys., 2008, 227, 6351–6371.
60. Y. Sui, X. Chen, Y. Chew, P. Roy and H. Low, Comput. Fluids, 2010, 39, 242–250.
61. Z. Wang, Y. Sui, A.-V. Salsac, D. Barthès-Biesel and W. Wang, J. Fluid Mech., 2016, 806, 603–626.
62. Z. Wang, Y. Sui, A.-V. Salsac, D. Barthès-Biesel and W. Wang, J. Fluid Mech., 2018, 849, 136–162.
63. R. Lu, Z. Wang, A.-V. Salsac, D. Barthès-Biesel, W. Wang and Y. Sui, J. Fluid Mech., 2021, 923, A11.
64. Z. Wang, R. Lu, W. Wang, F. Tian, J. Feng and Y. Sui, Biomech. Model. Mechanobiol., 2023, 1–15.
65. M. Bouzidi, M. Firdaouss and P. Lallemand, Phys. Fluids, 2001, 13, 3452–3459.
66. Z. Guo, C. Zheng and B. Shi, Phys. Fluids, 2002, 14, 2007–2010.
67. C. S. Peskin, Acta Numer., 2002, 11, 479–517.
68. S. Ramanujan and C. Pozrikidis, J. Fluid Mech., 1998, 361, 117–143.
69. R. V. Garimella and B. K. Swartz, Curvature estimation for unstructured triangulations of surfaces, Los Alamos National Laboratory Technical Report LA-UR-03-8240, 2003.
70. A. Yazdani and P. Bagchi, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 056308.
71. J. Huang, Z. Hu and W.-A. Yong, J. Comput. Phys., 2016, 310, 26–44.
72. F. A. Amiri and J. Zhang, Int. Commun. Heat Mass Transf., 2021, 128, 105601.
73. T. L. Bergman, A. S. Lavine, F. P. Incropera and D. P. DeWitt, Introduction to heat transfer, John Wiley & Sons, 2011, pp. 303–304.
74. D. R. Gossett, H. T. Tse, S. A. Lee, Y. Ying, A. G. Lindgren, O. O. Yang, J. Rao, A. T. Clark and D. Di Carlo, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 7630–7635.
75. F. J. Armistead, J. G. De Pablo, H. Gadêlha, S. A. Peyman and S. D. Evans, Biophys. J., 2019, 116, 1127–1135.
76. C. Pozrikidis, Introduction to theoretical and computational fluid dynamics, Oxford University Press, 2011.
77. Z. Zhang, Z. Li and Y. Wu, Front. Phys., 2022, 10, 875628.

---