diff --git a/biobalm/_pint_reachability.py b/biobalm/_pint_reachability.py
index 2bb2aab..7973211 100644
--- a/biobalm/_pint_reachability.py
+++ b/biobalm/_pint_reachability.py
@@ -16,32 +16,19 @@
 
 from biobalm.petri_net_translation import (
     place_to_variable,
-    _optimized_recursive_dnf_generator,
+    optimized_recursive_dnf_generator,
 )
-from biobalm.types import BooleanSpace
-import biobalm
+from biobalm.types import BooleanSpace, SuccessionDiagramConfiguration
 
 if TYPE_CHECKING:
     from biodivine_aeon import Bdd
 
-GOAL_SIZE_LIMIT = 8192
-"""
-Pint is called using command line and can only accept a limited number of arguments.
-This limit is applied when constructing goals to avoid reaching this argument limit.
-
-The limit currently applies to the total number of literals that can be used to
-represent a goal.
-
-The default value was empirically tested as safe on Debian linux, but other operating
-systems may need a different limit to stay safe. Nevertheless, this should not be
-an issue on smaller/simpler networks.
-"""
-
 
-def _pint_reachability(
+def pint_reachability(
     petri_net: DiGraph,
     initial_state: BooleanSpace,
     target_states: Bdd,
+    config: SuccessionDiagramConfiguration,
 ) -> bool:
     """
     Use Pint to check if a given `initial_state` can possibly reach some state
@@ -59,15 +46,17 @@ def _pint_reachability(
     for var, level in initial_state.items():
         pint_model.initial_state[var] = level
 
-    goal = _pint_build_symbolic_goal(target_states)
+    goal = _pint_build_symbolic_goal(target_states, config)
 
-    def failed(_model, _goal):
+    def failed(*_):
         raise RuntimeError("Cannot verify.")
 
     return pint_model.reachability(goal=goal, fallback=failed)  # type: ignore
 
 
-def _pint_build_symbolic_goal(states: Bdd) -> Goal:
+def _pint_build_symbolic_goal(
+    states: Bdd, config: SuccessionDiagramConfiguration
+) -> Goal:
     """
     A helper method which (very explicitly) converts a set of states
     represented through a BDD into a Pint `Goal`.
@@ -80,8 +69,8 @@ def _pint_build_symbolic_goal(states: Bdd) -> Goal:
     assert not states.is_false()
 
     goals: list[Goal] = []
-    limit = GOAL_SIZE_LIMIT
-    for clause in _optimized_recursive_dnf_generator(states):
+    limit = config["pint_goal_size_limit"]
+    for clause in optimized_recursive_dnf_generator(states):
         named_clause = {
             states.__ctx__().get_variable_name(var): int(value)
             for var, value in clause.items()
@@ -92,7 +81,7 @@ def _pint_build_symbolic_goal(states: Bdd) -> Goal:
             # If the goal is too large to be passed as a command line argument,
             # break here and don't continue. This is not ideal but I'm not sure
             # how to fix this other than modifying `pint` itself.
-            if biobalm.succession_diagram.DEBUG:
+            if config["debug"]:
                 print(
                     "WARNING: `pint` goal size limit exceeded. A partial goal is used."
                 )
@@ -127,8 +116,8 @@ def _petri_net_as_automata_network(petri_net: DiGraph) -> str:
         if kind != "transition":
             continue
 
-        predecessors = set(cast(Iterable[str], petri_net.predecessors(transition)))
-        successors = set(cast(Iterable[str], petri_net.successors(transition)))
+        predecessors = set(cast(Iterable[str], petri_net.predecessors(transition)))  # type: ignore
+        successors = set(cast(Iterable[str], petri_net.successors(transition)))  # type: ignore
 
         # The value under modification is the only
         # value that is different between successors and predecessors.
diff --git a/biobalm/_sd_algorithms/expand_attractor_seeds.py b/biobalm/_sd_algorithms/expand_attractor_seeds.py
index 65c4cc4..b0d85dc 100644
--- a/biobalm/_sd_algorithms/expand_attractor_seeds.py
+++ b/biobalm/_sd_algorithms/expand_attractor_seeds.py
@@ -5,8 +5,7 @@
 if TYPE_CHECKING:
     from biobalm.succession_diagram import SuccessionDiagram
 
-import biobalm
-import biobalm.succession_diagram
+from biobalm.types import BooleanSpace
 from biobalm.space_utils import intersect
 from biobalm.trappist_core import compute_fixed_point_reduced_STG
 from biobalm._sd_attractors.attractor_candidates import make_heuristic_retained_set
@@ -25,7 +24,7 @@ def expand_attractor_seeds(sd: SuccessionDiagram, size_limit: int | None = None)
     # motif-avoidant attractors.
     sd.expand_minimal_spaces(size_limit)
 
-    if biobalm.succession_diagram.DEBUG:
+    if sd.config["debug"]:
         print(
             "Minimal trap space expansion finished. Proceeding to attractor expansion."
         )
@@ -82,9 +81,11 @@ def expand_attractor_seeds(sd: SuccessionDiagram, size_limit: int | None = None)
                 intersect(successor_space, child) for child in expanded_motifs
             ]
             avoid = [x for x in avoid_or_none if x is not None]
-            avoid_restricted = []
+            avoid_restricted: list[BooleanSpace] = []
             for x in avoid:
-                y = {var: val for (var, val) in x.items() if var not in successor_space}
+                y: BooleanSpace = {
+                    var: val for (var, val) in x.items() if var not in successor_space
+                }
                 avoid_restricted.append(y)
 
             retained_set = make_heuristic_retained_set(
@@ -105,7 +106,7 @@ def expand_attractor_seeds(sd: SuccessionDiagram, size_limit: int | None = None)
                 successors.pop()
                 continue
 
-            if biobalm.succession_diagram.DEBUG:
+            if sd.config["debug"]:
                 print(
                     f"[{node}] Found successor with new attractor candidate seeds. Expand node {successors[-1]}."
                 )
@@ -114,7 +115,7 @@ def expand_attractor_seeds(sd: SuccessionDiagram, size_limit: int | None = None)
 
         if len(successors) == 0:
             # Everything is done for this `node` and we can continue to the next one.
-            if biobalm.succession_diagram.DEBUG:
+            if sd.config["debug"]:
                 print(f"[{node}] Finished node attractor expansion.")
             continue
 
diff --git a/biobalm/_sd_attractors/attractor_candidates.py b/biobalm/_sd_attractors/attractor_candidates.py
index 05979ec..a0de3c9 100644
--- a/biobalm/_sd_attractors/attractor_candidates.py
+++ b/biobalm/_sd_attractors/attractor_candidates.py
@@ -16,29 +16,15 @@
 
 import random
 import biobalm
-from biodivine_aeon import Bdd, AsynchronousGraph
+from biodivine_aeon import Bdd, AsynchronousGraph, BddVariable
 from biobalm.trappist_core import compute_fixed_point_reduced_STG
 from biobalm.symbolic_utils import state_list_to_bdd, valuation_to_state, state_to_bdd
 
 try:
-    PINT_AVAILABLE = True
-    import biobalm._pint_reachability as pint
+    pint_available = True
+    import biobalm._pint_reachability
 except ModuleNotFoundError:
-    PINT_AVAILABLE = False
-
-
-CANDIDATES_HARD_LIMIT = 100_000
-"""
-If there are more than this amount of attractor candidates, then
-attractor detection will fail with a runtime error. This is mainly to
-avoid out-of-memory errors or crashing clingo.
-"""
-
-RETAINED_SET_OPTIMIZATION_THRESHOLD = 100
-"""
-If there are more than this amount of attractor candidates, we will
-try to optimize the retained set using ASP (if enabled).
-"""
+    pint_available = False
 
 
 def compute_attractor_candidates(
@@ -103,7 +89,7 @@ def compute_attractor_candidates(
         The list of attractor candidate states.
     """
 
-    if biobalm.succession_diagram.DEBUG:
+    if sd.config["debug"]:
         print(f"[{node_id}] Start computing attractor candidates.")
 
     node_data = sd.node_data(node_id)
@@ -111,13 +97,13 @@ def compute_attractor_candidates(
     node_space = node_data["space"]
 
     if len(node_space) == sd.network.variable_count():
-        if biobalm.succession_diagram.DEBUG:
+        if sd.config["debug"]:
             print(f"[{node_id}] > Attractor candidates done: node is a fixed-point.")
         return [node_space]
 
     node_nfvs = sd.node_percolated_nfvs(node_id, compute=True)
 
-    if biobalm.succession_diagram.DEBUG:
+    if sd.config["debug"]:
         root_nfvs = sd.node_percolated_nfvs(sd.root(), compute=True)
         print(
             f"[{node_id}] > Percolated node NFVS contains {len(node_nfvs)} nodes (vs {len(root_nfvs)} in the root)."
@@ -140,7 +126,7 @@ def compute_attractor_candidates(
     # In either case, the space must contain at least one attractor.
     node_is_pseudo_minimal = len(child_motifs_reduced) == 0
 
-    if biobalm.succession_diagram.DEBUG and node_is_pseudo_minimal:
+    if sd.config["debug"] and node_is_pseudo_minimal:
         print(
             f"[{node_id}] > The node has no children (i.e. it is minimal or unexpanded)."
         )
@@ -162,7 +148,7 @@ def compute_attractor_candidates(
         assert not sd.node_is_minimal(node_id)
 
         if not node_is_pseudo_minimal:
-            if biobalm.succession_diagram.DEBUG:
+            if sd.config["debug"]:
                 print(
                     f"[{node_id}] > Attractor candidates done: empty NFVS in a non-minimal space."
                 )
@@ -197,7 +183,7 @@ def compute_attractor_candidates(
     if len(retained_set) == sd.network.variable_count() and node_is_pseudo_minimal:
         # If the retained set describes a fixed point, then only one attractor
         # is present in this space and it must contain the state described by the retained set.
-        if biobalm.succession_diagram.DEBUG:
+        if sd.config["debug"]:
             print(
                 f"[{node_id}] > Singular attractor found through fixed-point retained set. Done."
             )
@@ -208,13 +194,13 @@ def compute_attractor_candidates(
             pn_reduced,
             retained_set,
             avoid_subspaces=child_motifs_reduced,
-            solution_limit=CANDIDATES_HARD_LIMIT,
+            solution_limit=sd.config["attractor_candidates_limit"],
         )
-        if len(candidate_states) == CANDIDATES_HARD_LIMIT:
+        if len(candidate_states) == sd.config["attractor_candidates_limit"]:
             raise RuntimeError(
-                f"Exceeded the maximum amount of attractor candidates ({CANDIDATES_HARD_LIMIT=})."
+                f"Exceeded the maximum amount of attractor candidates ({sd.config['attractor_candidates_limit']}; see `SuccessionDiagramConfiguation.attractor_candidates_limit`)."
             )
-        if biobalm.succession_diagram.DEBUG:
+        if sd.config["debug"]:
             print(
                 f"[{node_id}] Computed {len(candidate_states)} candidate states without retained set optimization."
             )
@@ -223,21 +209,22 @@ def compute_attractor_candidates(
             pn_reduced,
             retained_set,
             avoid_subspaces=child_motifs_reduced,
-            solution_limit=RETAINED_SET_OPTIMIZATION_THRESHOLD,
+            solution_limit=sd.config["retained_set_optimization_threshold"],
         )
 
-        if len(candidate_states) < RETAINED_SET_OPTIMIZATION_THRESHOLD:
+        if len(candidate_states) < sd.config["retained_set_optimization_threshold"]:
             # The candidate set is small. No need to overthink it. The original
             # retained set is probably quite good. However, we can still try to
             # optimize it further if it makes sense.
             if len(candidate_states) > 1 or (
                 not node_is_pseudo_minimal and len(candidate_states) > 0
             ):
-                if biobalm.succession_diagram.DEBUG:
+                if sd.config["debug"]:
                     print(
                         f"[{node_id}] Initial retained set generated {len(candidate_states)} candidates. Optimizing..."
                     )
                 optimized = asp_greedy_retained_set_optimization(
+                    sd,
                     node_id,
                     petri_net=pn_reduced,
                     retained_set=retained_set,
@@ -249,9 +236,9 @@ def compute_attractor_candidates(
         else:
             # There seem to be many candidates, in which case it might be better
             # to optimize the retained set dynamically.
-            if biobalm.succession_diagram.DEBUG:
+            if sd.config["debug"]:
                 print(
-                    f"[{node_id}] Initial retained set generated >{RETAINED_SET_OPTIMIZATION_THRESHOLD} candidates. Regenerate retained set."
+                    f"[{node_id}] Initial retained set generated >{sd.config['retained_set_optimization_threshold']} candidates. Regenerate retained set."
                 )
             retained_set = {}
             candidate_states = []
@@ -261,11 +248,11 @@ def compute_attractor_candidates(
                     pn_reduced,
                     retained_set,
                     avoid_subspaces=child_motifs_reduced,
-                    solution_limit=CANDIDATES_HARD_LIMIT,
+                    solution_limit=sd.config["attractor_candidates_limit"],
                 )
 
                 if len(candidate_states_zero) <= len(candidate_states):
-                    if biobalm.succession_diagram.DEBUG:
+                    if sd.config["debug"]:
                         print(
                             f"[{node_id}] Chosen {var}=0 without increasing candidate count ({len(candidate_states_zero)}). {len(retained_set)}/{len(node_nfvs)} variables chosen."
                         )
@@ -281,15 +268,17 @@ def compute_attractor_candidates(
                 )
 
                 if (
-                    len(candidate_states_zero) == CANDIDATES_HARD_LIMIT
-                    and len(candidate_states_one) == CANDIDATES_HARD_LIMIT
+                    len(candidate_states_zero)
+                    == sd.config["attractor_candidates_limit"]
+                    and len(candidate_states_one)
+                    == sd.config["attractor_candidates_limit"]
                 ):
                     raise RuntimeError(
-                        f"Exceeded the maximum amount of attractor candidates ({CANDIDATES_HARD_LIMIT=})."
+                        f"Exceeded the maximum amount of attractor candidates ({sd.config['attractor_candidates_limit']}; see `SuccessionDiagramConfiguation.attractor_candidates_limit`)."
                     )
 
                 if len(candidate_states_one) <= len(candidate_states):
-                    if biobalm.succession_diagram.DEBUG:
+                    if sd.config["debug"]:
                         print(
                             f"[{node_id}] Chosen {var}=1 without increasing candidate count ({len(candidate_states_one)}). {len(retained_set)}/{len(node_nfvs)} variables chosen."
                         )
@@ -297,14 +286,14 @@ def compute_attractor_candidates(
                     continue
 
                 if len(candidate_states_zero) < len(candidate_states_one):
-                    if biobalm.succession_diagram.DEBUG:
+                    if sd.config["debug"]:
                         print(
                             f"[{node_id}] Chosen {var}=0 with better candidate count ({len(candidate_states_zero)}). {len(retained_set)}/{len(node_nfvs)} variables chosen."
                         )
                     candidate_states = candidate_states_zero
                     retained_set[var] = 0
                 else:
-                    if biobalm.succession_diagram.DEBUG:
+                    if sd.config["debug"]:
                         print(
                             f"[{node_id}] Chosen {var}=1 with better candidate count ({len(candidate_states_one)}). {len(retained_set)}/{len(node_nfvs)} variables chosen."
                         )
@@ -313,9 +302,10 @@ def compute_attractor_candidates(
 
                 # At this point, we know the candidate count increased and so we should
                 # try to bring it back down.
-                if biobalm.succession_diagram.DEBUG:
+                if sd.config["debug"]:
                     print(f"[{node_id}] Optimizing partial retained set...")
                 optimized = asp_greedy_retained_set_optimization(
+                    sd,
                     node_id,
                     petri_net=pn_reduced,
                     retained_set=retained_set,
@@ -327,23 +317,23 @@ def compute_attractor_candidates(
 
     # Terminate if done.
     if len(candidate_states) == 0:
-        if biobalm.succession_diagram.DEBUG:
+        if sd.config["debug"]:
             print(f"[{node_id}] > Initial candidate set empty. Done.")
         return []
     if node_is_pseudo_minimal and len(candidate_states) == 1:
-        if biobalm.succession_diagram.DEBUG:
+        if sd.config["debug"]:
             print(
                 f"[{node_id}] > Single candidate found in (pseudo) minimal trap space. Done."
             )
         return [candidate_states[0] | node_space]
 
-    if biobalm.succession_diagram.DEBUG:
+    if sd.config["debug"]:
         print(
             f"[{node_id}] > Attractor candidates from retained set: {len(candidate_states)}."
         )
 
     if simulation_minification:
-        if biobalm.succession_diagram.DEBUG:
+        if sd.config["debug"]:
             print(f"[{node_id}] Start simulation minification...")
         avoid_children_symbolic = state_list_to_bdd(
             graph_reduced.symbolic_context(), child_motifs_reduced
@@ -354,11 +344,12 @@ def compute_attractor_candidates(
         # cannot reduce any further states, we are done.
         iterations = 1024
         while len(candidate_states) > 0:
-            if biobalm.succession_diagram.DEBUG:
+            if sd.config["debug"]:
                 print(
                     f"[{node_id}] > Start simulation with {len(candidate_states)} states and simulation limit {iterations}."
                 )
             reduced = run_simulation_minification(
+                sd,
                 node_id,
                 graph_reduced,
                 candidate_states,
@@ -377,25 +368,25 @@ def compute_attractor_candidates(
             if len(candidate_states) == 1 and avoid_children_symbolic.is_false():
                 break
 
-        if biobalm.succession_diagram.DEBUG:
+        if sd.config["debug"]:
             print(f"[{node_id}] > Candidates after simulation: {len(candidate_states)}")
 
     # Terminate if done.
     if len(candidate_states) == 0:
-        if biobalm.succession_diagram.DEBUG:
+        if sd.config["debug"]:
             print(f"[{node_id}] > Candidate set empty. Done.")
         return []
     if node_is_pseudo_minimal and len(candidate_states) == 1:
-        if biobalm.succession_diagram.DEBUG:
+        if sd.config["debug"]:
             print(
                 f"[{node_id}] > Single candidate found in (pseudo) minimal trap space. Done."
             )
         return [candidate_states[0] | node_space]
 
-    if pint_minification and not PINT_AVAILABLE:
+    if pint_minification and not pint_available:
         print("WARNING: Using `pint`, but `pint` is not installed. Skipping.")
-    elif pint_minification and PINT_AVAILABLE:
-        if biobalm.succession_diagram.DEBUG:
+    elif pint_minification and pint_available:
+        if sd.config["debug"]:
             print(f"[{node_id}] Start `pint` minification...")
 
         children_bdd = state_list_to_bdd(
@@ -406,7 +397,7 @@ def compute_attractor_candidates(
         )
         avoid_bdd = children_bdd.l_or(candidates_bdd)
 
-        filtered_states = []
+        filtered_states: list[BooleanSpace] = []
         for i, state in enumerate(candidate_states):
             state_bdd = state_to_bdd(graph_reduced.symbolic_context(), state)
 
@@ -414,7 +405,9 @@ def compute_attractor_candidates(
 
             keep = True
             try:
-                if pint._pint_reachability(pn_reduced, state, avoid_bdd):
+                if biobalm._pint_reachability.pint_reachability(
+                    pn_reduced, state, avoid_bdd, sd.config
+                ):
                     keep = False
             except RuntimeError as e:
                 assert str(e) == "Cannot verify."
@@ -423,14 +416,14 @@ def compute_attractor_candidates(
                 avoid_bdd = avoid_bdd.l_or(state_bdd)
                 filtered_states.append(state)
 
-            if biobalm.succession_diagram.DEBUG:
+            if sd.config["debug"]:
                 print(
                     f"[{node_id}] > `pint` {i + 1}/{len(candidate_states)}: eliminated: {not keep}, retained: {len(filtered_states)}."
                 )
 
         candidate_states = filtered_states
 
-        if biobalm.succession_diagram.DEBUG:
+        if sd.config["debug"]:
             print(f"[{node_id}] > Candidates after `pint`: {len(candidate_states)}")
 
     # Finally, we need to augment each candidate state with the
@@ -442,6 +435,7 @@ def compute_attractor_candidates(
 
 
 def run_simulation_minification(
+    sd: SuccessionDiagram,
     node_id: int,
     graph: AsynchronousGraph,
     candidate_states: list[BooleanSpace],
@@ -484,7 +478,7 @@ def run_simulation_minification(
     # associate them with the corresponding network's update functions.
     symbolic_ctx = graph.symbolic_context()
     network_vars = graph.network_variables()
-    symbolic_vars = []
+    symbolic_vars: list[BddVariable] = []
     for var in network_vars:
         s_var = symbolic_ctx.find_network_bdd_variable(var)
         assert s_var is not None
@@ -502,7 +496,7 @@ def run_simulation_minification(
         filtered_candidates: list[BooleanSpace] = []
 
         for i, state in enumerate(candidate_states):
-            if i % 100 == 99 and biobalm.succession_diagram.DEBUG:
+            if i % 100 == 99 and sd.config["debug"]:
                 print(
                     f"[{node_id}] > Simulation progress: {i + 1}/{len(candidate_states)}"
                 )
@@ -566,12 +560,12 @@ def run_simulation_minification(
         # and (b) we can stop with one candidate instead of zero.
 
         candidates_bdd = state_list_to_bdd(symbolic_ctx, candidate_states)
-        printed = set()
+        printed: set[int] = set()
         for i in range(max_iterations):
             progress = int((i * len(candidate_states)) / max_iterations)
             if (
                 progress % 100 == 99
-                and biobalm.succession_diagram.DEBUG
+                and sd.config["debug"]
                 and (progress not in printed)
             ):
                 printed.add(progress)
@@ -613,6 +607,7 @@ def run_simulation_minification(
 
 
 def asp_greedy_retained_set_optimization(
+    sd: SuccessionDiagram,
     node_id: int,
     petri_net: DiGraph,
     retained_set: BooleanSpace,
@@ -674,7 +669,7 @@ def asp_greedy_retained_set_optimization(
                 retained_set = retained_set_2
                 candidate_states = candidate_states_2
                 done = False
-                if biobalm.succession_diagram.DEBUG:
+                if sd.config["debug"]:
                     print(
                         f"[{node_id}] > Candidate states optimized to {len(candidate_states)}."
                     )
@@ -723,7 +718,7 @@ def make_heuristic_retained_set(
         A :class:`BooleanSpace<biobalm.types.BooleanSpace>` object describing the
         retained set.
     """
-    retained_set = {}
+    retained_set: BooleanSpace = {}
 
     # First, if there are any child spaces present, we extend the retained set
     # with the values from the one that has the least amount of fixed variables
diff --git a/biobalm/_sd_attractors/attractor_symbolic.py b/biobalm/_sd_attractors/attractor_symbolic.py
index e713830..12cd0bb 100644
--- a/biobalm/_sd_attractors/attractor_symbolic.py
+++ b/biobalm/_sd_attractors/attractor_symbolic.py
@@ -4,10 +4,9 @@
 
 if TYPE_CHECKING:
     from biobalm.succession_diagram import SuccessionDiagram
-    from biodivine_aeon import VariableId, VertexSet
+    from biodivine_aeon import VariableId
 
-import biobalm
-from biodivine_aeon import AsynchronousGraph, ColoredVertexSet
+from biodivine_aeon import AsynchronousGraph, ColoredVertexSet, VertexSet
 from biobalm.symbolic_utils import state_list_to_bdd
 from biobalm.types import BooleanSpace
 
@@ -50,10 +49,12 @@ def compute_attractors_symbolic(
 
     # Ignore variables that are already fixed in the nodes space.
     # These are not used in the reduced network.
-    candidate_states_reduced = []
+    candidate_states_reduced: list[BooleanSpace] = []
     for candidate in candidate_states:
-        candidate = {k: v for (k, v) in candidate.items() if k not in node_space}
-        candidate_states_reduced.append(candidate)
+        candidate_reduced: BooleanSpace = {
+            k: v for (k, v) in candidate.items() if k not in node_space
+        }
+        candidate_states_reduced.append(candidate_reduced)
 
     child_motifs_reduced = []
     if node_data["expanded"]:
@@ -70,13 +71,13 @@ def compute_attractors_symbolic(
 
     avoid = candidate_set.union(children_set)
 
-    if biobalm.succession_diagram.DEBUG:
+    if sd.config["debug"]:
         print(
             f"[{node_id}] > Start symbolic seed state identification with {len(candidate_states)} candidates and {avoid} avoid states."
         )
 
     seeds: list[BooleanSpace] = []
-    sets = []
+    sets: list[ColoredVertexSet] = []
     for i, candidate in enumerate(candidate_states_reduced):
         is_last = i == len(candidate_states_reduced) - 1
         is_minimal = len(child_motifs_reduced) == 0
@@ -84,7 +85,7 @@ def compute_attractors_symbolic(
             # This node is pseudo-minimal, so it must contain at least
             # one attractor seed. If we only care about the attractor seeds,
             # we can thus return the last candidate without checking it.
-            if biobalm.succession_diagram.DEBUG:
+            if sd.config["debug"]:
                 print(
                     f"[{node_id}] > Single seed remaining in a (pseduo) minimal space. Done."
                 )
@@ -93,7 +94,7 @@ def compute_attractors_symbolic(
 
         avoid = avoid.minus(candidate_singleton)
 
-        closure = symbolic_attractor_test(node_id, graph_reduced, candidate, avoid)
+        closure = symbolic_attractor_test(sd, node_id, graph_reduced, candidate, avoid)
 
         if closure is None:
             # This candidate can reach someone else in the candidate set,
@@ -107,11 +108,11 @@ def compute_attractors_symbolic(
         seeds.append(candidate | node_space)
         sets.append(closure)
 
-    if biobalm.succession_diagram.DEBUG:
+    if sd.config["debug"]:
         print(f"[{node_id}] > Finished identification with {len(seeds)} seed states.")
 
     space_symbolic = sd.symbolic.mk_subspace(node_space).vertices()
-    sets_converted = []
+    sets_converted: list[VertexSet] = []
     for s in sets:
         # Extend the attractor set with fixed nodes from the node space.
         vertices = s.vertices()
@@ -123,6 +124,7 @@ def compute_attractors_symbolic(
 
 
 def symbolic_attractor_test(
+    sd: SuccessionDiagram,
     node_id: int,
     graph: AsynchronousGraph,
     pivot: BooleanSpace,
@@ -173,7 +175,7 @@ def symbolic_attractor_test(
     ]
     other_vars = sort_variable_list(other_vars)
 
-    if biobalm.succession_diagram.DEBUG:
+    if sd.config["debug"]:
         print(
             f"[{node_id}] > Start symbolic reachability with {len(conflict_vars)} conflict variables and {len(other_vars)} other variables."
         )
@@ -186,7 +188,7 @@ def symbolic_attractor_test(
         saturation_done = False
         while not saturation_done:
             if avoid is not None and not avoid.intersect(reach_set).is_empty():
-                if biobalm.succession_diagram.DEBUG:
+                if sd.config["debug"]:
                     print(f"[{node_id}] > Discovered avoid state. Done.")
                 return None
 
@@ -196,10 +198,7 @@ def symbolic_attractor_test(
                 if not successors.is_empty():
                     reach_set = reach_set.union(successors)
                     saturation_done = False
-                    if (
-                        reach_set.symbolic_size() > 100_000
-                        and biobalm.succession_diagram.DEBUG
-                    ):
+                    if reach_set.symbolic_size() > 100_000 and sd.config["debug"]:
                         print(
                             f"[{node_id}] > Saturation({len(saturated_vars)}) Expanded reach_set: {reach_set}"
                         )
@@ -227,13 +226,13 @@ def symbolic_attractor_test(
                 saturated_vars.append(var)
                 saturated_vars = sort_variable_list(saturated_vars)
 
-                if biobalm.succession_diagram.DEBUG:
+                if sd.config["debug"]:
                     print(
                         f"[{node_id}] > Saturation({len(saturated_vars)}) Added saturation variable. {len(conflict_vars)} conflict and {len(other_vars)} other variables remaining."
                     )
                 break
 
-    if biobalm.succession_diagram.DEBUG:
+    if sd.config["debug"]:
         print(f"[{node_id}] > Reachability completed with {reach_set}.")
 
     return reach_set
diff --git a/biobalm/petri_net_translation.py b/biobalm/petri_net_translation.py
index 261daee..0c7cda4 100644
--- a/biobalm/petri_net_translation.py
+++ b/biobalm/petri_net_translation.py
@@ -359,14 +359,14 @@ def network_to_petrinet(
     return pn
 
 
-def _optimized_recursive_dnf_generator(
+def optimized_recursive_dnf_generator(
     bdd: Bdd,
 ) -> Generator[BddPartialValuation, None, None]:
     """
     Yields a generator of `BddPartialValuation` objects, similar to
-    `bdd.clause_iterator`, but uses a recursive optimization strategy to return
-    a smaller result than the default method `Bdd` clause sequence. Note that
-    this is still not the "optimal" DNF, but is often close enough.
+    `Bdd.clause_iterator` in AEON, but uses a recursive optimization strategy
+    to return a smaller result than the default method `Bdd` clause sequence.
+    Note that this is still not the "optimal" DNF, but is often close enough.
 
     This is technically slower for BDDs that already have a small clause count,
     but can be much better in the long-term when the clause count is high.
@@ -398,11 +398,11 @@ def _optimized_recursive_dnf_generator(
             best_size = t_size + f_size
             best_var = var
 
-    for t_val in _optimized_recursive_dnf_generator(bdd.r_restrict({best_var: True})):
+    for t_val in optimized_recursive_dnf_generator(bdd.r_restrict({best_var: True})):
         t_val[best_var] = True
         yield t_val
 
-    for f_val in _optimized_recursive_dnf_generator(bdd.r_restrict({best_var: False})):
+    for f_val in optimized_recursive_dnf_generator(bdd.r_restrict({best_var: False})):
         f_val[best_var] = False
         yield f_val
 
@@ -421,7 +421,7 @@ def _create_transitions(
     """
     dir_str = "up" if go_up else "down"
     total = 0
-    for t_id, implicant in enumerate(_optimized_recursive_dnf_generator(implicant_bdd)):
+    for t_id, implicant in enumerate(optimized_recursive_dnf_generator(implicant_bdd)):
         total += 1
         t_name = f"tr_{var_name}_{dir_str}_{t_id + 1}"
         pn.add_node(  # type: ignore[reportUnknownMemberType]
diff --git a/biobalm/succession_diagram.py b/biobalm/succession_diagram.py
index 4f12d79..f370564 100644
--- a/biobalm/succession_diagram.py
+++ b/biobalm/succession_diagram.py
@@ -32,12 +32,12 @@
 )
 from biobalm.space_utils import percolate_network, percolate_space, space_unique_key
 from biobalm.trappist_core import trappist
-from biobalm.types import BooleanSpace, NodeData, SuccessionDiagramState
-
-DEBUG: bool = False
-"""
-If `True`, print debug and progress messages.
-"""
+from biobalm.types import (
+    BooleanSpace,
+    NodeData,
+    SuccessionDiagramState,
+    SuccessionDiagramConfiguration,
+)
 
 
 class SuccessionDiagram:
@@ -79,26 +79,6 @@ class SuccessionDiagram:
     <BLANKLINE>
     """
 
-    NFVS_NODE_THRESHOLD: int = 2_000
-    """
-    If the number of nodes in the (redunced) network is greater than this
-    threshold, we find a feedback vertex set (without consideration of cycle
-    sign) for reachability preprocessing. There is a trade-off between the speed
-    gains from a smaller node set to consider and the cost of determining which
-    FVS nodes only intersect negative cycles to find an NFVS subset. Typically,
-    for smaller networks, the trade-off is in favor of computing a smaller NFVS.
-    """
-
-    MAX_MOTIFS_PER_NODE = 100_000
-    """
-    An artificial limit on the number of stable motifs that can be
-    considered during the expansion a single node.
-
-    This limit is in place mainly to avoid surprising out of memory errors,
-    because currently there is no logging mechanism that would report the
-    number of stable motifs gradually.
-    """
-
     __slots__ = (
         "network",
         "symbolic",
@@ -106,9 +86,18 @@ class SuccessionDiagram:
         "nfvs",
         "dag",
         "node_indices",
+        "config",
     )
 
-    def __init__(self, network: BooleanNetwork):
+    def __init__(
+        self,
+        network: BooleanNetwork,
+        config: SuccessionDiagramConfiguration | None = None,
+    ):
+        if config is None:
+            config = SuccessionDiagram.default_config()
+        self.config = config
+
         # Original Boolean network.
         self.network: BooleanNetwork = cleanup_network(network)
         """
@@ -125,7 +114,7 @@ def __init__(self, network: BooleanNetwork):
         The Petri net representation of the network (see :mod:`petri_net_translation<biobalm.petri_net_translation>`).
         """
 
-        if DEBUG:
+        if self.config["debug"]:
             print(
                 f"Generated global Petri net with {len(self.petri_net.nodes)} nodes and {len(self.petri_net.edges)} edges."
             )
@@ -157,6 +146,7 @@ def __getstate__(self) -> SuccessionDiagramState:
             "nfvs": self.nfvs,
             "dag": self.dag,
             "node_indices": self.node_indices,
+            "config": self.config,
         }
 
     def __setstate__(self, state: SuccessionDiagramState):
@@ -167,6 +157,7 @@ def __setstate__(self, state: SuccessionDiagramState):
         self.nfvs = state["nfvs"]
         self.dag = state["dag"]
         self.node_indices = state["node_indices"]
+        self.config = state["config"]
 
     def __len__(self) -> int:
         """
@@ -174,6 +165,17 @@ def __len__(self) -> int:
         """
         return self.dag.number_of_nodes()
 
+    @staticmethod
+    def default_config() -> SuccessionDiagramConfiguration:
+        return {
+            "debug": False,
+            "max_motifs_per_node": 100_000,
+            "nfvs_size_threshold": 2_000,
+            "pint_goal_size_limit": 8_192,
+            "attractor_candidates_limit": 100_000,
+            "retained_set_optimization_threshold": 100,
+        }
+
     @staticmethod
     def from_rules(
         rules: str,
@@ -614,7 +616,7 @@ def reclaim_node_data(self):
         """
 
         for node_id in self.node_ids():
-            data = self.node_data[node_id]
+            data = self.node_data(node_id)
             data["percolated_network"] = None
             data["percolated_petri_net"] = None
             data["percolated_nfvs"] = None
@@ -910,10 +912,8 @@ def node_percolated_nfvs(self, node_id: int, compute: bool = False) -> list[str]
 
         if node["percolated_nfvs"] is None:
             percolated_network = self.node_percolated_network(node_id, compute)
-            if (
-                percolated_network.variable_count()
-                < SuccessionDiagram.NFVS_NODE_THRESHOLD
-            ):
+            percolated_size = percolated_network.variable_count()
+            if percolated_size < self.config["nfvs_size_threshold"]:
                 # Computing the *negative* variant of the FVS is surprisingly costly.
                 # Hence it mostly makes sense for the smaller networks only.
                 nfvs = feedback_vertex_set(percolated_network, parity="negative")
@@ -967,7 +967,7 @@ def node_percolated_network(
             network = percolate_network(
                 self.network, node_space, self.symbolic, remove_constants=True
             )
-            if DEBUG:
+            if self.config["debug"]:
                 print(
                     f"[{node_id}] Computed percolated network with {network.variable_count()} variables (vs {self.network.variable_count()})."
                 )
@@ -1032,7 +1032,7 @@ def node_percolated_petri_net(
 
             percolated_pn = restrict_petrinet_to_subspace(base_pn, percolate_space)
 
-            if DEBUG:
+            if self.config["debug"]:
                 print(
                     f"[{node_id}] Generated Petri net restriction with {len(percolated_pn.nodes)} nodes and {len(percolated_pn.edges)} edges."
                 )
@@ -1298,8 +1298,6 @@ def _expand_one_node(self, node_id: int):
         if node["expanded"]:
             return
 
-        node["expanded"] = True
-
         # If the node had any attractor data computed as unexpanded, these are
         # no longer valid and need to be erased.
         node["attractor_seeds"] = None
@@ -1308,7 +1306,7 @@ def _expand_one_node(self, node_id: int):
 
         current_space = node["space"]
 
-        if DEBUG:
+        if self.config["debug"]:
             print(
                 f"[{node_id}] Expanding: {len(self.node_data(node_id)['space'])} fixed vars."
             )
@@ -1316,8 +1314,9 @@ def _expand_one_node(self, node_id: int):
         if len(current_space) == self.network.variable_count():
             # This node is a fixed-point. Trappist would just
             # return this fixed-point again. No need to continue.
-            if DEBUG:
+            if self.config["debug"]:
                 print(f"[{node_id}] Found fixed-point: {current_space}.")
+            node["expanded"] = True
             return
 
         # We use the non-propagated Petri net for backwards-compatibility reasons here.
@@ -1338,7 +1337,7 @@ def _expand_one_node(self, node_id: int):
                 pn,
                 problem="max",
                 optimize_source_variables=source_nodes,
-                solution_limit=SuccessionDiagram.MAX_MOTIFS_PER_NODE,
+                solution_limit=self.config["max_motifs_per_node"],
             )
             sub_spaces = [(s | current_space) for s in partial_sub_spaces]
         else:
@@ -1349,7 +1348,7 @@ def _expand_one_node(self, node_id: int):
                 problem="max",
                 ensure_subspace=current_space,
                 optimize_source_variables=source_nodes,
-                solution_limit=SuccessionDiagram.MAX_MOTIFS_PER_NODE,
+                solution_limit=self.config["max_motifs_per_node"],
             )
 
         # Release the Petri net once the sub_spaces are computed.
@@ -1358,9 +1357,9 @@ def _expand_one_node(self, node_id: int):
         # in memory.
         node["percolated_petri_net"] = None
 
-        if len(sub_spaces) == SuccessionDiagram.MAX_MOTIFS_PER_NODE:
+        if len(sub_spaces) == self.config["max_motifs_per_node"]:
             raise RuntimeError(
-                f"Exceeded the maximum amount of stable motifs per node (SuccessionDiagram.MAX_MOTIFS_PER_NODE={SuccessionDiagram.MAX_MOTIFS_PER_NODE})."
+                f"Exceeded the maximum amount of stable motifs per node ({self.config['max_motifs_per_node']}; see `SuccessionDiagramConfiguration.max_motifs_per_node`)."
             )
 
         # Sort the spaces based on a unique key in case trappist is not always
@@ -1370,19 +1369,23 @@ def _expand_one_node(self, node_id: int):
         )
 
         if len(sub_spaces) == 0:
-            if DEBUG:
+            if self.config["debug"]:
                 print(f"[{node_id}] Found minimum trap space: {current_space}.")
+            node["expanded"] = True
             return
 
-        if DEBUG:
+        if self.config["debug"]:
             print(f"[{node_id}] Found sub-spaces: {len(sub_spaces)}")
 
         for sub_space in sub_spaces:
             child_id = self._ensure_node(node_id, sub_space)
 
-            if DEBUG:
+            if self.config["debug"]:
                 print(f"[{node_id}] Created edge into node {child_id}.")
 
+        # If everything else worked out, we can mark the node as expanded.
+        node["expanded"] = True
+
     def _ensure_node(self, parent_id: int | None, stable_motif: BooleanSpace) -> int:
         """
         Internal method that ensures the provided node is present in this
diff --git a/biobalm/types.py b/biobalm/types.py
index 9f5692f..01eb09c 100644
--- a/biobalm/types.py
+++ b/biobalm/types.py
@@ -45,6 +45,11 @@ class SuccessionDiagramState(TypedDict):
     diagram (see :func:`biobalm.space_utils.space_unique_key`).
     """
 
+    config: SuccessionDiagramConfiguration
+    """
+    "Global" configuration of a succession diagram.
+    """
+
 
 class NodeData(TypedDict):
     """
@@ -153,3 +158,73 @@ class NodeData(TypedDict):
     If `None`, these have not been computed, and the node may or may not
     have associated attractors.
     """
+
+
+class SuccessionDiagramConfiguration(TypedDict):
+    """
+    Describes the configuration options of a `SuccessionDiagram`.
+
+    Use :meth:`SuccessionDiagram.default_config` to create a
+    configuration dictionary pre-populated with default values.
+    """
+
+    debug: bool
+    """
+    If `True`, the `SuccessionDiagram` will print messages
+    describing the progress of the running operations.
+
+    [Default: False]
+    """
+
+    max_motifs_per_node: int
+    """
+    Limit on the number of stable motifs explored for one node of a succession
+    diagram. If this limit is exceeded during node expansion, a `RuntimeError`
+    is raised and the node remains unexpanded.
+
+    This limit is in place mainly to avoid surprising out of memory errors,
+    because currently there is no logging mechanism that would report the
+    number of stable motifs gradually.
+
+    [Default: 100_000]
+    """
+
+    nfvs_size_threshold: int
+    """
+    For networks larger than this threshold, we only run FVS detection
+    instead of NFVS detection. This is still correct, but can produce
+    a larger node set.
+
+
+    There is a trade-off between the speed gains from a smaller node set
+    to consider and the cost of determining which FVS nodes only
+    intersect negative cycles to find an NFVS subset. Typically,
+    for smaller networks, the trade-off is in favor of
+    computing a smaller NFVS.
+    """
+
+    pint_goal_size_limit: int
+    """
+    Pint is called using command line and can only accept a limited number of arguments.
+    This limit is applied when constructing goals to avoid reaching this argument limit.
+
+    The limit currently applies to the total number of literals that can be used to
+    represent a goal.
+
+    The default value was empirically tested as safe on Debian linux, but other operating
+    systems may need a different limit to stay safe. Nevertheless, this should not be
+    an issue on smaller/simpler networks.
+    """
+
+    attractor_candidates_limit: int
+    """
+    If more than `attractor_candidates_limit` states are produced during the
+    attractor detection process, then the process fails with a `RuntimeError`.
+    This is mainly to avoid out-of-memory errors or crashing `clingo`.
+    """
+
+    retained_set_optimization_threshold: int
+    """
+    If there are more than this amount of attractor candidates, the attractor
+    detection process will try to optimize the retained set using ASP (if enabled).
+    """