diff --git a/docs/examples/fastcore.md b/docs/examples/fastcore.md
index b1c99db..5d53aff 100644
--- a/docs/examples/fastcore.md
+++ b/docs/examples/fastcore.md
@@ -10,13 +10,11 @@ Fastcore defines a set of core reactions that is forced to be active (carry a no
 The Fastcore algorithm is greedy approximation of the exact problem which can be modelled as a MIP problem. With MIOM, the exact problem can be defined in a few lines, using the `opt_tol` parameter to control the level of optimality required:
 
 ```python
-from miom import miom, Solvers
-from miom.mio import load_gem
-from miom.miom import Comparator, ExtractionMode
+import miom
 import numpy as np
 
 # Use the flux-consistent subnetwork (fcm) of the Human1 GEM model 
-m = load_gem('https://github.com/pablormier/miom-gems/raw/main/gems/homo_sapiens_human1_fcm.miom')
+m = miom.mio.load_gem('https://github.com/pablormier/miom-gems/raw/main/gems/homo_sapiens_human1_fcm.miom')
 # Select reactions from the cholesterol metabolism as the core reactions to keep
 core_rxn = m.find_reactions_from_pathway("Cholesterol metabolism")
 print(sum(core_rxn))
@@ -26,18 +24,19 @@ weights = -1 * np.ones(m.num_reactions)
 weights[core_rxn == 1] = 1
 
 # Exact-Fastcore
-fmc = (miom(m, solver=Solvers.GUROBI_PYMIP)
-          .setup(opt_tol=0.01)
-          .steady_state()
-          .subset_selection(weights)
-          .keep(np.where(core_rxn == 1)[0])
-          .solve(verbosity=1)
-          .select_subnetwork(
-               mode=ExtractionMode.ABSOLUTE_FLUX_VALUE,
-               comparator=Comparator.GREATER_OR_EQUAL,
-               value=1e-8
-          )
-          .network
+fmc = (miom
+        .load(m, solver=miom.Solvers.GUROBI_PYMIP)
+        .setup(opt_tol=0.01)
+        .steady_state()
+        .subset_selection(weights)
+        .keep(core_rxn == 1)
+        .solve(verbosity=1)
+        .select_subnetwork(
+            mode=miom.ExtractionMode.ABSOLUTE_FLUX_VALUE,
+            comparator=miom.Comparator.GREATER_OR_EQUAL,
+            value=1e-8
+        )
+        .network
 )
 print(fmc.num_reactions)
 ```
diff --git a/docs/examples/fba.md b/docs/examples/fba.md
index abe5960..b0d2bfd 100644
--- a/docs/examples/fba.md
+++ b/docs/examples/fba.md
@@ -3,10 +3,11 @@
 A simple Flux Balance Analysis can be easily defined and solved with any of the commercial or open-source solvers available:
 
 ```python
-from miom import miom, load_gem
-network = load_gem("https://github.com/pablormier/miom-gems/raw/main/gems/mus_musculus_iMM1865.miom")
+import miom
+network = miom.mio.load_gem("https://github.com/pablormier/miom-gems/raw/main/gems/mus_musculus_iMM1865.miom")
 target = "BIOMASS_reaction"
-flux = (miom(network)
+flux = (miom
+        .load(network)
         .steady_state()
         .set_rxn_objective(target)
         .solve(verbosity=1)
diff --git a/docs/examples/imat.md b/docs/examples/imat.md
index 6299a69..f26ccec 100644
--- a/docs/examples/imat.md
+++ b/docs/examples/imat.md
@@ -6,12 +6,10 @@ Example implementation of iMAT using MIOM. Note that this implementation support
     Shlomi, T., Cabili, M. N., Herrgård, M. J., Palsson, B. Ø., & Ruppin, E. (2008). Network-based prediction of human tissue-specific metabolism. [Nature biotechnology, 26(9), 1003-1010](https://www.nature.com/articles/nbt.1487).
 
 ```python
-from miom import miom, Solvers
-from miom.mio import load_gem
-from miom.miom import Comparator, ExtractionMode
+import miom
 
 # Use the iHuman-GEM model
-m = load_gem('https://github.com/pablormier/miom-gems/raw/main/gems/homo_sapiens_human1.miom')
+m = miom.mio.load_gem('https://github.com/pablormier/miom-gems/raw/main/gems/homo_sapiens_human1.miom')
 
 # Add all the reactions from the Cholesterol pathway to the highly expressed set
 RH = m.find_reactions_from_pathway("Cholesterol metabolism")
@@ -20,14 +18,15 @@ RL = -1 * m.find_reactions_from_pathway("Pyruvate metabolism")
 w = RH + RL
 print("RH:", sum(RH), "RL:", sum(abs(RL)))
 
-m = (miom(m, solver=Solvers.GUROBI)
+m = (miom
+     .load(m, solver=miom.Solvers.GUROBI)
      .setup(int_tol=1e-8, opt_tol=0.01, verbosity=1)
      .steady_state()
      .subset_selection(w)
      .solve(max_seconds=30)
      .select_subnetwork(
-          mode=ExtractionMode.ABSOLUTE_FLUX_VALUE,
-          comparator=Comparator.GREATER_OR_EQUAL,
+          mode=miom.ExtractionMode.ABSOLUTE_FLUX_VALUE,
+          comparator=miom.Comparator.GREATER_OR_EQUAL,
           value=1e-8
      )
      .network)
diff --git a/docs/examples/sparse_fba.md b/docs/examples/sparse_fba.md
index a8af6ee..2f80ffb 100644
--- a/docs/examples/sparse_fba.md
+++ b/docs/examples/sparse_fba.md
@@ -7,18 +7,19 @@ Toolbox includes different LP heuristics to minimize different approximations of
 Here is the implementation of sparse-FBA with MIOM:
 
 ```python
-from miom import miom, load_gem, Solvers
+import miom
 
 # Load a genome-scale metabolic network. You can load SMBL or Matlab metabolic networks
 # as well using the same method, but it requires to have the cobratoolbox python library
 # installed.
-network = load_gem("https://github.com/pablormier/miom-gems/raw/main/gems/mus_musculus_iMM1865.miom")
+network = miom.mio.load_gem("https://github.com/pablormier/miom-gems/raw/main/gems/mus_musculus_iMM1865.miom")
 
 # Create the sparse FBA problem to get a solution that maximizes
 # the optimal flux through the BIOMASS_reaction minimizing the
 # number of active reactions. The solution should be not more than
 # 5% of the optimal solution (opt_tol = 0.05).
-V, X = (miom(network, solver=Solvers.GUROBI_PYMIP)
+V, X = (miom
+        .load(network, solver=miom.Solvers.GUROBI_PYMIP)
         # Set-up the solver options
         .setup(int_tol=1e-8, opt_tol=0.05, verbosity=1)
         # Add the steady-state constraints (S*V = 0)
@@ -51,7 +52,7 @@ V, X = (miom(network, solver=Solvers.GUROBI_PYMIP)
         .get_values())
 
 # Show reactions with a flux > 1e-7
-print("Number of reactions with flux above +/- 1e-7:", sum(abs(V)>1e-7))
+print("Number of reactions with flux above +/- 1e-8:", sum(abs(V) > 1e-8))
 
 # Count reactions with an indicator value of 0 (active). Note that since
 # the weights of the reactions are negative (for all rxns), an indicator
diff --git a/docs/index.md b/docs/index.md
index 771c5ea..24a752d 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -5,12 +5,10 @@
 
 
 
-__MIOM__ (Mixed Integer Optimization for Metabolism) is a python library for creating and solving complex optimization problems using genome-scale metabolic networks, in just a few lines. 
+__MIOM__ (Mixed Integer Optimization for Metabolism) is a python library for creating and solving complex constraint-based optimization problems for metabolism, in just a few lines. 
 
 By leveraging the power of modern Mixed Integer Optimization (MIO) solvers, it can transform any simple constraint-based problem, such as Flux Balance Analysis (FBA) into a more complex optimization problem, such as sparse FBA or context-specific reconstruction problems very easily, and solve it with the __required level of optimality__.
 
-But what is even better is that most of the time, algorithms formulated as Mixed Integer Optimization problems with MIOM can be solved faster and with better quality than currently existing alternatives that are approximations of the original problem. By using the MIO formulation, you can get also an estimation of how close to optimality a solution is, so you don't need to waste more time than needed.
-
 MIOM uses the [PICOS](https://picos-api.gitlab.io/picos/) and the [Python-MIP](https://www.python-mip.com/) libraries to build and solve the optimization problems using many commercial, academic and free solvers.
 
 !!! warning
@@ -37,17 +35,18 @@ CPLEX is also supported, but requires a license. To install MIOM with CPLEX supp
 Here is an example of how to load a metabolic network and maximize the flux through a target reaction using FBA, and then how to modify the original problem to implement the sparse FBA problem adding only a few lines to the original problem:
 
 ```python
-from miom import miom, load_gem, Solvers
+import miom
 
 # Load a genome-scale metabolic network using the miom format. 
 # You can load SMBL or Matlab metabolic networks as well using 
 # the same method, but it requires to have the cobratoolbox python library
 # installed (and scipy for mat files). To install these dependencies, run:
 # $ pip install cobra scipy
-network = load_gem("https://github.com/pablormier/miom-gems/raw/main/gems/mus_musculus_iMM1865.miom")
+network = miom.mio.load_gem("https://github.com/pablormier/miom-gems/raw/main/gems/mus_musculus_iMM1865.miom")
 target_rxn = "BIOMASS_reaction"
 # Create the optimization problem with miom and solve
-model = (miom(network)
+model = (miom
+        .load(network)
         .steady_state()
         .set_rxn_objective(target_rxn)
         .solve(verbosity=1))
@@ -109,7 +108,8 @@ approximate solution, controlled by the `opt_tol` parameter.
 To use other solvers, you only need to provide the specific solver to the `miom` method, for example:
 
 ```python
-model = (miom(network, solver=Solvers.GLPK)
+model = (miom
+        .load(network, solver=Solvers.GLPK)
         .steady_state()
         .set_rxn_objective(target_rxn)
         .solve(verbosity=1))
diff --git a/examples/exact_fastcore_example.py b/examples/exact_fastcore_example.py
index 4cb00c2..9a3ee81 100644
--- a/examples/exact_fastcore_example.py
+++ b/examples/exact_fastcore_example.py
@@ -1,27 +1,27 @@
-from miom import miom, Solvers
-from miom.mio import load_gem
-from miom.miom import Comparator, ExtractionMode
+import miom
 import numpy as np
 
 # Load a model
-m = load_gem('https://github.com/pablormier/miom-gems/raw/main/gems/homo_sapiens_human1_fcm.miom')
+m = miom.mio.load_gem('https://github.com/pablormier/miom-gems/raw/main/gems/homo_sapiens_human1_fcm.miom')
 print("Reactions in the consistent iHuman GEM:", m.num_reactions)
 core_rxn = m.find_reactions_from_pathway("Cholesterol metabolism")
 print("Num. of core reactions:", sum(core_rxn))
-# Assign a negative weight for reactions not in the core
+
+# Assign a negative weight for reactions not in the core and a positive value
+# for reactions in the core.
 weights = -1 * np.ones(m.num_reactions)
 weights[core_rxn == 1] = 1
 
-# Fastcore
-fmc = (miom(m, solver=Solvers.GUROBI_PYMIP)
+# Weighted Exact Fastcore algorithm with MIOM:
+fmc = (miom.load(m, solver=miom.Solvers.COIN_OR_CBC)
           .setup(opt_tol=0.05)
           .steady_state()
           .subset_selection(weights)
-          .keep(np.where(core_rxn == 1)[0])
+          .keep(core_rxn == 1)
           .solve(verbosity=1)
           .select_subnetwork(
-               mode=ExtractionMode.ABSOLUTE_FLUX_VALUE,
-               comparator=Comparator.GREATER_OR_EQUAL,
+               mode=miom.ExtractionMode.ABSOLUTE_FLUX_VALUE,
+               comparator=miom.Comparator.GREATER_OR_EQUAL,
                value=1e-8
           )
           .network
diff --git a/examples/fba_example.py b/examples/fba_example.py
index 67328b7..bfa96f1 100644
--- a/examples/fba_example.py
+++ b/examples/fba_example.py
@@ -1,11 +1,10 @@
-from miom import load_gem, miom, Solvers
-import numpy as np
+import miom
 
 # Load a model
-m = load_gem('https://github.com/pablormier/miom-gems/raw/main/gems/mus_musculus_iMM1865.miom')
+m = miom.mio.load_gem('https://github.com/pablormier/miom-gems/raw/main/gems/mus_musculus_iMM1865.miom')
 
 target = 'BIOMASS_reaction'
-opt_flux = (miom(network=m, solver=Solvers.GLPK)
+opt_flux = (miom.load(network=m, solver=miom.Solvers.GLPK)
                 .steady_state()
                 .set_rxn_objective(target)
                 .solve(verbosity=1)
diff --git a/examples/imat_example.py b/examples/imat_example.py
index 4b78cd0..48e0c77 100644
--- a/examples/imat_example.py
+++ b/examples/imat_example.py
@@ -18,7 +18,7 @@
 w = RH + RL
 print("RH:", sum(RH), "RL:", sum(abs(RL)))
 
-m = (miom(m, solver=Solvers.GUROBI)
+m = (miom(m, solver=Solvers.COIN_OR_CBC)
      .setup(int_tol=1e-8, opt_tol=0.01, verbosity=1)
      .steady_state()
      .subset_selection(w)
@@ -29,5 +29,5 @@
           value=1e-8
      )
      .network)
-print(m.num_reactions)
+print("Number of reactions in the subnetwork with active flux:", m.num_reactions)
 
diff --git a/examples/introduction_example.py b/examples/introduction_example.py
index c2cbbb0..ca32be7 100644
--- a/examples/introduction_example.py
+++ b/examples/introduction_example.py
@@ -1,4 +1,5 @@
-from miom import miom, load_gem, Solvers
+from miom import miom, Solvers
+from miom.mio import load_gem
 
 # Load a genome-scale metabolic network. You can load SMBL or Matlab metabolic networks
 # as well using the same method, but it requires to have the cobratoolbox python library
@@ -11,8 +12,9 @@
         .set_rxn_objective(target_rxn)
         .solve(verbosity=1))
 print("Optimal flux:", model.get_fluxes(target_rxn), "mmol/(h·gDW)")
+
 # Show reactions with non-zero flux
-V, _ = model.get_values()
+V = model.get_fluxes()
 print("Number of reactions with flux above +/- 1e-8:", sum(abs(V) > 1e-8))
 
 V, X = (model
@@ -31,6 +33,6 @@
         .get_values())
 
 print("Number of reactions with an absolute flux value above 1e-8:", sum(abs(V) > 1e-8))
-print("Active reactions:", sum(1 if activity == 1 else 0 for activity in model.variables.indicator_rxn_activity))
-print("Inconsistencies:", sum(1 if activity != activity else 0 for activity in model.variables.indicator_rxn_activity))
+print("Active reactions:", sum(1 if activity == 1 else 0 for activity in model.variables.reaction_activity))
+print("Inconsistencies:", sum(1 if activity != activity else 0 for activity in model.variables.reaction_activity))
 print("Solver status:", model.get_solver_status())
\ No newline at end of file
diff --git a/examples/sparse_fba_example.py b/examples/sparse_fba_example.py
index a20619d..5208671 100644
--- a/examples/sparse_fba_example.py
+++ b/examples/sparse_fba_example.py
@@ -1,11 +1,12 @@
-from miom import miom, load_gem, Solvers
+import miom
 
 # Load a model
-m = load_gem('https://github.com/pablormier/miom-gems/raw/main/gems/mus_musculus_iMM1865.miom')
-print(m.num_reactions)
+m = miom.mio.load_gem('https://github.com/pablormier/miom-gems/raw/main/gems/mus_musculus_iMM1865.miom')
+print("Num. of reactions in the network:", m.num_reactions)
 
 # Solve with Gurobi, CPLEX or CBC (other MIP solvers struggle with mediudm/large networks)
-V, X = (miom(network=m, solver=Solvers.GUROBI_PYMIP)
+V, X = (miom
+        .load(network=m, solver=miom.Solvers.COIN_OR_CBC)
         # Config the solver (e.g. set the optimality tolerance for MIP problems)
         .setup(int_tol=1e-8, opt_tol=0.05)
         # Add steady-state constraints to the model (S * V = 0)
@@ -30,7 +31,7 @@
         .get_values())
 
 # Show reactions with a flux > 1e-7
-print("Number of reactions with flux above +/- 1e-7:", sum(abs(V)>1e-7))
+print("Number of reactions with flux above +/- 1e-8:", sum(abs(V)>1e-8))
 
 # Count reactions with an indicator value of 0 (active). Note that since
 # the weights of the reactions are negative (for all rxns), an indicator
diff --git a/miom/__init__.py b/miom/__init__.py
index b844e2b..298c9fa 100644
--- a/miom/__init__.py
+++ b/miom/__init__.py
@@ -1,11 +1,9 @@
-from miom.miom import (
-    miom,
-    PicosModel, 
-    PythonMipModel, 
-    Solvers
+from miom.miom import ( 
+    load,
+    Solvers,
+    Comparator,
+    ExtractionMode
 )
 
-from miom.mio import (
-    load_gem, 
-    export_gem
-)
+__all__ = ["load", "Solvers", "Comparator", "ExtractionMode"]
+__version__ = "0.9.0-beta.0"
\ No newline at end of file
diff --git a/miom/__main__.py b/miom/__main__.py
new file mode 100644
index 0000000..f5e666a
--- /dev/null
+++ b/miom/__main__.py
@@ -0,0 +1,35 @@
+import argparse
+import miom
+
+def convert_gem(args):
+    print(f"MIOM v{miom.__version__}")
+    m = miom.mio.load_gem(args.input)
+    print(f"Loaded network with {m.num_reactions} reactions")
+    print(f"Exporting to {args.output}...")
+    miom.mio.export_gem(m, args.output)
+    print("Done.")
+
+def get_args():
+    parser = argparse.ArgumentParser(description=f"MIOM: Mixed Integer Optimization for Metabolism")
+    parser.add_argument(
+        "--version",
+        action="version", 
+        version=f"MIOM v{miom.__version__}"
+    )
+    subparsers = parser.add_subparsers(help="sub-command help")
+    convert = subparsers.add_parser("convert", help="Convert a model to miom format")
+    convert.add_argument(
+        "input",
+        help="Input model file (if cobra is installed, any format supported by cobra is allowed)"
+    )
+    convert.add_argument(
+        "output",
+        help="Output GEM file in MIOM format"
+    )
+    convert.set_defaults(func=convert_gem)
+    return parser.parse_args()
+
+if __name__ == '__main__':
+    args = get_args()
+    if args.func:
+        args.func(args)
\ No newline at end of file
diff --git a/miom/miom.py b/miom/miom.py
index e277d09..31dc697 100644
--- a/miom/miom.py
+++ b/miom/miom.py
@@ -148,7 +148,7 @@ class Comparator(str, Enum):
     ])
 
 
-def miom(network, solver=Solvers.COIN_OR_CBC):
+def load(network, solver=Solvers.COIN_OR_CBC):
     """
     Create a MIOM optimization model for a given solver.
     If the solver is Coin-OR CBC, an instance of PythonMipModel is used
@@ -157,12 +157,13 @@ def miom(network, solver=Solvers.COIN_OR_CBC):
 
     Example:
         Example of how to perform FBA to maximize flux through the
-        BIOMASS_reaction in the iMM1865 model:
+        `BIOMASS_reaction` in the iMM1865 model:
 
         ```python
-        >>> from miom import miom, load_gem
-        >>> network = load_gem("https://github.com/pablormier/miom-gems/raw/main/gems/mus_musculus_iMM1865.miom")
-        >>> V, _ = miom(network)
+        >>> import miom
+        >>> network = miom.mio.load_gem("https://github.com/pablormier/miom-gems/raw/main/gems/mus_musculus_iMM1865.miom")
+        >>> V, X = (miom
+                    .load(network)
                     .steady_state()
                     .set_rxn_objective("BIOMASS_reaction")
                     .solve(verbosity=1)
@@ -480,6 +481,9 @@ def keep(self, reactions):
             return False
         if isinstance(reactions, str):
             reactions = [reactions]
+        # Check if it's a binary vector
+        if len(reactions) == self.network.num_reactions and max(reactions)==1:
+            reactions = np.where(reactions==1)[0]
         # Check if the reactions have an indicator variable
         reactions = set(self.network.find_reaction(rxn)[0] for rxn in reactions)
         available = set(rxn.index for rxn in self.variables.assigned_reactions if rxn.cost > 0)
diff --git a/pyproject.toml b/pyproject.toml
index f5b4531..9979696 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "miom"
-version = "0.9.0-alpha.4"
+version = "0.9.0-beta.0"
 description = "Mixed Integer Optimization for Metabolism"
 authors = ["Pablo R. Mier <pablo.rodriguez.mier@gmail.com>"]
 keywords = ["optimization", "LP", "MIP", "metabolism", "metabolic-networks"]