ABSTRACT A significative fraction of all massive stars in the Milky Way move supersonically through their local interstellar medium (ISM), producing bow shock nebulae by wind-ISM interaction. The stability of these observed astrospheres around cool massive stars challenges precedent 2D (magneto-)hydrodynamical (MHD) simulations of their surroundings. We present 3D MHD simulations of the circumstellar medium of runaway M-type red supergiant stars moving with velocity $v_{\star }=50\, \rm km\, \rm s^{-1}$. We treat the stellar wind with a Parker spiral and assume a $7\, \rm \mu G$ magnetization of the ISM. Our free parameter is the angle θmag between ISM flow and magnetization, taken to 0°, 45°, and 90°. It is found that simulation dimension, coordinate systems, and grid effects can greatly affect the development of the modelled astrospheres. Nevertheless, as soon as the ISM flow and magnetization directions differs by more than a few degrees (θmag ≥ 5°), the bow shock is stabilized, most clumpiness and ragged structures vanishing. The complex shape of the bow shocks induce important projection effects, e.g. at optical H α line, producing complex of astrospheric morphologies. We speculate that those effects are also at work around earlier-type massive stars, which would explain their diversity of their observed arc-like nebula around runaway OB stars. Our 3D MHD models are fitting well observations of the astrospheres of several runaway red supergiant stars. The results interpret the smoothed astrosphere of IRC-10414 and Betelgeuse (αOri) are stabilized by an organized non-parallel ambient magnetic field. Our findings suggest that IRC-10414 is currently in a steady state of its evolution, and that Betelgeuse’s bar is of interstellar origin.