Parameter estimation:

local identifiability

JOHN A. JACQUEZ AND TIMOTHY PERRY Departments of Physiology and Biostatistics, University

JACQUEZ,

JOHN

A.,

AND

TIMOTHY PERRY. Parameter of parameters. Am. J. Physiol.

esti-

258 (Endocrinol. Metab. 21): E727-E736,1990.-For biological systems one often cannot set up experiments to measure all of the state variables. If only a subset of the state variables can be measured, it is possible that some of the system parameters cannot influence the measured state variables or that they do so in combinations that do not define the parameters’ effects separately. Such parameters are unidentifiable and are in theory unestimable. Given a model of the system, linear or nonlinear, and initial estimates of the values of all parameters, we exhibit a simple theory and describe a program for checking the local identifiability of the parameters at the initial estimates for given experiments on the model. The program, IDENT, is available from the authors.

mation:

local identifiability

of parameters

of Michigan,

Ann Arbor, Michigan

48109-0622

FIABILITY AND ESTIMABILITY: PRACTICAL CONSIDERATIONS, we consider some of the problems in interpretation of results on identifiability and estimation. In ESTIMATION OF NONIDENTIFIABLE PARAMETERS, we

examine what happens when one uses standard parameter estimation routines to “estimate” nonidentifiable parameters.In DESCRIPTION OF PROGRAMS, we describe the program, IDENT, and in EXAMPLES, we give examples to illustrate the use of the program and expand on the material in IDENTIFIABILITY AND ESTIMABILITY: PRACTICAL

CONSIDERATIONS.

Glossary

B

Input matrix. Inputs to a linear compartmental matrix are written as linear combinations, Bu, of the input vector,

C

Observation matrix, which gives observations on a compartmental system as linear combinations, Cq, of components of compartments vector q. In practice, units are often volume-‘, in which case observations have units of concentration Function describing how ri: depends on x, O*, u, and t. F(x, O*, u, t) has same units as jz Sensitivity matrix. Entry g;j is a partial derivative of response function, G(t, &, with respect to Oj at ith sample point in time and is evaluated at nominal value 0 Sensitivity matrix for identifiable parameters Sensitivity matrix for sensible parameters Matrix of constant fractional transfer coefficients; kG is fraction of compartment j transferred to i per unit time, i.e., a first-order rate constant Vector of compartment sizes. Components have units of mass, usually moles Time A vector of inputs to a system; units are those of jz A volume of distribution; units are liters-’

model

U

of physiological systems are used more and more in the analysis of experiments. That involves estimation of parameters of one or more models of the system. Although there is a large literature on parameter estimation in the biological sciences, it is only in the last two decades that attention has been paid to a problem that properly precedes parameter estimation, the identifiability problem. A rough first description follows. In biological systems one often cannot set up experiments to measure all of the state variables. If only a subset of the state variables can be observed (measured), then it is possible that some of the parameters cannot influence the observed state variables or, if they can, that they do so in combinations that do not permit us to estimate those parameters separately. Such parameters are said to be unidentifiable or nonidentifiable. The work reported here was undertaken to develop a *program to check local identifiability of parameters at a point in parameter space. It was done in the context of compartmental models but can be used generally for models described by systems of ordinary differential equations, linear or nonlinear. In SYSTEM, MODEL, EXPERIMENT and in IDENTIFIABILITY, we introduce notation and the terminology of identifiability in more detail. In CHECKING LOCAL IDENTIFIABILITY AT A POINT, the theory behind the method of checking local identifiability is presented. In IDENTIMATHEMATICAL

MODELS

0193-1849/90

$1.50

Copyright

Fb, o*,u,t) g

0

g1 gs

K

q t U v 1

0 1990 the American

Physiological

Society

E727

Downloaded from www.physiology.org/journal/ajpendo by ${individualUser.givenNames} ${individualUser.surname} (142.103.160.110) on September 25, 2018. Copyright © 1990 American Physiological Society. All rights reserved.

E728

LOCAL

A vector of state variables. Units depend on specific system; for physiological systems, units are commonly those of mass or concentration (mass/volume) Observation vector. A function, G(x, 0’), of state variables and experimental parameters. When solution for x is substituted, it becomes a function of basic parameters and time, the response function. Units depend on system and model Vector of observational parameters. These are parameters that are functions of basic parameters and are identifiable from observations in an experiment. A vector of parameters that contains components of 8* and 8’; vector of basic parameters Model parameters, i.e., vector of parameters that appear in a model of system. Units depend on model and system Experimental parameters. A vector of parameters not in the model but introduced by experimental procedure. Units depend on system, model, and experiment

X

Y

4

8 e* 6’

SYSTEM,

IDENTIFIABILITY

MODEL,

EXPERIMENT

We do our experiments on systems in the real world but analyze the results in terms of a model or models. Let x be the vector of state variables of a model. For given initial conditions and input to the model there is an equation that describes the time course of changes in the state variables jr= wx, 8”, u, t), x(0) = x() (I) Here 8* is a vector of parameters of the model, the model parameters. The inputs are u(t). An experiment specifies initial conditions, inputs, and the observations, y y = G(x, 0’)

(2 a ) In EQ. Za, G is the observation function that specifies the observations in terms of the state variables and new parameters 8’, which may be introduced in the observations and which we call experimental parameters. We note that experimental parameters can also be introduced in the initial conditions and the inputs. For simplicity we do not introduce such experimental parameters explicitly but note that they are like model parameters in that they influence the observations by way of the solutions for the state variables. Nonetheless, in physiology it is useful to distinguish between parameters inherent to the model and those introduced by the experiment. Thus the basic parameters are the components of a vector 8, which contains 8* and 0’. We should point out that there is some variation in usage for the term model. In some circumstances the term model has been used to include the initial conditions, inputs, any conditions of constraint, and the observations. i.e.. the exneriment (14). Others have kent

OF

PARAMETERS

the idea of experiment separate from that of model of the system (18, 20). In physiology it is more common to use the term model for the model of the system so that we can speak of doing a particular experiment or a number of different experiments on a model, and that is the usage we follow. To be precise one should distinguish a model of a system from a model of an experiment on a system. If one solves Eq. 1 to obtain the state variables as functions of time and the basic parameters, then one can write the observations as functions of time and the parameters, as

Y = w(t,

a 0’1= w, 4)

(2b)

In the form of Eq. Zb, G is often called the response function and is the form of G that is referred to in what follows (See CHECKING LOCAL IDENTIFIABILITY AT A POINT). In Eq. Zb, G is written as a function of time and the basic parameters and also as a function of time and a vector of parameters, 6, which we call observational parameters. By definition, the components of + are identifiable. The basic parameters may or may not be identifiable, but they appear in combinations in the observations that are identifiable, the +i. If the model is a compartmental model with constant fractional transfer coefficients, the equation corresponding to Eq. 1 is q(O) = 90 (3) q = Kq + Bu; Here B is a matrix that gives the inputs as linear combinations of the components of an input vector, u. A number of notations are in use for compartmental systems (17, 20); here we use the notation required by the Modeling Methodology Forum. In Eq. 3, K is the matrix of fractional transfer coefficients; its elements have the units of time? The components of q are the compartment sizes with units of mass, usually in moles. Sometimes it is convenient to recast Eq. 3 in terms of concentrations rather than the mass of the compartments (20). The components of Bu have the same units as 9, mass per’ unit time. A complete summary of the notation used in this paper is given in the GLOSSARY. The model parameters appear in the components of K. If the observations are linear combinations of the compartments, the observations are given by Eq. 4, in which C is the observation matrix

Y = cs

(4) . Experimental parameters could be introduced by way of the initial conditions, qo; the inputs, Bu; and the observation matrix, C. Now we can state the identifiability problem more exactly. For a given model and experiment, are the components of 8 uniquely specified by the components of @? The problem is even simpler if the experiment introduces no parameters: Do the observational -parameters uniquely determine the model parameters? IDENTIFIABILITY

The term identifiability has been used in economics, statistics. and svstems engineering and more recentlv in

Downloaded from www.physiology.org/journal/ajpendo by ${individualUser.givenNames} ${individualUser.surname} (142.103.160.110) on September 25, 2018. Copyright © 1990 American Physiological Society. All rights reserved.

LOCAL

IDENTIFIABILITY

mathematical biology (4, 10, 13,20, 29). There are in use a large number of related terms and concepts, many of which are overlapping. It is not our intention to review what is now an extensive literature. For a general introduction, see the relevant chapters in Refs. 3, 10, 17, 20, 21, and 29; for more details see the papers in Ref. 30. We first define the terms we use. In identifiability one is concerned with the question of theoretical uniqueness of solutions for a given model and experiment. Observational errors play no role; they enter subsequently in the theory of estimation of the identifiable parameters. Hence, in what follows, by observations we mean theoretically exact observations, i.e., without error. The reader should think of observations as the values of the response functions, Eq. Zb, calculated for given values of the 0; for a given experiment. Definitions Parameter identifiability. Suppose we generate observations for a nominal set of values, 8”, of the basic parameters. If a parameter sj is defined by the observations so that there is only one solution, ej”, that parameter is globally identifiable or uniquely identifiable for that model, experiment, and 0’. Another way to state this is that the solutions of the equations 4(e) = @(0”) give only one value in tij and that is 0;. If the observations determine a finite set of values for Bj, then 0j is locally identifiable. The latter includes cases of symmetry in models in which two or more parameters play equivalent roles and so can be interchanged. Note that by this definition local identifiability includes global identifiability as a special subset. Structural identifiability. A property of a parameter is structural if it holds for almost all values of the parameters, i.e., almost everywhere in parameter space. The qualification “almost all values of’ or “almost everywhere” is used to indicate that the property may not hold on a special subset of values, a set of measure zero. Thus structural global or local identifiability are generic properties that are not dependent on particular values of the parameters (29). In this paper we are concerned with local identifiability at a given point or a finite set of points and so we are not examining structural identifiability. Model identifiability. If for 8 = 8’ all of the parameters of a model are globally identifiable, then the model is globally identifiable. If all of the parameters are locally identifiable and at least one is not globally identifiable, then the model is only locally identifiable. If these properties hold for almost all values of 0”, then we have the corresponding structural properties for the model. Conditional identifiability. If for a model and experiment a parameter is not identifiable but setting the values of one or more other parameters makes it identifiable, then the parameter is identifiable conditioned on the values assigned to the parameters that are set. Note that by “setting a parameter” we mean we treat it as a constant of known value, i.e., we remove it from the parameter set. Interval identifiability and quasi identifiability. Berman and Schoenfeld (8) nointed out that even if the data do

OF

E729

PARAMETERS

not uniquely define a parameter, by taking all of the constraints into account such parameters usually have to fall in restricted regions of parameter space. Chau (9) calculated the intervals in which the parameters of a few simple compartmental models must fall when only blood concentration curves are available. That work was done independent of developments in identifiability, and the term was not used. DiStefano (13) has introduced the term interval identifiability to describe the restriction of a nonidentifiable parameter to an interval. It is possible for the interval to be small enough for such a parameter to be “identifiable for practical purposes,” and DiStefano calls this quasi identifiability. Classification

of the Parameters

We classify the basic parameters into three sets. If the values of some of the parameters are fixed, then the classification holds for the remaining parameters but in a conditional sense. Nonsensible or insensible parameters. If a parameter does not influence the observations of an experiment, the sensitivity of the observations to that parameter is zero. In that case the observational parameters are not functions of that parameter and the parameter is said to be nonsensible. Sensible parameters. If the observations change when a parameter changes in value, then that parameter is a sensible parameter; a parameter may be sensible and not be identifiable. In previous papers, one of us (21,22) has used the terms observable and nonobservable for sensible and nonsensible. However, in systems theory the term observable is used to mean measurable or calculable from observations and only in reference to state variables so there is a real possibility of confusion with this usage. Dickinson and Gelinas (11) use the terms sensitive and nonsensitive parameters. That too seems inappropriate, because it is the observations that are sensitive or insensitive to variation in the values of the parameters. We propose the term sensible in the dictionary sense of “capable of being perceived” (3&l),i.e., the observations are sensitive to the sensible parameters and nonsensitive to the nonsensible parameters. LOCALLY IDENTIFIABLE. For the locally identifiable parameters &, . . *,, 6,, the observational parameters & = f&t))have a finite number of solutions for &, . . . , 8,. NONIDENTIFIABLE. For the nonidentifiable parameters, the number of solutions for 8 is infinite. CHECKING

LOCAL

IDENTIFIABILITYATA

POINT

A large part of the literature on identifiability is concerned with structural global and structural local identifiability of models. That is a problem of great interest from the viewpoint of basic theory but also very difficult. A less difficult problem but one that gives considerable useful information is that of testing local identifiability at a given point in parameter space. In physiology we are concerned with testing hypotheses and need to focus on those parameters that have to be estimated to test a hypothesis, whereas the rest need not necessarily be identifiable. Furthermore, we often have auxilliarv infor-

Downloaded from www.physiology.org/journal/ajpendo by ${individualUser.givenNames} ${individualUser.surname} (142.103.160.110) on September 25, 2018. Copyright © 1990 American Physiological Society. All rights reserved.

E730

LOCAL

IDENTIFIABILITY

mation and results from other experiments that provide initial estimates of the parameter values. Even if our estimates are rough, tests of local identifiability at a number of points often provide sufficient information for practical applications. In fact, for compartmental models with constant fractional transfer coefficients, we are assured that if a parameter is locally identifiable at one value, it is locally identifiable for almost all values of the parameter (15, 29), i.e., structurally locally identifiable. This work was started (22) by noting that when the observational parameters are available as functions of the basic parameters, the tests for structural identifiability are done by examining the nature of the solutions of the functional relations, & = fi(B1, . . . , &). We note that the & is generally a nonlinear set of functions of the 0j whether the system model and observations are linear or nonlinear functions of the state variables. For problems of low dimensionality it is easy to generate the +i explicitly as functions of the ej and check structural identifiability on the functional relations, +i = fi(e,, . . . , &). A number of catalogs give structural local identifiability results for linear, constant coefficient models of two and three compartments (12, 27). To test for local identifiability of the sj, all the +i were linearized around the initial estimates of the ej. Local identifiability of a parameter, Sj, in the linearized equations implies local identifiability of the parameter in the nonlinear relations, +i = fi(e,, . . . , &). However, the reverse is not true. It is possible for a parameter to be locally identifiable in the nonlinear equations and to be nonidentifiable in the linearized equations; an example is given in DESCRIPTION OF PROGRAMS, Output. For problems of even moderate magnitude the algebraic work involved in finding the @i becomes limiting. An important finding is that if one has initial estimates of the basic parameters, then one can determine local identifiability numerically at the initial estimates directly without having to generate the observational parameters as explicit functions of the basic parameters. That approach involves linearizing the observations, Eq. Zb, in the ej and uses least squares to estimate the deviations of the basic parameters around the initial estimates (see chapt. 3 in Ref. 29). Again it is true that local identifiability of a parameter in the linearized observations implies local identifiability in the observations if they are nonlinear in the 0j, and again the reverse is not true. Least Squares

We provide the derivation for scalar y; the results extend directly to multiple observations, i.e., a vector y. Let there be p basic parameters and assume we have initial estimates, 0:, . . . , 0;. For the 0j set at the initial estimates, a large number of observations can be calculated, yl, . . . , yn, at n points in time for n > p. Then one can linearize the observations in the parameters Sj around the initial estimates that now play the role of known

OF

PARAMETERS

values p aG’ C TAOj + ei

Yi =G;+

j=l

(5)

j

In Eq. 5, G is the response function, given by Eq. 2b. One can then form the sum of squared deviations (S) over the simulated data set as in S=

i

aG

2

+Adj

yi-G:-C

(6)

i i j ) Then one can find the least-squares estimate of AOj.We let zi = yi - GF. If we let Ai = (A&, . . . , Agp)be the leastsquares estimates of A0 and g be the sensitivity matrix, then i=l

g -

aG0 dG0 - 1 ... - 1 de.l . de.p . . . 0 0 l aG aG -

n

de1

(7)

n

. . . -

deP

Note that zi = yi - GP = 0. The normal equations for the estimates Agj are given by g’giG = g’z = 0 (8) Note that if the determinant of g’g # 0 except possibly for a subspace of lower dimensions of the parameter space, then the model is structurally locally identifiable. Local Identifiability Insensible parameters. One should now examine the columns of g. For the insensible parameters the corresponding columns have only zero entries. Sensible parameters. Suppose we delete the columns of g corresponding to the insensible parameters to form the matrix gs. Then the normal equations corresponding to the sensible parameters are given by

g@ii

= 0

(9) For the identifiable parameters, 0j, the only possible solutions of Eq. 9 are Aij = 0. For the nonidentifiable parameters there should not be such unique solutions. To check that, we use row reduction of gsgs to row echelon form, with pivoting on maximum elements. When completed, the equations should have the following form. 1) There should be rows with only one nonzero entry and that on the diagonal. The column indexes of the diagonal elements of these rows are the indexes for the locally identifiable parameters. 2) There should be rows with nonzero entries on the diagonal and elsewhere. The column indexes of the nonzero elements of these rows are the indexes of nonidentifiable parameters. These parameters are thus shown to be related to one another locally. 3) There should be rows with all zero entries. These are evidence of redundancy in the equations. Correlations

Between Identifiable

Parameters

The final step is to calculate the pairwise correlations for the identifiable parameters. Even if identifiable it

Downloaded from www.physiology.org/journal/ajpendo by ${individualUser.givenNames} ${individualUser.surname} (142.103.160.110) on September 25, 2018. Copyright © 1990 American Physiological Society. All rights reserved.

LOCAL

IDENTIFIABILITY

may be difficult to estimate two parameters separately in the presence of measurement error if they are highly correlated. To check that, we reduce the matrix gsgs by eliminating rows and columns corresponding to the nonidentifiable parameters to obtain gig,. Because gigi is the matrix corresponding to the identifiable parameters, it is nonsingular, and so we invert it to obtain (gig&‘. The correlation matrix is obtained by dividing the i, j element of (gig&’ by the square root of the product of the ith and jth diagonal elements. Note that the values of the correlation coefficients depend on the parameter values and on the choice of points in time where the sensitivities are calculated. IDENTIFIABILITY PRACTICAL

AND ESTIMABILITY: CONSIDERATIONS

An important step in the development of our ideas was the recognition that the two following problems are separable: Are the observations sufficient to specify the parameters? and What are the effects of the errors of measurement on the estimates of the parameters?. The first is the identifiability problem; the second is the estimation problem. In theory one tests identifiability first and then proceeds to estimate the identifiable parameters. However, for practical purposes there are areas of overlap. Here we briefly discuss some of the problems so as to emphasize that the ideas and all of the programs available for identifiability and for estimation cannot be applied in a rote manner. Nonidentifiable

but Estimable for Practical Purposes

Even though a parameter may be nonidentifiable, the constraints in most problems restrict the physically permissible solution values to a subspace, an interval in each of the nonidentifiable parameters (13). If the interval is small, then that may be close enough for many practical purposes. So if a parameter turns out to be nonidentifiable, it may be worthwhile to check interval identifiability (13). Identifiable

but Poorly Estimable

A model or a set of parameters may be identifiable but the relevant g’g matrix may be near singular, i.e., the determinant of g’g is close to zero, resulting in poor estimability. That can occur in a number of ways. It can occur as the parameter vector approaches a subspace (a space of lower dimension) on which some parameters are not identifiable. An example is given in DESCRIPTION OF PROGRAMS. It can also hold in a region of the parameter space (not just a space of lower dimensions) where there are large differences in values of some of the parameters so that det(g’g) is near zero throughout the region. In such a region the errors in the observations will have such a large effect on the estimates that the estimates have high variances. Increasing the sample size and optimizing the sample set gives only small improvements. One must recognize the problem by checking the det(g’g) and seek other experiments and/or other information to constrain the estimation. Another wav for this to occur

OF

E731

PARAMETERS

is if one chooses a set of sample points that is far from optimal. That can be a difficult problem with real data sets. For checking identifiability, we have evaded it so far by choosing many points over the time period that covers the major range of variation of the observations; however, we expect to address that issue in the next development in IDENT. ESTIMATION

OF NONIDENTIFIABLE

PARAMETERS

A large number of parameter-fitting programs are used routinely, and none of the standard programs check identifiability. What happens if one tries to estimate a nonidentifiable parameter? We will call the subspace of solutions for nonidentifiable parameters the subspace of feasible solutions. What a particular fitting routine does depends on the search technique used. Derivative methods, such as the unmodified Newton-Raphson and Gauss-Newton methods (2, 3), calculate the Hessian of the criterion function (Newton-Raphson) or the equivalent of our g’g (Gauss-Newton) to take a step in the search. Both of these blow up if there are unidentifiable parame parameters, ters, because the Hessian and g’g are singular. However, most modern derivative methods, such as the Marquardt method (2,3), modify the Gauss-Newton g’g matrix to make it nonsingular, and such modified methods walk into the feasible subspace and stop when their termination test is satisfied. The same is true of derivative-free methods, such as the direct-search methods, that involve calculation of the values of the criterion function, e.g., sum of squares, directly in the search. Where such routines end up in the feasible subspace depends on the starting point and on the details of the search method. In what follows we refer specifically to the compartmental system-fitting software, SAAM (7, 16), that uses a modified Gauss-Newton method. A finding that should make one suspicious is a small overall sum of squares, indicating one is getting a good fit to the observations, and coupled with that, a very large standard deviation of the estimates of some parameters. That finding should make you suspect some of the parameters are unidentifiable or identifiable but very poorly estimable. The following simple examples illustrate what happens. Example 1

First we will consider the compartmental system, the connection diagram of which is shown in Fig Fig. , 1. We start with nothing in the system, inject 100 100 units into compartment 1 at t = 0, and then observe compartment 1. k21

1

k 01

1. Model partment 1, and nonidentifiable. FIG.

for example compartment

1. Input is an impulsive 1 is observed. Both

input into comb1 and kzl are

Downloaded from www.physiology.org/journal/ajpendo by ${individualUser.givenNames} ${individualUser.surname} (142.103.160.110) on September 25, 2018. Copyright © 1990 American Physiological Society. All rights reserved.

E732

LOCAL

IDENTIFIABILITY

OF

PARAMETERS

Example 2

FIG.

starting points.

2. Subspace of feasibl .e solutions for example points to estimates given by SAAM for

FIG. 3. in Fig. 2. fractional observations. choices of

1 Arrows 3 different

go from initial

3 kO1 Estimates given by SAAM for same 3 initial points shown In this case errors taken from a normal distribution with a standard deviation equal to 0.3162 were added to calculated For each initial point, estimates are shown for 3 different error sets.

The second example is the same as example 1 with the addition of an excretion from compartment 2 so that 12*2 is an insensible parameter (Fig. 4). The same data sets were run. For zero added error the search again goes to different points in the feasible subspace of hl + kZ1, but it makes no change in the starting value of the insensible parameter. In this case the insensible parameter had a very large FSD, but one of the other parameters had fairly small FSD. These examples illustrate some of the clues one can get from a SAAM output. One uses initial estimates of the parameters to generate a data set of observations on the model without added error; then one turns around and tries to fit the observations using a different starting point in the parameters. If during a fit, a parameter does not change in value and its FSD is high, then one should repeat the run starting at a number of different points. If the same thing happens, the parameter is likely to be insensible. If one starts at a number of different initial points and the estimation routine ends up in different final points, then a number of possibilities need to be explored. 1) The points are part of a feasible subspace for unidentifiable parameters. This is usually indicated by a good fit to the data, but the parameters have high FSDs at each of the final points. 2) The final points are different points of local identifiability. All will have the same sum of squares, zero, and the FSDs for the parameters will be very small. 3) The search got hung up in local minima in the sum of squares surface. The sum of squares will generally differ at the different final points. Further exploration may bring out a pattern of local minima. If one runs into such problems, a check of local identifiability at a number of points will often clarify the situation. DESCRIPTION

for example 2. Input is an impulsive partment 1, and compartment I is observed. Parameters nonidentifiable; koz is insensible. FIG.

4. Model

input into comkol and kzl are

For this experiment only, kol + k,, is identifiable. If koI + kzl = c, then that is the subspace of feasible solutions in parameter space for positive values of the parameters (Fig. 2). We show what the parameter estimation routine SAAM (7,16) does when presented with this problem. It walks into the feasible subspace and stops. Using c = 3, we generated data at 10 sample points with and without added error: in case 1 no error was added, and in case 2 normal deviates with a fractional standard deviation (FSD) of 0.3162 were added. Figure 2 shows where the dG -= de i

OF

PROGRAMS

Our first program for the method given in CHECKING IDENTIFIABILITY AT A POINT, IDENTl, was written in FORTRAN 66 and run on the University of Michigan computer center IBM 3090-400. In the next version, IDENT?, a front end was added that prompts the user through the input phase. In all versions of IDENT, the calculations are in double-precision arithmetic. In IDENTZ, to obtain the required accuracy in the calculation of the partial derivatives of the observations with respect to the parameters the higher order symmetrical central difference formula of Eq. 10 is used with h; = O.OOl&. The error in this is of order h4 LOCAL

8G(Oi + hi) - 8G(Bi - hi) - G(6i + Zhi) + G(Oi - Zhi) 12hi

search ended for three different starting points for no added error; for case 2 the search ended in points off but near the feasible subspace (Fig. 3). The exact point where the search ends now also depends on the pattern of errors in a particular data set.

This worked very well, but in the next version, IDENT3, we introduced a method that has been reported to be a better, direct integration of the differential equations for the partial derivatives (11, 26). It gives an improvement in accuracy of the calculated partials of less than one

Downloaded from www.physiology.org/journal/ajpendo by ${individualUser.givenNames} ${individualUser.surname} (142.103.160.110) on September 25, 2018. Copyright © 1990 American Physiological Society. All rights reserved.

LOCAL

IDENTIFIABILITY

significant figure, which is not really enough to compensate for the increased work of setting up even moderately sized problems. The user provides the differential equations for the state variables as an external function call. The experiments are assumed to be impulsive inputs into one or more compartments of a zero state model; the impulsive inputs are specified by the initial conditions. For linear models, impulsive inputs are sufficient, but for nonlinear models, more information may be gained from continuous inputs. The following is a brief description of our latest version IDENT2C, written in C. Inputs

The user provides values for the following, either in the subroutines or from program prompts. Global variables and subroutines. The following are requested: 1) num-states, the number of state variables (compartments); 2) num-pars, the number of basic parameters; 3) num-obs, the number of observations (for example, the number of rows in the C matrix); 4) modelname, a string containing a name for the model; 5) derivs, a C subroutine that contains the differential equations for the state variables; and 6) observs, a C subroutine that contains the equations for the observations. Program prompts. The following are requested: 1) name for the experiment (used to recall the experiment at a later date); 2) name for the output file; 3) initial conditions for each state variable; 4) initial estimates for each parameter; 5) number of sample points, in time; 6) initial time point; 7) time points at which the observations are made; 8) option to print the partials matrix (partials of the observations with respect to the parameters); 9) option to print the observations at each time point; 10) option to change any of the above program inputs; and 11) option to continue with the above experiment. We have used 50-200 points spread from t = 0 out to a time when all observations have fallen to ~5% of their maximum values. We are working on incorporating optimal sampling theory in our next version, IDENT4, so as to choose a minimum size, D-optimal (24) point set. The decrease in number of points at which sensitivities are calculated should increase the accuracy of the computations, and the optimal placement of the points should reduce the correlations between the identifiable parameters. Program and Subroutines

The program is made up of a number of modules. 1) getq: this solves the model equations for the values of the state variables at the time points. 2) dfdq: this solves the model differential equations repeatedly and calculates the partials matrix, i.e., the sensitivity matrix. Both getq and dfdq call a differential equation-integrating routine to integrate the user-provided derivs. 3) getgpg: this reads in the observation matrix (C matrix) and calculates the g’g matrix. The g’g matrix is scaled by dividing each row by its largest element to give a stand-

OF

PARAMETERS

E733

ardized g’g matrix. 4) rowred: this row reduces the standardized g’g matrix and classifies the parameters as insensible, nonidentifiable, or locally identifiable. 5) This subroutine calls subma, which sets up the identifiable submatrix of g’g and calculates the correlation matrix for the locally identifiable parameters. output

Several pages of output are generated in the file named as the output file. In addition, the main identifiability results and correlation matrix are printed on the screen. The data contained in the output file is organized as follows: 1) the experiment is displayed, including initial values and the sample time points; 2) the values of the observations at each time point [if you respond Y(es) to the prompt “. . . Observations printed?“]; 3) values of the partial derivatives for each observation with respect to the parameters [if you respond Y(es) to the prompt “. . . Partials printed?“]; 4) the g’g matrix and the standardized g’g matrix; 5) the row-reduced g’g matrix, the parameter classification, and the gf g, matrix; and 6) the correlation matrix for the identifiable parameters. Using the Program

The basic theory assumes one can do error-free calculations. For models with many parameters the many arithmetic operations involved lead to an accumulation of round-off and truncation errors in the final rowreduced matrix. For this reason we carry out the calculations in double-precision arithmetic and then truncate the final row-reduced matrix to eight significant figures. For larger problems it may be necessary to go to higher precision operations. Some judgement is needed in interpreting results when the final row-reduced matrix ends up with a very small number in the diagonal position of a row and zeros in the other positions in the row. One must always wonder whether the parameter involved (indicated by the column index of the diagonal position) is really identifiable or whether one is seeing the effect of error accumulation. One thing to do is to run the program for a series of nearby parameter values to see if there are any trends in the results. For a set of parameters that are in general identifiable, there may be particular combinations (a subspace) on which the parameters are nonidentifiable. Small entries appear in some of the diagonal elements of the rowreduced matrix as one approaches such a subspace; these rows become zero rows when the subspace is reached, indicating the equations have become redundant. An example will illustrate this. We now consider again the simple model shown in Fig. 4, for the experiment in which the model starts at zero initial conditions, and an impulse is injected into compartment 1 and compartment 2 is observed. For this experiment & is structurally locally identifiable, as are kol and ho2 in the nonlinear equations that relate the 6 to the 4. However, the program checks the linearized observations and then kol and kO, are not identifiable if they fall in the subspace bo2 = kO, + k,,. For rFzZ1 = 1.5, Kol = 1.25, and ko2 = 2.75, the

Downloaded from www.physiology.org/journal/ajpendo by ${individualUser.givenNames} ${individualUser.surname} (142.103.160.110) on September 25, 2018. Copyright © 1990 American Physiological Society. All rights reserved.

E734

LOCAL

IDENTIFIABILITY

row reduction gives the expected results: rlzZ1is locally identifiable and rFzoland ko2 are not. To describe the results, the columns corresponding to &, &, and & are 1, 2, and 3, respectively. In row 1 only position 1,l is not zero ( jlzZ1is locally identifiable). Row 2 has a 1.0 in positions 2,2 and 2,3. Row 3 is all zeros. As one moves off of the subspace & = rlzol+ &, the 1.0in position 2,3 is replaced by zero, and a small entry appears in position 3,3. For example, for ko2 = 2.76 or 2.74, the entry in position 3,3 becomes 7.3 x 10N7,and the row-reduced matrix is then diagonal. As one moves away from the condition ko2 = kol + hZ1, that small diagonal element increases rapidly. Thus if suspicious results are obtained, one should explore a series of parameter values going away from that point. Constraints on Parameters

Equality constraints are easily handled. For that one uses the equations of constraint to reduce the parameter set. To do that, one solves for parameters that are of least interest and substitutes the expressions for those parameters in the original differential equations. Inequality constraints are handled differently. For that one should check identifiability for a number of sets of values of the parameters, making sure that the values chosen cover the range of physiologically reasonable values and satisfy the inequality constraints. EXAMPLES

The program has been checked against a large number of two- and three-compartment models with single and multiple inputs and single and multiple observations as well as with a few large systems (10 compartments) that have connectivities and/or constraints such that we

FIG.

5. Two-compartment

model

for

distribution

nephrine.

FIG.

6. Model

for distribution

of tracer

norepinephrine.

of tracee

norepi-

OF

PARAMETERS

could solve for the observational parameters, $i, explicitly in terms of the & and determine structural local identifiability. We present two examples. One is a two-compartment model that has been used in the analysis of the distribution of norepinephrine in humans (25). The other is a nonlinear, nonautonomous system of equations used to model the intravenous glucose tolerance test (5, 6, 28). Two-Compartment

Model of Norepinephrine

Kinetics

Under constant conditions plasma norepinephrine levels are stable at a balance between production, metabolism, and reuptake at sympathetic nerve terminals. If tritiated norepinephrine is infused to constant plasma levels and the infusion is stopped, then the time course of the plasma tritiated norepinephrine is closely fitted by a two-exponential decay. Thus the minimal model for the unlabeled norepinephrine must be a two-compartment model, such as that shown in Fig. 5. In that model, compartment 1 is the plasma norepinephrine plus any that equilibrates rapidly with plasma, such as that in extracellular space of liver or bound in vessels. Compartment 2 is the norepinephrine in other, more slowly turning over extracellular spaces; the input into compartment 2 is primarily from sympathetic nerve terminals. In Fig. 5, klz is the fraction of compartment 2 transferred to compartment 1 per unit time, and & is the fraction removed by uptake into tissues and by metabolic change. Although &, &, &, and klz may be functions of concentrations of norepinephrine, in the steady state they are constant. For steady-state conditions, if labeled material is introduced it must distribute according to linear kinetics (20). The compartment model for the distribution of tracer is Fig. 6. Note that V1 (in liters-‘), the volume of distribution for compartment 1, can be estimated from the initial rate of decrease of plasma tritiated norepinephrine concentration after the infusion is stopped. For input into compartment 1 and observation of compartment 1, none of the parameters are identifiable. However, if we have independent knowledge of one of the parameters or an equation of constraint between the parameters, they become identifiable. In fact, estimates of urinary excretion that gave estimates of kol and an independent constraint relating Ixzo2 and klz gave similar results (25). Local identifiability was checked with IDENT for a number of values of the parameters in the range of values found (25). At all points IDENT gave the correct results. For another example of the identifiability problem in the distribution of a radioactive compound but involving a three-compartment model, see Jennrich and Bright (23) and the subsequent analysis by Griffiths (19). A Model of the Glucose-Insulin to Intravenous Glucose

Response

A model for the intravenous 6, 28) can be written as

glucose tolerance test (5,

kl = -(e, + x2)x1 + ~1x10

(IO

k2 = -&x2 + 63x3

(12)

23 =

-&lx3

+

e*tbl

-

x10)

(13)

Downloaded from www.physiology.org/journal/ajpendo by ${individualUser.givenNames} ${individualUser.surname} (142.103.160.110) on September 25, 2018. Copyright © 1990 American Physiological Society. All rights reserved.

LOCAL

IDENTIFIABILITY

Here x1 is the plasma glucose concentration, with xl0 being the steady ,-state fasting vvalue. Both x2 and x3 represent insulin deviati .ons from steady-state va lues; x3 1s the plasma insulin concentration, but x2 is “active” in .sulin in a remote (tissue) compartment (6). Of the state variables, x1 and x3 are potentially measureable, whereas x2 is not. The time dependence of Eq. 13 is an unusual feature that we consider to be unrealistic, but it was the combination of a nonlinear set of differential equations with time dependence that led us to choose this as an example set. In the base-line steady state, x1 = xlo, x2 = 0, and x3 = 0. The glucose tolerance test is then an impulsive input into x1 at t = 0. What parameters are identifiable if x1 and x3 are observed or if x1 alone or x3 alone are observed? We chose xl0 = 90 mg/dl for the steady-state glucose concentration and o1 = 0.045, ti2 = 0.025, e3 = 0.00018, & = 0.002, and & = 0.24, all of which are values falling in the ranges reported in Bergman et al. (6). If x1 and x3 are observed, then all parameters are locally identifiable. If only x1 is observed, then only &, ti2, and & are locally identifiable, whereas if only x3 is observed, then all parameters are again locally identifiable. We have repeated these runs but with treating xl0 as a parameter, &. All parameters are again locally identifiable from measurement of x1 and xi. If x1 is observed, only 81, 02, &, and 06 are locally identifiable, and if only x3 is observed, only t92and 0; are locally identifiable. The results for x1 and x3 observed are in accord with experience on estimating the parameters (3) CONCLUSION

It would be ideal to have simple methods for checking structural identifiability of parameters. However, that is extremely difficult to do for large models with many parameters and particularly so for nonlinear systems. We have developed a simple method to test local identifiability at a point in parameter space that works for linear and nonlinear systems, and we now have a number of programs to do that. Although this is not as powerful as testing structural identifiability, it is a very useful step in checking local identifiability. The following programs are currently available from the authors on request. 1) IDENT2: the source code is written in FORTRAN. The program is available in two versions, one for the VAX and the other for the ATT Unix PC. 2) IDENT2C: an updated and modified version of IDENT with the source code in C. This is our current operational version, run on a SUN SPARCstation 1. 3) IDENT3: essentially the same as IDENT but with the numerical calculation of partials replaced by direct integration of the differential equations for the sensitivities. IDENT is also available in versions for the VAX and the ATT Unix PC. We thank the referees of the prior versions of this paper for their constructive and very useful critiques. This work was supported in part by Division of Research Resources Grant RR-02176.

E735

OF PARAMETERS

Address for reprint requests: J. A. Jacquez, Dept. of Physiology, Univ. of Michigan, 7808 Medical Science II, Ann Arbor, MI 481090622. Received 23 November 1987; accepted in final form 8 December 1989. REFERENCES D. H. Lecture

1. ANDERSON, Modeling

and

Tracer

Kinetics.

Notes

in Biomathetics.

Compartment

New York: Springer-Verlag,

1983,

vol. 50. BARD,

Y. Nonlinear

Parameter

Estimation.

New York: Academic,

1974.

6. 7.

BATES, D. M., AND D. G. WATTS. Nonlinear Regression Analysis and Its Applications. New York: Wiley, 1988. BELLMAN, R., AND K. J. ASTROM. On structural identifiability. Math. Biosci. 7: 329-339, 1970. BERGMAN, R. N., Y. Z. IDER, C. R. BOWDEN, AND C. COBELLI. Quantitative estimation of insulin sensitivity. Am. J. Physiol. 236 (Endocrinol. Metab. Gastrointest. Physiol. 5): E667-E677, 1979. BERGMAN, R. N., L. S. PHILLIPS, AND C. COBELLI. Physiologic evaluation of factors controlling glucose tolerance in man. J. Clin. Invest. 68: 1456-1467, 1981. BERMAN, M., W. F. BELTZ, P. C. GREIF, R. CHABAY, AND R. C. BOSTON. CONSAM Users Guide. Bethesda, MD: Natl. Inst. Health

Lab. Math. Biology, 1983. 8. BERMAN, M., AND R. L. SCHOENFELD. Invariants in experimental data on linear kinetics and the formulation of models. J. Appl. Phys. 27: 1361-1370,1956. 9. CHAU, N. P. Linear pharmacokinetic models: geometric construction to determine transfer and elimination rate constants. J. Pharmacokinet. Biopharm. 5: 147-159, 1977. 10. COBELLI, C., AND J. J. DISTEFANO III. Parameter and structural identifiability concepts and ambiguities: a critical review and analysis. Am. J. Physiol. 239 (Regulatory Integrative Comp. Physiol. 8): R7-R24,1980. 11. DICKINSON, R. P., AND R. J. GELINAS. Sensitivity analysis of ordinary differential equation systems-a direct method. J. Comput. Phys. 21: 123-143,1976. 12. DISTEFANO, J. J., III. Design and optimization of tracer experiments in physiology and medicine. Federation Proc. 39: 84-90, 1980. 13. DISTEFANO, J. J., III. Complete parameter bounds and quasiidentifiability conditions for a class of unidentifiable linear systems. Math. Biosci. 65: 51-68, 1983. 14. DISTEFANO, J. J., III, AND C. COBELLI. On parameter and structural identifiability: non-unique observability/reconstructability for identifiable systems, other ambiguities, and new definitions. IEEE

Trans.

Autom.

Control

25: 830-833,198O.

15. EISENFELD, J. A simple solution to the compartmental structural identifiability problem. Math. Biosci. 79: 209-220, 1986. 16. FOSTER, D. M., AND R. C. BOSTON. The use of computers in compartmental analysis: the SAAM and CONSAM programs. In: Compartmental Distribution of Radiotracers, edited by J. R. Robertson. Boca Raton, FL: CRC, 1983. 17. GODFREY, K. Compartmental Models and Their Applications. London: Academic, 1983. 18. GREWAL, M. S., AND K. GLOVER. Identifiability of linear and nonlinear dynamical systems. IEEE Trans. Autom. Control 21: 833-837,1976. D. Structural identifiability for compartmental models. 19. GRIFFITHS, Technometrics 21: 257-259, 1979. Analysis in Biology and Medicine 20. JACQUEZ, J. A. Compartmental (2nd ed.). Ann Arbor: Univ. of Michigan Press, 1985. 21. JACQUEZ, J. A. Identifiability: the first step in parameter estimation. Federation Proc. 46: 2477-2480, 1987. 22. JACQUEZ, J. A., AND P. GREIF. Numerical parameter identifiability and estimability: integrating identifiability, estimability and optimal sampling design. Math. Biosci. 77: 201-227, 1985. 23. JENNRICH, R. I., AND P. B. BRIGHT. Fitting systems of linear differential equations using computer generated exact derivatives. Technometrics 18: 385-399, 1976. 24. LANDAW, E. M. Optimal multicompartmental sampling designs for parameter estimation: practical aspects of the identification problem. Math. Comput. Simul. 24: 525-530, 1982. 25. LINARES, 0. A., J. A. JACQUEZ, L. A. ZECH, M. J. SMITH, J. A. SANFIELD, L. A. MORROW, S. G. ROSEN, AND J. B. HALTER.

Downloaded from www.physiology.org/journal/ajpendo by ${individualUser.givenNames} ${individualUser.surname} (142.103.160.110) on September 25, 2018. Copyright © 1990 American Physiological Society. All rights reserved.

E736

LOCAL

IDENTIFIABILITY

Norepinephrine metabolism in humans. Kinetic analysis and model. J. Clin. Inuest. 80: 1332-1341, 1987. 26. MARKUSSEN, C. R., AND J. J. DISTEFANO III. Evaluation of four methods for computing parameter sensitivities in dynamic system models. Math. Biosci. 61: 135-148, 1982. 27. NORTON, J. P. An investigation of the sources of nonuniqueness in deterministic identifiability. Math. Biosci. 60: 89-108, 1982. 28. TOFFOLO, G., R. N. BERGMAN, D. T. FINEGOOD, C. R. BOWDEN,

OF PARAMETERS AND C. COBELLI. Quantitative estimation of beta cell sensitivity to glucose in the intact organism. Diabetes 29: 979-990,198O. 29. WALTER, E. Lecture Notes in Biomathematics. Identifiability of State Space Models. New York: Springer-Verlag, 1982, vol. 46. 30. WALTER, E. (Editor). Identifiability of Parametric Models. Oxford, UK: Pergamon, 1987. 31. WEBSTER’S NEW COLLEGIATE DICTIONARY. Springfield, MA: Merriam, 1961.

Downloaded from www.physiology.org/journal/ajpendo by ${individualUser.givenNames} ${individualUser.surname} (142.103.160.110) on September 25, 2018. Copyright © 1990 American Physiological Society. All rights reserved.