Shape Transitions of Red Blood Cell under Oscillatory Flows in Microchannels

This paper aims to examine the ability to control Red Blood Cell (RBCs) dynamics and the associated extracellular flow patterns in microfluidic channels via oscillatory flows. Our computational approach employs a hybrid continuum-particle coupling, in which the cell membrane and cytosol fluid are modeled using the Dissipative Particle Dynamics (DPD) method. The blood plasma is modeled as an incompressible fluid via the Immersed Boundary Method (IBM). This coupling is novel because it provides an accurate description of RBC dynamics while the extracellular flow patterns around the RBCs are also captured in detail. Our coupling methodology is validated with available experimental and computational data in the literature and shows excellent agreement. We explore the controlling regimes by varying the shape of the oscillatory flow waveform at the channel inlet. Our simulation results show that a host of RBC morphological dynamics emerges depending on the channel geometry, the incoming flow waveform, and the RBC initial location. Complex dynamics of RBC are induced by the flow waveform. Our results show that the RBC shape is strongly dependent on its initial location. Our results suggest that the controlling of oscillatory flows can be used to induce specific morphological shapes of RBCs and the surrounding fluid patterns in bio-engineering applications.


Introduction
Extensive research has been conducted in the last few decades on the morphological change of Red Blood Cells (RBCs) in fluid flows due to their importance in blood pathology (Kaul et al, 1983;Barabino et al, 2010;Secomb, 2017). It has been shown that the response of the RBC membrane to blood plasma dynamics can affect the overall patterns of microvascular blood flows (Tomaiulo et al, 2009;Guckenberger et al, 2018;Reichel et al, 2019). Despite a substantial body of literature, the dynamics of RBCs remain a significant challenge to be studied due to the complexity of various response modes, which result from the interaction of the suspended cellular membrane with shear flow (Vlahovska et al, 2011). Several factors can affect the dynamics of RBCs, such as stiffness of the membrane (Czaja et al, 2020), shear rate γ (Lanotte et al, 2016;Mauer et al, 2018), and viscosity contrast (λ) between the blood plasma and the cytosol (Mauer et al, 2018), among other factors. As a result, the RBC deformation process in shear flow is not well understood, especially under time-dependent shear rates (Schmidt et al, 2022;Krauss et al, 2022).
In free shear flows with a constant shear rate (γ 0 ), the shear strength is the controlling parameter of the RBC dynamics, in which the shape of RBCs becomes increasingly complex (more lobes) as the shear rate increases (Dupire et al, 2012). In particular, for the range of shear rate (γ 0 ) from 10 s −1 to 2,000 s −1 , the dynamics of RBCs can be classified into three main regions (Lanotte et al, 2016): (i) tumbling at weak shear rate (γ 0 < 10 s −1 ); (ii) circular/elliptical stomatocytes (10 s −1 <γ 0 < 400 s −1 ); and (iii) multilobes (400 s −1 <γ 0 < 2, 000 s −1 ). In the tumbling region, the RBC deformation is minimal and reversible, which allows the RBCs to maintain their biconcave discoid shape. As the shear rate increases to 400 s −1 , the percentage of discocytes decreases and stomatocytes start to dominate. The rolling and tumbling stomatocytes appear atγ 0 = 150 s −1 and 250 s −1 , respectively. This pattern persists up until γ 0 = 400 s −1 when the stomatocytes assume a shape with an elliptical rim. In the range 400 s −1 <γ 0 < 2, 000 s −1 , RBCs with large lobes on their surface, which are referred to as trilobes or hexalobes, emerge.
Studies of RBC dynamics in microchannels have shown that the RBC can transition from its biconcave discoid shape to different morphologies (Noguchi † and Gompper, 2005;Fedosov et al, 2014;Guckenberger et al, 2018) under specific combinations (state diagram) of viscosity contrast, shear rate (the capillary number -C a ) (Kaoui et al, 2011a) and channel confinement (χ) (Mauer et al, 2018;Fedosov et al, 2014;Reichel et al, 2019). The state diagram has revealed two main categories of RBC morphological shapes: (i) symmetrical; and (ii) asymmetrical types (Tomaiulo et al, 2009;Quint et al, 2017;Kihm et al, 2018). The simple type contains three modes (Coupier et al, 2012): (a) bullet; (b) croissant (in rectangular channels); and (c) parachute (in circular channels) shapes, while the complex type includes (Kaoui et al, 2009;Reichel et al, 2019): (a) slipper; (b) multilobes; (c) trilobes; and (d) hexalobes shapes. The shape transition in the simple type has been shown to reach a terminal shape (either bullet, croissant, or parachute). However, it is still not fully clear whether or not the complex shapes are stable or they are just transient states (Noguchi † and Gompper, 2005;Kaoui et al, 2011b). In the complex shape mode (Kaoui et al, 2009), the shape transition mostly depends on the flow lag, which is the difference between the translation velocity of the RBC and the velocity of blood plasma. In brief, it is unclear how the complex shapes emerge from the biconcave discoid shape.
Two shapes are the most frequently observed (Guckenberger et al, 2018): (i) the croissant shape (symmetrical); and (ii) the slipper shape (asymmetrical). In particular, the slipper shape is characterized by the tank-treading motion of the cell membrane, which is essentially a self-rotation of the membrane around its center of mass during the RBC propagation (McWhirter et al, 2009;Guckenberger et al, 2018). Experimental and computational studies have shown that these morphological shapes might result in distinct flow structures of blood plasma in the vicinity of the RBC (Guckenberger et al, 2018). For instance, there exists a closed vortex downstream of the RBC when the slipper shape emerges (Yaya et al, 2021). Such a vortex is absent during the croissant shape. To our knowledge, there has been no systematic effort in understanding the emergence of the extracellular flow patterns as the morphological shape of the RBC changes.
Recently, oscillatory flow (time-dependent shear rate -γ f ) has been shown to be a promising technique for cell separation. Because cell deformation is irreversible under time-dependent shear rates (Mutlu et al, 2018;Schmidt et al, 2022), oscillatory flows have been suggested to sort RBCs based on their size and deformability (Krauss et al, 2022). Oscillatory flows can reduce the required travel distance of cells because it induces the lateral migration of cells in a short axial distance. This feature simplifies the design of microfluidic channels and thus improves the cell separation process (Lafzi et al, 2020). However, it is unclear on the process of morphological transition as the RBC responds to the time-dependent shear rate (γ f ) during this lateral migration.
Therefore, it is necessary to investigate this process in detail.
In this work, we utilized our hybrid continuum-particle simulation methodology (Akerkouch and Le, 2021) to study the response of the RBC to time-dependent shear rates. Our paper is organized as follows. First, a brief description of the numerical methods for simulating the blood plasma and the RBC is presented. Second, the obtained RBC dynamics are validated with experimental data under (i) stretching force; (ii) constant shear rates (croissant and slipper shapes); (iii) oscillatory shear rates. Third, we perform a parametric study where the shear rate waveform, the peak flow rate, and the initial position of the RBC are varied to induce a host of RBC morphological changes. Finally, the relationship between the RBC's shape and the extracellular flow patterns is reported as a basis for cell manipulation in future applications.

The idealized shape of the RBC
The idealized shape of the RBC membrane is given by a set of points with coordinates (x, y, z) in 3D space with the analytical equation (Fedosov et al, 2010) as shown in Figure 1: the parameters are chosen in this work as D 0 = 7.82 µm (the equilibrium diameter), a 0 = 0.00518, a 1 = 2.0026, and a 2 = −4.491. Note that the idealized shape will be used as the initial shape of the RBC membrane for all simulation cases. The membrane mechanics that govern the cellular deformation under loadings will be described in the following sections.

RBC membrane model
As the idealized surface of the RBC membrane is known precisely according to Equation (1), a triangulation procedure is carried out to mimic the distribution of the spectrin links on the membrane as edges of each triangular element (links) (Fedosov et al, 2010). A network of non-linear springs is generated for each edge to model the dynamics of the spectrin links (Pivkin and Karniadakis, 2008;Fedosov et al, 2010;Akerkouch and Le, 2021). At each vertex i, the dynamics of the links are determined by the membrane force F membrane i , which is linked to Helmholtz's free energy V i at the same vertex i through the following relationship: with r i is the position vector of the vertex i.
The potential V ({r i }) incorporates the physical properties of the lipid bilayer: (a) in-plane stretching; (b) bending stiffness; (c) area and volume conservation; and (d)

Action potential models
The in-plane free energy term V in−plane includes the elastic energy stored in the membrane modeled using the nonlinear Wormlike Chain and Power (W LC −P OW ) spring model. Here, the W LC − P OW potential is computed for each link j formed by two vertices as, with N s as the total number of links forming the triangulated mesh.
The W LC attractive potentials U W LC (l j ) for individual link l j is expressed as: where the value x = lj lmax represents the spring deformation, in which l j , l max , p, k B , and T are the length of the spring j, the maximum allowable length for the links, the persistence length, Boltzmann's constant, and the temperature, respectively. The repulsive force, described by the energy potential U P OW (l j ), takes the form of a power function (P OW ). The separation distance l j is a determining factor in the calculation of U P OW , which is given by: where k p is the P OW force coefficient. In this work, the value m = 2 is used for the exponent (Fedosov et al, 2010).
The bending energy V bending accounting for the membrane resistance to bending is defined as, with k b , θ 0 and θ j are the bending rigidity, the spontaneous angle and the instantaneous angle between the normal vectors of two adjacent triangles sharing a common edge (link) j, respectively.
The area and volume conservation constraints account for the incompressibility of the lipid bilayer and the inner cytosol, respectively. They are defined as: with N t as the total number of triangles. k a , k d , and k v are the global area, local area and volume constraint coefficients, respectively. A k and A 0 are the instantaneous area of the k th triangle (element) and the initial value of the average area per element.
A tot 0 and V tot 0 are the RBC's equilibrium total area and volume, respectively. A and V are the instantaneous total area and total volume of the RBC. The detailed procedure to evaluate the values of A and V for individual elements was reported in our previous work (Akerkouch and Le, 2021).
Equation (2) is used to calculate the precise nodal forces for each potential energy V in Equations (4) -(8) (Akerkouch and Le, 2021;Fedosov et al, 2010). The internal force F membrane i contribution from i th vertex can be computed by summing all the nodal forces as:

Cellular membrane/cytoskeleton interaction
To account for the interaction between the cytoskeleton and the lipid bilayer, the bilayer-cytoskeletal interactions force F E is incorporated into the total RBC membrane forces (Peng et al, 2013a). The membrane thickness of a healthy RBC is typically 10 nm, composed of an outer layer (lipid bilayer) and the inner layer (skeleton) (Li et al, 2018). In this study, F E is applied when the distance between two elements with opposite normal vectors is at least 20 times the thickness of the membrane. The force F E is applied equally to all the vertices (i = 1, 2 and 3) for each of the two elements.
The bilayer-cytoskeletal interactions force is given by: with the stiffness of the bilayer-cytoskeletal (k bs = 4.12 pN/µm) is assumed to be in the same order as one of the membrane spectrin networks (Peng et al, 2013a). n is the normal vector of the triangle.

Modeling membrane-cytosol interactions
The interaction between the membrane and the cytosol 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 (Fedosov et al, 2010). 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 (Pivkin and Karniadakis, 2008;Fedosov et al, 2010) as shown in Figure 1.
Due to the different nature of the interactions, the component of the total force of each particle F i depends on the nature of the particle i (either membrane or cytosol particle). In general, each DPD particle i interact with surrounding particles j within a cut-off radius r c through three pairwise additive forces: (a) the conservative force F C ij ; (b) the dissipative force F D ij ; and (c) the random force F R ij . 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 vectorr ij = r ij rij . Also, v i,j = v i − v j is the relative velocity between the particles i and j with velocities v i and v j .
For a DPD particle i of the cytosol fluid, the total force F i is: For membrane particles, the total force F i acting on each membrane particle is given by the sum of the membrane force F membrane i , the bilayer-cytoskeletal interactions force F E and the contributing forces from the surrounding DPD fluid particles (the cytosol): The mathematical formulation of the conservative force F C ij , the dissipative force F D ij , and the random force F R ij for the membrane and the cytosol fluid particles are explained below.

The conservative force
The conservative force F C ij is given by : where a ij = 20 is the conservative force coefficient between particles i and j. Note that the particles i and j can be either membrane or cytosol fluid particles. Thus, there are two types of interactions: i) cytosol fluid/fluid; and ii) membrane/fluid-particle interactions (Fedosov et al, 2010).

The dissipative force
The dissipative force F D ij for the membrane particles is computed as: The membrane viscosity is a function of both dissipative parameters, Γ T and Γ C .
The superscripts T and C denote the translational and central components. Here, Γ T is responsible for a large portion of the membrane viscosity in comparison to Γ C . In addition, Γ C is assumed to be equal to one-third of Γ T in Equations (14) (Fedosov et al, 2010). Consequently, these parameters relate to the physical viscosity of the membrane η m as: Hence, the dissipative force F D ij among the membrane particles can be expressed as: The dissipative force F D ij for the cytosol fluid particles is defined as: the quantity γ is a constant coefficient defining the strength of the dissipative force.
The weight functions, ω D (r ij ) and ω R (r ij ) are given by: with s = 1 following the original DPD method (Fedosov et al, 2010). The particle i represents the fluid particle, while the particle j can be either a fluid or membrane particle within the cut-off radius r c .

The random force
Using the assumptions in Equation (15), the random force for membrane particles can be described as: where tr(dW ij ) is the trace of the random matrix of independent Wiener increments dW ij , and dW S ij = dW S ij − tr(dW S ij ) 3 is the traceless symmetric part.
The random force F R ij for the cytosol fluid are defined as: where σ is a constant coefficient defining the strength of the random force, dt is the physical time step, ϑ is a normally distributed random variable with zero mean and unit variance and ϑ ij = ϑ ji . Note that the particles i and j must be both cytosol fluid particles.

Plasma and cytosol viscosity contrast
At the physiological condition, the viscosity ratio between the blood plasma and the RBC cytosol is equal to 5.0 (λ = µ cytosol µ plasma = 5.0) (Wells and Schmid-Schönbein, 1969). To ensure that this condition is met, the dynamic viscosity of the plasma is set to be µ plasma = 1.2 mP a s. The viscosity condition is enforced on the cytosol fluid by calibrating the parameters of the dissipative and the random forces (Mai-Duy et al, 2020) (γ and σ). Specifically, the dynamic properties of the DPD particles of the cytosol fluid are given in the dimensionless DPD unit as (Fan et al, 2006): with ρ as the density. The DPD dimensionless parameters and physical units are linked (Ghoufi et al, 2013) in order to compute the coefficients of the dissipative and randoms forces for the cytosol dynamic viscosity of µ cytosol = 6 mP a s. The details of the conversion procedure are summarized in Table 1.

Scaling of model and physical units
One challenge in DPD modeling is establishing a relationship between the modeled quantities and the physical values (Ye et al, 2019). Since this relationship is not explicit, it is necessary to use a scaling argument to recover this relationship (Fedosov et al, 2010;Peng et al, 2013b). For each parameter, the superscript M and P correspond to the model and physical units, respectively. The details of the physical parameters for the RBC are shown in the Tables 2.
The length scale r M is defined as: The energy per unit mass k B T and the force N scaling values are given by: with Y is the membrane Young's modulus.
The timescale τ is defined as follows: with α = 1 is the scaling exponent.

Coarse-graining procedure
A full-scale model of the RBC typically consists of millions of particles, which are required to accurately simulate protein dynamics (Tang et al, 2017). However, it is not feasible to use such a full-scale model in a Fluid-Structure Interaction (FSI) simulation due to the high computational cost. We followed the coarse-graining procedure of Pivkin and Karniadakis (2008) to represent the RBC membrane by a smaller number of particles (coarse-grained model). This procedure does not allow a detailed simulation of separate proteins, but it is versatile enough to capture the overall dynamics of the RBC membrane. The parameters of the coarse-grained model (c) are computed from the ones of the fine-scaled model (f ) by a scaling procedure. Examples of such parameters are explained below.
Based on the equilibrium condition, Pivkin and Karniadakis (2008) proposed a coarse-graining procedure based on the area/volume constraint for the spring equilibrium l 0 and maximum l max lengths as follows: the role of l 0 and l max is critical in determining the response from the WLC model as seen in Equation (5), with l max = 2.2 l 0 in the fine-scaled model. Due to the scaling in Equation (26), the value of x 0 = l0 lmax = 1 2.2 does not change as the model is coarsegrained from the number of vertices N f v to N c v . In this work, a range of N c v values are considered as shown in Table 3 with N f v set to be 27, 472.
Furthermore, as the number of vertices reduces, the average angle between the pairs of adjacent triangles increases. Therefore, the spontaneous angle θ is adjusted accordingly in the coarse-grained model as: To maintain the shear and area-compression moduli, the parameters p and k p are adjusted as: The details of the parameters for the coarse-graining process are shown in Table 3.

Time integration
We implemented the modified Velocity-Verlet algorithm (Groot and Warren, 1997), which consists of two primary steps. The first step involves determining the new position of the particle i (r i ) while predicting the velocity (ṽ i ), and the second step involves correcting the velocity by utilizing the computed force (F i ) based on the predicted velocity and the new position as follows.
whereṽ i (t + dt) is the predictive velocity at time t + dt and Λ is the variable which accounts for the effects of the stochastic processes. The value of Λ is chosen to be the optimal value (Groot and Warren, 1997) Λ = 0.65.

Fluid-Structure Interaction simulation of RBC in flows
The blood plasma is considered as an incompressible Newtonian fluid modeled using the incompressible three-dimensional unsteady Navier-Stockes equations, with density ρ and kinematic viscosity ν = µ plasma ρ . The governing equations (continuity and momentum) are read in Cartesian tensor notation as follows (i = 1, 2, 3 and repeated indices imply summation): In the above equations, u i is the i th component of the velocity vector u; t is time; x i is the i th spatial coordinate; p is the pressure divided by ρ. The characteristic velocity scale is chosen as U 0 . The length scale L s is set to equal 8 µm for all cases. Note that this length scale is chosen to reflect the diameter of the RBC at the equilibrium The fluid solver is based on the sharp-interface curvilinear-immersed boundary (CURVIB) method in a background curvilinear domain that contains the RBC model (Ge and Sotiropoulos, 2007). The CURVIB method used here has been applied and validated in various FSI problems across different biological engineering areas (Le et al, 2010;Le and Sotiropoulos, 2012;Le et al, 2019). In our previous work (Akerkouch and Le, 2021), we utilize the capabilities of the CURVIB method to capture accurately the complex cellular dynamics of the RBC in fluid flows.

Computational setups
Fluid-Structure Interaction simulations are performed to determine the dynamics of RBC in a confined micro-channel (Guckenberger et al, 2018). The computation domain is defined as a rectangular channel containing a single RBC as illustrated in Figure 2a.
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 details of the channels used in the simulations are listed in Table 4.
The RBC locates initially at t = 0 in an axial distance of x 0 from the inlet. The transverse location of the RBC is placed along the bisector of the first quadrant with a radial shift (r). Thus, the transverse coordinates of the RBC are y 0 = r and z 0 = r, respectively as shown in Figure 2b (r is the radial shift). With this configuration, the RBC confinement is defined as the ratio between the effective RBC diameter π and the domain height L z : The initial shape of the RBC is first set to be the idealized shape (Equation 1) for all simulation cases at the initial time t = 0. A short period of relaxation t relax is allowed for the RBC under no external load (no flows) so that the internal forces of the RBC membrane balance. A uniform flow velocity U (t) is then applied at the channel inlet at t > t relax to induce the RBC's deformation depending on the controlling strategy.
The average shear rate across the channel height is defined as the ratio between the bulk velocity U (t) and the domain's height:

Constant shear rate condition (I 0 )
Following the experimental study of Guckenberger et al (2018) (Channel-1, Table 4), FSI simulations of the RBC in channel flow with a constant flow rate are carried out with x 0 = 22.5 µm. To highlight the constant flow rate, the notation I 0 is introduced to emphasize this condition. As shown in Table 5, a constant inflow velocity U (t) = ψ 0 is required at the inlet of the computational domain. Two values of ψ 0 are considered: (i) ψ 0 = U 3 = 2 mm/s; and (ii) ψ 0 = U 4 = 6 mm/s. In these cases, two values of the radial shift are also investigated: r 1 = 0 and r 3 = 0.7 µm. To simplify the discussions, the numerical values for the bulk velocity ψ 0 will not be explicitly referred to. Instead, only the acronyms (U 3 and U 4 ) will be used for reasons that will be evident in the subsequent texts.
Using these notations, the FSI simulation cases are named using the convention for each type of inflow waveform (I), the bulk velocity (U ), the radial shift (r), and the channel type, respectively. The first case (I 0 U 3 r 1 χ 1 ) is configured with (ψ 0 = U 3 = 2 mm/s) and r = r 1 = 0 µm. The second case (I 0 U 4 r 3 χ 1 ) is carried out with ψ 0 = U 4 = 6 mm/s and r = r 3 = 0.7 µm. Here, the Reynolds number is defined as Re = U Ls ν . The kinematic fluid viscosity of blood plasma is chosen as ν = µ plasma ρ = 1.2 × 10 −6 m 2 /s. The summary of the parameters for each simulation case is shown in Table 5.
First, t relax = 10 ms and 7.0 ms are set for I 0 U 3 r 1 χ 1 (croissant) and I 0 U 4 r 3 χ 1 (slipper) simulations, respectively. After the relaxation period, a linear ramping period is set for each simulation case t ramp = 30 ms and 20 ms are set for I 0 U 3 r 1 χ 1 and I 0 U 4 r 3 χ 1 , respectively. During this ramping period, the bulk velocity U (t) is linearly increased. The value of U (t) reaches ψ 0 at the end of the ramping period.

Stepwise oscillatory flows (I s )
To further validate our FSI model in oscillatory flows, the propulsion of the RBC in square channels is investigated (Schmidt et al, 2022). Two square channels (Channel-2 and Channel-3 in Table 4) with side lengths L z = 16 µm and 21 µm are used for the simulations, resulting in confinements χ 2 = 0.4 and χ 3 = 0.3, respectively. The initial location of the RBC is on the channel axis (x 0 = 16 µm, r = r 1 = 0). The computational configuration including the grid spacing, RBC surface meshes, and boundary conditions are shown in Figure 2a and Table 4. A stepwise asymmetric oscillatory waveform I s is used with two phases: (i) forward T f ; and (ii) backward T b periods T b T f = 4 as shown in Figure 3a. The velocities during the forward and backward phases are ψ f and ψ b − ψ f ψ b = 4 , respectively. The formula for the waveform is defined as: Following this formula, the flow has a forward phase (ψ f > 0) and a backward flow Lz . The Capillary number Ca f in the forward flow phase is given as: and µ 0 is the shear elastic modulus described in Table 2.
There are 6 values of ψ f are examined as ψ 1 = 1.05, ψ 2 = 1.58, ψ 3 = 2.1, ψ 4 = 1.9, ψ 5 = 2.8, and ψ 6 = 3.7 mm/s, respectively. Following the naming convention of the simulations, six cases are formed with the respective parameters: I s ψ 1 r 1 χ 2 ; I s ψ 2 r 1 χ 2 ; I s ψ 3 r 1 χ 2 ; I s ψ 4 r 1 χ 3 , I s ψ 5 r 1 χ 3 , I s ψ 6 r 1 χ 3 as shown in Table 5. As the waveform applied is of a stepwise nature, there is no relaxation time taken into consideration for these 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 = T 5 ), as a function of the forward (peak) capillary number Ca f for the chosen shear rates (Schmidt et al, 2022). Thus ∆x c is defined as the displacement of the RBC's centroid (C) at the end of the forward phase (t = T 5 ):

Sinusoidal flow simulations
To study the effect of the pulsatile flow on the propulsion and the behaviour of the cellular response (morphology changes) of the RBC, we considered time-periodic flow The symmetric waveform (I 1 ) is created with T f = T b (completely symmetric). The asymmetric waveforms (I 2 , I 3 , and I 4 ) are formed by progressively reducing the period of T b . Four distinct inflow types were generated with symmetry and asymmetric waveforms (I 1 , I 2 , I 3 , and I 4 ) as seen in Figure 4 and Tables 6 and 7. For each of these waveforms, three different velocity magnitudes (A = U 1 , U 2 and U 3 ) and three different radial shifts (r 1 , r 2 and r 3 ) are considered as shown in Tables 8. In total, the combinatoric arrangements lead to a total of 36 distinct simulation cases with the notation I m U n r p χ 1 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 8. In addition, the RBC shapes are recorded over a time period of two cycles 2T as exemplified in Figure   4b, in which the initial location of the RBC is set at x 0 = 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: 3 Results

Coarse-graining validation
To first validate the coarse-graining procedure employed in our study, a stretching test is carried out and aimed to replicate the experimental test of Mills et al (2004). In our simulations, the parameters to describe the physical characteristics of the RBC are listed in Table 2. 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 3. (ii) the slipper shape (I 0 U 4 r 3 χ 1 -γ 0 = 600 s −1 ) as shown in Figure 5 and Figure 6.
Under low shear rate (I 0 U 3 r 1 χ 1 ), the RBC is initially placed along the centerline of the microchannel (discocyte shape). As the RBC interacts with the incoming flow, deforms, and eventually transitions to a croissant shape as shown in Figure 6a. 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.
Under high shear rate (I 0 U 4 r 3 χ 1 ), the RBC transitions to the slipper shape as shown in Figure 6b, which exhibits a bistability mode with tank-treading behaviour.
Note that the RBC is placed at a radial shift r 3 = 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 Figure 6b. A counter-clockwise rotation is observed as indicated by the locations of two membrane particles (Lagrangian points -V 1 and V 2 ) at different time instances (t 1 = 22 ms and t 2 = 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 Figure 6. These transitions agree well with the corresponding experimental data of Guckenberger et al (2018) as well as described in recent experiments on RBC transient dynamics (Recktenwald et al, 2022;Prado et al, 2015). Furthermore, these transitions are in good agreement with the shape diagram produced by Agarwal and Biros (2022) for different Capillary numbers and confinements as seen in Figure 5.
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 of Guckenberger et al (2018). The extracellular flow pattern can be visualized by reconstructing the relative flow velocity field (Yaya et al, 2021) (co-moving frame). The relative velocity is defined as the difference between the flow velocity and the RBC's centroid velocity as shown in Figure 7. In the croissant shape (I 0 U 3 r 1 χ 1 ), the velocity streamlines closely resemble an axi-symmetrical flow pattern ( Figure 7a). The downstream side of the RBC membrane deforms significantly whereas the upstream side barely changes as depicted in Figure 7b. In the slipper shape (I 0 U 4 r 3 χ 1 ), 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. In all simulation cases (I s ψ 1 r 1 χ 3 , I s ψ 2 r 1 χ 3 , I s ψ 3 r 1 χ 3 , I s ψ 4 r 1 χ 3 , I s ψ 5 r 1 χ 3 , and I s ψ 6 r 1 χ 3 ), the RBC transitions from the idealized discocyte to the biconcave shape during the forward phase (0 < t < T 5 ) with all values of the peak forward flow (ψ f = 1.05 mm/s to ψ f = 4.34 mm/s) as shown in Figure 6c. 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 ψ 6 r 1 χ 3 in Figure 6c. 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 placed initially at the channel center (r = r 1 = 0). Therefore, it is possible to induce complex dynamics of RBCs in a confined channel by changing the inflow waveform only.

The emergence of RBC shapes
The sinusoidal flow waveform (I 1 , I 2 , I 3 , and I 4 ) 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 Figure 8 and Tables 9-12. The emergence of each type will be discussed as follows.
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 ( t T ≈ 0.25) as shown in Tables 9-12. This process is demonstrated for the cases I 1 U 1 r 1 χ 1 , I 3 U 2 r 1 χ 1 , I 4 U 2 r 1 χ 1 , and I 4 U 3 r 3 χ 1 as seen in Figure 8 (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 9-12, the S mode appears only when the RBC is initially placed not exactly at the cross-sectional center (r > 0) (for example, I 4 U 3 r 3 χ 1 ). 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 either I 1 , I 2 , I 3 , or I 4 as seen in cases I 1 U 1 r 1 χ 1 , I 3 U 2 r 1 χ 1 , and I 4 U 2 r 1 χ 1 . 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 ( t T > 0.5) as shown in Tables 9-12 and Figure 8 (middle column). The shape transformation is initiated by the buckling of the RBC membrane, which takes place in the resting interval (T r ) (see Figure 4) 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 Figure 8. 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 9-12. Thus the shape transition process is irreversible even if the inflow waveform is symmetric (I 1 ) as shown in Table 9.

The impacts of the initial position (the radial shiftr)
Our results in Figure 8 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 1 = 0 under both the symmetric (I 1 U 1 r 1 χ 1 ) and asymmetric (I 3 U 2 r 1 χ 1 and I 4 U 2 r 1 χ 1 ) waveforms in Figure 8a, Figure 8b, and Figure 8c, respectively. In these cases (I 1 U 1 r 1 χ 1 , I 3 U 2 r 1 χ 1 , and I 4 U 2 r 1 χ 1 ), 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 3 > 0 in the case I 4 U 3 r 3 χ 1 , the RBC remains mostly the slipper shape during the forward phase (t < T f ) as seen in Figure 8d. 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 3 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 below.
When the RBC is positioned at the center line (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 1 ) as shown in Figure 9a, 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 2 , I 3 , and I 4 ) 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 Figure 9a (left). The value of the axial displacement (∆x c ) (t) increases sequentially from I 1 , I 2 , I 3 to I 4 . ∆x c reaches a maximum value of approximately 4 × L s at the end of the second cycle for the case I 4 .
In the lateral direction, ∆y c reaches a value of approximately 0.16 × L s at the end of the first cycle for the symmetric waveform (I 1 ) as shown in Figure 9b (left). For the cases I 2 , I 3 , and I 4 the values of ∆y c remains comparably low (≈ 0.08 × L s ) during the cycles. Furthermore, the vertical displacement ∆z c as seen in Figure 9c follows a monotonically upward trend throughout the second cycle. This results in a maximum vertical displacement ∆z c of approximately 0.25 × L s for I 1 . For other waveforms (I 2 , I 3 , and I 4 ), 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 1 ) leads to a maximum displacement in the transverse directions while the asymmetric waveforms (I 2 , I 3 , and I 4 ) result in the maximum axial displacement for r = 0.
As the RBC is placed at the off-centered location (r = r 2 = 0.4 µm) as seen in The values of ∆y c and ∆z c reach 0.14 × L s for I 1 and I 2 at the end of the second cycle. Surprisingly, the impact of the applied waveform seems to diminish as the value of T f increases. Comparing the case I 3 U 2 r 2 χ 1 and I 4 U 2 r 2 χ 1 (I 3 and I 4 ), the values of ∆y c and ∆z c in Figure 9e and Figure 9f are nearly identical, resulting in a vertical displacement of approximately 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.

The impacts of the waveform (I)
It is striking to observe the irreversible dynamics of RBC in Figure 10. When the RBC is subjected to the symmetric waveform (I 1 ) with different radial shifts r = r 1 , r 2 and r 3 , the RBC oscillates around its initial position with a minimal displacement along the axial direction as shown in the cases I 1 U 1 r 1 χ 1 , I 1 U 1 r 2 χ 1 , and I 1 U 1 r 3 χ 1 . Despite the inflow waveform being completely symmetrical (a sine function -I 1 ), the axial position of the RBC in Figure 10a Figure 10c, 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 1 ). In conclusion, a symmetrical flow waveform (I 1 ) results in minimal propulsion along the axial direction but a significant lateral migration.
Under asymmetric waveform I 4 , the RBC propels along the channel direction with an axial displacement of approximately 2 × L s in each cycle as shown in Figure 10d.
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 to ones in the symmetric case (I 1 ). The values of (∆y c , ∆z c ) are within 0.15 × L s for all cases I 4 U 1 r 1 χ 1 , I 4 U 1 r 2 χ 1 , and I 4 U 1 r 3 χ 1 as shown in Figure 10e  The case I 1 U 1 r 1 χ 1 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 Figure 11a and Figure 11b. This case is chosen because the temporal variation of the waveform is completely symmetrical (I 1 ). Moreover, the RBC is placed initially at the channel axis (r = r 1 = 0) with the lowest forward velocity ψ f = U 1 = 1 mm/s. In the case I 1 U 1 r 1 χ 1 , Figure 11a 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 (Figure 11b). 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 1 U 3 r 3 χ 1 was selected to visualize the flow patterns. As shown in Figure 11c, 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 0 U 4 r 3 χ 1 with U 4 = 6 mm/s) in Figure   7b. This is remarkable since the peak flow ψ f is rather three times lower in this case ψ f = U 3 = 2 mm/s. The case I 1 U 3 r 1 χ 1 (Figures 11d, 11e, 11f, and 11g) is selected to illustrate further the impact of the peak forward flow ψ f . In this case, the peak velocity is ψ f = U 3 = 2 mm/s. In comparison to the case I 1 U 1 r 1 χ 1 (Figure 11a and Figure 11b), 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 1 U 1 r 1 χ 1 . During the forward phase, the RBC transitions quickly to the croissant shape in Figure 11d at (t = 0.28T ), just after the peak forward flow. The flow patterns are similar to those observed under case I 0 U 3 r 1 χ 1 in Figure 7a. 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 Figure 11e. The flow pattern is perturbed minimally surrounding the RBC as its shape turns to the trilobes shape as depicted in Figure 11f. 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 Table 9 to Table 12, 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 4 U 2 r 1 χ 1 , 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 12. Its corresponding flow patterns are shown in Figure   11h, 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 is the most commonly found as shown in Table 9 to Table 12. The flow pattern around a discocyte is exemplified in Figure 11i for the case I 1 U 1 r 2 χ 1 at t = 1.25T . Here, the flow streamlines show a distinct separation of upstream and downstream regions.

Discussion
Due to the membrane flexibility, RBC deforms swiftly under the shear flow (Lanotte et al, 2016). This morphological feature can be exploited to understand the mechanical properties of the RBC membrane (Prado et al, 2015) and thus it has the potential to identify the pathological changes (Recktenwald et al, 2022) of RBC's membrane.
However, the exact mechanism of this response is not yet fully understood. In this work, we explore the impacts of the unsteady shear rate to control cell deformation and migration in microchannels.
Our numerical method is based on the continuum-particle coupling (Akerkouch and Le, 2021), which allows the simulations of RBC dynamics under physiological conditions. Our numerical results show excellent agreements with available in vitro and computational data both in cellular mechanics and extracellular flow pattern of the blood plasma (Mills et al, 2004;Guckenberger et al, 2018;Yaya et al, 2021). While most previous studies (Fedosov et al, 2010;Pivkin and Karniadakis, 2008) have only focused on the impact of constant shear rate on the dynamics of the RBCs, our results show that the unsteady shear rate can induce complex RBC's morphology in confined channels as discussed below.

The emergence of the croissant shape and the slipper shape under a constant shear rateγ 0 )
In micro-channel flows with constant shear rate (γ 0 ), three common dynamics of RBCs are frequently observed: (i) tumbling; (ii) croissant/parachute; and (iii) slipper shapes as shown in Figure 5. In unconfined flows (Dupire et al, 2012), the RBC dynamics depends on only the shear rate (γ or Ca) and viscosity contrast (λ). However, the confinement of micro-channel flows imposes an additional condition for shape transition via the confinement ratio χ. Our results in Figure 5 further confirm that the combination of C a and χ dictates the emergence of RBC shape in either the croissant or slipper shapes.
Recent works (Guckenberger et al, 2018;Yaya et al, 2021) in rectangular microchannels, which are identical to our channels as shown in Figure 2 and Table 4, also suggest that the emergence of RBC shape is dependent on the radial shift (r -see Figure 2 for its definition). In previous works (Guckenberger et al, 2018;Yaya et al, 2021), the croissant shape emerged atγ 0 < 300 s −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 (γ 0 ≈ 500 s −1 ) and off-centered placement (r > 0) (Guckenberger et al, 2018;Yaya et al, 2021). In cylindrical micro-channels (Fedosov et al, 2014), similar observations were confirmed albeit at lower shear rates (0 <γ 0 < 80 s −1 ). These works point to the importance of the radial shift in regulating RBC dynamics. Our results in Figure 5 agree with these observations. Our results show the appearance of croissant-to-slipper transition as the Capillary number (and thusγ 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 −1 as illustrated in Figure 6. Therefore, our results further confirm the importance of the radial shift in the croissant-to-slipper transition.

The impact of time-varying shear rateγ(t) on RBC shape
When the inflow varies in a stepwise manner as seen in Figure 3, the shear rate (γ) changes as a function of timeγ(t) with distinct forward (T f ) and backward (T b ) time phases. In all cases (I s ψ 1 r 1 χ 2 , I s ψ 2 r 1 χ 2 , I s ψ 3 r 1 χ 2 , I s ψ 4 r 1 χ 3 , I s ψ 5 r 1 χ 3 , and I s ψ 6 r 1 χ 3 ), the RBC is placed exactly at the channel axis (r = r 1 = 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 Figure 6. 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 Figure   6 at the end of T b . This is remarkable given that the maximum shear rate during the backward phase can be sufficiently large (γ f = 200 s −1 ) but this symmetry is still maintained. Comparing the case I 0 U 3 r 1 χ 1 and I 0 U 4 r 3 χ 1 in Table 5, our results suggest that the complete break of symmetry (Recktenwald et al, 2022) (slipper shape) is observed only when the radial shift exists (r > 0).
When applying different sinusoidal waveforms (I 1 , I 2 , I 3 and I 4 ) shown in Figure   4, our results show the ubiquitous presence of croissant and slipper shapes across all shear rates (γ f = 100, 150, and 200 s −1 ). The slipper shape appears at (t ≈ 0.3T ) whenever the RBC is placed off the channel's axis (r > 0) as shown in Tables 9-12.
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 4 . As shown in Table 9, the slipper shape is observed even when the waveform is completely symmetric (I 1 ) given that r > 0 as in: (a) I 1 U 1 r 2 χ 1 ; (b) I 1 U 1 r 3 χ 1 ; (c) I 1 U 2 r 2 χ 1 ; (d) I 1 U 2 r 3 χ 1 ; (e) I 1 U 3 r 2 χ 1 ; and (f) I 1 U 3 r 3 χ 1 . 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 (Lanotte et al, 2016) that the presence of the discocyte shape is correlated with weak shear rates γ 0 . Underγ 0 < 10 s −1 , the RBC typically maintains its discocyte shape with an 80% probability (Lanotte et al, 2016). However, as the shear rate gradually rises (10 s −1 ≤γ 0 ≤ 400 s −1 ) the likelihood of a discocyte shape decreases to 30%. This observation has been validated even when different viscosity ratios (λ) (Mauer et al, 2018) are considered. In our work, the discocyte shape is ubiquitously observed during the second flow cycle across all applied waveforms (I 1 ,I 2 , I 3 , and I 4 ) and shear rates (γ f = 100 s −1 , 150 s −1 , and 200 s −1 ) with (r > 0) as shown in Tables 9-12. Therefore, our results further confirm that discocyte is the most common shape with the range shear rate less than 200 s −1 .
In our study, complex shapes evolve from simpler shapes under the influence of time-dependent waveforms. In the experimental work of Lanotte et al (2016) under constant shear rates, stomatocytes shape was observed to dominate the RBC population (65%) when the shear rate was between 10 s −1 <γ 0 < 400 s −1 . In our study, the elliptical-rim-shaped stomatocytes only emerge under symmetric waveform (I 1 ) in the case I 1 U 1 r 1 χ 1 under the peak shear rate ofγ f = 100 s −1 . At high constant shear rates 400 s −1 <γ 0 < 2, 000 s −1 , experimental data (Lanotte et al, 2016) 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 1 = 0) even at weak shear rates (γ f ≤ 200 s −1 ) as in Figure 8 and Tables 9-12. For example, the multilobes are observed with all waveforms. The trilobes are both observed in the symmetric waveform (I 1 U 2 r 1 χ 1 and I 1 U 3 r 1 χ 1 - Table 9) or the asymmetric waveform (I 3 U 2 r 1 χ 1 - Table 11). The hexalobes shape only appears under the most asymmetric waveform with (I 4 U 2 r 1 χ 1 ) −γ f = 150 s −1 as shown in Table 12. In the previous work of Li et al (2014), the RBC shape has been reported to deform further into elongated shapes as the shear rate increases. Our results support this observation as shown in Figure   11g as the elongated multilobes appear during the backward flow phase. Examining Table 9 to Table 12, 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γ f = 150 s −1 . Specifically, the elongated multilobes are observed under I 1 U 2 r 1 χ 1 , I 1 U 3 r 1 χ 1 , I 2 U 3 r 2 χ 1 , I 2 U 3 r 3 χ 1 , I 3 U 3 r 2 χ 1 , and I 3 U 3 r 3 χ 1 . 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.

Controlling lateral migration of cells with oscillatory flows
Microfluidic devices are typically used to isolate and separate cells (Gossett et al, 2010).
While these devices are promising for many cell-sorting applications (Suresh, 2007;Brandao et al, 2003), 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 (Amirouche et al, 2020) has been proposed as one efficient way.
Recent works have shown that varying the shear rates in time (Schmidt et al, 2022;Krauss et al, 2022) can reduce the required length based on the concept of velocity lift (Qi and Shaqfeh, 2017), which is the factor that drives the RBC's migration towards 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 by Krauss et al (2022) 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 Figures 10 and 9 for the symmetric waveform (I 1 ) 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 Figure 9. Our results in Figures 10b and 10e indicate that the RBC undergoes lateral migration. This migration is particularly evident in the cases I 1 U 1 r 3 χ 1 and I 4 U 1 r 3 χ 1 , where the cell migrates in the lateral direction for a distance of approximately 0.25L s and 0.13L s , respectively. Comparing the shapes of the waveform I 1 , I 2 , I 3 , and I 4 in Figure 4, our results suggest that significant lateral migration can be induced by adjusting the forward time interval (T f ). With our sinusoidal waveforms in Figure 4, T f should be less than three times the backward flow phase T b ( 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 Figure 7 and Figure 11 show that the extracellular flow patterns are directly influenced by the dynamics of the RBC. Under a stationary condition in Figure 7, the extracellular flow dynamics in the croissant shape is distinctively different from the ones of the slipper shape. In particular, the flow around the steady croissant shape was found similar to that of a rigid sphere (Lee et al, 2010), in which the flow streamlines move nearly symmetrically inwards 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 (Guckenberger et al, 2018). Our results in Figure 11 suggest that the extracellular flow patterns are far more complicated and highly dependent on the type of inflow waveform in timedependent shear rates. It is important to note that the extracellular flow has been found to play an important role in drug delivery strategies (Yaya et al, 2021) 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.x therapeutic nano-particles) to the RBCs.

Conclusion
Transient dynamics of Red Blood Cells (RBC) in confined channels under oscillatory flows are investigated using our continuum-particle approach (Akerkouch and Le, 2021). Our results reveal that the dynamics of RBCs are complex with different shape modes that are beyond the usually observed croissant and slipper modes. Our results indicate that the extracellular flow pattern around the RBC is dependent on the RBC shape. Our results suggest that the oscillatory flow can be used to control and manipulate the dynamics of RBC by adapting the appropriate flow waveform. Our specific conclusions are: • The RBC can transform into a variety of shapes such as multilobes, trilobes and hexalobes by varying the sinusoidal waveform even when it is subjected to a relatively weak shear rate (γ f ≤ 200 s −1 ).
• Simple shapes such as croissant, slipper, and rolling discocyte appear when the RBC is subjected to all waveforms. However, complex shapes such as rolling stomatocytes, trilobes, and hexalobes appeared only under specific conditions. In our study, the RBC transitions into 8 shapes under the symmetric waveform (I 1 ), and into 5 shapes under the asymmetric waveform (I 2 ). Therefore, it is possible to attain a certain shape using an appropriate waveform.
• Under the symmetric flow waveform, the axial displacement of the RBC is rather minimal. However, the lateral displacements are significantly large. Under the asymmetric flow waveform, the RBC experiences a large axial displacement but small lateral displacements.
• The maximum lateral displacement of the RBC during its propagation depends on the initial radial shift (r). This maximum value is also dependent on the asymmetry of the flow waveform (I).
• The extracellular flow surrounding the RBC depends on its morphological shape.
The flow pattern is thus distinct and unique for each shape.  Table 1 The relationship between DPD parameters and the physical units for viscosity ratio λ = 5. Nm, m, δt and µ cytosol correspond to the number of molecules in one bead, mass, time step and dynamic viscosity of the cytosol, respectively. V is the volume of the water molecule (30Å), M is the molar weight of water (18 g mol −1 ) and N A = 6.0221415×10 23 is the Avogardo's constant. The definitions of the parameters k B T , γ, ρ and rc are explained in Table 2 Table 4 Different channel geometries and their associated computational grids to simulate the dynamics of RBC in fluid flows. The channels have rectangular cross-sections of size Lx, Ly, and Lz along the axial, spanwise, and vertical directions, respectively. N i , N j and N k are respectively the number of grid points in x, y, z directions. χ is the channel confinement, which is defined in section 2.8.
Case  Table 5 Summary of the validation cases under constant shear rates (I 0 U 3 r 1 χ 1 and I 0 U 4 r 3 χ 1 ), and stepwise oscillatory flows (Isψ 1 r 1 χ 2 , Isψ 2 r 1 χ 2 , Isψ 3 r 1 χ 2 and Isψ 4 r 1 χ 3 , Isψ 5 r 1 χ 3 , Isψ 6 r 1 χ 3 ). The stepwise oscillatory flows with the forward (ψ f ), backward (ψ b ) velocities and the forward Capillary number are defined in section 2.8.2. The maximum shear rate (γ f ) and the maximum Reynolds number (Re f ) are defined in section 2.8.1. The definition of the RBC's radial shift r is shown in Figure 2b. Table 6 The controlling parameters of the pulsatile waveforms (I 1 , I 2 , I 3 and I 4 ). The waveforms are characterized by the intervals of the forward (T f ), rest (Tr) and backward (T b ) periods. The shapes of the waveforms are shown in Figure 4.
Tr ( Inflow waveform ψ f (mm/s) Radial shift r(µm) 1 I 1 U 1 = 1 r 1 = 0 2 I 2 U 2 = 1.5 r 2 = 0.4 3 I 3 U 3 = 2 r 3 = 0.7 4 I 4 U 4 = 6 - Table 7 The summary of the combinatoric configurations for steady and pulsatile flow simulations. The combination of the waveform type, the forward flow velocity, the radial shift, and the channel confinement results in the simulation configurations of ImUnrpχs. Here m = 1, 2, 3, and 4, n = 1, 2, 3 and 4, p = 1, 2 and 3 and s = 1, 2 and 3. The profile of the inflow waveforms (I 1 , I 2 , I 3 and I 4 ) are shown in Figure 4. The peak forward flow velocity U ( Figure 4a) varies from 1 to 6 mm/s. The radial shift of the RBC centroid along the bisector of the y − z plane at the initial time is defined in Figure 2b. 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 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  ), (b), (c), or (d) each consist of 9 possible combinations between the peak forward flow U and the radial shift r for each type of waveform I 1 , I 2 , I 3 and I 4 , respectively. The exact numeric value of U 1 , U 2 , U 3 and r 1 , r 2 , r 3 are shown in Table 7.  Table 9 Summary of the RBC morphology transition sequences recorded at different time instances t T under I 1 waveform, and different flow velocities U 1 , U 2 , U 3 and radial shift r 1 , r 2 , r 3 . Here, the time instances represent the first time the RBC deformed shape appeared. The acronyms C, S, CM , M , T , RS, EM and RD represent the croissant, slipper, complex multilobes, multilobes, trilobes, rolling stomatocytes, elongated multilobes and rolling discocyte, respectively. The exact numeric value of U 1 , U 2 , U 3 and r 1 , r 2 , r 3 are shown in Table 7 Table 10 Summary of the RBC morphology transition sequences recorded at different time instances t T under I 2 waveform, and the different flow velocities U 1 , U 2 , U 3 and initial placements r 1 , r 2 , r 3 . Here, the time instances represent the first time the RBC deformed shape appeared. The acronyms C, S, CM , EM and RD represent the croissant, slipper, complex multilobes, elongated multilobes and rolling discocyte, respectively. The exact numeric value of U 1 , U 2 , U 3 and r 1 , r 2 , r 3 are shown in Table 7 Table 11 Summary of the RBC morphology transition sequences recorded at different time instances t T under I 3 waveform, and the different flow velocities U 1 , U 2 , U 3 and initial placements r 1 , r 2 , r 3 . Here, the time instances represent the first time the RBC deformed shape appeared. The acronyms C, S, CM , T , EM and RD represent the croissant, slipper, complex multilobes, trilobes, elongated multilobes and rolling discocyte, respectively. The exact numeric value of U 1 , U 2 , U 3 and r 1 , r 2 , r 3 are shown in Table 7 Table 12 Summary of the RBC morphology transition sequences recorded at different time instances t T under I 4 waveform, and the different flow velocities U 1 , U 2 , U 3 and initial placements r 1 , r 2 , r 3 . Here, the time instances represent the first time the RBC deformed shape appeared. The acronyms C, S, CM , EM , RD and HX represent the croissant, slipper, complex multilobes, elongated multilobes, rolling discocyte and hexalobes, respectively. The exact numeric value of U 1 , U 2 , U 3 and r 1 , r 2 , r 3 are shown in Table 7. The dashed line shows that the RBC is placed at the channel's center line. The solid line depicts how the cell is transversely shifted from the cross-sectional center along the bisector of the first quadrant by a radial shift (r) in the y − z plane (Table 7).   The stable shapes are attained under the impact of a constant shear rate (I 0 ) in: (a) croissant shape (I 0 U 3 r 1 χ 1 ) and (b) slipper shape (I 0 U 4 r 3 χ 1 ). The RBC membrane only exhibits the tank-treading effect in the slipper shape (I 0 U 4 r 3 χ 1 ), which is characterized by the motions of two Lagrangian markers V 1 and V 2 . The slipper shape is maintained by the counter-clockwise rotation (the green arrow) of the cellular membrane around the RBC's centroid. The multilobe shape appears (c) under the oscillatory flow (Isψ 6 r 1 χ 3 ) during the backward phase (0.7T ).
!"# !$# Fig. 7 The extracellular flow patterns for: (a) the croissant shape (I 0 U 3 r 1 χ 1 ) and (b) the slipper shape (I 0 U 4 r 3 χ 1 ). The flow streamlines are reconstructed using the co-moving frame method as discussed in section 3.1.2. The tank-treading effect induces a closed vortex to form on the upstream side of the RBC.
(d) Fig. 8 The emergence of complex shapes induced by different inlet sinusoidal waveforms at different time instances in the flow cycle t 1 (end of forward phase), t 2 (during the resting period) and t 3 (backward flow phase). (a) I 1 U 1 r 1 χ 1 , (b) I 3 U 2 r 1 χ 1 , (c) I 4 U 2 r 1 χ 1 and (d) I 4 U 3 r 3 χ 1 .
(e) (f) r! r" Fig. 9 The impacts of the radial shift r on the time evolution of the RBC's centroid displacement (∆xc, ∆yc, ∆zc). The instantaneous evolution of the RBC's centroid position (xc, yc, zc)(t) is recorded as the RBC propagates along the channel. The displacements of the RBC from its initial location along three directions (∆xc, ∆yc, ∆zc) are measured in units of the length scale Ls. The evolution of the centroid position is examined under two conditions: (i) centred initial position r = r 1 = 0 (left column-(a − c) for the cases I 1 U 1 r 1 χ 1 , I 2 U 1 r 1 χ 1 , I 3 U 1 r 1 χ 1 , and I 4 U 1 r 1 χ 1 ; and (ii) off-centered initial position r = r 2 = 0.4 µm (right column -(d − f ) for the cases I 1 U 2 r 2 χ 1 , I 2 U 2 r 2 χ 1 , I 3 U 2 r 2 χ 1 , and I 4 U 2 r 2 χ 1 as described in Table 7. ); and (ii) the asymmetric I 4 (right column -(d − f )) waveforms at different values of the radial shift r 1 , r 2 , and r 3 . The symmetric flow cases (left column) include I 1 U 1 r 1 χ 1 , I 1 U 1 r 2 χ 1 , and I 1 U 1 r 3 χ 1 .
The asymmetric flow cases (right column) include I 4 U 1 r 1 χ 1 , I 4 U 1 r 2 χ 1 , and I 4 U 1 r 3 χ 1 cases. The exact values of r 1 , r 2 , and r 3 are described in Table 7.