The empirical harmonic potential function of elastic network models (ENMs) is augmented by three- and four-body interactions as well as by a parameter-free connection rule. In the new bend-twist-stretch (BTS) model the complexity of the parametrization is shifted from the spatial level of detail to the potential function, enabling an arbitrary coarse graining of the network. Compared to distance cutoff-based Hookean springs, the approach yields a more stable parametrization of coarse-grained ENMs for biomolecular dynamics. Traditional ENMs give rise to unbounded zero-frequency vibrations when (pseudo)atoms are connected to fewer than three neighbors. A large cutoff is therefore chosen in an ENM (about twice the average nearest-neighbor distance), resulting in many false-positive connections that reduce the spatial detail that can be resolved. More importantly, the required three-neighbor connectedness also limits the coarse graining, i.e., the network must be dense, even in the case of low-resolution structures that exhibit few spatial features. The new BTS model achieves such coarse graining by extending the ENM potential to include three-and four-atom interactions (bending and twisting, respectively) in addition to the traditional two-atom stretching. Thus, the BTS model enables reliable modeling of any three-dimensional graph irrespective of the atom connectedness. The additional potential terms were parametrized using continuum elastic theory of elastic rods, and the distance cutoff was replaced by a competitive Hebb connection rule, setting all free parameters in the model. We validate the approach on a carbon-alpha representation of adenylate kinase and illustrate its use with electron microscopy maps of E. coli RNA polymerase, E. coli ribosome, and eukaryotic chaperonin containing T-complex polypeptide 1, which were difficult to model with traditional ENMs. For adenylate kinase, we find excellent reproduction (>90% overlap) of the ENM modes and B factors when BTS is applied to the carbon-alpha representation as well as to coarser descriptions. For the volumetric maps, coarse BTS yields similar motions (70%-90% overlap) to those obtained from significantly denser representations with ENM. Our Python-based algorithms of ENM and BTS implementations are freely available.