How can neural networks approximate any continuous function but have $\mathcal{VC}$ dimension only proportional to their number of parameters? Neural networks typically have $\mathcal{VC}$ dimension that is proportional to their number of parameters and inputs. For example, see the papers Vapnik-Chervonenkis dimension of recurrent neural networks (1998) by Pascal Koirana and Eduardo D. Sontag and VC Dimension of Neural Networks (1998) by Eduardo D. Sontag for more details.

On the other hand, the universal approximation theorem (UAT) tells us that neural networks can approximate any continuous function. See Approximation by Superpositions of a Sigmoidal Function (1989) by G. Cybenko for more details. 

Although I realize that the typical UAT only applies to continuous functions, the UAT and the results about the $\mathcal{VC}$ dimension of neural networks seem to be a little bit contradictory, but this is only if you don't know the definition of $\mathcal{VC}$ dimension and the implications of the UAT.

So, how come that neural networks approximate any continuous function, but, at the same time, they usually have a $\mathcal{VC}$ dimension that is only proportional to their number of parameters? What is the relationship between the two?
I think DuttaA is right in his comment. VC Dimension of Neural Networks establishes VC bounds depending on the number of weights, whereas the UAT refers to a class of neural networks in which the number of weights a particular network can have is not bounded, although it needs to be finite.
I think that we can show, from theorem 2 and the observations below theorem 3 in Approximation by Superpositions of a Sigmoidal Function, that the VC dimension of
$$S=\left\{\sum_{i=1}^N \alpha_i\sigma(y_i^T x + \theta_i) : N\in\mathbb N, \alpha_i, \theta_i \in\mathbb R, y_i\in\mathbb{R}^n \right\}$$
is infinite.
Let $\{(x_i, y_i)\}_{i=1}^k$ be a sample of arbitrary size $k\in\mathbb N$, and let us see that there is a function in $S$ which can correctly classify it, i.e., $S$ shatters $\{x_i\}_{i=1}^k$.
We note $B(x, \varepsilon) := \{ y\in\mathbb{R}^n : d(x,y) &lt; \varepsilon \}$ (this is just standard notation to denote a ball).
First, let $\varepsilon &gt; 0$ be such that $B(x_i, \varepsilon)\cap B(x_j, \varepsilon) = \emptyset$ every time that $i \ne j$.
Now define $D = \cup_{y_i=1} B(x_i, \varepsilon)$. Define $f_{\varepsilon}(x)$ as in the observations below theorem 3 of Cybenko's paper, and use theorem 2 to find a function $G(x)$ in $S$ that classifies correctly all points at least $\varepsilon$ away from the boundary of $D$, i.e., all points in the sample.