diff --git a/src/vertical.jl b/src/vertical.jl
index 2ac1e69..860bbdc 100644
--- a/src/vertical.jl
+++ b/src/vertical.jl
@@ -2,6 +2,7 @@
 const VF = Union{
     MOI.VectorAffineFunction,
     MOI.VectorQuadraticFunction,
+    MOI.VectorNonlinearFunction,
 }
 
 function has_complementarity(model::MOI.ModelLike)
@@ -29,28 +30,13 @@ function _is_single_variable(func::MOI.ScalarQuadraticFunction)
         false
     end
 end
+function _is_single_variable(func::MOI.ScalarNonlinearFunction)
+    return func.head == :+ && length(func.args) == 1 && isa(func.args[1], MOI.VariableIndex)
+end
 _get_variable(func::MOI.ScalarAffineFunction) = func.terms[1].variable
 _get_variable(func::MOI.ScalarQuadraticFunction) = func.affine_terms[1].variable
+_get_variable(func::MOI.ScalarNonlinearFunction) = func.args[1]
 
-# Reformulate
-function _vertical_formulation!(model, terms::Vector{MOI.ScalarAffineTerm{T}}, c::T) where T
-    x = MOI.add_variable(model)
-    push!(terms, MOI.ScalarAffineTerm{T}(-one(T), x))
-    func = MOI.ScalarAffineFunction{T}(terms, zero(T))
-    MOI.add_constraint(model, func, MOI.EqualTo{T}(-c))
-    return x
-end
-
-function _add_slack!(model, func::MOI.ScalarAffineFunction{T}) where T
-    x = MOI.add_variable(model)
-    push!(func.terms, MOI.ScalarAffineTerm{T}(-one(T), x))
-    return x
-end
-function _add_slack!(model, func::MOI.ScalarQuadraticFunction{T}) where T
-    x = MOI.add_variable(model)
-    push!(func.affine_terms, MOI.ScalarAffineTerm{T}(-one(T), x))
-    return x
-end
 
 # TODO: add support for ScalarNonlinearTerm
 function _parse_complementarity_constraint(fun::MOI.AbstractVectorFunction, n_comp)
@@ -129,8 +115,9 @@ function reformulate_to_vertical!(model::MOI.ModelLike)
                             push!(ind_cc2, x2)
                         else
                             # Else, reformulate LHS using vertical form
-                            x1 = _add_slack!(model, lhs)
-                            MOI.add_constraint(model, lhs, MOI.EqualTo{Float64}(0))
+                            x1 = MOI.add_variable(model)
+                            new_lhs = MOIU.operate!(-, Float64, lhs, x1)
+                            MOI.add_constraint(model, new_lhs, MOI.EqualTo{Float64}(0))
                             push!(ind_cc1, x1)
                             push!(ind_cc2, x2)
                         end
diff --git a/test/runtests.jl b/test/runtests.jl
index 8efd349..0722628 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -21,13 +21,34 @@ function fletcher_leyffer_ex1_nonlinear_model()
     return model
 end
 
+function nonlinear_test_model()
+    model = Model()
+    @variable(model, x >= 0.0)
+    @variable(model, y >= 0.0)
+    @objective(model, Min, x^2 + y^2 - 4*x*y)
+    # Build complementarity constraints with nonlinear expression
+    @constraint(model, [sin(x), y] ∈ MOI.Complements(2))
+    return model
+end
+
+function nonlinear_test_reformulated_model()
+    model = Model()
+    @variable(model, x >= 0.0)
+    @variable(model, y >= 0.0)
+    @variable(model, slack)
+    @objective(model, Min, x^2 + y^2 - 4*x*y)
+    # Build complementarity constraints with nonlinear expression
+    @constraint(model, sin(x) == slack)
+    @constraint(model, slack * y <= 0.0)
+    return model
+end
+
 expected_models = Dict(
     Instances.fletcher_leyffer_ex1_model => fletcher_leyffer_ex1_nonlinear_model,
 )
 
 function test_model(model_func)
     model = model_func()
-    model = Instances.fletcher_leyffer_ex1_model()
     inner = MOI.Utilities.Model{Float64}()
     set_optimizer(model, () -> ComplementOpt.Optimizer(inner))
     MOI.Utilities.attach_optimizer(model)
@@ -42,6 +63,20 @@ function test_model(model_func)
     end
 end
 
+function test_nonlinear_expr()
+    model = nonlinear_test_model()
+    inner = MOI.Utilities.Model{Float64}()
+    set_optimizer(model, () -> ComplementOpt.Optimizer(inner))
+    MOI.Utilities.attach_optimizer(model)
+    MOI.Utilities.attach_optimizer(backend(model).optimizer.model)
+
+    expected = nonlinear_test_reformulated_model()
+    MOI.Bridges._test_structural_identical(
+        unsafe_backend(model).optimizer,
+        backend(expected),
+    )
+end
+
 instances = filter(names(Instances; all = true)) do name
     # The function types start with `#`
     s = String(name)
@@ -51,3 +86,9 @@ end
 @testset "$name" for name in instances
     test_model(getfield(Instances, name))
 end
+
+@testset "Unit-tests" begin
+    test_nonlinear_expr()
+end
+
+