# A Parameter-Free Model Comparison Test Using Differential Algebra

Authors: Heather A. Harrington; Kenneth L. Ho; Nicolette Meshkat
Journal: Complexity
 (2019)

---

## Abstract

We present a method for rejecting competing models from noisy time-course data that does not rely on parameter inference. First we characterize ordinary differential equation models in only measurable variables using differential-algebra elimination. This procedure gives input-output equations, which serve as invariants for time series data. We develop a model comparison test using linear algebra and statistics to reject incorrect models from their invariants. This algorithm exploits the dynamic properties that are encoded in the structure of the model equations without recourse to parameter values, and, in this sense, the approach is parameter-free. We demonstrate this method by discriminating between different models from mathematical biology.


---

## Body

## 1. Introduction

Given competing mathematical models to describe a process, we wish to know whether our data are compatible with the candidate models. Often comparing models requires optimization and fitting time-course data to estimate parameter values and then apply an information criterion to select a “best” model [1]. However sometimes it is not feasible to estimate the value of these unknown parameters (e.g., large parameter space, nonlinear objective function, nonidentifiability, etc.). In this paper, we compare candidate models with time-course data while avoiding the parameter estimation problem by considering a “parameter-free” approach.The parameter problem has motivated the growth of fields that embrace a parameter-free flavor such as chemical reaction network theory and stoichiometric theory [2–4]. However many of these approaches are limited to comparing the behavior of models at steady-state [5–7]. Inspired by techniques commonly used in applied algebraic geometry [8] and algebraic statistics [9], methods for discriminating between possible models without estimating parameters have been developed for steady-state data [10, 11]. These approaches characterize a model in only observable variables—called asteady-state invariant [5]—using techniques from computational algebraic geometry and determine whether the noisy steady-state data are compatible with this steady-state invariant via a statistical test. However, unlike other Bayesian and parameter estimation approaches, it does not select models; it can only rule them out. Notably the method does not require parameter estimation, hence there is the term parameter-free.Extending the method developed in [10], we present a method for comparing models using time-course data instead of steady-state data. In this approach we compute input-output equations, which we refer to as input-output invariants for time series data. We consider state-space ordinary differential equations (ODE) models of the form x˙(t)=f(x(t),u(t),p) and y(t)=g(x(t),p) where xk(t) are species variables, k=1,…,N,ui(t) is a known input into the system, i=1,…,L, yj(t) is a known output (measurement) from the system, j=1,…,M, p is the unknown R-dimensional parameter vector, and the functions f,g are rational functions of their arguments. The dynamics of the model can be observed in terms of a time series where u(t) is the input at discrete points and y(t) is the output.In this setting, we aim to characterize our ODE models by eliminating variables that we cannot measure using differential elimination from differential algebra [12]. From the elimination, we obtain a system of equations in 0, 1, and higher order derivatives forming theinput-output invariants: Fj(u,u˙,u¨,u⃛,…,y,y˙,y¨,y⃛,…)=0, j=1,…,M. Importantly, the coefficients of these equations are rational functions of the parameters. We will see shortly that, in the linear case, Fj is a linear differential equation. For nonlinear models, Fj is nonlinear. Computing input-output invariants is described in Section 2.In order to test model compatibility, we substitute the data into the input-output invariant, which is given in the form ofu,u˙,u¨,u⃛,…,y,y˙,y¨,y⃛,… evaluated at given time points. This results in a linear system of equations, Aκ=b, where each row of A and b corresponds to the input-output invariant evaluated at a different time point. The components of κ are the coefficient functions of the parameters in the input-output invariants. The set-up of the model compatibility test is given in Section 3.Then we ask: does there exist aκ such that Aκ=b. If b=0, of course we are guaranteed a zero trivial solution and the nontrivial solution can be determined via a rank test (e.g., singular value decomposition, or SVD). Since data may be imperfect, we can perform the statistical criterion developed in [10] with the bound improved in [11] to determine whether or not to reject the model. One of the key differences in adapting this method to time-course data is considering Aκ=b when b≠0. For Aκ=b, there may be no solutions. Thus, we must check if the linear system of equations Aκ=b is consistent, i.e., has one or infinitely many solutions. We present a rank test, based on the SVD, for determining the compatibility of the data with input-output invariants from various (potentially incorrect) models. The linear solvability test is described in Section 4. We assume our data have Gaussian measurement noise. In Section 5, we derive a statistical cut-off for when the model is incompatible with the data.Another key difference in this approach than previous parameter-free model discrimination methods is the occurrence of higher order derivatives of the input and output variables in the input-output invariants, requiring them to be known at various time instances. Often one does not have data points for the higher order derivatives, then these need to be estimated. Unlike numerical estimation or splines, which assume a specific functional form, one can use Gaussian Process Regression (GPR) to estimate the higher order derivatives from time-course data. In Section6, we present such a method, which has previously been done for first and second derivatives of biological data [13]. Bounding error of derivative estimates is a difficult problem, which requires us to remove certain data points; however, the advantage of GPR is that one can consider a family of functions, which [13] points out to be able to capture many more temporal trends in the data than any one equation. This enables us to substitute the newly estimated derivative data into the input-output invariant and test model compatibility using the solvability test with the statistical cut-off that we present in Sections 4 and 5.In Sections4 and 7, we showcase our method with examples from linear and nonlinear models. Finally we discuss special cases and other related topics in Section 8, before concluding in Section 9.


## 2. Differential Elimination

We now give some background on differential algebra since a crucial step in our algorithm is to perform differential elimination to obtain equations purely in terms of input variables, output variables, and parameters. For this reason, we will only give background on the ideas from differential algebra required to understand the differential elimination process. For a more detailed description of differential algebra and the algorithms listed below, see [12, 14, 15]. In what follows, we assume the reader is familiar with concepts such asrings andideals, which are covered in great detail in [8].Definition 1.
A ringS is said to be adifferential ring if there is a derivative defined on S and S is closed under differentiation. Adifferential ideal is an ideal which is closed under differentiation.Let our differential ideal be equipped with aranking, i.e., a total ordering, denoted <, among the variables and their derivatives. Let zi(μ) and zj(ν) be arbitrary derivatives. Then the ranking should be such that, for arbitrary positive integer k:(1)ziν<ziν+k,ziμ<zjν⇒ziμ+k<zjν+kLetuj be the leader of a polynomial Aj, which is the highest ranking derivative of the variables appearing in that polynomial. A polynomial Ai is said to be oflower rank than Aj if the order of ui is less than the order of uj or, whenever ui=uj, the highest algebraic degree of any term containing the leader of Ai is less than the highest algebraic degree of any term containing the leader of Aj. A polynomial Ai isreduced with respect to a polynomialAj if Ai contains neither the leader of Aj with equal or greater algebraic degree, nor its derivatives. If Ai is not reduced with respect to Aj, it can be reduced by using the pseudodivision algorithm in Section 2.1. A set of differential polynomials A=A1,A2,…,Ar that are all reduced with respect to each other is called anauto-reduced set.Two auto-reduced sets,A=A1,A2,…,Ar and B=B1,B2,…,Bs ordered in increasing rank so that A1<A2<…<Ar,B1<B2<…<Bs, are ranked according to the following principle: if there is an integer k, k≤min⁡s,r such that rank Ai = rank Bi, i=1,…,k-1, rank Ak< rank Bk, then A is said to be lower rank than B. If r<s and rank Ai = rank Bi, i=1,…,r, then A is also said to be of lower rank than B.A useful description of a differential ideal is called adifferential characteristic set, which is a finite description of a possibly infinite set of differential polynomials. We give the technical definition from [12].Definition 2.
LetΣ be a set of differential polynomials, not necessarily finite. If A⊂Σ is an auto-reduced set, such that no lower ranked auto-reduced set can be formed in Σ, then A is called adifferential characteristic set.A well-known fact in differential algebra is that differential ideals need not be finitely generated [12, 15]. However, a radical differential ideal is finitely generated by theRitt-Raudenbush basis theorem [16]. This result gives rise to Ritt’s pseudodivision algorithm (see below), allowing us to compute the differential characteristic set of a radical differential ideal. We now describe various methods to find a differential characteristic set and other related notions, and we describe why they are relevant to our problem; namely, they can be used to find theinput-output equations.In what follows, we will be considering the differential ringR(p)[u,y,x], where R(p) is the field of rational functions in the parameter vector p. The variables in this differential ring are the states, the inputs, the outputs, and possibly their derivatives.Consider an ODE system of the formx˙(t)=f(x(t),p,u(t)) and yj(t)=gj(x(t),p) for j=1,…,M with f and g rational functions of their arguments. Let our differential ideal be generated by the differential polynomials obtained by subtracting the right-hand-side from the ODE system to obtain x˙(t)-f(x(t),p,u(t)) and yj(t)-gj(x(t),p) for j=1,…,M. In what follows, we use the ranking in [17], which is given by(2)u<u˙<u¨<⋯<y<y˙<y¨<⋯<x1<x2<⋯<x˙1<x˙2<⋯Note that the notation reflects the fact that the ordering among the components of u and y is immaterial, since these are known variables, whereas different ordering of the components of x may lead to different characteristic sets [17]. With respect to this ordering, a differential characteristic set is of the form [17]:(3)A1u,y,…,AMu,yAM+1u,y,x1AM+2u,y,x1,x2⋮AM+Nu,y,x1,…,xNwhere Ai are differential polynomials. Note that the resulting system is not necessarily auto-reduced in R(p)[u,y,x], namely, A1(u,y),…,AM(u,y) may not be auto-reduced. The first M terms of the differential characteristic set, A1(u,y),…,AM(u,y), are those terms independent of the state variables and when set to zero form theinput-output equations:(4)Fu,u˙,u¨,u⃛,…,y,y˙,y¨,y⃛,…=0.Specifically, the M input-output equations F(u,u˙,u¨,u⃛,…,y,y˙,y¨,y⃛,…)=0 are polynomial equations in the variables u,u˙,u¨,u⃛,…,y,y˙,y¨,y⃛,… with rational coefficients in the parameter vector p. Note that the differential characteristic set is in general non-unique, but the coefficients of the input-output equations can be fixed uniquely by normalizing the equations to make them monic.We now discuss several methods to find the input-output equations. The first method (Ritt’s pseudodivision algorithm) can be used to find a differential characteristic set for a radical differential ideal. The second method (RosenfeldGroebner) gives a representation of the radical of the differential ideal as an intersection of regular differential ideals and can also be used to find a differential characteristic set under certain conditions [18, 19]. Finally, we discuss Gröbner basis methods to find theinput-output equations.


### 2.1. Ritt’s Pseudodivision Algorithm

An algorithm to find a differential characteristic set of a radical (in particular, prime) differential ideal generated by a finite set of differential polynomals is called Ritt’s pseudodivision algorithm. We describe the process in detail below, which comes from the description in [17]. Note that our differential ideal as described above is a prime differential ideal [12, 20]. Let Ai and Aj be differential polynomials.(1)
IfAi contains the kth derivative uj(k) of the leader of Aj, Aj is differentiated k times so its leader becomes uj(k).(2)
Multiply the polynomialAi by the coefficient of the highest power of uj(k); let R be the remainder of the division of this new polynomial by Aj(k) with respect to the variable uj(k). Then R is reduced with respect to Aj(k). The polynomial R is called thepseudoremainder of the pseudodivision.(3)
The polynomialAi is replaced by the pseudoremainder R and the process is iterated using Aj(k-1) in place of Aj(k) and so on, until the pseudoremainder is reduced with respect to Aj.This algorithm is applied to a set of differential polynomials, such that each polynomial is reduced with respect to each other, to form an auto-reduced set. The result is a differential characteristic set. Note that the multiplication mentioned in Step (2) above may yield a nonequivalent system if that coefficient happens to belong to the ideal. However, in practice, this does not occur for the ODE systems studied [17].


### 2.2. RosenfeldGroebner

Using theDifferentialAlgebra package in Maple, one can find a representation of the radical of a differential ideal generated by some equations, as an intersection of radical differential ideals with respect to a given ranking [21]. Specifically, the RosenfeldGroebner command in Maple takes two arguments: sys and R, where sys is a list of set of differential equations or inequations which are all rational in the independent and dependent variables and their derivatives and R is a differential polynomial ring built by the commandDifferentialRing specifying the independent and dependent variables and a ranking for them [21]. Then RosenfeldGroebner returns a representation of the radical of the differential ideal generated by sys, as an intersection of radical differential ideals saturated by the multiplicative family generated by the inequations found in sys. This representation consists of a list of regular differential chains with respect to the ranking of R. Note that RosenfeldGroebner returns a differential characteristic set if the differential ideal is prime [18].


### 2.3. Gröbner Basis Methods

Finally, both algebraic and differential Gröbner bases can be employed to find the input-output equations. To use an algebraic Gröbner basis, one can take a sufficient number of derivatives of the model equations and then treat the derivatives of the variables as indeterminates in the polynomial ring inx, x˙, x¨,..., u, u˙, u¨,..., y, y˙, y¨,..., etc. Then a Gröbner basis of the ideal generated by this full system of (differential) equations with an elimination ordering where the state variables and their derivatives are eliminated first can be found. Details of this approach can be found in [22]. Differential Gröbner bases have been developed by Carrà Ferro [23], Ollivier [24], and Mansfield [25], but currently there are no implementations in computer algebra systems [14].


## 2.1. Ritt’s Pseudodivision Algorithm

An algorithm to find a differential characteristic set of a radical (in particular, prime) differential ideal generated by a finite set of differential polynomals is called Ritt’s pseudodivision algorithm. We describe the process in detail below, which comes from the description in [17]. Note that our differential ideal as described above is a prime differential ideal [12, 20]. Let Ai and Aj be differential polynomials.(1)
IfAi contains the kth derivative uj(k) of the leader of Aj, Aj is differentiated k times so its leader becomes uj(k).(2)
Multiply the polynomialAi by the coefficient of the highest power of uj(k); let R be the remainder of the division of this new polynomial by Aj(k) with respect to the variable uj(k). Then R is reduced with respect to Aj(k). The polynomial R is called thepseudoremainder of the pseudodivision.(3)
The polynomialAi is replaced by the pseudoremainder R and the process is iterated using Aj(k-1) in place of Aj(k) and so on, until the pseudoremainder is reduced with respect to Aj.This algorithm is applied to a set of differential polynomials, such that each polynomial is reduced with respect to each other, to form an auto-reduced set. The result is a differential characteristic set. Note that the multiplication mentioned in Step (2) above may yield a nonequivalent system if that coefficient happens to belong to the ideal. However, in practice, this does not occur for the ODE systems studied [17].


## 2.2. RosenfeldGroebner

Using theDifferentialAlgebra package in Maple, one can find a representation of the radical of a differential ideal generated by some equations, as an intersection of radical differential ideals with respect to a given ranking [21]. Specifically, the RosenfeldGroebner command in Maple takes two arguments: sys and R, where sys is a list of set of differential equations or inequations which are all rational in the independent and dependent variables and their derivatives and R is a differential polynomial ring built by the commandDifferentialRing specifying the independent and dependent variables and a ranking for them [21]. Then RosenfeldGroebner returns a representation of the radical of the differential ideal generated by sys, as an intersection of radical differential ideals saturated by the multiplicative family generated by the inequations found in sys. This representation consists of a list of regular differential chains with respect to the ranking of R. Note that RosenfeldGroebner returns a differential characteristic set if the differential ideal is prime [18].


## 2.3. Gröbner Basis Methods

Finally, both algebraic and differential Gröbner bases can be employed to find the input-output equations. To use an algebraic Gröbner basis, one can take a sufficient number of derivatives of the model equations and then treat the derivatives of the variables as indeterminates in the polynomial ring inx, x˙, x¨,..., u, u˙, u¨,..., y, y˙, y¨,..., etc. Then a Gröbner basis of the ideal generated by this full system of (differential) equations with an elimination ordering where the state variables and their derivatives are eliminated first can be found. Details of this approach can be found in [22]. Differential Gröbner bases have been developed by Carrà Ferro [23], Ollivier [24], and Mansfield [25], but currently there are no implementations in computer algebra systems [14].


## 3. Model Rejection Using Input-Output Invariants

We now discuss how to use the input-output invariants obtained from differential elimination (using Ritt’s pseudodivision, differential Gröbner bases, or some other method) for model selection/rejection.We can write our input-output relations in (4), or input-output invariants, in the form:(5)∑icipψiu,y=0The functions ψi(u,y) are differential monomials, i.e., monomials in the input/output variables u,u˙,u¨,u⃛, …, y,y˙,y¨,y⃛,…, etc., and the functions ci(p) are rational functions in the unknown parameter vector p. In order to uniquely fix the rational coefficients ci(p) to the differential monomials ψi(u,y), we normalize each input/output equation to make it monic. In other words, we can rewrite our input-output relations as(6)∑ic~ipψiu,y=ξu,yHere ξ(u,y) is a differential monomial in the input/output variables u,u˙,u¨,u⃛, …, y,y˙,y¨,y⃛,…, etc. If the values of u,u˙,u¨,u⃛,…, y,y˙,y¨,y⃛,…, etc., were known at a sufficient number of time instances t1,t2,…,tm, then one could substitute in values of ψi(u,y) and ξ(u,y) at each of these time instances to obtain a linear system of equations in the variables c~i(p).First consider the case of a single input-output equation. If there aren unknown coefficients c~i(p), we obtain the system: (7)c~1pψ1ut1,yt1+⋯+c~npψnut1,yt1=ξut1,yt1⋮c~1pψ1utm,ytm+⋯+c~npψnutm,ytm=ξutm,ytmWe write this linear system asAκ=b, where A is an m by n matrix of the form:(8)ψ1ut1,yt1⋯ψnut1,yt1⋮⋮⋮ψ1utm,ytm⋯ψnutm,ytmκ is the vector of unknown coefficients [c~1(p),…,c~n(p)]T, and b is of the form [ξ(u(t1),y(t1)),…,ξ(u(tm),y(tm))]T.For the case of multiple input-output equations, we get the following block diagonal system of equationsAκ=b:(9)A100…00A20…0⋮⋮⋱⋮⋮000…AMκ1κ2⋮κM=b1b2⋮bMwhere A is a m=m1+⋯+mM by n=n1+⋯+nM matrix.In the symbolic setting, given a general input functionu and generic initial conditions and parameters, this system Aκ=b should have a unique solution for κ, due to the persistence of excitation conditions [26]. In other words, we assume that the vectors of differential monomials ψ1,…,ψn at various time points are linearly independent. This means the coefficients c~i(p) of the input-output equations can be uniquely determined in the generic setting [26]. Note that we have assumed that the parameters are all unknown and we have not taken any possible algebraic dependencies among the coefficients into account.The main idea of this paper is to translate the symbolic setting to the numerical setting and can be described as follows. Assume we have perfect data; i.e., we know values ofu,u˙,u¨,u⃛,…,y,y˙,y¨,y⃛,…, etc., at many time instances t1,…,tm, perfectly. Given a set of candidate models, we find their associated input-output invariants and then substitute in our values of u,u˙,u¨,u⃛,…,y,y˙,y¨,y⃛,…, etc., at time instances t1,…,tm, thus setting up the linear system Aκ=b for each model. With perfect data and assuming the persistence of excitation conditions mentioned above, the solution to Aκ=b should be unique for the correct model, but there should, in theory, be no solution for each of the incorrect models. Thus under ideal circumstances, one should be able to select the correct model since the input/output data corresponding to that model should satisfy its input-output invariant. Likewise, one should be able to reject the incorrect models since the input/output data should not satisfy their input-output invariants.However, with imperfect data, there could be no solution toAκ=b even for the correct model, and likewise there may or may not be a solution to Aκ=b for an incorrect model. Thus, with imperfect data, one may be unable to select the correct model. On the other hand, if there is no solution to Aκ=b for each of the candidate models, then the goal is to determine how “badly” each of the models fails and rejects models accordingly.A subtle point regarding this approach is that this model rejection technique works best if the models under consideration are in the simplest possible form. This means that, ideally, redundant parameters have been eliminated from the model so that the input-output equations are as reduced as possible; i.e., there are not extra columns inA when considering the linear system Aκ=b. Extra columns generically mean more possible solutions, which can make it harder for our algorithm to reject incorrect models. However, redundant parameters, and the related notion of parameter unidentifiability, do not necessarily yield more coefficients in the input-output equations, as can be seen from the structure of input-output equations for linear compartment models as discussed in Section 8.A related question to model compatibility is that of structural indistinguishability. Two models are structurally indistinguishable if for any choice of parameters in the first model there is a choice of parameters in the second model that will yield the same dynamics in both models, and vice versa [27]. One way to test for structural indistinguishability of two models is to find the associated input-output equations and then equate their coefficient functions and attempt to solve for one set of parameters in terms of the other set of parameters, and vice versa [27]. A necessary condition for models to be structurally indistinguishable is to have input-output equations with the same differential monomial terms. Since our approach only considers the structure of the input-output equations and not the specific coefficient functions, it is possible to have several different models, all with the same structure of their input-output equations, to be compatible using our model compatibility test. Thus, if a given model is found to be compatible, then any model that is structurally indistinguishable from that model is also compatible and thus our approach and structural indistinguishability testing can be applied in parallel. For more on structural indistinguishability, see [28–30]. The specific form of the coefficients of the input-output equations is considered in Section 8.We now describe criteria to reject models.


## 4. Linear Solvability

LetA∈Rm×n and consider the linear system(10)Aκ=b.Here, we study the solvability of (10) under noisy perturbation of both A and b. Let A~ and b~ denote the perturbed versions of A and b, respectively, and assume that A~-A and b~-b depend only on A~ and b~, respectively (see Section 5). Our goal is to infer theunsolvability of the unperturbed system (10) from observation of A~ and b~ only.Our method is based on detecting the rank of an augmented matrix, but first let us introduce some notation. The singular values of a matrixA∈Rm×n will be denoted by (11)σ1A≥⋯≥σlA≥σl+1A=⋯=σnA=0,l=min⁡m,n.(Note that we have trivially extended the number of singular values of A from l to n.) The rank of A is written rank(A). The range of A is denoted R(A). Throughout, · refers to the Euclidean norm.The basic strategy will be to assume as a null hypothesis that (10) has a solution, i.e., b∈R(A), and then to derive its consequences in terms of A~ and b~. If these consequences are not met, then we conclude by contradiction that (10) is unsolvable. In other words, we will providesufficient but not necessary conditions for (10) to have no solution; i.e., we can only reject (but not confirm) the null hypothesis. We will refer to this procedure astesting the null hypothesis.


### 4.1. Preliminaries

We first collect some useful results.Theorem 3 (Weyl’s inequality).
LetA,A~∈Rm×n. Then (12)σkA~-σkA≤A~-A,k=1,…,n.Corollary 4.
LetA,A~∈Rm×n and assume that rank(A)<k. Then(13)σkA~≤A~-A.Therefore, if (13) is not satisfied, then rank(A)≥k.


### 4.2. Augmented Matrix

Assume the null hypothesis. Thenb∈R(A), so rank([A,b])=rank(A)≤min(m,n). Therefore, σn+1([A,b])=0. But we do not have access to [A,b] and so must consider instead the perturbed augmented matrix [A~,b~].Theorem 5.
Under the null hypothesis,(14)σn+1A~,b~≤A~-A,b~-b≤A~-A+b~-b.Proof.
Apply Corollary4.In other words, if (14) does not hold, then (10) has no solution.Remark 6.
This approach can fail to correctly reject the null hypothesis ifA is (numerically) low-rank.Remark 7.
In principle, we should test directly the assertion thatrank([A,b])=rank(A). However, we can only establish lower bounds on the matrix rank (we can only tell if a singular value is “too large”), so this is not feasible in practice. An alternative approach is to consider onlynumerical ranks obtained by thresholding. How to choose such a threshold, however, is not at all clear and can be a very delicate matter especially if the data have high dynamic range.Remark 8.
The theorem is uninformative ifm≤n since then σn+1([A,b])=σn+1([A~,b~])=0 trivially. However, this is not a significant disadvantage beyond that described above since if A is full-rank, then it must be true that (10) is solvable.


### 4.3. Example: Perfect Data

As a proof of principle, we first apply Theorem5 to a simple linear model. We start by taking perfect input and output data and then add a specific amount of noise to the output data and attempt to reject the incorrect model. In the subsequent sections, we will see how to interpret Theorem 5 statistically under a particular “noise” model for the perturbations.Here, we take data from a linear 3-compartment model, add noise, and try to reject the general form of the linear 2-compartment model with the same input/output compartments. Linear compartment models are defined in Section8. In practice, one would like to compare models with the same input and output compartments, as the number of compartments involved may not be known, but the injection and measurement compartments would be known. Thus, in this particular example, we are assuming the linear 3-compartment model is the “true model” and want to reject a competing model, but as our method concerns model rejection and not model selection, this notion of a “true model” is not a requirement for our method to work.Example 9.
Let our model be a 3-compartment model of the following form:(15)x˙1x˙2x˙3=-2101-3101-2x1x2x3+2e-3t+12e-5t00,y=x1x10=1,x20=7,x30=9Here we have an input to the first compartment of the form u1=2e-3t+12e-5t and the first compartment is measured, so that y=x1 represents the output. Note that we have chosen a smooth, persistently exciting input function [26] so that derivatives can be taken and the coefficients of the input-output equation can be uniquely determined, as required. The solution to this system of ODEs can be easily found of the form:(16)x1x2x3=7111e-t+-101e-2t+1-21e-4t+-1-11e-3t+-53-1e-5tso thaty=7e-t-e-2t+e-4t-e-3t-5e-5t.The input-output equation for a 3-compartment model with a single input/output to the first compartment has the form:(17)y⃛+c1y¨+c2y˙+c3y=u¨1+c4u˙1+c5u1where c1,c2,c3 are the coefficients of the characteristic polynomial of the matrix A and c4,c5 are the coefficients of the characteristic polynomial of the matrix A1 which has the first row and first column of A removed [31].We now substitute values ofu1,u˙1,u¨1,y,y˙,y¨,y⃛ at time instances t=0,0.2,0.4,0.6,0.8,1 into our input-output equation and solve the resulting linear system of equations for c1,c2,c3,c4,c5. We get that c1=7,c2=14,c3=8,c4=5,c5=5, which agrees with the coefficients of the characteristic polynomials of A and A1.We now attempt to reject the 2-compartment model using 3-compartment model data. We find the input-output equations for a 2-compartment model with a single input/output to the first compartment, which has the form:(18)y¨+C1y˙+C2y=u˙1+C3u1where again C1,C2 are the coefficients of the characteristic polynomial of the matrix A and C3 is the coefficient of the characteristic polynomial of the matrix A1 which has the first row and first column of A removed.We substitute values ofu1,u˙1,y,y˙,y¨ at time instances t=0,0.2,0.4,0.6,0.8,1 into our input-output equation and attempt to solve the resulting linear system of equations for C1,C2,C3.The singular values for the matrixA with the substituted values of u1,y,y˙ at time instances t=0,0.2,0.4,0.6,0.8,1 are(19)24.8133,7.64917,0.0626589The singular values of the matrix [A,b] with the substituted values of u1,u˙1,y,y˙,y¨ at time instances t=0,0.2,0.4,0.6,0.8,1 are(20)57.174,7.69381,0.326204,0.00596031Since the smallest singular value is greater than zero (or order machine precision), it is evident that the 2-compartment model can be rejected.We now add noise to our matrixA in the following way. To each entry y˙, and y, we add ϵkij where kij is a random real number between 0 and 1, and ϵ=0.001. Then the noisy matrix A~ has the following singular values:(21)24.8134,7.64949,0.0627827We add noise to our vector b in a similar way. To each entry u˙1-y¨, we add ϵkij where kij is again a random real number between 0 and 1. Then the noisy matrix [A~,b~] has the following singular values:(22)57.1747,7.69409,0.326141,0.00579117To determine whether the noisy data are compatible, we need to compute [A~-A,b~-b]. Due to the specific noise model chosen, this can be bounded independently of the true unobservable data [A,b] as ϵ1=0.00489898, where 1 is a matrix of all ones of the appropriate size (the actual norm is 0.00207018). Since this norm isless than the smallest singular value 0.00579117, we can reject this model. Thus, using noisy 3-compartment model data, we are able to reject the 2-compartment model.


## 4.1. Preliminaries

We first collect some useful results.Theorem 3 (Weyl’s inequality).
LetA,A~∈Rm×n. Then (12)σkA~-σkA≤A~-A,k=1,…,n.Corollary 4.
LetA,A~∈Rm×n and assume that rank(A)<k. Then(13)σkA~≤A~-A.Therefore, if (13) is not satisfied, then rank(A)≥k.


## 4.2. Augmented Matrix

Assume the null hypothesis. Thenb∈R(A), so rank([A,b])=rank(A)≤min(m,n). Therefore, σn+1([A,b])=0. But we do not have access to [A,b] and so must consider instead the perturbed augmented matrix [A~,b~].Theorem 5.
Under the null hypothesis,(14)σn+1A~,b~≤A~-A,b~-b≤A~-A+b~-b.Proof.
Apply Corollary4.In other words, if (14) does not hold, then (10) has no solution.Remark 6.
This approach can fail to correctly reject the null hypothesis ifA is (numerically) low-rank.Remark 7.
In principle, we should test directly the assertion thatrank([A,b])=rank(A). However, we can only establish lower bounds on the matrix rank (we can only tell if a singular value is “too large”), so this is not feasible in practice. An alternative approach is to consider onlynumerical ranks obtained by thresholding. How to choose such a threshold, however, is not at all clear and can be a very delicate matter especially if the data have high dynamic range.Remark 8.
The theorem is uninformative ifm≤n since then σn+1([A,b])=σn+1([A~,b~])=0 trivially. However, this is not a significant disadvantage beyond that described above since if A is full-rank, then it must be true that (10) is solvable.


## 4.3. Example: Perfect Data

As a proof of principle, we first apply Theorem5 to a simple linear model. We start by taking perfect input and output data and then add a specific amount of noise to the output data and attempt to reject the incorrect model. In the subsequent sections, we will see how to interpret Theorem 5 statistically under a particular “noise” model for the perturbations.Here, we take data from a linear 3-compartment model, add noise, and try to reject the general form of the linear 2-compartment model with the same input/output compartments. Linear compartment models are defined in Section8. In practice, one would like to compare models with the same input and output compartments, as the number of compartments involved may not be known, but the injection and measurement compartments would be known. Thus, in this particular example, we are assuming the linear 3-compartment model is the “true model” and want to reject a competing model, but as our method concerns model rejection and not model selection, this notion of a “true model” is not a requirement for our method to work.Example 9.
Let our model be a 3-compartment model of the following form:(15)x˙1x˙2x˙3=-2101-3101-2x1x2x3+2e-3t+12e-5t00,y=x1x10=1,x20=7,x30=9Here we have an input to the first compartment of the form u1=2e-3t+12e-5t and the first compartment is measured, so that y=x1 represents the output. Note that we have chosen a smooth, persistently exciting input function [26] so that derivatives can be taken and the coefficients of the input-output equation can be uniquely determined, as required. The solution to this system of ODEs can be easily found of the form:(16)x1x2x3=7111e-t+-101e-2t+1-21e-4t+-1-11e-3t+-53-1e-5tso thaty=7e-t-e-2t+e-4t-e-3t-5e-5t.The input-output equation for a 3-compartment model with a single input/output to the first compartment has the form:(17)y⃛+c1y¨+c2y˙+c3y=u¨1+c4u˙1+c5u1where c1,c2,c3 are the coefficients of the characteristic polynomial of the matrix A and c4,c5 are the coefficients of the characteristic polynomial of the matrix A1 which has the first row and first column of A removed [31].We now substitute values ofu1,u˙1,u¨1,y,y˙,y¨,y⃛ at time instances t=0,0.2,0.4,0.6,0.8,1 into our input-output equation and solve the resulting linear system of equations for c1,c2,c3,c4,c5. We get that c1=7,c2=14,c3=8,c4=5,c5=5, which agrees with the coefficients of the characteristic polynomials of A and A1.We now attempt to reject the 2-compartment model using 3-compartment model data. We find the input-output equations for a 2-compartment model with a single input/output to the first compartment, which has the form:(18)y¨+C1y˙+C2y=u˙1+C3u1where again C1,C2 are the coefficients of the characteristic polynomial of the matrix A and C3 is the coefficient of the characteristic polynomial of the matrix A1 which has the first row and first column of A removed.We substitute values ofu1,u˙1,y,y˙,y¨ at time instances t=0,0.2,0.4,0.6,0.8,1 into our input-output equation and attempt to solve the resulting linear system of equations for C1,C2,C3.The singular values for the matrixA with the substituted values of u1,y,y˙ at time instances t=0,0.2,0.4,0.6,0.8,1 are(19)24.8133,7.64917,0.0626589The singular values of the matrix [A,b] with the substituted values of u1,u˙1,y,y˙,y¨ at time instances t=0,0.2,0.4,0.6,0.8,1 are(20)57.174,7.69381,0.326204,0.00596031Since the smallest singular value is greater than zero (or order machine precision), it is evident that the 2-compartment model can be rejected.We now add noise to our matrixA in the following way. To each entry y˙, and y, we add ϵkij where kij is a random real number between 0 and 1, and ϵ=0.001. Then the noisy matrix A~ has the following singular values:(21)24.8134,7.64949,0.0627827We add noise to our vector b in a similar way. To each entry u˙1-y¨, we add ϵkij where kij is again a random real number between 0 and 1. Then the noisy matrix [A~,b~] has the following singular values:(22)57.1747,7.69409,0.326141,0.00579117To determine whether the noisy data are compatible, we need to compute [A~-A,b~-b]. Due to the specific noise model chosen, this can be bounded independently of the true unobservable data [A,b] as ϵ1=0.00489898, where 1 is a matrix of all ones of the appropriate size (the actual norm is 0.00207018). Since this norm isless than the smallest singular value 0.00579117, we can reject this model. Thus, using noisy 3-compartment model data, we are able to reject the 2-compartment model.


## 5. Statistical Inference

We now consider the statistical inference of the solvability of (10). First, we need a noise model.


### 5.1. Noise Model

If the perturbationsA~-A and b~-b are bounded, e.g., A~-A≤ϵA~ and b~-b≤ϵb~ for some ϵ>0 (representing a relative accuracy of ϵ in the “measurements” A~ and b~), then Theorem 5 can be used at once. However, it is customary to model such perturbations as normal random variables, which are not bounded. Here, we will assume a noise model of the form (23)A~-A=CA~∘Z,b~-b=Cb~∘Z,where CA~ is a (computable) matrix that depends only on A~ and similarly with Cb~, A∘B denotes the Hadamard (entrywise) matrix product (A∘B)ij=AijBij, and Z is a matrix-valued random variable whose entries Zij~N(0,1) are independent standard normals.In our application of interest, the entries ofCA~ depend on those of A~ as follows. Let Aij=ϕij(v) for some input vector v, but suppose that we can only observe the “noisy” vector v~=(1+ϵZ)∘v. Then the corresponding perturbed matrix entries are (24)A~ij=ϕijv~=ϕijv+ϵ∑k∇ϕijvkvkZk+Oϵ2,Zk~N0,1.By the additivity formula(25)∑kakZk=∑kak2Z=aZfor standard Gaussians, (26)∑k∇ϕijvkvkZk=∑k∇ϕijv∘vkZk=∇ϕijv∘vZ.Therefore, (27)A~ij=Aij+ϵ∇ϕijv∘vZ+Oϵ2=Aij+ϵ∇ϕijv~∘v~Z+Oϵ2,so, to first order in ϵ, (28)CA~ij=ϵ∇ϕijv~∘v~.An analogous derivation holds for Cb~.The basic strategy is now as follows. Letτ be a test statistic, i.e., σn+1([A~,b~]) in Section 4.2. Then since (29)τω≤CA~∘Z+Cb~∘Zω,where we have made explicit the dependence of both sides on the same underlying random mechanism ω, the (cumulative) distribution function of τ must dominate that of CA~∘Z+Cb~∘Z, i.e., (30)Prτ≤v≥PrCA~∘Z+Cb~∘Z≤v.Thus,(31a)Prτ≥v≤PrCA~∘Z+Cb~∘Z≥v(31b)=∫0∞PrCA~∘Z=tPrCb~∘Z≥v-tdt(31c)=∫0vPr⁡CA~∘Z=tPr⁡Cb~∘Z≥v-tdt+∫v∞Pr⁡CA~∘Z=tdt(31d)≤∫0vPr⁡CA~∘Z≥tPr⁡Cb~∘Z≥v-tdt+Pr⁡CA~∘Z≥v.Using (31a)–(31d), we can associate a p-value to any given realization of τ by referencing upper tail bounds for quantities of the form C∘Z. Recall that τ=0 under the null hypothesis. In a classical statistical hypothesis testing framework, we may therefore reject the null hypothesis if (31d) is at most α, where α is the desired significance level (e.g., α=0.05).


### 5.2. Hadamard Tail Bounds

We now turn to boundingPr(C∘Z≥v), where we will assume that C,Z∈Rm×n. This can be done in several ways.One easy way is to recognize that(32)C∘Z≤C∘ZF≤CFZF,where ·F is the Frobenius norm, so (33)PrC∘Z≥v≤PrC∘ZF≥v≤PrZF≥vCF.But ZF~χmn has a chi distribution with mn degrees of freedom. Therefore, (34)PrZF≥vCF=PrV≥vCF,V~χmn.However, each inequality in (32) can be quite loose; a slightly better approach is to use the inequality [32] (35)C∘Z≤∥min⁡maxiCi,:,maxjC:,jZ,where Ci,: and C:,j denote the ith row and jth column, respectively, of C. The Z term can then be handled using a chi distribution via Z≤ZF as above or directly using a concentration bound (see below). Variations on this undoubtedly exist.Here, we will appeal to a result by Tropp [33]. The following is from Section 4.3 in [33].Theorem 10.
LetC,Z∈Rm×n, where each Zij~N(0,1). Then for any v≥0, (36)PrC∘Z≥v≤m+nexp⁡-v22σ2,σ2=max⁡maxiCi,:2,maxjC:,j2.


### 5.3. Test Statistic Tail Bounds

The bound (31d) for Pr(τ≥v) can then be computed as follows. Let (37)P1v=∫0vPrCA~∘Z≥tPrCb~∘Z≥v-tdt,P2v=PrCA~∘Z≥vso that Pr(τ≥v)≤P1(v)+P2(v). Then by Theorem 10, (38)P1v≤mm+n∫0vexp⁡-12t2σA2+v-t2σb2dt,where σA2 and σb2 are the “variance” parameters in the theorem for CA~ and Cb~, respectively. This simplifies to (39)P1v≤mm+nexp⁡-12v2σA2+σb2∫0vexp⁡-12σA2+σb2σA2σb2t-σA2σA2+σb2v2dton completing the square. Now set (40)σ2=σA2σb2σA2+σb2,α=σA2σA2+σb2so that the integral becomes (41)∫0vexp⁡-12σA2+σb2σA2σb2t-σA2σA2+σb2v2dt=∫0vexp⁡-t-αv22σ2dt.The variable substitution u=(t-αv)/σ then gives (42)∫0vexp⁡-t-αv22σ2dt=σ∫-αv/σ1-αv/σe-u2/2du=2πσΦ1-αvσ-Φ-αvσ,where (43)Φv=12π∫-∞ve-t2/2dtis the standard normal distribution function. Thus,(44)P1v≤2πσmm+nΦ1-αvσ-Φ-αvσexp⁡-12v2σA2+σB2.A similar analysis yields(45)P2v≤m+nexp⁡-v22σA2.Equations (44) and (45) together comprise the probability bound on the null hypothesis that we will use hereafter.


## 5.1. Noise Model

If the perturbationsA~-A and b~-b are bounded, e.g., A~-A≤ϵA~ and b~-b≤ϵb~ for some ϵ>0 (representing a relative accuracy of ϵ in the “measurements” A~ and b~), then Theorem 5 can be used at once. However, it is customary to model such perturbations as normal random variables, which are not bounded. Here, we will assume a noise model of the form (23)A~-A=CA~∘Z,b~-b=Cb~∘Z,where CA~ is a (computable) matrix that depends only on A~ and similarly with Cb~, A∘B denotes the Hadamard (entrywise) matrix product (A∘B)ij=AijBij, and Z is a matrix-valued random variable whose entries Zij~N(0,1) are independent standard normals.In our application of interest, the entries ofCA~ depend on those of A~ as follows. Let Aij=ϕij(v) for some input vector v, but suppose that we can only observe the “noisy” vector v~=(1+ϵZ)∘v. Then the corresponding perturbed matrix entries are (24)A~ij=ϕijv~=ϕijv+ϵ∑k∇ϕijvkvkZk+Oϵ2,Zk~N0,1.By the additivity formula(25)∑kakZk=∑kak2Z=aZfor standard Gaussians, (26)∑k∇ϕijvkvkZk=∑k∇ϕijv∘vkZk=∇ϕijv∘vZ.Therefore, (27)A~ij=Aij+ϵ∇ϕijv∘vZ+Oϵ2=Aij+ϵ∇ϕijv~∘v~Z+Oϵ2,so, to first order in ϵ, (28)CA~ij=ϵ∇ϕijv~∘v~.An analogous derivation holds for Cb~.The basic strategy is now as follows. Letτ be a test statistic, i.e., σn+1([A~,b~]) in Section 4.2. Then since (29)τω≤CA~∘Z+Cb~∘Zω,where we have made explicit the dependence of both sides on the same underlying random mechanism ω, the (cumulative) distribution function of τ must dominate that of CA~∘Z+Cb~∘Z, i.e., (30)Prτ≤v≥PrCA~∘Z+Cb~∘Z≤v.Thus,(31a)Prτ≥v≤PrCA~∘Z+Cb~∘Z≥v(31b)=∫0∞PrCA~∘Z=tPrCb~∘Z≥v-tdt(31c)=∫0vPr⁡CA~∘Z=tPr⁡Cb~∘Z≥v-tdt+∫v∞Pr⁡CA~∘Z=tdt(31d)≤∫0vPr⁡CA~∘Z≥tPr⁡Cb~∘Z≥v-tdt+Pr⁡CA~∘Z≥v.Using (31a)–(31d), we can associate a p-value to any given realization of τ by referencing upper tail bounds for quantities of the form C∘Z. Recall that τ=0 under the null hypothesis. In a classical statistical hypothesis testing framework, we may therefore reject the null hypothesis if (31d) is at most α, where α is the desired significance level (e.g., α=0.05).


## 5.2. Hadamard Tail Bounds

We now turn to boundingPr(C∘Z≥v), where we will assume that C,Z∈Rm×n. This can be done in several ways.One easy way is to recognize that(32)C∘Z≤C∘ZF≤CFZF,where ·F is the Frobenius norm, so (33)PrC∘Z≥v≤PrC∘ZF≥v≤PrZF≥vCF.But ZF~χmn has a chi distribution with mn degrees of freedom. Therefore, (34)PrZF≥vCF=PrV≥vCF,V~χmn.However, each inequality in (32) can be quite loose; a slightly better approach is to use the inequality [32] (35)C∘Z≤∥min⁡maxiCi,:,maxjC:,jZ,where Ci,: and C:,j denote the ith row and jth column, respectively, of C. The Z term can then be handled using a chi distribution via Z≤ZF as above or directly using a concentration bound (see below). Variations on this undoubtedly exist.Here, we will appeal to a result by Tropp [33]. The following is from Section 4.3 in [33].Theorem 10.
LetC,Z∈Rm×n, where each Zij~N(0,1). Then for any v≥0, (36)PrC∘Z≥v≤m+nexp⁡-v22σ2,σ2=max⁡maxiCi,:2,maxjC:,j2.


## 5.3. Test Statistic Tail Bounds

The bound (31d) for Pr(τ≥v) can then be computed as follows. Let (37)P1v=∫0vPrCA~∘Z≥tPrCb~∘Z≥v-tdt,P2v=PrCA~∘Z≥vso that Pr(τ≥v)≤P1(v)+P2(v). Then by Theorem 10, (38)P1v≤mm+n∫0vexp⁡-12t2σA2+v-t2σb2dt,where σA2 and σb2 are the “variance” parameters in the theorem for CA~ and Cb~, respectively. This simplifies to (39)P1v≤mm+nexp⁡-12v2σA2+σb2∫0vexp⁡-12σA2+σb2σA2σb2t-σA2σA2+σb2v2dton completing the square. Now set (40)σ2=σA2σb2σA2+σb2,α=σA2σA2+σb2so that the integral becomes (41)∫0vexp⁡-12σA2+σb2σA2σb2t-σA2σA2+σb2v2dt=∫0vexp⁡-t-αv22σ2dt.The variable substitution u=(t-αv)/σ then gives (42)∫0vexp⁡-t-αv22σ2dt=σ∫-αv/σ1-αv/σe-u2/2du=2πσΦ1-αvσ-Φ-αvσ,where (43)Φv=12π∫-∞ve-t2/2dtis the standard normal distribution function. Thus,(44)P1v≤2πσmm+nΦ1-αvσ-Φ-αvσexp⁡-12v2σA2+σB2.A similar analysis yields(45)P2v≤m+nexp⁡-v22σA2.Equations (44) and (45) together comprise the probability bound on the null hypothesis that we will use hereafter.


## 6. Gaussian Processes to Estimate Derivatives

We next present a method for estimating higher order derivatives and the estimation error using Gaussian Process Regression and then apply the input-output invariant method to both linear and nonlinear models in the subsequent sections.A Gaussian process (GP) is a stochastic processW(t)~N(μ(t),Σ(t,t′)), where μ(t) is a mean function and Σ(t,t′) a covariance function. GPs are often used for regression/prediction as follows.Suppose that there is an underlying deterministic functionw(t) that we can only observe with some measurement noise as w^(t)=w(t)+ϵ(t), where ϵ(t)~N(0,σ2(t)δ(t,t′)) for (46)δt,t′=1ift=t′0ift≠t′the Dirac delta. We consider the problem of finding w(t) in a Bayesian setting by assuming it to be a GP with prior mean and covariance functions μprior and Σprior, respectively. Then the joint distribution of w^(t)=[w^t1,…,w^tp]⊤ at the observation points t=[t1,…,tp]⊤ and w(s)=[w(s1),…,w(sq)]⊤ at the prediction points s=[s1,…,sq]⊤ is(47)w^tws~Nμpriortμpriors,Σpriort,t+σ2tIΣprior⊤s,tΣpriors,tΣpriors,sThe conditional distribution of w(s) given w(t)=w^(t) is also Gaussian:(48)ws∣wt=w^t~Nμpost,Σpost,where (49)μpost=μpriors+Σpriors,tΣpriort,t+σ2tI-1w^t-μpriort,Σpost=Σpriors,s-Σpriors,tΣpriort,t+σ2tI-1)Σprior⊤s,tare the posterior mean and covariance, respectively. This allows us to infer w(s) on the basis of observing w^(t). The diagonal entries of Σpost are the posterior variances and quantify the uncertainty associated with this inference procedure. In particular, the square roots of these variances give us estimates on the ϵw term in the assumed noise model w^=(1+ϵZ)∘w in Section 5.


### 6.1. Estimating Derivatives

Equation (48) provides an estimate for the function values w(s). What if we want to estimate its derivatives? Let cov(w(t),w(t′))=k(t,t′) for some covariance function k. Then cov(wm(t),wn(t′))=∂tm∂t′nk(t,t′) by linearity of differentiation. Thus, (50)
w
^
t
w
s
w
′
s
⋮
w
n
s
~
N
μ
p
r
i
o
r
t
μ
p
r
i
o
r
s
μ
p
r
i
o
r
1
s
⋮
μ
p
r
i
o
r
n
s
,
Σ
p
r
i
o
r
t
,
t
+
σ
2
t
I
Σ
p
r
i
o
r
⊤
s
,
t
Σ
p
r
i
o
r
1,0
,
⊤
s
,
t
⋯
Σ
p
r
i
o
r
n
,
0
,
⊤
s
,
t
Σ
p
r
i
o
r
s
,
t
Σ
p
r
i
o
r
s
,
s
Σ
p
r
i
o
r
1,0
,
⊤
s
,
s
⋯
Σ
p
r
i
o
r
n
,
0
,
⊤
s
,
s
Σ
p
r
i
o
r
1,0
s
,
t
Σ
p
r
i
o
r
1,0
s
,
s
Σ
p
r
i
o
r
1,1
s
,
s
⋯
Σ
p
r
i
o
r
n
,
1
,
⊤
s
,
s
⋮
⋮
⋮
⋱
⋮
Σ
p
r
i
o
r
n
,
0
s
,
t
Σ
p
r
i
o
r
n
,
0
s
,
s
Σ
p
r
i
o
r
n
,
1
s
,
s
⋯
Σ
n
,
n
s
,
s
,where μprior(i)(t) is the prior mean for w(i)(t) and Σpriori,j(t,t′)=∂ti∂t′jΣprior(t,t′). This joint distribution is exactly of the form (47). An analogous application of (48) then yields the posterior estimate of w(i)(s)∣(w(t)=w^(t)) for all i=0,1,…,n.Alternatively, if we are interested only in the posterior variances of eachw(i)(s), then it suffices to consider each 2×2 block independently: (51)w^twis~Nμpriortμprioris,Σpriort,t+σ2tIΣpriori,0,⊤s,tΣpriori,0s,tΣpriori,is,s.The cost of computing (Σprior(t,t)+σ2(t)I)-1 can clearly be amortized over all i.


### 6.2. Formulae for Squared Exponential Covariance Functions

We now consider the specific case of the squared exponential (SE) covariance function(52)kt,t′=θ2exp⁡-t-t′22l2,where θ2 is the signal variance and l is a length scale. The SE function is one of the most widely used covariance functions in practice. Its derivatives can be expressed in terms of the (probabilists’) Hermite polynomials (53)Hnw=-1new2/2dndwne-w2/2(these are also sometimes denoted Hen(w)). The first few Hermite polynomials are H0(w)=1, H1(w)=w, and H2(w)=w2-1.We need to compute the derivatives∂tm∂t′nk(t,t′). Let v=(t-t′)/l so that k(t,t′)=k(v)=θ2e-v2/2. Then ∂tmf(v)=(1/l)mf(m)(v) and ∂t′nf(v)=(-1/l)nf(n)(v). Therefore, (54)∂m∂tm∂n∂t′nkt,t′=-1nlm+nkm+nv=-1mlm+nHm+nvkv=-1mlm+nHm+nt-t′lkt,t′.The GP regression requires us to have the values of the hyperparametersσ2, θ2, and l. In practice, however, these are hardly ever known. In the examples below, we deal with this by estimating the hyperparameters from the data by maximizing the likelihood. We do this by using a nonlinear conjugate gradient algorithm, which can be quite sensitive to the initial starting point, so we initialize multiple runs over a small grid in hyperparameter space and return the best estimate found. This increases the quality of the estimated hyperparameters but can still sometimes fail.


## 6.1. Estimating Derivatives

Equation (48) provides an estimate for the function values w(s). What if we want to estimate its derivatives? Let cov(w(t),w(t′))=k(t,t′) for some covariance function k. Then cov(wm(t),wn(t′))=∂tm∂t′nk(t,t′) by linearity of differentiation. Thus, (50)
w
^
t
w
s
w
′
s
⋮
w
n
s
~
N
μ
p
r
i
o
r
t
μ
p
r
i
o
r
s
μ
p
r
i
o
r
1
s
⋮
μ
p
r
i
o
r
n
s
,
Σ
p
r
i
o
r
t
,
t
+
σ
2
t
I
Σ
p
r
i
o
r
⊤
s
,
t
Σ
p
r
i
o
r
1,0
,
⊤
s
,
t
⋯
Σ
p
r
i
o
r
n
,
0
,
⊤
s
,
t
Σ
p
r
i
o
r
s
,
t
Σ
p
r
i
o
r
s
,
s
Σ
p
r
i
o
r
1,0
,
⊤
s
,
s
⋯
Σ
p
r
i
o
r
n
,
0
,
⊤
s
,
s
Σ
p
r
i
o
r
1,0
s
,
t
Σ
p
r
i
o
r
1,0
s
,
s
Σ
p
r
i
o
r
1,1
s
,
s
⋯
Σ
p
r
i
o
r
n
,
1
,
⊤
s
,
s
⋮
⋮
⋮
⋱
⋮
Σ
p
r
i
o
r
n
,
0
s
,
t
Σ
p
r
i
o
r
n
,
0
s
,
s
Σ
p
r
i
o
r
n
,
1
s
,
s
⋯
Σ
n
,
n
s
,
s
,where μprior(i)(t) is the prior mean for w(i)(t) and Σpriori,j(t,t′)=∂ti∂t′jΣprior(t,t′). This joint distribution is exactly of the form (47). An analogous application of (48) then yields the posterior estimate of w(i)(s)∣(w(t)=w^(t)) for all i=0,1,…,n.Alternatively, if we are interested only in the posterior variances of eachw(i)(s), then it suffices to consider each 2×2 block independently: (51)w^twis~Nμpriortμprioris,Σpriort,t+σ2tIΣpriori,0,⊤s,tΣpriori,0s,tΣpriori,is,s.The cost of computing (Σprior(t,t)+σ2(t)I)-1 can clearly be amortized over all i.


## 6.2. Formulae for Squared Exponential Covariance Functions

We now consider the specific case of the squared exponential (SE) covariance function(52)kt,t′=θ2exp⁡-t-t′22l2,where θ2 is the signal variance and l is a length scale. The SE function is one of the most widely used covariance functions in practice. Its derivatives can be expressed in terms of the (probabilists’) Hermite polynomials (53)Hnw=-1new2/2dndwne-w2/2(these are also sometimes denoted Hen(w)). The first few Hermite polynomials are H0(w)=1, H1(w)=w, and H2(w)=w2-1.We need to compute the derivatives∂tm∂t′nk(t,t′). Let v=(t-t′)/l so that k(t,t′)=k(v)=θ2e-v2/2. Then ∂tmf(v)=(1/l)mf(m)(v) and ∂t′nf(v)=(-1/l)nf(n)(v). Therefore, (54)∂m∂tm∂n∂t′nkt,t′=-1nlm+nkm+nv=-1mlm+nHm+nvkv=-1mlm+nHm+nt-t′lkt,t′.The GP regression requires us to have the values of the hyperparametersσ2, θ2, and l. In practice, however, these are hardly ever known. In the examples below, we deal with this by estimating the hyperparameters from the data by maximizing the likelihood. We do this by using a nonlinear conjugate gradient algorithm, which can be quite sensitive to the initial starting point, so we initialize multiple runs over a small grid in hyperparameter space and return the best estimate found. This increases the quality of the estimated hyperparameters but can still sometimes fail.


## 7. Results

We showcase our method on competing models: linear compartment models (2 and 3 species), Lotka-Volterra models (2 and 3 species) and Lorenz. We compute the input-output invariants of the Lotka-Volterra and Lorenz usingRosenfeldGroebner. The method to compute the linear compartment input-output invariants is presented in the following section. We simulate each of these models to generate time-course data, add varying levels of noise, and estimate the necessary higher order derivatives using GP regression. Using the estimated GP regression data, we test each of the models using the input-output invariant method on other models.Example 11.
The two species Lotka-Volterra model is(55)x˙1=p1x1-p2x1x2,x˙2=-p3x2+p4x1x2,where x1 and x2 are variables and p1,p2,p3,p4 are parameters. We assume only x1 is observable and perform differential elimination and obtain our input-output invariant in terms of only y=x1(t):(56)p4y˙y2-p3y˙y-p1p4y3+p1p3y2=y¨y-y˙2.Example 12.
By including an additional variablez, the three species Lotka-Volterra model is (57)x˙1=p1x1-p2x1x2,x˙2=-p3x2+p4x1x2-p5x2x3,x˙3=-p6x3+p7x2x3,assuming only y=x1 is observable. After differential elimination, the input-output invariant is(58)p12p4p6-p13p4p7p2y5+p13p3p7p2-p12p3p6y4+3p12p4p7p2+p12p4+2p1p4p6y˙y4+2p1p3p6-3p12p3p7p2y˙y3+p4p6-2p1p4-3p1p4p7p2y˙2y3+p12p7+3p1p3p7p2-p3p6-p1p6y˙2y2+p4p7p2+p4y˙3y2+2p1+p6-2p1p7+p3p7p2y˙3y+p7p2y˙4+p1p6-p12p7p2y¨y3+2p1p7p2-3p1-p6y¨y˙y2-p7p2y¨y˙2y+p1y⃛y3=-y¨2y2+y⃛y˙y2-y¨y˙2y+y˙4.Example 13.
Another three species model, the Lorenz model, is described by the system of equations:(59)x˙1=p1x2-x1,x˙2=x1p2-x3-x2,x˙3=x1x2-p3x3,We assume only y=x1 is observable, perform differential elimination, and obtain the following invariant: (60)-p1+p3y¨y+p1y˙2-p1p3+p3y˙y-p1y4+p1p2p3-p1p3y2=y⃛y-y¨y˙+y¨y-y˙2+y˙y3.Example 14.
A linear 2-compartment model without input can be written as(61)x˙1=p11x1+p12x2,x˙2=p21x1+p22x2,where x1 and x2 are variables and p11,p12,p21,p22 are parameters. We assume only x1 is observable and perform differential elimination and obtain our input-output invariant in terms of only y=x1(t):(62)y¨-p11+p22y˙+p11p22-p12p21y=0Example 15.
A linear 3-compartment model without input is(63)x˙1=p11x1+p12x2+p13x3,x˙2=p21x1+p22x2+p23x3,x˙3=p31x1+p32x2+p33x3,where x1,x2,x3 are variables and p11.p12,p13,p21,p22,p23,p31,p32,p33 are parameters. We assume only x1 is observable and perform differential elimination and obtain our input-output invariant in terms of only y=x1(t): (64)y⃛-p11+p22+p33y¨+p12p21-p11p22+p13p31+p23p32-p11p33-p22p33y¨--p13p22p31+p12p23p31+p13p21p32-p11p23p32-p12p21p33+p11p22p33y=0By assumingy=x1 in Examples 6.1–6.5 representing the same observable variable, we apply our method to data simulated from each model and perform model comparison. The models are simulated and 100 time points are obtained for the variable x in each model. We add different levels of Gaussian noise to the simulated data and then estimate the higher order derivatives from the data. Accurate estimation of the derivatives was not always possible. For example, during our study we found that for some parameters of the Lotka-Volterra three species model, e.g., p1,p2,p3,p4,p5,p6,p7=[1.24;1.68;3.26;0.38;1.50;0.15;1.14], the data could not adequately fit with a GP, as indicated by a small likelihood. Furthermore, even when a good fit is achieved, the derivative estimates themselves could be poor as reflected in high posterior variances. This is a notoriously difficult problem and we offer only some pragmatic guidance here. In particular, we err on the side of being overly conservative by keeping only “good” time points defined as follows. Let σpostk(t) be the posterior standard deviation of the estimate of y(k)(t). Then a time point ti is considered good only if(65)σpostkti≤meanσpostktj+stdσpostktjfor all k, where mean(·) and std(·) give the mean and standard deviations, respectively, of a set. In this way, we adaptively filter out potentially problematic inputs to the ensuing model rejection framework. Note that filtering out data is equivalent to removing constraints so that we can only decrease the discriminatory power, i.e., models that are flagged as incompatible after data filtering would have been incompatible as well without filtering.Once the data are obtained and derivative data are estimated through the GP regression, each model data set is tested against the other input-output invariants. Results are shown in Figure1, which gives a probability bound that the data are compatible with a given model (i.e., ~0 means model rejected) at a variety of noise levels. We find that we can reject the three species Lotka-Volterra model and Lorenz model for data simulated from the Lotka-Volterra two species; however, both linear compartment models are compatible. For data from the three species Lotka-Volterra model, the linear compartment models and two species Lotka-Volterra can be rejected until the noise increases and then the method can no longer reject any models. Finally data generated from the Lorenz model can only reject the two species linear compartment and two species Lotka-Volterra model.Figure 1
Data simulated from model specified and differential algebraic statistics model rejection applied to five model invariants. Gaussian noise is added to data in factors of10 at levels shown in the figure. (a) Data simulated from two species Lotka-Volterra model with parameter values p1,p2,p3,p4=[1.24;1.68;3.26;0.38] and initial condition [10,1]. (b) Data simulated from three species Lotka-Volterra model with parameter values p1,p2,p3,p4,p5,p6,p7=[0.178;0.12;0.99;0.17;0.03;0.56;0.88] and initial condition [2,1,1]. (c) Data simulated from the Lorenz model with parameter values p1,p2,p3=[3.5,.3,2.8] and initial condition [2,1,1]. (d) Data simulated from the linear compartment three species model with parameter values p11,p12,p13,p21,p22,p23,p31,p32,p33=[-2,1,0,1,-3,1,0,1,-2] and initial condition [3,1,5].
(a)
Data from Lotka-Volterra two species (b)
Data from Lotka-Volterra three species (c)
Data from Lorenz model (d)
Data from linear compartment three species


## 8. Other Considerations: Known Parameter Values and Algebraic Dependencies

We have demonstrated our model discrimination algorithm on various models. In this section, we consider some other theoretical points regarding input-output invariants.As mentioned earlier, we have assumed that the parameters are all unknown and we have not taken any possible algebraic dependencies among the coefficients into account. This latter point is another reason our algorithm only concerns model rejection and not model selection. Thus, each unknown coefficient is essentially treated as an independent unknown variable in our linear system of equations. However, there may be instances where we would like to consider incorporating this additional information.To analyze the effects of incorporating known parameter values and algebraic dependencies, we will examine a particularly nice class of models, linear compartment models, whose input-output equations can be found using linear algebra techniques [31]; i.e., computation of the input-output equations does not rely on more computationally intensive approaches such as RosenfeldGroebner and Gröbner bases. Since we will now be considering the explicit form of the coefficients of the input-output equations, we describe the set-up of linear compartment models below.LetG=(V,E) be a directed graph with vertex set V and set of directed edges E. Each vertex i∈V corresponds to a compartment in our model and each edge j→i corresponds to a direct flow of material from the jth compartment to the ith compartment. Let In,Out,Leak⊆V be three sets of compartments: the set of input compartments, output compartments, and leak compartments, respectively. To each edge j→i we associate an independent parameter pij, the rate of flow from compartment j to compartment i. To each leak node i∈Leak, we associate an independent parameter p0i, the rate of flow from compartment i leaving the system.To such a graphG and set of leaks Leak we associate the matrix A in the following way:(66)Aij=-p0i-∑k:i→k∈Epkiifi=jandi∈Leak-∑k:i→k∈Epkiifi=jandi∉Leakpijifj→iisanedgeofG0otherwiseThen we construct a system of linear ODEs with inputs and outputs associated with the quadruple (G,In,Out,Leak) as follows:(67)x˙t=Axt+utyit=xitfori∈Outwhere ui(t)≡0 for i∉In. The coordinate functions xi(t) are the state variables, the functions yi(t) are the output variables, and the nonzero functions ui(t) are the inputs. The resulting model is called alinear compartment model.In [31], an explicit formula for the input-output equations for linear models was derived. In particular, it was shown that all linear n-compartment models corresponding to strongly connected graphs with at least one leak and having the same input and output compartments will have the same differential polynomial form of the input-output equations. For example, a linear 2-compartment model with a single input and output in the same compartment and corresponding to a strongly connected graph with at least one leak has the form:(68)y¨+c1y˙+c2y=u˙+c3uThus, our model discrimination method would not work for two distinct linear 2-compartment models with the above-mentioned form. In order to discriminate between two such models, we need to take other information into account, e.g., known parameter values.Example 16.
Consider the following two linear 2-compartment models:(69)x˙1x˙2=-p01-p21p12p21-p12x1x2+u0,y=x1x˙1x˙2=-p21p12p21-p02-p12x1x2+u0,y=x1
whose corresponding input-output equations are of the form:(70)y¨+p01+p21+p12y˙+p01p12y=u˙+p12uy¨+p21+p12+p02y˙+p02p21y=u˙+p02+p12u
Notice that both of these equations are of the above-mentioned form, i.e., both 2-compartment models have a single input and output in the same compartment and correspond to strongly connected graphs with at least one leak. In the first model, there is a leak from the first compartment and an exchange between compartments1 and 2. In the second model, there is a leak from the second compartment and an exchange between compartments 1 and 2. Assume that the parameter p12 is known. In the first model, this changes our invariant to(71)p01+p21y˙+p01p12y=u˙+p12u-y¨-p12y˙or,c1y˙+c2p12y=u˙+p12u-y¨-p12y˙
In the second model, our invariant is(72)p21+p02y˙+p02p21y-p02u=u˙+p12u-y¨-p12y˙or,c1y˙+c2y+c3u=u˙+p12u-y¨-p12y˙
In this case, the right-hand sides of the two equations are the same, but the first equation has two variables (coefficients) while the second equation has three variables (coefficients). Thus, if we had data from the second model, we could try to reject the first model (much like the 3-compartment versus 2-compartment model discrimination in the examples above). In other words, a vector in the span ofy˙,y, and u for t1,t2,t3 may not be in the span of y˙ and y only.We next consider the effect of incorporating coefficient dependency relationships. While we cannot incorporate the polynomial algebraic dependency relationships among the coefficients in our linear algebraic approach to model rejection, we can include certain dependency conditions, such as certain coefficients becoming known constants. We have already seen one way in which this can happen in the previous example (from known nonzero parameter values). We now explore the case where certain coefficients go to zero. From the explicit formula for input-output equations from [31], we get that a linear model without any leaks has a zero term for the coefficient of y. Thus a linear 2-compartment model with a single input and output in the same compartment and corresponding to a strongly connected graph without any leaks has the form:(73)y¨+c1y˙=u˙+c2uThus to discriminate between two distinct linear 2-compartment models, one with leaks and one without any leaks, we should incorporate this zero coefficient into our invariant.Example 17.
Consider the following two linear 2-compartment models:(74)x˙1x˙2=-p01-p21p12p21-p12x1x2+u0,y=x1x˙1x˙2=-p21p12p21-p12x1x2+u0,y=x1
whose corresponding input-output equations are of the form:(75)y¨+p01+p21+p12y˙+p01p12y=u˙+p12uy¨+p21+p12y˙=u˙+p12uIn the first model, there is a leak from the first compartment and an exchange between compartments 1 and 2. In the second model, there is an exchange between compartments 1 and 2 and no leaks. Thus, our invariants can be written as (76)c1y˙+c2y+c3u=u˙-y¨c1y˙+c2u=u˙-y¨
Again, the right-hand sides of the two equations are the same, but the first equation has three variables (coefficients) while the second equation has two variables (coefficients). Thus, if we had data from the first model, we could try to reject the second model. In other words, a vector in the span ofy˙,y, and u for t1,t2,t3 may not be in the span of y˙ and u only.


## 9. Conclusion

After performing this differential algebraic and statistical model rejection, one has already obtained the input-output equations and thus can test structural identifiability [17, 26, 34]. In a sense, our method extends the current spectrum of potential approaches for comparing models with time-course data, in that one can first reject incompatible models, then test structural identifiability of compatible models using input-output equations obtained from the differential elimination, infer parameter values of the admissible models, and apply an information criterion model selection method to assert the best model.Notably the presented differential algebraic and statistical method does not penalize for model complexity, unlike traditional model selection techniques. Rather, we reject when a model cannot, for any parameter values, be compatible with the given data. We found that simpler models, such as the linear 2-compartment model, could be rejected when data were generated from a more complex model, such as the three species Lotka-Volterra model, which elicits a wider range of behavior. On the other hand, more complex models, such as the Lorenz model, were often not rejected, from data simulated from less complex models. In the future it would be helpful to better understand the relationship between input-output invariants and dynamics. It would be useful to develop numerical algorithms in differential algebra (similar to that in numerical algebraic geometry); a natural extension, if such algorithms were available, would be to analyze models with data, although not parameter-free, similar to that done in [35, 36]. Another future direction is creating an algorithm that takes a probabilistic or randomized approach for eliminating variables of larger differential-algebra models [37].We believe there is large scope for additional parameter-free coplanarity model comparison methods. It would be beneficial to explore whether algorithms for differential elimination can handle larger systems and whether this area could be extended.



---

*Source: 6041981-2019-02-17.xml*