diff --git a/dwave/optimization/model.pxd b/dwave/optimization/_model.pxd
similarity index 91%
rename from dwave/optimization/model.pxd
rename to dwave/optimization/_model.pxd
index 6a795f7c..57675d03 100644
--- a/dwave/optimization/model.pxd
+++ b/dwave/optimization/_model.pxd
@@ -26,7 +26,7 @@ from dwave.optimization.libcpp.state cimport State as cppState
 __all__ = ["Model"]
 
 
-cdef class _Model:
+cdef class _GraphManager:
     cpdef bool is_locked(self) noexcept
     cpdef Py_ssize_t num_nodes(self) noexcept
 
@@ -59,9 +59,9 @@ cdef class _Model:
     cdef object _data_sources
 
 
-cdef class Model(_Model):
-    cpdef Py_ssize_t num_constraints(self) noexcept
-    cpdef Py_ssize_t num_decisions(self) noexcept
+cdef class _Model(_GraphManager):
+    cpdef Py_ssize_t _num_constraints(self) noexcept
+    cpdef Py_ssize_t _num_decisions(self) noexcept
 
     cdef readonly object objective  # todo: cdef ArraySymbol?
     """Objective to be minimized.
@@ -94,7 +94,7 @@ cdef class States:
     cpdef resolve(self)
     cpdef Py_ssize_t size(self) except -1
 
-    cdef _Model _model(self)
+    cdef _GraphManager _model(self)
 
     # In order to not create a circular reference, we only hold a weakref
     # to the model from the states. This introduces some overhead, but it
@@ -117,7 +117,7 @@ cdef class States:
 
 cdef class Symbol:
     # Inheriting nodes must call this method from their __init__()
-    cdef void initialize_node(self, _Model model, cppNode* node_ptr) noexcept
+    cdef void initialize_node(self, _GraphManager model, cppNode* node_ptr) noexcept
 
     cpdef uintptr_t id(self) noexcept
 
@@ -125,12 +125,12 @@ cdef class Symbol:
     cpdef bool expired(self) noexcept
 
     @staticmethod
-    cdef Symbol from_ptr(_Model model, cppNode* ptr)
+    cdef Symbol from_ptr(_GraphManager model, cppNode* ptr)
 
     # Hold on to a reference to the Model, both for access but also, importantly,
     # to ensure that the model doesn't get garbage collected unless all of
     # the observers have also been garbage collected.
-    cdef readonly _Model model
+    cdef readonly _GraphManager model
 
     # Hold Node* pointer. This is redundant as most observers will also hold
     # a pointer to their observed node with the correct type. But the cost
@@ -148,7 +148,7 @@ cdef class Symbol:
 # also Symbols (probably a fair assumption)
 cdef class ArraySymbol(Symbol):
     # Inheriting symbols must call this method from their __init__()
-    cdef void initialize_arraynode(self, _Model model, cppArrayNode* array_ptr) noexcept
+    cdef void initialize_arraynode(self, _GraphManager model, cppArrayNode* array_ptr) noexcept
 
     # Hold ArrayNode* pointer. Again this is redundant, because we're also holding
     # a pointer to Node* and we can theoretically dynamic cast each time.
diff --git a/dwave/optimization/model.pyi b/dwave/optimization/_model.pyi
similarity index 75%
rename from dwave/optimization/model.pyi
rename to dwave/optimization/_model.pyi
index 8cf11dcd..c7e43ea4 100644
--- a/dwave/optimization/model.pyi
+++ b/dwave/optimization/_model.pyi
@@ -26,7 +26,7 @@ from dwave.optimization.symbols import *
 _ShapeLike: typing.TypeAlias = typing.Union[int, collections.abc.Sequence[int]]
 
 
-class _Model:
+class _GraphManager:
     @property
     def states(self) -> States: ...
 
@@ -71,50 +71,18 @@ class _Model:
     def unlock(self): ...
 
 
-class Model(_Model):
+class _Model(_GraphManager):
     def __init__(self): ...
 
     @property
     def objective(self) -> ArraySymbol: ...
-    def add_constraint(self, value: ArraySymbol) -> ArraySymbol: ...
-    def binary(self, shape: typing.Optional[_ShapeLike] = None) -> BinaryVariable: ...
-    def decision_state_size(self) -> int: ...
 
-    def disjoint_bit_sets(
-        self, primary_set_size: int, num_disjoint_sets: int,
-        ) -> tuple[DisjointBitSets, tuple[DisjointBitSet, ...]]: ...
-
-    def disjoint_lists(
-        self, primary_set_size: int, num_disjoint_lists: int,
-        ) -> tuple[DisjointLists, tuple[DisjointList, ...]]: ...
-
-    def feasible(self, index: int = 0) -> bool: ...
-
-    def integer(
-        self,
-        shape: typing.Optional[_ShapeLike] = None,
-        lower_bound: typing.Optional[int] = None,
-        upper_bound: typing.Optional[int] = None,
-        ) -> IntegerVariable: ...
-
-    def iter_constraints(self) -> collections.abc.Iterator[ArraySymbol]: ...
-    def iter_decisions(self) -> collections.abc.Iterator[Symbol]: ...
-    def list(self, n: int) -> ListVariable: ...
-    def lock(self) -> contextlib.AbstractContextManager: ...
-    def minimize(self, value: ArraySymbol): ...
-    def num_constraints(self) -> int: ...
-    def num_decisions(self) -> int: ...
-
-    # dev note: this is underspecified, but it would be quite complex to fully
-    # specify the linear/quadratic so let's leave it alone for now.
-    def quadratic_model(self, x: ArraySymbol, quadratic, linear=None) -> QuadraticModel: ...
-
-    def set(
-        self,
-        n: int,
-        min_size: int = 0,
-        max_size: typing.Optional[int] = None,
-        ) -> SetVariable: ...
+    def _add_constraint(self, value: ArraySymbol) -> ArraySymbol: ...
+    def _iter_constraints(self) -> collections.abc.Iterator[ArraySymbol]: ...
+    def _iter_decisions(self) -> collections.abc.Iterator[Symbol]: ...
+    def _minimize(self, value: ArraySymbol): ...
+    def _num_constraints(self) -> int: ...
+    def _num_decisions(self) -> int: ...
 
 
 class States:
diff --git a/dwave/optimization/model.pyx b/dwave/optimization/_model.pyx
similarity index 84%
rename from dwave/optimization/model.pyx
rename to dwave/optimization/_model.pyx
index 0f48f2a2..01d22d1c 100644
--- a/dwave/optimization/model.pyx
+++ b/dwave/optimization/_model.pyx
@@ -54,7 +54,7 @@ __all__ = ["Model"]
 
 
 ctypedef fused ExpressionOrModel:
-    Model
+    _Model
     Expression
 
 
@@ -67,7 +67,7 @@ def locked(model):
         model.unlock()
 
 
-cdef class _Model:
+cdef class _GraphManager:
     def constant(self, array_like):
         r"""Create a constant symbol.
 
@@ -128,7 +128,7 @@ cdef class _Model:
             header_len = struct.unpack('<I', file.read(4))[0]
             header_data = json.loads(file.read(header_len).decode('ascii'))
 
-        cdef Model model = cls()
+        cdef _Model model = cls()
 
         with zipfile.ZipFile(file, mode="r") as zf:
             model_info = json.loads(zf.read("info.json"))
@@ -634,241 +634,20 @@ cdef class _Model:
                 self.states._states[i].resize(self.num_decisions())
 
 
-cdef class Model(_Model):
-    """Nonlinear model.
-
-    The nonlinear model represents a general optimization problem with an
-    :term:`objective function` and/or constraints over variables of various
-    types.
-
-    The :class:`.Model` class can contain this model and its methods provide
-    convenient utilities for working with representations of a problem.
-
-    Examples:
-        This example creates a model for a
-        :class:`flow-shop-scheduling <dwave.optimization.generators.flow_shop_scheduling>`
-        problem with two jobs on three machines.
-
-        >>> from dwave.optimization.generators import flow_shop_scheduling
-        ...
-        >>> processing_times = [[10, 5, 7], [20, 10, 15]]
-        >>> model = flow_shop_scheduling(processing_times=processing_times)
-    """
+cdef class _Model(_GraphManager):
     def __init__(self):
         self.states = States(self)
 
         self._data_sources = []
 
-    def add_constraint(self, ArraySymbol value):
-        """Add a constraint to the model.
-
-        Args:
-            value: Value that must evaluate to True for the state
-                of the model to be feasible.
-
-        Returns:
-            The constraint symbol.
-
-        Examples:
-            This example adds a single constraint to a model.
-
-            >>> from dwave.optimization.model import Model
-            >>> model = Model()
-            >>> i = model.integer()
-            >>> c = model.constant(5)
-            >>> constraint_sym = model.add_constraint(i <= c)
-
-            The returned constraint symbol can be assigned and evaluated
-            for a model state:
-
-            >>> with model.lock():
-            ...     model.states.resize(1)
-            ...     i.set_state(0, 1) # Feasible state
-            ...     print(constraint_sym.state(0))
-            1.0
-            >>> with model.lock():
-            ...     i.set_state(0, 6) # Infeasible state
-            ...     print(constraint_sym.state(0))
-            0.0
-        """
+    def _add_constraint(self, ArraySymbol value):
         if value is None:
             raise ValueError("value cannot be None")
         # TODO: shall we accept array valued constraints?
         self._graph.add_constraint(value.array_ptr)
         return value
 
-    def binary(self, shape=None):
-        r"""Create a binary symbol as a decision variable.
-
-        Args:
-            shape: Shape of the binary array to create.
-
-        Returns:
-            A binary symbol.
-
-        Examples:
-            This example creates a :math:`1 \times 20`-sized binary symbol.
-
-            >>> from dwave.optimization.model import Model
-            >>> model = Model()
-            >>> x = model.binary((1,20))
-        """
-        from dwave.optimization.symbols import BinaryVariable #avoid circular import
-        return BinaryVariable(self, shape)
-
-    def decision_state_size(self):
-        r"""An estimated size, in bytes, of the model's decision states.
-
-        Examples:
-            This example checks the size of a model with one
-            :math:`10 \times 10`-sized integer symbol.
-
-            >>> from dwave.optimization.model import Model
-            >>> model = Model()
-            >>> visit_site = model.integer((10, 10))
-            >>> model.decision_state_size()
-            800
-        """
-        return sum(sym.state_size() for sym in self.iter_decisions())
-
-    def disjoint_bit_sets(self, Py_ssize_t primary_set_size, Py_ssize_t num_disjoint_sets):
-        """Create a disjoint-sets symbol as a decision variable.
-
-        Divides a set of the elements of ``range(primary_set_size)`` into
-        ``num_disjoint_sets`` ordered partitions, stored as bit sets (arrays
-        of length ``primary_set_size``, with ones at the indices of elements
-        currently in the set, and zeros elsewhere). The ordering of a set is
-        not semantically meaningful.
-
-        Also creates from the symbol ``num_disjoint_sets`` extra successors
-        that output the disjoint sets as arrays.
-
-        Args:
-            primary_set_size: Number of elements in the primary set that are
-                partitioned into disjoint sets. Must be non-negative.
-            num_disjoint_sets: Number of disjoint sets. Must be positive.
-
-        Returns:
-            A tuple where the first element is the disjoint-sets symbol and
-            the second is a set of its newly added successors.
-
-        Examples:
-            This example creates a symbol of 10 elements that is divided
-            into 4 sets.
-
-            >>> from dwave.optimization.model import Model
-            >>> model = Model()
-            >>> parts_set, parts_subsets = model.disjoint_bit_sets(10, 4)
-        """
-
-        from dwave.optimization.symbols import DisjointBitSets, DisjointBitSet  # avoid circular import
-        main = DisjointBitSets(self, primary_set_size, num_disjoint_sets)
-        sets = tuple(DisjointBitSet(main, i) for i in range(num_disjoint_sets))
-        return main, sets
-
-    def disjoint_lists(self, Py_ssize_t primary_set_size, Py_ssize_t num_disjoint_lists):
-        """Create a disjoint-lists symbol as a decision variable.
-
-        Divides a set of the elements of ``range(primary_set_size)`` into
-        ``num_disjoint_lists`` ordered partitions.
-
-        Also creates ``num_disjoint_lists`` extra successors from the
-        symbol that output the disjoint lists as arrays.
-
-        Args:
-            primary_set_size: Number of elements in the primary set to
-                be partitioned into disjoint lists.
-            num_disjoint_lists: Number of disjoint lists.
-
-        Returns:
-            A tuple where the first element is the disjoint-lists symbol
-            and the second is a list of its newly added successor nodes.
-
-        Examples:
-            This example creates a symbol of 10 elements that is divided
-            into 4 lists.
-
-            >>> from dwave.optimization.model import Model
-            >>> model = Model()
-            >>> destinations, routes = model.disjoint_lists(10, 4)
-        """
-        from dwave.optimization.symbols import DisjointLists, DisjointList  # avoid circular import
-        main = DisjointLists(self, primary_set_size, num_disjoint_lists)
-        lists = [DisjointList(main, i) for i in range(num_disjoint_lists)]
-        return main, lists
-
-    def feasible(self, int index = 0):
-        """Check the feasibility of the state at the input index.
-
-        Args:
-            index: index of the state to check for feasibility.
-
-        Returns:
-            Feasibility of the state.
-
-        Examples:
-            This example demonstrates checking the feasibility of a simple model with
-            feasible and infeasible states.
-
-            >>> from dwave.optimization.model import Model
-            >>> model = Model()
-            >>> b = model.binary()
-            >>> model.add_constraint(b) # doctest: +ELLIPSIS
-            <dwave.optimization.BinaryVariable at ...>
-            >>> model.states.resize(2)
-            >>> b.set_state(0, 1) # Feasible
-            >>> b.set_state(1, 0) # Infeasible
-            >>> with model.lock():
-            ...     model.feasible(0)
-            True
-            >>> with model.lock():
-            ...     model.feasible(1)
-            False
-        """
-        return all(sym.state(index) for sym in self.iter_constraints())
-
-    def integer(self, shape=None, lower_bound=None, upper_bound=None):
-        r"""Create an integer symbol as a decision variable.
-
-        Args:
-            shape: Shape of the integer array to create.
-
-            lower_bound: Lower bound for the symbol, which is the
-                smallest allowed integer value. If None, the default
-                value is used.
-            upper_bound: Upper bound for the symbol, which is the
-                largest allowed integer value. If None, the default
-                value is used.
-
-        Returns:
-            An integer symbol.
-
-        Examples:
-            This example creates a :math:`20 \times 20`-sized integer symbol.
-
-            >>> from dwave.optimization.model import Model
-            >>> model = Model()
-            >>> i = model.integer((20,20), lower_bound=-100, upper_bound=100)
-        """
-        from dwave.optimization.symbols import IntegerVariable #avoid circular import
-        return IntegerVariable(self, shape, lower_bound, upper_bound)
-
-    def _header_data(self, *, only_decision, max_num_states=float('inf')):
-        """The header data associated with the model (but not the states)."""
-        num_nodes = self.num_decisions() if only_decision else self.num_nodes()
-        num_states = max(0, min(self.states.size(), max_num_states))
-
-        decision_state_size = self.decision_state_size()
-        state_size = decision_state_size if only_decision else self.state_size()
-
-        return dict(
-            decision_state_size=decision_state_size,
-            num_nodes=num_nodes,
-            state_size=state_size,
-            num_states=num_states,
-        )
-
-    def iter_constraints(self):
+    def _iter_constraints(self):
         """Iterate over all constraints in the model.
 
         Examples:
@@ -885,7 +664,7 @@ cdef class Model(_Model):
         for ptr in self._graph.constraints():
             yield symbol_from_ptr(self, ptr)
 
-    def iter_decisions(self):
+    def _iter_decisions(self):
         """Iterate over all decision variables in the model.
 
         Examples:
@@ -902,48 +681,7 @@ cdef class Model(_Model):
         for ptr in self._graph.decisions():
             yield symbol_from_ptr(self, ptr)
 
-    def list(self, n : int):
-        """Create a list symbol as a decision variable.
-
-        Args:
-            n: Values in the list are permutations of ``range(n)``.
-
-        Returns:
-            A list symbol.
-
-        Examples:
-            This example creates a list symbol of 200 elements.
-
-            >>> from dwave.optimization.model import Model
-            >>> model = Model()
-            >>> routes = model.list(200)
-        """
-        from dwave.optimization.symbols import ListVariable  # avoid circular import
-        return ListVariable(self, n)
-
-    def minimize(self, ArraySymbol value):
-        """Set the objective value to minimize.
-
-        Optimization problems have an objective and/or constraints. The objective
-        expresses one or more aspects of the problem that should be minimized
-        (equivalent to maximization when multiplied by a minus sign). For example,
-        an optimized itinerary might minimize the value of distance traveled or
-        cost of transportation or travel time.
-
-        Args:
-            value: Value for which to minimize the cost function.
-
-        Examples:
-            This example minimizes a simple polynomial, :math:`y = i^2 - 4i`,
-            within bounds.
-
-            >>> from dwave.optimization import Model
-            >>> model = Model()
-            >>> i = model.integer(lower_bound=-5, upper_bound=5)
-            >>> c = model.constant(4)
-            >>> y = i*i - c*i
-            >>> model.minimize(y)
-        """
+    def _minimize(self, ArraySymbol value):
         if value is None:
             raise ValueError("value cannot be None")
         if value.size() < 1:
@@ -953,7 +691,7 @@ cdef class Model(_Model):
         self._graph.set_objective(value.array_ptr)
         self.objective = value
 
-    cpdef Py_ssize_t num_constraints(self) noexcept:
+    cpdef Py_ssize_t _num_constraints(self) noexcept:
         """Number of constraints in the model.
 
         Examples:
@@ -973,7 +711,7 @@ cdef class Model(_Model):
         """
         return self._graph.num_constraints()
 
-    cpdef Py_ssize_t num_decisions(self) noexcept:
+    cpdef Py_ssize_t _num_decisions(self) noexcept:
         """Number of independent decision nodes in the model.
 
         An array-of-integers symbol, for example, counts as a single
@@ -992,54 +730,6 @@ cdef class Model(_Model):
         """
         return self._graph.num_decisions()
 
-    def quadratic_model(self, ArraySymbol x, quadratic, linear=None):
-        """Create a quadratic model from an array and a quadratic model.
-
-        Args:
-            x: An array.
-
-            quadratic: Quadratic values for the quadratic model.
-
-            linear: Linear values for the quadratic model.
-
-        Returns:
-            A quadratic model.
-
-        Examples:
-            This example creates a quadratic model.
-
-            >>> from dwave.optimization.model import Model
-            >>> model = Model()
-            >>> x = model.binary(3)
-            >>> Q = {(0, 0): 0, (0, 1): 1, (0, 2): 2, (1, 1): 1, (1, 2): 3, (2, 2): 2}
-            >>> qm = model.quadratic_model(x, Q)
-
-        """
-        from dwave.optimization.symbols import QuadraticModel
-        return QuadraticModel(x, quadratic, linear)
-
-    def set(self, Py_ssize_t n, Py_ssize_t min_size = 0, max_size = None):
-        """Create a set symbol as a decision variable.
-
-        Args:
-            n: Values in the set are subsets of ``range(n)``.
-            min_size: Minimum set size. Defaults to ``0``.
-            max_size: Maximum set size. Defaults to ``n``.
-
-        Returns:
-            A set symbol.
-
-        Examples:
-            This example creates a set symbol of up to 4 elements
-            with values between 0 to 99.
-
-            >>> from dwave.optimization.model import Model
-            >>> model = Model()
-            >>> destinations = model.set(100, max_size=4)
-        """
-        from dwave.optimization.symbols import SetVariable  # avoid circular import
-        return SetVariable(self, n, min_size, n if max_size is None else max_size)
-
 
 def _States_init(States self, ExpressionOrModel model):
     self._model_ref = weakref.ref(model)
@@ -1177,6 +867,10 @@ cdef class States:
         Returns:
             A model.
         """
+
+        # avoid circular import
+        from dwave.optimization.model import Model
+
         self.resolve()
 
         if not replace:
@@ -1184,7 +878,7 @@ cdef class States:
 
         # todo: we don't need to actually construct a model, but this is nice and
         # abstract. We should performance test and then potentially re-implement
-        cdef Model model = Model.from_file(file, check_header=check_header)
+        cdef _Model model = Model.from_file(file, check_header=check_header)
 
         # Check that the model is compatible
         for n0, n1 in zip(model.iter_symbols(), self._model().iter_symbols()):
@@ -1215,7 +909,7 @@ cdef class States:
         """Initialize any uninitialized states."""
         self.resolve()
 
-        cdef Model model = self._model()
+        cdef _Model model = self._model()
 
         if not model.is_locked():
             raise ValueError("Cannot initialize states of an unlocked model")
@@ -1238,9 +932,9 @@ cdef class States:
         return self._model().into_file(file, only_decision=True, max_num_states=self.size())
 
 
-    cdef _Model _model(self):
+    cdef _GraphManager _model(self):
         """Get a ref-counted Model object."""
-        cdef Model m = self._model_ref()
+        cdef _Model m = self._model_ref()
         if m is None:
             raise ReferenceError("accessing the states of a garbage collected model")
         return m
@@ -1334,7 +1028,7 @@ cdef class Symbol:
         cls = type(self)
         return f"<{cls.__module__}.{cls.__qualname__} at {self.id():#x}>"
 
-    cdef void initialize_node(self, _Model model, cppNode* node_ptr) noexcept:
+    cdef void initialize_node(self, _GraphManager model, cppNode* node_ptr) noexcept:
         self.model = model
 
         self.node_ptr = node_ptr
@@ -1365,7 +1059,7 @@ cdef class Symbol:
         return deref(self.expired_ptr)
 
     @staticmethod
-    cdef Symbol from_ptr(Model model, cppNode* ptr):
+    cdef Symbol from_ptr(_Model model, cppNode* ptr):
         """Construct a Symbol from a C++ Node pointer.
 
         There are times when a Node* needs to be passed through the Python layer
@@ -1386,7 +1080,7 @@ cdef class Symbol:
         raise ValueError("Symbols cannot be constructed directly")
 
     @classmethod
-    def _from_zipfile(cls, zf, directory, Model model, predecessors):
+    def _from_zipfile(cls, zf, directory, _Model model, predecessors):
         """Construct a symbol from a compressed file.
 
         Args:
@@ -1766,7 +1460,7 @@ cdef class ArraySymbol(Symbol):
         # via their subclasses.
         raise ValueError("ArraySymbols cannot be constructed directly")
 
-    cdef void initialize_arraynode(self, _Model model, cppArrayNode* array_ptr) noexcept:
+    cdef void initialize_arraynode(self, _GraphManager model, cppArrayNode* array_ptr) noexcept:
         self.array_ptr = array_ptr
         self.initialize_node(model, array_ptr)
 
diff --git a/dwave/optimization/expression/expression.pxd b/dwave/optimization/expression/expression.pxd
index 534437a5..ff5abbde 100644
--- a/dwave/optimization/expression/expression.pxd
+++ b/dwave/optimization/expression/expression.pxd
@@ -14,13 +14,13 @@
 #    See the License for the specific language governing permissions and
 #    limitations under the License.
 
-from dwave.optimization.model cimport ArraySymbol, _Model
+from dwave.optimization._model cimport ArraySymbol, _GraphManager
 
 
 __all__ = ["Expression"]
 
 
-cdef class Expression(_Model):
+cdef class Expression(_GraphManager):
     cdef readonly ArraySymbol output
 
     cpdef Py_ssize_t num_inputs(self) noexcept
diff --git a/dwave/optimization/expression/expression.pyi b/dwave/optimization/expression/expression.pyi
index 5ad7a07e..edf836ea 100644
--- a/dwave/optimization/expression/expression.pyi
+++ b/dwave/optimization/expression/expression.pyi
@@ -14,13 +14,13 @@
 
 from typing import Optional
 
-from dwave.optimization.model import _Model
+from dwave.optimization._model import _GraphManager
 
 
 _ShapeLike: typing.TypeAlias = typing.Union[int, collections.abc.Sequence[int]]
 
 
-class Expression(_Model):
+class Expression(_GraphManager):
     def __init__(
         self,
         num_inputs: int = 0,
diff --git a/dwave/optimization/expression/expression.pyx b/dwave/optimization/expression/expression.pyx
index 0ebae6ab..4210dff7 100644
--- a/dwave/optimization/expression/expression.pyx
+++ b/dwave/optimization/expression/expression.pyx
@@ -21,7 +21,7 @@ from libcpp.cast cimport dynamic_cast
 from dwave.optimization.libcpp.array cimport Array as cppArray
 from dwave.optimization.libcpp.graph cimport Node as cppNode
 from dwave.optimization.libcpp.nodes cimport InputNode as cppInputNode
-from dwave.optimization.model cimport ArraySymbol, _Model, States
+from dwave.optimization._model cimport ArraySymbol, _GraphManager, States
 from dwave.optimization.symbols cimport symbol_from_ptr
 
 ctypedef cppNode* cppNodePtr
@@ -30,7 +30,7 @@ ctypedef cppNode* cppNodePtr
 __all__ = ["Expression"]
 
 
-cdef class Expression(_Model):
+cdef class Expression(_GraphManager):
     def __init__(
         self,
         num_inputs: int = 0,
diff --git a/dwave/optimization/model.py b/dwave/optimization/model.py
new file mode 100644
index 00000000..6ec4a665
--- /dev/null
+++ b/dwave/optimization/model.py
@@ -0,0 +1,395 @@
+from typing import Optional
+from dwave.optimization._model import _Model, ArraySymbol, States
+
+
+class Model(_Model):
+    """Nonlinear model.
+
+    The nonlinear model represents a general optimization problem with an
+    :term:`objective function` and/or constraints over variables of various
+    types.
+
+    The :class:`.Model` class can contain this model and its methods provide
+    convenient utilities for working with representations of a problem.
+
+    Examples:
+        This example creates a model for a
+        :class:`flow-shop-scheduling <dwave.optimization.generators.flow_shop_scheduling>`
+        problem with two jobs on three machines.
+
+        >>> from dwave.optimization.generators import flow_shop_scheduling
+        ...
+        >>> processing_times = [[10, 5, 7], [20, 10, 15]]
+        >>> model = flow_shop_scheduling(processing_times=processing_times)
+    """
+
+    def __init__(self):
+        super().__init__()
+
+    def add_constraint(self, value: ArraySymbol):
+        """Add a constraint to the model.
+
+        Args:
+            value: Value that must evaluate to True for the state
+                of the model to be feasible.
+
+        Returns:
+            The constraint symbol.
+
+        Examples:
+            This example adds a single constraint to a model.
+
+            >>> from dwave.optimization.model import Model
+            >>> model = Model()
+            >>> i = model.integer()
+            >>> c = model.constant(5)
+            >>> constraint_sym = model.add_constraint(i <= c)
+
+            The returned constraint symbol can be assigned and evaluated
+            for a model state:
+
+            >>> with model.lock():
+            ...     model.states.resize(1)
+            ...     i.set_state(0, 1) # Feasible state
+            ...     print(constraint_sym.state(0))
+            1.0
+            >>> with model.lock():
+            ...     i.set_state(0, 6) # Infeasible state
+            ...     print(constraint_sym.state(0))
+            0.0
+        """
+        return super()._add_constraint(value)
+
+    def binary(self, shape=None):
+        r"""Create a binary symbol as a decision variable.
+
+        Args:
+            shape: Shape of the binary array to create.
+
+        Returns:
+            A binary symbol.
+
+        Examples:
+            This example creates a :math:`1 \times 20`-sized binary symbol.
+
+            >>> from dwave.optimization.model import Model
+            >>> model = Model()
+            >>> x = model.binary((1,20))
+        """
+        from dwave.optimization.symbols import BinaryVariable #avoid circular import
+        return BinaryVariable(self, shape)
+
+    def decision_state_size(self):
+        r"""An estimated size, in bytes, of the model's decision states.
+
+        Examples:
+            This example checks the size of a model with one
+            :math:`10 \times 10`-sized integer symbol.
+
+            >>> from dwave.optimization.model import Model
+            >>> model = Model()
+            >>> visit_site = model.integer((10, 10))
+            >>> model.decision_state_size()
+            800
+        """
+        return sum(sym.state_size() for sym in self.iter_decisions())
+
+    def disjoint_bit_sets(self, primary_set_size: int, num_disjoint_sets: int):
+        """Create a disjoint-sets symbol as a decision variable.
+
+        Divides a set of the elements of ``range(primary_set_size)`` into
+        ``num_disjoint_sets`` ordered partitions, stored as bit sets (arrays
+        of length ``primary_set_size``, with ones at the indices of elements
+        currently in the set, and zeros elsewhere). The ordering of a set is
+        not semantically meaningful.
+
+        Also creates from the symbol ``num_disjoint_sets`` extra successors
+        that output the disjoint sets as arrays.
+
+        Args:
+            primary_set_size: Number of elements in the primary set that are
+                partitioned into disjoint sets. Must be non-negative.
+            num_disjoint_sets: Number of disjoint sets. Must be positive.
+
+        Returns:
+            A tuple where the first element is the disjoint-sets symbol and
+            the second is a set of its newly added successors.
+
+        Examples:
+            This example creates a symbol of 10 elements that is divided
+            into 4 sets.
+
+            >>> from dwave.optimization.model import Model
+            >>> model = Model()
+            >>> parts_set, parts_subsets = model.disjoint_bit_sets(10, 4)
+        """
+
+        from dwave.optimization.symbols import DisjointBitSets, DisjointBitSet  # avoid circular import
+        main = DisjointBitSets(self, primary_set_size, num_disjoint_sets)
+        sets = tuple(DisjointBitSet(main, i) for i in range(num_disjoint_sets))
+        return main, sets
+
+    def disjoint_lists(self, primary_set_size: int, num_disjoint_lists: int):
+        """Create a disjoint-lists symbol as a decision variable.
+
+        Divides a set of the elements of ``range(primary_set_size)`` into
+        ``num_disjoint_lists`` ordered partitions.
+
+        Also creates ``num_disjoint_lists`` extra successors from the
+        symbol that output the disjoint lists as arrays.
+
+        Args:
+            primary_set_size: Number of elements in the primary set to
+                be partitioned into disjoint lists.
+            num_disjoint_lists: Number of disjoint lists.
+
+        Returns:
+            A tuple where the first element is the disjoint-lists symbol
+            and the second is a list of its newly added successor nodes.
+
+        Examples:
+            This example creates a symbol of 10 elements that is divided
+            into 4 lists.
+
+            >>> from dwave.optimization.model import Model
+            >>> model = Model()
+            >>> destinations, routes = model.disjoint_lists(10, 4)
+        """
+        from dwave.optimization.symbols import DisjointLists, DisjointList  # avoid circular import
+        main = DisjointLists(self, primary_set_size, num_disjoint_lists)
+        lists = [DisjointList(main, i) for i in range(num_disjoint_lists)]
+        return main, lists
+
+    def feasible(self, index: int = 0):
+        """Check the feasibility of the state at the input index.
+
+        Args:
+            index: index of the state to check for feasibility.
+
+        Returns:
+            Feasibility of the state.
+
+        Examples:
+            This example demonstrates checking the feasibility of a simple model with
+            feasible and infeasible states.
+
+            >>> from dwave.optimization.model import Model
+            >>> model = Model()
+            >>> b = model.binary()
+            >>> model.add_constraint(b) # doctest: +ELLIPSIS
+            <dwave.optimization.BinaryVariable at ...>
+            >>> model.states.resize(2)
+            >>> b.set_state(0, 1) # Feasible
+            >>> b.set_state(1, 0) # Infeasible
+            >>> with model.lock():
+            ...     model.feasible(0)
+            True
+            >>> with model.lock():
+            ...     model.feasible(1)
+            False
+        """
+        return all(sym.state(index) for sym in self.iter_constraints())
+
+    def integer(self, shape=None, lower_bound=None, upper_bound=None):
+        r"""Create an integer symbol as a decision variable.
+
+        Args:
+            shape: Shape of the integer array to create.
+
+            lower_bound: Lower bound for the symbol, which is the
+                smallest allowed integer value. If None, the default
+                value is used.
+            upper_bound: Upper bound for the symbol, which is the
+                largest allowed integer value. If None, the default
+                value is used.
+
+        Returns:
+            An integer symbol.
+
+        Examples:
+            This example creates a :math:`20 \times 20`-sized integer symbol.
+
+            >>> from dwave.optimization.model import Model
+            >>> model = Model()
+            >>> i = model.integer((20,20), lower_bound=-100, upper_bound=100)
+        """
+        from dwave.optimization.symbols import IntegerVariable #avoid circular import
+        return IntegerVariable(self, shape, lower_bound, upper_bound)
+
+    def _header_data(self, *, only_decision, max_num_states=float('inf')):
+        """The header data associated with the model (but not the states)."""
+        num_nodes = self.num_decisions() if only_decision else self.num_nodes()
+        num_states = max(0, min(self.states.size(), max_num_states))
+
+        decision_state_size = self.decision_state_size()
+        state_size = decision_state_size if only_decision else self.state_size()
+
+        return dict(
+            decision_state_size=decision_state_size,
+            num_nodes=num_nodes,
+            state_size=state_size,
+            num_states=num_states,
+        )
+
+    def iter_constraints(self):
+        """Iterate over all constraints in the model.
+
+        Examples:
+            This example adds a single constraint to a model and iterates over it.
+
+            >>> from dwave.optimization.model import Model
+            >>> model = Model()
+            >>> i = model.integer()
+            >>> c = model.constant(5)
+            >>> model.add_constraint(i <= c) # doctest: +ELLIPSIS
+            <dwave.optimization.symbols.LessEqual at ...>
+            >>> constraints = next(model.iter_constraints())
+        """
+        return super()._iter_constraints()
+
+    def iter_decisions(self):
+        """Iterate over all decision variables in the model.
+
+        Examples:
+            This example adds a single decision symbol to a model and iterates over it.
+
+            >>> from dwave.optimization.model import Model
+            >>> model = Model()
+            >>> i = model.integer()
+            >>> c = model.constant(5)
+            >>> model.add_constraint(i <= c) # doctest: +ELLIPSIS
+            <dwave.optimization.symbols.LessEqual at ...>
+            >>> decisions = next(model.iter_decisions())
+        """
+        return super()._iter_decisions()
+
+    def list(self, n: int):
+        """Create a list symbol as a decision variable.
+
+        Args:
+            n: Values in the list are permutations of ``range(n)``.
+
+        Returns:
+            A list symbol.
+
+        Examples:
+            This example creates a list symbol of 200 elements.
+
+            >>> from dwave.optimization.model import Model
+            >>> model = Model()
+            >>> routes = model.list(200)
+        """
+        from dwave.optimization.symbols import ListVariable  # avoid circular import
+        return ListVariable(self, n)
+
+    def minimize(self, value: ArraySymbol):
+        """Set the objective value to minimize.
+
+        Optimization problems have an objective and/or constraints. The objective
+        expresses one or more aspects of the problem that should be minimized
+        (equivalent to maximization when multiplied by a minus sign). For example,
+        an optimized itinerary might minimize the value of distance traveled or
+        cost of transportation or travel time.
+
+        Args:
+            value: Value for which to minimize the cost function.
+
+        Examples:
+            This example minimizes a simple polynomial, :math:`y = i^2 - 4i`,
+            within bounds.
+
+            >>> from dwave.optimization import Model
+            >>> model = Model()
+            >>> i = model.integer(lower_bound=-5, upper_bound=5)
+            >>> c = model.constant(4)
+            >>> y = i*i - c*i
+            >>> model.minimize(y)
+        """
+        super()._minimize(value)
+
+    def num_constraints(self):
+        """Number of constraints in the model.
+
+        Examples:
+            This example checks the number of constraints in the model after
+            adding a couple of constraints.
+
+            >>> from dwave.optimization.model import Model
+            >>> model = Model()
+            >>> i = model.integer()
+            >>> c = model.constant([5, -14])
+            >>> model.add_constraint(i <= c[0]) # doctest: +ELLIPSIS
+            <dwave.optimization.symbols.LessEqual at ...>
+            >>> model.add_constraint(c[1] <= i) # doctest: +ELLIPSIS
+            <dwave.optimization.symbols.LessEqual at ...>
+            >>> model.num_constraints()
+            2
+        """
+        return super()._num_constraints()
+
+    def num_decisions(self):
+        """Number of independent decision nodes in the model.
+
+        An array-of-integers symbol, for example, counts as a single
+        decision node.
+
+        Examples:
+            This example checks the number of decisions in a model after
+            adding a single (size 20) decision symbol.
+
+            >>> from dwave.optimization.model import Model
+            >>> model = Model()
+            >>> c = model.constant([1, 5, 8.4])
+            >>> i = model.integer(20, upper_bound=100)
+            >>> model.num_decisions()
+            1
+        """
+        return super()._num_decisions()
+
+    def quadratic_model(self, x: ArraySymbol, quadratic, linear=None):
+        """Create a quadratic model from an array and a quadratic model.
+
+        Args:
+            x: An array.
+
+            quadratic: Quadratic values for the quadratic model.
+
+            linear: Linear values for the quadratic model.
+
+        Returns:
+            A quadratic model.
+
+        Examples:
+            This example creates a quadratic model.
+
+            >>> from dwave.optimization.model import Model
+            >>> model = Model()
+            >>> x = model.binary(3)
+            >>> Q = {(0, 0): 0, (0, 1): 1, (0, 2): 2, (1, 1): 1, (1, 2): 3, (2, 2): 2}
+            >>> qm = model.quadratic_model(x, Q)
+
+        """
+        from dwave.optimization.symbols import QuadraticModel
+        return QuadraticModel(x, quadratic, linear)
+
+    def set(self, n: int, min_size: int = 0, max_size: Optional[int] = None):
+        """Create a set symbol as a decision variable.
+
+        Args:
+            n: Values in the set are subsets of ``range(n)``.
+            min_size: Minimum set size. Defaults to ``0``.
+            max_size: Maximum set size. Defaults to ``n``.
+
+        Returns:
+            A set symbol.
+
+        Examples:
+            This example creates a set symbol of up to 4 elements
+            with values between 0 to 99.
+
+            >>> from dwave.optimization.model import Model
+            >>> model = Model()
+            >>> destinations = model.set(100, max_size=4)
+        """
+        from dwave.optimization.symbols import SetVariable  # avoid circular import
+        return SetVariable(self, n, min_size, n if max_size is None else max_size)
diff --git a/dwave/optimization/symbols.pxd b/dwave/optimization/symbols.pxd
index 64c78e87..b2bbbe43 100644
--- a/dwave/optimization/symbols.pxd
+++ b/dwave/optimization/symbols.pxd
@@ -16,10 +16,10 @@
 
 from libcpp.typeinfo cimport type_info
 
-from dwave.optimization.model cimport _Model
+from dwave.optimization._model cimport _GraphManager
 from dwave.optimization.libcpp.graph cimport Array as cppArray
 from dwave.optimization.libcpp.graph cimport Node as cppNode
 
 cdef void _register(object cls, const type_info& typeinfo)
 
-cdef object symbol_from_ptr(_Model model, cppNode* ptr)
+cdef object symbol_from_ptr(_GraphManager model, cppNode* ptr)
diff --git a/dwave/optimization/symbols.pyi b/dwave/optimization/symbols.pyi
index b42656f7..00d8480e 100644
--- a/dwave/optimization/symbols.pyi
+++ b/dwave/optimization/symbols.pyi
@@ -16,7 +16,7 @@ import typing
 
 import numpy.typing
 
-from dwave.optimization.model import Symbol, ArraySymbol
+from dwave.optimization._model import Symbol, ArraySymbol
 
 
 class Absolute(ArraySymbol):
diff --git a/dwave/optimization/symbols.pyx b/dwave/optimization/symbols.pyx
index c07a95ec..5fae3741 100644
--- a/dwave/optimization/symbols.pyx
+++ b/dwave/optimization/symbols.pyx
@@ -99,11 +99,11 @@ from dwave.optimization.libcpp.nodes cimport (
     WhereNode as cppWhereNode,
     XorNode as cppXorNode,
     )
-from dwave.optimization.model cimport ArraySymbol, _Model, Model, Symbol
+from dwave.optimization._model cimport ArraySymbol, _GraphManager, _Model, Symbol
 from dwave.optimization.expression cimport Expression
 
 ctypedef fused ExpressionOrModel:
-    Model
+    _Model
     Expression
 
 
@@ -185,7 +185,7 @@ cdef void _register(object cls, const type_info& typeinfo):
 
 
 # TODO: should this use ExpressionOrModel?
-cdef object symbol_from_ptr(_Model model, cppNode* node_ptr):
+cdef object symbol_from_ptr(_GraphManager model, cppNode* node_ptr):
     """Create a Python/Cython symbol from a C++ Node*."""
 
     # If it's null, either after the cast of just as given, then we can't get a symbol from it
@@ -245,7 +245,7 @@ cdef class Absolute(ArraySymbol):
         <class 'dwave.optimization.symbols.Absolute'>
     """
     def __init__(self, ArraySymbol x):
-        cdef Model model = x.model
+        cdef _Model model = x.model
 
         self.ptr = model._graph.emplace_node[cppAbsoluteNode](x.array_ptr)
         self.initialize_arraynode(model, self.ptr)
@@ -283,7 +283,7 @@ cdef class Add(ArraySymbol):
         if lhs.model is not rhs.model:
             raise ValueError("lhs and rhs do not share the same underlying model")
 
-        cdef _Model model = lhs.model
+        cdef _GraphManager model = lhs.model
 
         self.ptr = model._graph.emplace_node[cppAddNode](lhs.array_ptr, rhs.array_ptr)
         self.initialize_arraynode(model, self.ptr)
@@ -317,7 +317,7 @@ cdef class All(ArraySymbol):
         <class 'dwave.optimization.symbols.All'>
     """
     def __init__(self, ArraySymbol array):
-        cdef Model model = array.model
+        cdef _Model model = array.model
         self.ptr = model._graph.emplace_node[cppAllNode](array.array_ptr)
         self.initialize_arraynode(model, self.ptr)
 
@@ -346,7 +346,7 @@ cdef class And(ArraySymbol):
         if lhs.model is not rhs.model:
             raise ValueError("lhs and rhs do not share the same underlying model")
 
-        cdef Model model = lhs.model
+        cdef _Model model = lhs.model
 
         self.ptr = model._graph.emplace_node[cppAndNode](lhs.array_ptr, rhs.array_ptr)
         self.initialize_arraynode(model, self.ptr)
@@ -388,7 +388,7 @@ cdef class Any(ArraySymbol):
     .. versionadded:: 0.4.1
     """
     def __init__(self, ArraySymbol array):
-        cdef Model model = array.model
+        cdef _Model model = array.model
         self.ptr = model._graph.emplace_node[cppAnyNode](array.array_ptr)
         self.initialize_arraynode(model, self.ptr)
 
@@ -409,7 +409,7 @@ _register(Any, typeid(cppAnyNode))
 
 cdef class _ArrayValidation(Symbol):
     def __init__(self, ArraySymbol array_node):
-        cdef Model model = array_node.model
+        cdef _Model model = array_node.model
 
         self.ptr = model._graph.emplace_node[cppArrayValidationNode](array_node.array_ptr)
 
@@ -446,7 +446,7 @@ cdef class AdvancedIndexing(ArraySymbol):
         <class 'dwave.optimization.symbols.AdvancedIndexing'>
     """
     def __init__(self, ArraySymbol array, *indices):
-        cdef Model model = array.model
+        cdef _Model model = array.model
 
         cdef vector[cppAdvancedIndexingNode.array_or_slice] cppindices
 
@@ -517,7 +517,7 @@ cdef class AdvancedIndexing(ArraySymbol):
         return sym
 
     @classmethod
-    def _from_zipfile(cls, zf, directory, Model model, predecessors):
+    def _from_zipfile(cls, zf, directory, _Model model, predecessors):
         cdef cppNode* ptr
 
         indices = []
@@ -575,7 +575,7 @@ cdef class BasicIndexing(ArraySymbol):
     """
     def __init__(self, ArraySymbol array, *indices):
 
-        cdef Model model = array.model
+        cdef _Model model = array.model
 
         cdef vector[cppBasicIndexingNode.slice_or_int] cppindices
         for index in indices:
@@ -616,7 +616,7 @@ cdef class BasicIndexing(ArraySymbol):
         return sym
 
     @classmethod
-    def _from_zipfile(cls, zf, directory, Model model, predecessors):
+    def _from_zipfile(cls, zf, directory, _Model model, predecessors):
         if len(predecessors) != 1:
             raise ValueError(f"`BasicIndexing` should have exactly one predecessor")
 
@@ -675,7 +675,7 @@ cdef class BinaryVariable(ArraySymbol):
         >>> type(x)
         <class 'dwave.optimization.symbols.BinaryVariable'>
     """
-    def __init__(self, Model model, shape=None):
+    def __init__(self, _Model model, shape=None):
         # Get an observing pointer to the node
         cdef vector[Py_ssize_t] vshape = _as_cppshape(tuple() if shape is None else shape)
 
@@ -695,7 +695,7 @@ cdef class BinaryVariable(ArraySymbol):
         return x
 
     @classmethod
-    def _from_zipfile(cls, zf, directory, Model model, predecessors):
+    def _from_zipfile(cls, zf, directory, _Model model, predecessors):
         """Construct a binary symbol from a compressed file.
 
         Args:
@@ -807,7 +807,7 @@ cdef class Constant(ArraySymbol):
         >>> type(a)
         <class 'dwave.optimization.symbols.Constant'>
     """
-    def __init__(self, Model model, array_like):
+    def __init__(self, _Model model, array_like):
         # In the future we won't need to be contiguous, but we do need to be right now
         array = np.asarray(array_like, dtype=np.double, order="C")
 
@@ -908,7 +908,7 @@ cdef class Constant(ArraySymbol):
         return constant
 
     @classmethod
-    def _from_zipfile(cls, zf, directory, Model model, predecessors):
+    def _from_zipfile(cls, zf, directory, _Model model, predecessors):
         """Construct a constant symbol from a compressed file.
 
         Args:
@@ -1005,7 +1005,7 @@ cdef class DisjointBitSets(Symbol):
         <class 'dwave.optimization.symbols.DisjointBitSets'>
     """
     def __init__(
-        self, Model model, Py_ssize_t primary_set_size, Py_ssize_t num_disjoint_sets
+        self, _Model model, Py_ssize_t primary_set_size, Py_ssize_t num_disjoint_sets
     ):
         # Get an observing pointer to the node
         self.ptr = model._graph.emplace_node[cppDisjointBitSetsNode](
@@ -1026,7 +1026,7 @@ cdef class DisjointBitSets(Symbol):
         return x
 
     @classmethod
-    def _from_zipfile(cls, zf, directory, Model model, predecessors):
+    def _from_zipfile(cls, zf, directory, _Model model, predecessors):
         """Construct a disjoint-sets symbol from a compressed file.
 
         Args:
@@ -1166,7 +1166,7 @@ cdef class DisjointBitSet(ArraySymbol):
         if set_index > <Py_ssize_t>(parent.ptr.successors().size()):
             raise ValueError("`DisjointBitSet`s must be created successively")
 
-        cdef Model model = parent.model
+        cdef _Model model = parent.model
         if set_index == <Py_ssize_t>(parent.ptr.successors().size()):
             # The DisjointBitSet has not been added to the model yet, so add it
             self.ptr = model._graph.emplace_node[cppDisjointBitSetNode](parent.ptr)
@@ -1190,7 +1190,7 @@ cdef class DisjointBitSet(ArraySymbol):
         return x
 
     @classmethod
-    def _from_zipfile(cls, zf, directory, Model model, predecessors):
+    def _from_zipfile(cls, zf, directory, _Model model, predecessors):
         """Construct a disjoint-set symbol from a compressed file.
 
         Args:
@@ -1262,7 +1262,7 @@ cdef class DisjointLists(Symbol):
         <class 'dwave.optimization.symbols.DisjointLists'>
     """
     def __init__(
-        self, Model model, Py_ssize_t primary_set_size, Py_ssize_t num_disjoint_lists
+        self, _Model model, Py_ssize_t primary_set_size, Py_ssize_t num_disjoint_lists
     ):
         # Get an observing pointer to the node
         self.ptr = model._graph.emplace_node[cppDisjointListsNode](
@@ -1282,7 +1282,7 @@ cdef class DisjointLists(Symbol):
         return x
 
     @classmethod
-    def _from_zipfile(cls, zf, directory, Model model, predecessors):
+    def _from_zipfile(cls, zf, directory, _Model model, predecessors):
         """Construct a disjoint-lists symbol from a compressed file.
 
         Args:
@@ -1419,7 +1419,7 @@ cdef class DisjointList(ArraySymbol):
         if list_index > <Py_ssize_t>(parent.ptr.successors().size()):
             raise ValueError("`DisjointList`s must be created successively")
 
-        cdef Model model = parent.model
+        cdef _Model model = parent.model
         if list_index == <Py_ssize_t>(parent.ptr.successors().size()):
             # The DisjointListNode has not been added to the model yet, so add it
             self.ptr = model._graph.emplace_node[cppDisjointListNode](parent.ptr)
@@ -1443,7 +1443,7 @@ cdef class DisjointList(ArraySymbol):
         return x
 
     @classmethod
-    def _from_zipfile(cls, zf, directory, Model model, predecessors):
+    def _from_zipfile(cls, zf, directory, _Model model, predecessors):
         """Construct a disjoint-list symbol from a compressed file.
 
         Args:
@@ -1520,7 +1520,7 @@ cdef class Equal(ArraySymbol):
         if lhs.model is not rhs.model:
             raise ValueError("lhs and rhs do not share the same underlying model")
 
-        cdef Model model = lhs.model
+        cdef _Model model = lhs.model
 
         self.ptr = model._graph.emplace_node[cppEqualNode](lhs.array_ptr, rhs.array_ptr)
         self.initialize_arraynode(model, self.ptr)
@@ -1581,7 +1581,7 @@ cdef class IntegerVariable(ArraySymbol):
         >>> type(i)
         <class 'dwave.optimization.symbols.IntegerVariable'>
     """
-    def __init__(self, Model model, shape=None, lower_bound=None, upper_bound=None):
+    def __init__(self, _Model model, shape=None, lower_bound=None, upper_bound=None):
         cdef vector[Py_ssize_t] vshape = _as_cppshape(tuple() if shape is None else shape )
 
         if lower_bound is None and upper_bound is None:
@@ -1607,7 +1607,7 @@ cdef class IntegerVariable(ArraySymbol):
         return x
 
     @classmethod
-    def _from_zipfile(cls, zf, directory, Model model, predecessors):
+    def _from_zipfile(cls, zf, directory, _Model model, predecessors):
         if predecessors:
             raise ValueError(f"{cls.__name__} cannot have predecessors")
 
@@ -1689,7 +1689,7 @@ cdef class LessEqual(ArraySymbol):
         if lhs.model is not rhs.model:
             raise ValueError("lhs and rhs do not share the same underlying model")
 
-        cdef Model model = lhs.model
+        cdef _Model model = lhs.model
 
         self.ptr = model._graph.emplace_node[cppLessEqualNode](lhs.array_ptr, rhs.array_ptr)
         self.initialize_arraynode(model, self.ptr)
@@ -1721,7 +1721,7 @@ cdef class ListVariable(ArraySymbol):
         >>> type(l)
         <class 'dwave.optimization.symbols.ListVariable'>
     """
-    def __init__(self, Model model, Py_ssize_t n):
+    def __init__(self, _Model model, Py_ssize_t n):
         # Get an observing pointer to the node
         self.ptr = model._graph.emplace_node[cppListNode](n)
 
@@ -1739,7 +1739,7 @@ cdef class ListVariable(ArraySymbol):
         return x
 
     @classmethod
-    def _from_zipfile(cls, zf, directory, Model model, predecessors):
+    def _from_zipfile(cls, zf, directory, _Model model, predecessors):
         if predecessors:
             raise ValueError(f"{cls.__name__} cannot have predecessors")
 
@@ -1799,7 +1799,7 @@ cdef class Logical(ArraySymbol):
         :func:`~dwave.optimization.mathematical.logical`: equivalent function.
     """
     def __init__(self, ArraySymbol x):
-        cdef Model model = x.model
+        cdef _Model model = x.model
 
         self.ptr = model._graph.emplace_node[cppLogicalNode](x.array_ptr)
         self.initialize_arraynode(model, self.ptr)
@@ -1834,7 +1834,7 @@ cdef class Max(ArraySymbol):
         <class 'dwave.optimization.symbols.Max'>
     """
     def __init__(self, ArraySymbol node):
-        cdef Model model = node.model
+        cdef _Model model = node.model
 
         self.ptr = model._graph.emplace_node[cppMaxNode](node.array_ptr)
 
@@ -1876,7 +1876,7 @@ cdef class Maximum(ArraySymbol):
         if lhs.model is not rhs.model:
             raise ValueError("lhs and rhs do not share the same underlying model")
 
-        cdef Model model = lhs.model
+        cdef _Model model = lhs.model
 
         self.ptr = model._graph.emplace_node[cppMaximumNode](lhs.array_ptr, rhs.array_ptr)
 
@@ -1912,7 +1912,7 @@ cdef class Min(ArraySymbol):
         <class 'dwave.optimization.symbols.Min'>
     """
     def __init__(self, ArraySymbol node):
-        cdef Model model = node.model
+        cdef _Model model = node.model
 
         self.ptr = model._graph.emplace_node[cppMinNode](node.array_ptr)
 
@@ -1954,7 +1954,7 @@ cdef class Minimum(ArraySymbol):
         if lhs.model is not rhs.model:
             raise ValueError("lhs and rhs do not share the same underlying model")
 
-        cdef Model model = lhs.model
+        cdef _Model model = lhs.model
 
         self.ptr = model._graph.emplace_node[cppMinimumNode](lhs.array_ptr, rhs.array_ptr)
 
@@ -1993,7 +1993,7 @@ cdef class Modulus(ArraySymbol):
         if lhs.model is not rhs.model:
             raise ValueError("lhs and rhs do not share the same underlying model")
 
-        cdef Model model = lhs.model
+        cdef _Model model = lhs.model
 
         self.ptr = model._graph.emplace_node[cppModulusNode](lhs.array_ptr, rhs.array_ptr)
         self.initialize_arraynode(model, self.ptr)
@@ -2031,7 +2031,7 @@ cdef class Multiply(ArraySymbol):
         if lhs.model is not rhs.model:
             raise ValueError("lhs and rhs do not share the same underlying model")
 
-        cdef Model model = lhs.model
+        cdef _Model model = lhs.model
 
         self.ptr = model._graph.emplace_node[cppMultiplyNode](lhs.array_ptr, rhs.array_ptr)
         self.initialize_arraynode(model, self.ptr)
@@ -2072,7 +2072,7 @@ cdef class NaryAdd(ArraySymbol):
         if len(inputs) == 0:
             raise TypeError("must have at least one predecessor node")
 
-        cdef Model model = inputs[0].model
+        cdef _Model model = inputs[0].model
         cdef vector[cppArrayNode*] cppinputs
 
         cdef ArraySymbol array
@@ -2129,7 +2129,7 @@ cdef class NaryMaximum(ArraySymbol):
         if len(inputs) == 0:
             raise TypeError("must have at least one predecessor node")
 
-        cdef Model model = inputs[0].model
+        cdef _Model model = inputs[0].model
         cdef vector[cppArrayNode*] cppinputs
 
         cdef ArraySymbol array
@@ -2179,7 +2179,7 @@ cdef class NaryMinimum(ArraySymbol):
         if len(inputs) == 0:
             raise TypeError("must have at least one predecessor node")
 
-        cdef Model model = inputs[0].model
+        cdef _Model model = inputs[0].model
         cdef vector[cppArrayNode*] cppinputs
 
         cdef ArraySymbol array
@@ -2228,7 +2228,7 @@ cdef class NaryMultiply(ArraySymbol):
         if len(inputs) == 0:
             raise TypeError("must have at least one predecessor node")
 
-        cdef Model model = inputs[0].model
+        cdef _Model model = inputs[0].model
         cdef vector[cppArrayNode*] cppinputs
 
         cdef ArraySymbol array
@@ -2295,7 +2295,7 @@ cdef class NaryReduce(ArraySymbol):
         if len(initial_values) != expression.num_inputs():
             raise ValueError("must have same number of initial values as inputs")
 
-        cdef Model model = operands[0].model
+        cdef _Model model = operands[0].model
         cdef cppArrayNode* output = expression.output.array_ptr
         cdef vector[double] cppinitial_values
         cdef cppInputNode* cppinput
@@ -2369,7 +2369,7 @@ cdef class Negative(ArraySymbol):
         <class 'dwave.optimization.symbols.Negative'>
     """
     def __init__(self, ArraySymbol x):
-        cdef Model model = x.model
+        cdef _Model model = x.model
 
         self.ptr = model._graph.emplace_node[cppNegativeNode](x.array_ptr)
         self.initialize_arraynode(model, self.ptr)
@@ -2396,7 +2396,7 @@ cdef class Not(ArraySymbol):
         :func:`~dwave.optimization.mathematical.logical_not`: equivalent function.
     """
     def __init__(self, ArraySymbol x):
-        cdef Model model = x.model
+        cdef _Model model = x.model
 
         self.ptr = model._graph.emplace_node[cppNotNode](x.array_ptr)
         self.initialize_arraynode(model, self.ptr)
@@ -2426,7 +2426,7 @@ cdef class Or(ArraySymbol):
         if lhs.model is not rhs.model:
             raise ValueError("lhs and rhs do not share the same underlying model")
 
-        cdef Model model = lhs.model
+        cdef _Model model = lhs.model
 
         self.ptr = model._graph.emplace_node[cppOrNode](lhs.array_ptr, rhs.array_ptr)
         self.initialize_arraynode(model, self.ptr)
@@ -2461,7 +2461,7 @@ cdef class PartialSum(ArraySymbol):
         <class 'dwave.optimization.symbols.PartialSum'>
     """
     def __init__(self, ArraySymbol array, int axis):
-        cdef Model model = array.model
+        cdef _Model model = array.model
         self.ptr = model._graph.emplace_node[cppPartialSumNode](array.array_ptr, axis)
         self.initialize_arraynode(model, self.ptr)
 
@@ -2480,7 +2480,7 @@ cdef class PartialSum(ArraySymbol):
         return ps
 
     @classmethod
-    def _from_zipfile(cls, zf, directory, Model model, predecessors):
+    def _from_zipfile(cls, zf, directory, _Model model, predecessors):
         if len(predecessors) != 1:
             raise ValueError("PartialSum must have exactly one predecessor")
 
@@ -2550,7 +2550,7 @@ cdef class Prod(ArraySymbol):
         <class 'dwave.optimization.symbols.Prod'>
     """
     def __init__(self, ArraySymbol node):
-        cdef Model model = node.model
+        cdef _Model model = node.model
 
         self.ptr = model._graph.emplace_node[cppProdNode](node.array_ptr)
 
@@ -2735,7 +2735,7 @@ cdef class QuadraticModel(ArraySymbol):
         return self.ptr.get_quadratic_model().get_quadratic(u, v)
 
     @classmethod
-    def _from_zipfile(cls, zf, directory, Model model, predecessors):
+    def _from_zipfile(cls, zf, directory, _Model model, predecessors):
         """Construct a QuadraticModel from a zipfile."""
         if len(predecessors) != 1:
             raise ValueError("Reshape must have exactly one predecessor")
@@ -2816,7 +2816,7 @@ cdef class Reshape(ArraySymbol):
         <class 'dwave.optimization.symbols.Reshape'>
     """
     def __init__(self, ArraySymbol node, shape):
-        cdef Model model = node.model
+        cdef _Model model = node.model
 
         self.ptr = model._graph.emplace_node[cppReshapeNode](
             node.array_ptr,
@@ -2837,7 +2837,7 @@ cdef class Reshape(ArraySymbol):
         return m
 
     @classmethod
-    def _from_zipfile(cls, zf, directory, Model model, predecessors):
+    def _from_zipfile(cls, zf, directory, _Model model, predecessors):
         if len(predecessors) != 1:
             raise ValueError("Reshape must have exactly one predecessor")
 
@@ -2873,7 +2873,7 @@ cdef class SetVariable(ArraySymbol):
         >>> type(s)
         <class 'dwave.optimization.symbols.SetVariable'>
     """
-    def __init__(self, Model model, Py_ssize_t n, Py_ssize_t min_size, Py_ssize_t max_size):
+    def __init__(self, _Model model, Py_ssize_t n, Py_ssize_t min_size, Py_ssize_t max_size):
         self.ptr = model._graph.emplace_node[cppSetNode](n, min_size, max_size)
         self.initialize_arraynode(model, self.ptr)
 
@@ -2889,7 +2889,7 @@ cdef class SetVariable(ArraySymbol):
         return x
 
     @classmethod
-    def _from_zipfile(cls, zf, directory, Model model, predecessors):
+    def _from_zipfile(cls, zf, directory, _Model model, predecessors):
         if predecessors:
             raise ValueError(f"{cls.__name__} cannot have predecessors")
 
@@ -2952,7 +2952,7 @@ _register(SetVariable, typeid(cppSetNode))
 
 cdef class Size(ArraySymbol):
     def __init__(self, ArraySymbol array):
-        cdef Model model = array.model
+        cdef _Model model = array.model
 
         self.ptr = model._graph.emplace_node[cppSizeNode](array.array_ptr)
         self.initialize_arraynode(array.model, self.ptr)
@@ -2989,7 +2989,7 @@ cdef class Square(ArraySymbol):
         <class 'dwave.optimization.symbols.Square'>
     """
     def __init__(self, ArraySymbol x):
-        cdef Model model = x.model
+        cdef _Model model = x.model
 
         self.ptr = model._graph.emplace_node[cppSquareNode](x.array_ptr)
         self.initialize_arraynode(model, self.ptr)
@@ -3024,7 +3024,7 @@ cdef class SquareRoot(ArraySymbol):
         <class 'dwave.optimization.symbols.SquareRoot'>
     """
     def __init__(self, ArraySymbol x):
-        cdef Model model = x.model
+        cdef _Model model = x.model
 
         self.ptr = model._graph.emplace_node[cppSquareRootNode](x.array_ptr)
         self.initialize_arraynode(model, self.ptr)
@@ -3063,7 +3063,7 @@ cdef class Subtract(ArraySymbol):
         if lhs.model is not rhs.model:
             raise ValueError("lhs and rhs do not share the same underlying model")
 
-        cdef Model model = lhs.model
+        cdef _Model model = lhs.model
 
         self.ptr = model._graph.emplace_node[cppSubtractNode](lhs.array_ptr, rhs.array_ptr)
         self.initialize_arraynode(model, self.ptr)
@@ -3098,7 +3098,7 @@ cdef class Sum(ArraySymbol):
         <class 'dwave.optimization.symbols.Sum'>
     """
     def __init__(self, ArraySymbol array):
-        cdef Model model = array.model
+        cdef _Model model = array.model
         self.ptr = model._graph.emplace_node[cppSumNode](array.array_ptr)
         self.initialize_arraynode(model, self.ptr)
 
@@ -3124,7 +3124,7 @@ cdef class Where(ArraySymbol):
         :func:`~dwave.optimization.mathematical.where`: equivalent function.
     """
     def __init__(self, ArraySymbol condition, ArraySymbol x, ArraySymbol y):
-        cdef Model model = condition.model
+        cdef _Model model = condition.model
 
         if condition.model is not x.model:
             raise ValueError("condition and x do not share the same underlying model")
@@ -3162,7 +3162,7 @@ cdef class Xor(ArraySymbol):
         if lhs.model is not rhs.model:
             raise ValueError("lhs and rhs do not share the same underlying model")
 
-        cdef Model model = lhs.model
+        cdef _Model model = lhs.model
 
         self.ptr = model._graph.emplace_node[cppXorNode](lhs.array_ptr, rhs.array_ptr)
         self.initialize_arraynode(model, self.ptr)
diff --git a/meson.build b/meson.build
index 5e3e8336..dc807f03 100644
--- a/meson.build
+++ b/meson.build
@@ -66,8 +66,8 @@ else
 endif
 
 py.extension_module(
-    'model',
-    'dwave/optimization/model.pyx',
+    '_model',
+    'dwave/optimization/_model.pyx',
     dependencies: libdwave_optimization,
     gnu_symbol_visibility: 'default',
     install: true,
diff --git a/tests/test_model.py b/tests/test_model.py
index bfa8353e..476dbd9c 100644
--- a/tests/test_model.py
+++ b/tests/test_model.py
@@ -22,18 +22,19 @@
 
 import numpy as np
 
+import dwave.optimization._model
 import dwave.optimization.symbols
 from dwave.optimization import Model
 
 
 class TestArraySymbol(unittest.TestCase):
     def test_abstract(self):
-        from dwave.optimization.model import ArraySymbol
+        from dwave.optimization._model import ArraySymbol
         with self.assertRaisesRegex(ValueError, "ArraySymbols cannot be constructed directly"):
             ArraySymbol()
 
     def test_bool(self):
-        from dwave.optimization.model import ArraySymbol
+        from dwave.optimization._model import ArraySymbol
 
         # bypass the init, this should be done very carefully as it can lead to
         # segfaults dependig on what methods are accessed!
@@ -107,7 +108,7 @@ class UnknownType():
 
         with self.subTest("__pow__"):
             self.assertIsInstance(x ** 1, type(x))
-            self.assertIsInstance(x ** 1, dwave.optimization.model.ArraySymbol)
+            self.assertIsInstance(x ** 1, dwave.optimization._model.ArraySymbol)
             self.assertIsInstance(x ** 2, dwave.optimization.symbols.Square)
             self.assertIsInstance(x ** 3, dwave.optimization.symbols.NaryMultiply)
             self.assertIsInstance(x ** 4, dwave.optimization.symbols.NaryMultiply)
@@ -121,7 +122,7 @@ def __init__(self, array):
         def __getitem__(self, index):
             if not isinstance(index, tuple):
                 return self[(index,)]
-            i0, i1 = dwave.optimization.model._split_indices(index)
+            i0, i1 = dwave.optimization._model._split_indices(index)
             np.testing.assert_array_equal(self.array[index], self.array[i0][i1])
 
     def test_split_indices(self):
@@ -677,7 +678,7 @@ def test_serialization_bad(self):
 
 class TestSymbol(unittest.TestCase):
     def test_abstract(self):
-        from dwave.optimization.model import Symbol
+        from dwave.optimization._model import Symbol
         with self.assertRaisesRegex(ValueError, "Symbols cannot be constructed directly"):
             Symbol()