diff --git a/.github/workflows/test_sbml_semantic_test_suite.yml b/.github/workflows/test_sbml_semantic_test_suite.yml
index b3394d54ea..891ddb91d4 100644
--- a/.github/workflows/test_sbml_semantic_test_suite.yml
+++ b/.github/workflows/test_sbml_semantic_test_suite.yml
@@ -8,8 +8,7 @@ on:
     paths:
       - .github/workflows/test_sbml_semantic_test_suite.yml
       - python/sdist/amici/exporters/sundials/de_export.py
-      - python/sdist/amici/de_model_components.py
-      - python/sdist/amici/de_model.py
+      - python/sdist/amici/_symbolic/**
       - python/sdist/amici/importers/sbml/**
       - python/sdist/amici/importers/utils.py
       - scripts/run-SBMLTestsuite.sh
diff --git a/.gitignore b/.gitignore
index 06d5fa7197..4a3df8cbdd 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,6 +7,7 @@ swig/python/build/*
 tests/cpp/build/*
 tests/cpp/build_xcode/*
 tests/cpp/Testing/*
+.cache/
 
 doc-venv/*
 fonts/*
diff --git a/doc/python_modules.rst b/doc/python_modules.rst
index 3df9216bbc..b750d9adad 100644
--- a/doc/python_modules.rst
+++ b/doc/python_modules.rst
@@ -29,8 +29,6 @@ AMICI Python API
    amici.importers.petab.v1.simulator
    amici.jax
    amici.exporters.sundials.de_export
-   amici.de_model
-   amici.de_model_components
    amici.plotting
    amici.pandas
    amici.logging
diff --git a/python/sdist/amici/_symbolic/__init__.py b/python/sdist/amici/_symbolic/__init__.py
new file mode 100644
index 0000000000..fe86dd0a2c
--- /dev/null
+++ b/python/sdist/amici/_symbolic/__init__.py
@@ -0,0 +1,13 @@
+"""
+Private functionality for symbolic model generation and manipulation.
+
+This module provides tools to create, modify, and analyze symbolic
+representations of differential equation models used in AMICI.
+Most symbolic functionality is provided by `SymPy <https://www.sympy.org/>`_.
+"""
+
+from .de_model import *  # noqa: F403
+from .de_model_components import *  # noqa: F403
+from .sympy_utils import *  # noqa: F403
+
+__all__ = de_model.__all__ + de_model_components.__all__ + sympy_utils.__all__
diff --git a/python/sdist/amici/de_model.py b/python/sdist/amici/_symbolic/de_model.py
similarity index 99%
rename from python/sdist/amici/de_model.py
rename to python/sdist/amici/_symbolic/de_model.py
index 31a9650938..542ee8284d 100644
--- a/python/sdist/amici/de_model.py
+++ b/python/sdist/amici/_symbolic/de_model.py
@@ -16,6 +16,22 @@
 import sympy as sp
 from sympy import ImmutableDenseMatrix, MutableDenseMatrix
 
+from amici.exporters.sundials.cxx_functions import (
+    nobody_functions,
+    sensi_functions,
+    sparse_functions,
+    var_in_function_signature,
+)
+from amici.exporters.sundials.cxxcodeprinter import csc_matrix
+from amici.importers.utils import (
+    ObservableTransformation,
+    _default_simplify,
+    amici_time_symbol,
+    toposort_symbols,
+    unique_preserve_order,
+)
+from amici.logging import get_logger, log_execution_time, set_log_level
+
 from .de_model_components import (
     AlgebraicEquation,
     AlgebraicState,
@@ -39,21 +55,6 @@
     SigmaZ,
     State,
 )
-from .exporters.sundials.cxx_functions import (
-    nobody_functions,
-    sensi_functions,
-    sparse_functions,
-    var_in_function_signature,
-)
-from .exporters.sundials.cxxcodeprinter import csc_matrix
-from .importers.utils import (
-    ObservableTransformation,
-    _default_simplify,
-    amici_time_symbol,
-    toposort_symbols,
-    unique_preserve_order,
-)
-from .logging import get_logger, log_execution_time, set_log_level
 from .sympy_utils import (
     _parallel_applyfunc,
     smart_is_zero_matrix,
diff --git a/python/sdist/amici/de_model_components.py b/python/sdist/amici/_symbolic/de_model_components.py
similarity index 99%
rename from python/sdist/amici/de_model_components.py
rename to python/sdist/amici/_symbolic/de_model_components.py
index c48e2d28c6..23ec0cac8d 100644
--- a/python/sdist/amici/de_model_components.py
+++ b/python/sdist/amici/_symbolic/de_model_components.py
@@ -6,8 +6,8 @@
 
 import sympy as sp
 
-from .constants import SymbolId
-from .importers.utils import (
+from amici.constants import SymbolId
+from amici.importers.utils import (
     RESERVED_SYMBOLS,
     ObservableTransformation,
     amici_time_symbol,
diff --git a/python/sdist/amici/sympy_utils.py b/python/sdist/amici/_symbolic/sympy_utils.py
similarity index 96%
rename from python/sdist/amici/sympy_utils.py
rename to python/sdist/amici/_symbolic/sympy_utils.py
index 85a0e986a8..8e2d75d977 100644
--- a/python/sdist/amici/sympy_utils.py
+++ b/python/sdist/amici/_symbolic/sympy_utils.py
@@ -10,10 +10,19 @@
 
 import sympy as sp
 
-from .logging import get_logger, log_execution_time
+from amici.logging import get_logger, log_execution_time
 
 logger = get_logger(__name__, logging.ERROR)
 
+__all__ = [
+    "smart_jacobian",
+    "smart_multiply",
+    "smart_is_zero_matrix",
+    "_monkeypatch_sympy",
+    "_parallel_applyfunc",
+    "_piecewise_to_minmax",
+]
+
 
 def _custom_pow_eval_derivative(self, s):
     """
diff --git a/python/sdist/amici/exporters/sundials/de_export.py b/python/sdist/amici/exporters/sundials/de_export.py
index 06f5ecca68..4fdfc43feb 100644
--- a/python/sdist/amici/exporters/sundials/de_export.py
+++ b/python/sdist/amici/exporters/sundials/de_export.py
@@ -24,20 +24,15 @@
 
 import sympy as sp
 
+from amici._symbolic import *
+
 from ... import (
     __commit__,
     __version__,
     get_model_dir,
 )
-from ...de_model import DEModel
-from ...de_model_components import *
 from ...importers.sbml import splines
 from ...logging import get_logger, log_execution_time, set_log_level
-from ...sympy_utils import (
-    _custom_pow_eval_derivative,
-    _monkeypatched,
-    smart_is_zero_matrix,
-)
 from ..template import apply_template
 from .compile import build_model_extension
 from .cxx_functions import (
@@ -236,16 +231,14 @@ def __init__(
         self.generate_sensitivity_code: bool = generate_sensitivity_code
         self.hybridisation = hybridization
 
+    @_monkeypatch_sympy
     @log_execution_time("generating cpp code", logger)
     def generate_model_code(self) -> None:
         """
         Generates the model code (Python package + C++ extension).
         """
-        with _monkeypatched(
-            sp.Pow, "_eval_derivative", _custom_pow_eval_derivative
-        ):
-            self._prepare_model_folder()
-            self._generate_c_code()
+        self._prepare_model_folder()
+        self._generate_c_code()
 
     @log_execution_time("compiling cpp code", logger)
     def compile_model(self) -> None:
diff --git a/python/sdist/amici/importers/petab/_petab_importer.py b/python/sdist/amici/importers/petab/_petab_importer.py
index 34ca4df1c9..263f18b0eb 100644
--- a/python/sdist/amici/importers/petab/_petab_importer.py
+++ b/python/sdist/amici/importers/petab/_petab_importer.py
@@ -22,7 +22,7 @@
 
 import amici
 from amici import SensitivityOrder, get_model_dir
-from amici.de_model import DEModel
+from amici._symbolic import DEModel, Event
 from amici.importers.utils import MeasurementChannel, amici_time_symbol
 from amici.logging import get_logger
 
@@ -32,8 +32,6 @@
 if TYPE_CHECKING:
     import pysb
 
-    import amici.de_model_components
-
 __all__ = [
     "PetabImporter",
     "ExperimentManager",
@@ -1926,7 +1924,7 @@ def __init__(self, petab_problem: v2.Problem):
                     )
         #: The PEtab problem to convert. Not modified by this class.
         self._petab_problem = petab_problem
-        self._events: list[amici.de_model_components.Event] = []
+        self._events: list[Event] = []
         self._new_problem: v2.Problem | None = None
 
     @staticmethod
@@ -1936,7 +1934,7 @@ def _get_experiment_indicator_condition_id(experiment_id: str) -> str:
 
     def convert(
         self,
-    ) -> tuple[v2.Problem, list[amici.de_model_components.Event]]:
+    ) -> tuple[v2.Problem, list[Event]]:
         """Convert PEtab experiments to amici events and pysb initials.
 
         Generate events, add Initials, or convert Parameters to Expressions
@@ -1952,7 +1950,7 @@ def convert(
             to be passed to `pysb2amici`.
         """
         self._new_problem = copy.deepcopy(self._petab_problem)
-        self._events: list[amici.de_model_components.Event] = []
+        self._events: list[Event] = []
 
         self._add_preequilibration_indicator()
 
@@ -2177,7 +2175,7 @@ def _create_period_start_event(
             else:
                 raise AssertionError(change)
 
-        event = amici.de_model_components.Event(
+        event = Event(
             identifier=sp.Symbol(event_id),
             name=event_id,
             value=root_fun,
diff --git a/python/sdist/amici/importers/pysb/__init__.py b/python/sdist/amici/importers/pysb/__init__.py
index 235ae8b72c..573a556e2d 100644
--- a/python/sdist/amici/importers/pysb/__init__.py
+++ b/python/sdist/amici/importers/pysb/__init__.py
@@ -23,8 +23,8 @@
 import sympy as sp
 
 import amici
-from amici.de_model import DEModel
-from amici.de_model_components import (
+from amici._symbolic.de_model import DEModel
+from amici._symbolic.de_model_components import (
     DifferentialState,
     Event,
     Expression,
@@ -486,7 +486,7 @@ def get_cached_index(symbol, sarray, index_cache):
     dflux_dw = sp.ImmutableSparseMatrix(n_r, n_w, dflux_dw_dict)
     dflux_dp = sp.ImmutableSparseMatrix(n_r, n_p, dflux_dp_dict)
 
-    # use dok format to convert numeric csc to sparse symbolic
+    # use dok format to convert numeric csc to sparse _symbolic
     S = sp.ImmutableSparseMatrix(
         n_x,
         n_r,  # don't use shape here as we are eliminating rows
@@ -649,7 +649,7 @@ def _add_expression(
         name of the expression
 
     :param expr:
-        symbolic expression that the symbol refers to
+        _symbolic expression that the symbol refers to
 
     :param pysb_model:
         see :py:func:`_process_pysb_expressions`
@@ -681,10 +681,10 @@ def _add_expression(
 
         y = sp.Symbol(name)
         trafo = noise_distribution_to_observable_transformation(noise_dist)
-        # note that this is a bit iffy since we are potentially using the same symbolic identifier in expressions (w)
+        # note that this is a bit iffy since we are potentially using the same _symbolic identifier in expressions (w)
         # and observables (y). This is not a problem as there currently are no model functions that use both. If this
         # changes, I would expect symbol redefinition warnings in CPP models and overwriting in JAX models, but as both
-        # symbols refer to the same symbolic entity, this should not be a problem (untested)
+        # symbols refer to the same _symbolic entity, this should not be a problem (untested)
         obs = Observable(
             y, name, _parse_special_functions(expr), transformation=trafo
         )
@@ -723,7 +723,7 @@ def _get_sigma(
     pysb_model: pysb.Model, obs_name: str, sigma_name: str | None
 ) -> sp.Symbol:
     """
-    Tries to extract standard deviation symbolic identifier and formula
+    Tries to extract standard deviation _symbolic identifier and formula
     for a given observable name from the pysb model and if no specification is
     available sets default values
 
@@ -738,7 +738,7 @@ def _get_sigma(
         sigma or ``None``.
 
     :return:
-        symbolic variable representing the standard deviation of the observable
+        _symbolic variable representing the standard deviation of the observable
     """
     if sigma_name is None:
         return sp.Symbol(f"sigma_{obs_name}")
@@ -933,7 +933,7 @@ def _compute_possible_indices(
                 )
                 # TODO: implement this, multiply species by the volume of
                 # their respective compartment and allow total_cl to depend
-                # on parameters + constants and update the respective symbolic
+                # on parameters + constants and update the respective _symbolic
                 # derivative accordingly
 
             prototype = dict()
@@ -1161,7 +1161,7 @@ def _is_in_cycle(
 def _greedy_target_index_update(cl_prototypes: CL_Prototype) -> None:
     """
     Computes unique target indices for conservation laws from possible
-    indices  such that expected fill in in symbolic derivatives is minimized
+    indices  such that expected fill in in _symbolic derivatives is minimized
 
     :param cl_prototypes:
         dict that contains possible indices and non-unique target indices
diff --git a/python/sdist/amici/importers/sbml/__init__.py b/python/sdist/amici/importers/sbml/__init__.py
index b67b4fcb43..346d3d0bfb 100644
--- a/python/sdist/amici/importers/sbml/__init__.py
+++ b/python/sdist/amici/importers/sbml/__init__.py
@@ -30,14 +30,20 @@
 
 import amici
 from amici import get_model_dir, has_clibs
-from amici.constants import SymbolId
-from amici.de_model import DEModel
-from amici.de_model_components import (
+from amici._symbolic.de_model import DEModel
+from amici._symbolic.de_model_components import (
     DifferentialState,
     Event,
     Expression,
     symbol_to_type,
 )
+from amici._symbolic.sympy_utils import (
+    _monkeypatch_sympy,
+    _piecewise_to_minmax,
+    smart_is_zero_matrix,
+    smart_multiply,
+)
+from amici.constants import SymbolId
 from amici.importers.sbml.splines import AbstractSpline
 from amici.importers.sbml.utils import SBMLException
 from amici.importers.utils import (
@@ -64,12 +70,6 @@
     toposort_symbols,
 )
 from amici.logging import get_logger, log_execution_time, set_log_level
-from amici.sympy_utils import (
-    _monkeypatch_sympy,
-    _piecewise_to_minmax,
-    smart_is_zero_matrix,
-    smart_multiply,
-)
 
 SymbolicFormula = dict[sp.Symbol, sp.Expr]
 
@@ -91,7 +91,7 @@ class SbmlImporter:
         displayed
 
     :ivar symbols:
-        dict carrying symbolic definitions
+        dict carrying _symbolic definitions
 
     :ivar sbml_reader:
 
@@ -127,11 +127,11 @@ class SbmlImporter:
 
     :ivar species_assignment_rules:
         Assignment rules for species.
-        Key is symbolic identifier and value is assignment value
+        Key is _symbolic identifier and value is assignment value
 
     :ivar compartment_assignment_rules:
         Assignment rules for compartments.
-        Key is symbolic identifier and value is assignment value
+        Key is _symbolic identifier and value is assignment value
 
     :ivar parameter_assignment_rules:
         assignment rules for parameters, these parameters are not permissible
@@ -879,7 +879,7 @@ def _gather_locals(self, hardcode_symbols: Sequence[str] = None) -> None:
 
         This is later used during sympifications to avoid sympy builtins
         shadowing model entities as well as to avoid possibly costly
-        symbolic substitutions
+        _symbolic substitutions
         """
         self._gather_base_locals(hardcode_symbols=hardcode_symbols)
         self._gather_dependent_locals()
@@ -1291,7 +1291,7 @@ def _process_parameters(
                 }
 
         # Parameters that need to be turned into expressions or species
-        #  so far, this concerns parameters with symbolic initial assignments
+        #  so far, this concerns parameters with _symbolic initial assignments
         #  (those have been skipped above) that are not rate rule targets
 
         # Set of symbols in initial assignments that still allows handling them
@@ -1824,7 +1824,7 @@ def _process_events(self) -> None:
                 except TypeError:
                     raise SBMLException(
                         "Could not process event assignment for "
-                        f"{str(variable_sym)}. AMICI only allows symbolic "
+                        f"{str(variable_sym)}. AMICI only allows _symbolic "
                         "expressions as event assignments."
                     )
                 if variable_sym in concentration_species_by_compartment:
@@ -1898,7 +1898,7 @@ def _process_observation_model(
         observation_model: list[MeasurementChannel] = None,
     ) -> None:
         """
-        Perform symbolic computations required for observable and objective
+        Perform _symbolic computations required for observable and objective
         function evaluation.
 
         Add user-provided observables or make all species, and compartments
@@ -2084,7 +2084,7 @@ def _process_log_likelihood(
         event_reg: bool = False,
     ):
         """
-        Perform symbolic computations required for objective function
+        Perform _symbolic computations required for objective function
         evaluation.
 
         :param sigmas:
@@ -2255,7 +2255,7 @@ def _make_initial(
         replacing species by their initial values.
 
         :param sym_math:
-            symbolic expression
+            _symbolic expression
         :return:
             transformed expression
         """
@@ -2549,7 +2549,7 @@ def _add_conservation_for_non_constant_species(
                 compartment_size = self.compartments[symbol["compartment"]]
             all_compartment_sizes.append(compartment_size)
 
-        # iterate over list of conservation laws, create symbolic expressions,
+        # iterate over list of conservation laws, create _symbolic expressions,
         for target_state_model_idx, state_idxs, coefficients in raw_cls:
             if all_state_ids[target_state_model_idx] in eliminated_state_ids:
                 # constants state, already eliminated
@@ -2607,13 +2607,13 @@ def _replace_in_all_expressions(
         self, old: sp.Symbol, new: sp.Expr, replace_identifiers=False
     ) -> None:
         """
-        Replace 'old' by 'new' in all symbolic expressions.
+        Replace 'old' by 'new' in all _symbolic expressions.
 
         :param old:
-            symbolic variables to be replaced
+            _symbolic variables to be replaced
 
         :param new:
-            replacement symbolic variables
+            replacement _symbolic variables
         """
         fields = [
             "stoichiometric_matrix",
@@ -2918,7 +2918,7 @@ def _get_element_stoichiometry(self, ele: libsbml.SBase) -> sp.Expr:
         :param ele:
             reactant or product
         :return:
-            symbolic variable that defines stoichiometry
+            _symbolic variable that defines stoichiometry
         """
         if ele.isSetId():
             sym = self._get_element_initial_assignment(ele.getId())
diff --git a/python/sdist/amici/jax/ode_export.py b/python/sdist/amici/jax/ode_export.py
index d855372b23..f1c6d60a95 100644
--- a/python/sdist/amici/jax/ode_export.py
+++ b/python/sdist/amici/jax/ode_export.py
@@ -21,16 +21,16 @@
 from amici import (
     amiciModulePath,
 )
-from amici.de_model import DEModel
+from amici._symbolic.de_model import DEModel
+from amici._symbolic.sympy_utils import (
+    _monkeypatch_sympy,
+)
 from amici.exporters.sundials.de_export import is_valid_identifier
 from amici.exporters.template import apply_template
 from amici.jax.jaxcodeprinter import AmiciJaxCodePrinter, _jnp_array_str
 from amici.jax.model import JAXModel
 from amici.jax.nn import generate_equinox
 from amici.logging import get_logger, log_execution_time, set_log_level
-from amici.sympy_utils import (
-    _monkeypatch_sympy,
-)
 
 #: python log manager
 logger = get_logger(__name__, logging.ERROR)
@@ -96,7 +96,7 @@ def _jax_variable_ids(model: DEModel, sym_names: tuple[str, ...]) -> dict:
 class ODEExporter:
     """
     The ODEExporter class generates AMICI jax files for a model as
-    defined in symbolic expressions.
+    defined in _symbolic expressions.
 
     :ivar model:
         DE definition
diff --git a/python/tests/test_de_model.py b/python/tests/test_de_model.py
index 4a569af03b..9a78b100d7 100644
--- a/python/tests/test_de_model.py
+++ b/python/tests/test_de_model.py
@@ -1,5 +1,5 @@
 import sympy as sp
-from amici.de_model_components import Event
+from amici._symbolic.de_model_components import Event
 from amici.importers.utils import amici_time_symbol
 from amici.testing import skip_on_valgrind
 
diff --git a/python/tests/test_ode_export.py b/python/tests/test_ode_export.py
index b11425a98f..90153e6967 100644
--- a/python/tests/test_ode_export.py
+++ b/python/tests/test_ode_export.py
@@ -100,7 +100,7 @@ def test_csc_matrix_vector():
 
 @skip_on_valgrind
 def test_match_deriv():
-    from amici.de_model import DERIVATIVE_PATTERN as pat
+    from amici._symbolic.de_model import DERIVATIVE_PATTERN as pat
 
     def check(str, out1, out2):
         match = pat.match(str)
diff --git a/python/tests/test_pysb.py b/python/tests/test_pysb.py
index ed796ebe68..de946f5bd6 100644
--- a/python/tests/test_pysb.py
+++ b/python/tests/test_pysb.py
@@ -17,7 +17,7 @@
 import pytest
 import sympy as sp
 from amici import ParameterScaling, parameter_scaling_from_int_vector
-from amici.de_model_components import Event
+from amici._symbolic.de_model_components import Event
 from amici.gradient_check import check_derivatives
 from amici.importers.pysb import pysb2amici
 from amici.importers.utils import amici_time_symbol
diff --git a/python/tests/test_sciml.py b/python/tests/test_sciml.py
index 24a9bca30a..84af564bfe 100644
--- a/python/tests/test_sciml.py
+++ b/python/tests/test_sciml.py
@@ -375,8 +375,8 @@ class TestProcessHybridizationErrors:
     def mock_de_model(self):
         """Create a mock DEModel instance for testing."""
         import sympy as sp
-        from amici.de_model import DEModel
-        from amici.de_model_components import (
+        from amici._symbolic.de_model import DEModel
+        from amici._symbolic.de_model_components import (
             DifferentialState,
             Expression,
             FreeParameter,
diff --git a/python/tests/test_sympy_utils.py b/python/tests/test_sympy_utils.py
index a868cd8f8a..c60d356868 100644
--- a/python/tests/test_sympy_utils.py
+++ b/python/tests/test_sympy_utils.py
@@ -1,7 +1,7 @@
 """Tests related to the sympy_utils module."""
 
 import sympy as sp
-from amici.sympy_utils import (
+from amici._symbolic.sympy_utils import (
     _custom_pow_eval_derivative,
     _monkeypatched,
     _piecewise_to_minmax,