`` grouped by the number of acetylated lysines (of K5, K8, K12, K16) are:
+
+- zero acetylated lysines: ``x_0ac``
+- one acetylated lysine: ``x_k05``, ``x_k08``, ``x_k12``, ``x_k16``
+- two acetylated lysines: ``x_k05k08``, ``x_k05k12``, ``x_k05k16``, ``x_k08k12``, ``x_k08k16``, ``x_k12k16``
+- three acetylated lysines: ``x_k05k08k12``, ``x_k05k08k16``, ``x_k05k12k16``, ``x_k08k12k16``
+- four acetylated lysines: ``x_4ac``
+
+The 32 acetylation reactions
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Each acetylation reaction converts one motif :math:`p` into another motif :math:`q`, :math:`p \to q`, at rate :math:`a_{p\to q}\, a_b\, x_p`,
+where :math:`a_b` is the basal (shared) acetylation rate and :math:`a_{p\to q} \equiv` ``a__`` is a factor that represents the research question. This factor is either fixed to one for reactions that occur with the shared basal rate constant, or the factor is estimated, to change the rate constant to be motif-specific.
+Each
+acetylation reaction also has an associated deacetylation reaction :math:`q \to p` with rate
+:math:`d\, x_q`, where :math:`d \equiv` ``da_b`` is a single, shared basal deacetylation rate constant (fixed to ``1``).
+
+The acetylation reactions are explicitly:
+
+.. code-block:: text
+
+ x_0ac -> x_k05 (rate a_0ac_k05 * a_b * x_0ac)
+ x_0ac -> x_k08 (rate a_0ac_k08 * a_b * x_0ac)
+ x_0ac -> x_k12 (rate a_0ac_k12 * a_b * x_0ac)
+ x_0ac -> x_k16 (rate a_0ac_k16 * a_b * x_0ac)
+ x_k05 -> x_k05k08 (rate a_k05_k05k08 * a_b * x_k05)
+ x_k05 -> x_k05k12 (rate a_k05_k05k12 * a_b * x_k05)
+ x_k05 -> x_k05k16 (rate a_k05_k05k16 * a_b * x_k05)
+ x_k08 -> x_k05k08 (rate a_k08_k05k08 * a_b * x_k08)
+ x_k08 -> x_k08k12 (rate a_k08_k08k12 * a_b * x_k08)
+ x_k08 -> x_k08k16 (rate a_k08_k08k16 * a_b * x_k08)
+ x_k12 -> x_k05k12 (rate a_k12_k05k12 * a_b * x_k12)
+ x_k12 -> x_k08k12 (rate a_k12_k08k12 * a_b * x_k12)
+ x_k12 -> x_k12k16 (rate a_k12_k12k16 * a_b * x_k12)
+ x_k16 -> x_k05k16 (rate a_k16_k05k16 * a_b * x_k16)
+ x_k16 -> x_k08k16 (rate a_k16_k08k16 * a_b * x_k16)
+ x_k16 -> x_k12k16 (rate a_k16_k12k16 * a_b * x_k16)
+ x_k05k08 -> x_k05k08k12 (rate a_k05k08_k05k08k12 * a_b * x_k05k08)
+ x_k05k08 -> x_k05k08k16 (rate a_k05k08_k05k08k16 * a_b * x_k05k08)
+ x_k05k12 -> x_k05k08k12 (rate a_k05k12_k05k08k12 * a_b * x_k05k12)
+ x_k05k12 -> x_k05k12k16 (rate a_k05k12_k05k12k16 * a_b * x_k05k12)
+ x_k05k16 -> x_k05k08k16 (rate a_k05k16_k05k08k16 * a_b * x_k05k16)
+ x_k05k16 -> x_k05k12k16 (rate a_k05k16_k05k12k16 * a_b * x_k05k16)
+ x_k08k12 -> x_k05k08k12 (rate a_k08k12_k05k08k12 * a_b * x_k08k12)
+ x_k08k12 -> x_k08k12k16 (rate a_k08k12_k08k12k16 * a_b * x_k08k12)
+ x_k08k16 -> x_k05k08k16 (rate a_k08k16_k05k08k16 * a_b * x_k08k16)
+ x_k08k16 -> x_k08k12k16 (rate a_k08k16_k08k12k16 * a_b * x_k08k16)
+ x_k12k16 -> x_k05k12k16 (rate a_k12k16_k05k12k16 * a_b * x_k12k16)
+ x_k12k16 -> x_k08k12k16 (rate a_k12k16_k08k12k16 * a_b * x_k12k16)
+ x_k05k08k12 -> x_4ac (rate a_k05k08k12_4ac * a_b * x_k05k08k12)
+ x_k05k08k16 -> x_4ac (rate a_k05k08k16_4ac * a_b * x_k05k08k16)
+ x_k05k12k16 -> x_4ac (rate a_k05k12k16_4ac * a_b * x_k05k12k16)
+ x_k08k12k16 -> x_4ac (rate a_k08k12k16_4ac * a_b * x_k08k12k16)
+
+The 32 reverse deacetylations occur at rate e.g. ``da_b * x_k05`` for ``k05 -> 0ac``. These are not listed explicitly but are present for every acetylation reaction above.
+
+The ordinary differential equations
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Based on these reactions and rates, the ODE system is therefore
+
+.. math::
+
+ \frac{\mathrm{d}x_{\mathrm{0ac}}}{\mathrm{d}t} &= d\left(x_{\mathrm{k05}} + x_{\mathrm{k08}} + x_{\mathrm{k12}} + x_{\mathrm{k16}}\right) - a_b\left(a_{\mathrm{0ac}\to\mathrm{k05}} + a_{\mathrm{0ac}\to\mathrm{k08}} + a_{\mathrm{0ac}\to\mathrm{k12}} + a_{\mathrm{0ac}\to\mathrm{k16}}\right)x_{\mathrm{0ac}} \\
+ \frac{\mathrm{d}x_{\mathrm{k05}}}{\mathrm{d}t} &= a_b\,a_{\mathrm{0ac}\to\mathrm{k05}} x_{\mathrm{0ac}} + d\left(x_{\mathrm{k05k08}} + x_{\mathrm{k05k12}} + x_{\mathrm{k05k16}}\right) - a_b\left(a_{\mathrm{k05}\to\mathrm{k05k08}} + a_{\mathrm{k05}\to\mathrm{k05k12}} + a_{\mathrm{k05}\to\mathrm{k05k16}}\right)x_{\mathrm{k05}} - d\,x_{\mathrm{k05}} \\
+ \frac{\mathrm{d}x_{\mathrm{k08}}}{\mathrm{d}t} &= a_b\,a_{\mathrm{0ac}\to\mathrm{k08}} x_{\mathrm{0ac}} + d\left(x_{\mathrm{k05k08}} + x_{\mathrm{k08k12}} + x_{\mathrm{k08k16}}\right) - a_b\left(a_{\mathrm{k08}\to\mathrm{k05k08}} + a_{\mathrm{k08}\to\mathrm{k08k12}} + a_{\mathrm{k08}\to\mathrm{k08k16}}\right)x_{\mathrm{k08}} - d\,x_{\mathrm{k08}} \\
+ \frac{\mathrm{d}x_{\mathrm{k12}}}{\mathrm{d}t} &= a_b\,a_{\mathrm{0ac}\to\mathrm{k12}} x_{\mathrm{0ac}} + d\left(x_{\mathrm{k05k12}} + x_{\mathrm{k08k12}} + x_{\mathrm{k12k16}}\right) - a_b\left(a_{\mathrm{k12}\to\mathrm{k05k12}} + a_{\mathrm{k12}\to\mathrm{k08k12}} + a_{\mathrm{k12}\to\mathrm{k12k16}}\right)x_{\mathrm{k12}} - d\,x_{\mathrm{k12}} \\
+ \frac{\mathrm{d}x_{\mathrm{k16}}}{\mathrm{d}t} &= a_b\,a_{\mathrm{0ac}\to\mathrm{k16}} x_{\mathrm{0ac}} + d\left(x_{\mathrm{k05k16}} + x_{\mathrm{k08k16}} + x_{\mathrm{k12k16}}\right) - a_b\left(a_{\mathrm{k16}\to\mathrm{k05k16}} + a_{\mathrm{k16}\to\mathrm{k08k16}} + a_{\mathrm{k16}\to\mathrm{k12k16}}\right)x_{\mathrm{k16}} - d\,x_{\mathrm{k16}} \\
+ \frac{\mathrm{d}x_{\mathrm{k05k08}}}{\mathrm{d}t} &= a_b\left(a_{\mathrm{k05}\to\mathrm{k05k08}} x_{\mathrm{k05}} + a_{\mathrm{k08}\to\mathrm{k05k08}} x_{\mathrm{k08}}\right) + d\left(x_{\mathrm{k05k08k12}} + x_{\mathrm{k05k08k16}}\right) - a_b\left(a_{\mathrm{k05k08}\to\mathrm{k05k08k12}} + a_{\mathrm{k05k08}\to\mathrm{k05k08k16}}\right)x_{\mathrm{k05k08}} - 2\,d\,x_{\mathrm{k05k08}} \\
+ \frac{\mathrm{d}x_{\mathrm{k05k12}}}{\mathrm{d}t} &= a_b\left(a_{\mathrm{k05}\to\mathrm{k05k12}} x_{\mathrm{k05}} + a_{\mathrm{k12}\to\mathrm{k05k12}} x_{\mathrm{k12}}\right) + d\left(x_{\mathrm{k05k08k12}} + x_{\mathrm{k05k12k16}}\right) - a_b\left(a_{\mathrm{k05k12}\to\mathrm{k05k08k12}} + a_{\mathrm{k05k12}\to\mathrm{k05k12k16}}\right)x_{\mathrm{k05k12}} - 2\,d\,x_{\mathrm{k05k12}} \\
+ \frac{\mathrm{d}x_{\mathrm{k05k16}}}{\mathrm{d}t} &= a_b\left(a_{\mathrm{k05}\to\mathrm{k05k16}} x_{\mathrm{k05}} + a_{\mathrm{k16}\to\mathrm{k05k16}} x_{\mathrm{k16}}\right) + d\left(x_{\mathrm{k05k08k16}} + x_{\mathrm{k05k12k16}}\right) - a_b\left(a_{\mathrm{k05k16}\to\mathrm{k05k08k16}} + a_{\mathrm{k05k16}\to\mathrm{k05k12k16}}\right)x_{\mathrm{k05k16}} - 2\,d\,x_{\mathrm{k05k16}} \\
+ \frac{\mathrm{d}x_{\mathrm{k08k12}}}{\mathrm{d}t} &= a_b\left(a_{\mathrm{k08}\to\mathrm{k08k12}} x_{\mathrm{k08}} + a_{\mathrm{k12}\to\mathrm{k08k12}} x_{\mathrm{k12}}\right) + d\left(x_{\mathrm{k05k08k12}} + x_{\mathrm{k08k12k16}}\right) - a_b\left(a_{\mathrm{k08k12}\to\mathrm{k05k08k12}} + a_{\mathrm{k08k12}\to\mathrm{k08k12k16}}\right)x_{\mathrm{k08k12}} - 2\,d\,x_{\mathrm{k08k12}} \\
+ \frac{\mathrm{d}x_{\mathrm{k08k16}}}{\mathrm{d}t} &= a_b\left(a_{\mathrm{k08}\to\mathrm{k08k16}} x_{\mathrm{k08}} + a_{\mathrm{k16}\to\mathrm{k08k16}} x_{\mathrm{k16}}\right) + d\left(x_{\mathrm{k05k08k16}} + x_{\mathrm{k08k12k16}}\right) - a_b\left(a_{\mathrm{k08k16}\to\mathrm{k05k08k16}} + a_{\mathrm{k08k16}\to\mathrm{k08k12k16}}\right)x_{\mathrm{k08k16}} - 2\,d\,x_{\mathrm{k08k16}} \\
+ \frac{\mathrm{d}x_{\mathrm{k12k16}}}{\mathrm{d}t} &= a_b\left(a_{\mathrm{k12}\to\mathrm{k12k16}} x_{\mathrm{k12}} + a_{\mathrm{k16}\to\mathrm{k12k16}} x_{\mathrm{k16}}\right) + d\left(x_{\mathrm{k05k12k16}} + x_{\mathrm{k08k12k16}}\right) - a_b\left(a_{\mathrm{k12k16}\to\mathrm{k05k12k16}} + a_{\mathrm{k12k16}\to\mathrm{k08k12k16}}\right)x_{\mathrm{k12k16}} - 2\,d\,x_{\mathrm{k12k16}} \\
+ \frac{\mathrm{d}x_{\mathrm{k05k08k12}}}{\mathrm{d}t} &= a_b\left(a_{\mathrm{k05k08}\to\mathrm{k05k08k12}} x_{\mathrm{k05k08}} + a_{\mathrm{k05k12}\to\mathrm{k05k08k12}} x_{\mathrm{k05k12}} + a_{\mathrm{k08k12}\to\mathrm{k05k08k12}} x_{\mathrm{k08k12}}\right) + d\,x_{\mathrm{4ac}} - a_b\,a_{\mathrm{k05k08k12}\to\mathrm{4ac}}x_{\mathrm{k05k08k12}} - 3\,d\,x_{\mathrm{k05k08k12}} \\
+ \frac{\mathrm{d}x_{\mathrm{k05k08k16}}}{\mathrm{d}t} &= a_b\left(a_{\mathrm{k05k08}\to\mathrm{k05k08k16}} x_{\mathrm{k05k08}} + a_{\mathrm{k05k16}\to\mathrm{k05k08k16}} x_{\mathrm{k05k16}} + a_{\mathrm{k08k16}\to\mathrm{k05k08k16}} x_{\mathrm{k08k16}}\right) + d\,x_{\mathrm{4ac}} - a_b\,a_{\mathrm{k05k08k16}\to\mathrm{4ac}}x_{\mathrm{k05k08k16}} - 3\,d\,x_{\mathrm{k05k08k16}} \\
+ \frac{\mathrm{d}x_{\mathrm{k05k12k16}}}{\mathrm{d}t} &= a_b\left(a_{\mathrm{k05k12}\to\mathrm{k05k12k16}} x_{\mathrm{k05k12}} + a_{\mathrm{k05k16}\to\mathrm{k05k12k16}} x_{\mathrm{k05k16}} + a_{\mathrm{k12k16}\to\mathrm{k05k12k16}} x_{\mathrm{k12k16}}\right) + d\,x_{\mathrm{4ac}} - a_b\,a_{\mathrm{k05k12k16}\to\mathrm{4ac}}x_{\mathrm{k05k12k16}} - 3\,d\,x_{\mathrm{k05k12k16}} \\
+ \frac{\mathrm{d}x_{\mathrm{k08k12k16}}}{\mathrm{d}t} &= a_b\left(a_{\mathrm{k08k12}\to\mathrm{k08k12k16}} x_{\mathrm{k08k12}} + a_{\mathrm{k08k16}\to\mathrm{k08k12k16}} x_{\mathrm{k08k16}} + a_{\mathrm{k12k16}\to\mathrm{k08k12k16}} x_{\mathrm{k12k16}}\right) + d\,x_{\mathrm{4ac}} - a_b\,a_{\mathrm{k08k12k16}\to\mathrm{4ac}}x_{\mathrm{k08k12k16}} - 3\,d\,x_{\mathrm{k08k12k16}} \\
+ \frac{\mathrm{d}x_{\mathrm{4ac}}}{\mathrm{d}t} &= a_b\left(a_{\mathrm{k05k08k12}\to\mathrm{4ac}} x_{\mathrm{k05k08k12}} + a_{\mathrm{k05k08k16}\to\mathrm{4ac}} x_{\mathrm{k05k08k16}} + a_{\mathrm{k05k12k16}\to\mathrm{4ac}} x_{\mathrm{k05k12k16}} + a_{\mathrm{k08k12k16}\to\mathrm{4ac}} x_{\mathrm{k08k12k16}}\right) - 4\,d\,x_{\mathrm{4ac}}
+
+The observation model
+^^^^^^^^^^^^^^^^^^^^^^
+
+Each observable is the steady-state abundance of one motif,
+
+.. math::
+
+ y_m(\theta) = x_m^{\mathrm{ss}}(\theta),
+
+for the 15 observed motifs, :math:`m \in \mathcal{O}`. The observed set :math:`\mathcal{O}` excludes the ``x_k05k16`` motif because it was below the quantification limit of the experiment.
+The measured are assumed to follow a log-normal (multiplicative) noise distribution:
+
+.. math::
+
+ \ln \bar{y}_{m,r} = \ln y_m(\theta) + \varepsilon_{m,r},
+ \qquad \varepsilon_{m,r} \sim \mathcal{N}(0, \sigma^2),
+
+where :math:`\bar{y}_{m,r}` is replicate :math:`r` of the measurement of motif
+:math:`m`, and :math:`\sigma =` ``sigma_`` is the noise parameter (fixed to
+``1``).
+
+The likelihood function
+^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The likelihood of the estimated parameters :math:`\theta` (the basal rate
+:math:`a_b` and the estimated motif-specific factors) is
+
+.. math::
+
+ \mathcal{L}(\theta) = \prod_{m \in \mathcal{O}} \prod_{r=1}^{R_m}
+ \frac{1}{\bar{y}_{m,r}\, \sigma \sqrt{2\pi}}
+ \exp\!\left( -\frac{\left(\ln \bar{y}_{m,r} - \ln y_m(\theta)\right)^2}{2\sigma^2} \right),
+
+and the corresponding negative log-likelihood, which PEtab Select uses to
+compute model selection criteria, is
+
+.. math::
+
+ \mathrm{NLLH}(\theta) = \sum_{m \in \mathcal{O}} \sum_{r=1}^{R_m}
+ \left[ \ln\!\left(\bar{y}_{m,r}\, \sigma \sqrt{2\pi}\right)
+ + \frac{\left(\ln \bar{y}_{m,r} - \ln y_m(\theta)\right)^2}{2\sigma^2} \right],
+
+where :math:`R_m` is the number of replicates of motif :math:`m`.
+
+Experimental data
+-----------------
+
+The model is fitted to the relative abundances of the H4 acetylation motifs
+measured in *Drosophila melanogaster* Kc cells. The data come from the
+quantitative mass-spectrometry study of Feller *et al.* (2015), as used by
+Blasi *et al.* [Blasi2016]_:
+
+Briefly, histone H4 was extracted from wild-type Kc cells. The relative abundances of the H4 N-terminal acetylation motifs were quantified by liquid chromatography–mass spectrometry (LC–MS). The measurements represent the steady-state distribution of motifs.
+
+Of the 16 motifs, one (``K5K16``) lies below the quantification limit and is
+not used; hence, the PEtab problem has 15 observables (one per measured motif)
+and 252 measurements in total.
+
+The model selection problem
+---------------------------
+
+There is one hypothesis per acetylation reaction: its rate is either the shared
+basal rate (:math:`a_{p\to q}` fixed to ``1``) or an estimated
+motif-specific rate (:math:`a_{p\to q}` estimated). With 32
+reactions, the model space therefore contains
+:math:`2^{32} \approx 4.3` billion models. The task is to identify the subset of
+reactions that require a motif-specific rate constant to explain the data. The published best model
+had seven motif-specific reactions [Blasi2016]_.
+
+This PEtab Select formulation differs from the original publication in three
+ways in that it omits 11 additional models that were considered in the original publication. These 11 models add negligible computational cost to the model selection problem and were not amongst the best models in the original publication, so are ignored here. Furthermore, we use the FAMoS search method as a general strategy, instead of the highly-tailored problem-specific approach used in the original publication that enabled them to use the brute-force method.
+
+The PEtab Select files
+----------------------
+
+PEtab Select problem YAML
+^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The problem file specifies the criterion, search method, model space file and additional arguments for the search method.
+
+.. literalinclude:: ../../test_cases/0009/petab_select_problem.yaml
+ :language: yaml
+
+The ``candidate_space_arguments`` configure FAMoS:
+
+- ``predecessor_model``: the initial model the search starts from (see below).
+- ``critical_parameter_sets``: empty here — no reaction is forced to always be
+ motif-specific.
+- ``swap_parameter_sets``: a single set containing all 32 reaction parameters.
+ FAMoS *lateral* (swap) moves exchange an estimated parameter for a fixed one;
+ restricting swaps to within a set means any motif-specific reaction may be
+ swapped for any other.
+- ``consecutive_laterals: true``: keep performing lateral moves while they keep
+ improving the model.
+- ``summary_tsv``: where to write a summary of the search status and history.
+
+Model space
+^^^^^^^^^^^
+
+The model space is a single subspace ``M`` with 32 parameter columns, one
+per acetylation reaction. Every column has the value ``1.0;estimate``, meaning
+each reaction can be either fixed to the basal rate (``1.0``) or estimated
+(motif-specific). This concisely encodes all ~4.3 billion
+models in a single row (the table is wide, with one column per reaction):
+
+.. csv-table::
+ :file: ../../test_cases/0009/model_space.tsv
+ :delim: tab
+ :header-rows: 1
+
+The PEtab problem (``petab/``)
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The model space references a standard PEtab problem that defines the superset
+model (all reactions present):
+
+- ``model.xml``: the SBML model with the 16 motif species and the acetylation /
+ deacetylation reactions.
+- ``parameters.tsv``: the 32 reaction-rate factors ``a_`` (``log10``
+ scale, bounds ``[1e-3, 1e3]``), the basal acetylation rate ``a_b``, the
+ deacetylation reference rate ``da_b`` (fixed to ``1``), and the noise
+ parameter ``sigma_``.
+- ``observables.tsv``: 15 observables, one per measured motif, and the log-normal noise distribution.
+- ``measurements.tsv``: the 252 steady-state (``time = inf``) abundance
+ measurements.
+- ``conditions.tsv``: a single dummy condition.
+
+Predecessor (initial) model
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+``predecessor_model.yaml`` is the model the FAMoS search is initialised from, in the PEtab Select model YAML file format. It
+is a specific starting model, with 15 motif-specific reactions, that was found
+to reproducibly lead to the best model.
+Without it, solving the problem from scratch requires on the order of 100
+randomly initialised FAMoS searches.
+
+.. literalinclude:: ../../test_cases/0009/predecessor_model.yaml
+ :language: yaml
+
+Expected results
+^^^^^^^^^^^^^^^^
+
+``expected.yaml`` is the expected selected model, in the PEtab Select model YAML file format. It has the seven
+motif-specific reactions of the published best model
+(``a_0ac_k08``, ``a_k05_k05k12``, ``a_k12_k05k12``, ``a_k16_k12k16``,
+``a_k05k12_k05k08k12``, ``a_k12k16_k08k12k16``, ``a_k08k12k16_4ac``) plus the
+basal rate ``a_b`` — eight estimated parameters in total — with
+``AICc ≈ -1708.1``.
+
+.. literalinclude:: ../../test_cases/0009/expected.yaml
+ :language: yaml
+
+``expected_summary.tsv`` is the FAMoS search trajectory, one row per iteration,
+showing how the method switches between ``forward``, ``backward``, and
+``lateral`` moves and how the criterion improves at each step:
+
+.. csv-table::
+ :file: ../../test_cases/0009/expected_summary.tsv
+ :delim: tab
+ :header-rows: 1
+
+Why this problem is challenging
+-------------------------------
+
+This problem is difficult to solve. The search space is large (~4.3 billion models); hence, many models have criterion values that differ by less than numerical noise, so the
+ranking of near-optimal models is effectively non-deterministic across
+machines and tolerances. Plain forward or backward selection mostly fails to
+reach the optimum.
+
+The FAMoS method reproducibly
+converges to models with markedly better criterion values.
+Multi-start FAMoS searches recover the published best model while calibrating
+only a small fraction (~0.002 %) of the full model space, making this
+large-scale problem computationally feasible.
+
+Because of the numerical-noise sensitivity noted above, when running this test
+case you should expect to obtain a similar (but not necessarily identical)
+``expected_summary.tsv`` (a few rows may be reordered, or the path through model space
+may differ). However, the improvement in criterion value over consecutive rows should be conserved, and the same select model should be found.
+
+References
+----------
+
+.. [Blasi2016] Blasi T, Feller C, Feigelman J, Hasenauer J, Imhof A, Theis FJ,
+ Becker PB, Marr C. *Combinatorial Histone Acetylation Patterns Are Generated
+ by Motif-Specific Reactions.* Cell Systems, 2016, 2(1):49–58.
+ https://doi.org/10.1016/j.cels.2016.01.002
diff --git a/doc/index.rst b/doc/index.rst
index b9c73e63..40de99d3 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -41,7 +41,8 @@ Python package or command line interface. An example of this is provided for the
For users, in most
cases, model selection is performed entirely via one of the calibration tools,
-rather than the PEtab Select package directly. After model selection, the
+rather than the PEtab Select package directly. Hence, we recommend reading the `documentation from your calibration tool on model selection `_.
+After model selection, the
analysis and visualization methods in the PEtab Select package can be used
directly with the results from your calibration tool.
diff --git a/doc/test_suite.rst b/doc/test_suite.rst
index 4684963a..015e903f 100644
--- a/doc/test_suite.rst
+++ b/doc/test_suite.rst
@@ -75,4 +75,4 @@ the model format.
.. [#f2] Noise parameter is removed, noise is fixed to ``1``.
-.. [#f3] This is a computationally expensive problem to solve. Developers can try a model selection initialized with the provided predecessor model, which is a model start that reproducibly finds the expected model. To solve the problem reproducibly *ab initio*, on the order of 100 random model starts are required. This test case reproduces the model selection problem presented in https://doi.org/10.1016/j.cels.2016.01.002.
+.. [#f3] This is a computationally expensive problem to solve. Developers can try a model selection initialized with the provided predecessor model, which is a model start that reproducibly finds the expected model. To solve the problem reproducibly *ab initio*, on the order of 100 random model starts are required. This test case reproduces the model selection problem presented in https://doi.org/10.1016/j.cels.2016.01.002. A description of this model selection problem is also provided at :ref:`test_case_0009`.