We study the problem of how efficiently, in terms of the number of parameters, deep neural networks with the ReLU activation function can approximate functions in the Sobolev space $W^s(L_q(\Omega))$ on a bounded domain $\Omega$, where the error is measured in $L_p(\Omega)$. This problem is important for studying the application of neural networks in scientific computing and has previously been solved only in the case $p=q=\infty$. Our contribution is to provide a solution for all $1\leq p,q\leq \infty$ and $s > 0$. Our results show that deep ReLU networks significantly outperform classical methods of approximation, but that this comes at the cost of parameters which are not encodable.