Algebra, Geometry and Topology of ERK Kinetics

The MEK/ERK signalling pathway is involved in cell division, cell specialisation, survival and cell death (Shaul and Seger in Biochim Biophys Acta (BBA)-Mol Cell Res 1773(8):1213–1226, 2007). Here we study a polynomial dynamical system describing the dynamics of MEK/ERK proposed by Yeung et al. (Curr Biol 2019, 10.1016/j.cub.2019.12.052) with their experimental setup, data and known biological information. The experimental dataset is a time-course of ERK measurements in different phosphorylation states following activation of either wild-type MEK or MEK mutations associated with cancer or developmental defects. We demonstrate how methods from computational algebraic geometry, differential algebra, Bayesian statistics and computational algebraic topology can inform the model reduction, identification and parameter inference of MEK variants, respectively. Throughout, we show how this algebraic viewpoint offers a rigorous and systematic analysis of such models.

Michalis-Menten constant associated to compound C i QSSA Quasi-steady-state-approximation f (x, θ) Rates of change of concentrations given x and θ V( f ) Vanishing ideal of a polynomial f Y θ , V θ Parameter/QSSA-variety of θ κ i , π, γ i Parameters of reduced models φ t 1 ,...,t l (θ ) The model prediction map at times t 1 , . . . , t l and parameters θ K A differentially closed field I The differential ideal studied to determine structural identifiability C(θ ){x, y} The differential ring in indeterminates x, y over the fraction field C(θ ) k The The bandwidth of a kernel KDE Kernel density estimation VR b (v) A Vietoris-Rips complex on vertices v at resolution b B A barcodê π A p-value estimate

Introduction
In systems biology, dynamics play a crucial role in cellular decision making (e.g. whether a cell responds appropriately to a particular signal) (Voit 2017). Molecular interactions can be modelled as systems of chemical reactions with a choice of kinetics, such as the law of mass action, which assumes that the rate at which a chemical reaction proceeds is proportional to the product of the concentrations of its reactants. From a finite set of reactions, the mass-action modelling assumption gives rise to a system of polynomial ordinary differential equations (ODEs), which are sums of monomials in which each term includes concentrations of molecular species as variables and coefficients as rates of reaction. Chemical reaction network theory (CRNT) is a mathematical field developed by Horn and Jackson, and independently by Bykov, Gorban, Volpert and Yablonsky, for analysing such reactions, and the mathematical techniques employed extend beyond dynamical systems theory to include algebraic geometry, differential algebra, algebraic statistics and discrete mathematics (Dickenstein 2016). CRNT often focuses on steady-state analysis through the lens of computational and real algebraic geometry, asking questions about the capability or preclusion of multiple positive real steady-states (i.e. multistationarity) or more complex dynamics, often without requiring specialised parameter values (Banaji and Craciun 2009;Craciun and Feinberg 2005;Millán et al. 2012;Angeli 2009;Wang and Sontag 2008;Feliu and Wiuf 2012;Müller 2016;Conradi and Pantea 2019). Multi-site protein phosphorylation systems, such as the ERK/MEK signalling pathway, can be translated into such chemical reactions and their multistationarity, corresponding to different biological cellular decisions, has attracted much attention (Thomson and Gunawardena 2009b;Gunawardena 2007;Aoki 2011;Takahashi et al. 2010;Markevich et al. 2004). Algebraic analyses and invariants of multi-site phosphorylation have revealed geometric information of steady-state varieties, informed experimental design and enabled model comparison using steady-state data (Manrai and Gunawardena 2008; Thomson and Gunawardena 2009a;Harrington et al. 2012;Gross et al. 2016;MacLean et al. 2015). However, such systems have also been shown to exhibit nontrivial transient dynamics and oscillations Qiao et al. 2007). In recent years, the fields of systems biology and CRNT have extended the repertoire of techniques to assert other dynamics (Banaji 2020;Domijan and Kirkilionis 2009;Mincheva and Roussel 2007;Kay 2017;Errami 2015;Angeli et al. 2013), reduce models systematically (Pantea et al. 2014;Goeke et al. 2017;Feliu et al. 2019;Sweeney 2017;Boulier et al. 2011;Hubert and Labahn 2013), and assess identifiability (Ljung and Glad 1994;Ollivier 1990;Meshkat et al. 2009;Hong et al. 2020;Bellu et al. 2007). Furthermore, combinatorial structures, such as simplicial complexes, and techniques from computational algebraic topology have enabled comparison of chemical reaction network models and their parameters (Vittadello and Stumpf 2020;Nardini et al. 2020).
A previous algebraic systems biology case study (Gross et al. 2016) analysed a chemical reaction network model at steady state, by studying the steady-state ideal, chamber complex and algebraic matroids of the model. Here we present a sequel of such analysis to study the dynamics of chemical reaction networks with time-course data, which relies on studying the QSS variety (Sect. 3), the model prediction map (Sect. 4) and the topology of a parameter inference (Sect. 5).
We perform a detailed mathematical analysis of recently published models and experimental data (Yeung et al. 2019). The Full ERK model describes dual phosphorylation of ERK by MEK, two molecular species whose activation regulates cell division, cell specialisation, survival and cell death (Shaul and Seger 2007). The dynamics of the six ERK/MEK molecular species x ∈ R n=6 in the Full ERK model are governed by a polynomial dynamical systemẋ(t) = f (x(t), θ ), where θ ∈ R m=6 is the vector of parameters and there are two conservation relations between the species. The Full ERK model is presented in Sect. 2. Analysing the kinetic parameters of a model depends on the available data. The accompanying time-course experimental observations include measurements of ERK in 3 different states, at 7 time points following activation by its activated enzyme kinase MEK, which is either wild-type (WT) or mutated MEK. Mutations of MEK are known to be involved in human cancer and embryonic developmental defects; therefore, understanding their kinetics and differences between wild-type and 4 mutants (e.g. Y130C, F53S, E203K or SSDD) may increase fundamental biological understanding of the pathway and contribute to the development of potential therapies. The experimental data and relevant biological information are presented in Sect. 2.
Using algebraic approaches first presented by Goeke, Walcher and Zerz in Goeke et al. (2017), we decrease the number of variables and parameters in the Full ERK model. We derive two model reductions: the Rational ERK model and the Linear ERK model. We show, with known biological information (see Sect. 2), that the reduction to the Linear ERK model by Yeung et al. (2019) is mathematically sound. We note that the Rational ERK model was not analysed in Yeung et al. (2019), although it can be derived from the Full ERK model using singular perturbation methods. A natural question is whether a quasi-steady-state approximation is justified given the experimental setup, which equates to solving an algebraic problem (Goeke et al. 2017). We identify algebraic varieties V θ that are (analytic) invariant sets of the ODE system and characterise neighbourhoods in parameter space for which the ODE solutions remain close to these varieties. This systematic analysis allows us to simplify the model equations such that the dynamics of both reduced models are good approximations to the Full ERK model. Algebraic model reduction and derivation of the reduced ERK models are given in Sect. 3.
Before estimating the parameters of a model from observations, one must determine its identifiability. Identifiability is concerned with asking whether it is possible to recover values of the model parameters given data. A model is structurally identifiable if parameter recovery is possible with perfect data. Mathematically, this task is equivalent to asking whether the model prediction map is injective. The model prediction map, defined precisely in Sect. 4.1, is a map that takes a parameter to the corresponding predicted noise-free data point(s) (Dufresne et al. 2018). Real data is often noisy; testing whether parameter recovery is possible with imperfect data is the problem of practical identifiability (Raue et al. 2009;Dufresne et al. 2018). Mathematically, measurement noise induces a probability distribution in data space. Assuming that the model prediction map is injective (at least generically), practical identifiability can be defined in terms of the boundedness (with respect to a reference metric in parame-ter space) of the confidence regions of a likelihood test. Under our assumptions, this translates to asking whether the preimages of small bounded regions in data space are bounded in parameter space. We prove the following:

Theorem 1 The Linear ERK/MEK model, with the given experimental setup (number of species, number of replicates, number of measurement time-points and initial conditions), is structurally and practically identifiable.
We provide a definition of practical identifiability that improves a previous definition (Dufresne et al. 2018), and which is an alternative to that of Raue et al. (2009). We also propose a computable algorithm for practical identifiability, implement it and apply it to the ERK models. We prove Theorem 1 in Sect. 4. We use the differential algebra method to show that the Full ERK model and the Rational ERK model are generically structurally identifiable. These results are guaranteed to be valid if we have at least 2m + 1 generic time points by Sontag's result (Sontag 2002), but can be valid with fewer generic time points. Indeed, as the Linear ERK model admits analytic solutions, we can prove that it is globally structurally identifiable for any choice of three distinct time points. Determining structural identifiability for specific time points in the absence of analytic solutions is an open problem.
We numerically show that the Full ERK model and Rational ERK model are not practically identifiable; however, the source of this practical non-identifiability is not completely clear (see Sect. 4).
Finally, for a model that is structurally and practically identifiable, one would like to infer parameters, i.e. what parameter values are consistent with the observations? We perform Bayesian inference, as done in Yeung et al. (2019), and extend this to the Rational ERK model. The result of the parameter inference on the Linear ERK model is a sample point cloud of posterior densities of inferred ERK parameter kinetics that are consistent with the data. We obtain five different sample densities corresponding to the five MEK variants.
In Sect. 5, we compare the geometry of the admissible regions of parameter space of the five MEK variants. The computational field of topological data analysis quantifies the shape and connectivity of data through computation of topological invariants across resolutions (or threshold values) of metric data. In recent years, topological methods have dramatically improved in computational speed as well as theoretical advancements that facilitate the analysis of scientific datasets. We implement a theoretical framework originally proposed by Taylor (2019a), in order to quantify the shape of the resulting posterior distributions of kinetic parameters and facilitate a comparison between mutants. Specifically, their theorem provides a direction for hypothesis testing of two densities using distances between topological summaries. The framework relies on approximating the persistent homology of super-level sets of posterior densities by simplicial complexes. We perform these measurements on the distributions obtained from Bayesian parameter inference for the 5 MEK variants and compare them via a topological bottleneck distance.

Biological Result
The topological data analysis quantifies that the Linear ERK model parameter posteriors are most different between the WT and SSDD mutant data. The kinetics of the SSDD mutant, which mimics phosphorylated MEK, has the largest topological distance from all other MEK/ERK mutants. This biological result raises the question of whether the SSDD variant is a suitable approximation for wild-type MEK activated by Raf, and suggests further experimental studies are needed. While the previous analysis by Yeung et al. (2019) compared the variants by the inferred kinetics of each parameter, here we complement that analysis by comparing the three parameters together as a point cloud.
Our aim is to showcase how systematic algebraic, geometric and topological approaches can be applied to a biologically relevant model with state-of-the-art experimental time-course data. Each of these approaches incorporates the structure of the mathematical model, experimental observations, and experimental setup and observations (e.g. experimental initial condition, observable species, number of experimental replicates, number of time points collected, etc.), as well as known biological information (e.g. published parameter values). Due to the multiple disciplines and different notation conventions (as well as standard abbreviations), we include a glossary of symbols at the start of the paper. The framework we present is not limited to this case study and may enhance the analysis of similar models in systems and synthetic biology.

From ERK Biochemical Reactions to a Polynomial Dynamical System
Protein phosphorylation alters protein function in signalling pathways and plays a crucial role in cellular decisions and homeostasis. Phosphorylation is the addition of a phosphate group by an enzyme known as a kinase, and dephosphorylation is the removal of a phosphate group by an enzyme known as a phosphatase. Multisite phosphorylation is the process of having multiple possible locations on a protein phosphorylated, which increases the number of potential ways protein function can be altered. The algebra, geometry, combinatorics and dynamics of multisite phosphorylation has been a source of interesting mathematical problems (Dickenstein 2016;Manrai and Gunawardena 2008;Conradi and Pantea 2019). A protein with q phosphorylation sites has been shown to have 2 q phospho-states; the sites on the protein can be phosphorylated in q! possible ways (Thomson and Gunawardena 2009b). One of the simplest multisite phosphorylation systems is when a protein has two phosphorylation sites. We focus on the sequential dual phosphorylation of the extracellular signal regulated kinase (ERK) by its kinase activated (dually phosphorylated) MEK. The model developed by Yeung et al. (2019) encodes a mixed phosphorylation mechanism (i.e. distributive and processive) by changes in parameter values rather than separate models (see, for example, Conradi and Shiu 2015;Gunawardena 2007 and references therein). This enabled them to quantify the extent to which a MEK variant is processive or distributive. We remark that the model presented by Yeung et al. does not include dephosphorylation mechanisms, since the experimental setup omitted the addition of phosphatases.
Next, we introduce the model and the experimental data published by Yeung et al. (2019).

The Model
The protein substrate ERK, is activated through dual phosphorylation by its activated enzyme kinase MEK. As shown in the chemical reaction network (see Fig. 1), unphosphorylated ERK (S 0 ) binds reversibly with its kinase MEK (E) to form an intermediate complex C 1 . The complex becomes C 2 when a phosphate group is added. Complex C 2 can then disassociate to form MEK (E) and ERK phosphorylated on the tyrosine site (S 1 ), or a second phosphate group is added to C 2 , resulting in product reactants E + S 2 . The six species and six rate constants are given in the following chemical reaction network (Fig. 1).
We can translate this reaction network into a dynamical systemẋ = f (x, θ). Here, f is a vector-valued function of the vectors of species concentrations x = {S 0 , C 1 , C 2 , S 1 , S 2 , E} and rate constants, referred to as parameters θ = {k f 1 , k r 1 , k c 1 , k f 2 , k r 2 , k c 2 }. The kinetics assumption for f is a modelling choice; here we assume that the law of mass action holds (Klipp et al. 2016, §2.1.1), as for the original model (Yeung et al. 2019). The resulting dynamical system of ODEs is given in Eqs. (1).
We assume that initially all species are zero, except for S 0 (t = 0) = S tot and E(t = 0) = E tot . Equations (2)-(3) define two conserved quantities that constitute a basis for the linear space of conservation relations of the model: where the total amounts of substrate ERK (S tot ) and enzyme MEK (E tot ) are constant and known from the initial conditions. We aim to study the relationship between the species x, parameters θ , conserved quantities and available biological information (previous knowledge and experimental observations). The emphasis in this paper is not to analyse the steady-state variety as in Gross et al. (2016), rather here we focus on the transient dynamics of the model and algebraic approaches to analyse ERK kinetics in light of the available biological information.

The Data
The data is published. Details on measurement techniques and experimental methods can be found in Yeung et al. (2019). We present the experimental setup for the data we analyse.

Known Biological Information
The relationship between some kinetic rate constants is known. When a substrate binds reversibly to an enzyme to form an enzyme-substrate complex, which then reacts irreversibly to form a product and release the enzyme, one can define the Michaelis-Menten constant k M . In the reaction network given by Eqs. (1), there are two Michaelis-Menten constants k M i = (k c i + k r i )/k f i for i = 1, 2. Measurements show that in our experimental setup k M i ≈ 25μM for i = 1, 2 (Taylor 2019b). While the reaction rates k c i and k r i for i = 1, 2 cannot be measured directly, they have been shown to be of the same order of magnitude (Bar-Even 2011). We will use these insights to assume, henceforth, that S 0 , S 1 and S 2 were measured (without added compound variables). We justify this mathematically in Sect. 3.2.4.

Algebraic Model Reduction
The first step to studying most models typically involves model reduction, which reduces the number of dependent variables and constant parameters. For many chemical reactions, there are time scales on which the rate of change of some variables is negligible and their dynamics is dominated by those of the remaining variables. This observation motivates the Quasi-Steady-State-Approximation (QSSA).
In recent years, algebraic approaches to reduce polynomial ODE models have been extended by Walcher and colleagues. In 2013, Pantea et al. (2014) used Galois theory to characterise chemical reaction networks for which no explicit QSSA reduction is possible. Furthermore, they provided computational tools for determining the feasibility of an explicit reduction. Subsequently, Sweeney (2017) proved that the nonsolvability of polynomials poses no issue to the CRNs most commonly encountered in practice and derived a more efficient algorithm for determining explicit reducibility by translating algebraic structures into graphs. Goeke and Walcher (2014) provide an explicit formula for obtaining a reduced QSSA model using a subset of an algebraic variety defined by the slow manifold. Subsequently, Goeke et al. (2017) characterised parameter values at which QSSA reduction is accurate using algebraic varieties and bounds on the polynomials governing the ODE system on a bounded parameter-and variabledomain. Most recently, Feliu et al. (2019) derived necessary and sufficient conditions for purely algebraic reductions of a CRN model to agree with model reductions derived via classical singular-perturbation theory (Keener and Sneyd 2011;Segel 1988).
In this section, we briefly review QSSA using classical singular-perturbation theory as well as the algebraic approaches developed by Goeke et al. (2017); Goeke and Walcher (2014). We then apply both methods to the full ERK model (Eqs. (1)). We show both approaches can generate the same QSSA-reduction of our ERK model, which we will call the Rational ERK model. Additionally, the algebraic method can yield a linear QSSA-reduction of our ERK model in a single step (which we call the Linear ERK model). By contrast, the singular-perturbation-theory approach requires additional assumptions on parameter values to arrive at the Linear ERK model (see Sect. 3.2). We show that the Linear ERK model approximates the Full ERK model (Eqs. (1)) with similar accuracy as the Rational ERK model in the context of the experimental setup, data and known biological information (see Sect. 3.2; Appendix A.2 for details).
With the algebraic method, we provide a rigorous mathematical justification of the Linear ERK model presented by Yeung et al. (2019). By comparing the singular perturbation method with the algebraic method and the two resulting model reductions, we illustrate how the algebraic methods form a well-structured approach for arriving at a QSS reduction and for assessing the accuracy of such reductions systematically.

Notation for Model Reduction
Throughout, we will assume we have an ODE system in variables x ∈ R n and parameters θ ∈ R m . If the system dynamics are governed by f , a vector of polynomials in R[x, θ] n , then our ODE system is given by For 1 ≤ q < n, we may define We wish to retain the variables x [1] in the reduced model and seek to eliminate variables x [2] as part of our model reduction.

Remark
In the current section, we treat the (non-zero) initial conditions of the ODE systems as parameters (and include them in the parameter count m), as they are central to determining the goodness of a model reduction. In Sect. 4 (Identifiability) and Sect. 5 (Inference & Comparison), we will not include the initial conditions in the set of parameters, as they are given by the experimental setup and, as such, do not need to be identified or inferred.

The Algebraic QSSA Approach
The algebraic approach to QSSA, as presented by Goeke, Walcher and Zerz in Goeke et al. (2017), differs from the classical approach in several ways. Most notably, an a priori separation of time scales is not needed. On the other hand, we require a choice of fast and slow variables (i.e., a choice of which variables we eliminate, and which we retain in the reduced model).

Remark
To the best of our knowledge, all existing algebraic approaches to QSSA, including (Goeke et al. 2011;Goeke and Walcher 2014;Boulier et al. 2011;Goeke et al. 2017), require a choice, explicit or implicit, of slow and fast variables. In Goeke et al. (2011) the relevant choice is made by expanding f (in Goeke et al. (2011): h) at different orders of ε.
First, we characterise points in parameter space, i.e. parameter values, where the fast variables are exactly determined by the slow variables, which yields a reduced model. This set of parameter values is defined as the vanishing set of the polynomials governing the ODEs of the fast variables. This defines an algebraic variety in the parameter space. Typically, the ODE system will be degenerate at these values. Secondly, we characterise neighbourhoods of these values in parameter space, as well as time scales, for which the reduction is a good approximation to the original model.
To describe the characterisation from Goeke et al. (2017), we use x [1] , x [2] , f [1] and f [2] as before. In addition, we denote the partial derivative with respect to x [i] by D i . For a fixed θ * ∈ R m , we let Y θ * denote the algebraic variety defined by f [2] ( · , θ * ).
Definition 2 Let y ∈ Y θ * be such that the (n −q)×(n −q) matrix D 2 f [2] has full rank at (y, θ * ). Then we denote by V θ * ⊆ Y θ * a relatively Zariski-open neighbourhood of y in which this rank is maximal. We call V θ * a quasi-steady-state (QSS) variety in the sense of Goeke et al. (2017) and may assume without loss of generality that it is irreducible.
If, furthermore, V θ * is an invariant set of the ODE system dx [1] /dt = f (x, θ * ), then we call θ * a QSS parameter value. Recall that in dynamical systems theory, V θ * is an invariant set of R q if whenever the initial condition of an ODE at t = 0 is in V θ * , then the corresponding trajectories of the ODE remain in V θ * for all t > 0.
Remark Note that the steady-state variety (see Gross et al. 2016) and the QSS variety at a parameter value θ * are not as closely related as one may first think. Indeed with our notation, the steady state variety is the zero set θ)), but the steady-state variety and V θ * are not contained in one another in general.
To apply the theory of Goeke, Walcher and Zerz in Goeke et al. (2017), we assume that the initial condition of our ODE system (Eq. (4)) lies in V θ * . As D 2 f [2] has full rank on V θ * , we have that x [2] = x [1] for some continuous by the Implicit Function Theorem. Hence, writing x = (x [1] , x [2] ), we obtain a reduced model: on some open neighbourhood in R j that naturally includes V θ * . This corresponds to determining the fast variables in terms of the slow variables. We do this by setting their time rates of change equal to zero on the short timescale in classical QSSA, with the addition that on V θ * the above yields an exact solution rather than an approximation. As a caveat, we note that, in both settings, it may not be possible to find an algebraic expression for ; this was pointed out and completely characterised by Pantea et al. (2014) in terms of Galois theory. Because of the possible non-solvability issue with Eq. (5), we require a more general methodology (Proposition 3) to study the accuracy of a model reduction (Proposition 4). Goeke, Walcher and Zerz showed that locally, in the variable x [1] , the reduced system given by Eq. (5) has the same solution as the following ODE system

Proposition 3 (Lemma 1 & Proposition 1 in Goeke et al. (2017)) Let V θ * be a QSSvariety. Then V θ * is an invariant set of Eq. (5). Moreover, any solution of Eq. (6) with initial condition in V θ * locally solves Eq. (5). Conversely, any solution of Eq. (5) with initial condition in V θ * locally solves Eq. (6). In addition, V θ * is an invariant set of Eq. (4) if and only if the solutions of Eqs. (4) and (6) are equal for all initial conditions
This proposition equips us with a method to obtain a solution for x [1] in an algebraic QSSA without explicitly determining . In Sects. 4 and 5, we will use Eq. (5) as our model reduction.
First, however, we assess the accuracy of Eq. (6) as an approximation to the full system, for parameter-values θ in some neighbourhood of θ * . For convenience, we abbreviate system (6) as dx/dt = f red (x, θ * ).
Proposition 4 (Outline of Proposition 2 in Goeke et al. (2017)) Let K * ⊆ R n + × R m + be a compact domain in the product of the variable and parameter spaces which satisfies a number of conditions (we refer the interested reader to Appendix A.2 for details). Let θ * be given such that V θ * × {θ * } has non-empty intersection with int K * , let (y, θ * ) be a point in this intersection, and let V θ * be an open neighbourhood of y such that Then there exists a compact neighbourhood A θ * ⊆ V θ * of y such that:

with initial condition z, exists and remains in
For every ε > 0, there exists a δ ∈ (0, δ 1 ] such that, for any z ∈ V θ * ∩ A θ * , the difference between the solutions of Eqs. (6) and (4), with initial condition z, is at most ε for t ∈ [0, t * ] whenever f − f red < δ on V θ * . Here, · denotes the infinity-norm over the interval [0, t * ] for a fixed parameter value.
In summary, given some technical assumptions on the variables and the domain K * , we can bound the difference between the solutions of Eqs. (4) and (6) in terms of f − f red up to some time t * > 0. The full statement of this proposition also includes lower bounds on this difference. Note that we do not assume that θ * is a QSS-parameter value, but the assumptions on K * (as detailed in Appendix A.2) require it to be close to some QSS-parameter value.

Reducing the ERK Model Algebraically
We now apply the theory from Sect. 3.1 to the Full ERK model (Eqs. (1)) in two different ways to derive two reduced models (the Linear ERK and Rational ERK models). The full details of the derivations can be found in Appendix A.1. We also give a brief biological explanation of why both systems explain the phenomena underlying the given experimental data equally well.

Reduction via Conservation Laws
We can exploit the conservation laws (2) and (3) to eliminate a variable before using the analytic or algebraic QSSA approach. First, we choose to eliminate E and note that there are two choices: Subsequently, we choose to eliminate the variables C 1 and C 2 via (algebraic) QSSA. For the Rational ERK model, using (7) to eliminate E, we obtain while for the Linear ERK model, employing substitution (8), we have

Reduction via an Algebraic QSSA
To reduce the model further, we apply an algebraic QSSA, as described in Sect. 3.1. We start by identifying QSS-parameter-values. For f [2] rat , we have In both cases, assuming that (k r i + k c i ) > 0 for i = 1, 2 (otherwise, the reaction network would be degenerate, meaning some or all variables would remain constant), and given that S 0 and S 1 are non-constant, we deduce that these matrices are invertible. Hence, both substitutions (7) and (8) are good candidates for an algebraic QSSA reduction. We note that the assumption E tot = 0 is required to ensure that the initial condition lies in V θ * . This is not physically realistic, as the absence of free enzyme makes the reaction rates negligible, however, in parameter space this assumption is close to the experimental setup (E tot ≈ 0.65μM). In fact, unlike the rate parameters, we know the value of E tot and can, therefore, bound the error associated with such an idealisation (cf. Appendix A.2). The assumption that E tot = 0 is similar to the classical singular-perturbation theory approach, where a typical choice of short timescale is (t S = E tot k f 1 ) and one then subsequently assumes ε = E tot /S tot → 0.
As E tot = 0 will yield a stationary model and ensure that V θ * contains the initial condition, we find that any parameter value θ * satisfying (k r i + k c i ) > 0 for i = 1, 2 and E tot = 0 is a QSS-parameter-value for both the Rational and Linear ERK model.
For both models, we have For the Linear ERK model, we can show that Y lin θ * is irreducible (at generic parameter values) and thus its QSS-variety is V lin . At generic parameter values, only the first irreducible component will contain the initial condition. Hence, the natural choice for the QSS-variety is The substitution (7) yields the Rational ERK model given by while the substitution (8) gives the Linear ERK model: Here, for i = 1, 2, we use the newly introduced quantities Both models are reductions obtained via the ODE system (5). The processivity parameter, which is the probability that both phosphorylations are carried out by the same enzyme, is represented by π in the reduced models. The κ i represents the kinetic efficiencies of the first and second phosphorylation steps, respectively (Yeung et al. 2019).
It should be noted that the Rational ERK model is the system we would obtain via the classical singular perturbation approach (Keener and Sneyd 2011).

Assessing Accuracy
We can use the algebraic framework of Goeke, Walcher and Zerz and, in particular, Proposition 4 to bound the error of the Linear ERK model reduction to the full model. Given the measurements of the Michaelis-Menten constants k M i , we can derive simple expressions which bound the approximation error (see Appendix A.2 for both the Rational & Linear ERK model). Unfortunately, the bound on the approximation error depends on parameters with unknown values. However, we can compare the bounds derived for the Linear ERK model to those for the Rational ERK model and show that in the regime where k M i ≈ 25μM, both approximate the full model equally well (see Appendix A.2).
Recall that we can also derive the Rational ERK model via singular perturbation theory. When using perturbation theory, it is uncommon to bound the approximation error as explicitly as we do via the algebraic methods of Goeke et al. (2017). However, we can still show that the Linear ERK model is a good approximation of the Rational ERK model when 0 ≤ γ 1 , γ 2 1. Again, we can use knowledge of the Michaelis-Menten constant to show that in our experimental setup, γ 1 and γ 2 are small. Indeed, we can rewrite Since k M i ≈ 25μM and the parameters k c i and k r i are of similar magnitude (see Bar-Even 2011), we conclude that γ 1 ≈ 1/25 (1/μM). We reiterate that by employing an algebraic approach, we can derive a reduced model (without taking further limits) that approximates the Full ERK model as well as that obtained via singular perturbation theory, but has several advantages: it has fewer parameters, is interpretable as a chemical reaction network, and identifiable, as discussed in the next section.

Choice of Output Variables
Recall from Sect. 2.2, the experimental measurements correspond to the following linear combinations of variables: S 0 + C 1 , S 1 + C 2 and S 2 . Here we argue that in the context of available data, S 0 , S 1 and S 2 are sufficient approximations of the output variables, which simplifies both the identifiability analysis and the parameter inference.
We argued in Sect. 3.2.3 that in the context of experimental data the Linear ERK model is as good of an approximation to the Full ERK model as the Rational ERK model. On the long timescale, substitutions for C 1 and C 2 from the Linear ERK model give approximately Recall that k M i ≈ 25μM and E tot = 0.65μM. We then find that the measurements of S i + C i+1 will be dominated by S i . Henceforth we will use S i interchangeably with our measurements S i + C i+1 .

Identifiability
One of the goals of this ERK study is to determine the kinetic parameters of the models given the data. Each model and experimental setup induces a map from the space of model parameters to observable model solutions (here, this is the measurement of the 3 species at the 7 time points over the course of r experimental replicates, i.e. a subset of R 21r ). We call this map φ t 1 ,...,t 7 : → R 21r the model prediction map (see Dufresne et al. 2018). Here, the parameter space is a subset of the positive octant R 6 ≥0 for the Full ERK model, R 5 ≥0 for the Rational ERK model, and R 3

≥0
for the Linear ERK model. One can think of the data as being a point z * in the space of observable model solutions, i.e. R 21r , and parameter estimation corresponds to attempting to compute the inverse image φ −1 t 1 ,...,t 7 (z * ) of this map at that point. Structural identifiability generally corresponds to the model prediction map φ t 1 ,...,t 7 being injective. Real-world observations are noisy, hence the data point z * may not be in the image of the map φ t 1 ,...,t 7 . Thus, when performing parameter estimation, we instead search for parameters yielding model predictions close to the data point z * . Practical identifiability broadly corresponds to having the set of parameters with model predictions close to the data point z * being bounded. In Sect. 4.1 we show that the Linear ERK model is structurally identifiable on its whole parameter space, while the Rational ERK model and the Full ERK model are structurally identifiable on some open dense subset of their parameter space. In Sect. 4.2 we show that the Linear ERK model is practically identifiable for our experimental data, providing the proof of Theorem 1. By contrast, we provide evidence that the Rational ERK model and Full ERK model are not practically identifiable.

Structural Identifiability
First, we study the structural identifiability of our ODE models, that is whether the model prediction map φ t 1 ,...,t 7 : → R 21r is one-to-one, or at least locally one-to-one. We start by providing a formal definition of structural identifiability for models given by ODE systems with specific time points. Suppose we have a rational ODE system in variables x ∈ R n and parameters θ ∈ R m , given by where f is a vector of rational functions in R(x, θ) n . We assume that the measurable output is y = g(x, θ) where g is also a vector of rational functions. Letx(θ, t) be a solution of (12) for the parameter value θ ∈ and then letŷ(θ, t) = g(x(θ, t), θ ) be the observable solution for the same parameter value. Then, supposing that there are r replicates of the experiment, for the specific time points t 1 , . . . , t l the model prediction map is given by The model prediction map then induces an equivalence relation ∼ t 1 ,...,t l on the parameter space via for any θ, θ ∈ .
Definition 5 (c.f. Definition 2.8 in Dufresne et al. (2018)) Suppose we have a model given by a system of rational ODEs (as above) with parameter space and model prediction map φ t 1 ,...,t l . We say a model is: • globally identifiable if every equivalence class of ∼ t 1 ,...,t l on has size exactly 1.
• generically identifiable if for almost all θ ∈ the equivalence class of θ has size exactly 1. • locally identifiable if for almost all θ ∈ the equivalence class of θ is finite.
• generically non-identifiable if for almost all θ ∈ the equivalence class of θ is infinite.
Here "almost all" means everywhere except possibly in a closed subvariety (i.e. the set of common zeroes of some polynomials).
There are several approaches to assess structural identifiability. All identifiability methods involve a certain number of assumptions of genericity, but not always explicitly (see for example discussions in Ovchinnikov et al. (2021) )). First, all methods assume that one has access to the whole trajectory of the observable output, and so are looking at the size of the equivalence classes of the equivalence relation ∼ ∞ on defined as θ ∼ ∞ θ if and only ifŷ(θ, t) =ŷ(θ , t) for all t ≥ 0.
For rational ODE models with time series data as considered here, a result of Sontag (2002) proves if at least 2m + 1 generic time points are observed, where m is the dimension of the parameter space, then the equivalence relation ∼ t 1 ,...,t 2m+1 coincides with the equivalence relation ∼ ∞ . If there are fewer time points or they are not generic, it could be that almost all equivalence classes of ∼ ∞ have size 1 but those of ∼ t 1 ,...,t l are larger. For the Linear ERK model, the parameter space has dimension 3, so we have enough time points, although we do not know a priori if they are generic. In fact, this model admits analytic solutions (see Sect. 4.2), so we can build the model prediction map explicitly and determine its identifiability directly. By a straightforward computation, we can show that for any choice of three distinct non-zero time points, the model prediction map φ t 1 ,t 2 ,t 3 of the Linear ERK model is injective and so the model is globally structurally identifiable (see Appendix A.4 for details). In particular, it follows that any choice of three distinct time points is generic. For the Rational ERK model and the Full ERK model, the parameter space has dimensions 5 and 6, respectively; hence, we may not have enough time points, and we cannot determine the validity of any structural identifiability results for these specific model prediction maps. Indeed, these two models are non-linear and do not admit analytic solutions that would allow us to make the same argument as for the Linear ERK model. This is an instance of a more general open problem: Open Problem 6 Find and implement an algorithm to determine structural identifiability of a rational ODE model with time series data at specific given time points Methods to assess the structural identifiability of ODE models include the classical approach via Taylor series (Pohjanpalo 1978) and generating series (Grewal and Glover 1976), and, more recently, approaches based on differential algebra (Audoly 2001;Saccomani et al. 2003;Hong et al. 2020). In this paper, we use SIAN (Hong et al. 2019), an approach based on differential algebra implemented in Maple (2019).
Similar to other methods based on differential algebra (for example, the method implemented in DAISY (Bellu et al. 2007)), SIAN is based on the differential Nullstellensatz (Ritt 1950, Chapter 1) or (Seidenberg 1952, Section 4). For a differentially closed field K,this theorem establishes a correspondence between radical differential ideals and differentially closed subsets of K n . In the context of an ODE system, this implies that the solutions of the ODE system are completely determined by a prime differential ideal in a differential ring (see below). Criteria for identifiability can then be extracted from the ideal (or the quotient ring). The requirement that K is differentially closed then means that the solutions in question are possibly complex-valued, and the identifiability results will be about complex parameters, whether this is stated explicitly or not. For this reason, Hong et al. (2020) state their definition for complex parameters.
Remark As mentioned above, the first difference between our definition of identifiability and Hong et al.'s is that their parameter space is a subset of C n instead of R n . A second difference to note is that what Hong et al. (2020) call "globally identifiable" corresponds to what we call generically identifiable. Finally, Hong et al.'s (2020) definition is written for components of the parameters and makes the notion of "almost all" more precise.
The starting point is an ODE system of the same form as in Eq. (12) together with the initial condition x(0) = x 0 . Let Q be the least common multiple of all the polynomials appearing in the denominators in f and g. Then we have f = F/Q and g = G/Q where F and G are polynomial functions. Note that SIAN usually views the initial conditions as additional unknown components of the parameter that one may want to identify. The differential ring of interest is the differential ring C(θ ){x, y} (the differential ring in indeterminates x and y over the fraction field C(θ ), i.e. the field of complex rational functions in the parameters). We can think of this ring as a polynomial ring in infinitely many indeterminates: θ , x, y and the infinitely many higher derivatives of x and y (i.e. x (i) and y (i) for i ≥ 1). We are interested in the differential ideal I of C(θ ){x, y} given by where for non-empty subsets T , S of a ring R, the set T : S ∞ is defined as follows: T : S ∞ := {r ∈ R | there exist s ∈ S, n ∈ Z ≥0 such that s n r ∈ T }.
Note that for polynomial systems like the Full ERK model and the Linear ERK model, we have Q = 1, and so the column operation is not needed and the ideal I is simply the differential ideal generated by the equations defining the ODE system and their derivatives. The ideal I is the ideal of all differential polynomials in C(θ ){x, y} that vanish on the solutions of the system of ODE system (12) (Saccomani et al. 2003;Hong et al. 2020).
The ideal I is prime (Hong et al. 2020) and so the quotient ring C(θ ){x, y}/I is an integral domain. Let K := Q(C[θ ]{x, y}/I ) be the field of fractions of the domain C(θ ){x, y}/I , and let k be the subfield of K generated by the image of C{y}, that is, the subfield generated by the elements of the form y i + I . We can now state the non-constructive algebraic criterion for structural identifiability: Proposition 7 (c.f. Proposition 3.4 in Hong et al. 2020) Suppose we have a model given by a system of rational ODEs as described above.
• If the fields k and k(θ ) coincide, then the model is generically identifiable.
• If the field extension k ⊆ k(θ ) is algebraic, then the model is locally identifiable.
Remark Note that Proposition 3.4 in Hong et al. (2020) implies that k and k(θ ) coincide (respectively the field extension k ⊆ k(θ ) is algebraic) if and only if the model is globally identifiable (respectively, locally identifiable) in the sense of Hong et al. (2020). We are interested in something weaker; we only wish to identify parameters in the parameter space , which is a subset of the real positive octant.
The criterion provided by the proposition above is not constructive, as it involves the field of rational functions of an infinitely generated C-algebra. Hong et al. (2020) go on to provide a constructive version of the criterion (Hong et al. 2020, Section 3). The software SIAN (Hong et al. 2019), which we use here, is in turn based on a probabilistic version of the criterion (Hong et al. 2020, Section 4). Note that local identifiability is determined via the Taylor series approach.
We now consider the issue of initial conditions. As mentioned above, by default, SIAN considers the initial conditions as parameters that one may wish to identify. Other methods, like the differential algebra method as implemented in DAISY (Bellu et al. 2007), do not explicitly address initial conditions. Ovchinnikov et al. show in (Ovchinnikov et al. 2021, Theorem 19) that input-output identifiability corresponds to what they call multiple experiment identifiability, that is, identifiability from sufficiently many generic initial conditions. DAISY and COMBOS verify input-output identifiability (Meshkat et al. 2014).
Using SIAN (Hong et al. 2019), we verify that all three models are generically identifiable. In particular, in all three models all parameters are generically globally identifiable. Recall that this result is valid under the assumption that we have measurements at sufficiently many generic time-points, and for generic initial conditions. Inspired by the discussion in Saccomani et al. (2003), in Appendix A.3 we show that the set of differential polynomials in C(θ ){x, y} vanishing on those solutions of the system 12 with initial conditions S 0 (0) = 5μM and S 1 (0) = S 2 (0) = 0μM for all three models, as well as C 1 (0) = C 2 (0) = 0μM and E(0) = 0.65μM for the Full ERK model coincides with the ideal I . This means that the set of solutions with initial conditions corresponding to our experimental setup is dense in the set of all solutions for the Kolchin topology (induced by the differential ideals of C(θ ){x, y}). We can, therefore, conclude that the initial conditions specific to the experimental setup are indeed generic. Therefore, our structural identifiability results hold for the initial conditions specific to the experimental setup.
Remark Using SIAN we can show that the Full ERK model is also generically identifiable with measurable outputs S 0 + C 1 , S 1 + C 2 and S 2 which is what was actually measured experimentally (see Sect. 3.2.4).

Practical Identifiability
Suppose a model is generically identifiable, then, generically, distinct parameters produce distinct data points. However, if there are parameter values that are arbitrarily far from one another but produce data points close to each other, parameter estimation would not be meaningful in practice. Practical (non-)identifiability aims to categorise models exhibiting such undesirable behaviour. For example, sloppiness (Gutenkunst 2007), uncertainty quantification (Smith 2013) and filtering problems (Shi et al. 1999) study mathematical models with a similar aim. We use a definition of practical identifiability introduced in Dufresne et al. (2018), which was adapted from the definition given in Raue et al. (2009).
Practical identifiability depends on more than the defining equations and specification of input and output of the model. Practical identifiability will be influenced by the precise choice of time points, the method used for parameter estimation, the assumption on measurement noise of the data, and the way we measure distances in parameter space. It may also vary on the area in the data space. A data point z * is an experimental observation in the form of an N -dimensional vector whose entries are the observed values of the measured variables at each of the specific time points for each replicate of the experiment. We focus on practical identifiability for maximum likelihood estimation (MLE), one of the most widely used methods for parameter estimation (see, for example Ljung et al. 1987). Accordingly, in the remainder of this section, we consider models (M, φ t 1 ,...,t s , ψ, d ) with a precise choice of model prediction map φ t 1 ,...,t s with specific time points t 1 , . . . , t s , a specific assumption for the probability distribution ψ of measurement noise and a choice of reference metric d on parameter space . We will also assume that the model considered is at least generically identifiable, so that MLE exist and are unique for generic data (see (Dufresne et al. 2018, Proposition 4.15)). We writeθ(z * ) to denote the MLE for z * , that is,θ(z * ) := max θ∈ ψ(θ, z * ).
We define an δ-confidence region U δ (z * ) as follows: The set U δ (z * ), often known as a likelihood-based confidence region (Vajda et al. 1989;Casella and Berger 2002), is intimately connected with the likelihood ratio test. Specifically, suppose we had a null hypothesis H 0 that data point z * has true parameter θ * , and we wished to test the alternative hypothesis H 1 that z * 's true parameter is something else. By definition, a likelihood ratio test would reject the null hypothesis when where k * is a critical value, with the significance level α equal to the probability Pr( (z * ) ≤ k * |H 0 ) of rejecting the null hypothesis when it is in fact true. The set of parameters such that the null hypothesis is not rejected at significance level α is Definition 8 (Dufresne et al. 2018, Definition 4.17) The model (M, φ t 1 ,...,t s , ψ, d ) is practically identifiable for a data point z * ∈ R N at significance level α if and only if the confidence region U δ (z * ) is bounded with respect to d , where For our analysis, we make the common assumption that the measurement noise is additive Gaussian with covariance matrix equal to a multiple of the identity matrix. The assumption is implicit when performing a least-squares fit computation for MLE. In our setup, we are measuring 3 substances at 7 time-points and there were r replicates, so our assumption on the measurement noise means that the probability distribution of the data is given by where σ 2 I 21 is the covariance. It then follows that Therefore, we have that where B ρ (z * ) is the Euclidean open ball of radius ρ := ( z * − φ t 1 ,...,t 7 (θ(z * )) 2 2 − 2σ 2 log k * ) around the data point z * . It follows that under our assumptions, determining whether the various models we study are practically identifiable corresponds to determining whether the preimages under the model prediction map of small open balls around data points are bounded in parameter space. The size of the balls will depend on the data point and the significance level α (or equivalently the critical value k * ).

Algorithm for Testing Practical Identifiability
The Rational ERK model and the Full ERK model do not admit analytic solutions, hence we do not have access to an explicit model prediction map φ t 1 ,...,t l . Therefore, we must approximate φ t 1 ,...,t l and thus also U δ using numerical methods and repeated sampling.
First, we assume that our measurements have been corrupted with some Gaussian noise with mean 0 and variance σ 2 . This variance is identical across measurement quantities, time points and trials. The noise distributions are independent across measurements.
As we have assumed that measurement noise is additive Gaussian with covariance matrix equal to a multiple of the identity matrix, we can obtain an MLE, given some data z * , by solving a least squares problem. This gives usθ(z * ). We use this parameter to calculate the sample variance, assuming that the mean of each quantity is the model trajectory at each time point. This gives us an estimate of the covariance σ 2 .
Recall that δ is defined to be − log ψ(θ(z * ), z * ) − log k * . The log-likelihood is easy to compute, as we already know z * andθ(z * ), and can estimate φ t 1 ,...,t 7 using a numerical solution to the ODE system. We use the following procedure to approximate − log k * : Algorithm 1: Computing approximate − log k * Data: z * ,p(z * ), σ 2 , α, n iter (number of iterations) 1 {z i } i=1,..,n iter ← n iter corruptions of z * by adding i.i.d. random samples from N (0, σ 2 ) to each measurement 2 LogLik ← [] (create an empty list) 3 for i = 1, ..., n iter do 4 Calculatep(z i ) by least-squares solving This simply follows the definition of k * in Eq. (14) and approximates − log k * by repeatedly sampling likelihood-ratios under our given noise assumptions and then taking a (1 − α)-quantile (as − log( · ) is a monotonically decreasing function).

Remark
In a situation where the number of replicates r is large, an approximate δ can be computed from α that depends primarily on the distance between the data point z * and the predicted data point φ t 1 ,...,t 7 (θ(z * )) corresponding to the MLE.
Wilk's theorem (Fan et al. 2000) implies that 2(log ψ(θ(ẑ),ẑ) − log ψ(θ(z * ),ẑ)) is asymptotically χ 2 with three degrees of freedom. If F(ẑ) is the asymptotic cumulative distribution function of 2(log ψ(θ z),ẑ) − log ψ(θ(z * ),ẑ)), then α is approximately equal to Therefore, asymptotically we have that Unfortunately, this is not applicable here, as the number of experiments is 5, 6 or 11, which are not large numbers. Indeed, the δ obtained by applying Wilks' Theorem and the δ obtained via Algorithm 1 are notably different. For example, for the wildtype and the Linear model, we approximate − log k * as 0.477 while Wilks' theorem approximates it as 3.907.
In order to demonstrate practical non-identifiability for the Full and Rational ERK models, we pick two parameters from each model, based on which we can illustrate non-identifiability well by presenting confidence areas marginalised to these two parameters. This choice of parameters is informed by performing a (illposed) Bayesian parameter inference first (see next section). This procedure is described here for the Rational ERK model, but works similarly for the full model: Find parameters (κ 2 , π , γ 2 ) minimising ψ((κ 1 , κ 2 , π ,γ 1 , γ 2 ), z * ) by least-squares solving 9 if − log ψ((κ 1 , κ 2 , π ,γ 1 , γ 2 ), z * ) < δ then 10 Append (κ 1 ,γ 1 ) to C A 11 return C A While we do not know the values of κ 1 and γ 1 , previous experimental work has provided bounds for κ 1 and γ 1 , which we pass to the algorithm above. The list returned by the algorithm is a discrete approximation of the confidence area, marginalised to the pair of parameters κ 1 and γ 1 . We plot these points for visual inspection, which can be seen in Fig. 2. The blue area reaching the upper and leftmost boundary of the plot indicates that the confidence region is very unlikely to be bounded and that this model is very unlikely to be practically identifiable.
The source of this practical non-identifiability of the Full ERK model and the Rational ERK model is not completely clear. One possible source of non-identifiability could be the choice of time points. Indeed, as mentioned in Sect. 4.1, in both cases we do not know if the time points are sufficiently generic. There are reasons to believe that not all practical non-identifiability can be explained by having an insufficient number of time points. Indeed, as part of earlier work during the preparation of Yeung et al. (2019), additional time point data were simulated for the Full ERK model, but confidence regions still appeared unbounded. Another possible source of non-identifiability could be that for the given experimental data there is a valid quasi-steady-state approximation resulting in a smaller dimensional parameter space. At quasi-steady-state parameter values, the reduction is exact and so for these parameters, the equivalence class of ∼ t 1 ,...,t 7 is positive dimensional. Intuitively, since the solutions of the Full ERK model and the Rational ERK model are close to those of the Linear ERK model near quasi-steady-state parameter values, the confidence regions should contain the equivalence class of the nearby quasi-steady-state parameter value, which in this case, was unbounded. This might be an example of more widespread phenomena.

The Practical Identifiability of the Linear ERK Model
We now consider the practical identifiability of the Linear ERK model. What distinguishes the Linear ERK model from the Full ERK model and the Rational ERK model is that an analytic solution to the ODE system is available and so we can construct an explicit model prediction map. The solution to the ODE system (10) with initial conditions S 0 (0) = 5μM and S 1 (0) = S 2 (0) = 0 is given by: As we did for the Rational ERK model in Sect. 4.3, for a given data point z * , we obtain an MLEθ(z * ) by solving a least-squares problem. We then use Algorithm 1 to approximate − log k * , and then δ, using the explicit model prediction map we construct based on the analytic solutions. In Fig. 3 we plot the boundary of the confidence regions at significance level α = 0.05 for the data points corresponding to the wild-type and each mutant. All five confidence regions are seen to be bounded, and we conclude that the model is practically identifiable for those data points.
the Bayesian approach for inferring parameters of the Linear ERK model, as already computed by Yeung et al. (2019). We then introduce topological data analysis (TDA) and present previous results that enable us to analyse the parameters sampled from the posteriors of the wild-type (WT) and four mutants with topological data analysis. Specifically we exploit a theorem by Bobrowski et al. (2017) for hypothesis testing of topological distances in noisy settings. We implement their theoretical result and compare the topological distances between WT and mutants.

Bayesian Inference
Given experimental data and a mathematical model, we seek to infer parameters for which the model accurately fits the data. We choose to do this via Bayesian inference. The theory of Bayesian statistics captures how our belief in the true values of these parameters changes when we make observations (in this case: measurements) in the language of probability theory. Most importantly, Bayesian inference does not infer a single value for each parameter, as would a frequentist approach; rather, it infers a probability distribution of parameter values expressing how strongly we believe a certain set of parameter values is correct.
Formally, we are given a parameter space and observations x from some sample space X . Combining the mathematical model with noise assumptions on available measurements, we obtain an expression for p(x|θ), the likelihood of observing x assuming that the parameter of the model is θ ∈ . In addition, we need to specify a measure of belief in the parameter values before we observe any data, expressed through a probability density p(θ ), called the prior distribution. Theoretically, we want to inform a Bayesian inference only through observations. Consequently, we do not want to inform the inference by placing strong prior beliefs on certain parameter values. In practice, however, a trade-off between neutral prior beliefs (which should only account for substantive prior knowledge and possibly scientific conjectures), analytical convenience and computational tractability is commonplace (Gelman and Shalizi 2013, 11-12).
Having selected a mathematical model and a prior distribution, our formal belief in parameter values becomes by making observations x ∈ X . The probability density p(θ |x) is called the posterior distribution. The proportionality in the above equation indicates that we omitted a normalisation which is independent of θ . As one can approximately sample from p(θ |x) without normalising, the normalisation factor is not necessary for our application.
For the Linear ERK model (Eqs. (10)), the parameter is θ = (κ 1 , κ 2 , π, σ ) ∈ R 4 = . Here, the first three components come from the parameter of the Linear ERK model, while σ , the variance of the distribution of the data, which must be inferred in order to construct a Bayesian model, and will be subsequently marginalised (i.e. integrated out). The observations are measurements of S 0 , S 1 and S 2 . As measurements of each MEK type are taken from r replicates, at 7 different times, for 3 phosphorylation states of substrate, we formally have X = R r ·3·7 = R r ·21 . We have r = 11 for the wild-type, r = 6 for SSDD and r = 5 for all other variants.
Given samples S * 0 , S * 1 and S * 2 , we assume that where t denotes the respective measurement time and i indexes the sample. Here, S j (κ 1 , κ 2 , π, t) is a solution to the ODE system at time t for parameters κ 1 , κ 2 and π . For the Linear ERK model, we can construct an analytic solution to the governing equations, but generally, a numerical solution suffices. Such ODE solutions give rise to an expression for the likelihood p(x|θ). We note that in the above Bayesian model, some standard simplifying assumptions were made. First, in the given setup, negative values of measurements of S 0 , S 1 and S 2 have strictly positive likelihoods, which is not true in reality. Second, we assume that (S * 0 ) t,i , (S * 1 ) t,i and (S * 2 ) t,i are independent random variables for all t and i and that they have the same standard deviation. Despite these assumptions, we obtained good fits to the data. For example, performing an inference with three different standard deviation parameters σ 0 , σ 1 and σ 2 for S 0 , S 1 and S 2 , respectively, did not significantly improve the fits to the data. This Bayesian inference framework can also be applied to other ODE models describing the measurements, including the Rational ERK model (Eqs. (9)) and the Full ERK model (Eqs. (1)). In these cases, we employ numerical solutions and adapt priors to the larger parameter spaces.
We note that for the Full ERK model and the Rational ERK model, the choice of prior distributions significantly changes both the location and prominence of modes of the posterior distributions. In particular, they tend to be near the endpoints of the prior distributions. This is linked to the practical non-identifiability of these models and prevents us from interpreting parameter modes, and also from conducting a sensible topological comparison that is not highly dependent on the choice of prior distribution.
In order to compute posterior distributions of the involved parameters, we used PyStan, the Python version of the statistical software STAN (Carpenter et al. 2017). While analytical expressions for the posterior distributions are too complex to be feasible for interpretation, PyStan enables us to approximately sample from them via Hamiltonian MCMC. The resulting samples (visualised in Fig. 4) form the basis of our further analysis.

Topological Analysis
To analyse the topology of the samples of the resulting posterior distributions, we introduce notation and methodology from Topological Data Analysis (TDA).

Definition 9
Let v be a finite set of vertices. A subset of the power-set of v, K ⊆ P(v), is called a simplicial complex if for any τ ∈ K the relation τ ⊆ τ implies τ ∈ K.
We write K i = {τ ∈ K | |τ | = i + 1} and call the elements of K i the i-simplices. A map h : v → v which extends to a map h : K → K by h(τ ) := {h(v) | v ∈ τ } for each τ ∈ K is called a simplicial map. We can view a simplicial complex as a combinatorial description of a topological space. Given a simplicial complex K, we can investigate its geometric realisation where cvx denotes the convex hull in the real free vector space generated by the vertices V . The realisation | K | is endowed with the subspace topology in R v . An example of a simplicial complex and its geometric realisation can be found in Fig. 5a. Since K is a discrete and combinatorial entity, one can compute meaningful topological information from topological spaces (or datasets) described by simplicial complexes.

Homology
One topological invariant we can compute from simplicial complexes is homology. In each dimension k, the dimension of the k-th homology group can be thought of as the number of voids in a simplicial complex enclosed by a k-dimensional boundary. We restrict our definition of homology over the field of two elements, F 2 , which is the setting for our computations. For a simplicial complex, the homology groups coincide with those of its geometric realisation (viewed as a topological space).
Definition 10 Let K be a simplicial complex. We define its chain complex C • (K) over F 2 to be the collection of vector spaces C i = F 2 K i , together with the collection of linear maps ∂ i : C i → C i−1 induced by We observe that ∂ i • ∂ i+1 = 0 for all i. Furthermore, we note that any simplicial map h : K → K induces a collection of maps on corresponding chain complexes C • and C • , denoted {ĥ i : We call such a collection of maps a chain map from Definition 11 Let K be a simplicial complex and let C • (K) be its associated chain complex over F 2 . Then the k-th homology group of K is defined to be the quotient of vector spaces Note that forĥ : where c ∈ ker ∂ k and the brackets denote equivalence up to translation by im ∂ k and im ∂ k respectively, is well defined for all k (Otter 2017). Moreover, for simplicial maps h : K → K and h : This property is called the functorality of homology and will be used when we introduce persistence.

Persistence
We view point clouds as a discrete subset of a continuous geometric object embedded in Euclidean space. The underlying continuous space is the primary subject of interest. In order to obtain information about this geometric object, we wish to inflate our discrete points to a continuous space, or to capture a relative offset between points in this space. In practice, we usually do not know the adequate inflation resolution. Persistence theory offers an elegant way to overcome this caveat by scaling the resolution from fine to coarse, and tracking how the homology of these spaces evolves by considering their canonical inclusion relations.

Definition 12
Let K be a simplicial complex and let g : K → R be a function such that τ ⊆ τ implies g(τ ) ≤ g(τ ) for any τ, τ ∈ K. A filtration of the simplicial complex K by g is then defined to be the sequence of simplicial complexes {K L } L∈R , where together with the canonical inclusions ι L L : K L → K L whenever L ≤ L . An example of a filtration is visualised in Fig. 5b. In the same spirit, let T be a topological space and g : T → R be a continuous function. A filtration of the topological space T is then defined to be the sequence of topological spaces {T L } L∈R , where together with the canonical inclusions ι L L : T L → T L whenever L ≤ L . A common way of constructing a filtration from a point cloud v ⊂ R d is to set K = P(X ) and g(τ ) = max{d(x, y) | x, y ∈ τ }. This is called the Vietoris-Rips filtration, and K L is a good approximation to an inflation of v by placing balls of radius L/2 at each point (Oudot 2015). We will consider the following alternative filtration. For a fixed L ∈ R and map p : R d → R, we set K := K L in the Vietoris-Rips sense and consider the filtration by the map g : Definition 13 Let F 2 [t] be the ring of polynomials in the indeterminate t with coefficients in F 2 . Let {K L } L∈R be a filtration of a simplicial complex. Moreover, define Crit L := {L ∈ R |ι L L−ε = id ∀ε > 0}, the set of all L at which K L changes (which is a finite set at K is finite). Define the function c : N 0 → Crit L ∪ {inf Crit L − 1} by mapping 0 to inf Crit L − 1 and n > 0 to the n-th smallest element of Crit L (without loss of generality, we map integers bigger than the cardinality of Crit L to the largest element of Crit L ).
For a fixed integer k, let H k ( · ) denote the k-th simplicial homology with coefficients in F 2 . Define together with the action of F 2 [t] on M k induced by t a · x = ι c(n) c(n+a) (x) * ∈ H k (K c(n+a) ) for x ∈ H k (K c(n) ) and non-negative integer a. Then M k is a (graded) F 2 [t]-module, called the persistence module of the filtration.
The definition works analogously for a filtration of a topological space (assuming that the homology of the spaces changes at only finitely many filtration values). It can be shown that the operation of taking a persistence module of a filtration of a simplicial complex (or a topological space) is functorial. Hence, persistence modules are algebraic invariants of filtrations.
Since K is finite, the persistence module M k is finitely generated as a F 2 [t]-module. As F 2 [t] is a principal ideal domain, M k decomposes into summands generated by a single object uniquely up to (graded) isomorphism and permutation of summands. Hence, we can write where G F is the subset of chosen generators that are free and G T is the subset of generators that are torsion. In particular, any element in G F or G T will have a non-zero entry in exactly one summand of the decomposition in Eq. (15). We call the integer n indexing this entry the degree of that element.

Definition 14
Let M k be a persistence module that decomposes as above. Let deg : G F ∪ G T → N 0 be the function mapping each element to its degree. The barcode of M k is defined to be the multiset deg(a)), c(deg(a) + d a )) | a ∈ G T }.
We call the elements of B bars, the first coordinate of each bar its birth-value, the latter coordinate its death-value and the absolute difference of the coordinates its persistence.
A matching of barcodes B and B is a partial injection : B → B . The bottleneck distance between B and B is defined to be where the infimum is taken over all possible matchings and elements of a barcode are viewed as elements of R 2 (we assume ∞ − ∞ = 0). Here, dom is the domain of , i.e. the set of inputs at which is defined.
The bottleneck distance defines a metric on the space of barcodes (Oudot 2015). This metric is stable in the following sense: Theorem 15 (e.g. Corollary 3.6 in Oudot 2015) Let K be a simplicial complex and let g, g : K → R be functions defining filtrations of K, and subsequently persistence modules M k and M k , and barcodes B and B . Then Henceforth, we write PH k (g) to denote the k-dimensional persistent homology (which can equivalently be summarised by a barcode or a persistence module) of a simplicial complex or a topological space filtered by a function g.

Persistent Homology of Random Data
In this section, we study the persistent homology of the posterior distributions of the parameter inferences of Sect. 5.1. Note that simplicial complexes, filtrations and persistent homology can also be employed to compare biological models a priori (i.e. with no dependence on measurement data) (Vittadello and Stumpf 2020).
We demonstrate that filtering a Vietoris-Rips complex for a fixed value L by a function g , as described at the beginning of this section, yields more discriminative power. Here, we pick g to be an estimated probability density function. These filtrations turn out to be highly discriminative between the mutants and offer novel insight at the biological level. While a Vietoris-Rips filtration is entirely based on distances, the construction we employ, using a Vietoris-Rips complex at a fixed parameter value and then filtering it by a probability density function (pdf), places an emphasis on density. The information encoded is directly related to the probability distribution and the resulting barcodes will stabilise as the sample size increases (Theorem 3.5.1 in Rabadan and Blumberg 2020). Furthermore, the chosen construction is stable with respect to outliers. By contrast, in a Vietoris-Rips filtration, bars in the resulting barcodes will converge towards zero length when increasing the sample size and a single outlier, even in a large sample, can change a barcode drastically.
Initially, we assume that we are given a probability density function p : R m → R. This pdf defines a filtration of the graph by − p, T say, via Such a filtration is visualised for the case m = 1 in Fig. 6. By analogy with filtrations of simplicial complexes, we can theoretically compute a barcode for each such topological filtration and investigate the resulting bottleneck distances.
For each (homological) dimension, these barcodes provide a topological signature of a posterior distribution. We point out that although this signature is not a sufficient statistic, it is effective at distinguishing between posteriors corresponding to distinct mutants in our application. In particular, for any pdf p 1 : R d → R, the pdf p 2 (x) = p 1 (x − x 0 ) gives rise to the same topological signature for any constant x 0 ∈ R d . Thus, rather than comparing the location of probability density in parameter space, in the context of a Bayesian inference, this topological signature captures the quality of the certainty we have in parameter values, irrespective of their location.
For example, bars in the H 0 -barcode encode the density (as negative of the birthvalue) and the prominence (as the persistence) of the modes of a pdf. Similarly, Morse Theory tells us that for a (smooth) pdf on R d , the (d − 1)th barcode captures local minima by their density (as death-value) and the depth of their basin of attraction (as persistence).
In order to conduct such a topological analysis, two questions must be addressed: (1) How can we approximate the topology of a graph of a probability density combinatorially (i.e. in a manner amenable to the application of discrete computational methods) if only point samples are available? (2) Can we test the statistical significance of the resulting bottleneck distances?
To resolve the first question, we will employ a result from Bobrowski et al. (2017) that relies on the concept of kernel density estimation (KDE). In order to test the significance of the resulting bottleneck distance, we will use an empirical p-value estimate.

Definition 16
Let v = {v 1 , . . . , v N } ⊆ R m be a set of N samples drawn independently from a probability distribution governed by the density function p : R m → R. Let K : R m → R be smooth, unimodal, symmetric probability density function whose support is contained in the unit ball centred at 0. Then is called a kernel density estimate (KDE) of p with bandwidth b. On each sample v i , we place a pdf and average it, where b controls the width of each pdf, that is, how much of the probability mass is centred around v i . Loosely speaking, if b is too large, then the resulting function underfits a histogram given by the data, while if it is too small, then the bandwidth overfits the histograms (see Fig. 7). The bandwidth is negatively correlated with the sample size and there are standardised ways of picking optimal bandwidths for the case where p is unknown (Henderson and Parmeter 2012).
Given such an i.i.d. sample v = {v 1 , . . . , v N } ⊆ R m from our probability density function p and an optimal bandwidth b, we can construct a Vietoris-Rips complex with fixed parameter b (equalling the bandwidth) For the sake of brevity, let K = VR b (v). The KDEp b of p based on v then extends to a function on K viâ In return, the extended functionp r defines a filtration {K L } L∈R of K by We seek to relate the persistent homology of the filtration of simplicial complexes K L to the persistent homology of the filtration of topological spaces T L .
In order to use results from Bobrowski et al. (2017), we introduce some notation.
Theoretically, the above theorem can be exploited for testing the null hypothesis H 0 : PH k ( p) = PH k p for two distributions P and P with associated densities p and p , as the result enables us to establish a bound on how large a bottleneck distance can be explained by sampling noise at a given significance level. However, we estimate that to use this theorem for showing that the bottleneck distances between posterior distributions associated with the wild-type and the four mutants are significant, we must sample at least 1.5×10 7 points per distribution. This makes persistent homology computation infeasible.
At the same time, we observe that there is little change in the bottleneck distances between the barcodes resulting from the wild-type's and the four mutants' posterior distributions when resampling point clouds containing as few as 200 points. This leads us to think that the true p-value associated with the null hypothesis H 0 : PH k ( p) = PH k p , where p and p are posterior densities corresponding to the wild-type and a mutant is possibly much lower than the upper bound derived by appealing to Theorem 17. One factor that may explain this discrepancy is that while our distributions are technically distributions on R 3 , they have compact support. Similarly, major sources of instability for KDE, and subsequently for the filtration of density functions, are modes linked to outliers, while repeated simulations suggest that in our case all density functions are unimodal. Together, these aspects imply that the computed barcodes could converge to the barcode obtained by filtering the unknown density function at a faster rate than in the general setting of Theorem 17. Henceforth, we use the method of constructing a filtration based on a point cloud proposed in Bobrowski et al. (2017), which is provably well-behaved asymptotically but uses a different approach to estimate significance. To do this we opt for a Monte Carlo p-value estimate, also known as the empirical p-value (e.g. see Davison and Hinkley 1997).
is a p-value estimate for a hypothesis test H 0 : PH 1 ( p) = PH 1 p . The resulting p-value estimates, for each pair of mutants and wild-type, can be found in Table 1. It is likely that these p-value estimates over-estimate the actual value, but they allow us to reject all null hypotheses at a significance level of 0.05 (North et al. 2002).
The results (Table 1) of the topological data analysis quantify the differences between the Linear ERK model parameter posteriors for WT and mutants and find SSDD mutant kinetics are most different from WT and other mutants. This biological result raises the suitability for using the SSDD variant as a replacement for wild-type MEK activated by Raf. We suggest this should be investigated with further experimen-tal studies. The previous work by Yeung et al. (2019) found that π , the processivity parameter, of E203K differed the most from WT MEK. Here we extended and complemented their analysis by comparing the three parameters together as a point cloud.
It remains to address the practical computability of all the constructions involved. As mentioned in the previous section, we use the statistical software STAN (in particular, PyStan) to sample from the posterior distributions. This sampling is approximate via Hamiltonian MCMC (Carpenter et al. 2017), but we can verify via output summaries and trace plots of the Markov chains involved that all chains have converged close to their stationary distribution during the warm-up phase.
In order to construct the KDE, we used the KernelDensity method of the Python package sklearn. We used the Epanechnikov kernel, which satisfies the hypothesis of the kernel in Theorem 17. As a guess for the bandwidth, this package uses a rule of thumb proportional to Silverman's method, which we then cross-validate and plot against a histogram of our samples for each marginal distribution. Given experimental data, we construct a Vietoris-Rips complex with a radius b, equalling the bandwidth from the KDE, using the Python package dionysus (version 2). We compute the resulting bottleneck distances using the package persim.

Conclusion
We presented an exhaustive mathematical analysis that supports the three main findings presented in Yeung et al. (2019): model reduction, analysis of the model parameters and comparing mutation kinetics. Yeung et al. observed that certain values of parameter combinations from the Full ERK model fit the data, which in turn motivated the creation of a reduced model, the Linear ERK model. We confirmed the derivation of the Linear ERK model using algebraic QSS and the validity of the QSS approximation using the QSS variety. We performed systematic identifiability analyses on all three models, which is a prerequisite for meaningful parameter estimation. We found the Full, Rational and Linear ERK models are structurally identifiable. We then improved a previous definition of practical identifiability and showed that the Linear ERK model is practically identifiability but Rational and Full ERK models are not, which is consistent with (Yeung et al. 2019). Hitherto, testing structural identifiability has been limited to small models due to computational costs; however, recent work significantly improves computing structural identifiability, enabling analysis of larger models (Dong et al. 2021;Villaverde et al. 2019). We remark there are many realistic models, such as this ERK study or those by the group of Marisa Eisenberg, that benefit from existing methods and motivate the development of new identifiability tools.
We reproduced the parameter inference for wild-type and mutant MEK experiments. While Yeung et al visually inspected samples of the posteriors, here we quantified these point clouds with computational algebraic topology. In future, it would be interesting to further explore the relationship between topological analysis and practical identifiability and how they may be used to inform experimental design (Apgar et al. 2010;Hagen et al. 2013). Throughout we showcase the potential role of algebra, geometry and topology in systems and synthetic biology. Complementary to the analysis here is an inference of models in systems and single-cell biology that relies on algebra and topology (Wang et al. 2019;Vittadello and Stumpf 2021;Rizvi 2017). We believe that topological data analysis in combination with modelling and parameter estimation is a promising area for the sciences (Thorne et al. 2022;Carriere et al. 2018;Suzuki 2021). We hope our analysis of this ERK case study will motivate other systems biologists to partner with algebraists and topologists to analyse dynamical systems together with their experimental setup and data.
Observe that the Linear ERK model obeys the conservation law S 0 + S 1 + S 2 =S tot . Here,S tot is a constant close to S tot . Hence, E = E tot − (S tot −S tot ) is a constant close to E tot . The implication that E = E tot at all times t is physically infeasible. WhetherS tot = S tot depends on whether one updates the initial conditions for S j through an inner solution, which is typically assumed in singular-perturbation-theory approaches. The algebraic approach does not require computing an inner solution.
It is important to reiterate at this point that this reduction is not the result of a singular perturbation analysis. Indeed, if we had followed a singular perturbation analysis with substitution (8) for E, we could not factor out E tot /S tot in the non-dimensionalised differential equations for C 1 and C 2 as nicely as we were able to in the previous subsection. In this instance, we would have to leave factors of ε −1 in the algebraic equations, which would be ambiguous when taking a limit ε → 0.

A.1.2 Deriving the Rational ERK Model
We substitute E using Eq. (7) into the Full ERK model (Eqs. (1)). We then follow the same steps as in the algebraic QSSA model reduction for the Linear model to arrive at the Rational model.
Alternatively, after substituting for E using Eq. (7), we can perform classical QSSA to arrive at the same model reduction (Keener and Sneyd 2011).

A.2 Accuracy of Algebraic QSSA
We start by providing the full statement of Proposition 2 of Goeke et al. (2017) (i.e. the full version of Proposition 4 of this manuscript): Let K * ⊂ R n + × R m + satisfy the following: • There exist y 0 ∈ R n and r > 0 with the following property: Whenever (x, θ) ∈ K * for some x ∈ R n and some θ ∈ R m These conditions imply that every V θ * , with θ * nearθ , is a submanifold. Note that every (y, θ * ) with y in the interior of R n + is contained in some K * that satisfies the last three of the above conditions. Proposition 18 Assume that the above conditions are satisfied for K * .
(a) Let θ * be given such that V θ * × {θ * } has non-empty intersection with int K * , let (y, θ * ) be a point in this intersection and V θ * ⊂ R n be some open neighbourhood of y such that (V θ * ∩ V θ * ) × {θ * } ⊂ K * . Moreover let t * > 0 such that the solution of (4) with initial value y exists and remains in V θ * for all t ∈ [0, t * ]. Then there exists a compact neighbourhood A θ * ⊆ V θ * of y with the following properties: (i) For every z ∈ A θ * the solution of (4) with initial value z exists and remains in V θ * for all t ∈ [0, t * ]. (ii) For every ε > 0 there is a δ 1 > 0 such that the solution of (6) with initial value z ∈ A θ * ∩ V θ * exists and remains in V θ * for t ∈ [0, t * ] whenever f − f red < δ 1 on V θ * . (iii) For every ε > 0 there is a δ ∈ (0, δ 1 ] such that the difference of the solutions of (4) resp. of (6) with initial value z ∈ A θ * ∩ V θ * has norm less than ε for all t ∈ [0, t * ] whenever f − f red < δ on V θ * . (b) Let y ∈ V θ * and let ρ 0 > 0 such that Let ρ ≤ ρ 0 such that f (y, θ * ) − f red (y, θ * ) ≥ 2ρ. Then for t := ρ/(2L R) the solutions of (4) resp. of (6) with initial value y exist and remain in B ρ 0 /2L (y) for 0 ≤ t ≤ t , and their difference has norm at least ρ 2 /(2L R) at t = t .
Returning to our two model reductions of interest, we use the notation of Sects. 2 and 3 and label the two possible substitutions for the variable E as (7) and (8), as before. Moreover, we writeṠ 0 andṠ 1 as shorthand for the algebraic expressions for the time derivatives of S 0 and S 1 respectively. Recall that K i := k c i + k r i . Throughout, we will view our set of polynomials f , which govern the ODE system, as a smooth function from R 5 to R. Any norm in this subsection will refer to the · ∞ -norm on the respective space.
Define a domain K * such that k f i K i for i = 1, 2 and such that it satisfies all hypotheses in the list at the start of this section (i.e. from Proposition 2 of Goeke et al. 2017).
Note that for both substitutions, E is bounded above by E tot ≈ 0.65μE. Hence, the norms on the slow variables are going to be approximately the same for both substitutions.
For the case (7), we get that f red (showing only the components that have changed compared to f ) is given by on variables C 1 and C 2 . Here, the quantity det is given by For the quantity f − f red , for substitution (7), we get Similarly, for (8) we get Both of the above norms are O(1) (asṠ 0 andṠ 1 contain terms K 1 and K 2 respectively). Hence, all quantities use to bound accuracy in Proposition 18 (above, by part (a), or below, by part (b)) are of a similar order. We conclude that the accuracy of both reductions is approximately the same when assuming 0 ≤ k f i K i . This follows from the measurements k M i ≈ 25 μM (Bar-Even 2011). Moreover, we have seen previously (in the main text) that K i > 0 for i = 1, 2 is a sufficient condition for a parameter value θ to be a QSS-parameter value in the sense of Goeke et al. (2017).

A.3 Deriving Equality of Ideals
In this section of the appendix, we derive that I (V 0 ) = I for the ERK models with given initial condition. In Saccomani et al. (2003), V 0 is defined as all trajectories of a given ODE system with initial condition x 0 . That is, if x(t) satisfies S(0) = x 0 anḋ S = f (S, θ), then V 0 = {x(t) | t ≥ 0} ⊂ C(θ ) n . Then I (V 0 ) is defined to be the ideal of all polynomials in C(θ )[S] vanishing on V 0 .

Proposition 19
Assume we are given a complex-valued ODE model, including an initial condition. Let θ denote the set of its parameters and S the set of its variables. Consider the affine space A given by its variables and their derivative terms (of all orders), viewed as an affine space over the fraction field C(θ ). Define V 0 to be the manifold in A given by the ODE trajectories under the given initial condition.
Then for the Full, Rational and Linear ODE model considered in this manuscript, we have that I (V 0 ), the differential ideal of all polynomials in the differential ring C(θ )[S] vanishing on V 0 (c.f. Saccomani et al. 2003), is exactly I .
Proof By definition, I (V 0 ) ⊃ I . Hence, we need only show the inclusion I ⊂ I (V 0 ). By way of contradiction, suppose there are elements of I (V 0 ) which are not elements of I . Without loss of generality, we may assume that such polynomials contain no derivatives of variables. Indeed, starting with any polynomial in I (V 0 ), we can use the differential equations in I to cancel out all terms containing derivatives of variables and so obtain an element of I (V 0 ) without such terms.
We need to take some extra care for the Rational ERK model: strictly speaking, the polynomial A(γ 1 S 0 + γ 2 S 1 + 1) − 1 is in I (V 0 ) but not in I . However, we know this relation a priori and it is merely a result of us converting a rational ODE model into a polynomial one. In particular, it is independent of an initial condition. Hence, in the context of this section, we will assume that the Rational ERK model contains an output variable y 3 := A(γ 1 S 0 + γ 2 S 1 + 1) (in addition to y 1 = S 0 and y 2 = S 1 ). Whilst y 3 adds polynomials to I , we will omit it from any polynomials we study henceforth as it is equivalent to the constant 1 (given the definition of A).
For the Rational ERK model, assume there exists a nonzero p ∈ C(θ )[S 0 , S 1 , A] such that p(S 0 , S 1 , A) = 0 for all t ≥ 0 given our initial condition. By using A(γ 1 S 0 + γ 2 S 1 +1) = 1, we can turn p into a rational function in variables S 0 , S 1 which vanishes on the ODE trajectories. Hence, there must exist a polynomial q ∈ C(θ )[S 0 , S 1 ] such that q(S 0 , S 1 ) = 0 for all t ≥ 0 given our initial condition. Any properties we will use in the following argument are satisfied by both the Rational and Linear ERK models unless stated otherwise. As we need to disprove the existence of such polynomial q when studying the Linear ERK model, the following conclusion will hold for both the Linear and Rational ERK models.
Without loss of generality, assume that q is a polynomial of the smallest non-zero degree satisfying q(S 0 , S 1 ) = 0 on the ODE trajectories and that q is irreducible. If it were reducible, simply replace q with an irreducible factor of q such that q(S 0 , S 1 ) = 0 for all 0 ≤ t < t 1 for some t 1 > 0. The following argument will still hold.
As this change of basis is an invertible affine transformation, we may assume that q is a polynomial in variables X 0 and X 1 . We note that this diagonalisation of the ODE models would not have been possible if we had not removed S 2 from our systems, as S 2 decouples.
We will write (X 0 , X 1 ) = i, j≥0 a i j X i 0 X j 1 = 0.
Taking the derivative of this equation with respect to t gives i, j≥0 a i j (iẊ 0 X i−1 0 X j 1 + jẊ 1 X i 0 X j−1 1 ) = 0, which, after substitution for the derivative variables, yields q (X 0 , X 1 ) := i, j≥0 a i j (iκ 1 + jκ 2 )X i 0 X j 1 = 0 (for the Rational ERK model, we need to divide by A to obtain such q ). Note that, given our initial condition, both X 0 and X 1 vary, hence the intersection of V (q) and V (q ) must contain a smooth point, as both contain the ODE trajectory. This would imply that q is a constant multiple of q, which is not true at generic parameter values. Hence, no such q can exist.
For the Full ERK model, assume that there is a non-zero polynomial p ∈ C(θ )[S 0 , S 1 , C 1 , C 2 ] such that p(S 0 , S 1 , C 1 , C 2 ) = 0 for all t ≥ 0 given our initial condition. Let p have degree n. Again, we may assume without loss of generality that n is minimal. We know that there are no conserved quantities in the Full ERK model (after removing E and S 2 ) and thus n ≥ 2. Let p(S 0 , S 1 , C 1 , C 2 ) = i, j,k,l≥0 a i jkl S i 0 S j 1 C k 1 C l 2 = 0.
By taking the derivative with respect to t, we get i, j,k,l≥0 For S 0 , this rearranges tȯ (assuming there are S 0 terms in p) and we can derive similar expressions for the other variables. Note that the denominator is of a smaller degree than p, hence it will be non-zero at almost all points of our ODE trajectory. In addition, for our given initial condition, all of our variables vary and hence the numerator must be a non-zero polynomial and must not vanish on the ODE trajectory. Denote i * the highest power of S 0 in p. Then a i * jkl = 0 implies k = 0, as otherwise the highest power of S 0 in the numerator is i * + 1 and i * − 1 in the denominator, implying that we can find an S 2 0 -term inṠ 0 , a contradiction. We can derive similar statements for S 1 , C 1 , and C 2 . Then, in the above equation, the highest power of C 1 in the numerator is k * + 1, while the highest power in the denominator is k * − 1 (at generic parameter values), implying thatṠ 0 is quadratic in C 1 , a contradiction. If p does not contain terms in S 0 , we can apply a similar argument to S 1 , C 1 or C 2 .
In conclusion, no such p can exist.

A.4 The Linear ERK Model is Globally Structurally Identifiable
Proposition 20 For any choice of three time points 0 < t 1 < t 2 < t 3 the model prediction map φ t 1 ,t 2 ,t 3 is injective and so the Linear ERK model is globally structurally identifiable.
As the derivative of Eq. (17) has exactly one solution in t, Eq. (17) has at most two solutions in t by Rolle's theorem. Therefore, the second, fifth and eighth components of φ t 1 ,t 2 ,t 3 cannot take the same value. This is a contradiction, and so we must have (κ 1 , κ 2 , π) = (κ 1 , κ 2 , π ).