Revisiting: Shape transitions of red blood cell under oscillatory flows in microchannels

A. Problem formulation

The red blood cell (RBC) is considered as a deformable body, which is immersed in the fluid plasma Ω_f, as shown in Fig. 1(a). The RBC consists of two components: (i) the cellular membrane (Γ_RBC) and (ii) the cytosol fluid (Ω_c). The outer surface of the RBC’s membrane (Γ_RBC) is the interface between the RBC and the fluid plasma, which is considered as a Newtonian incompressible fluid. Inside the RBC, the cytosol fluid occupies the inner part of the domain Ω_c. The viscosity of the cytosol fluid is different from the one of the plasma fluid (viscosity contrast- $\frac{{\nu }_c}{{\nu }_f}\ne 1$). The RBC deforms under the loading from the fluid domain Ω_f onto the interface Γ_RBC. In our hybrid continuum–particle approach, different numerical methods are used to model the dynamics of Γ_RBC, Ω_c, and Ω_f, as explained in the following sections.

B. The idealized shape of the RBC—Γ_RBC

C. RBC membrane model—Γ_RBC

In this work, the model for membrane mechanics follows closely standard dissipative particle dynamics (DPD) methodologies.^30–33 The details of the respective potentials are listed in Appendix A 1.

2. Cellular membrane/cytoskeleton self-interaction

Due to its mechanical properties, the RBC membrane is highly deformable and it can come into contact with its own surface if the external stretching force is sufficiently large locally. To account for this interaction, the bilayer–cytoskeletal interaction force F^E is incorporated into the total RBC membrane forces.³⁶ The membrane thickness of a healthy RBC is typically 10 nm, composed of an outer layer (lipid bilayer) and the inner layer (skeleton).³⁷ In this study, F^E is applied when the distance between two elements of Γ_RBC with opposite normal vectors is less than d_a. Note that d_a is a parameter of the model, which plays its role only when the cell membrane comes into contact with itself. In this work, d_a is chosen as 20 times the thickness of the membrane (d_a = 0.2 µm). A penalty approach is adopted to compute F^E, which is applied equally to all the vertices (i = 1, 2 and 3) for each of the two elements. We assume that the cytoplasmic pressure increases locally as two adjacent elements come significantly closer together. We further assume that the maximum intracellular pressure reaches 4.12 Pa, as shown in Fig. 2. The bilayer–cytoskeletal interaction force is given by

$${\mathbf{F}}_i^E=p_{bs}{A}_e{\mathbf{n}}_{\mathbf{e}},$$(4)

with a penalty coefficient for the contact traction of the bilayer–cytoskeletal p_bs = 4.12 Pa. A_e and n_e are the area and the normal vector of the element, respectively.

D. Modeling the interaction $(\mathfrak{S})$ between RBC’s membrane (Γ_RBC) and the cytosol fluid (Ω_c)

The interaction between the membrane and the cytosol fluid leads to the emergence of the overall RBC dynamics.³⁸ Differentiating between membrane and cytosol particles is essential for accurately capturing RBC mechanics, as it enables realistic modeling of membrane elasticity, internal viscosity contrast, and deformation behaviors under flow.³⁹ This interaction is modeled using the dissipative particles dynamics (DPD method). DPD is a microscopic simulation technique widely used to model the flow of complex fluids. Here, the flow is described as a group of clustered interacting particles according to the Lagrangian approach.³⁰ In this work, the cytosol within the RBC is modeled using a set of randomly distributed DPD particles (N_f) that fill the internal volume of the cell,^30,31 as shown in Fig. 1(a). Therefore, both the membrane and the cytosol fluid are considered collections of DPD particles.

Due to the different nature of the interactions, the components of the force of each particle depend on the nature of the particle i (either membrane or cytosol particle). In general, each DPD particle i interact with surrounding particles j (either a membrane or cytosol fluid particle) within a cutoff radius r_c through three pairwise additive forces: (a) the conservative force ${\mathbf{F}}_{ij}^C$; (b) the dissipative force ${\mathbf{F}}_{ij}^D$; (c) the random force ${\mathbf{F}}_{ij}^R$, and the penalty force ${\mathbf{F}}_{ij}^P$. The relative position vector between the particles i and j and related terms are given by r_ij = r_i − r_j, the distance r_ij = |r_ij|, and the unit vector ${\widehat{\mathbf{r}}}_{ij}=\frac{{\mathbf{r}}_{\mathbf{ij}}}{r_{ij}}$. In addition, v_i,j = v_i − v_j is the relative velocity between the particles i and j with velocities v_i and v_j.

The mathematical formulation of the conservative force ${\mathbf{F}}_{ij}^C$, the dissipative force ${\mathbf{F}}_{ij}^D$, the random force ${\mathbf{F}}_{ij}^R$, and the penalty force ${\mathbf{F}}_{ij}^P$ for the membrane and the cytosol fluid particles are explained in Appendix A 2.

Here, the computed force is the sum of the external and internal forces as ${\mathbf{F}}_i={\mathbf{F}}_i^{\mathit{internal}}+{\mathbf{F}}_i^{\mathit{external}}$. The external force ${\mathbf{F}}_i^{\mathit{external}}$ is either known (applied directly on the cell membrane) or computed from the extracellular fluid flows (fluid-structure interaction), as explained in the following.

E. Continuum approach for the blood plasma $(\mathfrak{F})$

The fluid solver is based on the sharp-interface curvilinear-immersed boundary (CURVIB) method (https://www.osti.gov/biblio/1312901) in a background structured-grid that contains the RBC model.²⁹ The computational grid is formulated using the hybrid staggered/non-staggered methodology.⁴² The Navier–Stokes equations are discretized using second-order accurate differencing for all spatial derivatives. The discrete equations are integrated in time via a second-order accurate fractional step method. The momentum equation is solved with an implicit Runge–Kutta solver. The total mass flux of the entire domain is enforced strictly to machine zero to satisfy the incompressibility constraint of the blood plasma. The flexible generalized Minimal residual (FGMRES) method with multigrid pre-conditioner is used to solve the Poisson equation, which results from the discrete continuity equation.⁴³

The complex surface of Γ_RBC is handled by the sharp-interface immersed boundary method, which requires only the velocity reconstruction at the immersed surface along the local normal to the body.⁴³ The presence of Γ_RBC is represented by a layer of computational grid points (immersed surface—IB nodes), which is distributed around the surface Γ_RBC. The boundary conditions reconstructed at the IB nodes are the three velocity components, which are reconstructed by interpolation.⁴⁴ In this work, a linear velocity profile along the normal direction of the immersed surface is assumed near RBC’s boundary.⁴⁵ Therefore, the presence of the IB nodes in our algorithm can be considered as an interface where the fluid and the cellular structure exchanges. The CURVIB method used here has been applied and validated in various FSI problems across different biological engineering areas.^46–49

To enable the continuum–particle coupling of the blood plasma and the cellular structures in our FSI algorithm, the flow imparted loads $({\mathbf{F}}_i^{\mathit{external}})$ must be calculated on Γ_RBC as required by Eq. (9).

F. Fluid–structure interaction simulation of RBC in flows

G. Computational setups

Fluid–structure interaction simulations are performed to determine the dynamics of RBC in a confined micro-channel.⁴ The physical properties of the RBC are described in Table I. The corresponding DPD parameters are shown for the corresponding number of vertices N_v in Table II. The computation domain is defined as a rectangular channel containing a single RBC, as illustrated in Fig. 3(a). The dimensions of the domain along the x, y, and z are L_x (the length), L_y (the width), and L_z (the height), respectively. The computational domain is discretized as a structured grid of size N_i × N_j × N_k with the spatial resolution in three directions (i, j, k) are Δx ×Δy ×Δz, respectively. The no-slip condition is applied to all channel walls. The details of the channels used in the simulations are listed in Table III.

2. Stepwise oscillatory flows (I_s)

To further validate our FSI model in oscillatory flows, the propulsion of the RBC in square channels is investigated.¹³ Two square channels (Channel-2 and Channel-3 in Table III) with side lengths L_z = 16 and 21 µm are used for the simulations, resulting in confinements χ₂ = 0.4 and χ₃ = 0.3, respectively. The initial location of the RBC is on the channel axis (x₀ = 16 µm, r = r₁ = 0). The computational configuration including the grid spacing, RBC surface meshes, and boundary conditions are shown in Fig. 3(a) and Table III. A stepwise asymmetric oscillatory waveform I_s is used with two phases: (i) forward T_f and (ii) backward T_b periods $\left(\frac{T_b}{T_f}=4\right)$, as shown in Fig. 4(a). The velocities during the forward and backward phases are ψ_f and ${\psi }_b\left(-\frac{{\psi }_f}{{\psi }_b}=4\right)$, respectively. The formula for the waveform is defined as

$$U(t)=\left\{\begin{array}{cc}{\psi }_f\hfill & \text{for}0\le t<\frac{T}{5},\hfill \\ {\psi }_b\hfill & \text{for}\frac{T}{5}\le t<T.\hfill \end{array}\right.$$(16)

Following this formula, the flow has a forward phase (ψ_f > 0) and a backward flow phase (ψ_b < 0). The maximum shear rate is defined as $\overline{\dot{{\gamma }_f}}=\frac{{\psi }_f}{L_z}$.

The Capillary number Ca^f in the forward flow phase is given as

$$C{a}^f=\frac{4{\psi }_f{L}_s{t}_R}{L_z^2},$$(17)

with $t_R=\frac{L_s{\mu }_{\mathit{plasma}}}{2{\mu }_0}$ and μ₀ is the shear elastic modulus described in Table I.

There are six values of ψ_f examined as ψ₁ = 1.05, ψ₂ = 1.58, ψ₃ = 2.1, ψ₄ = 1.9, ψ₅ = 2.8, and ψ₆ = 3.7 mm s⁻¹, respectively. Following the naming convention of the simulations, six cases are formed with the respective parameters: I_sψ₁r₁χ₂; I_sψ₂r₁χ₂; I_sψ₃r₁χ₂; and I_sψ₄r₁χ₃, I_sψ₅r₁χ₃, and I_sψ₆r₁χ₃, as shown in Table IV. As the waveform applied is of a stepwise nature, there is no relaxation time taken into consideration for these cases (t_relax = 0).

Under these oscillatory conditions, the axial propulsion step (Δx_c) is recorded at the end of the forward time interval of the asymmetric oscillating flow $(t=T_f=\frac{T}{5})$, as a function of the forward (peak) capillary number Ca^f for the chosen shear rates.¹³ Thus, Δx_c is defined as the displacement of the RBC’s centroid (C) at the end of the forward phase $(t=\frac{T}{5})$,

$$\mathrm{\Delta }{x}_c=x_c\left(t=\frac{T}{5}\right)-x_c(t=0).$$(18)

3. Sinusoidal flow simulations

To study the effect of the pulsatile flow on the propulsion and the behavior of the cellular response (morphology changes) of the RBC, we considered time-periodic flow U(t). The flow time period consists of three separate phases: (i) the forward (T_f); (ii) the resting (T_r), and the backward (T_b) periods. Following experimental protocols,¹⁴ the resting period T_r was introduced in the bulk flow waveform to allow the RBC to restore its initial equilibrium shape. The oscillatory frequency is chosen to be f = 20 Hz²⁷ and thus T = T_f + T_r + T_b = 50 ms. The asymmetry of the waveform is adjusted by changing the values of T_f, T_r, and T_b. The formula for the waveform is

$$U(t)=\left\{\begin{array}{cc}A\sin \left(2\pi \frac{t}{T_f}\right)\hfill & \text{for}0\le t\le T_f,\hfill \\ 0\hfill & \text{for}{T}_f\le t\le T_f+T_r,\hfill \\ A\sin \left(2\pi \frac{(t-T_f-T_r)}{T_b}\right)\hfill & \text{for}{T}_f+T_r\le t\le T.\hfill \end{array}\right.$$(19)

The symmetric waveform (I₁) is created with T_f = T_b (completely symmetric). The asymmetric waveforms (I₂, I₃, and I₄) are formed by progressively reducing the period of T_b. Four distinct inflow types were generated with symmetry and asymmetric waveforms (I₁, I₂, I₃, and I₄), as seen in Fig. 5 and Tables V and VI. For each of these waveforms, three different velocity magnitudes (A = U₁, U₂, and U₃) and three different radial shifts (r₁, r₂, and r₃) are considered, as shown in Table VII. In total, the combinatoric arrangements lead to a total of 36 distinct simulation cases with the notation I_mU_nr_pχ₁ with the corresponding values of the indices m = 1, 2, 3, 4; n = 1, 2, 3; and p = 1, 2, 3. The outline of the simulation cases is shown in Table VII. In addition, the RBC shapes are recorded over a time period of two cycles 2T as exemplified in Fig. 5(b), in which the initial location of the RBC is set at x₀ = 22.5 µm. Due to the nature of the sinusoidal waveform applied, there was no relaxation time for all of these cases (t_relax = 0). The centroid’s displacement is monitored continuously as the function of time,

$$\mathrm{\Delta }{x}_c(t)=x_c(t)-x_c(t=0).$$(20)

Time
Waveforms	Tf (ms)	Tr (ms)	Tb (ms)
I 1	20	10	20
I 2	25	12.5	12.5
I 3	27.3	9	13.7
I 4	28.6	7.1	14.3

Subscript	Inflow waveform	ψf (mm s−1)	Radial shift r (μm)
1	I 1	U1 = 1	r1 = 0
2	I 2	U2 = 1.5	r2 = 0.4
3	I 3	U3 = 2	r3 = 0.7
4	I 4	U4 = 6	⋯

(a)		r 1	r 2	r 3
I 1	U 1	I 1 U 1 r 1 χ 1	I 1 U 1 r 2 χ 1	I 1 U 1 r 3 χ 1
	U 2	I 1 U 2 r 1 χ 1	I 1 U 2 r 2 χ 1	I 1 U 2 r 3 χ 1
	U 3	I 1 U 3 r 1 χ 1	I 1 U 3 r 2 χ 1	I 1 U 3 r 3 χ 1

(b)		r 1	r 2	r 3
I 2	U 1	I 2 U 1 r 1 χ 1	I 2 U 1 r 2 χ 1	I 2 U 1 r 3 χ 1
	U 2	I 2 U 2 r 1 χ 1	I 2 U 2 r 2 χ 1	I 2 U 2 r 3 χ 1
	U 3	I 2 U 3 r 1 χ 1	I 2 U 3 r 2 χ 1	I 2 U 3 r 3 χ 1

(c)		r 1	r 2	r 3
I 3	U 1	I 3 U 1 r 1 χ 1	I 3 U 1 r 2 χ 1	I 3 U 1 r 3 χ 1
	U 2	I 3 U 2 r 1 χ 1	I 3 U 2 r 2 χ 1	I 3 U 2 r 3 χ 1
	U 3	I 3 U 3 r 1 χ 1	I 3 U 3 r 2 χ 1	I 3 U 3 r 3 χ 1

(d)		r 1	r 2	r 3
I 4	U 1	I 4 U 1 r 1 χ 1	I 4 U 1 r 2 χ 1	I 4 U 1 r 3 χ 1
	U 2	I 4 U 2 r 1 χ 1	I 4 U 2 r 2 χ 1	I 4 U 2 r 3 χ 1
	U 3	I 4 U 3 r 1 χ 1	I 4 U 3 r 2 χ 1	I 4 U 3 r 3 χ 1

A. Coarse-graining validation

A stretching test is carried out and aimed to replicate the experimental test of Ref. 51. In our simulations, the parameters used to describe the physical characteristics of the RBC are listed in Table I. Following the coarse-graining procedure, the model parameters for the cell membrane such as the equilibrium length, the persistence length, the spring stiffness, and the spontaneous angle are computed for each value of N_v, as in Table II. The cytosol fluid is modeled by a set of particles N_f = 100, placed within the interior volume of the cell membrane, as shown in Fig. 6(b).

In this experiment, two external forces F_stretch and −F_stretch with opposite directions are applied on both sides of the RBC. The magnitude of the force F_stretch is increased in a stepwise manner from 0 to 200 pN (a total of 16 steps). The axial diameter (D_a) and transverse diameter (Dt) are measured for every step. The definitions of D_a and D_t are shown in Fig. 6(a). D_a is calculated as |x_max − x_min| representing the axial distance between the maximum x position, denoted as x_max, and the minimum x position, denoted as x_min of the membrane vertices. D_t is calculated using the formula $2\times {\max }_{i=1\dots N_v}\sqrt{{(y_i-c_y)}^2+{(z_i-c_z)}^2}$, where c_y and c_z represent the y and z coordinates of the RBC’s center of mass, respectively. The simulations are performed systematically with different RBC surface mesh resolutions by changing the number of vertices (N_v) (also see Table II).

The current RBC model accurately replicates the elastic response of the RBC under stretching forces, as revealed by the results of D_a and D_t in Fig. 6. During stretching, the dynamic response of cytosol particles is visible, indicating the coupling between the membrane and the cytosol fluid. The shapes of the RBC under loading conditions agree with ones from experimental data of Ref. 51. The computed values of the axial (D_a) and transverse (D_t) diameters agree well with the experimental values, as seen in Fig. 6(a). In particular, the values of D_a and D_t are consistent across the different values of N_v, which indicates a robust performance of the coarse-graining procedure. There is a disagreement between the simulated results and the experimental value of D_t. Examining the shapes of the RBC in the simulations [Fig. 6(b)], it is revealed that the RBC tends to rotate around the stretching direction. This rotation leads to the difference between the experimental and numerical results of D_t. In brief, the mechanics of RBC is well-replicated by the computational model across different levels of coarse-graining. Thus, the value of N_v = 1000 is chosen to report the dynamics of the RBC.

B. Model validation

1. Deformation of the RBC under a constant shear rates $\overline{\dot{{\gamma }_0}}$

Under constant shear rate conditions (I₀), as described in Sec. II G 1, two districts of the RBC shape are observed: (i) the croissant shape (I₀U₃r₁χ₁‐$\overline{\dot{{\gamma }_0}}=200s^{-1}$) and (ii) the slipper shape (I₀U₄r₃χ₁—$\overline{\dot{{\gamma }_0}}=600s^{-1}$), as shown in Fig. 7.

Under low shear rate (I₀U₃r₁χ₁), the RBC is initially placed along the centerline of the microchannel (discocyte shape). As the RBC interacts with the incoming flow, it deforms and eventually transitions to a croissant shape, as shown in Fig. 7(a). The terminal shape (croissant) is attained as the RBC continues to propagate along the channel’s symmetry axis. Note that the croissant shape, in this case, is not fully axis-symmetric. In addition, the croissant shape observed in this study is comparable with the parachute shape from our previous work,²⁹ although not identical due to differences in channel geometry and flow conditions.

Under high shear rate (I₀U₄r₃χ₁), the RBC transitions to the slipper shape, as shown in Fig. 7(b), which exhibits a bistability mode with tank-treading behavior. Note that the RBC is placed at a radial shift r₃ = 0.7 µm. Thus, the initial location of the RBC is not at the channel’s symmetry axis. The tank-treading effect is a complex dynamic in which the RBC membrane propagates axially along the channel while it rotates around its own center of mass, as seen in Fig. 7(b). A counterclockwise rotation is observed as indicated by the locations of two membrane particles (Lagrangian points—V₁ and V₂) at different time instances (t₁ = 22 ms and t₂ = 25 ms).

In both the croissant or slipper shapes, the transition from the initial shape (discocyte) to the terminal shape (either croissant or slipper) occurs within around 30 ms, as seen in Fig. 7. Our simulation results agree well with the corresponding experimental data of Ref. 4 as well as those described in recent experiments on RBC transient dynamics.^52,53 Furthermore, these transitions are in good agreement with the shape diagram⁵⁴ for different capillary numbers and confinements. In conclusion, our simulations are able to replicate the dynamics of the croissant and slipper shapes excellently.

The extracellular patterns of the croissant and slipper shapes also agree well with the experimental data.^4,26 The extracellular flow pattern can be visualized by reconstructing the relative flow velocity field²⁶ (co-moving frame). The relative velocity is defined as the difference between the flow velocity and the RBC’s centroid velocity, as shown in Fig. 8. In the croissant shape (I₀U₃r₁χ₁), the velocity streamlines closely resemble an axi-symmetrical flow pattern [Fig. 8(a)]. The downstream side of the RBC membrane deforms significantly, whereas the upstream side barely changes, as depicted in Fig. 8(b). In the slipper shape (I₀U₄r₃χ₁), there exists an asymmetrical vortical structure in the vicinity of the RBC membrane. As the slipper shape emerges, a fully closed vortex ring is created by a separated flow region, which is close to the channel wall. In short, the emergence of the RBC shape dictates the extracellular flow pattern.

2. Propulsion of RBC under stepwise oscillatory flows (I_s)

Under stepwise flow waveform (I_s), our simulation results agree well with the propulsion step map (Δx_c, Ca^f), which was developed¹³ for both channels χ₂ = 0.5 and χ₃ = 0.38. In both cases, the propulsion step (Δx_c) is observed to monotonically increase with the values of Ca^f. However, Δx_c is higher in the lower confinement channel (χ₃), which indicates the importance of channel confinement.

In all simulation cases (I_sψ₁r₁χ₃, I_sψ₂r₁χ₃, I_sψ₃r₁χ₃, I_sψ₄r₁χ₃, I_sψ₅r₁χ₃, and I_sψ₆r₁χ₃), the RBC transitions from the idealized discocyte to the biconcave shape during the forward phase $(0<t<\frac{T}{5})$ with all values of the peak forward flow (ψ_f = 1.05 mm s⁻¹ to ψ_f = 4.34 mm s⁻¹), as shown in Fig. 7(c). Strikingly, the complex multilobe shape emerges during the backward phase T_b. The elastic response of the RBC membrane to the oscillatory flow during the cycle T is depicted for the case I_sψ₆r₁χ₃ in Fig. 7(c). The reversal of the flow direction during T_b results in membrane buckling and stretching, which give rise to the multilobe shape. Note that these cases are initially placed at the channel center (r = r₁ = 0). Therefore, it is possible to induce complex dynamics of RBCs in a confined channel by changing the inflow waveform only.

C. The impact of sinusoidal flows on RBC dynamics

1. The emergence of RBC shapes

The sinusoidal flow waveform (I₁, I₂, I₃, and I₄) further adds complexity to the membrane dynamics as the shape of RBC is highly sensitive to the extracellular flow condition. As a result of the pulsatile flow condition, the RBC shape continuously responds to the applied flow in the channel. Our simulations show that the RBC alternates its shapes in one of the following types: (1) croissant; (2) slipper; (3) trilobes; (4) simple/complex/elongated multilobes; (5) rolling stomatocytes; (6) hexalobes; and (7) rolling discocyte, as shown in Fig. 9 and Tables VIII–XI. The emergence of each type will be discussed in the following.

Waveform (I2)	C	S		CM	EM		RD
	r 1	r 2	r 3	r 1	r 2	r 3	r 1	r 2	r 3
U 1	0.27	0.29	0.29	0.9	⋯	⋯	1.55	1.33	1.33
U 2	0.2	0.22	0.22	1.06	⋯	⋯	1.44	1.3	1.3
U 3	0.2	0.22	0.22	1.06	1.16	1.16	1.8	1.3	1.3

Waveform (I4)	C	S		CM	EM		RD			HX
	r 1	r 2	r 3	r 1	r 2	r 3	r 1	r 2	r 3	r 1
U 1	0.32	0.31	0.31	0.6	⋯	⋯	1.38	1.27	1.27	⋯
U 2	0.27	0.24	0.24	0.56	⋯	⋯	1.5	1.24	1.24	1.15
U 3	0.27	0.24	0.24	0.56	1.15	1.15	1.1	1.24	1.24	⋯

In all cases, the RBC evolves to the croissant (C) or the slipper (S) shape during the forward phase (0 < t < T_f) of the flow cycle $(\frac{t}{T}\approx 0.25)$ , as shown in Tables VIII–XI. This process is demonstrated for the cases I₁U₁r₁χ₁, I₃U₂r₁χ₁, I₄U₂r₁χ₁, and I₄U₃r₃χ₁, as seen in Fig. 9 (first column). Note that the transition to C or S mode from the biconcave shape is dependent on the value of the radial shift (r). As shown in Tables VIII–XI, the S mode appears only when the RBC is initially placed not exactly at the cross-sectional center (r > 0) (for example, I₄U₃r₃χ₁). The RBC remains in C mode during the forward phase if it is initially placed at the cross section center (r = 0) regardless of the bulk flow waveform I₁, I₂, I₃, or I₄ as seen in cases I₁U₁r₁χ₁, I₃U₂r₁χ₁, and I₄U₂r₁χ₁. In brief, the croissant and the slipper shapes exist during the forward phase and their emergence depends on the radial shift (r).

The RBC transitions from the simple shapes (croissant and slipper) toward more complex shapes (trilobes, simple/complex/elongated multilobes, rolling stomatocytes, hexalobes, and rolling discocyte) later in the flow cycle during the resting/reverse periods $(\frac{t}{T}>0.5)$ , as shown in Tables VIII–XI and Fig. 9 (middle column). The shape transformation is initiated by the buckling of the RBC membrane, which takes place in the resting interval (T_r) (see Fig. 5) of the flow. As a result of the change in flow direction in T_b, the RBC experiences considerable stretching and compression, leading to significant alterations in its membrane shape, as seen in Fig. 9. Note that the RBC does not return to either C or S mode at the beginning of the second cycle T < t < 2T, as shown in Tables VIII–XI. Thus, the shape transition process is irreversible even if the inflow waveform is symmetric (I₁), as shown in Table VIII.

	C	S		CM	M	T	RS	EM	RD
Waveform (I1)	r 1	r 2	r 3	r 1	r 1	r 1	r 1	r 1	r 1	r 2	r 3
U 1	0.32	0.3	0.3	⋯	0.75	⋯	0.9	⋯	1.8	1.25	1.25
U 2	0.21	0.2	0.2	0.5	⋯	0.8	⋯	1.2	⋯	1.21	1.21
U 3	0.21	0.2	0.2	0.5	⋯	0.8	⋯	1.2	⋯	1.21	1.21

2. The impacts of the initial position (the radial shift—r)

Our results in Fig. 9 show that the initial position (r) plays a critical role in the emergence of RBC shapes. The RBC is placed initially at the channel axis r = r₁ = 0 under both the symmetric (I₁U₁r₁χ₁) and asymmetric (I₃U₂r₁χ₁ and I₄U₂r₁χ₁) waveforms shown in Figs. 9(a)–9(c), respectively. In these cases (I₁U₁r₁χ₁, I₃U₂r₁χ₁, and I₄U₂r₁χ₁), the RBCs all transition sequentially from the croissant shape toward the multilobes/complex multilobes and finally one of the poly-lobes (trilobes, rolling stomatocytes, or hexalobes). When r = r₃ > 0 in the case I₄U₃r₃χ₁, the RBC mostly remains in the slipper shape during the forward phase (t < T_f), as seen in Fig. 9(d). It transitions toward the elongated multilobe during the backflow phase (T_f + T_r < t < T). Finally, the RBC becomes a rolling discocyte in the second cycle (t ≈ 1.2T). Note that r₃ is a small value (0.7 µm) and thus the shape transition process is strongly sensitive to the initial placement of the RBC. The dependence of RBC dynamics on the value of r is explained in the following.

When the RBC is positioned at the centerline (r = 0) of the channel, its migration is highly dependent on the type of inflow waveform. When the RBC is subjected to a symmetric waveform (I₁) as shown in Fig. 10(a), it oscillates around its initial position with a minimal migration along the channel (axial) direction. However, as the inflow profile transitions to an asymmetric waveform (I₂, I₃, and I₄) with an increasing T_f, the RBC gains more momentum during the forward phase and propels far away from its initial position, as shown in Fig. 10(a) (left). The value of the axial displacement (Δx_c) (t) increases sequentially from I₁, I₂, I₃ to I₄. Δx_c reaches a maximum value of ∼4 × L_s at the end of the second cycle for the case I₄. In the lateral direction, Δy_c reaches a value of ∼0.16 × L_s at the end of the first cycle for the symmetric waveform (I₁), as shown in Fig. 10(b) (left). For the cases I₂, I₃, and I₄, the values of Δy_c remain comparably low $(\approx 0.08\times L_s)$ during the cycles. Furthermore, the vertical displacement Δz_c seen in Fig. 10(c) follows a monotonically upward trend throughout the second cycle. This results in a maximum vertical displacement Δz_c of ∼0.25 × L_s for I₁. For other waveforms (I₂, I₃, and I₄), a smaller upward trend is observed with a vertical displacement Δz_c ≈ 0.08 × L_s at the end of the first cycle. During the second cycle, the cell is observed to migrate downward. In summary, the symmetric waveform (I₁) leads to a maximum displacement in the transverse directions while the asymmetric waveforms (I₂, I₃, and I₄) result in the maximum axial displacement for r = 0.

As the RBC is placed at the off-centered location (r = r₂ = 0.4 µm), as seen in Fig. 10(d), its axial migration is similar to ones with r = 0 [Fig. 10(a)]. The axial displacement (Δx_c) is dependent on the applied inflow waveform and also reaches 4 × L_s for the I₄ waveform. These values indicate that the radial shift (r) does not significantly affect the axial migration of the RBC. However, the lateral migration of RBC in Figs. 10(e) and 10(f) does show a dependence of the applied waveform. The values of Δy_c and Δz_c reach 0.14 × L_s for I₁ and I₂ at the end of the second cycle, respectively. Surprisingly, the impact of the applied waveform seems to diminish as the value of T_f increases. Comparing the case I₃U₂r₂χ₁ and I₄U₂r₂χ₁ (I₃ and I₄), the values of Δy_c and Δz_c shown in Figs. 10(e) and 10(f) are nearly identical, resulting in a vertical displacement of ∼0.04 × L_s. To summarize, the impact of the applied waveform is significant on the axial displacement of the RBC but it is milder on the lateral displacements when r > 0. The impact of the applied waveform diminishes as the forward time T_f increases.

3. The impacts of the waveform (I)

It is striking to observe the irreversible dynamics of RBC in Fig. 11. When the RBC is subjected to the symmetric waveform (I₁) with different radial shifts r = r₁, r₂, and r₃, the RBC oscillates around its initial position with a minimal displacement along the axial direction as shown in the cases I₁U₁r₁χ₁, I₁U₁r₂χ₁, and I₁U₁r₃χ₁. Despite the inflow waveform being completely symmetrical (a sine function—I₁), the axial position of the RBC in Fig. 11(a) (left column) shows a net positive value of the displacement Δx_c at the end of the first (t = T) and second cycle (t = 2T) even when there is no radial shift (r = 0). Although small, this positive value of Δx_c indicates that the RBC does not go back exactly to its initial location, which is Δx_c(t = 0) = 0. At all values of the radial shift of r = 0, 0.4, and 0.7 µm, this irreversible dynamics is even more evident, as shown in the lateral displacements in Figs. 11(b) and 11(c). The magnitudes of Δy_c and Δz_c are comparable for all values of r during the cycles. For the case I₁U₁r₁χ₁ (r = 0), the value of Δy_c reaches a value of ∼0.16 × L_s at the end of the first cycle. For the cases I₁U₁r₂χ₁ and I₁U₁r₃χ₁, the values of Δy_c and Δz_c reach ∼0.25 × L_s at the end of the second cycle. In the vertical direction (z_c) shown in Fig. 11(c), the well-centered RBC (r = 0) is influenced by the change of flow direction, which is depicted by the upward and downward trends in the first cycle. However, the cell follows a dominant upward trend during the entire second cycle, resulting in a lateral migration of around 0.25 × L_s. Therefore, there exists a significant lateral migration of the RBC during its propagation regardless of its initial position under the symmetrical waveform condition (I₁). In conclusion, a symmetrical flow waveform (I₁) results in minimal propulsion along the axial direction but a significant lateral migration.

Under asymmetric waveform I₄, the RBC propels along the channel direction with an axial displacement of ∼2 × L_s in each cycle, as shown in Fig. 11(d). As the waveform becomes asymmetric with a longer forward phase, the RBC does not go back significantly during the reverse phase. It rather remains at a displacement value of Δx_c ≈ 1.9 × L_s at the end of the first cycle. It continues to propel in the second cycle up to Δx_c ≈ 4.0 × L_s. Surprisingly, the lateral displacements of the RBC (Δy_c, Δz_c) are smaller in comparison with ones in the symmetric case (I₁). The values of (Δy_c, Δz_c) are within 0.15 × L_s for all cases I₄U₁r₁χ₁, I₄U₁r₂χ₁, and I₄U₁r₃χ₁, as shown in Figs. 11(e) and 11(f). In brief, the RBC propels significantly under the impact of the asymmetrical flow waveform I₄ along the axial direction but it does not migrate significantly in the lateral directions.

4. Extracellular flow dynamics at the vicinity of the RBC under oscillatory flows

The emergence of the RBC shape has a close relationship with the flow pattern of the surrounding fluid (extracellular flow). Under the impact of the channel confinement, the deformation of RBC is well regulated by the flow waveform, which results in distinct extracellular flow patterns surrounding the RBC, as shown in Figs. 9 and 12. To highlight the impact of the RBC motion, the flow pattern is visualized in the co-moving frame with the RBC’s centroid (see Sec. III B 1). Thus, the flow streamlines are represented from the perspective of the RBC.

The case I₁U₁r₁χ₁ is selected to illustrate the evolution of flow pattern as the RBC deforms from a relatively simple shape to a more complicated shape, as depicted in Figs. 12(a) and 12(b). This case is chosen because the temporal variation of the waveform is completely symmetrical (I₁). Moreover, the RBC is placed initially at the channel axis (r = r₁ = 0) with the lowest forward velocity ψ_f = U₁ = 1 mm s⁻¹. In the case I₁U₁r₁χ₁, Fig. 12(a) revealed that the RBC has a multilobes shape at the end of the forward phase. The presence of the large lobes results in a more convoluted streamlines pattern during the resting phase. As the RBC undergoes a morphological transition to rolling stomatocytes at the end of the first cycle (t = 0.9T), the streamlines exhibit changes [Fig. 12(b)]. During the second cycle, the RBC gradually transforms into a rolling discocyte by the end of the second cycle (t = 1.8T).

The impact of the radial shift (r) on the RBC shape and the resulting flow pattern is significant. To highlight the impact of the initial location, the case I₁U₃r₃χ₁ was selected to visualize the flow patterns. As shown in Fig. 12(c), due to the off-centered initial location (r > 0), the slipper shape emerges during the forward phase. A closed vortex ring is also observed downstream of the RBC as the flow velocity reaches its maximum magnitude in the forward phase. This phenomenon is similar to the one observed in the constant shear rate case (I₀U₄r₃χ₁ with U₄ = 6 mm s⁻¹) in Fig. 8(b). This is remarkable since the peak flow ψ_f is rather three times lower in this case ψ_f = U₃ = 2 mm/s.

The case I₁U₃r₁χ₁ [Figs. 12(d)–12(g)] is selected to illustrate further the impact of the peak forward flow ψ_f. In this case, the peak velocity is ψ_f = U₃ = 2 mm s⁻¹. In comparison with the case I₁U₁r₁χ₁ [Figs. 12(a) and 12(b)], only the value of the peak velocity ψ_f is increased. However, the RBC shape evolves in a completely different sequence, as opposed to the case I₁U₁r₁χ₁. During the forward phase, the RBC transitions quickly to the croissant shape in Fig. 12(d) at (t = 0.28T), just after the peak forward flow. The flow patterns are similar to those observed under case I₀U₃r₁χ₁ in Fig. 8(a). During the resting period (T_f < T < T_f + T_r), the flow velocity surrounding the cell decreased notably and the complex multilobes shape emerges, as seen in Fig. 12(e). The flow pattern is perturbed minimally surrounding the RBC as its shape turns to the trilobes shape, as depicted in Fig. 12(f). During the forward flow phase of the second cycle, the elongated multilobes appear (t = 1.2T).

The overall dynamics of the RBC over the cycles depend significantly on the radial shift. As demonstrated in Tables VIII–XI, the multilobes shape appears at the beginning of each cycle if the RBC is located initially at the channel center r = 0. Under the specific condition of the case I₄U₂r₁χ₁, the multilobes shape can further transform into the hexalobes shape during the forward flow phase of the second cycle (t = 1.15T), as shown in Table XI. Its corresponding flow patterns are shown in Fig. 12(h), in which the extracellular flow was observed to exhibit a minimal disturbance around the hexalobes shape as the RBC. For the majority of the off-centered cases (r > 0), the rolling discocyte case is found the most commonly, as shown in Tables VIII–XI. The flow pattern around a discocyte is exemplified in Fig. 12(i) for the case I₁U₁r₂χ₁ at t = 1.25T. Here, the flow streamlines show a distinct separation of upstream and downstream regions.

A. The emergence of the croissant shape and the slipper shape under a constant shear rate $\overline{\dot{{\gamma }_0}}$)

Recent studies^4,26 in rectangular microchannels, which are identical to our channels as shown in Fig. 3 and Table III, suggest that the emergence of the RBC shape is dependent on the radial shift (r—see Fig. 3 for its definition). In previous studies,^4,26 the croissant shape emerged at $\overline{\dot{{\gamma }_0}}<300s^{-1}$ if the RBC is placed exactly at the channel’s center (r = 0). On the other hand, the slipper shape emerged whenever the RBC was not placed exactly at the centerline (r > 0). The RBC was found to exhibit a (tank-treading) slipper shape at sufficiently high shear rate $(\overline{\dot{{\gamma }_0}}\approx 500s^{-1})$ and off-centered placement (r > 0).^4,26 In cylindrical microchannels,¹⁷ similar observations were confirmed, albeit at lower shear rates $(0<\overline{\dot{{\gamma }_0}}<80s^{-1})$. These studies point to the importance of the radial shift in regulating RBC dynamics. Our results shown in Fig. 7 agree with these observations (also see Tables III and IV). Our results show the appearance of croissant-to-slipper transition as the capillary number (and thus $\overline{\dot{{\gamma }_0}}$) increases from 0.1 to 0.37 for a confinement of χ = 0.65. The croissant shape emerges when the initial position of the RBC is placed exactly at the channel centerline at a sufficiently low shear rate (Ca = 0.1). When the shear rate is increased to Ca = 0.37, the slipper shape emerges. Furthermore, our model is able to capture the intricate dynamics of the tank-treading motion, which is characterized by the rotation of the membrane at the shear rate of 600 s⁻¹, as illustrated in Fig. 7. Therefore, our results further confirm the importance of the radial shift in the croissant-to-slipper transition.

B. The impact of time-varying shear rate $\overline{\dot{\gamma (t)}}$ on RBC shape

When the inflow varies in a stepwise manner, as seen in Fig. 4, the shear rate $(\overline{\dot{\gamma }})$ changes as a function of time $\overline{\dot{\gamma }}(t)$ with distinct forward (T_f) and backward (T_b) time phases. In all cases (I_sψ₁r₁χ₂, I_sψ₂r₁χ₂, I_sψ₃r₁χ₂, I_sψ₄r₁χ₃, I_sψ₅r₁χ₃, and I_sψ₆r₁χ₃), the RBC is placed exactly at the channel axis (r = r₁ = 0). The shape transitions are accomplished through consistent transient stretching and compression of the membrane. This occurs as the RBC experiences forward and backward flow phases during the flow cycles. The RBC transitions from a discocyte shape toward a croissant shape during its forward propulsion, as shown in Fig. 7. Although the backward phase induces the buckling of the cellular membrane, the RBC shape remains relatively symmetrical with respect to the channel axis (multilobes) as shown in Fig. 7 at the end of T_b. This is remarkable given that the maximum shear rate during the backward phase can be sufficiently large $({\overline{\dot{\gamma }}}_f=200s^{-1})$ but this symmetry is still maintained. Comparing the case I₀U₃r₁χ₁ and I₀U₄r₃χ₁ in Table IV, our results suggest that the complete break of symmetry⁵² (slipper shape) is observed only when the radial shift exists (r > 0).

When applying different sinusoidal waveforms (I₁, I₂, I₃, and I₄) shown in Fig. 5, our results show the ubiquitous presence of croissant and slipper shapes across all shear rates (${\overline{\dot{\gamma }}}_f=100,150$, and 200 s⁻¹). The slipper shape appears at (t ≈ 0.3T) whenever the RBC is placed off the channel’s axis (r > 0), as shown in Tables VIII–XI. Note that these waveforms are different in terms of the forward (T_f) and backward (T_b) phases, with the backward phase being the shortest in I₄. As shown in Table VIII, the slipper shape is observed even when the waveform is completely symmetric (I₁) given that r > 0, as in (a) I₁U₁r₂χ₁, (b) I₁U₁r₃χ₁, (c) I₁U₂r₂χ₁, (d) I₁U₂r₃χ₁, (e) I₁U₃r₂χ₁, and (f) I₁U₃r₃χ₁. Hence, our results indicate that the flow waveform does not affect the emergence of the RBC shape during the forward phase. Instead, the radial shift plays an essential role in this process.

It has been demonstrated¹¹ that the presence of the discocyte shape is correlated with weak shear rates $\overline{{\gamma }_0}$. Under ${\overline{\dot{\gamma }}}_0<10s^{-1}$, the RBC typically maintains its discocyte shape with an 80% probability.¹¹ However, as the shear rate gradually rises $(10s^{-1}\le {\overline{\dot{\gamma }}}_0\le 400s^{-1})$, the likelihood of a discocyte shape decreases to 30%. This observation has been validated even when different viscosity ratios (λ)¹² are considered. In our work, the discocyte shape is ubiquitously observed during the second flow cycle across all applied waveforms (I₁, I₂, I₃, and I₄) and shear rates (${\overline{\dot{\gamma }}}_f=100$, 150, and 200 s⁻¹) with (r > 0) as shown in Tables VIII–XI. Therefore, our results further confirm that discocyte is the most common shape with the range shear rate less than 200 s⁻¹.

In our study, complex shapes evolve from simpler shapes under the influence of time-dependent waveforms. In the experimental work¹¹ under constant shear rates, stomatocytes shape was observed to dominate the RBC population (65%) when the shear rate was between $10s^{-1}<\overline{\dot{{\gamma }_0}}<400s^{-1}$. In our study, the elliptical-rim-shaped stomatocytes only emerge under symmetric waveform (I₁) in the case I₁U₁r₁χ₁ under the peak shear rate of ${\overline{\dot{\gamma }}}_f=100s^{-1}$. At high constant shear rates $400s^{-1}<\overline{\dot{{\gamma }_0}}<2000s^{-1}$, experimental data¹¹ showed that polylobes shape could emerge. This polylobes shape is characterized by a large number of lobes on the RBC surface, known as multilobes, trilobes, and hexalobes. In the current study, polylobes are also observed across all applied waveforms, given that the cell is placed initially at the channel axis (r = r₁ = 0) even at weak shear rates $({\overline{\dot{\gamma }}}_f\le 200s^{-1})$ as in Fig. 9 and Tables VIII–XI. For example, the multilobes are observed with all waveforms. The trilobes are both observed in the symmetric waveform (I₁U₂r₁χ₁ and I₁U₃r₁χ₁—Table VIII) or the asymmetric waveform (I₃U₂r₁χ₁—Table X). The hexalobes shape only appears under the most asymmetric waveform with $(I_4{U}_2{r}_1{\chi }_1)-\overline{\dot{{\gamma }_f}}=150s^{-1}$, as shown in Table XI. In the previous work,⁵⁵ the RBC shape has been reported to deform further into elongated shapes as the shear rate increases. Our results support this observation shown in Fig. 12(g) as the elongated multilobes appear during the backward flow phase. Examining Tables VIII–XI, our results suggest that this elongated shape is generally present regardless of the applied waveform but it only manifests under higher shear rates of at least ${\overline{\dot{\gamma }}}_f=150s^{-1}$. In particular, the elongated multilobes are observed under I₁U₂r₁χ₁, I₁U₃r₁χ₁, I₂U₃r₂χ₁, I₂U₃r₃χ₁, I₃U₃r₂χ₁, and I₃U₃r₃χ₁. In essence, our results strongly suggest that the complex shape can appear in a micro-channel even at weak shear rates given that the inflow is time-dependent.

While the radial shift is important for the emergence of the RBC shape, it is less important in controlling the cell displacement (the center of mass), as shown in Fig. 11. The background flow velocities vary largely over the cycle in the forward and reverse phases, depending on the flow waveform, as shown in Fig. 5. This variation dominates the plasma flow surrounding the RBC. Therefore, small radial shifts (0–0.7 µm) play a minimal role in regulating the overall displacement of the RBC’s center of mass.

C. Controlling lateral migration of cells with oscillatory flows

Microfluidic devices are typically used to isolate and separate cells.⁵⁶ While these devices are promising for many cell-sorting applications,^57,58 the main challenge is the difficulty in obtaining high-throughputs due to the required length of the microfluidic channels. Changing the geometrical design of channels⁵⁹ has been proposed as one efficient way. Recent studies have shown that varying the shear rates in time^13,14 can reduce the required length based on the concept of velocity lift,⁶⁰ which is the factor that drives the RBC’s migration toward the center of the channel.

As the inertial effect is negligible at a very low Reynolds number (R_e ≈ 0.01), the flow is reversible for a rigid body. Thus, a rigid body will return to its initial position if the inflow conditions in the backward phase are reversed in the exact opposite way of its own during the forward phase. However, the RBC is not a rigid body and its membrane is highly flexible. Experiments¹⁴ have shown that the responses of RBCs are different from the ones of stiff beads in oscillatory flows. The RBC migration is observed to have a net actuation in oscillating flows whereas a stiff bead does not.

Our results in Figs. 10 and 11 for the symmetric waveform (I₁) show that the RBC does not go back to its initial position at the end of the cycle. There is an axial shift of the RBC from its original position (Δx_c ≠ 0) at the end of the cycle. Moreover, the RBC migrates significantly in the lateral cross section (Δy_c ≫ 0 and (Δz_c ≫ 0), as shown in Fig. 10. Our results in Figs. 11(b) and 11(e) indicate that the RBC undergoes lateral migration. This migration is particularly evident in the cases I₁U₁r₃χ₁ and I₄U₁r₃χ₁, where the cell migrates in the lateral direction for a distance of ∼0.25L_s and 0.13L_s, respectively. Comparing the shapes of the waveform I₁, I₂, I₃, and I₄ in Fig. 5, our results suggest that significant lateral migration can be induced by adjusting the forward time interval (T_f). With our sinusoidal waveforms shown in Fig. 5, T_f should be less than three times the backward flow phase T_b ($\frac{T_f}{T_b}<3$) to be able to induce significant lateral migration. This lateral migration in oscillatory flows might be used to separate cells selectively based on their mechanical properties.

Our findings in Figs. 8 and 12 show that the extracellular flow patterns are directly influenced by the dynamics of the RBC. Under a stationary condition in Fig. 8, the extracellular flow dynamics in the croissant shape are distinctively different from the ones of the slipper shape. In particular, the flow around the steady croissant shape was similar to that of a rigid sphere,⁶¹ in which the flow streamlines move nearly symmetrically inward and outward from the cell in the upstream and downstream sides, respectively. In contrast, a fully closed vortex ring (“bolus”) was observed downstream the cell for the slipper shape.⁴ Our results in Fig. 12 suggest that the extracellular flow patterns are far more complicated and highly dependent on the type of inflow waveform in time-dependent shear rates. It is important to note that the extracellular flow has been found to play an important role in drug delivery strategies²⁶ due to its potential use of particle trapping. Therefore, our results suggest that controlling the inflow waveform either by adjusting the peak flow ψ_f or the shape of the waveform (specifically T_f) might lead to the desired effects in delivering small particles (e.g., therapeutic nano-particles) to the RBCs.

D. Limitations

While the current model has successfully simulated various RBC morphologies under different flow conditions, it is crucial to acknowledge its limitations. First, the current model assumes that the impact of the collision force of the cytosol particles [Eq. (8)] is negligible in determining the shape of the cellular membrane, as shown in many validations of Sec. III A with high shear rates. However, it is important to recognize that this collision force might have impacts on the dynamics of the membrane and cytosol under very low shear rates. Therefore, future work will need to develop new algorithms for the collision scheme.^30,40 Second, analyzing the distribution of cytosol fluid particles and their conformational temperature could offer valuable insights into cytosol dynamics. In this work, the number of cytosol particles (N_f) is set to be 100, which is much less than the number of membrane particles (N_v). The coarse resolution of the cytosol fluid may lead to significant under-representation of hydrodynamic interactions and spatial gradients in the cytoplasm. Therefore, the intracellular flow cannot be reliably reproduced. In future studies, a larger number of N_f should be implemented to capture the intracellular flow patterns accurately. Third, previous work⁶² has shown that Eq. (A11) does not satisfy angular momentum conservation. This violation might induce inaccurate behaviors in simulating the dynamics of fluids such as the cytosol fluid. This issue needs to be investigated in future work. Fourth, the intracellular pressure coefficient p_bs = 4.12 Pa is chosen empirically, which is based on the range of RBC intracellular pressure.⁶³ Further investigations might be needed to identify the optimal value of p_bs. In addition, clustering of red blood cells is an important phenomenon, which determines the overall rheology of blood flow. Clusters of RBCs are formed by the adherence of individual cells, which have been shown to be dependent not only on the characteristics of the RBC membrane but also on the hydrodynamic condition of the plasma solution.⁶⁴ While this clustering process is beyond the scope of this work, our methodology can be extended to study this phenomenon by incorporating the interaction forces between adjacent cell membranes. This issue will be addressed in our future work. Finally, our coupling methodology utilizes the sharp-interface immersed boundary method,⁴⁵ which requires only the velocities on the cellular membrane as boundary conditions. In this work, the physical time step is chosen in the order of microseconds. Therefore, the thermal fluctuations of the velocity fields on the cellular membrane are fully resolved via the implicit momentum solver of the fractional step method (Sec. II E). Future studies need to be carried out to accelerate the convergence of the fractional step method for the Navier–Stokes equations and allow larger physical time steps.

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Lahcen Akerkouch: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Trung Bao Le: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Methodology (equal); Resources (equal); Software (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).