As a reduced representation of the nonlinear spectral fluxes of ideal invariants in incompressible magnetohydrodynamics, we construct a gradient-diffusion network model that combines phenomenological considerations and geometrical analysis of the exact nonlinear energy transfer function. The reduced-order representation of the conservative spectral transport of energy and cross-helicity is of port-Hamiltonian form, which highlights the flexibility and modularity of this approach. Numerical experiments with Reynolds numbers up to $~10^6$ yield clear power-law signatures of inertial-range energy spectra. Depending on the dominant timescale of energy transfer, Kolmogorov (-5/3), weak-turbulence (-2), or Iroshnikov-Kraichnan-like (-3/2) scaling exponents are observed. Anisotropic turbulence in a mean magnetic field is successfully modelled as well. The characteristic exponents of turbulence decay and the observed influence of cross-helicity on the energy transfer are consistent with the literature and in agreement with theoretical results.