Let $f:2^{\{1,2\}} \rightarrow \mathbb{R}$ with
$$f(\{1\})=100$$
$$f(\{2\})=\frac{1}{2}$$
$$f(\{1,2\})= 101$$
$$f(\emptyset)=0$$
$f$ is supermodular. The solution of
\begin{equation}
\max_{S \subseteq V}\, f(S)\quad \text{s.t}\, \lvert S \rvert \leq k
\end{equation}
for $k=1$ is $\{1\}$ with $f(\{1\})=100$.
If we now minimize $-f$ and turn around the constraint direction, we obtain the following problem:
\begin{equation}
\min_{S \subseteq V}\, -f(S)\quad \text{s.t}\, \lvert S \rvert \geq 1
\end{equation}
Now, we obtain the optimum is at $S=\{1,2\}$. Thus, the two problems are not equivalent and just switching the constraint direction is not enough.
In general, I don't think it's very easy, but there are a few things you can do. Firstly $\|x\|_2 \leq \epsilon$ is equivalent to $\|x\|_2^2 \leq \epsilon^2$, that is
$$
\sum_{i=1}^{n}x_i^2 \leq \epsilon^2.\tag{1}\label{1}
$$
You can then introduce $n$ auxiliary scalar decision variables, $t_1,\ldots, t_n$ and replace \eqref{1} by
\begin{align}
|x_i| {}\leq{}& t_i,
\\
\sum_{i=1}^{n}t_i {}\leq{}& \epsilon^2.
\end{align}
In case you cannot handle $\sum_{i=1}^{n}t_i {}\leq{} \epsilon^2$ directly, you can associate with it the penalty function $\psi(t) = \max\{\sum_{i=1}^{n}t_i -\epsilon^2, 0\}$ and define the modified cost function
$$
\tilde{f}(x,t; \lambda) = f(x) +\lambda \psi(t).
$$
This is a Lipschitz continuous and you can apply a penalty method on top of your algorithm (or an augmented Lagrangian method - you have to see which one works better). Simply take $\lambda \to \infty$ and stop once $\psi(t)$ becomes sufficiently small.
Let me just clarify that for each $\lambda$ you will have to solve the following problem
\begin{align}
&\operatorname*{Minimize}_{x,t} \tilde{f}(x,t; \lambda)
\\
&\text{subject to } |x_i| \leq t_i, \text{ for all } i =1,\ldots, n
\end{align}
Update 1. Under the stronger assumption that $f$ is continuously differentiable with Lipschitz gradient, it turns out the the original problem is equivalent to the unconstrained minimization problem
$$
\operatorname*{Minimize}_{x} f(x) - \frac{\gamma}{2}\|\nabla f(x)\|^2 + \frac{\gamma}{2}\operatorname{dist}^2_D(x - \gamma \nabla f(x))^2,\tag{2}\label{2}
$$
where $D$ is the set of box constraints (I'm using the same notation as the paper you cited), and $\operatorname{dist}^2_D$ is the squared distance from $D$, which is easy to compute. The cost function in \eqref{2} is known as the forward-backward envelope of the original problem (see for example this paper).
If $f$ is just Lipschitz the above does not apply. Generally speaking, this is too weak an assumption in optimization (with the exception of global optimization where people resort to branch-and-bound-type methods).
Update 2. You can use a homeomorphism between $B_2 := \{x\in\mathbb{R}^n: \|x\|_2 \leq 1\}$ and $B_\infty := \{x\in\mathbb{R}^n: \|x\|_\infty \leq 1\}$, i.e., a function $\phi: B_2\to B_{\infty}$ such as
$$
\phi: B_{\infty} \ni x \mapsto \phi(x) = \frac{\|x\|_\infty}{\|x\|_2}x \in B_2
$$
for $x\neq 0$ and $\phi(0) = 0$
The inverse is
$$
\phi^{-1}: B_2\ni z \mapsto \phi^{-1}(z) = \frac{\|z\|_2}{\|z\|_\infty}z \in B_{\infty}
$$
for $z\neq 0$ and $\phi^{-1}(0) = 0$.
The constraints in your problem can be written as $\frac{x}{\epsilon} \in B_2$. Define $u = \phi^{-1}(x/\epsilon)$. Then $\frac{x}{\epsilon} = \phi(u) \Leftrightarrow x = \epsilon \phi(u)$. The constraints become
\begin{align}
\frac{x}{\epsilon} \in B_2
\Leftrightarrow{}&{}\frac{x}{\epsilon} \in \phi (B_\infty)
\\
\Leftrightarrow{}&{}\phi(u) \in \phi (B_\infty)
\\
\Leftrightarrow{}&{}u \in B_{\infty}
\end{align}
and the optimization problem becomes
\begin{align}
&\operatorname*{Minimize}_u f(\epsilon \phi(u))
\\
&\text{subject to } -1 \leq u \leq 1
\end{align}
Once you find an optimal solution $u^\star$ (if one exists), you can retrieve $x^\star = \epsilon \phi(u^\star)$.
Best Answer
This only works if $A$ is positive definite. In this case, just drop the quadratic part and obtain $$\begin{array}{rcl} \min\limits_{x,y} && y + a^Tx\\ st && Bx\leq b\\ && x^TAx \leq y \end{array}$$ Then, we play a long circuitous game with this last constraint where $$\begin{array}{rl} &y\geq x^TAx\\ \Longrightarrow& 0\geq x^TAx - y\\ \Longrightarrow& 0\geq 4x^TAx - 4y\\ \Longrightarrow& 0\geq 4x^TAx + (1-y)^2-(1+y)^2\\ \Longrightarrow& (1+y)^2\geq 4x^TAx + (1-y)^2\\ \Longrightarrow& (1+y)^2\geq 4x^TAx + (1-y)^2, 1+y\geq 0\\ \end{array}$$ where the extra constraint, $1+y\geq 0$, can be added since $y\geq 0$ since $x^TAx\geq 0$ since $A\succ 0$ and $y\geq x^TAx$. Anyway, then we have $$\begin{array}{rl} \Longrightarrow& 1+y\geq \sqrt{4x^TAx + (1-y)^2}, 1+y\geq 0\\ \Longrightarrow& 1+y\geq \sqrt{4x^TU^TUx + (1-y)^2}, 1+y\geq 0 \end{array}$$ where we use the Choleski factorization $A=U^TU$ since $A\succ 0$. Then,
$$\begin{array}{rl} \Longrightarrow& 1+y\geq \left\|\begin{bmatrix}2Ux\\1-y\end{bmatrix}\right\|, 1+y\geq 0\\ \Longrightarrow& \begin{bmatrix}1+y\\2Ux\\1-y\end{bmatrix}\succeq_Q 0\\ \end{array}$$ Putting this all together, we have $$\begin{array}{rcl} \min\limits_{x,y} && y + a^Tx\\ st && Bx\leq b\\ && \begin{bmatrix}1+y\\2Ux\\1-y\end{bmatrix}\succeq_Q 0 \end{array}$$