Dynamics of a Bilayer Membrane Coupled to a Two-dimensional Cytoskeleton: Scale Transfers of Membrane Deformations

TL;DR

This paper models the complex dynamics of a lipid bilayer membrane coupled with a cytoskeleton, revealing how microscopic interactions influence large-scale membrane deformations and relaxation behaviors.

Contribution

It introduces a theoretical framework accounting for lattice structure, intermonolayer friction, and prestress, highlighting their effects on membrane dynamics and mode coupling.

Findings

Long-wavelength deformations excite short-wavelength modes.

Slow relaxation is dominated by intermonolayer friction.

Long-term forces can induce large-scale membrane responses.

Abstract

We theoretically investigate the dynamics of a floating lipid bilayer membrane coupled with a two-dimensional cytoskeleton network, taking into explicitly account the intermonolayer friction, the discrete lattice structure of the cytoskeleton, and its prestress. The lattice structure breaks lateral continuous translational symmetry and couples Fourier modes with different wavevectors. It is shown that within a short time interval a long-wavelength deformation excites a collection of modes with wavelengths shorter than the lattice spacing. These modes relax slowly with a common renormalized rate originating from the long-wavelength mode. As a result, and because of the prestress, the slowest relaxation is governed by the intermonolayer friction. Reversely, and most interestingly, forces applied at the scale of the cytoskeleton for a sufficiently long time can cooperatively excite large-scale modes.

Figures

Tables

Table 1. Table 1: List of the parameter values used in numerical calculations.

$\sigma$	$\kappa$	$k$	$d$	$\nu$	$a$	$\eta$	$b$	$\mu +\zeta$
$N/m$	$J$	$N/m$	$m$	$N/m$	$m$	$Js/m^3$	$Js/m^4$	$Js/m^2$
${10}^{-11}$	$2\times {10}^{-20}$	$7\times {10}^{-2}$	${10}^{-9}$	${10}^{-6}$	${10}^{-7}$	${10}^{-3}$	$2\times {10}^8$	$2\times {10}^{-9}$

Table 2. Table 2: Approximate expressions of 0th order eigenvalues and eigenvectors for sufficiently small bare tension, σ ≪ σ c much-less-than 𝜎 subscript 𝜎 𝑐 \sigma\ll\sigma_{c} .

	${\gamma }_{+}^{\left(0\right)}\left(q\right)$	${\gamma }_{-}^{\left(0\right)}\left(q\right)$	$e_{+}\left(q\right)$	$e_{-}\left(q\right)$
(i) $b\sigma /\left(2k\eta \right)\ll q\ll q_c$	$k{q}^2/\left(2b\right)$	$\left(\sigma q+\kappa q^3\right)/\left(4\eta \right)$	$-\eta q/b$	${\left(qd\right)}^2$
(ii) $q\gg q_c$	$\tilde{\kappa }{q}^3/\left(4\eta \right)$	$k\kappa q^2/\left(2b\tilde{\kappa }\right)$	$-2k\eta d^{2}q/\left(b\tilde{\kappa }\right)$	$\tilde{\kappa }{q}^2/\left(2k\right)$

Equations

$$\displaystyle F_{0}=\int\mathrm{d}^{2}x\,\Big{[}$$

$$\displaystyle+\frac{k}{2}\sum_{\epsilon=\pm}\,[\rho^{\epsilon}+\epsilon d(\nabla^{2}h)]^{2}\Big{]},$$

$$F_{\mathrm{c}}=\frac{\nu }{4}\sum _{\ell }\sum _{\mathbf{n}}[h({\mathbf{R}}_{\ell })-h({\mathbf{R}}_{\ell }+\mathbf{n}){]}^2,$$

$$Q_{\mathbf{q}}=\{\mathbf{q}+2\pi m_{\alpha }{\mathbf{e}}^{\alpha }|m_{\alpha }\in \mathbb{Z}\}.$$

$$\eta {\hat{\nabla }}^2\mathbf{V}-\hat{\nabla }p=0,\hat{\nabla }\cdot \mathbf{V}=0,$$

$${\tau }_{ij}^{\pm }=\mu ({\partial }_i{v}_j^{\pm }+{\partial }_j{v}_i^{\pm })+(\zeta -\mu ){\delta }_{ij}\nabla \cdot {\mathbf{v}}^{\pm },$$

$$-{\partial }_i\left(\frac{\delta F}{\delta {\rho }^{\pm }}\right)+{\partial }_j{\tau }_{ij}^{\pm }\pm T_{iz}^{\pm }\mp b(v_i^{+}-v_i^{-})=0,$$

$$T_{zz}^{+}-T_{zz}^{-}=\frac{\delta F}{\delta h}.$$

$$\frac{\partial {\rho }^{\pm }}{\partial t}\simeq -\nabla \cdot {\mathbf{v}}^{\pm }.$$

$$4\eta d^{2}q\frac{\partial \hat{h}(\mathbf{q},t)}{\partial t}=-{\mathcal{F}}_{\mathbf{q}}\left[\frac{\delta F}{\delta \hat{h}}\right]+u_h(\mathbf{q},t),$$

$$\frac{2c(q)}{q^2}\frac{\partial \rho (\mathbf{q},t)}{\partial t}=-{\mathcal{F}}_{\mathbf{q}}\left[\frac{\delta F}{\delta \rho }\right]+u_{\rho }(\mathbf{q},t),$$

$$\displaystyle\gamma_{+}\simeq\frac{\sigma_{\rm eff}q+\kappa_{\rm eff}q^{3}}{4\eta},\quad\gamma_{-}\simeq\Big{(}\frac{k}{2b}+\frac{3^{1/2}d^{2}b}{4\eta^{2}}\nu\Big{)}q^{2},$$

$${\hat{h}}_{\mathrm{qe}}({\mathbf{q}}^{\prime };\hat{h}(\mathbf{q}))=\frac{{\rho }_{\mathrm{qe}}({\mathbf{q}}^{\prime };\hat{h}(\mathbf{q}))}{d^2{q}^{\prime 2}}\simeq -\frac{2\nu \hat{h}(\mathbf{q})K_{\mathbf{q}}}{\sqrt{3}\kappa a^2{q}^{\prime 4}},$$

$${\hat{h}}_{\mathrm{st}}(\mathbf{q})\simeq -\frac{\Delta {\sigma }_{\mathrm{eff}}+\Delta {\kappa }_{\mathrm{eff}}{q}^2}{{\sigma }_{\mathrm{eff}}+{\kappa }_{\mathrm{eff}}{q}^2}\frac{\bar{u}}{\kappa d^2{q}_1^4},$$

$$F_{\mathrm{c}}=\frac{1}{2}\int {\mathrm{d}}^{2}x\frac{\nu }{2}\sum _{\ell }\sum _{\mathbf{n}}\delta (\mathbf{x}-{\mathbf{R}}_{\ell })[h(\mathbf{x})-h(\mathbf{x}+\mathbf{n}){]}^2.$$

$$\displaystyle{\cal F}_{\bm{q}}\Big{[}\frac{\delta F}{\delta h}\Big{]}=$$

$$+\nu \sum _{\ell }\sum _{\mathbf{n}}{e}^{-i\mathbf{q}\cdot {\mathbf{R}}_{\ell }}[h({\mathbf{R}}_{\ell })-h({\mathbf{R}}_{\ell }+\mathbf{n})].$$

$$\sum _{\ell }{e}^{-i\mathbf{q}\cdot {\mathbf{R}}_{\ell }}=\frac{(2\pi {)}^2}{\sqrt{g}}\sum _{{\mathbf{q}}^{\prime }\in Q_{\mathbf{q}}}\delta ({\mathbf{q}}^{\prime }),$$

$$\displaystyle{\cal F}_{\bm{q}}\Big{[}\frac{\delta F}{\delta h}\Big{]}=$$

$$+\frac{\nu }{\sqrt{g}}\sum _{{\mathbf{q}}^{\prime }\in Q_{\mathbf{q}}}h({\mathbf{q}}^{\prime })K_{{\mathbf{q}}^{\prime }},$$

$$\eta ({\partial }_z^2-q^2)V_z-{\partial }_{z}p=0,$$

$$\eta ({\partial }_z^2-q^2)V_{\parallel }-iqp=0,$$

$$\eta ({\partial }_z^2-q^2)V_{\perp }=0,$$

$$iq{V}_{\parallel }+{\partial }_z{V}_z=0.$$

$$p=p^{\pm }{e}^{\mp qz},$$

$$V_{\perp }=v_{\perp }^{\pm }{e}^{\mp qz},$$

$$\displaystyle V_{z}=\Big{[}\dot{h}+\frac{p^{\pm}}{2\eta}z\Big{]}e^{\mp qz}=[\dot{h}\pm(\dot{h}\mp v^{\pm}_{\parallel})qz]e^{\mp qz},$$

$$\displaystyle V_{\parallel}=-\frac{1}{iq}\Big{[}\frac{p^{\pm}}{2\eta}\mp q\Big{(}\dot{h}+\frac{p^{\pm}}{2\eta}z\Big{)}\Big{]}e^{\mp qz}$$

$$=[v_{\parallel }^{\pm }-(i\dot{h}\pm v_{\parallel }^{\pm })qz]e^{\mp qz},$$

$$\displaystyle 4\eta q\frac{\partial h({\bm{q}})}{\partial t}=-{\cal F}_{\bm{q}}\Big{[}\frac{\delta F}{\delta h}\Big{]},$$

$$0=$$

$$\pm b[v_{\parallel }^{+}(\mathbf{q})-v_{\parallel }^{-}(\mathbf{q})].$$

$$\displaystyle 0=-\frac{2}{iq}c(q)\frac{\partial\rho({\bm{q}})}{\partial t}+iq{\cal F}_{\bm{q}}\Big{[}\frac{\delta F}{\delta\rho^{+}}-\frac{\delta F}{\delta\rho^{-}}\Big{]},$$

$$|\Psi (t){⟩}_{\mathbf{q}}\equiv \sum _{{\mathbf{q}}^{\prime }\in Q_{\mathbf{q}}}[\hat{h}({\mathbf{q}}^{\prime },t)|h;{\mathbf{q}}^{\prime }⟩+\rho ({\mathbf{q}}^{\prime },t)|\rho ;{\mathbf{q}}^{\prime }⟩],$$

$$\frac{\partial }{\partial t}|\Psi (t){⟩}_{\mathbf{q}}=-({\Gamma }_{\mathrm{SL}}+\nu {\Gamma }_h)|\Psi (t){⟩}_{\mathbf{q}},$$

$${\Gamma }_{\mathrm{SL}}(\mathbf{q})=$$

Full text

I Introduction

In biological materials, much attention has been paid to the dynamics from the viewpoint of nonequilibrium physics, because of the high complexity of composition, hydrodynamic interactions, active components, etc. Mizuno ; Weitz ; Mikhailov . In particular, shape relaxation and fluctuations of lipid bilayer membranes have intensively been studied both in, or near, equilibrium Brochard ; Sackmann ; Seifert and far from equilibrium Prost ; Garcia ; Turlier ; BetzPNAS ; ParkPNAS ; YasudaPRE . The dynamics of a bilayer membrane is determined by many factors, such as the viscosity of the surrounding fluid Brochard , the membrane bending rigidity, the inter-monolayer friction caused by relative lateral motions of two monolayers Seifert ; JBPRL2 , and possibly by active inclusions Prost . Red blood cell (RBC) membranes have further complexity because of the cytoskeleton which is attached to the lipid bilayer. It has been argued that the cytoskeleton plays crucial roles both in the statics and the dynamics, e.g., drastic effective tension increase in equilibrium GovPRL ; JBPRL ; JBEPL ; Popescu , tension decrease in the presence of ATP BetzPNAS ; Turlier , and enhanced non-equilibrium fluctuations on the scale of the cytoskeleton mesh size ParkPNAS .

In RBCs, the cytoskeleton consists of spectrin filaments forming a pre-stressed Turlier two-dimensional (2D) triangular lattice with a protein at each vertex embedded in the bilayer membrane. The lattice spacing $a\approx 100\,\mathrm{nm}$ is quite large, and therefore what matters in the membrane collective dynamics is not only the modes whose wavelengths are much larger than $a$, but also those having wavelengths smaller than $a$. To understand simultaneously the dynamics on such a wide range of spatial scales, we need to take into explicitly account the discrete nature of the lattice structure. The latter breaks lateral continuous translational symmetry, giving rise to a coupling between modes on different length scales and thus to a rich dynamical behavior.

This paper is organized as follows. In Sec. II, our model free energy is constructed on the basis of the previous theories for a bilayer membrane (without a cytoskeleton) Seifert and for a pre-stressed 2D cytoskeleton coupled to a membrane JBEPL . In Sec. III, following Ref. Seifert , the hydrodynamic equations are introduced, where we take into account the hydrodynamic flows of the surrounding fluid and of the monolayers, inter-monolayer friction between the monolayers. Then we obtain the coupled equations for the membrane variables by integrating out the flow velocity fields. In Sec. IV, we discuss how the cytoskeleton alters the dynamics of the membrane, when initially a large-scale deformation is imposed, and when force(s) are applied to small-scale mode(s) for a long time ($\gtrsim 10$ ms). Section V is devoted for discussion and summary. We also present our detailed calculations in the Appendices.

II Free energy

We consider out-of-plane deformations of a RBC membrane patch described by its height $h(x,y)$ above a reference plane $z=0$. Our model free energy is given by $F=F_{0}+F_{\rm c}$ as follows. We take into account the areal compression which is necessarily coupled with $h$ due to the finite thickness $d\approx 1\,{\rm nm}$ of the monolayers Seifert . Then the bilayer membrane free energy $F_{0}$ is given by

[TABLE]

where $\kappa$ is the bare bending rigidity, $\sigma$ the bare tension, $k$ the areal compression modulus and $\nabla=(\partial_{x},\partial_{y})$. In the above, $\rho^{+}$ (resp. $\rho^{-}$) denotes the dimensionless projected excess lipid density in the upper (resp. lower) monolayer Seifert .

Note that in our paper, the surface tension is not considered as a constant. The quantity $\sigma$ in Eq. (1) is only the background tension for a flat membrane with homogeneous, reference lipid density ($\rho^{\pm}=0$). The actual tension fluctuates about the zeroth-order tension $\sigma$ according to a model that is closely related to the area-difference-elasticity model Miao . Indeed, the term proportional to $k$ in Eq. (1), that involves the variables $\rho^{+}$ and $\rho^{-}$ is the excess energy associated with a local compression or dilation of the lipids in each monolayer. The actual tension (without the cytoskeleton) is $\sigma+k(\rho^{+}+d\nabla^{2}h)$ and $\sigma+k(\rho^{-}-d\nabla^{2}h)$ in the upper and lower monolayers, respectively.

The other contribution, $F_{\rm c}$, arises from the membrane–cytoskeleton coupling. We assume, to simplify, that the cytoskeleton network is a uniform triangular lattice without defects. An anchoring protein at each lattice site is embedded in the membrane and interacts with its nearest neighbor through an effective spring of relaxed length $a_{0}$ and stiffness $k_{\rm s}$ (Fig. 1a). In the ground state ($\rho=h=0$), the network forms a regular triangular lattice, and the lattice points $\{{\bm{R}}_{\ell}\}$ are expressed in terms of the primitive lattice vectors ${\bm{e}}_{\alpha}$ $(\alpha=1,2)$ as ${\bm{R}}_{\ell}=R_{\ell}^{\alpha}{\bm{e}}_{\alpha}$ with $R_{\ell}^{\alpha}\in\mathbb{Z}$ positive or negative integers (see Fig. 1b). The lattice spacing is $a=|{\bm{e}}_{\alpha}|$.

If the out-of-plane deformation of the membrane is sufficiently small, $F_{\rm c}$ is given by JBEPL

[TABLE]

where $\nu=k_{\rm s}(1-a_{0}/a)$ is the effective stiffness of the harmonic potentials associated with the out-of-plane deformations, and $\sum_{\bm{n}}$ denotes the sum over the nearest neighbor sites. Note that $\nu$ is nonzero only if the lattice is pre-stressed ($a\neq a_{0}$). In Ref. Turlier , as a result of fitting their experimental data, it has been shown that the cytoskeleton in healthy RBCs is naturally stretched (by about $4\%$) while the bare membrane tension can be negative. Let us introduce the in-plane Fourier transform as ${\cal F}_{\bm{q}}[\cdots]\equiv\int{\rm d}^{2}x\,(\cdots)e^{-i{\bm{q}}\cdot{\bm{x}}}$ and the reciprocal vectors ${{\bm{e}}^{\alpha}}$ satisfying ${\bm{e}}^{\alpha}\cdot{\bm{e}}_{\beta}=\delta^{\alpha}_{\beta}$, the Kronecker delta (see Fig. 1c). Because $F_{\rm c}$ breaks the lateral continuous translational symmetry, modes with different wavevectors are coupled to one another. As shown in Appendix A, the modes coupled to a given $\bm{q}$ belong to the subset

[TABLE]

III Dynamic equations

Following Seifert and Langer Seifert , we regard each monolayer as a compressible 2D fluid having the shear viscosity $\mu$ and the bulk viscosity $\zeta$. The upper and lower monolayers can have different fluid velocities, ${\bm{v}}^{+}$ and ${\bm{v}}^{-}$, respectively (Fig 1a). The full dynamic equations consist of (i) lateral force balance for each monolayer, (ii) force balance normal to the bilayer, and (iii) the continuity equation for lipids in each monolayer.

We use the Stokes equation for the solvent velocity field ${\bm{V}}$ and the pressure field $p$, with the shear viscosity $\eta$:

[TABLE]

where $\hat{\nabla}=(\partial_{x},\partial_{y},\partial_{z})$ is the 3D nabla operator. No-slip boundary condition is employed at the membrane surface, $v^{\pm}_{i}=V_{i}\ (i=x,y)$ and $V_{z}=\partial h/\partial t$ at $z\to 0^{\pm}$. We also impose ${\bm{V}}\to 0$ and $p\to p_{0}$ as $z\to\pm\infty$. The 2D viscous stress tensors in the monolayers are given by

[TABLE]

where the superscript “$+$” (resp. “$-$”) denote the upper (resp. lower) monolayer. Then the lateral force balance equation in each monolayer reads

[TABLE]

where $T_{ij}^{+}$ (resp. $T_{ij}^{-}$) is the stress tensor $T_{ij}=-p\delta_{ij}+\eta(\partial_{i}V_{j}+\partial_{j}V_{i})$ in the solvent fluid evaluated at $z\to 0^{+}$ (resp. $z\to 0^{-}$). The last term is due to the inter-monolayer friction, with the friction coefficient $b$ Seifert . In the normal direction, the force exerted by the surrounding fluid is balanced with the restoring force of the membrane,

[TABLE]

At linear order in ${\bm{v}}^{\pm}$ and $\rho^{\pm}$, which are both considered to be small, the continuity equation in each monolayer is given by

[TABLE]

The velocities, ${\bm{V}}$ and ${\bm{v}}^{\pm}$, can be eliminated from the dynamic equations by integrating the Stokes equations along $z$ for each mode $\bm{q}$ (see Appendix B). This yields coupled linear equations for $\hat{h}\equiv h/d$ and $\rho\equiv(\rho^{+}-\rho^{-})/2$:

[TABLE]

where $c(q)=2b+2\eta q+(\mu+\zeta)q^{2}$ with $q=|\bm{q}|$. We have added $u_{h}$ and $u_{\rho}$, representing external forces applied mechanically (e.g., by active molecules) which act on the variables $\hat{h}$ and $\rho$, respectively. Since ${\cal F}_{\bm{q}}[\delta F/\delta h]$ includes $\hat{h}({\bm{q}}^{\prime})$ for $\forall{\bm{q}}^{\prime}\in Q_{\bm{q}}$, these equations actually consist of sets of coupled equations for the variables $\{\hat{h}({\bm{q}}^{\prime}),\rho({\bm{q}}^{\prime})\}$ in each set $Q_{\bm{q}}$ (see Appendix A).

Without the cytoskeleton, the modes for different wavevectors are not coupled in Eqs. (9) and (10). Then, $\hat{h}({\bm{q}})$ and $\rho({\bm{q}})$ exhibit two relaxation rates, $\gamma_{+}^{(0)}(q)>\gamma_{-}^{(0)}(q)$, associated with some linear combinations of $\hat{h}({\bm{q}})$ and $\rho({\bm{q}})$. Seifert and Langer discussed these relaxation modes for vanishing tension Seifert . They found a crossover wavenumber $q_{\rm c}=2\eta k/(b\tilde{\kappa})\approx 4.4\times 10^{6}\,{\rm m}^{-1}$, at which the relaxation behavior of the membrane changes qualitatively. Here we set $\kappa=2\times 10^{-20}\,{\rm J}$ as in JBPRL (the value of $\kappa$ measured in experiments lies in quite a wide range, $1$ to $30\times 10^{-20}\,{\rm J}$ Brochard ; BetzPNAS ; Turlier ; Popescu , but the following results remain almost unchanged even with these different values). For large scales satisfying $q\ll q_{\rm c}$, the rates correspond to $\rho$ relaxing quickly followed by $h$ relaxing slowly with $\rho$ being slaved JBNLM . For small scales, $q\gg q_{\rm c}$, conversely, they correspond to $h$ relaxing quickly followed by $\rho$ relaxing slowly with $h$ being slaved. Hence, the dynamics on the small scales is dominated by the inter-monolayer friction, whereas that on the large scales is dominated by the solvent viscosity. In the presence of tension, their results hold for $\sigma\ll\sigma_{\rm c}\equiv(2\eta k)^{2}/(\tilde{\kappa}b^{2})$, except at very large scales (see Refs. JBNLM ; OkamotoEPJE ; Sachin and Appendix C). However, for $\sigma\gtrsim\sigma_{\rm c}$ the dynamics is dominated at all scales by the inter-monolayer friction, with $\gamma_{+}^{(0)}\simeq(\sigma q+\kappa q^{3})/(4\eta)>\gamma_{-}^{(0)}\approx kq^{2}/(2b)$ JBNLM .

IV results

IV.1 Relaxation of a large-scale deformation

The cytoskeleton shifts the mode relaxation rates by an amount that depends on the prestress $\sim\!\nu$, and at the same times it couples all the modes belonging to a common set $Q_{\bm{q}}$. Let us first discuss how the rates of the large scale modes, with $q\ll q_{\rm c}\ll 2\pi/a$, are shifted by the cytoskeleton. For such modes, the dependence on the direction of $\bm{q}$ is negligible. In the following, analytical expressions will be given systematically at first-order in a perturbative expansion in power series of $\nu$ (see Appendix C). The parameter values used in the following numerical calculations are summarized in Table 1.

From the dynamic equations (9)–(10), we find that the rates of the large-scale modes, $\gamma_{+}>\gamma_{-}$, are shifted according to

[TABLE]

where $\sigma_{\rm eff}=\sigma+3^{1/2}\nu$ and $\kappa_{\rm eff}=\kappa-3^{1/2}\nu a^{2}/16$ are the tension and the bending rigidity renormalized by the cytoskeleton JBEPL ; Turlier . Note that the fast and slow rates have been exchanged with respect to their bare value $\gamma_{+}^{(0)}\simeq kq^{2}/(2b)>\gamma_{-}^{(0)}\simeq(\sigma q+\kappa q^{3})/(4\eta)$ because we anticipate $\sigma_{\mathrm{eff}}\gtrsim\sigma_{\rm c}>\sigma$ for large enough $\nu$. The other rates $\gamma_{\pm}(q^{\prime})$ of the subset ${\bm{q}}^{\prime}\in Q_{\bm{q}}$ associated with ${\bm{q}}$ are also shifted from their bare values $\gamma_{\pm}^{(0)}(q^{\prime})$; they correspond to wavelengths comparable to or smaller than the cytoskeleton mesh size and are much faster than those in Eq. (11). Our detailed calculations show that the shifts of these rates are small (see Appendix C), but this does not mean that the small-scale modes are not affected by the cytoskeleton

The relaxation of a large scale mode $\bm{q}$ excites all the small scale modes in $Q_{\bm{q}}$ (Fig. 2). To investigate this effect, we set the initial condition $\hat{h}({\bm{x}})=e^{i{\bm{q}}\cdot{\bm{x}}}$ and $\rho({\bm{x}})=0$ with $q=10^{6}\,\mathrm{m}^{-1}$ in the direction ${\bm{q}}/q=(\sqrt{3}/2,1/2)$, and we integrate numerically the dynamical equations up to the cutoff $20\pi/a$, of the order of the inverse membrane thickness. Experiments indicate $\sigma_{\mathrm{eff}}\approx\nu\approx 10^{-7}$–$10^{-5}\,\mathrm{N/m}$ with $\sigma$ very small or even negative Turlier ; Sackmann ; Garcia ; Popescu . Accordingly, besides the values already given, we set $\sigma=10^{-11}\,\mathrm{N/m}$ and $\nu=10^{-6}\,\mathrm{N/m}$, yielding $\kappa_{\rm eff}=1.9\times 10^{-20}\,\mathrm{J}$ and $\sigma_{\rm eff}=1.73\times 10^{-6}\,\mathrm{N/m}$. We study the coupled evolution of $\hat{h}$ and $\rho$ for ${\bm{q}}$ and for small-scale modes ${\bm{q}}^{\prime}\in Q_{\bm{q}}$. In Fig. 2, we present as an example only ${\bm{q}}^{\prime}={\bm{q}}+2\pi{\bm{e}}^{2}\simeq 2\pi{\bm{e}}^{2}$, as we find the other small-scale modes in $Q_{\bm{q}}$ also exhibit a similar behavior. In the short time interval $0<t\ll 1/\gamma_{\pm}(q^{\prime})\approx 1/\gamma_{\pm}^{(0)}(q^{\prime})$, the small-scale modes $\hat{h}({\bm{q}}^{\prime})$ and $\rho({\bm{q}}^{\prime})$, that are initially zero, are excited, while $\hat{h}({\bm{q}})$ almost remains unchanged (Fig. 2a and c). All the excited small-scale modes rapidly approach their respective quasi-equilibrium states $\hat{h}_{\rm qe}$ and $\rho_{\rm qe}$, which minimize the free energy for a fixed value of $\hat{h}({\bm{q}})$, given by

[TABLE]

where $K_{\bm{q}}=\sum_{\bm{n}}(1-e^{i{\bm{q}}\cdot{\bm{n}}})$. Figure 2b illustrates the long time evolution of the system. For $1/\gamma_{\pm}^{(0)}(q^{\prime})\ll t\lesssim\gamma_{+}^{-1}$, $\hat{h}({\bm{q}})$, $\hat{h}({\bm{q}}^{\prime})$ and $\rho({\bm{q}}^{\prime})$ decay with the common rate $\gamma_{+}$. Then, around $t\approx\gamma_{+}^{-1}$, $\hat{h}({\bm{q}})$ follows the dynamical quasi-equilibrium value $\hat{h}_{\rm qe}^{\rho}({\bm{q}};\rho({\bm{q}},t))$ that minimizes the free energy at fixed $\rho({\bm{q}},t)$. Finally, for $t\gg\gamma_{+}^{-1}$, all the modes decay with the common rate $\gamma_{-}$, with $\hat{h}({\bm{q}}^{\prime})$, $\rho({\bm{q}}^{\prime})$ and $\hat{h}({\bm{q}})$ following their respective dynamical quasi-equilibrium values $\hat{h}({\bm{q}}^{\prime})\simeq\hat{h}_{\rm qe}({\bm{q}}^{\prime};\hat{h}_{\rm qe}^{\rho}({\bm{q}}))$, $\rho({\bm{q}}^{\prime})\simeq\rho_{\rm qe}({\bm{q}}^{\prime};\hat{h}_{\rm qe}^{\rho}({\bm{q}}))$ and $\hat{h}({\bm{q}})\simeq\hat{h}_{\rm qe}^{\rho}({\bm{q}};\rho({\bm{q}},t))$ (Fig. 2b and c). Suppose in Fig. 2 the initial amplitude of $h({\bm{q}})$ is comparable with the mode wavelength, $2\pi/q\approx 6.3\,\mu{\rm m}$. Then the amplitude of the excited small-scale mode $h({\bm{q}}^{\prime})$ is about $2.4\times 10^{-5}\,\mu{\rm m}$, which is much smaller than the mode wavelength $2\pi/q^{\prime}\approx 8.6\times 10^{-2}\,\mu{\rm m}$, and may not be observable in experiments. This is because an energy cost to make a deformation with amplitude $q^{\prime-1}$ at the small-scale $q^{\prime}$ is larger than to make a deformation with amplitude $q^{-1}$ at the large-scale $q$. Nevertheless, we notice that the cytoskeleton alters qualitatively the large-scale dynamics; because of the cytoskeleton that yields $\sigma_{\mathrm{eff}}\sim\sigma_{\rm c}\approx 3\times 10^{-6}\,\mathrm{N/m}$, the slowest relaxation process is dominated by the large-scale compression mode $\rho({\bm{q}})$ limited by the inter-monolayer friction.

IV.2 Large-scale deformation induced by small-scale deformation via the cytoskeleton

We have seen that large-scale deformations excite the modes whose scales are comparable or smaller than the cytoskeleton mesh. Now a question arises. Can small scale deformations excite large scale ones? If yes, is the amplitude of the excited modes large enough to be observable? The answer is no, if there are no applied forces. This is because the small-scale modes are much faster than the large-scale ones, so that the small-scale modes rapidly relax before large-scale modes are excited. However, if we keep applying forces only to the small-scale modes for a time longer than the relaxation time of the large-scale modes, the latter will be excited via the cytoskeleton. Furthermore, if forces are applied to many small-scale modes, the amplitude of the excited large-scale mode can be noticeably large. In fresh RBCs, active molecules could be the source of these forces, as their characteristic time is of the order of $1\,\mathrm{s}$ BetzPNAS ; Turlier which is much larger than the typical relaxation time of the modes for $q\sim 10^{6}\,\mathrm{m}^{-1}$. It was further proposed that the active force is particularly enhanced on the scales of the cytoskeleton mesh ParkPNAS .

To study this, we choose again $q=10^{6}\,\mathrm{m}^{-1}$, oriented as before, and we apply a constant force $u_{h}({\bm{q}}_{1})=\bar{u}$ only to $\hat{h}({\bm{q}}_{1})$ with ${\bm{q}}_{1}={\bm{q}}+2\pi{\bm{e}}^{1}\in Q_{\bm{q}}$ (Fig. 3). The corresponding wavelengths of $q$ and ${\bm{q}}_{1}$ are then $2\pi/q\approx 6.3\,\mu{\rm m}$ and $2\pi/|{\bm{q}}_{1}|\approx 0.87a$, respectively. We investigate the response of $\hat{h}$ at the large scale $\bm{q}$ but also at another small scale, ${\bm{q}}_{2}={\bm{q}}+2\pi{\bm{e}}^{2}\in Q_{\bm{q}}$. With the initial condition $\hat{h}=\rho=0$, all the modes will be proportional to $\bar{u}$. For $t\lesssim 1/\gamma_{\pm}^{(0)}(q_{1})$, the small scale deformation $\hat{h}({\bm{q}}_{1})$ is excited and reaches the stationary value $\hat{h}_{\rm st}({\bm{q}}_{1})\simeq\bar{u}/(\kappa d^{2}q_{1}^{4})$ minimizing $F-(2\pi)^{-2}\bar{u}\hat{h}({\bm{q}}_{1})$. Then, for $1/\gamma_{\pm}^{(0)}(q_{1})\lesssim t\lesssim\gamma_{\pm}^{-1}$, the large-scale mode $\hat{h}({\bm{q}})$ gets excited by $\hat{h}({\bm{q}}_{1})$ via the cytoskeleton deformation, and then for $t\gtrsim\gamma_{\pm}^{-1}$, $\hat{h}({\bm{q}})$ reaches the stationary value

[TABLE]

where $\Delta\sigma_{\rm eff}=\sigma_{\rm eff}-\sigma$ and $\Delta\kappa_{\rm eff}=\kappa_{\rm eff}-\kappa$. With our choice of parameters, $(\Delta\sigma_{\rm eff}+\Delta\kappa_{\rm eff}q^{2})/(\sigma_{\rm eff}+\kappa_{\rm eff}q^{2})\simeq 1$, and thus $\hat{h}_{\rm st}({\bm{q}})\simeq-\hat{h}_{\rm st}({\bm{q}}_{1})$, consistent with Fig. 3a. We find that the other small-scale modes, such as $\hat{h}(\bm{q}_{2})$, are also excited, but not significantly (Fig. 3a).

When a force distribution is applied to multiple small-scale modes, the magnitude of the excited large-scale mode $\hat{h}({\bm{q}})$ can become much larger than in the case examined above. To show this, let us consider at each lattice site of the cytoskeleton active forces inducing some local curvature. Such forces can formally be derived by adding a “fictitious” potential $U=-\sum_{\ell}w_{\ell}\int{\rm d}^{2}x\,\delta({\bm{x}}-{\bm{R}}_{\ell})\nabla^{2}h$ to the free energy. This yields $u({\bm{q}}^{\prime})=q^{\prime 2}w_{{\bm{q}}}$, $\forall{\bm{q}}^{\prime}\in Q_{\bm{q}}$ with $w_{\bm{q}}=\sum_{\ell}w_{\ell}e^{-i{\bm{q}}\cdot{\bm{R}}_{\ell}}$. Note that some force is applied also to the large-scale mode ${\bm{q}}$. By linearity, the effect of the small-scale modes $\bm{q}^{\prime}$ on the large-scale mode $\bm{q}$, denoted by $\delta\hat{h}(\bm{q})$, will be enhanced by a factor $r=\sum^{\prime}_{m_{1},m_{2}}q_{1}^{2}/|2\pi m_{\alpha}{\bm{e}}^{\alpha}|^{2}$ with respect to the case shown in Fig. 3 where only the mode $\bm{q}_{1}$ was excited. The sum is taken up to the high wavevector cutoff while excluding $m_{1}=m_{2}=0$. With the parameters given above, we find $r\approx 20$. Assuming that the microscopic forces can produce a deformation of amplitude comparable to the mesh size, i.e., $h_{1}\simeq 50\,\mathrm{nm}$ BetzPNAS ; ParkPNAS , we expect the large scale response $\delta h$ to be about $rh_{1}$, which is a sizeable deformation of the order of $1\,\mu\mathrm{m}$. Note that this scale transfer of membrane deformation requires applying the small scale forces for at least about $10\,\mathrm{ms}$ (Fig. 3).

V Discussion and Summary

In this paper, since we assume only vertical motion of the quasi-planer RBC patch, and equilibrium dynamics, we neglect the dissipation due to tangential motion of the cytoskeleton as well as the cytoskeleton activity that was studied in Ref. Turlier . Nevertheless, in Appendix D, we have considered the friction between the tangential monolayer flow and the anchored proteins, yielding an extra contribution to the lateral force balance equation. However, in our detailed calculation, it is shown to be negligible. As for the viscous drag of the spectrin filaments due to the surrounding fluid, it was also shown to be negligible Turlier . For simplicity, we have assumed a quasi-planer membrane, i.e., small deformations about a flat reference shape. However, real RBCs are intrinsically curved objects that fluctuate about a curved reference shape Turlier ; EvansPNAS . For such cases, not only $\rho=(\rho^{+}-\rho^{-})/2$ but also $\bar{\rho}=(\rho^{+}+\rho^{-})/2$ is coupled to the membrane deformation $h$ Sachin . Furthermore, for a curved membrane, the tangential deformation of the cytoskeleton is also coupled to $h$ Turlier . For the relaxation of $\bar{\rho}$, we can show that the inter-monolayer friction is not a dissipation source, while the friction between the anchored proteins and the monolayers is one of the major dissipation sources for large scales satisfying $2\eta q+(\mu+\zeta)q^{2}\ll\lambda/a^{2}$. However, in the relaxation of $\bar{\rho}$ for a bilayer without the cytoskeleton, the inertia effect of the surrounding fluid can not be neglected Seifert , so that a more careful study is necessary in the future.

In summary, we have studied the dynamics of RBC membranes modelled as bilayers coupled to a pre-stressed discrete elastic network, and subject to viscous dissipation in the solvent, in each monolayer and between the monolayers. Given the mesh size of the cytoskeleton ($\approx\!100\,\mathrm{nm}$), it is important from the biological point of view to address the dynamics at scales both larger and smaller than the cytoskeleton. Because the latter breaks lateral translational symmetry, each mode is coupled to all the modes that are congruent modulo a wavevector of the cytoskeleton’s reciprocal lattice. We have characterized how the small modes renormalize the relaxation rates of the large modes. We have found that, because of the large renormalized tension $\sigma_{\rm eff}$, the shape relaxation dynamics on the large-scales is dominated by the inter-monolayer friction that has regularly been neglected in the previous theories on RBC dynamics GovPRL ; Turlier ; EvansPNAS . It has been also shown that applying forces on the small scale modes for a sufficiently long time can excite large scale deformations.

To the best of our knowledge, however, the correlations between different Fourier modes, $\langle h({\bm{q}},t)h({\bm{q}}^{\prime},t^{\prime})\rangle$, has not been measured in previous experiments. It is informative to measure this quantity in order to know the precise dynamical processes where modes in different scales are coupled due to the cytoskeleton, and also to understand the behavior of active forces in fresh RBCs.

Acknowledgements.

R.O. and S.K. acknowledge support from the Grant-in-Aid for Scientific Research on Innovative Areas “Fluctuation and Structure” (Grant No. 25103010) from the Ministry of Education, Culture, Sports, Science, and Technology of Japan, and from the Grant-in-Aid for Scientific Research (C) (Grant No. 15K05250) from the Japan Society for the Promotion of Science (JSPS).

Appendix A Coupling between different Fourier modes

Notice that the Fourier modes between different wavevectors are coupled with each other because the cytoskeleton network breaks the continuous translational symmetry. More precisely, while ${\cal F}_{\bm{q}}[\delta F/\delta\rho]=2k[\rho({\bm{q}})-dq^{2}h({\bm{q}})]$ in Eq. (28) does not couple different Fourier modes, ${\cal F}_{\bm{q}}[\delta F/\delta h]$ in Eq. (26) couples the Fourier modes in the set $Q_{\bm{q}}=\{{\bm{q}}+2\pi m_{\alpha}{\bm{e}}^{\alpha}|m_{\alpha}\in\mathbb{Z}\}$. To see this, we rewrite Eq. (2) as

[TABLE]

Then we can calculate the Fourier-transformed functional derivative as

[TABLE]

Now we use the identity

[TABLE]

where $g=|{\bm{e}}_{\alpha}\cdot{\bm{e}}_{\beta}|$ is the determinant of the metric tensor in the primitive-vector-frame $\{{\bm{e}}_{\alpha}\}$, and is given by $g=3a^{4}/4$. Using this identity, we can rewrite Eq. (15) as

[TABLE]

where $K_{\bm{q}}=\sum_{\bm{n}}(1-e^{i{\bm{q}}\cdot{\bm{n}}})$. We can clearly see from Eq. (17) that in Eq. (9) the Fourier modes in the common set $Q_{\bm{q}}$ are coupled to each other. Note that without the cytoskeleton $\nu=0$ (and without the force terms $u_{h}$ and $u_{\rho}$), Eqs. (9) and (10) reduce to the equations studied by Seifert and Langer Seifert .

Appendix B Elimination of the velocity fields

In accordance with the in-plane Fourier transform with a wavevector ${\bm{q}}$, we introduce ${\bm{q}}_{\perp}={\bm{e}}_{z}\times{\bm{q}}$, where ${\bm{e}}_{z}$ is the unit vector pointing towards the $z$-direction. Then unit vectors $\hat{\bm{q}}$ and $\hat{\bm{q}}_{\perp}$ are defined as $\hat{\bm{q}}={\bm{q}}/q$ and $\hat{\bm{q}}_{\perp}={\bm{q}}_{\perp}/q={\bm{e}}_{z}\times\hat{\bm{q}}$, respectively. The logitudinal and transverse components of the Fourier transform of ${\bm{V}}$ and ${\bm{v}}^{\pm}$ are defined by $V_{\parallel}({\bm{q}},z)=\hat{\bm{q}}\cdot{\bm{V}}({\bm{q}},z)$, $V_{\perp}({\bm{q}},z)=\hat{\bm{q}}_{\perp}\cdot{\bm{V}}({\bm{q}},z)$, $v^{\pm}_{\parallel}({\bm{q}})=\hat{\bm{q}}\cdot{\bm{v}}^{\pm}({\bm{q}})$, and $v^{\pm}_{\perp}({\bm{q}})=\hat{\bm{q}}_{\perp}\cdot{\bm{v}}^{\pm}({\bm{q}})$. The Fourier transform of Eq. (4) is written as

[TABLE]

These equations are solved to obtain

[TABLE]

where the upper and the lower signs indicate the solutions for $z>0$ and for $z<0$, respectively.

Substituting Eqs. (22) and (24) into the Fourier transform of Eq. (7), we obtain

[TABLE]

which yields Eq. (9) (without the force term $u_{h}$). Next, we use Eq. (25) to eliminate $V_{\parallel}$ and $V_{z}$ from the longitudinal component of Eq. (6), and obtain

[TABLE]

The Fourier transform of Eq. (8) relates $\partial\rho/\partial t=(\partial/\partial t)(\rho^{+}-\rho^{-})/2$ with the longitudinal velocity $v_{\parallel}$ as $\partial\rho/\partial t=-iq(v_{\parallel}^{+}-v_{\parallel}^{-})/2$. Then Eq. (27) yields

[TABLE]

where $c(q)=2b+2\eta q+(\mu+\zeta)q^{2}$. Using the identity $\delta(\cdots)/\delta\rho=\delta(\cdots)/\delta\rho^{+}-\delta(\cdots)/\delta\rho^{-}$, we obtain Eq. (10) (without the force term $u_{\rho}$).

Appendix C Operator representation and perturbation expansion

We seek the relaxation rates and their associated eigen modes of Eqs. (9) and (10) as a power series of $\nu$, since for $\nu=0$ they can be obtained analytically. It is convenient to introduce for each $Q_{\bm{q}}$ a Hilbert space $S_{\bm{q}}$ spanned by an orthonormal set \big{\{}|h;{\bm{q}}^{\prime}\rangle,\ |\rho;{\bm{q}}^{\prime}\rangle\ |\ {\bm{q}}^{\prime}\in Q_{\bm{q}}\big{\}}, and its dual space $S_{\bm{q}}^{*}$ spanned by an orthonormal set \big{\{}\langle h;{\bm{q}}^{\prime}|,\ \langle\rho;{\bm{q}}^{\prime}|\ |\ {\bm{q}}^{\prime}\in Q_{\bm{q}}\big{\}}. We define the following state vector in $S_{\bm{q}}$,

[TABLE]

where $\hat{h}=h/d$. Then Eqs. (26) and (28) (or Eqs. (9) and (10) without the force terms $u_{h}$ and $u_{\rho}$) are written as

[TABLE]

where the linear operators $\Gamma_{\rm SL}$ and $\Gamma_{h}$ are defined by

[TABLE]

Notice that, in the absence of the cytoskeleton ($\nu=0$), Eq. (30) reduces to $0=|\dot{\Psi}\rangle_{\bm{q}}+\Gamma_{\rm SL}|\Psi\rangle_{\bm{q}}$ which was discussed by Seifert and Langer for vanishing tension $\sigma=0$ Seifert .

Let us consider the eigenvalue problem

[TABLE]

where $\gamma$ and $|e\rangle$ are the eigenvalue and the eigenvector of $\Gamma_{\rm SL}+\nu\Gamma_{h}$, respectively. We expand $\gamma$ and $|e\rangle$ in powers of $\nu$ as $\gamma=\gamma^{(0)}+\nu\gamma^{(1)}+\cdots$ and $|e\rangle=|e^{(0)}\rangle+\nu|e^{(1)}\rangle+\cdots$. The 0th and the 1st order equations read

[TABLE]

For the following perturbation calculation, it is convenient to introduce the 0th and 1st order left eigenvectors $\langle e^{(0)\dagger}|$ and $\langle e^{(1)\dagger}|$, respectively. These satisfy

[TABLE]

Here $\gamma^{(0)}$ and $\gamma^{(1)}$ are common to Eqs. (34) and (35), respectively. Note that $\langle e^{(0)\dagger}|$ (resp. $\langle e^{(1)\dagger}|$) is not the Hermitian conjugate of $|e^{(0)}\rangle$ (resp. $|e^{(1)}\rangle$), because neither $\Gamma_{\rm SL}$ nor $\Gamma_{h}$ is a Hermitian operator.

C.1 0th order

Since $\Gamma_{\rm SL}$ does not couple the Fourier modes of different wavevectors, we readily obtain the 0th order eigenvalues,

[TABLE]

with $g(c)=(\sigma+\tilde{\kappa}q^{2})c(q)-4\eta kq$. Their associated right eigenvectors are given by

[TABLE]

where $e_{\pm}(q)=[g(q)\mp\sqrt{g(q)^{2}+32\eta k^{2}d^{2}q^{3}c(q)}]/[4kc(q)]$. The corresponding left eigenvectors are

[TABLE]

with $e_{\pm}^{\dagger}(q)=c(q)e_{\pm}(q)/(2\eta d^{2}q^{3})$. In the above, the left eigenvector is normalized such that $\langle e_{\pm}^{(0)\dagger};{\bm{q}}|e_{\pm}^{(0)};{\bm{q}}\rangle=1$. Since a contraction of left and right eigenvectors associated with different eigenvalues vanishes, we obtain

[TABLE]

This yields $e_{+}^{\dagger}(q)e_{-}(q)=e_{-}^{\dagger}(q)e_{+}(q)=-1$, which can also be confirmed directly from the definition of $e_{\pm}$ and $e^{\dagger}_{\pm}$.

In Table 2, we present approximate expressions of $\gamma_{\pm}^{(0)}$ and $e_{\pm}$ for different length scales classified by the characteristic wavenumbers $b\sigma/(2k\eta)$ and $q_{\rm c}=2\eta k/(b\tilde{\kappa})$, where a small bare tension $\sigma\ll\sigma_{\rm c}$ is assumed JBNLM . Typical parameter values quoted in the main text yield $\sigma_{\rm c}\equiv(2\eta k)^{2}/(\tilde{\kappa}b^{2})\approx 3\times 10^{-6}\,\mathrm{J/m}^{2}$. The behavior in the two regimes (i) and (ii) in Table 2 is essentially the same as those in Ref. Seifert for $\sigma=0$. In the presence of tension, there appears another regime $q\ll b\sigma/(2k\eta)$ where the dynamics is again dominated by the inter-monolayer friction JBNLM . For $\sigma=10^{-11}{\rm N/m}$ chosen in the main text, the characteristic wavenumber corresponds to the length $4\pi k\eta/(b\sigma)\approx 0.44\,{\rm m}$, which is too large to be measured in experiments. However, with $\sigma=4\times 10^{-7}{\rm N/m}$, for instance, we have $4\pi k\eta/(b\sigma)\approx 10\mu{\rm m}$, which is relevant for giant unilamellar vesicles. For large tension $\sigma\gtrsim\sigma_{\rm c}$, on the other hand, the dynamics is dominated by the inter-monolayer friction in all length scales JBNLM .

C.2 1st order

We assume for simplicity that the 0th order eigenvalue $\gamma_{\pm}^{(0)}(q)$ is not degenerated in the Hilbert space $S_{\bm{q}}$. This assumption is always valid for scales much larger than the lattice spacing of the cytoskeleton, $q\ll a^{-1}$. Equation (35) yields the 1st order correction $\nu\gamma_{\pm}^{(1)}(q)$ to the 0th order eigenvalue $\gamma_{\pm}^{(0)}(q)$ as

[TABLE]

To obtain the shifted eigenvectors, we expand the 1st order eigenvectors in terms of the 0th order eigenvectors,

[TABLE]

Here we set $s^{\pm}_{\pm}({\bm{q}};{\bm{q}})=0$ (as for the perturbation theory in quantum mechanics). Operation of $\langle e_{\epsilon}^{(0)\dagger};{\bm{q}}^{\prime}|$ to the both sides of Eq. (35) yields

[TABLE]

Similarly, we can calculate the 1st order left eigenvector $\langle e_{\pm}^{(1)\dagger};{\bm{q}}|$ from Eq. (37). One can also show that Eq. (41) is generalized to the 1st order in $\nu$ as

[TABLE]

Case (i) $q\ll 2\eta k/(b\tilde{\kappa})$

Let us suppose a wavevector ${\bm{q}}_{\rm L}$ satisfies $q_{\rm L}=|{\bm{q}}_{\rm L}|\ll 2\eta k/(b\tilde{\kappa})$, i.e., case (i) in Table 2. We then study how the cytoskeleton alters the rates $\gamma_{\pm}^{(0)}(q_{\rm L})$. We notice that the 0th order eigenvalue $\gamma_{\pm}^{(0)}(q_{\rm L})$ is not degenerated in $S_{{\bm{q}}_{\rm L}}$, because, for such a small wavevector, there is no other wavevector ${\bm{q}}^{\prime}\in Q_{{\bm{q}}_{\rm L}}$ satisfying $|{\bm{q}}^{\prime}|=q_{\rm L}$. Using Eq. (42) and approximate expressions in Table 2, we find

[TABLE]

when $q_{\rm L}\ll 2\eta k/(b\tilde{\kappa})$. In general, $K({\bm{q}}_{\rm L})$ depends on the direction of the wavevector $\hat{\bm{q}}_{\rm L}$. However, it is approximated as

[TABLE]

for $q_{\rm L}\ll a^{-1}$, which is isotropic. We can always use this approximation in Eq. (46) which is valid for $q_{\rm L}\ll 2\eta k/(b\tilde{\kappa})$. This is because $2\eta k/(b\tilde{\kappa})=4.37\times 10^{6}$ ${\mathrm{m}}^{-1}$ for typical parameter values chosen in the main text, and it is smaller than the reciprocal of the lattice spacing $a\sim 10^{-7}$ ${\mathrm{m}}$ of the cytoskeleton. Substitution of Eq. (47) into Eq. (46) yields Eq. (11). Note for $\nu=10^{-6}\,{\rm N/m}$ chosen in the main text, $\gamma^{(0)}_{-}+\nu\gamma^{(1)}_{-}$ is larger than $\gamma^{(0)}_{+}+\nu\gamma^{(1)}_{+}$, and therefore by definition, $\gamma_{+}\simeq\gamma^{(0)}_{-}+\nu\gamma^{(1)}_{-}$ and $\gamma_{-}=\gamma^{(0)}_{+}+\nu\gamma^{(1)}_{+}$.

Next we discuss the shifted eigenvector $|e_{\pm}^{(1)};{\bm{q}}_{\rm L}\rangle$, $q_{\rm L}\ll 2\eta k/(b\tilde{\kappa})$. In the present case of $q_{\rm L}\ll 2\eta k/(b\tilde{\kappa})$, we can assume $q^{\prime}\gtrsim 2\pi/a\gg 2\eta k/(b\tilde{\kappa})$ for $\forall{\bm{q}}^{\prime}\in Q_{{\bm{q}}_{\rm L}}\setminus\{{\bm{q}}_{\rm L}\}$ (here “$\setminus$” indicates set difference). Thus, for ${\bm{q}}^{\prime}\neq{\bm{q}}_{\rm L}$, we can set $\gamma_{\pm}^{(0)}(q_{\rm L})-\gamma_{\epsilon}^{(0)}(q^{\prime})\simeq-\gamma_{\epsilon}^{(0)}(q^{\prime})$ in the denominator of Eq. (44). Using the expressions in (ii) of Table 2, we obtain for ${\bm{q}}^{\prime}\neq{\bm{q}}_{\rm L}$,

[TABLE]

For ${\bm{q}}^{\prime}={\bm{q}}_{\rm L}$, we can assume $\gamma_{+}^{(0)}(q_{\rm L})\gg\gamma_{-}^{(0)}(q_{\rm L})$ in the denominator of Eq. (44) and obtain

[TABLE]

where $E({\bm{q}}_{\rm L})=K_{{\bm{q}}_{\rm L}}/[4\sqrt{g}\eta q_{\rm L}\gamma_{+}^{(0)}(q_{\rm L})]\simeq bK_{{\bm{q}}_{\rm L}}/[2\sqrt{g}\eta kq_{\rm L}^{3}]$. To make the physical meaning of the shifted eigenvectors clear, we examine the time evolution of the vector

[TABLE]

We may suppose that the modes for $\forall{\bm{q}}^{\prime}\in Q_{{\bm{q}}_{\rm L}}\setminus\{{\bm{q}}_{\rm L}\}$ decay much faster than the modes for ${\bm{q}}_{\rm L}$, because in the 0th order $\gamma_{\pm}^{(0)}(q^{\prime})$ is much larger than $\gamma_{\pm}^{(0)}(q_{\rm L})$. Then, after sufficiently large time $t$ satisfying $t\gg\gamma_{\pm}^{(0)}(q^{\prime})^{-1}$ for $\forall{\bm{q}}^{\prime}\in Q_{{\bm{q}}_{\rm L}}\setminus\{{\bm{q}}_{\rm L}\}$, the modes for $\forall{\bm{q}}^{\prime}\in Q_{{\bm{q}}_{\rm L}}\setminus\{{\bm{q}}_{\rm L}\}$ are in the quasi-equilibrium state $\hat{h}_{\mathrm{qe}}({\bm{q}}^{\prime};\hat{h}({\bm{q}}_{\rm L}))$ and $\rho_{\mathrm{qe}}({\bm{q}}^{\prime};\hat{h}({\bm{q}}_{\rm L}))$ determined by ${\cal F}_{{\bm{q}}^{\prime}}[\delta F/\delta h]={\cal F}_{{\bm{q}}^{\prime}}[\delta F/\delta\rho^{\pm}]=0$ for a given value of $\hat{h}({\bm{q}}_{\rm L})$. Then we may approximate Eq. (50) as

[TABLE]

This will be justified later (see discussion around Eq. (58)). Using ${\cal F}_{{\bm{q}}^{\prime}}[\delta F/\delta\rho]=2k[\rho({\bm{q}}^{\prime})-dq^{\prime 2}h({\bm{q}}^{\prime})]$ and Eq. (17), we obtain Eq. (12) at linear order in $\nu$. We can then show that the following relation holds up to the first order in $\nu$:

[TABLE]

where

[TABLE]

and

[TABLE]

with $\bar{E}({\bm{q}}_{\rm L})=2\eta d^{2}q_{\rm L}^{3}E({\bm{q}}_{\rm L})/c(q_{\rm L})$ and $e_{\pm}=e_{\pm}(q_{\rm L})$. Therefore, after sufficiently large time $t$, $A_{+}$ and $A_{-}$ decay with the rates in Eq. (11). The above discussion indicates that the rates in Eq. (11) and their associated eigenvectors correspond to the relaxation of the (long wavelength) modes of ${\bm{q}}_{\rm L}$ accompanied by instantaneous relaxation of the other (short wavelength) modes of $\forall{\bm{q}}^{\prime}\in Q_{{\bm{q}}_{\rm L}}\setminus\{{\bm{q}}_{\rm L}\}$ to the quasi-equilibrium state. Note that even though at initial time we set $\hat{h}({\bm{q}}^{\prime})=\rho({\bm{q}}^{\prime})=0$ for ${\bm{q}}^{\prime}\in Q_{{\bm{q}}_{\rm L}}\setminus\{{\bm{q}}_{\rm L}\}$, after sufficiently large time ($t\gg\gamma_{\pm}^{(0)}(q^{\prime})^{-1}$) they are excited by non-zero $\hat{h}({\bm{q}}_{\rm L})$ to their quasi-equilibrium values in Eq. (12).

Case (ii) $q\gg 2\eta k/(b\tilde{\kappa})$

Next we discuss how the rate $\gamma_{\pm}^{(0)}(q_{\rm S})$ with $q_{\rm S}=|{\bm{q}}_{\rm S}|\gg 2\eta k/(b\tilde{\kappa})$ is altered by the cytoskeleton. In this regime, using the expressions in (ii) of Table 2, we find

[TABLE]

We also examine the magnitude of $\gamma_{\pm}^{(1)}$ compared with that of $\gamma_{\pm}^{(0)}$. Using the approximate expressions in (ii) of Table 2, we obtain

[TABLE]

Therefore the effect of cytoskeleton is negligible if $(aq_{\rm S})^{4}\gg 2a^{2}\nu K_{{\bm{q}}_{\rm S}}/(\sqrt{3}\tilde{\kappa})\approx 0.0722K_{{\bm{q}}_{\rm S}}$. Since $K_{\bm{q}}\lesssim 10$, the correction to the rates due to the cytoskeleton is not very large. However, the associated modes do not necessarily decay to zero for $t\gg 1/\gamma_{\pm}^{(0)}(q_{\rm S})$. In fact, if a large scale Fourier mode exists in $Q_{{\bm{q}}_{\rm S}}$, i.e., ${\bm{q}}_{\rm L}\in Q_{{\bm{q}}_{\rm S}}$ (or equivalently ${\bm{q}}_{\rm S}\in Q_{{\bm{q}}_{\rm L}}$), the Fourier modes for ${\bm{q}}_{\rm S}$ are rather excited to the quasi-equilibrium state by the large scale Fourier mode $\hat{h}({\bm{q}}_{\rm L})$, as discussed above. We can confirm this explicitly by operating $\langle e_{\pm}^{(0)\dagger};{\bm{q}}_{\rm S}|+\nu\langle e_{\pm}^{(1)\dagger};{\bm{q}}_{\rm S}|$ to Eq. (50). Owing to Eq. (45), the resultant quantity decays with the rate $\gamma_{\pm}^{(0)}(q_{\rm S})+\nu\gamma_{\pm}^{(1)}({\bm{q}}_{\rm S})\sim\gamma_{\pm}^{(0)}(q_{\rm S})$. Hence, for $t\gg\gamma_{\pm}^{(0)}(q_{\rm S})^{-1}$, we have

[TABLE]

Using the approximations in Table 2, we can show that Eq. (58) is equivalent to the quasi-equilibrium condition Eq. (12) with ${\bm{q}}^{\prime}={\bm{q}}_{\rm S}$ and ${\bm{q}}={\bm{q}}_{\rm L}$.

Appendix D Effects of friction between the bilayer and the proteins at the cytoskeleton vertices

We can take into account the friction between the bilayer and the proteins at the vertices of the cytoskeleton by modifying the lateral force balance Eq. (6) to

[TABLE]

In the above, the friction coefficient $\lambda$ is assumed to be common for both the upper and lower monolayers. The new term $\lambda\sum_{\ell}v_{i}^{\pm}\delta({\bm{x}}-{\bm{R}}_{\ell})$ also breaks the translational symmetry, and hence leads to the coupling between different Fourier modes. With this new term, Eq. (28) is modified to

[TABLE]

With the use of the identity Eq. (16), the new term is rewritten as

[TABLE]

where $w_{\perp}=-iq(v_{\perp}^{+}-v_{\perp}^{-})/2$. To eliminate $w_{\perp}$ from Eq. (60), we need the transverse part of Eq. (59). As in a similar way to derive Eq. (60), we obtain

[TABLE]

with $c_{w}(q)=2b+\eta q+\mu q^{2}$. Equations (26), (60) and (62) are the complete set of the relaxation equations for $\rho$ and $\hat{h}$. We can see from Eq. (61) that the Fourier modes for ${\bm{q}}$ are coupled to the Fourier modes for $\forall{\bm{q}}^{\prime}\in Q_{\bm{q}}$, as in the case without the friction at the vertices of the cytoskeleton.

The full equations are also represented in terms of operators and vectors in the Hilbert space $S_{\bm{q}}$, as in the previous section. To perform perturbation calculations, we regard both $\nu$ and $\lambda$ as small parameters. Then, to the first order in $\nu$ and $\lambda$, the governing equation is

[TABLE]

where the operator $\Gamma_{\rho}$ is defined as

[TABLE]

The correction to the rate $\gamma_{\pm}^{(0)}(q)$ due to the friction is then given by

[TABLE]

To measure the relevance of $D_{\pm}(q)$, we shall consider the ratio $|D_{\pm}/\gamma_{\pm}^{(0)}|$. Using the approximations in (i) $b\sigma/(2k\eta)\ll q\ll 2\eta k/(b\tilde{\kappa})$ of Table 2, we obtain

[TABLE]

Similarly, with the approximations in (ii) $q\gg 2\eta k/(b\tilde{\kappa})$ in Table 2, the ratio $|D_{\pm}/\gamma_{\pm}^{(0)}|$ is comparable or smaller than $\lambda/(ba^{2})$. Therefore the friction due to the network is negligible as long as the friction coefficient per area $\lambda/a^{2}$ is much smaller than the coefficient $b$ for the inter-monolayer friction. We use Saffman–Delbrück theory to estimate the value of $\lambda$ SD . In this theory two physical situations are examined; (i) membrane of finite size with surrounding fluid being neglected, and (ii) membrane of infinite size with surrounding fluid being taken into account. Since $\lambda$ is a “bare” friction constant between a protein and a monolayer with surrounding fluid being neglected, we may use the result for (i). Then we can set $\lambda\simeq 4\pi\mu/\ln(L/r_{0})$, with $L$ the membrane size and $r_{0}$ the protein size. Using $L/r_{0}\simeq 10^{2}$, $b\simeq 2\times 10^{8}\,\mathrm{J\,s}/\mathrm{m}^{4}$, $\mu\simeq 2\times 10^{-9}\,\mathrm{J\,s}/\mathrm{m}^{2}$ and $a\simeq 10^{-7}$ m, we estimate $\lambda/(a^{2}b)\simeq 5\times 10^{-3}$. We can thus neglect the effects of $\lambda$. Note that the estimation with (ii) also leads to the same conclusion; for (ii), we set $\lambda\simeq 4\pi\mu/\ln[\mu/(\eta r_{0})]$, and obtain $\lambda/(a^{2}b)\simeq 1.8\times 10^{-3}$ when $r_{0}\simeq 2\times 10^{-9}\,\mathrm{m}$.