Graphs provide a powerful abstraction for modeling real-world systems, such as social and transportation networks. Of particular significance is the study of diffusion processes on graphs, which is crucial to disciplines spanning biology, engineering and computer science, and to the understanding of phenomena such as disease propagation or the flux of goods and/or people through a transportation network. In this work, we formulate the solution of the diffusion equation on graphs as a parametric model, which gives us insight into the behavior of the diffusion processes by studying how the diffusivity parameters influence the model output. However, accurately simulating the model output of these diffusion processes can be computationally demanding since it involves solving large systems of ordinary differential equations. To address this challenge, we propose to construct surrogate models able to approximate the state of a graph at a given time from the knowledge of the diffusivity parameters. In particular, we consider recently introduced techniques from high-dimensional approximation based on sparse polynomial expansions, which are known to produce accurate and sample-efficient approximations when the function to be approximated has holomorphic regularity. Hence, to justify our methodology, we present theoretical results showing that solution maps resulting from a certain class of parametric graph diffusion processes are holomorphic. Then, we demonstrate numerically that it is possible to efficiently compute accurate sparse polynomial surrogate models from a few random samples, hence empirically showing the validity of our approach.