diff --git a/src/library/format/bqpjson/printer.jl b/src/library/format/bqpjson/printer.jl
index e7a56ea..ed358ac 100644
--- a/src/library/format/bqpjson/printer.jl
+++ b/src/library/format/bqpjson/printer.jl
@@ -1,4 +1,4 @@
-function write_model(io::IO, model::AbstractModel{V}, fmt::BQPJSON) where {V<:Integer}
+function write_model(io::IO, model::AbstractModel{V}, fmt::BQPJSON) where {V}
     if fmt.version === v"1.0.0"
         _print_bqpjson_model_v1_0_0(io, model, fmt)
     else
@@ -8,13 +8,13 @@ function write_model(io::IO, model::AbstractModel{V}, fmt::BQPJSON) where {V<:In
     return nothing
 end
 
-function _print_bqpjson_model_v1_0_0(io::IO, model::AbstractModel{V}, fmt::BQPJSON) where {V<:Integer}
+function _print_bqpjson_model_v1_0_0(io::IO, model::AbstractModel{V}, fmt::BQPJSON) where {V}
     json_data = Dict{String,Any}(
         "id"              => 0,
         "variable_domain" => _BQPJSON_VARIABLE_DOMAIN(domain(model)),
         "linear_terms"    => Dict{String,Any}[],
         "quadratic_terms" => Dict{String,Any}[],
-        "variable_ids"    => variables(model),
+        "variable_ids"    => indices(model),
         "scale"           => scale(model),
         "offset"          => offset(model),
         "metadata"        => Dict{String,Any}(),
@@ -24,7 +24,7 @@ function _print_bqpjson_model_v1_0_0(io::IO, model::AbstractModel{V}, fmt::BQPJS
     for (i, l) in linear_terms(model)
         push!(
             json_data["linear_terms"],
-            Dict{String,Any}("id" => variable(model, i), "coeff" => l),
+            Dict{String,Any}("id" => i, "coeff" => l),
         )
     end
 
@@ -32,8 +32,8 @@ function _print_bqpjson_model_v1_0_0(io::IO, model::AbstractModel{V}, fmt::BQPJS
         push!(
             json_data["quadratic_terms"],
             Dict{String,Any}(
-                "id_head" => variable(model, i),
-                "id_tail" => variable(model, j),
+                "id_head" => i,
+                "id_tail" => j,
                 "coeff"   => q,
             ),
         )
diff --git a/src/library/format/minizinc/printer.jl b/src/library/format/minizinc/printer.jl
index 7690f2b..f885bc8 100644
--- a/src/library/format/minizinc/printer.jl
+++ b/src/library/format/minizinc/printer.jl
@@ -1,4 +1,4 @@
-function write_model(io::IO, model::AbstractModel{V}, fmt::MiniZinc) where {V<:Integer}
+function write_model(io::IO, model::AbstractModel{V}, fmt::MiniZinc) where {V}
     _print_metadata(io, model, fmt)
     _print_domain(io, model, fmt)
     _print_variables(io, model, fmt)
@@ -32,9 +32,7 @@ end
 
 function _print_variables(io::IO, model::AbstractModel, ::MiniZinc)
     for i = indices(model)
-        k = variable(model, i)
-        
-        println(io, "var Domain: x$(k);")
+        println(io, "var Domain: x$(i);")
     end
 
     return nothing
@@ -47,16 +45,11 @@ function _print_objective(io::IO, model::AbstractModel, ::MiniZinc)
     println(io, "float: offset = $(offset(model));")
 
     for (i, v) in linear_terms(model)
-        xi = "x$(variable(model, i))"
-
-        push!(objective_terms, "$(v)*$(xi)")
+        push!(objective_terms, "$(v)*x$(i)")
     end
 
     for ((i, j), v) in quadratic_terms(model)
-        xi = "x$(variable(model, i))"
-        xj = "x$(variable(model, j))"
-
-        push!(objective_terms, "$(v)*$(xi)*$(xj)")
+        push!(objective_terms, "$(v)*x$(i)*x$(j)")
     end
 
     if !isempty(objective_terms)
diff --git a/src/library/format/qubo/parser.jl b/src/library/format/qubo/parser.jl
index 3886f8d..41391c0 100644
--- a/src/library/format/qubo/parser.jl
+++ b/src/library/format/qubo/parser.jl
@@ -17,7 +17,7 @@ function read_model(io::IO, fmt::QUBO)
     end
 
     return Model{Int,Float64,Int}(
-        Set{Int}(0:data[:dimension]-1),
+        Set{Int}(1:data[:dimension]),
         data[:linear_terms],
         data[:quadratic_terms];
         scale       = data[:scale],
@@ -55,8 +55,8 @@ function _parse_entry!(
         return false
     end
 
-    i = parse(Int, m[1])
-    j = parse(Int, m[2])
+    i = parse(Int, m[1]) + 1
+    j = parse(Int, m[2]) + 1
     c = parse(Float64, m[3])
 
     if i == j
@@ -78,8 +78,8 @@ function _parse_entry!(data::Dict{Symbol,Any}, line::AbstractString, ::QUBO, ::V
         return false
     end
 
-    i = parse(Int, m[1])
-    j = parse(Int, m[2])
+    i = parse(Int, m[1]) + 1
+    j = parse(Int, m[2]) + 1
     c = parse(Float64, m[3])
 
     if i == j
diff --git a/src/library/format/qubo/printer.jl b/src/library/format/qubo/printer.jl
index e69afbf..3cb9d0d 100644
--- a/src/library/format/qubo/printer.jl
+++ b/src/library/format/qubo/printer.jl
@@ -1,7 +1,7 @@
 function write_model(io::IO, model::AbstractModel{V,T,U}, fmt::QUBO) where {V,T,U}
     data = Dict{Symbol,Any}(
-        :linear_terms    => Dict{Int,T}(variable(model, i) => v for (i, v) in linear_terms(model)),
-        :quadratic_terms => Dict{Tuple{Int,Int},T}((variable(model, i), variable(model, j)) => v for ((i, j), v) in quadratic_terms(model)),
+        :linear_terms    => Dict{Int,T}(linear_terms(model)),
+        :quadratic_terms => Dict{Tuple{Int,Int},T}(quadratic_terms(model)),
         :linear_size     => linear_size(model),
         :quadratic_size  => quadratic_size(model),
         :scale           => scale(model),
@@ -87,11 +87,11 @@ end
 
 function _print_entries(io::IO, data::Dict{Symbol,Any}, ::QUBO, ::Val{_}) where {_}
     for (i, l) in data[:linear_terms]
-        println(io, "$(i) $(i) $(l)")
+        println(io, "$(i-1) $(i-1) $(l)")
     end
 
     for ((i, j), q) in data[:quadratic_terms]
-        println(io, "$(i) $(j) $(q)")
+        println(io, "$(i-1) $(j-1) $(q)")
     end
 
     return nothing
@@ -101,13 +101,13 @@ function _print_entries(io::IO, data::Dict{Symbol,Any}, ::QUBO, ::Val{:dwave})
     println(io, "c linear terms")
 
     for (i, l) in data[:linear_terms]
-        println(io, "$(i) $(i) $(l)")
+        println(io, "$(i-1) $(i-1) $(l)")
     end
 
     println(io, "c quadratic terms")
 
     for ((i, j), q) in data[:quadratic_terms]
-        println(io, "$(i) $(j) $(q)")
+        println(io, "$(i-1) $(j-1) $(q)")
     end
 
     return nothing
@@ -117,7 +117,7 @@ function _print_entries(io::IO, data::Dict{Symbol,Any}, ::QUBO, ::Val{:mqlib})
     println(io, "# linear terms")
 
     for (i, l) in data[:linear_terms]
-        println(io, "$(i) $(i) $(l)")
+        println(io, "$(i-1) $(i-1) $(l)")
     end
 
     println(io, "# quadratic terms")
@@ -125,7 +125,7 @@ function _print_entries(io::IO, data::Dict{Symbol,Any}, ::QUBO, ::Val{:mqlib})
     for ((i, j), q) in data[:quadratic_terms]
         # NOTE: in MQLib qubo files, quadratic coefficients
         # are halved when written to the file
-        println(io, "$(i) $(j) $(q/2)")
+        println(io, "$(i-1) $(j-1) $(q/2)")
     end
 
     return nothing
diff --git a/src/library/model/model.jl b/src/library/model/model.jl
index 5d403a9..9fe1526 100644
--- a/src/library/model/model.jl
+++ b/src/library/model/model.jl
@@ -344,3 +344,13 @@ function Model{V,T,U}(f::F; kws...) where {V,T,U,F<:PBO.AbstractFunction{V,T}}
 
     return Model{V,T,U}(L, Q; offset = β, sense = :min, domain = :bool, kws...)
 end
+
+function map_variables(::Type{V}, vm::Function, model::AbstractModel{_,T,U}) where {_,V,T,U}
+    new_model = copy(model)::AbstractModel{V,T,U}
+    new_model.variable_map = VariableMap{V}(Dict{Int,V}(i => vm(i)::V for i in indices(model)))
+
+    return new_model
+end
+
+map_variables(vm::Dict{Int,V}, model::AbstractModel) where {V}       = map_variables(V, i -> vm[i], model)
+map_variables(vm::AbstractVector{V}, model::AbstractModel) where {V} = map_variables(V, i -> vm[i], model)
diff --git a/src/library/model/variable_map.jl b/src/library/model/variable_map.jl
index 14af997..0441cb6 100644
--- a/src/library/model/variable_map.jl
+++ b/src/library/model/variable_map.jl
@@ -17,6 +17,24 @@ struct VariableMap{V}
     end
 end
 
+function VariableMap{V}(vm::Dict{Int,V}) where {V}
+    map = sizehint!(Dict{V,Int}(), length(vm))
+    inv = Vector{V}(undef, length(vm))
+
+    for i = 1:length(vm)
+        if !haskey(vm, i)
+            error("Invalid variable mapping: Mappings should contain values for all indices")
+        else
+            let v = vm[i]
+                map[v] = i
+                inv[i] = v
+            end
+        end
+    end
+
+    return VariableMap{V}(map, inv)
+end
+
 function VariableMap{V}(
     variables::X,
 ) where {V,X<:Union{AbstractVector{V},AbstractSet{V}}}
diff --git a/test/unit/library/formats.jl b/test/unit/library/formats.jl
index 5ea8322..7f63397 100644
--- a/test/unit/library/formats.jl
+++ b/test/unit/library/formats.jl
@@ -30,16 +30,19 @@ function test_bqpjson_format()
     @testset "⋅ BQPJSON" begin
         @testset "bool" begin
             for i = 0:2
-                file_path = joinpath(__TEST_PATH__, "data", Printf.@sprintf("%02d", i), "bool.json")
+                file_path =
+                    joinpath(__TEST_PATH__, "data", Printf.@sprintf("%02d", i), "bool.json")
                 temp_path = "$(tempname()).bool.json"
 
                 src_model = QUBOTools.read_model(file_path)
+                variables = QUBOTools.variables(src_model)
 
                 @test src_model isa QUBOTools.Model
 
                 QUBOTools.write_model(temp_path, src_model)
 
-                dst_model = QUBOTools.read_model(temp_path)
+                dst_model =
+                    QUBOTools.map_variables(variables, QUBOTools.read_model(temp_path))
 
                 @test dst_model isa QUBOTools.Model
 
@@ -49,16 +52,19 @@ function test_bqpjson_format()
 
         @testset "spin" begin
             for i = 0:2
-                file_path = joinpath(__TEST_PATH__, "data", Printf.@sprintf("%02d", i), "spin.json")
+                file_path =
+                    joinpath(__TEST_PATH__, "data", Printf.@sprintf("%02d", i), "spin.json")
                 temp_path = "$(tempname()).spin.json"
 
                 src_model = QUBOTools.read_model(file_path)
+                variables = QUBOTools.variables(src_model)
 
                 @test src_model isa QUBOTools.Model
 
                 QUBOTools.write_model(temp_path, src_model)
 
-                dst_model = QUBOTools.read_model(temp_path)
+                dst_model =
+                    QUBOTools.map_variables(variables, QUBOTools.read_model(temp_path))
 
                 @test dst_model isa QUBOTools.Model
 
@@ -75,17 +81,22 @@ function test_qubo_format()
         src_fmt = QUBOTools.QUBO(:dwave)
 
         for i = 0:2
-            file_path = joinpath(__TEST_PATH__, "data", Printf.@sprintf("%02d", i), "bool.qubo")
+            file_path =
+                joinpath(__TEST_PATH__, "data", Printf.@sprintf("%02d", i), "bool.qubo")
             temp_path = "$(tempname()).bool.qubo"
 
             src_model = QUBOTools.read_model(file_path, src_fmt)
+            variables = QUBOTools.variables(src_model)
 
             @test src_model isa QUBOTools.Model
 
             for dst_fmt in QUBOTools.QUBO.([:dwave, :mqlib])
                 QUBOTools.write_model(temp_path, src_model, dst_fmt)
 
-                dst_model = QUBOTools.read_model(temp_path, dst_fmt)
+                dst_model = QUBOTools.map_variables(
+                    variables,
+                    QUBOTools.read_model(temp_path, dst_fmt),
+                )
 
                 @test dst_model isa QUBOTools.Model
 
@@ -100,16 +111,18 @@ end
 function test_qubist_format()
     @testset "⋅ Qubist" begin
         for i = 0:2
-            file_path = joinpath(__TEST_PATH__, "data", Printf.@sprintf("%02d", i), "spin.qh")
+            file_path =
+                joinpath(__TEST_PATH__, "data", Printf.@sprintf("%02d", i), "spin.qh")
             temp_path = "$(tempname()).spin.qh"
 
             src_model = QUBOTools.read_model(file_path)
+            variables = QUBOTools.variables(src_model)
 
             @test src_model isa QUBOTools.Model
 
             QUBOTools.write_model(temp_path, src_model)
 
-            dst_model = QUBOTools.read_model(temp_path)
+            dst_model = QUBOTools.map_variables(variables, QUBOTools.read_model(temp_path))
 
             @test dst_model isa QUBOTools.Model
 
@@ -124,16 +137,19 @@ function test_qubin_format()
     @testset "⋅ QUBin" begin
         @testset "bool" begin
             for i = 0:2
-                file_path = joinpath(__TEST_PATH__, "data", Printf.@sprintf("%02d", i), "bool.qb")
+                file_path =
+                    joinpath(__TEST_PATH__, "data", Printf.@sprintf("%02d", i), "bool.qb")
                 temp_path = "$(tempname()).bool.qb"
 
                 src_model = QUBOTools.read_model(file_path)
+                variables = QUBOTools.variables(src_model)
 
                 @test src_model isa QUBOTools.Model
 
                 QUBOTools.write_model(temp_path, src_model)
 
-                dst_model = QUBOTools.read_model(temp_path)
+                dst_model =
+                    QUBOTools.map_variables(variables, QUBOTools.read_model(temp_path))
 
                 @test dst_model isa QUBOTools.Model
 
@@ -148,12 +164,14 @@ function test_qubin_format()
                 temp_path = "$(tempname()).spin.qb"
 
                 src_model = QUBOTools.read_model(file_path)
+                variables = QUBOTools.variables(src_model)
 
                 @test src_model isa QUBOTools.Model
 
                 QUBOTools.write_model(temp_path, src_model)
 
-                dst_model = QUBOTools.read_model(temp_path)
+                dst_model =
+                    QUBOTools.map_variables(variables, QUBOTools.read_model(temp_path))
 
                 @test dst_model isa QUBOTools.Model
 
@@ -175,7 +193,7 @@ function test_minizinc_format()
                         (1, 3) => -13.0,
                         (2, 3) => -23.0,
                     );
-                    scale  =  2.0,
+                    scale  = 2.0,
                     offset = -1.0,
                     sense  = :min,
                     domain = :bool,
@@ -206,7 +224,7 @@ function test_minizinc_format()
                         (1, 3) => -13.0,
                         (2, 3) => -23.0,
                     );
-                    scale  =  2.0,
+                    scale  = 2.0,
                     offset = -1.0,
                     sense  = :max,
                     domain = :spin,