Shell models of hydrodynamic turbulence originated in the seventies. Their main aim was to describe the statistics of homogeneous and isotropic turbulence in spectral space, using a simple set of ordinary differential equations. In the eighties, shell models of magnetohydrodynamic (MHD) turbulence emerged based on the same principles as their hydrodynamic counter-part but also incorporating interactions between magnetic and velocity fields. In recent years, significant improvements have been made such as the inclusion of non-local interactions and appropriate definitions for helicities. Though shell models cannot account for the spatial complexity of MHD turbulence, their dynamics are not over simplified and do reflect those of real MHD turbulence including intermittency or chaotic reversals of large-scale modes. Furthermore, these models use realistic values for dimensionless parameters (high kinetic and magnetic Reynolds numbers, low or high magnetic Prandtl number) allowing extended inertial range and accurate dissipation rate. Using modern computers it is difficult to attain an inertial range of three decades with direct numerical simulations, whereas eight are possible using shell models.In this review we set up a general mathematical framework allowing the description of any MHD shell model. The variety of the latter, with their advantages and weaknesses, is introduced. Finally we consider a number of applications, dealing with free-decaying MHD turbulence, dynamo action, Alfvén waves and the Hall effect. Summary and outlook 65Appendix A L1-models 66 2 Appendix B L2-models 67Appendix C N1-models 68 Appendix D N2-models 68Appendix E Numerical aspects 683 Variables u, b, p, z ± , a, u ± , b ±

We solve the Navier-Stokes equations with two simultaneous forcings. One forcing is applied at a given large scale and it injects energy. The other forcing is applied at all scales belonging to the inertial range and it injects helicity. In this way we can vary the degree of turbulence helicity from nonhelical to maximally helical. We find that increasing the rate of helicity injection does not change the energy flux. On the other hand, the level of total energy is strongly increased and the energy spectrum gets steeper. The energy spectrum spans from a Kolmogorov scaling law k^{-5/3} for a nonhelical turbulence, to a non-Kolmogorov scaling law k^{-7/3} for a maximally helical turbulence. In the latter case we find that the characteristic time of the turbulence is not the turnover time but a time based on the helicity injection rate. We also analyze the results in terms of helical modes decomposition. For a maximally helical turbulence one type of helical mode is found to be much more energetic than the other one, by several orders of magnitude. The energy cascade of the most energetic type of helical mode results from the sum of two fluxes. One flux is negative and can be understood in terms of a decimated model. This negative flux, however, is not sufficient to lead an inverse energy cascade. Indeed, the other flux involving the least energetic type of helical mode is positive and the largest. The least energetic type of helical mode is then essential and cannot be neglected.