Identifiability and numerical algebraic geometry.

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.

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: where f and g are vectors of rational functions, x(t) is the state variable vector, p is the parameter vector which is assumed to be constant, u(t) is the input vector, and y(t) is the output vector. In the following, only the input u(t) and output y(t) vectors are assumed to be known, i.e., the state variables x(t) and the parameters p are unknown.

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 jth one can be written as where each ψ(u, y) is a differential monomial. Each c(p) is a rational function in the parameters p, forming a collection c(p) called the coefficients of the input-output equations. The coefficients of each input-output equation can be determined uniquely by normalizing each input-output equation so that one of the coefficients is one.

Deciding identifiability

Let m1 denote the number of independent parameters p and m2 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 . 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 , let

ℓ = dim X ≥ 0, and . 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., , and unidentifiable if ℓ > 0, i.e., k = ∞.

When identifiable, the number 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, . The corank of c is corank c = m1 − rank c. The following, which is Lemma 13.4.1 of [7] (see also [17]), relates ℓ and corank c.

Proposition 2. For a general , ℓ = dim X (2) is equal to corank c where (1) is identifiable if and only if corank .

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: m2 input-output equation coefficients c(p), depending on parameters .

Output: Dimension of unidentifibility ℓ = corank c = dim c−1(c(q)) for general .

Choose random, complex values .

Return ℓ = corank J(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) = (x1(t), x2(t), x3(t)), and unknown parameters p = (k01, k02, k03, k12, k13, k21, k32), where k represents the rate of transfer from compartment j to compartment i and k0 represents a leak from compartment i to outside the system:

Fig 1 presents a pictorial representation of this model.

Fig 1

A 3-compartment model.

A 3-compartment model with input (represented by the arrowhead) and output (represented by the line segment with a circle at the end) in the first compartment and “leaks” from every compartment (represented by arrows leaving the compartments). Here, the input could represent a drug injection and the first compartment could represent blood, with the other two compartments representing organs, e.g., tissue and stomach. The unknown parameters represent rates of transfer from one compartment to another (drawn as arrows in the figure), or in the case of leaks, from one compartment to outside the system. The state variables represent drug concentration in the blood and organs, with output from the first compartment representing measured drug concentration in the blood.

The approach described in [17, 23] yields the input-output equation: where such that E(z1, …, z) is the kth elementary symmetric polynomial in z1, …, z. Thus, for c = (c1, …, c5), it is easy to verify that rank c = 5 and c = 2 so that this model is unidentifiable with 2 dimensions of unidentifiability.

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 . 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 be a polynomial system, be irreducible of multiplicity 1 with respect to F, and π(x1, …, x) = (x1, …, x) for some a ≤ N. For general z ∈ V, where JF(z) is the Jacobian matrix of F evaluated at z and corank M is the corank of the last N − j columns of M.

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 be the system obtained by taking coefficients of 1, t, t2, …, t. For this linear compartment model, the coefficients of t are

Based on the structure of each G, it is clear that the Jacobian matrix of F has full rank, namely 4(r + 1), at every point. In fact, F = 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, x0, and u0, …, u, and trivially computing the unique x and y sequentially for j = 0, …, r via G = 0.

Next, one treats the coefficients of the input u(t) and output y(t) as constants in F. Thus, we have that F depends upon N = 13 + 3r variables and apply Prop. 5 to compute since . We trivially know d ≥ d since F is a subset of F. Hence, is a sequence of nonincreasing nonnegative integers that stabilizes with

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 d7 = d8 = 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.

Table 1

Summary of computations showing 2 dimensions of unidentifiability.

r	Nr	corank0 JFr	corank7 JFr	dr
0	13	9	2	7
1	16	8	1	7
2	19	7	0	7
3	22	6	0	6
4	25	5	0	5
5	28	4	0	4
6	31	3	0	3
7	34	2	0	2
8	37	2	0	2
9	40	2	0	2
10	43	2	0	2

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(q, x, u, y) consisting of the coefficients of 1, t, t2, …, t and general point z such that F(z) = 0 where q consists of m1 parameters.

Output: Dimension of unidentifibility ℓ = corank c = dim c−1(c(p)) for general .

For r = 0, 1, 2, …

Compute

If either d = 0 or r > 0 with corank0JF(z) = corank0JF(z) and return d.

Such an approach naturally extends to problems when the parameters and initial conditions are restricted to an irreducible component by simply appending to F 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]: with state variables x(t) = (x1(t), x2(t), x3(t)), input u(t), output y(t), and unknown parameters p = (p12, p13, p21). Using a similar setup from Ex. 6 summarized in Method 2, Table 2 shows that the model (5) is identifiable.

Table 2

Summary of computations showing the model (5) is identifiable.

r	Nr	corank0 JFr	corank3 JFr	dr
0	9	5	2	3
1	12	4	1	3
2	15	3	0	3
3	18	2	0	2
4	21	1	0	1
5	24	0	0	0

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

Table 3

Summary of computations showing that (5) is unidentifiable when x3(0) = 0.

r	Nr	corank_0JFr′	corank_3JFr′	dr
0	9	4	1	3
1	12	3	0	3
2	15	2	0	2
3	18	1	0	1
4	21	1	0	1

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 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) = k01 + k21 is globally identifiable with f = ϕ ○ c where ϕ(x1, …, x5) = x4 − x1, i.e., f = c4 − c1. The function g(p) = k02 + k12 + k32 is locally identifiable with g2 + c4 g + c5 = 0, i.e., g = ξ ○ c where That is, for we have

The entry of this 2-valued function which is equal to g(p) is selected based on the sign of i.e., the “+” entry when k02 + k12 + k32 − k03 − k13 ≥ 0 and the “–” entry otherwise.

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 .

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 , 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 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 (a34, a43) or the pair (a33, a44) 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 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 k01 = k03 = 0 yielding a new model with only one leak parameter k02. The coefficients c of the input-output equation for this simplified model are which is easily seen to have rank 5, i.e., the model is identifiable. For general , the system consists of 5 equations (1 cubic, 2 quadratic, and 2 linear) in 5 variables. Using a total degree homotopy (see [6, 7] for more details), one tracks 3 ⋅ 22 ⋅ 12 = 12, i.e., the total degree of (7), solution paths. Tracking these paths with homotopy continuation, e.g., via Bertini [31], yields 2 solutions to (7). One can also use a Gröbner basis computation to see that (7) has 2 solutions. These computations show that the model (3) with k01 = k03 = 0 is locally identifiable with identifiability degree of 2.

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: m2 input-output equation coefficients c(q), depending on parameters for which corank c = 0, i.e., corresponding model is identifiable.

Output: Identifiability degree .

Choose random, complex values p of parameters q.

Use homotopy continuation to compute

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 and c is the set of coefficients of the input-output equations. Following the notation before Definition 1, let m1 be the number of independent parameters p and m2 be the number of entries in c so that . Assume that the model is identifiable so that corank c = 0 and rank c = dim X = m1 where . The continuity of c yields that X is irreducible. The graph of c, namely is also irreducible of dimension m1. In fact, for the projection map , we know that and π restricted to Graph(c) is generically a k-to-1 map.

One can compute k via a pseudowitness point set [27] for X. To that end, let be a general linear space of codimension m1. The finite set 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 of codimension m1 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 be a smooth path consisting of general linear spaces of codimension m1 in such that . Hence, this defines paths w(t) defined by 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 in (6), we know that , 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 has codimension 5 with . Consider the loop where s(t) = 1 − e2 and . Hence, is a loop with . For the path with w(1) = (p, c(p)), we have w(0) = (q, c(q)) where q = (5/6, −2, −6, −1, 37/6) and c(q) = c(p) showing {w(0), w(1)} ⊂ W and k = #W ≥ 2.

Running finitely many random monodromy loops necessarily yields a set that may fail to achieve the goal of equality. However, information about the model can be obtained even if . For example, if (p1, c(p1)) and (p2, c(p2)) are known points in W with c(p1) = c(p2) and p1 ≠ p2, then one knows the identifiability degree is larger than 1, i.e., the model is locally identifiable. A heuristic stopping criterion for when 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 . 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 : with . We take

The irreducible curve has multidegree (5, 2), which is verified using the multihomogeneous trace test applied to . This yields k = 2.

We summarize this computation in the following pseudocode.

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

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

Output: Identifiability degree 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 of codimension m1 that passes through c(p).

Set numUselessLoops = 0.

While numUselessLoops < maxUselessLoops

Increment numUselessLoops = numUselessLoops + 1.

Construct a general loop of linear spaces such that .

For each w ∈ W

Use homotopy continuation applied to the homotopy 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 which depends upon the parameters p, input U = {u0, …, u}, output Y = {y0, …, y}, and state variables X = {x0, …, x}, 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 F7 provides the model (8) is identifiable.

Table 4

Summary of computations showing (8) is identifiable.

r	corank0 JFr	corank0 JFr	dr
0	7	2	5
1	6	1	5
2	5	0	5
3	4	0	4
4	3	0	3
5	2	0	2
6	1	0	1
7	0	0	0

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 y8, 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 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 = 0 is irreducible, smooth, and parameterized by p, U, and x0. Thus, it is trivial to construct a generic point (p*, X*, U*, Y*) in the solution set of F = 0. From this point, we can use Prop. 5 to compute the dimension d ≥ 0 of the solution set of F(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 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 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 F7 = 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.

Table 5

12 possible values of parameters of system F7 from model (8).

p1	p2	p3	p4	p5	y8
±0.1253	−2.4825	4.4249	−0.9210	−0.2137	0.1706
±0.2602	−2.4825	4.4249	−0.2137	−0.9210	0.1706
0.3023 ± 0.0779i	−3.5234 ± 0.5105i	4.2201 ± 1.9168i	−1.3367 ∓ 0.0298i	−0.1080 ∓ 0.2292i	0.1107 ∓ 0.4040i
−0.3023 ± 0.0779i	−3.5234 ∓ 0.5105i	4.2201 ∓ 1.9168i	−1.3367 ± 0.0298i	−0.1080 ± 0.2292i	0.1107 ± 0.4040i
0.6847 ± 0.2133i	−3.5234 ∓ 0.5105i	4.2201 ∓ 1.9168i	−0.1080 ± 0.2292i	−1.3367 ± 0.0298i	0.1107 ± 0.4040i
−0.6847 ± 0.2133i	−3.5234 ± 0.5105i	4.2201 ± 1.9168i	−0.1080 ∓ 0.2292i	−1.3367 ∓ 0.0298i	0.1107 ∓ 0.4040i

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

The bidegree in p5 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 f1 = p4 and f2 = p5 are both identifiable functions. From the first two rows of Table 5, we see that both f1 and f2 are not globally identifiable since each of them take two different values. The functions g1 = p2, g2 = p3, and g3 = p4 + p5 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 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 be a general linear space of codimension d passing through p and 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.

One can easily compute other points in this same fiber V by repeating with a different linear space 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 , we take a general affine linear mapping so that has finitely many solutions for generic , say , i.e., the model with parameters is identifiable over 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 q1, …, q.

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 to be and where . Tracking the path defined by (9) yields the endpoint (to four decimal places):

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 defined by the resulting model is identifiable with identifiability degree 8 and the following 7 other points corresponding with b(1, 2, 3, 4, 5) = (1, 2, 3, 4, 5, 6, 7):

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 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 F7. Starting with parameter values p = (1, 2, 3, 4, 5, 6, 7), the structure of F7 makes it trivial to generate general input, output, and state variables satisfying F7 = 0, i.e., randomly selecting input U and initial conditions x0 for the state variables trivially yields the values of the other state variables x1, …, x8 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 F7 = 0 that deforms in (10) to 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 . Following Prop. 2, an identifiable function 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 as where . If p is a generic value of the parameters, using the sampling method above, we can compute a generic q in the same irreducible component V. Hence, the condition f(q) = 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. 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 .

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

Example 21. Let , p = (1, 2, 3, 4, 5, 6, 7), and q as in (12). Then, integer solutions to (q − p) ⋅ a = 0 computed using [40] correspond to:

Alternatively, one can sample V 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 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 (if available), else the truncated system F for some r and a basis f1, …, f for a linear space of polynomials of interest.

Output: Identifiable functions in .

Choose random, complex values p of parameters q.

Compute a point on each irreducible component of c−1(c(p)) using either c or F.

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 .

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 .

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

The selection of the test space 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 = (a11, a12, a21, a22, a23, a33, a34, a42, a43, a44), input u(t), state variables x1(t), x2(t), x3(t), x4(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 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 f1 = a11 and f2 = a22 are identifiable. Applying interpolation as above to the space of linear forms also yields the identifiable linear function f3 = a33 + a44.

Table 6

Two values of the parameters rounded to four decimal places.

a11	−0.6690 − 0.1758i	−0.6690 − 0.1758i
a12	−0.1669 + 0.3165i	1.3705 − 0.4117i
a21	2.3433 + 0.6225i	−0.5219 + 0.3086i
a22	−0.6286 − 0.1868i	−0.6286 − 0.1868i
a23	0.4005 − 0.5144i	2.5585 + 0.5746i
a33	2.1248 − 0.6011i	0.2095 − 0.4521i
a34	1.1295 − 0.8604i	0.8611 + 0.5272i
a42	−0.4210 + 0.6785i	0.2734 − 0.0567i
a43	−1.1126 − 0.0416i	−0.1132 − 0.7724i
a44	−0.6880 + 0.3317i	1.2273 + 0.1827i

Considering the space of polynomials of degree at most 2 which are algebraically independent of f1, f2, f3 yields f4 = a12a21 and f5 = a33a44−a34a43.

Finally, the space of polynomials of degree at most 3 which are algebraically independent of f1, …, f5 yields f6 = a23a34a42.

To show that f1, …, f6 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 and 22 parameters

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.

Table 7

Summary of computations showing model is identifiable.

r	corank0 JFr	corank22 JFr	dr
0	28	6	22
1	23	1	22
2	18	0	18
3	13	0	13
4	8	0	8
5	3	0	3
6	0	0	0

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.

Table 8

Summary of computations showing model is identifiable.

r	corank0 JFr	corank22 JFr	dr
0	32	10	22
1	30	8	22
2	28	6	22
3	26	4	22
4	24	2	22
5	22	0	22
6	20	0	20
7	18	0	18
8	16	0	16
9	14	0	14
10	12	0	12
11	10	0	10
12	8	0	8
13	6	0	6
14	4	0	4
15	2	0	2
16	0	0	0

Third, if we take the following 4 measurable outputs:

Table 9 shows that the resulting model is still identifiable.

Table 9

Summary of computations showing model is identifiable.

r	corank0 JFr	corank22 JFr	dr
0	30	8	22
1	26	4	22
2	22	0	22
3	18	0	18
4	14	0	14
5	10	0	10
6	6	0	6
7	2	0	2
8	0	0	0

Finally, we consider 10 new mixing parameters, namely with the following 4 measurable outputs:

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

Table 10

Summary of computations showing 1 dimension of unidentifiability.

r	corank0 JFr	corank32 JFr	dr
0	40	8	32
1	36	4	32
2	32	0	32
3	28	0	28
4	24	0	24
5	20	0	20
6	16	0	16
7	12	0	12
8	8	0	8
9	4	0	4
10	1	0	1
11	1	0	1

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: meaning a00, a01, a10, α01, α10, α11 and the 10 mixing parameters are unidentifiable. In fact, no nonconstant linear function in these 16 latter unidentifiable parameters is 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.

20 Sep 2019

PONE-D-19-15150

Identifiability and numerical algebraic geometry

PLOS ONE

Dear Prof. Hauenstein,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

We would appreciate receiving your revised manuscript by Nov 03 2019 11:59PM. When you are ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter.

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). This letter should be uploaded as separate file and labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. This file should be uploaded as separate file and labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. This file should be uploaded as separate file and labeled 'Manuscript'.

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

We look forward to receiving your revised manuscript.

Kind regards,

Fang-Bao Tian

Academic Editor

PLOS ONE

Journal Requirements:

When submitting your revision, we need you to address these additional requirements.

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

http://www.journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and http://www.journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. We note that you have stated that you will provide repository information for your data at acceptance. Should your manuscript be accepted for publication, we will hold it until you provide the relevant accession numbers or DOIs necessary to access your data. If you wish to make changes to your Data Availability statement, please describe these changes in your cover letter and we will update your Data Availability statement to reflect the information you provide.

Additional Editor Comments (if provided):

Thank you for submitting your work to the Plos One. It has been reviewed by two specialists. As you can see, their opinions are quit different. I would encourage you to revise the paper based on their comments, especially those raised by Reviewer 1.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Partly

Reviewer #2: Yes

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: No

Reviewer #2: Yes

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: No

Reviewer #2: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: The paper is well-written and the questions on identifiability posed at the beginning are important. However,

all systems considered in examples have constant coefficients, hence presumably analytical solutions for the main dependent variable are available explicitly in terms of the unknown constant parameters. In reality, systems have non-homogeneous and/or nonlinearities. Can the approach proposed deal with such practical generalities of inhomogeneous and nonlinear systems?

There is some exaggeration of formalism in the first parts of the paper, which is actually not needed.

Can the analysis be extended to PDE's?

The references should be written in the style of the journal and be consistent throughout.

Several references are on arXiv, i.e. not reviewed yet, and limit the credibility of the citations.

The word 'geometry' used is not fully justified and it could be misleading. Suggest removing and find a better formulation.

Given a set of parameters that are desired to estimate and a governing set of equations, can you design the input-output data for identifiability? One can als consult the early works of Kitamura and Nakagiri on the identifiability of distributed systems.

Reviewer #2: I strongly recommend publication after the revisions suggested in the attached PDF. The techniques you provide establish identifiability results for models that would be too large for established methods. As such, it is a substantial contribution to the field of parameter identifiability analysis. Please see attached document for more details.

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files to be viewed.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org. Please note that Supporting Information files do not need this step.

Submitted filename: bates-review.pdf

Click here for additional data file.

31 Oct 2019

Thank you for your careful review. We have incorporated all of your comments into our revised manuscript as described in the cover letter.

Submitted filename: Response to Reviewers.pdf

Click here for additional data file.

25 Nov 2019

Identifiability and numerical algebraic geometry

PONE-D-19-15150R1

Dear Dr. Hauenstein,

We are pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it complies with all outstanding technical requirements.

Within one week, you will receive an e-mail containing information on the amendments required prior to publication. When all required modifications have been addressed, you will receive a formal acceptance letter and your manuscript will proceed to our production department and be scheduled for publication.

Shortly after the formal acceptance letter is sent, an invoice for payment will follow. To ensure an efficient production and billing process, please log into Editorial Manager at https://www.editorialmanager.com/pone/, click the "Update My Information" link at the top of the page, and update your user information. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to enable them to help maximize its impact. If they will be preparing press materials for this manuscript, you must inform our press team as soon as possible and no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

With kind regards,

Fang-Bao Tian

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #2: All comments have been addressed

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #2: Yes

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #2: Yes

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #2: Yes

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #2: Yes

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #2: (No Response)

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #2: No

5 Dec 2019

PONE-D-19-15150R1

Identifiability and numerical algebraic geometry

Dear Dr. Hauenstein:

I am pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please notify them about your upcoming paper at this point, to enable them to help maximize its impact. If they will be preparing press materials for this manuscript, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

For any other questions or concerns, please email plosone@plos.org.

Thank you for submitting your work to PLOS ONE.

With kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Fang-Bao Tian

Academic Editor

PLOS ONE