diff --git a/README.md b/README.md
index b915b1c9..dcafb478 100644
--- a/README.md
+++ b/README.md
@@ -71,33 +71,29 @@ This model can be formulated using JuMP code as:
 ```julia
 using JuMP, EAGO
 
-m = Model(EAGO.Optimizer)
+model = Model(EAGO.Optimizer)
 
 # Define bounded variables
 xL = [10.0; 0.0; 0.0; 0.0; 0.0; 85.0; 90.0; 3.0; 1.2; 145.0]
 xU = [2000.0; 16000.0; 120.0; 5000.0; 2000.0; 93.0; 95.0; 12.0; 4.0; 162.0]
-@variable(m, xL[i] <= x[i=1:10] <= xU[i])
+@variable(model, xL[i] <= x[i=1:10] <= xU[i])
 
-# Define nonlinear constraints
-@NLconstraint(m, e1, -x[1]*(1.12 + 0.13167*x[8] - 0.00667*(x[8])^2) + x[4] == 0.0)
-@NLconstraint(m, e3, -0.001*x[4]*x[9]*x[6]/(98.0 - x[6]) + x[3] == 0.0)
-@NLconstraint(m, e4, -(1.098*x[8] - 0.038*(x[8])^2) - 0.325*x[6] + x[7] == 57.425)
-@NLconstraint(m, e5, -(x[2] + x[5])/x[1] + x[8] == 0.0)
+# Define constraints
+@constraint(model, e1, -x[1]*(1.12 + 0.13167*x[8] - 0.00667*(x[8])^2) + x[4] == 0.0)
+@constraint(model, e2, -x[1] + 1.22*x[4] - x[5] == 0.0)
+@constraint(model, e3, -0.001*x[4]*x[9]*x[6]/(98.0 - x[6]) + x[3] == 0.0)
+@constraint(model, e4, -(1.098*x[8] - 0.038*(x[8])^2) - 0.325*x[6] + x[7] == 57.425)
+@constraint(model, e5, -(x[2] + x[5])/x[1] + x[8] == 0.0)
+@constraint(model, e6, x[9] + 0.222*x[10] == 35.82)
+@constraint(model, e7, -3.0*x[7] + x[10] == -133.0)
 
-# Define linear constraints
-@constraint(m, e2, -x[1] + 1.22*x[4] - x[5] == 0.0)
-@constraint(m, e6, x[9] + 0.222*x[10] == 35.82)
-@constraint(m, e7, -3.0*x[7] + x[10] == -133.0)
-
-# Define nonlinear objective
-@NLobjective(m, Max, 0.063*x[4]*x[7] - 5.04*x[1] - 0.035*x[2] - 10*x[3] - 3.36*x[5])
+# Define objective
+@objective(model, Max, 0.063*x[4]*x[7] - 5.04*x[1] - 0.035*x[2] - 10*x[3] - 3.36*x[5])
 
 # Solve the optimization problem
-JuMP.optimize!(m)
+JuMP.optimize!(model)
 ```
 
-Special handling has been included for linear/quadratic functions defined using the `@constraint` macro in JuMP and these can generally be expected to perform better than specifying quadratic or linear terms with the `@NLconstraint` macro.
-
 ## A Cautionary Note on Global Optimization
 
 As a global optimization platform, EAGO's solvers can be used to find solutions of general nonconvex problems with a guaranteed certificate of optimality. However, global solvers suffer from the curse of dimensionality and therefore their performance is outstripped by convex/local solvers. For users interested in large-scale applications, be warned that problems generally larger than a few variables may prove challenging for certain types of global optimization problems.
diff --git a/src/eago_optimizer/moi_wrapper.jl b/src/eago_optimizer/moi_wrapper.jl
index 1b6301fd..6727a4e0 100644
--- a/src/eago_optimizer/moi_wrapper.jl
+++ b/src/eago_optimizer/moi_wrapper.jl
@@ -21,6 +21,7 @@ function check_inbounds!(m::Optimizer, vi::VI)
     end
     return nothing
 end
+check_inbounds!(m::Optimizer, f::Number) = nothing
 check_inbounds!(m::Optimizer, f::SAF) = foreach(x -> check_inbounds!(m, x.variable), f.terms)
 check_inbounds!(m::Optimizer, f::VECOFVAR) = foreach(x -> check_inbounds!(m, x), f.variables)
 function check_inbounds!(m::Optimizer, f::SQF)
@@ -31,19 +32,23 @@ function check_inbounds!(m::Optimizer, f::SQF)
     end
     return nothing
 end
+function check_inbounds!(m::Optimizer, f::MOI.ScalarNonlinearFunction)
+    foreach(x -> check_inbounds!(m, x), f.args)
+    return nothing
+end
 
 MOI.supports_constraint(::Optimizer, ::Type{VI}, ::Type{S}) where S <: VAR_SETS = true
-MOI.supports_constraint(::Optimizer,::Type{T},::Type{S}) where {T<:Union{SAF,SQF},S<:INEQ_SETS} = true
+MOI.supports_constraint(::Optimizer,::Type{T},::Type{S}) where {T<:Union{SAF,SQF,MOI.ScalarNonlinearFunction},S<:INEQ_SETS} = true
 
 MOI.is_valid(m::Optimizer, v::VI) = (1 <= v.value <= m._input_problem._variable_count)
 MOI.is_valid(m::Optimizer, c::CI{VI,S}) where S <: VAR_SETS = (1 <= c.value <= m._input_problem._variable_count)
-MOI.is_valid(m::Optimizer, c::CI{F,S}) where {F<:Union{SAF,SQF}, S<:INEQ_SETS} = (1 <= c.value <= length(_constraints(m,F,S)))
+MOI.is_valid(m::Optimizer, c::CI{F,S}) where {F<:Union{SAF,SQF,MOI.ScalarNonlinearFunction}, S<:INEQ_SETS} = (1 <= c.value <= length(_constraints(m,F,S)))
 
 MOI.get(m::Optimizer, ::MOI.NumberOfConstraints{VI,S}) where S<:VAR_SETS = length(_constraints(m,VI,S))
-MOI.get(m::Optimizer, ::MOI.NumberOfConstraints{F,S}) where {F<:Union{SAF,SQF},S<:INEQ_SETS} = length(_constraints(m,F,S))
+MOI.get(m::Optimizer, ::MOI.NumberOfConstraints{F,S}) where {F<:Union{SAF,SQF,MOI.ScalarNonlinearFunction},S<:INEQ_SETS} = length(_constraints(m,F,S))
 
 MOI.get(m::Optimizer, ::MOI.ListOfConstraintIndices{VI,S}) where S<:VAR_SETS = collect(keys(_constraints(m,VI,S)))
-MOI.get(m::Optimizer, ::MOI.ListOfConstraintIndices{F,S}) where {F<:Union{SAF,SQF},S<:INEQ_SETS} = collect(keys(_constraints(m,F,S)))
+MOI.get(m::Optimizer, ::MOI.ListOfConstraintIndices{F,S}) where {F<:Union{SAF,SQF,MOI.ScalarNonlinearFunction},S<:INEQ_SETS} = collect(keys(_constraints(m,F,S)))
 
 function MOI.add_variable(m::Optimizer)
     vi = VI(m._input_problem._variable_count += 1)
@@ -57,10 +62,31 @@ function MOI.add_variable(m::Optimizer, name::String)
     return vi
 end
 
-function MOI.add_constraint(m::Optimizer, f::F, s::S) where {F<:Union{SAF,SQF},S<:INEQ_SETS}
+function MOI.add_constraint(m::Optimizer, f::F, s::S) where {F<:Union{SAF,SQF,MOI.ScalarNonlinearFunction},S<:INEQ_SETS}
     check_inbounds!(m, f)
     ci = CI{F, S}(m._input_problem._constraint_count += 1)
-    _constraints(m, F, S)[ci] = (f, s)
+    if f isa MOI.ScalarNonlinearFunction
+        if isnothing(m._input_problem._nlp_data)
+            model = MOI.Nonlinear.Model()
+            backend = MOI.Nonlinear.SparseReverseMode()
+            vars = MOI.get(m, MOI.ListOfVariableIndices())
+            evaluator = MOI.Nonlinear.Evaluator(model, backend, vars) 
+            m._input_problem._nlp_data = MOI.NLPBlockData(evaluator)
+        end
+        MOI.Nonlinear.add_constraint(m._input_problem._nlp_data.evaluator.model, f, s)
+        constraint_bounds = m._input_problem._nlp_data.constraint_bounds
+        has_objective = m._input_problem._nlp_data.has_objective
+        if s isa MOI.LessThan
+            push!(constraint_bounds, MOI.NLPBoundsPair(-Inf, s.upper))
+        elseif s isa MOI.GreaterThan
+            push!(constraint_bounds, MOI.NLPBoundsPair(s.lower, Inf))
+        else
+            push!(constraint_bounds, MOI.NLPBoundsPair(s.value, s.value))
+        end
+        m._input_problem._nlp_data = MOI.NLPBlockData(constraint_bounds, m._input_problem._nlp_data.evaluator, has_objective)
+    else
+        _constraints(m, F, S)[ci] = (f, s)
+    end
     return ci
 end
 function MOI.add_constraint(m::Optimizer, f::VI, s::S) where S<:VAR_SETS
@@ -230,7 +256,7 @@ MOI.get(m::Optimizer, p::MOI.RawStatusString) = string(m._global_optimizer._end_
 #####
 MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
 MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true
-MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{F}) where {F <: Union{VI, SAF, SQF}} = true
+MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{F}) where {F <: Union{VI, SAF, SQF, MOI.ScalarNonlinearFunction}} = true
 
 function MOI.set(m::Optimizer, ::MOI.NLPBlock, nlp_data::MOI.NLPBlockData)
     if nlp_data.has_objective
@@ -239,9 +265,20 @@ function MOI.set(m::Optimizer, ::MOI.NLPBlock, nlp_data::MOI.NLPBlockData)
     m._input_problem._nlp_data = nlp_data
 end
 
-function MOI.set(m::Optimizer, ::MOI.ObjectiveFunction{T}, f::T) where T <: Union{VI,SAF,SQF}
+function MOI.set(m::Optimizer, ::MOI.ObjectiveFunction{T}, f::T) where T <: Union{VI,SAF,SQF,MOI.ScalarNonlinearFunction}
     check_inbounds!(m, f)
-    m._input_problem._objective = f
+    if f isa MOI.ScalarNonlinearFunction
+        if isnothing(m._input_problem._nlp_data)
+            model = MOI.Nonlinear.Model()
+            backend = MOI.Nonlinear.SparseReverseMode()
+            vars = MOI.get(m, MOI.ListOfVariableIndices())
+            evaluator = MOI.Nonlinear.Evaluator(model, backend, vars) 
+            m._input_problem._nlp_data = MOI.NLPBlockData(evaluator)
+        end
+        MOI.Nonlinear.set_objective(m._input_problem._nlp_data.evaluator.model, f)
+    else
+        m._input_problem._objective = f
+    end
 end
 MOI.get(m::Optimizer, ::MOI.ObjectiveFunction{T}) where T <: Union{VI,SAF,SQF} = m._input_problem._objective
 MOI.get(m::Optimizer, ::MOI.ObjectiveFunctionType) = typeof(m._input_problem._objective)
diff --git a/test/moit_tests.jl b/test/moit_tests.jl
index 7bf9f85a..5e888c80 100644
--- a/test/moit_tests.jl
+++ b/test/moit_tests.jl
@@ -73,9 +73,6 @@ function test_runtests()
             "test_conic_NormOneCone_VectorAffineFunction",
             "test_conic_NormOneCone_VectorOfVariables",
             "test_conic_linear_VectorAffineFunction",
-            "test_cpsat_Circuit",
-            "test_cpsat_CountAtLeast",
-            "test_cpsat_Table",
             "test_linear_Semicontinuous_integration",
             "test_linear_Semiinteger_integration",
             "test_linear_integer_solve_twice",
@@ -83,15 +80,39 @@ function test_runtests()
             "test_quadratic_duplicate_terms",
             "test_quadratic_integration",
             "test_quadratic_nonhomogeneous",
-            "test_quadratic_constraint_LessThan",
-            "test_quadratic_constraint_GreaterThan",
             "test_modification_affine_deletion_edge_cases",
-            "test_solve_SOS2_add_and_delete",
+            # MOI bridge issues
+            "test_basic_VectorNonlinearFunction_AllDifferent",
+            "test_basic_VectorNonlinearFunction_BinPacking",
+            "test_basic_VectorNonlinearFunction_Circuit",
+            "test_basic_VectorNonlinearFunction_Complements",
+            "test_basic_VectorNonlinearFunction_CountAtLeast",
+            "test_basic_VectorNonlinearFunction_CountBelongs",
+            "test_basic_VectorNonlinearFunction_CountDistinct",
+            "test_basic_VectorNonlinearFunction_CountGreaterThan",
+            "test_basic_VectorNonlinearFunction_GeometricMeanCone",
+            "test_basic_VectorNonlinearFunction_HyperRectangle",
+            "test_basic_VectorNonlinearFunction_Nonnegatives",
+            "test_basic_VectorNonlinearFunction_Nonpositives",
+            "test_basic_VectorNonlinearFunction_NormInfinityCone",
+            "test_basic_VectorNonlinearFunction_NormOneCone",
+            "test_basic_VectorNonlinearFunction_RotatedSecondOrderCone",
+            "test_basic_VectorNonlinearFunction_SecondOrderCone",
+            "test_basic_VectorNonlinearFunction_SOS1",
+            "test_basic_VectorNonlinearFunction_SOS2",
+            "test_basic_VectorNonlinearFunction_Table",
+            "test_basic_VectorNonlinearFunction_Zeros",
+            # TODO: bug related to unbounded_check
+            "test_nonlinear_expression_quartic",
+            # TODO: bug related to reform_epigraph_min
+            "test_nonlinear_expression_overrides_objective",
+            # Remove after MOI v1.31.3
+            "test_nonlinear_expression_hs110",
             # TODO: wrong error thrown. Likely (trivial) bug in MOI wrapper
             "test_model_LowerBoundAlreadySet",
             "test_model_UpperBoundAlreadySet",
         ],
-        exclude_tests_after = v"1.22.0",
+        exclude_tests_after = v"1.31.2",
     )
     return
 end