# Loading and relaxation dynamics
of a red blood cell †

Fabio Guglietta
Department of Physics & INFN, University of Rome “Tor Vergata”
Chair for Computational Analysis of Technical Systems (CATS), RWTH Aachen
University
Computation-based Science and Technology Research Center, The Cyprus Institute
<EMAIL_ADDRESS>
&Marek Behr
Chair for Computational Analysis of Technical Systems (CATS), RWTH Aachen
University†
&Giacomo Falcucci
Department of Enterprise Engineering “Mario Lucertini,” University of Rome
“Tor Vergata"
Mauro Sbragaglia
Department of Physics & INFN, University of Rome “Tor Vergata”∗
Via della Ricerca Scientifica 1, 00133, Rome, Italy52056 Aachen, Germany20
Konstantinou Kavafi Str., 2121 Nicosia, CyprusVia del Politecnico 1, 00133,
Rome, Italy

###### Abstract

We use mesoscale numerical simulations to investigate the unsteady dynamics of
a single red blood cell (RBC) subjected to an external mechanical load. We
carry out a detailed comparison between the loading (L) dynamics, following
the imposition of the mechanical load on the RBC at rest, and the relaxation
(R) dynamics, allowing the RBC to relax to its original shape after the sudden
arrest of the mechanical load. Such a comparison is carried out by analyzing
the characteristic times of the two corresponding dynamics, i.e., $t_{L}$ and
$t_{R}$. When the intensity of the mechanical load is small enough, the two
kinds of dynamics are symmetrical ($t_{L}\approx t_{R}$) and independent of
the typology of mechanical load (intrinsic dynamics); otherwise, in marked
contrast, an asymmetry is found, wherein the loading dynamics is typically
faster than the relaxation one. This asymmetry manifests itself with non-
universal characteristics, e.g., dependency on the applied load and/or on the
viscoelastic properties of the RBC membrane. To deepen such a non-universal
behaviour, we consider the viscosity of the erythrocyte membrane as a variable
parameter and focus on three different typologies of mechanical load
(mechanical stretching, shear flow, elongational flow): this allows to clarify
how non-universality builds up in terms of the deformation and rotational
contributions induced by the mechanical load on the membrane. Finally, we also
investigate the effect of the elastic shear modulus on the characteristic
times $t_{L}$ and $t_{R}$. Our results provide crucial and quantitative
information on the unsteady dynamics of RBC and its membrane response to the
imposition/cessation of external mechanical loads.

00footnotetext: † Electronic Supplementary Information (ESI) available: one
PDF containing further details regarding the simulations performed; four
videos showing the simulations performed. See DOI: 10.1039/D1SM00246E

##  1 Introduction

Red blood cells (RBCs) are biological cells made of a viscoelastic membrane
enclosing a viscous fluid (cytoplasm): their main features are the biconcave
shape and the absence of a nucleus and most organelles, that allow them to
carry oxygen even inside the smallest capillaries [1, 2, 3]. In fact, during
circulation, RBCs deform multiple times, rearranging their shape to adapt to
the physiological conditions of the blood flow. The mechanical properties of
RBC’s membrane have been deeply investigated, both numerically [3, 4, 5, 6, 7,
8, 9, 10, 11, 12, 13, 14] and experimentally [15, 16, 17, 18, 19, 20, 21, 22,
23]. Research interest on the mechanical response of RBC’s membrane was
prompted by several reasons: among the others, the link between its properties
and the erythrocyte’s health conditions [11, 16], or the role played by the
membrane dynamics in the design of biomedical devices[24, 25, 26, 27]. A huge
effort has been devoted to the characterisation of the time-independent
properties of the membrane, for the study of the corresponding steady-state
configurations. In recent years, also the dynamical behaviour of RBCs has been
investigated in several works, both numerically [28, 29, 30, 31, 8, 32, 33,
34] and experimentally [17, 19, 20, 35]. When dealing with time-dependent
properties of biological cells (or capsules, in general), the membrane
viscosity plays a crucial role [36, 37, 38, 33, 34, 39, 40, 41, 42, 43]. Evans
[44] showed that the RBC relaxation time is affected by both the membrane
viscosity and the dissipation in the adjacent aqueous phases (i.e., cytoplasm
and external solution); neglecting the membrane viscosity, i.e.,
$\mu_{\mbox{m}}=0$, he predicted a relaxation time $t_{\mbox{R}}\approx
1\times 10^{-3}$ s (also confirmed by numerical simulations [33]), a
remarkably lower value compared to other works in the literature [22, 19, 17,
45, 33]. Several works aimed at the precise estimation of the value of the
membrane viscosity [46, 18, 45, 47, 22, 48, 49, 17, 19, 32], finding that
$\mu_{m}$ roughly ranges between $10^{-7}\mbox{ m Pa s}$ and $10^{-6}\mbox{ m
Pa s}$. Such variability may be ascribed to different factors, e.g., the
different theoretical models used to infer $\mu_{\mbox{m}}$ [19, 46], the
different experimental apparatuses (such as micro-pipette aspiration [46, 18,
45, 17], microchannel deformation [49], or other setups [47, 22, 48]), etc. As
a matter of fact, although $\mu_{\mbox{m}}$ is an essential parameter to
quantitatively characterise the time dynamics of RBCs [34, 37], its precise
value has not been accurately determined so far, which warrants a parametric
investigation. Moreover, some earlier studies [50] proposed to use an
increased apparent viscosity ratio to account for the energy dissipation due
to the presence of a viscous membrane: even though this assumption provides a
qualitative description of the effects of the membrane viscosity, it does not
account for a quantitative characterisation, as shown by recent studies [51,
37, 34, 41, 40]. Our previous work [33] aimed to investigate the effect of
membrane viscosity $\mu_{\mbox{m}}$ on the relaxation dynamics of a single
RBC, and we found that increasing the value of $\mu_{\mbox{m}}$, as well as
increasing the intensity of the loading strength, leads to faster recovery
dynamics. Moreover, we simulated two experimental setups, i.e., the stretching
with optical tweezers and the deformation due to an imposed shear flow, and we
found a dependency on the kind of mechanical load when the strengths of load
are large enough.
The relaxation dynamics, however, gives only a partial characterisation of the
time-dependent response of RBCs to external forces: therefore, the loading
process should be considered as well, as already pointed out in earlier
literature papers. Chien et al. [18] experimentally studied both the loading
and the relaxation dynamics of RBC membrane through micro-pipette aspiration,
providing evidence that the two dynamics are not symmetrical in certain
conditions; however, a systematic study involving different stress values and
different typologies of mechanical loads was not performed. Diaz et al. [43]
studied the dynamics of a pure elastic capsule with a hyperelastic membrane
deformed by an elongational flow: they focused on both loading and relaxation,
finding an asymmetry. However, their model did not take into account the
membrane viscosity and the asymmetry was not studied for different typologies
of mechanical loads. Thus, although previous literature points to two distinct
dynamics for loading and relaxation [18, 43, 42], a comprehensive parametric
study on the effects of $\mu_{\mbox{m}}$ for different typologies of
mechanical loads and flow conditions (such as simple shear flow or
elongational flow) has never been attempted, so far. This paper aims at
filling this gap with the help of mesoscale numerical simulations. Indeed, for
this kind of characterisation, numerical simulations can be thought of as the
appropriate tool of analysis [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14], due to
the obvious experimental difficulties in carrying out such systematic
investigation [17, 18, 19, 21, 23]. We provide a quantitative characterisation
of loading and relaxation dynamics exploring three typologies of mechanical
loads. To do this, we built three different simulation setups: the stretching
simulation (STS), which simulates the deformation with optical tweezers [16]
(see Fig. 1, panel (a)); the shear simulation (SHS), i.e., the deformation in
simple shear flow (see Fig. 1, panel (c)); the four-roll mill simulation
(FRMS), where the deformation is induced by an elongational axisymmetric flow
made by the rotation of four cylinders (see Fig. 1, panel (e)). These three
numerical setups are chosen to inspect the different roles that the membrane
rotation and/or membrane deformation have in the time-dependent dynamics. This
information is summarised in Tab. 1. First, we systematically study the
characteristic times of both the loading ($t_{\mbox{L}}$) and the relaxation
($t_{\mbox{R}}$) processes and their ratio
$\tilde{t}=t_{\mbox{L}}/t_{\mbox{R}}$, as a function of the load strength and
membrane viscosity $\mu_{\mbox{m}}$. For small strengths, the two
characteristic times are essentially equal and set by the value of
$\mu_{\mbox{m}}$; however, for strengths large enough, the loading dynamics is
found to be faster than the relaxation dynamics, leading to a non-universal
behaviour while changing the typology of the mechanical load. Such non-
universal contributions are further characterised in terms of the importance
of rotation and deformation of the membrane, according to the different load
mechanisms. Some useful parametrizations for both $t_{\mbox{R}}$ and
$t_{\mbox{L}}$ as a function of the membrane viscosities are also provided.
Finally, since different pathologies that cause the reduction of RBC membrane
elasticity are known [11, 16, 52, 53, 54, 55, 56, 57, 58, 59], we also study
the dependency of both the loading ($t_{\mbox{L}}$) and the relaxation
($t_{\mbox{R}}$) times as well as their ratio $\tilde{t}$ on the elastic shear
modulus $k_{\mbox{S}}$, for a fixed value of membrane viscosity.

|
---|---
(a) STretching Simulation (STS) | (b) Stretching simulation (STS): deformation.
|
(c) SHear Simulation (SHS) | (d) Shear simulation (SHS): deformation.
|
(e) Four-Roll Mill Simulation (FRMS) | (f) Four-roll mill simulation (FRMS): deformation.
Figure 1: Loading-relaxation (L-R) simulations for red blood cell (RBC) at
changing the typology of mechanical load. Left panels: the three different L-R
simulations investigated in the paper are sketched: grey arrows refer to the
mechanical load, either an applied force $\mathbf{F}$ or an applied velocity
$\text{U}_{\mbox{w}}$, while the RBC membrane forces
($\mathbf{F}_{\mbox{mem}}$) are sketched with green arrows. In all
simulations, the deformation $D(t)$, i.e., the ratio between the difference
and the sum of the axial and transversal diameters (see Eq. (11)), is used to
fit the loading and relaxation times ($t_{\mbox{L}}$ and $t_{\mbox{R}}$,
respectively; see Eqs. (12) (13) (14)). Right panels: we report the
deformation $D(t)$ (see Eq. (11)) as a function of time for two values of
membrane viscosity $\mu_{\mbox{m}}$. Panels (a-b): we simulate the stretching
with optical tweezers [16] (STS), in which two forces with the same intensity
and opposite direction stretch the membrane in two areas at the ends of the
RBC (see Sec. 3.1); the deformation $D(t)$ is reported for $F=90\times
10^{-12}$ N. Panels (c-d): deformation induced by simple shear flow (SHS),
with $\text{U}_{\mbox{w}}=(\pm\dot{\gamma}/2,0,0)$, where $\dot{\gamma}$ is
the shear rate (see Sec. 3.2); the deformation $D(t)$ is reported for
$\dot{\gamma}=86\mbox{ s}^{-1}$. Panels (e-f): in the four-roll mill
simulation (FRMS) we simulate four rotating cylinders to reproduce an
elongational flow that deforms the membrane [60] (see Sec. 3.3); the
deformation $D(t)$ is reported for $\dot{\gamma}_{\mbox{FRMS}}=80\mbox{
s}^{-1}$. Four videos showing these simulations are available (see ESI†).

The paper is organised in the following way: in Sec. 2 we provide some details
on the numerical method used to simulate both the fluid and the membrane of
the RBC; in Sec. 3 we analyse the three simulated loading mechanisms (the
stretching simulation, Sec. 3.1; the shear simulation, Sec. 3.2; the four-roll
mill simulation, Sec. 3.3); a detailed discussion section with comparisons
between the loading mechanisms is provided in Sec. 4; finally, conclusions are
reported in Sec. 5.

Load Type | Rotation | Direct Forcing
---|---|---
STS | NO | YES
SHS | YES | NO
FRMS | NO | NO
Table 1: Summary of the main characteristics of the three kinds of applied
mechanical loads (see Fig. 1). For each load type, we specify if the rotation
is induced on the membrane while loading and if the forcing is directly
applied on the nodes of the mesh used to discretise the membrane (otherwise,
the membrane is forced indirectly via hydrodynamic flow).

##  2 Numerical method

We perform three-dimensional numerical simulations in the framework of the
Immersed Boundary – Lattice Boltzmann method (IB–LBM) [61, 62]. The
methodology, as well as the membrane model, are the same already used and
validated in [33]: here we report an essential summary.
The equation of motion for a fluid with viscosity $\mu$ is given by the
Navier-Stokes (NS) equations:

$\rho\left(\frac{\partial\mathbf{u}}{\partial
t}+(\mathbf{u}\cdot\boldsymbol{\nabla})\mathbf{u}\right)=-\boldsymbol{\nabla}p+\mu\boldsymbol{\nabla}^{2}\mathbf{u}+\mathbf{F}\;,$
(1)

where $\rho$ and $\mathbf{u}$ are the density and the velocity of the fluid,
respectively; $p$ is the isotropic pressure; $\mathbf{F}$ is an external body
force density. If the fluid is incompressible (as in the present work) the
condition $\boldsymbol{\nabla}\cdot\mathbf{u}=0$ holds.
In the LBM, instead of directly solving the NS equations by integrating Eq.
(1), the fluid is represented by the so-called populations
$f_{i}(\mathbf{x},t)$, that stand for the density of fluid molecules moving
with velocity $\mathbf{c}_{i}$ at position $\mathbf{x}$ and time $t$. The
populations evolve according to the Lattice Boltzmann equation:

$\begin{split}f_{i}(\mathbf{x}+\mathbf{c}_{i}\Delta t,t+\Delta
t)-f_{i}(\mathbf{x},t)=\\\ =-\frac{\Delta
t}{\tau}\left(f_{i}(\mathbf{x},t)-f_{i}^{(\mbox{eq})}(\mathbf{x},t)\right)&+f_{i}^{(F)}\;,\end{split}$
(2)

in which $\Delta t$ is the discrete time step, $\tau$ is the relaxation time,
$f_{i}^{(F)}$ is the source term that takes into account the force density (it
has been implemented according to the “Guo” scheme [63]), and
$f_{i}^{(\mbox{eq})}$ is the equilibrium distribution function (we refer back
to [61, 62] for the details). The fluid density $\rho$ and the velocity
$\mathbf{u}$ are given by:

$\begin{split}\rho(\mathbf{x},t)=\sum_{i}f_{i}(\mathbf{x},t)\;,\\\
\rho\mathbf{u}(\mathbf{x},t)=\sum_{i}\mathbf{c}_{i}f_{i}(\mathbf{x},t)\;.\end{split}$
(3)

The link between NS and LB equations (Eq. (1) and Eq. (2), respectively) is
given by the following relation:

$\mu=\rho c_{\mbox{s}}^{2}\left(\tau-\frac{\Delta t}{2}\right)\;,$ (4)

where $c_{\mbox{s}}=\Delta x/\Delta t\sqrt{3}$ is the speed of sound. In the
following, we considered both the lattice spacing $\Delta x$ and the time
interval $\Delta t$ equal to 1.
We simulate two fluids: one outside the membrane (the plasma, with viscosity
$\mu_{\mbox{out}}=1.2\times 10^{-3}$ Pa s) and one inside it (the cytosol,
with viscosity $\mu_{\mbox{in}}=6\times 10^{-3}$ Pa s). The viscosity ratio is
given by

$\lambda=\frac{\mu_{\mbox{in}}}{\mu_{\mbox{out}}}\;,$ (5)

providing $\lambda=5$. We implement the parallel Hoshen-Kopelman algorithm to
recognise which lattice sites are inside or outside the membrane (see [64] for
details).
The RBC membrane is described as a 3D triangular mesh of $\approx~{}4000$
elements, whose shape at rest is given by [65]

$\begin{split}z(x,y)=\pm\sqrt{1-\frac{x^{2}+y^{2}}{r^{2}}}\cdot\\\
\cdot\left(C_{0}+C_{1}\frac{x^{2}+y^{2}}{r^{2}}+C_{2}\left(\frac{x^{2}+y^{2}}{r^{2}}\right)^{2}\right)\;,\end{split}$
(6)

with $C_{0}=0.81\times 10^{-6}\mbox{ m}$, $C_{1}=7.83\times 10^{-6}\mbox{ m}$
and $C_{2}=-4.39\times 10^{-6}\mbox{ m}$; $r=3.91\times 10^{-6}\mbox{ m}$ is
the large radius.
The membrane is characterised by a resistance to shear deformation, area
dilation and bending; the viscoelastic behaviour is implemented, as well. The
first two terms form the strain energy $W_{\mbox{S}}$ are described by Skalak
model [66]:

$W_{\mbox{S}}=\sum_{j}A_{j}\left[\frac{k_{\mbox{S}}}{12}\left(I_{1}^{2}+2I_{1}-2I_{2}\right)+\frac{k_{\alpha}}{12}I_{2}^{2}\right]\;,$
(7)

where $k_{\mbox{S}}=5.3\times 10^{-6}\mbox{ N m}^{-1}$ [16] and
$k_{\alpha}=50k_{\mbox{S}}$ [3] are the surface elastic shear modulus and the
area dilation modulus, respectively; $I_{1}=\lambda_{1}^{2}+\lambda_{2}^{2}-2$
and $I_{2}=\lambda_{1}^{2}\lambda_{2}^{2}-1$ are the strain invariants for the
$j$-th element, while $\lambda_{1}$ and $\lambda_{2}$ are the principal
stretch ratios [66, 3]; $A_{j}$ is the surface are of the $j$-th element. We
adopt the Helfrich formulation to compute the free-energy $W_{\mbox{B}}$
related to the resistance to bending [67]. Following [3], we discretise the
bending energy as:

$W_{\mbox{B}}=\frac{k_{\mbox{B}}\sqrt{3}}{2}\sum_{\langle
i,j\rangle}\left(\theta_{ij}-\theta_{ij}^{(0)}\right)^{2}\;,$ (8)

where $k_{\mbox{B}}=2\times 10^{-19}\mbox{ N m}$ [68] is the bending modulus;
the sum runs over all the neighbouring triangular elements, and $\theta_{ij}$
is the angle between the normals of the $i$-th and $j$-th elements
($\theta_{ij}^{(0)}$ is the same angle in the unperturbed configuration). Once
we have the total free-energy $W=W_{\mbox{S}}+W_{\mbox{B}}$, we compute the
force acting on the $i-$th node by performing the derivative of $W$ with
respect to the coordinates of the node $\mathbf{x}_{i}$:

$\mathbf{F}_{i}=-\frac{\partial W(\mathbf{x}_{i})}{\partial\mathbf{x}_{i}}\;.$
(9)

Note that we are implementing neither area nor volume conservation: in fact,
as stated in [33], we checked that both area and volume were conserved, even
without an explicit area or volume conservation law (see Electronic
Supplementary Information in [33]).
Regarding the viscoelastic term, we implement the Standard Linear Solid (SLS)
model [69, 36]. The viscous stress tensor is given by

$\boldsymbol{\tau}^{\nu}=\mu_{\mbox{s}}\left(2\dot{\mathbf{E}}-\mbox{tr}(\dot{\mathbf{E}})\mathbb{1}\right)+\mu_{\mbox{d}}\mbox{tr}(\dot{\mathbf{E}})\mathbb{1}\;,$
(10)

where $\mathbf{E}$ is the strain tensor (see [69, 33]); $\mu_{\mbox{s}}$ and
$\mu_{\mbox{d}}$ are the shear and dilatational membrane viscosity: in this
work, we assume $\mu_{\mbox{s}}~{}=~{}\mu_{\mbox{d}}~{}=~{}\mu_{\mbox{m}}$
[33]. We refer to [69, 33] for the computation of the stress tensor
$\boldsymbol{\tau}^{\nu}$ (Eq. (10)) as well as for the nodal force
$\mathbf{F}_{i}$ (Eq. (9)).
Finally, once we have the nodal force $\mathbf{F}_{i}$ for each node $i$ of
the 3D mesh, we spread this force to the lattice nodes via the IBM (see [61]
for details). We adopt the same scheme as in [33].

Figure 2: Characteristic times $t_{\mbox{L}}$ (first column of panels) and
$t_{\mbox{R}}$ (second column of panels), as well as the ratio
$\tilde{t}=t_{\mbox{L}}/t_{\mbox{R}}$ (third column of panels) are reported
for the three simulations performed, i.e., stretching simulation (STS,  ,
panels (a-c), Sec. 3.1), shear simulation (SHS,  , panels (d-f), Sec. 3.2),
four-roll mill simulation (FRMS,  , panels (g-i), Sec. 3.3), for different
values of membrane viscosity $\mu_{\mbox{m}}$ (from lightest to darkest
color): $\mu_{\mbox{m}}=0$ m Pa s ( ,  ,  ), $\mu_{\mbox{m}}=0.64\times
10^{-7}$ m Pa s ( ,  ,  ), $\mu_{\mbox{m}}=1.59\times 10^{-7}$ m Pa s ( ,  ,
), $\mu_{\mbox{m}}=3.18\times 10^{-7}$ m Pa s ( ,  ,  ). The red dashed line
represents the reference value for the symmetric case, i.e., $\tilde{t}=1$.

##  3 Loading and relaxation time

In this section, we quantitatively study the loading time $t_{\mbox{L}}$ and
the relaxation time $t_{\mbox{R}}$ with three different loading mechanisms
(see Fig. 1): the stretching with optical tweezers (STS, see Sec. 3.1), the
deformation in simple shear flow (SHS, see Sec. 3.2) and the deformation in an
elongational flow (FRMS, see Sec. 3.3). These three simulations differ mainly
for two aspects (summarised in Tab. 1): the first one is that the membrane can
be deformed by an external force that acts directly on the membrane (like in
the STS) or by the viscous friction with the fluid (SHS and FRMS); moreover,
the membrane can rotate (like in the SHS) or not (STS and FRMS). The idea
underlying the choice of these three different setups is to catch which of the
aforementioned characteristics affects the loading and relaxation dynamics.
To quantify the loading time $t_{\mbox{L}}$ and the relaxation time
$t_{\mbox{R}}$, we define the deformation parameter

$D(t)=\frac{d_{\mbox{A}}(t)-d_{\mbox{T}}(t)}{d_{\mbox{A}}(t)+d_{\mbox{T}}(t)}\;,$
(11)

where $d_{\mbox{A}}$ and $d_{\mbox{T}}$ represent the length of the axial and
transversal diameters, i.e., the greater and medium eigenvalues of the inertia
tensor (see [33]). In our computational domain, $d_{\mbox{A}}$ and
$d_{\mbox{T}}$ lie in the $x-y$ plane. We define the average deformation
$D_{\mbox{av}}$, i.e., the value of the deformation $D$ such that
$\lim_{t\to\infty}D(t)=D_{\mbox{av}}$ in the loading simulation and
$D(0)=D_{\mbox{av}}$ in the relaxation simulation.
Qualitatively, the loading time $t_{\mbox{L}}$ is the characteristic time the
deformation $D(t)$ takes to reach $D_{\mbox{av}}$; the relaxation time
$t_{\mbox{R}}$ is the characteristic time needed to relax to the initial
shape, after the arrest of the mechanical load. Quantitatively, we can get
$t_{\mbox{L}}$ and $t_{\mbox{R}}$ via a fit of $D(t)/D_{\mbox{av}}$ with the
following functions:

$L_{\mbox{1}}(t)=1-\exp\left\\{-\left(\frac{t}{t_{\mbox{L}}}\right)^{\delta_{\mbox{L}}}\right\\}\;,$
(12)
$L_{\mbox{2}}(t)=1-\exp\left\\{-\left(\frac{t}{t_{\mbox{L}}}\right)^{\delta_{\mbox{L}}}\right\\}\
\cos\left(\frac{t}{t_{\mbox{L}}^{\mbox{cos}}}\right)\;,$ (13)
$R(t)=\exp\left\\{-\left(\frac{t}{t_{\mbox{R}}}\right)^{\delta_{\mbox{R}}}\right\\}\;,$
(14)

where $L_{\mbox{1}}$ is used to fit the loading time for the STS and the FRMS
(see Sec. 3.1-3.3, respectively); $L_{\mbox{2}}$ is used to fit the loading
time for the SHS (see Sec. 3.2); $R$ is used to fit the relaxation time
$t_{\mbox{R}}$ for all three simulations; $\delta_{\mbox{L}}$ and
$\delta_{\mbox{R}}$ are parameters introduced to improve the fit [32] (see
[33] for some more details) and will be characterised in Sec. 4.
Note that we propose two different functions to fit data during the loading
(i.e., Eq. (12) and Eq. (13)); this is due to the different deformation
process of the RBC: in the STS and FRMS, $D(t)$ is a monotonic increasing
function with an asymptote in $D=D_{\mbox{av}}$ (see Fig. 1, panels (b) and
(f)); in the SHS, $D(t)$ oscillates around $D_{\mbox{av}}$, and the amplitude
of the oscillations varies with the value of the membrane viscosity
$\mu_{\mbox{m}}$ (see Fig. 1, panel (d)). These oscillations have been also
observed for viscoelastic capsules[38, 69, 33, 70]. Before starting the
relaxation, for the STS and FRMS, we waited long enough to achieve the steady
value of the deformation $D_{\mbox{av}}$; for the SHS, we ensured the
oscillations were small enough if compared to $D_{\mbox{av}}$. Notice,
however, that depending on the importance of viscous effects with respect to
elastic ones, there may be cases where such oscillations are damped on very
long times [70, 71, 38].
Moreover, while the values of the external mechanical loads we simulated are
comparable to each other in terms of stress (see Sec. 4), the values of
$D_{\mbox{av}}$ for the SHS are much smaller than for the STS and FRMS (see
Fig. 4, panel (c)).
Since we want to compare the loading time $t_{\mbox{L}}$ and the relaxation
time $t_{\mbox{R}}$, we define the ratio

$\tilde{t}=\frac{t_{\mbox{L}}}{t_{\mbox{R}}}.$ (15)

In all the following simulations, the membrane viscosity ranges between
$\mu_{\mbox{m}}\in[0,3.18]\times 10^{-7}\mbox{ m Pa s}$ [46, 18, 45, 47, 22,
48, 49, 17, 19, 32]. In the STS, the applied force is in the range
$F\in[5,70]\times 10^{-12}\mbox{ N}$; in the SHS, we simulated shear rates in
the range $\dot{\gamma}\in[1.23,123]\mbox{ s}^{-1}$; finally, in the FRMS, we
simulated $\dot{\gamma}_{\mbox{FRMS}}\in[1,120]\mbox{ s}^{-1}$ (see Eq. (16)).
Values for all parameters used in the simulations in both physical and lattice
units are reported in Tab. 1 in ESI†.

###  3.1 Stretching simulation (STS)

In order to simulate the stretching with optical tweezers, we apply two forces
with the same intensity $F$ and opposite directions at the ends of the RBC
(see Fig. 1, panel (a)). Simulations are performed in a 3D box
$L_{\mbox{x}}\times L_{\mbox{y}}\times L_{\mbox{z}}=(28,12,12)\times 10^{-6}$
m. In Fig. 2, we report the loading time $t_{\mbox{L}}$ (panel (a)) and the
relaxation time $t_{\mbox{R}}$ (panel (b)) as a function of $F$, for different
values of membrane viscosity $\mu_{\mbox{m}}$. In both cases, increasing the
loading strength (as well as decreasing the value of membrane viscosity
$\mu_{\mbox{m}}$) results in faster dynamics. It is interesting to compare
$t_{\mbox{L}}$ and $t_{\mbox{R}}$: in Fig. 2, panel (c), we report the ratio
$\tilde{t}$ (see Eq. (15)).

###  3.2 Shear simulation (SHS)

In the shear simulation, the deformation is due to a linear shear flow with
intensity $\dot{\gamma}$. We set the wall velocity
$\text{U}_{\mbox{w}}=(\pm\dot{\gamma}L_{\mbox{z}}/2,0,0)$ on the two plane
walls at $z=\pm L_{\mbox{z}}/2$, and the RBC is oriented as reported in Fig.
1, panel (c). This choice is surely the optimal one for the purpose of the
present study, since we can focus on the relaxation/loading process without
any further complication. In real-world experiments, indeed, RBCs do not
necessarily orient in the shear plane and evolve into a complex dynamics with
many dynamical modes [72, 29]. In particular, for the values of shear rate
$\dot{\gamma}$ we are interested in, a rolling dynamics is expected if the RBC
is not oriented in the shear plane [73]. Further complications can be
introduced by polydispersity, i.e., that RBCs may have different sizes and
viscoelastic properties [16, 74]. We have preferred to avoid all these extra
complications which could not allow us for a fair assessment in the comparison
between a pure relaxation dynamics and a pure loading dynamics. Simulations
are performed in a 3D box $L_{\mbox{x}}\times L_{\mbox{y}}\times
L_{\mbox{z}}=(20,32,20)\times 10^{-6}$ m. In Fig. 2, the loading time
$t_{\mbox{L}}$ (panel (d)) and the relaxation time $t_{\mbox{R}}$ (panel (e))
as a function of $\dot{\gamma}$ for different values of membrane viscosity
$\mu_{\mbox{m}}$ are reported, as well as the ratio $\tilde{t}$ (panel (f)).
While the relaxation time $t_{\mbox{R}}$ shows a similar behaviour compared to
the STS (see Fig. 2, panel (b)), a few more words are needed regarding the
loading time $t_{\mbox{L}}$. Unlike the STS, now we have two characteristic
times, that are $t_{\mbox{L}}$ and $t_{\mbox{L}}^{\mbox{cos}}$ (see Eq. (13)):
$t_{\mbox{L}}$ measures the time the membrane takes to reach the average
deformation $D_{\mbox{av}}$, while $t_{\mbox{L}}^{\mbox{cos}}$ measures the
period of the oscillations. Data for $t_{\mbox{L}}^{\mbox{cos}}$ are reported
in ESI†, Fig. 1. In contrast to the STS, the loading time $t_{\mbox{L}}$ first
decreases, and then it does not change much with the intensity of the
mechanical load.

Figure 3: Comparison between the characteristic times $t_{\mbox{L}}$ (panel
(a)), $t_{\mbox{R}}$ (panel (b)) and $\tilde{t}=t_{\mbox{L}}/t_{\mbox{R}}$
(panel (c)) as a function of the stress $\sigma$ (see Sec. 4) for the three
simulations performed, i.e., stretching simulation (STS,  , Sec. 3.1), shear
simulation (SHS,  , Sec. 3.2), four-roll mill simulation (FRMS,  , Sec. 3.3),
for different values of membrane viscosity $\mu_{\mbox{m}}$ (from lightest to
darkest color): $\mu_{\mbox{m}}=0$ m Pa s ( ,  ,  ),
$\mu_{\mbox{m}}=0.64\times 10^{-7}$ m Pa s ( ,  ,  ),
$\mu_{\mbox{m}}=1.59\times 10^{-7}$ m Pa s ( ,  ,  ),
$\mu_{\mbox{m}}=3.18\times 10^{-7}$ m Pa s ( ,  ,  ). The red dashed line
represents the reference value for the symmetric case, i.e., $\tilde{t}=1$.

Figure 4: Panel (a): snapshots of RBC at selected times. A point is selected
on the membrane (green sphere) to perform a Lagrangian tracking and determine
the time dependency of the angle $\phi$ that the point direction forms with
the y axis in the deformation plane. Panel (b): we report the angle $\phi$
(see panel (a) above) as a function of time for the stretching simulation
(STS,  , Sec. 3.1), shear simulation (SHS,  , Sec. 3.2), four-roll mill
simulation (FRMS,  , Sec. 3.3), for $\mu_{\mbox{m}}=3.18\times 10^{-7}$ m Pa
s. The red and blue shades represent loading and relaxation regions,
respectively. Panel (c): average deformation $D_{\mbox{av}}$ (see text for
details) as a function of the stress $\sigma$ for the three simulations
performed (STS, SHS and FRMS). SHS data are displayed for different values of
membrane viscosity $\mu_{\mbox{m}}$ (from lightest to darkest color):
$\mu_{\mbox{m}}=0$ m Pa s ( ), $\mu_{\mbox{m}}=0.64\times 10^{-7}$ m Pa s ( ),
$\mu_{\mbox{m}}=1.59\times 10^{-7}$ m Pa s ( ), $\mu_{\mbox{m}}=3.18\times
10^{-7}$ m Pa s ( ). STS data ( ) and FRMS data ( ) data are only reported for
$\mu_{\mbox{m}}=3.18\times 10^{-7}$ m Pa s.

###  3.3 Four-roll mill simulation (FRMS)

In this case, we simulate the effect of four cylinders rotating [60], as shown
in Fig. 1, panel (e), in order to create a flow similar to a pure elongational
one. Simulations are performed in a 3D box $L_{\mbox{x}}\times
L_{\mbox{y}}\times L_{\mbox{z}}=(48,48,20)\times 10^{-6}$ m. The idea is to
simulate a loading mechanism that is a mixture of stretching with optical
tweezers and deformation in simple shear flow (see Tab. 1): in fact, in this
case, the membrane does not rotate (like in the STS) and the deformation is
caused by the flow (like in the SHS).
To create such a flow, we impose a force density [60]

$\mathbf{F}(x,y)=2k\mu\dot{\gamma}_{\mbox{FRMS}}\left(\begin{array}[]{c}\sin(kx)\cos(ky)\\\
-\cos(kx)\sin(ky)\\\ 0\end{array}\right)\;,$ (16)

where $k=2\pi/L_{\mbox{x}}$, $\mu$ is the local fluid viscosity, and
$\dot{\gamma}_{\mbox{FRMS}}$ is used to tune the load strength. We multiplied
Eq. (16) by $k$ to make the velocity gradient independent of the size of the
fluid domain111The following result is valid in a homogeneous fluid with
dynamics viscosity $\mu$.:

$\frac{\partial\mathbf{u}}{\partial\mathbf{x}}=\dot{\gamma}_{\mbox{FRMS}}\left(\begin{matrix}\cos(kx)\cos(ky)&-\sin(kx)\sin(ky)\\\
\sin(kx)\sin(ky)&-\cos(kx)\cos(ky)\end{matrix}\right)\;,$ (17)

where we have reported only $x$ and $y$ components, i.e., the components in
the plane of the shear. Note that Eq. (16) gives a pure elongational flow only
in $x=\pi/2,3\pi/2$ and in $y=\pi/2,3\pi/2$.
In Fig. 2 we report the loading time $t_{\mbox{L}}$ (panel (g)) and the
relaxation time $t_{\mbox{R}}$ (panel (h)) as a function of
$\dot{\gamma}_{\mbox{FRMS}}$. As for the STS and the SHS, both $t_{\mbox{L}}$
and $t_{\mbox{R}}$ decrease when the loading force increases or the membrane
viscosity $\mu_{\mbox{m}}$ decreases. In Fig. 2, panel (i), the ratio
$\tilde{t}$ is reported.

##  4 Discussion

In our simulations, the intensity of the three kinds of mechanical loads is
changed by varying different quantities, i.e., $F$ for the STS, $\dot{\gamma}$
for the SHS and $\dot{\gamma}_{\mbox{FRMS}}$ for the FRMS. To facilitate a
comparison between them, we first consider the characteristic times
$t_{\mbox{L}}$ and $t_{\mbox{R}}$ as well as the ratio $\tilde{t}$ as a
function of the characteristic simulation stress $\sigma$ (Fig. 3). To
evaluate the stress $\sigma$ for the STS, we computed the area $A$ at the end
of the RBC where the force $F$ is applied. Then, the stress is given by
$\sigma^{\mbox{STS}}=F/A$; for the SHS, we wrote the stress as
$\sigma^{\mbox{SHS}}=2\dot{\gamma}\mu_{\mbox{out}}$ [3]. Finally, for the
FRMS, the stress is given by the stress-peak
$\sigma^{\mbox{FRMS}}=\mu_{\mbox{out}}\dot{\gamma}_{\mbox{FRMS}}$. In all
three simulations, the loading and relaxation times ($t_{\mbox{L}}$ and
$t_{\mbox{R}}$, respectively) show qualitatively the same behaviour, i.e.,
they decrease when the loading strength increases or when the membrane
viscosity $\mu_{\mbox{m}}$ decreases (see Fig. 3, panels (a) and (b)); the
ratio $\tilde{t}=t_{\mbox{L}}/t_{\mbox{R}}$ is reported in Fig. 3, panel (c).
For small forces ($\sigma\to 0$) we observe a clear tendency towards symmetry
between loading and relaxation ($\tilde{t}(\sigma,\mu_{\mbox{m}})\rightarrow
1$), meaning that the characteristic times $t_{\mbox{L}}$ and $t_{\mbox{R}}$
tend to be equal. This is the limit where one expects to recover the intrinsic
dynamics of the membrane, which depends only on the value of membrane
viscosity $\mu_{\mbox{m}}$ [33, 19].
On the other hand, for force strengths large enough, loading and relaxation
dynamics are asymmetrical, i.e., $\tilde{t}\neq 1$. As already noticed
elsewhere [42], this asymmetry could be explained by energetic considerations:
in fact, during the loading phase, the deformation is driven by the external
load (i.e., an external source of energy), while during the relaxation, the
membrane provides the whole energy. Beyond these qualitative considerations,
results in Fig. 3 provide a systematic characterisation of the relaxation
times, as a function of either the stress $\sigma$ or the membrane viscosity
$\mu_{\mbox{m}}$: an important message conveyed by our analysis is that the
asymmetry is not universal, i.e., on equal values of membrane viscosity
$\mu_{\mbox{m}}$, the ratio $\tilde{t}$ depends on the kind of mechanical
load. Just to give some numbers, the difference between the values of
$\tilde{t}$ for the STS and FRMS is roughly constant ($\approx$ 30%) and it
goes to zero for small values of $\sigma$; for the SHS the situation is a bit
more complex because $\tilde{t}$ depends on $\mu_{\mbox{m}}$. However, if we
compare SHS against STS for $\mu_{\mbox{m}}=3.18\times 10^{-7}$ m Pa s, we
find a difference of less then 30% for small values of $\sigma$ (i.e.,
$\sigma<0.1$ Pa), while such a difference goes over the 50% for large values
of $\sigma$ (i.e., $\sigma>0.1$ Pa).
If we think that the asymmetry comes from the presence of a mechanical load
with load strength large enough [43, 42], it comes natural to expect a non-
universality and a dependency on the details of the loading mechanism. Thanks
to our analysis, we are in a condition to further characterise this non-
universality: indeed, we observe that while $\tilde{t}$ does not depend on
$\mu_{\mbox{m}}$ for the STS and FRMS
($\tilde{t}^{\mbox{{STS}}}=\tilde{t}^{\mbox{{STS}}}(\sigma)$,
$\tilde{t}^{\mbox{{FRMS}}}=\tilde{t}^{\mbox{{FRMS}}}(\sigma)$), it actually
does in the SHS
($\tilde{t}^{\mbox{{SHS}}}=\tilde{t}^{\mbox{{SHS}}}(\sigma,\mu_{\mbox{m}})$).
The collapse shown by $\tilde{t}$ in the STS and FRMS (see Fig. 3, panel (c))
suggests a factorisation of the loading and relaxation times in two
contributions: one depending on the membrane viscosity $\mu_{\mbox{m}}$ and
one on the load intensity $\sigma$:

$t_{\mbox{L}}^{K}(\sigma,\mu_{\mbox{m}})\approx
t_{\mbox{L}}^{*K}(\sigma)t_{0}^{K}(\mu_{\mbox{m}})\qquad\mbox{for
}K=\mbox{STS, FRMS}\;,$ (18) $t_{\mbox{R}}^{K}(\sigma,\mu_{\mbox{m}})\approx
t_{\mbox{R}}^{*K}(\sigma)t_{0}^{K}(\mu_{\mbox{m}})\qquad\mbox{for
}K=\mbox{STS, FRMS}\;,$ (19)

where the superscript $K$ stands for the kind of mechanical load. Given this
factorisation, we have

$\tilde{t}^{K}(\sigma)=\frac{t_{\mbox{L}}^{*K}(\sigma)}{t_{\mbox{R}}^{*K}(\sigma)}\qquad\mbox{for
}K=\mbox{STS, FRMS}\;.$ (20)

Figure 5: Characteristic times $t_{\mbox{L}}$ (panels (a-c)) and
$t_{\mbox{R}}$ (panels (d-f)) as a function of the membrane viscosity
$\mu_{\mbox{m}}$ (see Sec. 4) for the three performed simulations: stretching
simulation (STS,  , Sec. 3.1), shear simulation (SHS,  , Sec. 3.2), four-roll
mill simulation (FRMS,  , Sec. 3.3), for different values of stress $\sigma$
(from lightest to darkest color): $\sigma=0.001$ Pa ( ,  ,  ), $\sigma=0.01$
Pa ( ,  ,  ), $\sigma=0.1$ Pa ( ,  ,  ). Figure 6: Characteristic times
$t_{\mbox{L}}$ (first column of panels) and $t_{\mbox{R}}$ (second column of
panels) as well as the ratio $\tilde{t}=t_{\mbox{L}}/t_{\mbox{R}}$ (third
column of panels) are reported for the three simulations performed, i.e.,
stretching simulation (STS,  , panels (a-c), Sec. 3.1), shear simulation (SHS,
, panels (d-f), Sec. 3.2), four-roll mill simulation (FRMS,  , panels (g-i),
Sec. 3.3), for different values of surface elastic shear modulus
$k_{\mbox{S}}$ (from lightest to darkest color): $k_{\mbox{S}}=5.3\times
10^{-6}\mbox{ N m}^{-1}$ ( ,  ,  ), $k_{\mbox{S}}=53\times 10^{-6}\mbox{ N
m}^{-1}$ ( ,  ,  ), $k_{\mbox{S}}=530\times 10^{-6}\mbox{ N m}^{-1}$ ( ,  ,
). The value of membrane viscosity $\mu_{\mbox{m}}$ was kept fixed in all the
simulations ($\mu_{\mbox{m}}=1.59\times 10^{-7}\mbox{ m Pa s}$). The red
dashed line represents the reference value for the symmetric case, i.e.,
$\tilde{t}=1$.

To make progress, we investigated what is the physical ingredient at the core
of the factorisation given in Eqs. (18)-(19), or, alternatively, why the SHS
does not factorize as in Eqs. (18)-(19). We investigated if such non-
factorisation could be related to the oscillations of $D(t)$ that appear only
in the SHS (see Sec. 3): however, since $\tilde{t}^{\mbox{{SHS}}}$ does not
factorise even when the deformation $D(t)$ does not oscillate (e.g., for small
values of the membrane viscosity $\mu_{\mbox{m}}$ and/or the shear rate
$\dot{\gamma}$), we conclude that these oscillations cannot be fully
responsible for the non-factorisation. We rather think that the difference
between SHS and STS/FRMS is mainly due to the different dynamics that are
induced by the mechanical load (see Tab. 1). In fact, regarding the SHS, one
can split the velocity gradient
$\frac{\partial\mathbf{u}}{\partial\mathbf{x}}$ in the symmetric (rotational)
and antisymmetric (elongational) parts:

$\frac{\partial\mathbf{u}}{\partial\mathbf{x}}=\left(\begin{matrix}0&\dot{\gamma}\\\
0&0\end{matrix}\right)=\left(\begin{matrix}0&\frac{\dot{\gamma}}{2}\\\
\frac{\dot{\gamma}}{2}&0\end{matrix}\right)+\left(\begin{matrix}0&\frac{\dot{\gamma}}{2}\\\
-\frac{\dot{\gamma}}{2}&0\end{matrix}\right)\;,$ (21)

where the only two components in the shear plane are reported. The rotational
part causes the rolling motion of the membrane (see Fig. 4, panels (a) and
(b)), while the elongational one tends to deform the RBC and pushes the main
diameter to an angle of $\pi/4$ with respect to the shear direction; an
increase in the membrane viscosity causes an increase in the time needed for
the membrane to adapt to the flow and to deform; meanwhile, the rotational
component promotes a rotation of the main diameter. Overall, the increase in
membrane viscosity $\mu_{\mbox{m}}$ leads to a decrease of the average
deformation $D_{\mbox{av}}$ (see also Fig. 4, panel (c)). To make these
arguments clearer, we have made two videos available in the ESI†: in one we
show the simulation with $\dot{\gamma}=123$ s-1 and $\mu_{\mbox{m}}=3.18\times
10^{-7}\mbox{ m Pa s}$, while in the other one the simulation with
$\dot{\gamma}=123$ s-1 and $\mu_{\mbox{m}}=0\mbox{ m Pa s}$ is reported. In
both cases, the tank-treading motion of the membrane appears, but, for
$\mu_{\mbox{m}}=0\mbox{ m Pa s}$, the membrane deforms more than in the case
with $\mu_{\mbox{m}}=3.18\times 10^{-7}\mbox{ m Pa s}$. In Fig. 4, panel (c),
we report the average deformation $D_{\mbox{av}}$ as a function of the stress
$\sigma$ for all three mechanical loads at changing the membrane viscosity
$\mu_{\mbox{m}}$. As already observed in [33] for the STS, we found that
$D_{\mbox{av}}$ is not sensitive to the value of membrane viscosity
$\mu_{\mbox{m}}$. Moreover, in the FRMS we found that $D_{\mbox{av}}$ shows
very little dependency on $\mu_{\mbox{m}}$, at least in the range of
$\mu_{\mbox{m}}$ and $\sigma$ analysed. Hence, for STS and FRMS, we report
only points for $\mu_{\mbox{m}}=3.18\times 10^{-7}$ m Pa s. It emerges that,
in the SHS, the average deformation $D_{\mbox{av}}$ saturates at a constant
value: an increase in the shear rate $\dot{\gamma}$ causes an initial increase
of the average deformation $D_{\mbox{av}}$; then, $D_{\mbox{av}}$ reaches a
plateau and increasing the shear rate $\dot{\gamma}$ beyond a certain value
does not result in an increased deformation. The higher the membrane viscosity
$\mu_{\mbox{m}}$, the lower is the value of $\dot{\gamma}$ for which the
plateau is reached. Furthermore, when compared to the STS and FRMS, we can see
that on the same values of $\sigma$ the average deformation $D_{\mbox{av}}$ is
much smaller in the SHS. Again, this is due to the rotation of the membrane
during the loading. In both STS and FRMS, the membrane does not rotate, so
that the energy injected by the flow is used to deform the membrane. These
investigations reveal that it is impossible to predict the loading and
relaxation times if we only know the deformation and have no information about
the kind of mechanical load.
In view of the above considerations on the deformation, it appears also
natural to study the characteristic times as a function of the average
deformation $D_{\mbox{av}}$. We performed this analysis (see ESI†, Fig. 2),
confirming the picture displayed in Fig. 3: again, $\tilde{t}$ shows a
collapse for the STS and the FRMS and does not depend on the value of the
membrane viscosity $\mu_{\mbox{m}}$; for the SHS $\tilde{t}$ shows a
dependency on both $D_{\mbox{av}}$ and the membrane viscosity
$\mu_{\mbox{m}}$. The results on $t_{\mbox{L}}(D_{\mbox{av}})$ and
$t_{\mbox{R}}(D_{\mbox{av}})$ (Fig. 2 in ESI†, panels (a) and (b),
respectively) further confirm that in general there is no correlation between
the degree of deformation of the membrane and the characteristic times for
different kinds of mechanical loads.
In our previous work [33], we have already seen that, for small forces,
$t_{\mbox{R}}$ is linear in $\mu_{\mbox{m}}$, in agreement with literature
predictions [19]. Now we can go further, and we study the dependency of both
$t_{\mbox{L}}$ and $t_{\mbox{R}}$ as a function of $\mu_{\mbox{m}}$ for
different values of $\sigma$. This will help further to determine to what
degree these two kinds of dynamics can be regarded as different dynamics [43,
42]. In Fig. 5, we report both $t_{\mbox{L}}(\mu_{\mbox{m}})$ and
$t_{\mbox{R}}(\mu_{\mbox{m}})$ (first and second row of panels, respectively)
for three values of $\sigma$ spanning two orders of magnitude as well as their
linear fit (whose coefficients are reported in ESI†, Tab. 2) for all three
simulations (STS in panel (a) and (d); SHS in panel (b) and (e); FRMS in panel
(c) and (f)). In all three simulations, for a fixed value of $\sigma$, the
linear approximation is reasonably good. For small values of $\sigma$ (e.g.,
$\sigma=0.001$ Pa), both $t_{\mbox{L}}$ and $t_{\mbox{R}}$ are similar for all
three simulations; that is not surprising, since at small values of stress
$\sigma$ the intrinsic properties of the membrane arise. Regarding the
sensitivity of the linear trend with respect to a change in $\sigma$, we
observe different behaviours in the two dynamics. Regarding the loading
dynamics, we observe that for high values of $\sigma$, the three load
mechanisms provide similar linear fits, while in the intermediate region of
$\sigma$, the SHS shows a different behaviour than the STS and FRMS. Regarding
the relaxation dynamics, the variability of the linear trends with the value
of the stress is more pronounced in presence of hydrodynamical forces (i.e.,
SHS and FRMS), while in the STS the linear behaviour of $t_{\mbox{R}}$ with
respect to $\mu_{\mbox{m}}$ is only slightly perturbed by a change in the
stress if compared to the others. These quantitative observations are
summarized in Tab. 2 in ESI†.
A dimensionless analysis could be performed to try to gain a deeper insight
into the problem; however, if we try to make both the characteristic times
$t_{\mbox{L}}$ and $t_{\mbox{R}}$ and the shear rates $\dot{\gamma}$ and
$\dot{\gamma}_{\mbox{FRMS}}$ dimensionless by using the characteristic elastic
time $t_{\mbox{el}}=\mu_{\mbox{out}}r/k_{\mbox{S}}$, as well as the force $F$
by using the characteristic elastic force $F_{\mbox{el}}=rk_{\mbox{S}}$, we
get only a rescaling along the x and y axis. Making the membrane viscosity
$\mu_{\mbox{m}}$ dimensionless by introducing the Boussinesq number
$\text{Bq}=\mu_{\mbox{m}}/r\mu_{\mbox{out}}$ (see also [69, 38, 33]) results
again in a constant rescaling of all the values, without giving a new insight
into the problem. In principle, one can look for some more refined non-
dimensionalisation procedure combining membrane viscosity and
rotational/deformation contributions: in the case of the SHS, for example,
this would mean to find a shear time dependent on both Bq and Ca. This would
anyhow require a more precise knowledge (e.g. a phenomenological model [19])
on how the rotational/deformation contributions couple with the membrane
viscosity effects.
Next, we discuss the parameters $\delta_{\mbox{L}}$ and $\delta_{\mbox{R}}$
used to improve the fit (see Eqs. (12)-(13)-(14)). In all three setups,
$\delta_{\mbox{L}}$ and $\delta_{\mbox{R}}$ are close to one, especially in
the STS and FRMS (see Fig. 3 in ESI†). The biggest deviation can be found
during the loading in the SHS, where the parameter $\delta_{\mbox{L}}$ seems
to tend asymptotically to $\delta_{\mbox{L}}\approx 0.6$ at increasing values
of shear rate $\dot{\gamma}$: this deviation from 1 reflects the effect of the
kind of mechanical load also on $\delta_{\mbox{L}}$, showing that, during the
loading in the SHS, $D(t)$ is not that close to an exponential function, and
then multiple time scales arise [33]. Indeed, having fitting parameters
$\delta_{\mbox{L}}$ and $\delta_{\mbox{R}}$ different from 1 is symptomatic of
the presence of multiple loading and relaxation times, respectively. This
issue has already been investigated for both SHS and STS during the relaxation
dynamics in our earlier study [33]; here we go deeper with the investigation
during the loading dynamics. We have monitored the time evolution of the
deformation $\frac{D(t)}{D_{\mbox{av}}}$ in log-lin scale (see Figs. 4-5-6 in
ESI†). The initial stage of the loading/relaxation process is well
characterised by a single "dominant" time scale, and only later, when the
difference of $D(t)$ from $D_{\mbox{av}}$ (during the loading) or from 0
(during the relaxation) is less (or even much less) than about $10\%$, other
time scales appear. As already pointed out in [33], it is difficult to make
quantitative assessments on the "late" dynamics, because then the deformation
is close to its steady value (for the loading) and/or to the rest value (for
the relaxation), and in this situation, discretisation errors could have more
influence. One could make a deeper analysis of these multiple relaxation times
and fit data with two (or more) characteristic times, both during loading and
relaxation; however, we notice that the values of $t_{\mbox{L}}$ and
$t_{\mbox{R}}$ found by fitting with the stretched exponential (i.e., by using
$\delta_{\mbox{L}},\delta_{\mbox{R}}\neq 1$) are in good agreement with the
"dominant" time scale (see Figs. 4-5-6 in ESI†). On the other hand, the small
difference between the fitted characteristic times and the "dominant"
characteristic times also explains why the ratio $\tilde{t}$ becomes greater
than 1 for one case in the FRMS (see Fig. 2 panel (i) and Fig. 7 in ESI†).
Before closing this section, we also discuss the effect of the surface elastic
shear modulus $k_{\mbox{S}}$ (see Eq. (7)) on the loading and relaxation
times. As already stated in the introduction, some pathologies can affect the
value of membrane elasticity [11, 16, 52, 53, 54, 55, 56, 57, 58, 59]: for
example, for RBCs infected with the malaria parasite Plasmodium falciparum,
experiments with optical tweezers estimated values of elastic shear modulus
ranging from $k_{\mbox{S}}=5.3\times 10^{-6}\mbox{ N m}^{-1}$ (i.e., for the
healthy RBC) to $k_{\mbox{S}}=100\times 10^{-6}\mbox{ N m}^{-1}$ [16]. For
this purpose, we fixed the value of membrane viscosity
$\mu_{\mbox{m}}=1.59\times 10^{-7}\mbox{ m Pa s}$ and we varied the elastic
shear modulus $k_{\mbox{S}}$ in a range of two orders of magnitude: from
$k_{\mbox{S}}=5.3\times 10^{-6}\mbox{ N m}^{-1}$ to $k_{\mbox{S}}=530\times
10^{-6}\mbox{ N m}^{-1}$. Results are reported in Fig. 6. For all three kinds
of mechanical loads simulated, increasing the value of the elastic shear
modulus by a factor 10 or 100 results in a reduction of the characteristic
times by about the same factor, as expected [46, 19]. Moreover, we analysed
the ratio $\tilde{t}$, finding that it gets closer to 1 at increasing values
of $k_{\mbox{S}}$. Mechanical balance at the interface tells us that the
relative importance of viscous to elastic effects rescales as the ratio
between the shear rate $\dot{\gamma}$ developed in the fluid and the surface
elastic modulus $k_{\mbox{S}}$ of the membrane. The shear rate developed in
the fluid during both loading and relaxation is larger at increasing load
strength; thus, for a fixed $\dot{\gamma}$, by increasing $k_{\mbox{S}}$ we
should have the same results observed for $k_{\mbox{S}}=5.3\times
10^{-6}\mbox{ N m}^{-1}$ at shear rates that are smaller by a factor given by
the ratio of the $k_{\mbox{S}}$’s. In other words, we expect
$\tilde{t}\rightarrow 1$ when $k_{\mbox{S}}$ gets very large. This fact is
confirmed by our results. We hasten to remark that these are only preliminary
results, for at least two reasons: first, a proper dimensionless analysis of
the governing equations [38, 57] reveals that also the importance of the
bending modulus $k_{\mbox{B}}$ with respect to $k_{\mbox{S}}$ needs to be
taken into account via a suitable dimensionless number $k_{\mbox{B}}^{*}$.
Since it is not known whether certain blood-related pathologies, such as
Plasmodium falciparum malaria parasite infection, affect the value of the
bending modulus $k_{\mbox{B}}$ [56, 57, 58, 59], we have kept it fixed at the
value for the healthy RBC (see Sec. 2). Therefore, in our simulations, at
changing $k_{\mbox{S}}$, the dimensionless number
$k_{\mbox{B}}^{*}=\frac{k_{\mbox{B}}}{r^{2}\ k_{\mbox{S}}}$ is changing: a
more comprehensive analysis on the effects of $k_{\mbox{S}}$ should be done by
assessing also the impact of a variation in $k_{\mbox{B}}^{*}$ separately.
Second, the SLS model we implemented to take into account the viscoelastic
effects (see Sec. 2) contains an artificial elastic contribution $k^{\prime}$
that is proportional to $k_{\mbox{S}}$ and needs to be tuned in order to
recover physical results (see [33, 69] for further details and validations).
Varying the value of the elastic shear modulus $k_{\mbox{S}}$ modifies also
the value of $k^{\prime}$: if and how this affects the physical dynamics of
the RBC at the very large $k_{\mbox{S}}$ we considered requires a detailed
computational analysis by its own, which is out of the scope of the present
paper. All these considerations surely warrant dedicated studies in the
future.

##  5 Conclusions

A comprehensive characterisation of the viscoelastic properties of the RBC
membrane, as well as the way the membrane responds to an external force, is of
paramount interest in different fields, from the detection of pathologies [11,
16, 19], to the design of biomedical devices [24, 25, 26, 27]. A paradigmatic
example is provided by ventricular assist devices [75] where RBCs evolve in a
complex flow and their fate is closely linked to their residence time inside
the device: if the residence time inside the impeller (that is the region
where the RBCs experience a wide range of stress) is much shorter then the
loading time, RBCs deform without reaching a steady state configuration; on
the contrary, a higher residence time leads to a deformation that can cause
hemolysis (that is, the release of the cytoplasm into surrounding plasma due
to damage or rupture of the membrane).
In general, when an external force acts on a viscoelastic membrane, two main
kinds of dynamics arise: the loading and the relaxation dynamics with
associated times $t_{\mbox{L}}$ and $t_{\mbox{R}}$. Earlier investigations
pointed to the fact that these two kinds of dynamics are two distinct
processes, since during the relaxation there is no external force to drive the
membrane, in contrast to the loading [43, 42]. To the best of our knowledge,
however, an exhaustive comparative characterisation of these two kinds of
dynamics has never been conducted. This motivated our work to investigate
these two kinds of dynamics with different setups that involve different
typologies of mechanical loads (whose main features are summarised in Tab. 1)
while performing a parametric study on the values of membrane viscosity
$\mu_{\mbox{m}}$. The latter choice is motivated by the large variability of
membrane viscosity values reported in the literature [46, 18, 45, 47, 22, 48,
49, 17, 19, 32].
The two kinds of dynamics are symmetrical
($\tilde{t}=t_{\mbox{L}}/t_{\mbox{R}}\rightarrow 1$) in the limit of small
load strengths ($\sigma\rightarrow 0$), i.e., in the limit where the response
function of the RBC is dominated by the "intrinsic" properties of the
membrane; in marked contrast, we found an asymmetry in the two kinds of
dynamics for load strengths large enough
($\tilde{t}=t_{\mbox{L}}/t_{\mbox{R}}\neq 1$ for $\sigma>0$), meaning that the
loading dynamics is always faster than the relaxation one. We found that the
asymmetry profoundly depends on the kind of mechanical load and we have
demonstrated this non-universality via a quantitative study in terms of the
applied load strength $\sigma$ and the value of membrane viscosity
$\mu_{\mbox{m}}$. There are some realistic load mechanisms, like shear flows,
that make the membrane rotate during loading while leaving the membrane
relaxing to the shape at rest without rotation: in this case, the contribution
that the membrane viscosity gives to the characteristic times $t_{\mbox{L}}$
and $t_{\mbox{R}}$ differs, and then the ratio $\tilde{t}$ is a function of
both the stress $\sigma$ and the membrane viscosity $\mu_{\mbox{m}}$. From the
other side, there are other realistic load mechanisms, like the stretching
with optical tweezers or the deformation with an elongational flow, in which
the membrane deforms without rotating during both processes. In this case, the
contribution given by the membrane viscosity $\mu_{\mbox{m}}$ to the
characteristic times is the same during both loading and relaxation, and as a
consequence, the ratio $\tilde{t}$ is a function of the stress $\sigma$ only.
Even though we showed that both loading and relaxation dynamics are not
universal, we found that for a given value of the stress $\sigma$, a linear
increase of the characteristic times as a function of the membrane viscosity
$\mu_{\mbox{m}}$ is a fair approximation in all cases.
Finally, since some blood-related diseases [11, 16, 52, 53, 54, 55, 56, 57,
58, 59] can alter the values of the elastic shear modulus $k_{\mbox{S}}$, we
also investigated the loading and relaxation dynamics at changing
$k_{\mbox{S}}$ for a fixed value of $\mu_{\mbox{m}}$, finding that larger
values of $k_{\mbox{S}}$ promote symmetrization ($\tilde{t}\rightarrow 1$) of
the dynamics.
We argue our findings offer interesting physical and practical insights on the
response function and the unsteady dynamics of RBCs driven by realistic
mechanical loads.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

The authors acknowledge L. Biferale and G. Koutsou. This project has received
funding from the European Union Horizon 2020 Research and Innovation Program
under the Marie Skłodowska-Curie grant agreement No. 765048. We also
acknowledge support from the project “Detailed Simulation of Red blood Cell
Dynamics accounting for membRane viscoElastic propertieS” (SorCeReS, CUP No.
E84I19002470005) financed by the University of Rome “Tor Vergata” (“Beyond
Borders 2019” call).

## References

  * Popel and Johnson [2005] A. S. Popel and P. C. Johnson, _Annu. Rev. Fluid Mech._ , 2005, 37, 43–69.
  * Skalak and Branemark [1969] R. Skalak and P. Branemark, _Science_ , 1969, 164, 717–719.
  * Krüger [2012] T. Krüger, _Computer simulation study of collective phenomena in dense suspensions of red blood cells under shear_ , Springer Science & Business Media, 2012.
  * Mountrakis [2015] L. Mountrakis, _Transport of blood cells studied with fully resolved models_ , 9789462597754, 2015.
  * Janoschek [2013] F. Janoschek, _Mesoscopic simulation of blood and general suspensions in flow_ , Technische Universiteit Eindhoven, 2013.
  * Gross _et al._ [2014] M. Gross, T. Krüger and F. Varnik, _Soft Matter_ , 2014, 10, 4360–4372.
  * Krüger _et al._ [2013] T. Krüger, M. Gross, D. Raabe and F. Varnik, _Soft Matter_ , 2013, 9, 9008–9015.
  * Fedosov _et al._ [2010] D. A. Fedosov, B. Caswell and G. E. Karniadakis, _Biophysical journal_ , 2010, 98, 2215–2225.
  * Pan _et al._ [2011] W. Pan, D. A. Fedosov, B. Caswell and G. E. Karniadakis, _Microvascular research_ , 2011, 82, 163–170.
  * Fedosov _et al._ [2010] D. A. Fedosov, B. Caswell and G. E. Karniadakis, _Computer Methods in Applied Mechanics and Engineering_ , 2010, 199, 1937–1948.
  * Krüger _et al._ [2014] T. Krüger, D. Holmes and P. V. Coveney, _Biomicrofluidics_ , 2014, 8, 054114.
  * Fedosov _et al._ [2011] D. A. Fedosov, W. Pan, B. Caswell, G. Gompper and G. E. Karniadakis, _Proceedings of the National Academy of Sciences_ , 2011, 108, 11772–11777.
  * Guckenberger _et al._ [2018] A. Guckenberger, A. Kihm, T. John, C. Wagner and S. Gekle, _Soft Matter_ , 2018, 14, 2032–2043.
  * Fedosov _et al._ [2014] D. A. Fedosov, H. Noguchi and G. Gompper, _Biomechanics and modeling in mechanobiology_ , 2014, 13, 239–258.
  * Mills _et al._ [2004] J. Mills, L. Qie, M. Dao, C. Lim and S. Suresh, _MCB-TECH SCIENCE PRESS_ , 2004, 1, 169–180.
  * Suresh _et al._ [2005] S. Suresh, J. Spatz, J. P. Mills, A. Micoulet, M. Dao, C. Lim, M. Beil and T. Seufferlein, _Acta biomaterialia_ , 2005, 1, 15–30.
  * Braunmüller _et al._ [2012] S. Braunmüller, L. Schmid, E. Sackmann and T. Franke, _Soft Matter_ , 2012, 8, 11240–11248.
  * Chien _et al._ [1978] S. Chien, K. L. Sung, R. Skalak, S. Usami and A. Tözeren, _Biophysical Journal_ , 1978, 24, 463–487.
  * Prado _et al._ [2015] G. Prado, A. Farutin, C. Misbah and L. Bureau, _Biophysical Journal_ , 2015, 108, 2126–2136.
  * Hénon _et al._ [1999] S. Hénon, G. Lenormand, A. Richert and F. Gallet, _Biophysical Journal_ , 1999, 76, 1145 – 1151.
  * Hochmuth _et al._ [1979] R. Hochmuth, P. Worthy and E. Evans, _Biophysical journal_ , 1979, 26, 101–14.
  * Baskurt and Meiselman [1996] O. Baskurt and H. Meiselman, _Biorheology_ , 1996, 33, 489 – 503\.
  * Bronkhorst _et al._ [1995] P. Bronkhorst, G. Streekstra, J. Grimbergen, E. Nijhof, J. Sixma and G. Brakenhoff, _Biophysical journal_ , 1995, 69, 1666–1673.
  * Murakami _et al._ [1979] T. Murakami, L. R. Golding, G. B. Jacobs, S. Takatani, R. Sukalac, H. Harasaki and Y. NONE, _Jinko Zoki_ , 1979, 8, 636–639.
  * Nonaka _et al._ [2001] K. Nonaka, J. Linneweber, S. Ichikawa, M. Yoshikawa, S. Kawahito, M. Mikami, T. Motomura, H. Ishitoya, I. Nishimura, D. Oestmann _et al._ , _Artificial Organs_ , 2001, 25, 675–682.
  * Behbahani _et al._ [2009] M. Behbahani, M. Behr, M. Hormes, U. Steinseifer, D. Arora, O. Coronado and M. Pasquali, _European Journal of Applied Mathematics_ , 2009, 20, 363–397.
  * Arora _et al._ [2006] D. Arora, M. Behr and M. Pasquali, _Artificial organs_ , 2006, 30, 539–547.
  * Zhu and Asaro [2019] Q. Zhu and R. J. Asaro, _Physics of Fluids_ , 2019, 31, 051901\.
  * Cordasco and Bagchi [2017] D. Cordasco and P. Bagchi, _Physics of Fluids_ , 2017, 29, 041901\.
  * Cordasco and Bagchi [2016] D. Cordasco and P. Bagchi, _Journal of Fluid Mechanics_ , 2016, 800, 484.
  * Cordasco and Bagchi [2014] D. Cordasco and P. Bagchi, _Journal of Fluid Mechanics_ , 2014, 759, 472.
  * Fedosov [2010] D. A. Fedosov, _Multiscale modeling of blood flow and soft matter_ , Brown University, 2010.
  * Guglietta _et al._ [2020] F. Guglietta, M. Behr, L. Biferale, G. Falcucci and M. Sbragaglia, _Soft Matter_ , 2020, 16, 6191–6205.
  * Guglietta _et al._ [2020] F. Guglietta, M. Behr, L. Biferale, G. Falcucci and M. Sbragaglia, _arXiv preprint arXiv:2012.06190_ , 2020.
  * Dao _et al._ [2003] M. Dao, C. T. Lim and S. Suresh, _Journal of the Mechanics and Physics of Solids_ , 2003, 51, 2259–2280.
  * Li and Zhang [2020] P. Li and J. Zhang, _Biomechanics and Modeling in Mechanobiology_ , 2020, 1–15.
  * Li and Zhang [2021] P. Li and J. Zhang, _Cardiovascular Engineering and Technology_ , 2021\.
  * Yazdani and Bagchi [2013] A. Yazdani and P. Bagchi, _Journal of Fluid Mechanics_ , 2013, 718, 569–595.
  * Graessel _et al._ [2021] K. Graessel, C. Bächer and S. Gekle, _Journal of Fluid Mechanics_ , 2021, 910, A46.
  * Noguchi and Gompper [2005] H. Noguchi and G. Gompper, _Physical Review E_ , 2005, 72, 011901\.
  * Noguchi and Gompper [2007] H. Noguchi and G. Gompper, _Physical review letters_ , 2007, 98, 128103\.
  * Barthes-Biesel [2016] D. Barthes-Biesel, _Annual Review of fluid mechanics_ , 2016, 48, 25–52.
  * Diaz _et al._ [2000] A. Diaz, N. Pelekasis and D. Barthes-Biesel, _Physics of Fluids_ , 2000, 12, 948–957.
  * Evans [1989] E. A. Evans, _Methods in enzymology_ , Elsevier, 1989, vol. 173, pp. 3–35.
  * Hochmuth _et al._ [1979] R. M. Hochmuth, P. Worthy and E. A. Evans, _Biophysical journal_ , 1979, 26, 101–114.
  * Evans and Hochmuth [1976] E. Evans and R. Hochmuth, _Biophysical journal_ , 1976, 16, 1–11.
  * Tran-Son-Tay _et al._ [1984] R. Tran-Son-Tay, S. Sutera and P. Rao, _Biophysical journal_ , 1984, 46, 65–72.
  * Riquelme _et al._ [2000] B. D. Riquelme, J. R. Valverde and R. J. Rasia, Optical Diagnostics of Biological Fluids V, 2000, pp. 132–140.
  * Tomaiuolo and Guido [2011] G. Tomaiuolo and S. Guido, _Microvascular research_ , 2011, 82, 35–41.
  * Keller and Skalak [1982] S. R. Keller and R. Skalak, _Journal of Fluid Mechanics_ , 1982, 120, 27–47.
  * Matteoli _et al._ [2021] P. Matteoli, F. Nicoud and S. Mendez, _Physical Review Fluids_ , 2021, 6, 043602.
  * Suresh [2006] S. Suresh, _Journal of materials research_ , 2006, 21, 1871–1877.
  * Brandao _et al._ [2003] M. Brandao, A. Fontes, M. Barjas-Castro, L. Barbosa, F. Costa, C. Cesar and S. Saad, _European journal of haematology_ , 2003, 70, 207–211.
  * Briole _et al._ [2021] A. Briole, T. Podgorski and B. Abou, _Soft Matter_ , 2021.
  * Brandao _et al._ [2003] M. Brandao, S. Saad, C. Cezar, A. Fontes, F. Costa and M. Barjas-Castro, _Vox sanguinis_ , 2003, 85, 213–215.
  * Fedosov _et al._ [2011] D. Fedosov, B. Caswell, S. Suresh and G. Karniadakis, _Proceedings of the National Academy of Sciences_ , 2011, 108, 35–39.
  * Luo _et al._ [2013] Z. Y. Luo, S. Q. Wang, L. He, F. Xu and B. F. Bai, _Soft Matter_ , 2013, 9, 9651–9660.
  * Ye _et al._ [2013] T. Ye, N. Phan-Thien, B. C. Khoo and C. T. Lim, _Biophysical journal_ , 2013, 105, 1103–1109.
  * Hosseini and Feng [2012] S. M. Hosseini and J. J. Feng, _Biophysical journal_ , 2012, 103, 1–10.
  * Malaspinas _et al._ [2010] O. Malaspinas, N. Fiétier and M. Deville, _Journal of Non-Newtonian Fluid Mechanics_ , 2010, 165, 1637–1653.
  * Krüger _et al._ [2016] T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva and E. M. Viggen, _The Lattice Boltzmann Method - Principles and Practice_ , 2016.
  * Succi [2001] S. Succi, _The lattice Boltzmann equation: for fluid dynamics and beyond_ , Oxford university press, 2001.
  * Guo _et al._ [2002] Z. Guo, C. Zheng and B. Shi, _Phys. Rev. E_ , 2002, 65, 046308\.
  * Frijters _et al._ [2015] S. Frijters, T. Krüger and J. Harting, _Computer Physics Communications_ , 2015, 189, 92 – 98.
  * Evans and Fung [1972] E. Evans and Y.-C. Fung, _Microvascular research_ , 1972, 4, 335–347.
  * Skalak _et al._ [1973] R. Skalak, A. Tozeren, R. P. Zarda and S. Chien, _Biophysical journal_ , 1973, 13, 245–264.
  * Helfrich [1973] W. Helfrich, _Zeitschrift fur Naturforschung. Teil C: Biochemie, Biophysik, Biologie, Virologie_ , 1973, 28, 693–703.
  * Gompper and Schick [2008] G. Gompper and M. Schick, _Soft Matter: Lipid Bilayers and Red Blood Cells_ , Wiley-VCH, 2008.
  * Li and Zhang [2019] P. Li and J. Zhang, _International Journal for Numerical Methods in Biomedical Engineering_ , 2019, 35, e3200.
  * Barthes-Biesel and Sgaier [1985] D. Barthes-Biesel and H. Sgaier, _Journal of Fluid Mechanics_ , 1985, 160, 119–135.
  * Rallison [1980] J. Rallison, _Journal of Fluid Mechanics_ , 1980, 98, 625–633.
  * Minetti _et al._ [2019] C. Minetti, V. Audemar, T. Podgorski and G. Coupier, _Journal of fluid mechanics_ , 2019, 864, 408–448.
  * Mauer _et al._ [2018] J. Mauer, S. Mendez, L. Lanotte, F. Nicoud, M. Abkarian, G. Gompper and D. A. Fedosov, _Physical review letters_ , 2018, 121, 118103.
  * Alberts _et al._ [2002] B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts and P. Walter, _Molecular biology of the cell_ , New York: Garland, 4th edn, 2002.
  * Haßler _et al._ [2020] S. Haßler, A. M. Ranno and M. Behr, _Computer Methods in Applied Mechanics and Engineering_ , 2020, 369, 113171.