The fractional differential equation $L^βu = f$ posed on a compact metric graph is considered, where $β>0$ and $L = κ^2 - \nabla(a\nabla)$ is a second-order elliptic operator equipped with certain vertex conditions and sufficiently smooth and positive coefficients $κ, a$. We demonstrate the existence of a unique solution for a general class of vertex conditions and derive the regularity of the solution in the specific case of Kirchhoff vertex conditions. These results are extended to the stochastic setting when $f$ is replaced by Gaussian white noise. For the deterministic and stochastic settings under generalized Kirchhoff vertex conditions, we propose a numerical solution based on a finite element approximation combined with a rational approximation of the fractional power $L^{-β}$. For the resulting approximation, the strong error is analyzed in the deterministic case, and the strong mean squared error as well as the $L_2(Γ\times Γ)$-error of the covariance function of the solution are analyzed in the stochastic setting. Explicit rates of convergences are derived for all cases. Numerical experiments for ${L = κ^2 - Δ, κ>0}$ are performed to illustrate the results.