Here is a proof that adding the condition $\nabla I(\mathbf x) \ne 0$ if $I(\mathbf x) = c$ makes the statement true. In fact, a stronger statement will be proved:
Let $M \subseteq \mathbb R^n$ be a nonempty, $k$-dimensional differentiable manifold, and suppose that $F(\mathbf x)$ is tangent to $M$ at every point $\mathbf x \in M$. Then, $M$ is invariant.
The lemma stated in the question follows from the latter due to the level set $I^{-1}(c)$ being an $(n - 1)$-dimensional differentiable manifold (this follows from the implicit function theorem when $I$ has nonzero gradient).
Now to the proof. Fix $\mathbf x_0 \in M$ and let $\varphi: U \longrightarrow V$, where $U$ is a neighborhood of $0$ in $\mathbb R^n$ and $V$ is a neighborhood of $\mathbf x_0$, be a diffeomorphism such that $\varphi(0) = \mathbf x_0$ and $\varphi(E \cap U) = M \cap V$, where $E = \mathrm{Span}\{e_1, \, \dots, \, e_k\}$. Any trajectory $t \mapsto \mathbf x(t)$ whose support lies entirely inside $V$ can be mapped to $U$ via $\varphi^{-1}$, i.e. we can define $\mathbf y(t) = \varphi^{-1}(\mathbf x(t))$. Then, $\mathbf y$ satisfies $$\dot{\mathbf y} = J_{\varphi^{-1}}(\mathbf x)\dot{\mathbf x} = (J_{\varphi}(\mathbf y))^{-1}F(\varphi(\mathbf y)) =: G(\mathbf y).$$ Now, whenever $\mathbf y \in E$ (which means $\varphi(\mathbf y) \in M$), $F(\varphi(\mathbf y))$ belongs to the tangent space to $M$ at $\mathbf x
= \varphi(\textbf y)$, $T_{\mathbf x}M = J_{\varphi}(\mathbf y)(E)$. Thus $F(\mathbf x) = J_{\varphi}(\mathbf y)\mathbf v$ for some $\mathbf v \in E$, and $G(\mathbf y) = (J_{\varphi}(\mathbf y))^{-1}J_{\varphi}(\mathbf y)\mathbf v = \mathbf v$. If we manage to show that the trajctories of $\dot{\mathbf y} = G(\mathbf y)$ starting at any $\tilde{\mathbf y} \in E$ lie (locally) in $E$ (that is, $E$ is invariant for $G$) we are done, since these trajectories are one-to-one with the trajectories of $\dot{\mathbf x} = F(\mathbf x)$ starting at $\tilde{\mathbf x} \in M$ via $\varphi$, and $\varphi$ maps points of $E$ to points of $M$.
To do this, consider any solution of $\dot{\mathbf y} = G(\mathbf y)$ with initial value $\mathbf y(0) = \tilde{\mathbf y}$. If we let $\tilde{\mathbf y} = (\tilde{\mathbf z}, \, 0)$ (where the $0$ is $(n - k)$-dimensional) and, for $\mathbf z \in \mathbb R^k$, $G((\mathbf z, \, 0)) = (H(\mathbf z), \, 0)$ (recall that $G(E \cap U) \subseteq E$), we can consider the system $\begin{cases} \dot{\mathbf z} = H(\mathbf z) \\ \mathbf z(0) = \tilde{\mathbf z} \end{cases}$ in $\mathbb R^k$. Due to the regularity of $H$, local existence and uniqueness of the solution hold. But now, if $\mathbf z(t)$ is the local solution, $\mathbf y(t) = (\mathbf z(t), \, 0)$ is a solution to $\dot{\mathbf y} = F(\mathbf y)$ that lies entirely in $E$, and it is in fact the only one.
Best Answer
Hint: Check how vector field points along the boundary of your domain of interest. If it points inside, then Bony-Brezis theorem can be applied. Or you just can say that vector field along the boundary doesn't let trajectories go out of domain $x \leqslant 0$.