We establish summability results for coefficient sequences of Wiener-Hermite polynomial chaos expansions for countably-parametric solutions of linear elliptic and parabolic partial differential equations with Gaussian random field inputs. Our proof technique is based on the holomorphic extension of solution maps, and thus allows a unified ``differentiation-free'' sparsity analysis applicable also to posterior densities in Bayesian inverse problems. We present dimension-independent convergence rates of high-dimensional deterministic numerical approximation schemes such as single-level and multi-level versions of anisotropic sparse-grid Hermite-Smolyak interpolation and quadrature in both forward and inverse uncertainty quantification. Additionally, we provide expression rate bounds for neural network approximations of the response surfaces to the parametric PDEs.