Identifiability and numerical algebraic geometry

Summary

A common problem when analyzing models, such as mathematical modeling of a biological process, is to determine if the unknown parameters of the model can be determined from given input-output data. Identifiable models are models such that the unknown parameters can be determined to have a finite number of values given input-output data. The total number of such values over the complex numbers is called the identifiability degree of the model. Unidentifiable models are models such that the unknown parameters can have an infinite number of values given input-output data. For unidentifiable models, a set of identifiable functions of the parameters are sought so that the model can be reparametrized in terms of these functions yielding an identifiable model. In this work, we use numerical algebraic geometry to determine if a model given by polynomial or rational ordinary differential equations is identifiable or unidentifiable. For identifiable models, we present a novel approach to compute the identifiability degree. For unidentifiable models, we present a novel numerical differential algebra technique aimed at computing a set of algebraically independent identifiable functions. Several examples are used to demonstrate the new techniques.

Keywords:

Blood – Polynomials – Interpolation – Computing methods – Vector spaces – Algebraic geometry – Complex numbers – Compartment models

Introduction

Parameter identifiability analysis for dynamical system models consisting of ordinary differential equations (ODEs) addresses the question of which unknown parameters can be determined from given input-output data. In this paper, we address structural identifiability, which concerns whether the parameters of a model can be determined from perfect input-output data, i.e., noise-free and of any time duration required. This is a necessary condition for the practical or numerical identifiability problem, which involves parameter estimation with real, and often noisy, data. For this reason, structural identifiability is often referred to as a priori identifiability [1]. Even if a model fails to be structurally identifiable, some useful information about the parameters can still be determined, which is the main motivation for this paper.

There are two possible outcomes of the structural identifiability check of a mathematical model. If the parameters of the model have a unique or finite number of values given input-output data, then the model and its parameters are said to be identifiable. However, if some subset of the parameters can take on an infinite number of values and yet yield the same input-output data, then the model and this subset of parameters are called unidentifiable. In the latter case, we attempt to find a set of identifiable functions of the parameters. These can then be used to reparameterize the model and also to give additional insight into which parameters should be experimentally measured [2].

Several methods have been proposed to find identifiable functions. In linear models, this can be done using the transfer function method [3]. However, in nonlinear models, the problem has been more challenging with only ad hoc methods proposed, e.g., [2, 4, 5]. For example, the approach in [2] requires the calculation of many Gröbner bases and can thus be computationally expensive. It should be noted, however, that even in the linear case, the identifiable functions of parameters found using the transfer function method are not necessarily (and are usually not) the simplest identifiable functions of parameters. Since our goal is to reparametrize a model over identifiable functions of the parameters, simpler functions are preferred.

In this paper, we use techniques from numerical algebraic geometry (e.g., see [6, 7] for a general overview) to investigate both identifiable and unidentifiable models. For an identifiable model, we compute the finite number of values of the parameters given input-output data. The total number of such values over the complex numbers is called the identifiability degree which is computed in two ways. The first method relies on differential algebra tools to first generate the input-output equations while the second does not utilize these equations.

For unidentifiable models, we also introduce two novel approaches for finding identifiable functions of the parameters. The first method relies on knowing the input-output equations and uses them to find globally identifiable functions of parameters, as in [2]. In the case where these input-output equations cannot be calculated using conventional differential algebra techniques, we also introduce a method to compute locally identifiable functions of parameters. This combination of numerical algebraic geometry and differential algebra could be thought of as numerical differential algebra. We demonstrate our methods on various models.

Materials and methods

Identifiability

We consider ODE models of the form:

Input-output equations

One approach to determine identifiability properties of the model (1) using known input-output data is via the input-output equations, i.e., equations that relate the input u(t), output y(t), and parameters p. Thus, the input-output equations result from eliminating the state variables x(t). Several methods have already been proposed, e.g., [5, 8–16], to compute the input-output equations, including the so-called differential algebra approach to identifiability [11, 13, 15]. Using differential algebra, the state variables x(t) are eliminated using differential elimination techniques. If the number of outputs y(t) is m, this procedure produces m differential polynomial equations that are solely in input and output variables with rational coefficients in the parameters so that the j^th one can be written as

Deciding identifiability

Let m₁ denote the number of independent parameters p and m₂ denote the total number of non-constant coefficients taken from all m input-output equations. Thus, we can treat the coefficients of the input-output equations as a rational map c : C m 1 → C m 2. Identifiability refers to whether it is possible to recover the parameters of the model only by observing the relations among the input and output variables. In other words, assuming known input-output data for a sufficient number of time instances so that c can theoretically be computed, identfiability asks whether it is possible to recover the parameters p.

Definition 1. Let c be the coefficients of the input-output equations for a model (1). For general p ∈ C m 1, let

ℓ = dim X_p ≥ 0, and k = # X p ∈ N ∪ { ∞ }. That is, ℓ is the dimension of a general fiber of c and c is generically a k-to-one map when ℓ = 0. The model (1) is identifiable from c if ℓ = 0, i.e., k ∈ N, and unidentifiable if ℓ > 0, i.e., k = ∞.

When identifiable, the number k ∈ N is called the identifiability degree. If k = 1, the model (1) is called globally identifiable and called locally identifiable if 1 < k < ∞.

When unidentifiable, the number ℓ ≥ 1 is called the dimension of unidentifiability.

To distinguish between identifiable and unidentifiable models, one simply needs to compute the dimension ℓ of a general fiber of c. As defined in Section 13.4 of [7], the rank of c, denoted rank c, is the rank of the Jacobian matrix of c evaluated at a general, i.e., random, p ∈ C m 1. The corank of c is corank c = m₁ − rank c. The following, which is Lemma 13.4.1 of [7] (see also [17]), relates ℓ and corank c.

Proposition 2. For a general p ∈ C m 1 , ℓ = dim X_p as defined in (2) is equal to corank c where c is the set of coefficients of the input-output equations. In particular, the model (1) is identifiable if and only if c has full rank and the dimension of unidentifiability is equal to corank c.

In particular, Prop. 2 indicates a method to distinguish between identifiable and unidentifiable models provided that the coefficients c of the input-output equations can be computed, which is summarized in the following pseudocode.

Method 1: Computing dimension of unidentifiability from input-output equations

Input: m₂ input-output equation coefficients c(p), depending on parameters p = ( p 1 , … , p m 1 ).

Output: Dimension of unidentifibility ℓ = corank c = dim c⁻¹(c(q)) for general q ∈ C m 1.

Choose random, complex values q ∈ C m 1.

Return ℓ = corank J_c(q) where Jc(p) is the Jacobian matrix of c evaluated at p.

Example 3. Linear compartment models are frequently used models arising in pharmacokinetics, toxicology, cell biology, physiology, and ecology [18–22]. The following from [17] is an example of a linear three-compartment model with input u(t), output y(t), state variables x(t) = (x₁(t), x₂(t), x₃(t)), and unknown parameters p = (k₀₁, k₀₂, k₀₃, k₁₂, k₁₃, k₂₁, k₃₂), where k_ij represents the rate of transfer from compartment j to compartment i and k_0i represents a leak from compartment i to outside the system:

Fig 1 presents a pictorial representation of this model.

The approach described in [17, 23] yields the input-output equation:

For an identifiable model, one approach to distinguish between global and local identifiability is to solve the system of equations c(q) = c(p) given a general point p ∈ C m 1. If there is a unique solution, namely q = p, the model is globally identifiable. If there are a finite number of solutions, the model is locally identifiable. Such an approach, for example, is implemented in the software package DAISY [1, 24] which randomly selects a point p and uses Gröbner basis methods to count the number of solutions to c(q) = c(p) yielding the identifiability degree. Since such an approach can only be applied when c has first been computed, we will consider the following problem using numerical algebraic geometric methods.

Problem 4. Given a model (1), decide if it is identifiable or unidentifiable. If identifiable, determine its identifiability degree to decide if it is globally identifiable or locally identifiable.

One technique for determining whether a model is identifiable without computing c is via the Exact Arithmetic Rank (EAR) approach [25]. In particular, rather than eliminating to compute the corank of c, one considers projections of a system that still involves the state variables derived by replacing functions with Taylor series expansions and taking a finite-size system via the Cartan-Kuranishi Theorem that underlies differential elimination, e.g., see [26]. Projections yield contructible algebraic sets whose closure in both the Euclidean and Zariksi topologies are equal. The following, Lemma 3 from [27], is essential for computing corank c without first computing c.

Proposition 5. Let F : C N → C n be a polynomial system, V ⊂ { x ∈ C N | F ( x ) = 0 } ⊂ C N be irreducible of multiplicity 1 with respect to F, and π(x₁, …, x_N) = (x₁, …, x_a) for some a ≤ N. For general z ∈ V,

Example 6. With the setup from Ex. 3, write the function x(t), u(t), and y(t) using a Taylor series expansion centered at t = 0, namely

Since (3) holds for all t, one obtains equations by substituting these Taylor series expansions into (3) and taking coefficients with respect to t. For r ≥ 0, let F_r be the system obtained by taking coefficients of 1, t, t², …, t^r. For this linear compartment model, the coefficients of t^j are

Based on the structure of each G_j, it is clear that the Jacobian matrix of F_r has full rank, namely 4(r + 1), at every point. In fact, F_r = 0 defines an irreducible and smooth solution set of codimension 4(r + 1) (dim = 11 + r). We can compute a random point on this solution set by randomly selecting the following 11 + r values: p, x₀, and u₀, …, u_r, and trivially computing the unique x_j+1 and y_j sequentially for j = 0, …, r via G_j = 0.

Next, one treats the coefficients of the input u(t) and output y(t) as constants in F_r. Thus, we have that F_r depends upon N_r = 13 + 3r variables and apply Prop. 5 to compute

This limit is obtained at a finite value of r in accordance with the Cartan-Kuranishi Theorem and can be observed by checking for stabilization between the values obtained from r to r + 1 as demonstrated in Table 1. We see that d₇ = d₈ = 2 = corank c and provide the extra rows to show how the entries stabilize. In particular, this confirms that (3) is unidentifiable with 2 dimensions of unidentifiability.

We summarize this computation of the dimension of unidentifiability without first explicitly computing the input-output equations c in the following pseudocode.

Method 2: Computing dimension of unidentifiability without input-output equations

Input: For each r ≥ 0, system F_r(q, x, u, y) consisting of the coefficients of 1, t, t², …, t^r and general point z_r such that F_r(z_r) = 0 where q consists of m₁ parameters.

Output: Dimension of unidentifibility ℓ = corank c = dim c⁻¹(c(p)) for general p ∈ C m 1.

For r = 0, 1, 2, …

 Compute d r = corank 0 J F r ( z r ) - corank m 1 J F r ( z r )

 If either d_r = 0 or r > 0 with corank₀JF_r(z_r) = corank₀JF_r−1(z_r−1) and corank m 1 J F r ( z r ) = corank m 1 J F r - 1 ( z r - 1 ) , return d_r.

Such an approach naturally extends to problems when the parameters and initial conditions are restricted to an irreducible component by simply appending to F_r the requested constraints and taking the test points to be general on the corresponding irreducible component. The following demonstrates this on Example 1 from [28].

Example 7. Consider the following three-compartment model [29]:

Let F r ′ be the system obtained by adding the constraint x₃(0) = 0 to F_r. Table 3 shows that the model (5) is now unidentifiable with one dimension of unidentifiability.

Identifiable functions

When a model (1) is unidentifiable, one can ask for functions of the parameters p which are actually functions of the coefficients c(p) of the input-output equations. Such functions are called identifiable functions. For example, every element of c is itself an identifiable function. This is algebraically formalized in the following.

Definition 8. Let c be as above. A real-valued function f(p) is identifiable if the field extension R ( f , c ) / R ( c ) is an algebraic field extension.

One can also consider the global and local identifiability of functions.

Definition 9. Let c be as above and f be an identifiable function. The function f is called globally identifiable from c if there exists a function ϕ such that ϕ ○ c = f. The function f is called locally identifiable from c if there exists a multi-valued function ξ such that for every p, f(p) is equal to an entry of the multi-valued function ξ ○ c(p).

Example 10. With the setup from Ex. 3, the function f(p) = k₀₁ + k₂₁ is globally identifiable with f = ϕ ○ c where ϕ(x₁, …, x₅) = x₄ − x₁, i.e., f = c₄ − c₁. The function g(p) = k₀₂ + k₁₂ + k₃₂ is locally identifiable with g² + c₄ g + c₅ = 0, i.e., g = ξ ○ c where That is, for

The entry of this 2-valued function which is equal to g(p) is selected based on the sign of

When a model is unidentifiable with ℓ = corank c dimensions of unidentifiability, the goal is to compute d = rank c algebraically independent identifiable functions. The problem of finding d “nice” algebraically independent identifiable functions can be described in the following way, where “nice” could be taken to mean low degree, sparse, or are easy to interpret in terms of the model, depending on the context.

Problem 11. For rational functions c with d = rank c, compute a “nice” transcendence basis of the field extension R ( c ) / R.

One way to locate identifiable functions is by computing Gröbner bases with respect to various elimination orderings of the ideal 〈c(q) − c(p)〉. This approach is described in [2, 30] and has been implemented in the web application COMBOS [30]. In addition to requiring c, e.g., computed using differential elimination techniques, the biggest disadvantage of this method is that Gröbner basis computations can be computationally expensive. Thus, COMBOS can fail even for relatively simple models. Alternatively, the program DAISY [1, 24] can sometimes be used to find identifiable functions. Specifically, the DAISY program gives the solution to the system of equations c(q) = c(p) for a randomly chosen numerical point p. Sometimes one can algebraically manipulate the solution to obtain functions of the form f(q) = f(p), but there are many cases where this cannot be done [2, 30]. Nonetheless, if one is able to obtain such f, the following shows that they are indeed identifiable functions.

Proposition 12. If f(q) − f(p) is an element of the ideal I = 〈 c ( q ) - c ( p ) 〉 ⊂ R [ p , q ], then f is an identifiable function. That is, if f is constant on irreducible components of generic fibers of c, then f is an identifiable function.

proof. If f(q) − f(p) is contained in I, then the dimension of the image of the combined map (c, f) is equal to the dimension of the image of the map c. In other words, the field extension R ( f , c ) / R ( c ) is an algebraic field extension showing that f is identifiable.

Reparametrization and other uses of identifiable functions

If one can solve Problem 11, one then tries to use the new basis to reparametrize the model. In [23], a method to find identifiable scaling reparametrizations is given for a certain class of linear compartment models where the identifiable functions are monomials. Currently, there is no general approach to find identifiable reparametrizations and, for most models, the reparametrizations are found using ad hoc approaches which work on a case-by-case basis.

Even if a reparametrization cannot be found, identifiable functions have other important uses. From the identifiable functions, one can determine which parameters need to be known in order to render the entire model identifiable. This information can also be determined from the solution of the system of equations c(q) = c(p). However, identifiable functions give us additional information if only a subset of those parameters can be determined. In other words, we can obtain a simpler set of identifiable functions of parameters if a subset of the parameters are known and, perhaps, for this new set of identifiable functions, computing an identifiable reparametrization is possible. This is the case for Ex. 23 below where knowledge of either the pair (a₃₄, a₄₃) or the pair (a₃₃, a₄₄) renders all the identifiable functions to be monomials, in which case the method in [23] can be used to find an identifiable scaling reparametrization.

Computing identifiability degree

For a model (1) that is identifiable, Problem 4 can be solved by computing the identifiability degree k ∈ N in order to distinguish between globally identifiable (k = 1) and locally identifiable (k > 1) models. k is simply the number of solutions of c(q) = c(p) for general p, where c is the collection of coefficients of the input-output equations. As mentioned above, the software package DAISY [1, 24] uses such an approach with Gröbner basis methods to count the number of solutions. One could also use numerical homotopy methods, e.g., as summarized in [6, 7], to compute k, as illustrated in the following example.

Example 13 As shown in Ex. 3, the model (3) has 2 dimensions of unidentifiability. With the aim of constructing an identifiable model, we modify (3) by adding the extra constraints k₀₁ = k₀₃ = 0 yielding a new model with only one leak parameter k₀₂. The coefficients c of the input-output equation for this simplified model are

We summarize this most basic approach for computing the identifiability degree when the input-output equations are known in the following pseudocode.

Method 3: Computing identifiability degree from input-output equations (direct solving)

Input: m₂ input-output equation coefficients c(q), depending on parameters q ∈ C m 1 for which corank c = 0, i.e., corresponding model is identifiable.

Output: Identifiability degree k ∈ N.

Choose random, complex values p of parameters q.

Use homotopy continuation to compute Z = { q ∈ C m 1 | c ( q ) = c ( p ) } .

Return k = #Z.

Rather than using a direct global solving method which is based on knowing the coefficients c, we next consider an alternative approach based on monodromy computations in numerical algebraic geometry that also can be used without computing c. We first describe the approach when c is known and then extend to the case when c is not explicitly computed.

Identifiability degree from input-output equations

Suppose that (1) is identifiable with identifiability degree k ∈ N and c is the set of coefficients of the input-output equations. Following the notation before Definition 1, let m₁ be the number of independent parameters p and m₂ be the number of entries in c so that c : C m 1 → C m 2. Assume that the model is identifiable so that corank c = 0 and rank c = dim X = m₁ where X = c ( C m 1 ) ¯ ⊂ C m 2. The continuity of c yields that X is irreducible. The graph of c, namely

One can compute k via a pseudowitness point set [27] for X. To that end, let L 2 ⊂ C m 2 be a general linear space of codimension m₁. The finite set W = Graph ( c ) ∩ ( C m 1 × L 2 ) is a pseudowitness point set for X with respect to Graph(c) and π where #W = k ⋅ deg X and #π(W) = deg X, i.e., k = #W/#π(W). In order to compute W, we follow the approach in [32] using monodromy loops [33], as follows.

We first note that it is trivial to construct one point w ∈ W as follows. One first selects a general point (p, c(p)) ∈ Graph(c) and then constructs a general linear space L 2 ⊂ C m 2 of codimension m₁ that passes through c(p). Hence, w = (p, c(p)) ∈ W.

Next, the irreducibility of Graph(c) ensures that pairs of points in W are connected via smooth paths on Graph(c). We aim to discover such connecting paths using random monodromy loops. For t ∈ [0, 1], let L ( t ) be a smooth path consisting of general linear spaces of codimension m₁ in C m 2 such that L ( 0 ) = L ( 1 ) = L 2. Hence, this defines paths w(t) defined by Graph ( c ) ∩ ( C m 1 × L ( t ) ) where w(1) ∈ W is known. Homotopy continuation computes the endpoint w(0), which must also be a point in W. If w(0) ≠ w(1), the resulting loop has produced a nontrivial monodromy action and potentially yielded a previously unknown point in W.

Example 14. For c : C 5 → C 5 in (6), we know that X = c ( C 5 ) ¯ = C 5, i.e., deg X = 1. Hence, we have that the identifiability degree k = #W where W is a pseudowitness point set for X.

For illustrative purposes, consider p = (−1, −2, 5, −1, −3) with c(p) = (−2, −31, 5, −1, −30) so L 2 = { ( - 2 , - 31 , 5 , - 1 , - 30 ) } has codimension 5 with c ( p ) ∈ L 2. Consider the loop

Running finitely many random monodromy loops necessarily yields a set W ^ ⊂ W that may fail to achieve the goal of equality. However, information about the model can be obtained even if W ^ ⊊ W. For example, if (p₁, c(p₁)) and (p₂, c(p₂)) are known points in W with c(p₁) = c(p₂) and p₁ ≠ p₂, then one knows the identifiability degree is larger than 1, i.e., the model is locally identifiable. A heuristic stopping criterion for when W ^ = W provided in [32] is simply to have many different random monodromy loops yielding no new points.

We use trace tests [34, 35] to provide a stopping criterion to recognize when W ^ = W. These are described and illustrated well in [36]. To make these monodromy and trace test computations more efficient, see [37, 38].

Example 15. To show that Ex. 14 computed both points in W, i.e., the degree of identifiability k = 2, for illustrative purposes, we consider the following three linear spaces in C 5:

The irreducible curve C = Graph ( c ) ∩ ( C 5 × M 2 ) has multidegree (5, 2), which is verified using the multihomogeneous trace test applied to C ∩ H ( t ). This yields k = 2.

We summarize this computation in the following pseudocode.

Method 4: Computing identifiability degree from input-output equations (monodromy)

Input: m₂ input-output equation coefficients c(q), depending on parameters q ∈ C m 1 for which corank c = 0, i.e., model is identifiable, and an integer maxUselessLoops.

Output: Identifiability degree k ∈ N or error along with a lower bound on k if the number of loops in a row that do not yield any new points is more than maxUselessLoops.

Choose random, complex values p of parameters q and compute c(p).

Form w = (p, c(p)) and W = {w}.

Construct general linear space L 2 ⊂ C m 2 of codimension m₁ that passes through c(p).

Set numUselessLoops = 0.

While numUselessLoops < maxUselessLoops

 Increment numUselessLoops = numUselessLoops + 1.

 Construct a general loop of linear spaces L ( t ) such that L ( 0 ) = L ( 1 ) = L 2.

 For each w ∈ W

  Use homotopy continuation applied to the homotopy Graph ( c ) ∩ ( C m 1 × L ( t ) ) to

  track from w at t = 1 to t = 0 yielding w′.

 If w′ ∉ W

  Update w = {W, w′} and numUselessLoops = 0.

  If trace test passes, return k = #W.

Return error with k = #W.

The advantage to using such a monodromy approach is that the structure of c may be such that k is small but this structure is not known a priori meaning that a homotopy for solving c(q) = c(p) requires tracking many homotopy paths. The disadvantage is that many monodromy loops may be needed to find all points necessary for the trace test to validate completeness when k is large.

Identifiability degree without input-output equations

In the previous section, we computed the degree of a general fiber of a generically finite-to-one coefficient map. This is based on the fact that one has the same input-output equation if and only if the coefficients agree. However, when we are using a truncated system as described in Example 6, namely F_r which depends upon the parameters p, input U = {u₀, …, u_r}, output Y = {y₀, …, y_r}, and state variables X = {x₀, …, x_r+1}, it provides necessary conditions to have the same input-output as shown in the following example.

Example 16. The following model is a modification of an HIV model from [39]:

As in the previous section, Table 4 shows that the system F₇ provides the model (8) is identifiable.

For example, consider the sufficiently general truncated output

We know that there are finitely many values of the parameters p which yield this output. Monodromy yields the following 12 values of the parameters (listed to four decimal places):

This table shows that there are 3 distinct values of y₈, each of which is obtained by 4 values of the parameters indicating that the identifiability degree is 4.

This example shows that even though F_r is enough to show identifiability, we may only need to consider a subset of the corresponding parameter values which have the same input-output.

The structure of (1) clearly shows that the solution set of F_r = 0 is irreducible, smooth, and parameterized by p, U, and x₀. Thus, it is trivial to construct a generic point (p*, X*, U*, Y*) in the solution set of F_r = 0. From this point, we can use Prop. 5 to compute the dimension d ≥ 0 of the solution set of F_r(p*, X, U*, Y*) = 0, i.e., the dimension of the state variables. If d > 0, we can add d general linear slices in X to F_r to reduce to the case when d = 0.

With this reduction, we repeatedly apply random monodromy loops to compute all values of p such that there exists X with

By testing the finitely many values of the parameters p, the identifiability degree k is the number of points corresponding to the same input-output. To verify the completeness, we simply apply the multihomogeneous trace test via the parameter space and the input-output space.

To save space, we exclude pseudocode for this method as it is so similar to Method 1. The primary change is that the set of coefficients c is replaced by the truncated system F_r for some value of r along with an extra computation to test for the same input-output values.

Example 17. Reconsidering the model (8) in Ex. 16 which has no input, we first restrict the output space to, for illustration purposes, the sufficiently general line

Thus, we apply the multihomogeneous trace test by solving F₇ = 0 on this line intersected with the sufficiently general family of bilinear hyperplanes in the parameter and output space:

Monodromy followed by the trace test confirms that the bidegree is (222, 12). Hence, the number of elements in Table 5 is complete.

We can simplify this computation, for example, by instead taking the following family

The bidegree in p₅ and the output space is (60, 12) which again shows that Table 5 is complete.

Computing identifiable functions

A model (1) is identifiable if and only if every function of the parameters is an identifiable function. In particular, each irreducible component of a generic fiber of the coefficients c of the input-output equations is a singleton for an identifiable model. Since every function of the parameters is trivially constant on each singleton, Prop. 2 yields that every function is identifiable. To be a globally identifiable function, it must take the same constant value on all of the irreducible components of a general fiber.

Example 18. With the setup from Ex. 16, the model (8) is identifiable with identifiability degree 4. Hence, for example, we know that f₁ = p₄ and f₂ = p₅ are both identifiable functions. From the first two rows of Table 5, we see that both f₁ and f₂ are not globally identifiable since each of them take two different values. The functions g₁ = p₂, g₂ = p₃, and g₃ = p₄ + p₅ are all globally identifiable since each of them take the same value at all four points.

To compute identifiable functions, we will first use numerical algebraic geometry to sample points from fibers. Then, given a finite collection of terms, we will use exactness recovery techniques, e.g., [40], or interpolation to construct identifiable functions from the sample data. Computing globally identifiable functions simply requires computing points on all irreducible components and adding additional constraints.

Sampling

In the case that input-output equations have been computed, let c be the collection of coefficients of the input-output equations and suppose that d ≥ 0 is the dimension of unidentifiability. Thus, for a given generic point p, the point q = p is a smooth point on an irreducible component V_p of dimension d of the solution set defined by c(q) − c(p) = 0. Hence, when d > 0, we can sample other points in this irreducible component as follows. Let L p be a general linear space of codimension d passing through p and L be some other general linear space of codimension d. By using homotopy continuation, we can track the solution path q(t) defined by q(1) = p and

This yields the point q(0) which is also a generic point in V_p.

One can easily compute other points in this same fiber V_p by repeating with a different linear space L and sample other fibers by repeating the process with a different generic point p.

With the aim of computing globally identifiable functions, sample points in every irreducible component of c(q) − c(p) = 0 are needed. In this case, one simply constructs an identifiable system by restricting the parameters to a general linear space of codimension d and applying the techniques above to the resulting system. That is, if p ∈ C m 1, we take a general affine linear mapping b : C m 1 - d → C m 1 so that c ( b ( q ^ ) ) - c ( b ( p ^ ) ) = 0 has finitely many solutions for generic p ^, say q 1 = b ( q ^ 1 ) , … , q k = b ( q ^ k ), i.e., the model with parameters p = b ( p ^ ) is identifiable over C m 1 - d with identifiability degree k. Applying the slice moving described above, one can sample points in all components of the fiber over p using the points q₁, …, q_k.

Example 19. Reconsider (3) in Ex. 3 for which c shows the model has d = 2 dimensions of unidentifiability. For illustration, with p = (1, 2, 3, 4, 5, 6, 7), we can take L p to be

Hence, since all of the values of the parameters changed, we know that each parameter itself is an unindentifiable function.

If, for illustration, we take the affine linear mapping b : C 5 → C 7 defined by

Thus, we have computed at least one point in each irreducible component of the fiber over p.

Without input-output equations, one simply uses a truncated system F_r as described in Example 6 to perform the same computations. The only potential issues were addressed above, namely reduction to the case that the state variables are generically zero-dimensional over the parameter-input-output space and restricting to the irreducible components which have the same input-output. The latter is accomplished by simply ignoring the components which have different input-output values.

Example 20. To illustrate moving on an irreducible component, we describe the setup to yield the same corresponding endpoint in (12). To that end, following Ex. 6, we utilize F₇. Starting with parameter values p = (1, 2, 3, 4, 5, 6, 7), the structure of F₇ makes it trivial to generate general input, output, and state variables satisfying F₇ = 0, i.e., randomly selecting input U and initial conditions x₀ for the state variables trivially yields the values of the other state variables x₁, …, x₈ and output Y. Then, by holding the input U and output Y fixed, we track along the solution path where the variables consist of the model parameters and the state variables defined by F₇ = 0 that deforms L p in (10) to L in (11). The resulting endpoint corresponds with the endpoint in (12).

Functions from samples

From the ability to sample points described in the previous section, we can reconstruct identifiable functions in a given finite-dimensional vector space of functions, say F = span { f 1 , … , f j }. Following Prop. 2, an identifiable function f ∈ F is constant on irreducible components of generic fibers of c, which corresponds with computing null spaces of linear equations described as follows.

We can express every f ∈ F as f = ∑ i = 1 j a i f i where a = ( a 1 , … , a j ) ∈ C j. If p is a generic value of the parameters, using the sampling method above, we can compute a generic q_p in the same irreducible component V_p. Hence, the condition f(q_p) = f(p) imposes a linear constraint on a, namely

One option is to keep imposing more such conditions by selecting other general values of p with corresponding q_p. The dimension of the null space reduces by one monotonically with each new condition until it reaches the dimension of the linear span of the identifiable functions in F.

Alternatively, for computing identifiable functions with integer coefficients, i.e., a ∈ Z j, one general point is enough via exactness recovery methods [40].

Example 21. Let F = span { k 01 , k 02 , k 03 , k 12 , k 13 , k 21 , k 32 }, p = (1, 2, 3, 4, 5, 6, 7), and q_p as in (12). Then, integer solutions to (q_p − p) ⋅ a = 0 computed using [40] correspond to:

Alternatively, one can sample V_p for five general values of p and observe that the first four impose a new linear constraint on the coefficients a while the fifth one is redundant. This shows that there is a three-dimensional linear space of identifiable functions in F spanned by the three linear functions above.

We bring all methods of this section together in the following brief high-level pseudocode.

Method 5: Computing identifiable functions via sampling

Input: Input-output equation coefficients c(q), depending on parameters q ∈ C m 1 (if available), else the truncated system F_r for some r and a basis f₁, …, f_j for a linear space of polynomials F of interest.

Output: Identifiable functions in F.

Choose random, complex values p of parameters q.

Compute a point on each irreducible component of c⁻¹(c(p)) using either c or F_r.

Use homotopy sampling to find additional points on each irreducible component.

Use the sample points together with exactness recovery methods to find all identifiable

 functions in F.

Return all discovered identifiable functions.

Globally identifiable functions are computed by simply adding the condition that the function takes the same constant value on all irreducible components of general fibers which are sampled using the methods described above. Since globally identifiable functions are a subset of the identifiable functions, one need only search inside of the space of identifiable functions in F.

Example 22. From the seven points in (13) corresponding with p = (1, 2, 3, 4, 5, 6, 7), we see that k₀₁ + k₂₁ is globally identifiable (always taking the value 7 on these eight points) whereas k₀₃ + k₁₃ and k₀₂ + k₁₂ + k₃₂ are not globally identifiable. However, from the sample points, it is easy to see that their sum, namely k₀₂ + k₀₃ + k₁₂ + k₁₃ + k₃₂, is globally identifiable.

The selection of the test space F is a user-defined input and is based on the structure of the identifiable functions of interest, e.g., linear functions, polynomials of low degree, or linear span of rational monomials where the numerator and denominator have low degree.

Results

We now demonstrate our methods on two larger examples. Throughout the paper, for illustrative purposes, the examples presented typically selected small integer values for random numbers. In practice, including the following examples, we select random complex numbers. Data for computations available at http://dx.doi.org/10.7274/R03T9F91.

Example 23. The following is a 4-compartment model from Example 6.3 of [23]:

This model, which has parameters p = (a₁₁, a₁₂, a₂₁, a₂₂, a₂₃, a₃₃, a₃₄, a₄₂, a₄₃, a₄₄), input u(t), state variables x₁(t), x₂(t), x₃(t), x₄(t), and output y(t), does not fit the criteria presented in [23] for computing identifiable functions. Nonetheless, the method provided in [17, 23] is able to compute the input-output equations where the set c : C 10 → C 7 of coefficients is

Using Prop 2, the model is unidentifiable with 4 dimensions of unidentifiability. Therefore, to solve Problem 11, we need to compute 6 algebraically independent identifiable functions.

We utilize the sampling and interpolation methods above to sample and construct the identifiable functions. For example, sampling yields two values of the parameters, provided in Table 6 rounded to four decimal places, so that every identifiable function must take the same value on both. In particular, we immediately see that both f₁ = a₁₁ and f₂ = a₂₂ are identifiable. Applying interpolation as above to the space of linear forms also yields the identifiable linear function f₃ = a₃₃ + a₄₄.

Considering the space of polynomials of degree at most 2 which are algebraically independent of f₁, f₂, f₃ yields f₄ = a₁₂a₂₁ and f₅ = a₃₃a₄₄−a₃₄a₄₃.

Finally, the space of polynomials of degree at most 3 which are algebraically independent of f₁, …, f₅ yields f₆ = a₂₃a₃₄a₄₂.

To show that f₁, …, f₆ are actually globally identifiable, we use the approach above to sample points from every irreducible component. The result of this process is that a generic fiber only has one irreducible component thereby showing global identifiability. We could also have used Defn. 9 to show global identifiability. This is demonstrated by the following:

Example 24. The following is a model from biochemical reaction network theory for the mitogen-activated protein kinase (MAPK) pathway [41] which is part of a molecular signaling network that governs the growth, proliferation, differentiation, and survival of many cell types:

This model has 12 state variables

We will consider several different cases of what is measured as output. In all of our examples, we attempted to first compute input-output equations using differential elimination via the command RosenfeldGroebner in Maple [42]. In all of our attempts, the differential elimination failed to terminate meaning that we will just utilize the model equations in the following.

First, for taking the standard 6 measurable outputs:

Table 7, computed in about a minute on a single processor, shows that the resulting model is identifiable.

For comparison of methods, neither DAISY [1, 24] nor COMBOS [30] finished the identifiability computations for this model after running for 24 hours. To the best of our knowledge, this is the first successful implementation of a structural identifiability test for this model.

Second, if we adjust the model so that we only take the following 2 measurable outputs:

Table 8 shows that the resulting model is still identifiable.

Third, if we take the following 4 measurable outputs:

Table 9 shows that the resulting model is still identifiable.

Finally, we consider 10 new mixing parameters, namely

Table 10 shows that the resulting model, which has a total of 32 parameters, is unidentifiable with one dimension of unidentifiability.

Using the results above, we can observe from sampling that each irreducible component of a general fiber is simply a line and the following 16 parameters are all identifiable:

Conclusion

In this article, we considered the problems of determining the identifiability of an ODE model, computing the identifiability degree in the case that the model is identifiable and identifiable functions in the case that the model is unidentifiable. To summarize, the results of this article include numerical methods for the following:

• compute the dimension of unidentifiability with or without input-output equations;

• for identifiable models, compute the identifiability degree with or without input-output equations using basic homotopy continuation or monodromy loops;

• for unidentifiable models, compute identifiable and globally identifiable functions inside of a linear family of functions with or without input-output equations.

These methods were illustrated on several examples, including the first known structural identifiability result for MAPK in Example 24.

In the future, we hope to apply similar numerical algebraic geometry methods to other areas in biological modelling, such as controllability, observability, and indistinguishability.