To quantitatively understand the correlation between the molecular structure of an erythrocyte (red blood cell, RBC) and its mechanical response, and to predict mechanically induced structural remodelling in physiological conditions, we developed a computational model by coupling a multiscale approach of RBC membranes with a boundary element method (BEM) for surrounding Stokes flows. The membrane is depicted at three levels: in the whole cell level, a finite element method (FEM) is employed to model the lipid bilayer and the cytoskeleton as two distinct layers of continuum shells. The mechanical properties of the cytoskeleton are obtained from a molecular-detailed model of the junctional complex. The spectrin, a major protein of the cytoskeleton, is simulated using a molecular-based constitutive model. The BEM model is coupled with the FEM model through a staggered coupling algorithm. Using this technique, we first simulated RBC dynamics in capillary flow and found that the protein density variation and bilayer–skeleton interaction forces are much lower than those in micropipette aspiration, and the maximum interaction force occurs at the trailing edge. Then we investigated mechanical responses of RBCs in shear flow during tumbling, tank-treading and swinging motions. The dependencies of tank-treading frequency on the blood plasma viscosity and the membrane viscosity we found match well with benchmark data. The simulation results show that during tank-treading the protein density variation is insignificant for healthy erythrocytes, but significant for cells with a smaller bilayer–skeleton friction coefficient, which may be the case in hereditary spherocytosis.