diff --git a/examples/circle_packing.jl b/examples/circle_packing.jl
new file mode 100644
index 0000000..3c06275
--- /dev/null
+++ b/examples/circle_packing.jl
@@ -0,0 +1,438 @@
+# Circle Packing Optimization Example - Packing Circles into Circular Containers
+
+using JuMP
+using SCIP
+using MadNLP
+using Ipopt
+using LinearAlgebra
+using Random
+using JSON3
+using NablaMatheuristicKit
+
+const NMK = NablaMatheuristicKit
+
+"""
+    CirclePackingProblem
+
+A structure to hold the data for a circle packing problem.
+
+# Fields
+- `num_containers::Int`: Number of circular containers
+- `container_radii::Vector{Float64}`: Radius of each container
+- `num_items::Int`: Number of items (circles)
+- `item_radii::Vector{Float64}`: Radius of each item
+- `item_rewards::Vector{Float64}`: Reward value of each item
+"""
+struct CirclePackingProblem
+    num_containers::Int
+    container_radii::Vector{Float64}
+    num_items::Int
+    item_radii::Vector{Float64}
+    item_rewards::Vector{Float64}
+
+    function CirclePackingProblem(
+        num_containers::Int,
+        container_radii::Vector{Float64},
+        num_items::Int,
+        item_radii::Vector{Float64},
+        item_rewards::Vector{Float64}
+    )
+        # Validate inputs
+        @assert num_containers == length(container_radii) "Number of containers must match container radii array"
+        @assert num_items == length(item_radii) == length(item_rewards) "Number of items must match radii and rewards arrays"
+        @assert all(container_radii .> 0) "Container radii must be positive"
+        @assert all(item_radii .> 0) "Item radii must be positive"
+
+        new(num_containers, container_radii, num_items, item_radii, item_rewards)
+    end
+end
+
+"""
+    create_sample_problem()
+
+Create a sample circle packing problem with 3 circular containers and 15 items.
+
+# Returns
+- `CirclePackingProblem`: A sample problem instance
+"""
+function create_sample_problem()
+    # Define 3 circular containers with different radii
+    num_containers = 2
+    container_radii = [20.0, 25.0]
+    #container_radii = [20.0, 15.0]
+
+    # Define 15 items with different radii and rewards
+    num_items = 15
+    #num_items = 7
+    Random.seed!(42)  # For reproducibility
+
+    # Generate random radii between 5 and 15
+    item_radii = round.(5.0 .+ 10.0 .* rand(num_items))
+
+    @show container_radii
+    @show item_radii
+
+    # Generate rewards somewhat proportional to the area of the circle
+    # with some randomness
+    item_rewards = [π * r^2 * (0.8 + 0.4 * rand()) for r in item_radii]
+
+    return CirclePackingProblem(
+        num_containers,
+        container_radii,
+        num_items,
+        item_radii,
+        item_rewards
+    )
+end
+
+"""
+    build_jump_model(problem::CirclePackingProblem)
+
+Build a JuMP model for the circle packing problem with continuous variables.
+
+# Arguments
+- `problem::CirclePackingProblem`: The problem instance
+
+# Returns
+- `JuMP.Model`: The JuMP model
+"""
+function build_jump_model(problem::CirclePackingProblem)
+    model = Model(Ipopt.Optimizer)
+
+    # Extract problem data for convenience
+    num_containers = problem.num_containers
+    container_radii = problem.container_radii
+    num_items = problem.num_items
+    item_radii = problem.item_radii
+    item_rewards = problem.item_rewards
+
+    # Decision variables
+    items = 1:num_items
+    containers = 1:num_containers
+
+    # x[i,c], y[i,c]: Position of item i in container c (if placed)
+    # These are relative to the center of the container
+    @variable(model, x[items, containers])
+    @variable(model, y[items, containers])
+
+    # z[i,c]: 1 if item i is placed in container c, 0 otherwise
+    @variable(model, 0 <= α[items, containers] <= 1, Bin)
+
+    # Objective: Maximize total reward of placed items
+    @objective(model, Max, sum(item_rewards[i] * α[i, c] for i in items, c in containers))
+
+    # Constraints
+    @constraint(model, covering[i in items], sum(α[i, c] for c in containers) <= 1)
+
+    @constraint(model, item_fits_container[i in items, c in containers],
+        x[i, c]^2 + y[i, c]^2 <= α[i, c] * (container_radii[c] - item_radii[c])^2
+    )
+
+    item_pairs = [(i, j) for i in items for j in items if i < j]
+
+    @constraint(model, non_overlap[(i, j) in item_pairs, c in containers],
+        (x[i, c] - x[j, c])^2 + (y[i, c] - y[j, c])^2 >= α[i, c] * α[j, c] * (item_radii[i] + item_radii[j])^2
+    )
+
+    @constraint(model, pos_x_lb[i in items, c in containers],
+        x[i, c] >= -α[i, c] * (container_radii[c] - item_radii[i])
+    )
+
+    @constraint(model, pos_x_ub[i in items, c in containers],
+        x[i, c] <= α[i, c] * (container_radii[c] - item_radii[i])
+    )
+
+    @constraint(model, pos_y_lb[i in items, c in containers],
+        y[i, c] >= -α[i, c] * (container_radii[c] - item_radii[i])
+    )
+
+    @constraint(model, pos_y_ub[i in items, c in containers],
+        y[i, c] <= α[i, c] * (container_radii[c] - item_radii[i])
+    )
+
+    return model
+end
+
+
+################################################
+
+mutable struct CirclePackingSearchSpace <: NMK.TreeSearch.AbstractSearchSpace
+    node_limit::Int
+    nb_nodes_evaluated::Int
+    problem::CirclePackingProblem
+    model::JuMP.Model
+    helper::NMK.MathOptState.DomainChangeTrackerHelper
+    root_state
+    tracker
+    var_names::Dict{MOI.VariableIndex, String}
+    # For DOT file generation
+    dot_nodes::Dict{Int, String}  # node_id => node label
+    dot_edges::Vector{Tuple{Int, Int, String}}  # (parent_id, child_id, edge label)
+    next_node_id::Int
+end
+
+function CirclePackingSearchSpace(
+    problem::CirclePackingProblem, 
+    model::JuMP.Model, 
+    helper::NMK.MathOptState.DomainChangeTrackerHelper,
+    root_state,
+    tracker; 
+    node_limit=1000
+)
+    var_names = Dict{MOI.VariableIndex, String}()
+    for var in JuMP.all_variables(model)
+        var_names[JuMP.index(var)] = JuMP.name(var)
+    end
+    
+    # Initialize DOT file generation structures
+    dot_nodes = Dict{Int, String}()
+    dot_edges = Vector{Tuple{Int, Int, String}}()
+    
+    return CirclePackingSearchSpace(
+        node_limit, 0, problem, model, helper, root_state, tracker, 
+        var_names, dot_nodes, dot_edges, 1
+    )
+end
+
+struct CirclePackingNode
+    id::Int
+    depth::Int
+    value_guess::Float64
+    state
+end
+
+struct CirclePackingResult
+    node_id::Int
+    assignement::Matrix{Float64} # (container_id, item_id)
+    obj_value::Float64
+    integral::Bool
+end
+
+function CirclePackingResult(node_id::Int, model::JuMP.Model, problem::CirclePackingProblem)
+    assignment = zeros(problem.num_items, problem.num_containers)
+    α = model[:α]
+    for container_id in 1:problem.num_containers, item_id in 1:problem.num_items
+        assignment[item_id, container_id] = value(α[item_id, container_id])
+    end
+
+    obj_value = MOI.get(model, MOI.ObjectiveValue())
+    integral = all(val -> abs(val - round(val)) < 1e-5, assignment)
+
+    return CirclePackingResult(
+        node_id,
+        assignment,
+        obj_value,
+        integral
+    )
+end
+
+##########
+##########
+#########
+
+# 1 is the best value (meaning the value x = round(x) +/- 1/2)
+# 0 is the worst (meaning x = round(x))
+function most_fractional_value(space, var_id)
+    val = MOI.get(JuMP.backend(space.model), MOI.VariablePrimal(), var_id)
+    return 2 * abs(round(val) - val)
+end
+
+NMK.TreeSearch.new_root(space::CirclePackingSearchSpace) = CirclePackingNode(1, 0, 0.0, space.root_state)
+NMK.TreeSearch.stop(space::CirclePackingSearchSpace, untreated_nodes) = space.nb_nodes_evaluated >= space.node_limit || isempty(untreated_nodes)
+
+function NMK.TreeSearch.children(space::CirclePackingSearchSpace, current::CirclePackingNode)
+    println("--- node depth = $(current.depth)")
+    NMK.MathOptState.apply_change!(JuMP.backend(space.model), NMK.MathOptState.forward(current.state), space.helper)
+    space.nb_nodes_evaluated += 1
+    depth = current.depth + 1
+    
+    # Assign a unique ID to the current node if it doesn't have one already
+    current_node_id = current.id
+    if current_node_id == 0
+        current_node_id = space.next_node_id
+        space.next_node_id += 1
+    end
+    
+    # Create node label for DOT file
+
+    MOI.set(space.model, MOI.Silent(), true)
+    optimize!(space.model)
+    @show JuMP.termination_status(space.model)
+
+    obj_value = round(MOI.get(space.model, MOI.ObjectiveValue()), digits=2)
+    node_label = "Node $(current_node_id)\nDepth: $(current.depth)\nObj: $(obj_value)"
+    space.dot_nodes[current_node_id] = node_label
+
+
+    term_status = JuMP.termination_status(space.model)
+    if term_status == MOI.INFEASIBLE || term_status == MOI.LOCALLY_INFEASIBLE
+        # Update node label to indicate infeasibility
+        space.dot_nodes[current_node_id] = "Node $(current_node_id)\nDepth: $(current.depth)\nINFEASIBLE"
+        
+        NMK.MathOptState.apply_change!(JuMP.backend(space.model), NMK.MathOptState.backward(current.state), space.helper)
+        return CirclePackingNode[]
+    end
+
+    result = CirclePackingResult(current_node_id, space.model, space.problem)
+    branching_scores = [(var_id, most_fractional_value(space, var_id)) for var_id in space.helper.original_binary_vars]
+  
+    best_candidate_score, best_candidate_pos = findmax(Iterators.map(elem -> last(elem), branching_scores))
+    best_candidate_var_id = first(branching_scores[best_candidate_pos])
+
+    @show best_candidate_var_id, best_candidate_score
+    children = CirclePackingNode[]
+
+    if best_candidate_score > 1e-3
+        # Get variable name for better labeling
+        var_name = get(space.var_names, best_candidate_var_id, string(best_candidate_var_id))
+        
+        left_lb_forward_changes = [NMK.MathOptState.LowerBoundVarChange(best_candidate_var_id, 1)]
+        left_ub_forward_changes = NMK.MathOptState.UpperBoundVarChange[]
+        left_local_forward_change = NMK.MathOptState.DomainChangeDiff(left_lb_forward_changes, left_ub_forward_changes)
+        left_forward_diff = NMK.MathOptState.merge_forward_change_diff(NMK.MathOptState.forward(current.state), left_local_forward_change)
+
+        left_lb_backward_changes = [NMK.MathOptState.LowerBoundVarChange(best_candidate_var_id, 0)]
+        left_ub_backward_changes = NMK.MathOptState.UpperBoundVarChange[]
+        left_local_backward_change = NMK.MathOptState.DomainChangeDiff(left_lb_backward_changes, left_ub_backward_changes)
+        left_backward_diff = NMK.MathOptState.merge_backward_change_diff(NMK.MathOptState.backward(current.state), left_local_backward_change)
+
+        left_state = NMK.MathOptState.new_state(space.tracker, left_forward_diff, left_backward_diff)
+
+        right_lb_forward_changes = NMK.MathOptState.LowerBoundVarChange[]
+        right_ub_forward_changes = [NMK.MathOptState.UpperBoundVarChange(best_candidate_var_id, 0)]
+        right_local_forward_change = NMK.MathOptState.DomainChangeDiff(right_lb_forward_changes, right_ub_forward_changes)
+        right_forward_diff = NMK.MathOptState.merge_forward_change_diff(NMK.MathOptState.forward(current.state), right_local_forward_change)
+
+        right_lb_backward_changes = NMK.MathOptState.LowerBoundVarChange[]
+        right_ub_backward_changes = [NMK.MathOptState.UpperBoundVarChange(best_candidate_var_id, 1)]
+        right_local_backward_change = NMK.MathOptState.DomainChangeDiff(right_lb_backward_changes, right_ub_backward_changes)
+        right_backward_diff = NMK.MathOptState.merge_backward_change_diff(NMK.MathOptState.backward(current.state), right_local_backward_change)
+
+        right_state = NMK.MathOptState.new_state(space.tracker, right_forward_diff, right_backward_diff)
+
+        # Create child nodes with unique IDs
+        left_node_id = space.next_node_id
+        space.next_node_id += 1
+        
+        right_node_id = space.next_node_id
+        space.next_node_id += 1
+        
+        # Create children with a branching algorithm.
+        children = CirclePackingNode[
+            CirclePackingNode(left_node_id, depth, 0, left_state),
+            CirclePackingNode(right_node_id, depth, 0, right_state)
+        ]
+        
+        # Add edges to DOT file representation
+        push!(space.dot_edges, (current_node_id, left_node_id, "$(var_name) = 1"))
+        push!(space.dot_edges, (current_node_id, right_node_id, "$(var_name) = 0"))
+        
+        # Generate DOT file after each branching
+        generate_dot_file(space)
+    else
+        # Update node label to indicate leaf node (optimal solution)
+        space.dot_nodes[current_node_id] = "Node $(current_node_id)\nDepth: $(current.depth)\nObj: $(obj_value)\nOPTIMAL"
+        generate_dot_file(space)
+    end
+
+    NMK.MathOptState.apply_change!(JuMP.backend(space.model), NMK.MathOptState.backward(current.state), space.helper)
+
+    return children
+end
+
+function NMK.TreeSearch.output(space::CirclePackingSearchSpace, untreated_nodes)
+    # Generate final DOT file
+    generate_dot_file(space)
+    return nothing
+end
+
+# Function to generate DOT file representation of the decision tree
+function generate_dot_file(space::CirclePackingSearchSpace)
+    dot_file_path = joinpath(dirname(@__FILE__), "circle_packing_tree.dot")
+    open(dot_file_path, "w") do f
+        # Write DOT file header
+        write(f, "digraph CirclePackingDecisionTree {\n")
+        write(f, "    // Graph settings\n")
+        write(f, "    graph [rankdir=TB, fontname=\"Arial\", splines=true];\n")
+        write(f, "    node [shape=box, style=\"rounded,filled\", fillcolor=lightblue, fontname=\"Arial\"];\n")
+        write(f, "    edge [fontname=\"Arial\"];\n\n")
+        
+        # Write nodes
+        for (node_id, label) in space.dot_nodes
+            # Set different colors for different node types
+            if occursin("INFEASIBLE", label)
+                write(f, "    node$(node_id) [label=\"$(label)\", fillcolor=lightcoral];\n")
+            elseif occursin("OPTIMAL", label)
+                write(f, "    node$(node_id) [label=\"$(label)\", fillcolor=lightgreen];\n")
+            else
+                write(f, "    node$(node_id) [label=\"$(label)\"];\n")
+            end
+        end
+        
+        write(f, "\n")
+        
+        # Write edges
+        for (parent_id, child_id, label) in space.dot_edges
+            write(f, "    node$(parent_id) -> node$(child_id) [label=\"$(label)\"];\n")
+        end
+        
+        # Add note about the branching strategy
+        write(f, "\n    // Notes about the branching process\n")
+        write(f, "    note1 [shape=note, fillcolor=lightyellow, label=\"Branching Strategy:\nMost fractional binary variable\"];\n")
+        write(f, "    note2 [shape=note, fillcolor=lightyellow, label=\"Search Strategy:\nBest-first search based on\nnode.value_guess\"];\n")
+        
+        # Close the graph
+        write(f, "}\n")
+    end
+    
+    println("DOT file generated at: ", dot_file_path)
+    return dot_file_path
+end
+
+struct BestValueSearchStrategy <: NMK.TreeSearch.AbstractBestFirstSearchStrategy end
+
+function NMK.TreeSearch.get_priority(::BestValueSearchStrategy, space::CirclePackingSearchSpace, node::CirclePackingNode)
+    return node.value_guess
+end
+
+"""
+    solve_circle_packing(problem::CirclePackingProblem; time_limit_seconds::Float64=60.0)
+
+Solve the circle packing problem using JuMP and return the solution.
+
+# Arguments
+- `problem::CirclePackingProblem`: The problem instance
+- `time_limit_seconds::Float64=60.0`: Time limit for the solver in seconds
+
+# Returns
+- `Dict`: A dictionary containing the solution information
+"""
+function solve_circle_packing(problem::CirclePackingProblem; time_limit_seconds::Float64=600.0)
+    model = build_jump_model(problem)
+    tracker = NMK.MathOptState.DomainChangeTracker()
+    helper = NMK.MathOptState.transform_model!(tracker, JuMP.backend(model))
+
+    original_state = NMK.MathOptState.root_state(tracker, JuMP.backend(model))
+    relaxed_state = NMK.MathOptState.relax_integrality!(JuMP.backend(model), helper)
+    NMK.MathOptState.recover_state!(JuMP.backend(model), original_state, relaxed_state, helper)
+
+    root_state = NMK.MathOptState.root_state(tracker, JuMP.backend(model))
+
+    space = CirclePackingSearchSpace(problem, model, helper, root_state, tracker; node_limit=1000)
+
+    strategy = BestValueSearchStrategy()
+    result = NMK.TreeSearch.search(strategy, space)
+end
+
+# Main execution
+function main()
+    println("Creating circle packing problem...")
+    problem = create_sample_problem()
+
+    println("Solving circle packing problem...")
+    solution = solve_circle_packing(problem, time_limit_seconds=600.0)
+    @show solution
+end
+
+# Run the example
+main()
diff --git a/src/MathOptState/MathOptState.jl b/src/MathOptState/MathOptState.jl
index e73aec8..7fb99c2 100644
--- a/src/MathOptState/MathOptState.jl
+++ b/src/MathOptState/MathOptState.jl
@@ -198,5 +198,6 @@ function new_state end
 include("var_bounds_state.jl")
 include("cut_rhs_state.jl")
 include("fixed_var_state.jl")
+include("integrality_state.jl")
 
 end # module MathOptState
\ No newline at end of file
diff --git a/src/MathOptState/integrality_state.jl b/src/MathOptState/integrality_state.jl
new file mode 100644
index 0000000..ebce3ee
--- /dev/null
+++ b/src/MathOptState/integrality_state.jl
@@ -0,0 +1,129 @@
+struct IntegralityStateTracker <: AbstractMathOptStateTracker end
+
+function _apply_integrality_relaxation!(backend, var_id, helper)
+    ci = get(helper.map_integer, var_id, nothing)
+    if !isnothing(ci)
+        MOI.delete(backend, ci)
+        delete!(helper.map_integer, var_id)
+    else
+        @warn "Cannot relax a variable that is not integral."
+    end
+    return
+end
+
+function _apply_integrality_restriction!(backend, var_id, helper)
+    ci = get(helper.map_integer, var_id, nothing)
+    if isnothing(ci)
+        ci = MOI.add_constraint(backend, var_id, MOI.Integer())
+        helper.map_integer[var_id] = ci
+    else
+        @warn "Variable is already integral."
+    end
+    return
+end
+
+function _apply_zero_one_relaxation!(backend, var_id, helper)
+    ci = get(helper.map_binary, var_id, nothing)
+    if !isnothing(ci)
+        MOI.delete(backend, ci)
+        delete!(helper.map_binary, var_id)
+
+        # Check bounds
+        ci_lb = get(helper.map_lb, var_id, nothing)
+        if isnothing(ci_lb) || MOI.get(backend, MOI.ConstraintSet(), ci_lb).lower < 0
+            !isnothing(ci_lb) && MOI.delete(backend, ci_lb)
+            new_ci_lb = MOI.add_constraint(backend, var_id, MOI.GreaterThan(0.0))
+            helper.map_lb[var_id] = new_ci_lb
+        end
+
+        ci_ub = get(helper.map_ub, var_id, nothing)
+        if isnothing(ci_ub) || MOI.get(backend, MOI.ConstraintSet(), ci_ub).upper > 1
+            !isnothing(ci_ub) && MOI.delete(backend, ci_ub)
+            new_ci_ub = MOI.add_constraint(backend, var_id, MOI.LessThan(1.0))
+            helper.map_ub[var_id] = new_ci_ub
+        end
+    else
+        @warn "Cannot relax a variable that is not binary."
+    end
+    return
+end
+
+function _apply_zero_one_restriction!(backend, var_id, helper)
+    ci = get(helper.map_binary, var_id, nothing)
+    if isnothing(ci)
+        ci = MOI.add_constraint(backend, var_id, MOI.ZeroOne())
+        helper.map_binary[var_id] = ci
+    else
+        @warn "Variable is already binary."
+    end
+    return
+end
+
+@enum IntegralityChangeType relax_integrality restrict_integrality relax_zero_one restrict_zero_one
+
+struct IntegralityChange
+    var_id::MOI.VariableIndex
+    change::IntegralityChangeType
+end
+
+function apply_change!(backend, change::IntegralityChange, helper::DomainChangeTrackerHelper)
+    if change.change == relax_integrality
+        _apply_integrality_relaxation!(backend, change.var_id, helper)
+    elseif change.change == restrict_integrality
+        _apply_integrality_restriction!(backend, change.var_id, helper)
+    elseif change.change == relax_zero_one
+        _apply_zero_one_relaxation!(backend, change.var_id, helper)
+    elseif change.change == restrict_zero_one
+        _apply_zero_one_restriction!(backend, change.var_id, helper)
+    end
+end
+
+struct IntegralityChangeDiff <: AbstractMathOptStateDiff
+    integrality_changes::Vector{IntegralityChange}
+end
+
+IntegralityChangeDiff() = IntegralityChangeDiff(
+    IntegralityChange[]
+)
+
+function apply_change!(backend, change::IntegralityChangeDiff, helper::DomainChangeTrackerHelper)
+    for change in change.integrality_changes
+        apply_change!(backend, change, helper)
+    end
+    return
+end
+
+# merge_forward_change_diff(
+#     parent_forward_diff::IntegralityChangeDiff,
+#     local_forward_change::IntegralityChangeDiff
+# ) = IntegralityChangeDiff(
+#     vcat(parent_forward_diff.integrality_changes, local_forward_change.integrality_changes)
+# )
+
+# merge_backward_change_diff(
+#     parent_backward_diff::IntegralityChangeDiff,
+#     local_backward_change::IntegralityChangeDiff
+# ) = IntegralityChangeDiff(
+#     vcat(parent_backward_diff.integrality_changes, local_backward_change.integrality_changes)
+# )
+
+new_state(::IntegralityStateTracker, backward::IntegralityChangeDiff, forward::IntegralityChangeDiff) = ModelState(backward, forward)
+
+function relax_integrality!(backend, helper::DomainChangeTrackerHelper)
+    # Relax integrality constraints (Integer and ZeroOne)
+    # Add bounds for the binary variables.
+
+    forward_changes_itr = Iterators.flatten((
+        Iterators.map(var_id -> IntegralityChange(var_id, relax_integrality), keys(helper.map_integer)),
+        Iterators.map(var_id -> IntegralityChange(var_id, relax_zero_one), keys(helper.map_binary))
+    ))
+    forward_changes = IntegralityChangeDiff(collect(forward_changes_itr))
+
+    backward_changes_itr = Iterators.flatten((
+        Iterators.map(var_id -> IntegralityChange(var_id, restrict_integrality), keys(helper.map_integer)),
+        Iterators.map(var_id -> IntegralityChange(var_id, restrict_zero_one), keys(helper.map_binary))
+    ))
+    backward_changes = IntegralityChangeDiff(collect(backward_changes_itr))
+
+    return new_state(IntegralityStateTracker(), forward_changes, backward_changes)
+end
\ No newline at end of file
diff --git a/src/MathOptState/var_bounds_state.jl b/src/MathOptState/var_bounds_state.jl
index ef0b189..8f38ede 100644
--- a/src/MathOptState/var_bounds_state.jl
+++ b/src/MathOptState/var_bounds_state.jl
@@ -9,16 +9,29 @@ This is used to efficiently apply and track changes to variable bounds in an opt
 - `map_lb`: Maps variable indices to their lower bound constraints
 - `map_ub`: Maps variable indices to their upper bound constraints
 - `map_eq`: Maps variable indices to their equality constraints
+- `is_integer`: Boolean indicating whether the variable is integer
+- `is_binary`: Boolean indicating whether the variable is binary
 """
 struct DomainChangeTrackerHelper
     map_lb::Dict{MOI.VariableIndex, MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{Float64}}}
     map_ub::Dict{MOI.VariableIndex, MOI.ConstraintIndex{MOI.VariableIndex, MOI.LessThan{Float64}}}
     map_eq::Dict{MOI.VariableIndex, MOI.ConstraintIndex{MOI.VariableIndex, MOI.EqualTo{Float64}}}
+    map_integer::Dict{MOI.VariableIndex, MOI.ConstraintIndex{MOI.VariableIndex, MOI.Integer}}
+    map_binary::Dict{MOI.VariableIndex, MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}
+
+    # Store variable index of original integer and binary variables.
+    # Required to branch on a relaxed model.
+    original_integer_vars::Set{MOI.VariableIndex}
+    original_binary_vars::Set{MOI.VariableIndex}
     function DomainChangeTrackerHelper()
         return new(
             Dict{MOI.VariableIndex,MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{Float64}}}(),
             Dict{MOI.VariableIndex,MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{Float64}}}(),
-            Dict{MOI.VariableIndex,MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Float64}}}()
+            Dict{MOI.VariableIndex,MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Float64}}}(),
+            Dict{MOI.VariableIndex,MOI.ConstraintIndex{MOI.VariableIndex,MOI.Integer}}(),
+            Dict{MOI.VariableIndex,MOI.ConstraintIndex{MOI.VariableIndex,MOI.ZeroOne}}(),
+            Set{MOI.VariableIndex}(),
+            Set{MOI.VariableIndex}()
         )
     end
 end
@@ -26,60 +39,39 @@ end
 """
     _register_constraints!(helper, vi, ci)
 
-Register a constraint in the appropriate mapping in the helper.
-This is a fallback method that does nothing for constraint types not handled.
+Register a constraint on a single variable in the appropriate mapping in the helper.
 
 # Arguments
 - `helper`: The DomainChangeTrackerHelper to update
 - `vi`: Variable index
-- `ci`: Constraint index
+- `ci`: Constraint index involving only the variable `vi`.
 """
 _register_constraints!(helper, vi, ci) = nothing
 
-"""
-    _register_constraints!(helper, vi::F, ci::MOI.ConstraintIndex{F,S}) where {F<:MOI.VariableIndex,S<:MOI.GreaterThan}
-
-Register a lower bound constraint in the helper.
-
-# Arguments
-- `helper`: The DomainChangeTrackerHelper to update
-- `vi`: Variable index
-- `ci`: Lower bound constraint index
-"""
 function _register_constraints!(helper, vi::F, ci::MOI.ConstraintIndex{F,S}) where {F<:MOI.VariableIndex,S<:MOI.GreaterThan}
     helper.map_lb[vi] = ci
 end
 
-"""
-    _register_constraints!(helper, vi::F, ci::MOI.ConstraintIndex{F,S}) where {F<:MOI.VariableIndex,S<:MOI.LessThan}
-
-Register an upper bound constraint in the helper.
-
-# Arguments
-- `helper`: The DomainChangeTrackerHelper to update
-- `vi`: Variable index
-- `ci`: Upper bound constraint index
-"""
 function _register_constraints!(helper, vi::F, ci::MOI.ConstraintIndex{F,S}) where {F<:MOI.VariableIndex,S<:MOI.LessThan}
     helper.map_ub[vi] = ci
 end
 
-"""
-    _register_constraints!(helper, vi::F, ci::MOI.ConstraintIndex{F,S}) where {F<:MOI.VariableIndex,S<:MOI.EqualTo}
-
-Register an equality constraint in the helper.
-
-# Arguments
-- `helper`: The DomainChangeTrackerHelper to update
-- `vi`: Variable index
-- `ci`: Equality constraint index
-"""
 function _register_constraints!(helper, vi::F, ci::MOI.ConstraintIndex{F,S}) where {F<:MOI.VariableIndex,S<:MOI.EqualTo}
     helper.map_eq[vi] = ci
 end
 
+function _register_constraints!(helper, vi::F, ci::MOI.ConstraintIndex{F,S}) where {F<:MOI.VariableIndex,S<:MOI.Integer}
+    helper.map_integer[vi] = ci
+    push!(helper.original_integer_vars, vi)
+end
+
+function _register_constraints!(helper, vi::F, ci::MOI.ConstraintIndex{F,S}) where {F<:MOI.VariableIndex,S<:MOI.ZeroOne}
+    helper.map_binary[vi] = ci
+    push!(helper.original_binary_vars, vi)
+end
+
 """
-    LowerBoundVarChange <: AbstractAtomicChange
+    LowerBoundVarChange <: AbstractAtomicChange@
 
 Represents a change to the lower bound of a variable.
 
@@ -112,8 +104,10 @@ function apply_change!(backend, change::LowerBoundVarChange, helper::DomainChang
     if isnothing(ci)
         new_ci = MOI.add_constraint(backend, change.var_id, MOI.GreaterThan(change.new_lb))
         helper.map_lb[change.var_id] = new_ci
+        @debug "add constraint $(change.var_id) => $(change.new_lb)"
     else
         MOI.set(backend, MOI.ConstraintSet(), ci, MOI.GreaterThan(change.new_lb))
+        @debug "set constraint to $(change.var_id) => $(change.new_lb)"
     end
     return
 end
@@ -152,8 +146,10 @@ function apply_change!(backend, change::UpperBoundVarChange, helper::DomainChang
     if isnothing(ci)
        new_ci = MOI.add_constraint(backend, change.var_id, MOI.LessThan(change.new_ub))
        helper.map_ub[change.var_id] = new_ci
+       @debug "add constraint $(change.var_id) <= $(change.new_ub)"
     else
         MOI.set(backend, MOI.ConstraintSet(), ci, MOI.LessThan(change.new_ub))
+        @debug "set constraint to $(change.var_id) <= $(change.new_ub)"
     end
     return
 end
@@ -172,6 +168,13 @@ struct DomainChangeDiff <: AbstractMathOptStateDiff
     upper_bounds::Dict{ColId,UpperBoundVarChange}
 end
 
+function DomainChangeDiff(lb_changes::Vector{LowerBoundVarChange}, ub_changes::Vector{UpperBoundVarChange})
+    return DomainChangeDiff(
+        Dict(change.var_id.value => change for change in lb_changes),
+        Dict(change.var_id.value => change for change in ub_changes)
+    )
+end
+
 """
     DomainChangeDiff()
 
@@ -182,7 +185,7 @@ A new empty `DomainChangeDiff`.
 """
 DomainChangeDiff() = DomainChangeDiff(
     Dict{ColId,LowerBoundVarChange}(),
-    Dict{ColId,UpperBoundVarChange}()
+    Dict{ColId,UpperBoundVarChange}(),
 )
 
 """
@@ -261,7 +264,6 @@ function apply_change!(backend, diff::DomainChangeDiff, helper)
     return
 end
 
-
 """
     DomainChangeTracker <: AbstractMathOptStateTracker
 
diff --git a/test/MathOptStateTests/MathOptStateTests.jl b/test/MathOptStateTests/MathOptStateTests.jl
index 5f7bd3d..51fc5e7 100644
--- a/test/MathOptStateTests/MathOptStateTests.jl
+++ b/test/MathOptStateTests/MathOptStateTests.jl
@@ -527,11 +527,210 @@ function test_math_opt_state3()
     end
 end
 
+function test_domain_change_tracker_continuous()
+    model = Model(GLPK.Optimizer)
+    @variable(model, x >= 0)
+    @variable(model, 0 <= y <= 3)
+
+    @objective(model, Min, x + y)
+
+    tracker = NMK.MathOptState.DomainChangeTracker()
+    helper = NMK.MathOptState.transform_model!(tracker, JuMP.backend(model))
+
+    @test JuMP.index(x) in keys(helper.map_lb)
+    @test JuMP.index(y) in keys(helper.map_lb)
+    @test length(helper.map_lb) == 2
+
+    @test JuMP.index(y) in keys(helper.map_ub)
+    @test length(helper.map_ub) == 1
+
+    @test isempty(helper.map_eq)
+    @test isempty(helper.map_binary)
+    @test isempty(helper.map_integer)
+    @test isempty(helper.original_binary_vars)
+    @test isempty(helper.original_integer_vars)
+    return
+end
+
+function test_domain_change_tracker_binary()
+    model = Model(GLPK.Optimizer)
+    @variable(model, x, Bin)
+    @variable(model, 0 <= y <= 1, Bin)
+
+    @objective(model, Min, x + y)
+
+    tracker = NMK.MathOptState.DomainChangeTracker()
+    helper = NMK.MathOptState.transform_model!(tracker, JuMP.backend(model))
+
+    @test JuMP.index(y) in keys(helper.map_lb)
+    @test length(helper.map_lb) == 1
+
+    @test JuMP.index(y) in keys(helper.map_ub)
+    @test length(helper.map_ub) == 1
+
+    @test isempty(helper.map_eq)
+    @test JuMP.index(x) in keys(helper.map_binary)
+    @test JuMP.index(y) in keys(helper.map_binary)
+    @test isempty(helper.map_integer)
+
+    @test JuMP.index(x) in helper.original_binary_vars
+    @test JuMP.index(y) in helper.original_binary_vars
+    @test isempty(helper.original_integer_vars)
+
+    # What happens when we relax ?
+    undo_relax = JuMP.relax_integrality(model)
+
+    @test !MOI.is_valid(JuMP.backend(model), helper.map_binary[JuMP.index(x)])
+    @test !MOI.is_valid(JuMP.backend(model), helper.map_binary[JuMP.index(y)])
+
+    @test MOI.is_valid(JuMP.backend(model), helper.map_lb[JuMP.index(y)])
+    @test MOI.is_valid(JuMP.backend(model), helper.map_ub[JuMP.index(y)])
+    return
+end
+
+function test_domain_change_tracker_integer()
+    model = Model(GLPK.Optimizer)
+    @variable(model, x >= 0, Int)
+    @variable(model, 0 <= y <= 9, Int)
+
+    @objective(model, Min, x + y)
+
+    tracker = NMK.MathOptState.DomainChangeTracker()
+    helper = NMK.MathOptState.transform_model!(tracker, JuMP.backend(model))
+
+    @test JuMP.index(x) in keys(helper.map_lb)
+    @test JuMP.index(y) in keys(helper.map_lb)
+    @test length(helper.map_lb) == 2
+
+    @test JuMP.index(y) in keys(helper.map_ub)
+    @test length(helper.map_ub) == 1
+
+    @test isempty(helper.map_eq)
+    @test isempty(helper.map_binary)
+    @test JuMP.index(x) in keys(helper.map_integer)
+    @test JuMP.index(y) in keys(helper.map_integer)
+
+    @test JuMP.index(x) in helper.original_integer_vars
+    @test JuMP.index(y) in helper.original_integer_vars
+    @test isempty(helper.original_binary_vars)
+
+    # What happens when we relax ?
+    undo_relax = JuMP.relax_integrality(model)
+
+    @test !MOI.is_valid(JuMP.backend(model), helper.map_integer[JuMP.index(x)])
+    @test !MOI.is_valid(JuMP.backend(model), helper.map_integer[JuMP.index(y)])
+    return
+end
+
+function test_relax_integrality()
+    model = JuMP.Model()
+    @variable(model, 1 <= w <= 5)
+    @variable(model, x, Bin)
+    @variable(model, 0 <= y <= 2, Bin)
+    @variable(model, z >= 2, Int)
+
+    tracker = NMK.MathOptState.DomainChangeTracker()
+    helper = NMK.MathOptState.transform_model!(tracker, JuMP.backend(model))
+
+    integral_state = NMK.MathOptState.root_state(tracker, JuMP.backend(model))
+    integrality_relaxation_state = NMK.MathOptState.relax_integrality!(JuMP.backend(model), helper)
+
+    relax_var_2 = false
+    relax_var_3 = false
+    relax_var_4 = false
+
+    for change in integrality_relaxation_state.forward_diff.integrality_changes
+        if change.var_id.value == 2
+            @test change.change == NMK.MathOptState.relax_zero_one
+            relax_var_2 = true
+        elseif change.var_id.value == 3
+            @test change.change == NMK.MathOptState.relax_zero_one
+            relax_var_3 = true
+        elseif change.var_id.value == 4
+            @test change.change == NMK.MathOptState.relax_integrality
+            relax_var_4 = true
+        end
+    end
+
+    @test relax_var_2
+    @test relax_var_3
+    @test relax_var_4
+
+    restrict_var_2 = false
+    restrict_var_3 = false
+    restrict_var_4 = false
+
+    for change in integrality_relaxation_state.backward_diff.integrality_changes
+        if change.var_id.value == 2
+            @test change.change == NMK.MathOptState.restrict_zero_one
+            restrict_var_2 = true
+        elseif change.var_id.value == 3
+            @test change.change == NMK.MathOptState.restrict_zero_one
+            restrict_var_3 = true
+        elseif change.var_id.value == 4
+            @test change.change == NMK.MathOptState.restrict_integrality
+            restrict_var_4 = true
+        end
+    end
+
+    @test restrict_var_2
+    @test restrict_var_3
+    @test restrict_var_4
+
+    NMK.MathOptState.recover_state!(JuMP.backend(model), integral_state, integrality_relaxation_state, helper)
+
+    @test !JuMP.is_binary(w)
+    @test !JuMP.is_integer(w)
+    @test JuMP.has_lower_bound(w) && JuMP.lower_bound(w) == 1
+    @test JuMP.has_upper_bound(w) && JuMP.upper_bound(w) == 5
+    
+    @test !JuMP.is_binary(x)
+    @test !JuMP.is_integer(x)
+    @test JuMP.has_lower_bound(x) && JuMP.lower_bound(x) == 0
+    @test JuMP.has_upper_bound(x) && JuMP.upper_bound(x) == 1
+
+    @test !JuMP.is_binary(y)
+    @test !JuMP.is_integer(y)
+    @test JuMP.has_lower_bound(y) && JuMP.lower_bound(y) == 0
+    @test JuMP.has_upper_bound(y) && JuMP.upper_bound(y) == 1
+
+    @test !JuMP.is_binary(z)
+    @test !JuMP.is_integer(z)
+    @test JuMP.has_lower_bound(z) && JuMP.lower_bound(z) == 2
+    @test !JuMP.has_upper_bound(z)
+
+    NMK.MathOptState.recover_state!(JuMP.backend(model), integrality_relaxation_state, integral_state, helper)
+
+    @test !JuMP.is_binary(w)
+    @test !JuMP.is_integer(w)
+    @test JuMP.has_lower_bound(w) && JuMP.lower_bound(w) == 1
+    @test JuMP.has_upper_bound(w) && JuMP.upper_bound(w) == 5
+    
+    @test JuMP.is_binary(x)
+    @test !JuMP.is_integer(x)
+    @test JuMP.has_lower_bound(x) && JuMP.lower_bound(x) == 0
+    @test JuMP.has_upper_bound(x) && JuMP.upper_bound(x) == 1
+
+    @test JuMP.is_binary(y)
+    @test !JuMP.is_integer(y)
+    @test JuMP.has_lower_bound(y) && JuMP.lower_bound(y) == 0
+    @test JuMP.has_upper_bound(y) && JuMP.upper_bound(y) == 1
+
+    @test !JuMP.is_binary(z)
+    @test JuMP.is_integer(z)
+    @test JuMP.has_lower_bound(z) && JuMP.lower_bound(z) == 2
+    @test !JuMP.has_upper_bound(z)
+end
+
 function run()
     @testset "MathOptStateTests" begin
         test_math_opt_state1()
         test_math_opt_state2()
         test_math_opt_state3()
+        test_domain_change_tracker_continuous()
+        test_domain_change_tracker_binary()
+        test_domain_change_tracker_integer()
+        test_relax_integrality()
     end
 end
 end # end module
\ No newline at end of file