I know that the divergence theorem admits a direct generalization to $\mathbb{R}^n$ (and thus $\mathbb{R}^4$). My question is :
Is there also a version of Kelvin-Stokes (curl) theorem for surfaces in $\mathbb{R}^4$ ? I'm guessing that Stokes theorem can be also applied in this case, but how to define the curl in $\mathbb{R}^4$ ?
You can get an understanding of the required structure of the Stokes generalization by considering the $\mathbb{R}^{3}$ coordinate expansion of Stokes theorem. Given a two parameter surface
$$\mathbf{x} = \mathbf{x}(a, b) = \mathbf{e}_k x_k(a, b),$$
we can expand the Stokes integral on that surface $$\begin{aligned}\int d\mathbf{A} \cdot \left( { \boldsymbol{\nabla} \times \mathbf{F} } \right)&=\int da db \left( { \frac{\partial {\mathbf{x}}}{\partial {a}} \times \frac{\partial {\mathbf{x}}}{\partial {b}} } \right)\cdot \left( { \boldsymbol{\nabla} \times \mathbf{F} } \right) \\ &=\int da db \epsilon_{rst}\frac{\partial {x_r}}{\partial {a}} \frac{\partial {x_s}}{\partial {b}} \epsilon_{ruv}\partial_u F_v \\ &=\int da db\delta_{st}^{[uv]}\frac{\partial {x_r}}{\partial {a}} \frac{\partial {x_s}}{\partial {b}} \partial_u F_v \\ &=\int da db\frac{\partial(x_u, x_v)}{\partial(a,b)}\partial_u F_v.\end{aligned}$$
You can see that the primary structure involved here is the antisymmetric volume element. With a bit of care to adjust upper and lower indexes, this can be used as a Stokes integral for a 4D space. Given a two parameter surface in a 4D space $ \mathbf{x}(a,b) = \gamma_\mu x^\mu(a, b) $, and a four vector $ \mathbf{F} = F^\mu \gamma_\mu $, the 4D equivalent to the standard $\mathbb{R}^{3}$ Stokes integral has the coordinate expansion
$$\begin{aligned}\int da db\frac{\partial(x^\mu, x^\nu)}{\partial(a,b)}\partial_\mu F_\nu&=\int da db\left( { \frac{\partial {x^\mu}}{\partial {a}}\frac{\partial {x^\nu}}{\partial {b}} -\frac{\partial {x^\mu}}{\partial {b}}\frac{\partial {x^\nu}}{\partial {a}} } \right)\partial_\mu F_\nu \\ &=\int db \frac{\partial {x^\nu}}{\partial {b}} \int da \frac{\partial {x^\mu}}{\partial {a}} \frac{\partial {F_\nu}}{\partial {x^\mu}}-\int da \frac{\partial {x^\nu}}{\partial {a}} \int db\frac{\partial {x^\mu}}{\partial {b}}\frac{\partial {F_\nu}}{\partial {x^\mu}} \\ &=\int db {\left.{{\left( {\frac{\partial {x^\nu}}{\partial {b}} F_\nu} \right)}}\right\vert}_{{\Delta a}}-\int da {\left.{{\left( {\frac{\partial {x^\nu}}{\partial {a}} F_\nu} \right)}}\right\vert}_{{\Delta b}} \\ &=\int db {\left.{{\left( {\frac{\partial {\mathbf{x}}}{\partial {b}} \cdot \mathbf{F}} \right)}}\right\vert}_{{\Delta a}}-\int da {\left.{{\left( {\frac{\partial {\mathbf{x}}}{\partial {a}} \cdot \mathbf{F}} \right)}}\right\vert}_{{\Delta b}}.\end{aligned}$$
This is an integral around the boundary of the surface parameterized by $ (a,b) $. A more careful approach would have to consider breaking up the surface into smaller regions and summing the contributions from each region, but this is the rough idea.
As mentioned in the comments, there are a number of integral formalisms that describe the generalized Stokes theorem. One can generalize Stokes theorem by generalizing the degree of the volume element (two, three, or four parameter volume elements), as well as operating on algebraic structures of different degrees (four-vectors, second-degree antisymmetric tensors, or third degree antisymmetric tensors). This can be described using tensor algebra, with differential forms, or using geometric calculus.
The geometric calculus generalization of Stokes theorem is probably closer to the standard 3D vector formalism than any of the other formalisms. In the GC formalism, Stokes theorem has the form
$$\int d^k x \cdot \left( { \partial \wedge F } \right) = \int d^{k-1} x \cdot F,$$
where $ \partial $ is the projection of the gradient onto the integration manifold (i.e. the subspace that is being integrated over), $ F $ is a ``blade'' of grade $ <= k -1 $, $\wedge$ is the wedge product, and $d^k x = dx_a \wedge dx_b \wedge \cdots $ is the k-volume element (a wedge product of the differentials along each of the parameterization directions). The wedge product is antisymmetric like the cross product, and in 3D the wedge product is dual to the cross product ($\mathbf{a} \wedge \mathbf{b} = \mathbf{e}_1 \mathbf{e}_2 \mathbf{e}_3 (\mathbf{a} \times \mathbf{b})$).
To give precise meaning to all of this, I'd recommend Alan Macdonald's excellent little book Vector and Geometric Calculus (a prerequisite for that book is some basic knowledge of geometric algebra).