INTRODUCTION

TOP ABSTRACT INTRODUCTION PROTOCOL AND DATA FORMALISM AND NUMERICAL METHODS PHYSIOLOGICAL BACKGROUND MODEL BUILDING DISCUSSION APPENDIX A APPENDIX B APPENDIX C REFERENCES

DESPITE THE NUMEROUS NONLINEARITIES involved in metabolic processes per se and in their regulation (e.g., hormonal), most compartmental and/or noncompartmental models of in vivo mineral metabolism are linear. On the basis of tracer kinetic data and with the assumption that the system is maintained in a constant steady state, at least over the experimental duration, these models have been widely used to obtain quantitative estimates of processes such as intestinal absorption and urinary excretion. These models have also been used to estimate internal distribution pools (IDP) and other transfer pathways of many minerals; they may be macronutrients, such as Ca (1) and Mg (2), or trace elements, such as Zn (11) and Se (28). In some respects, compartmental and noncompartmental models of in vivo radioactive Sr kinetics (9, 20, 30) differ from other investigations, in that the trace element Sr has been mainly used as a "tracer" for Ca in comparative metabolic studies.

Over the past few years, interest has increased in the mechanisms underlying the regulation of trace element metabolism in healthy subjects or animals given high, but nontoxic, doses of the mineral (18, 39, 44, 45) or elements that interact with it (23, 31, 35). Linear compartmental analysis has been used to identify the sites of regulation by comparing the parameters of a model built from initial steady-state tracer kinetic data with those obtained under other steady state(s) after changes in the mineral intake, i.e., long after the suspected regulatory mechanisms have been in full nonlinear operation.

Mathematical modeling may help describe the physiological characteristics of regulatory and/or adaptive processes of mineral metabolism, provided the study is carried out in a non-steady state.¹ Under these conditions, the nonlinear expressions result in variations that are observed under physiological conditions (37) or experimentally induced by displacing the metabolism or its related controlling system from the initial (steady) state. In any case, the compartmental formalism, with or without tracer data (8), must explicitly incorporate the nonlinearities postulated to be appropriate. Typical recent applications of such nonlinear models include the study of endocrine-metabolic systems in non-steady state (5, 24).

It is difficult to obtain information about the nonlinear processes that help regulate Ca metabolism in vivo under physiological conditions because of the powerful self-regulatory system governing extracellular Ca homeostasis (27, 38). The situation is quite different for Sr, a mineral trace element very similar to Ca in its physicochemical properties. A fairly high oral intake of Sr is well tolerated and results in a large increase in extracellular Sr concentration, without any apparent hormonal response.

The purpose of our study is to use compartmental modeling to analyze the kinetics of Sr metabolism that is greatly shifted from its physiological state. We used experimental plasma and urinary data from four groups of ~10 postmenopausal women given four oral doses of Sr (S-12911, molecule containing 2 atoms of nonradioactive Sr; Institut de Recherches Internationales Servier). The kinetic data cover the increase in plasma Sr concentration over the 25-day period for which Sr was given [administration period (AdP)] and its decrease over the 28 days after Sr administration [postadministration period (PAdP)]. This report describes the development of a nonlinear seven-compartment model of human Sr metabolism and the underlying assumptions. The building strategy led to the introduction of two kinds of nonlinearity. The first accounts for the saturable processes [Michaelis-Menten (M-M) and Langmuir-type equations] that behave as time-implicit nonlinear functions and are intrinsic to Sr metabolism, because they form part of the system of differential equations for Sr metabolism. The second is extrinsic to Sr metabolism and requires a specific time-explicit differential formulation that defines the kinetic behavior of biological component(s) other than Sr. Surprisingly, these last model properties may reflect self-regulatory processes specifically operating on Sr metabolism. Because there has been no evidence that Sr has any physiological role, an attractive alternative is that our modeling procedure has revealed some difficult-to-measure regulatory processes linked to Ca metabolism because of the analogies between Sr and Ca. At the stage of the present study, the building procedure is reported with the primary goal to show how pertinent information relative to nonlinear expressions can be extracted from the data. Therefore, only a minimal interpretation of the model with regard to Sr metabolism per se is made.

	PROTOCOL AND DATA

TOP ABSTRACT INTRODUCTION PROTOCOL AND DATA FORMALISM AND NUMERICAL METHODS PHYSIOLOGICAL BACKGROUND MODEL BUILDING DISCUSSION APPENDIX A APPENDIX B APPENDIX C REFERENCES

Total Sr plasma and urine concentrations were measured by inductively coupled plasma emission spectrophotometry. All data were obtained during a study of the acceptability of repeated oral doses of S-12911. Four groups, each composed of 9-10 healthy postmenopausal women, were followed for 54 days. All subjects gave informed consent. The women were 59.5 ± 4.9 yr old and of normal height (158 ± 5 cm) and weight (62 ± 10 kg). They consumed a normal diet (actual daily dietary Ca and Sr intake data are not available). They were given 0.5, 1, 2, or 4 g of S-12911 for 25 days (AdP) in two daily doses: the first was given 1 h after breakfast and the other 1 h before dinner, i.e., with intervals of 10 and 14 h between doses, respectively. Plasma samples were collected just before administration (time 0) and on days 1, 2, 7, 14, 21, 25, and 26-54, with detailed kinetics after the morning daily administration on days 1, 14, and 25 and with sampling frequency of 1-3 days during the 28 days after withdrawal of repeated administration (PAdP; Fig. 1). The number of measurement times over the experimental duration was 38 (26 and 12 during AdP and PAdP, respectively) per individual, with a total of ~350 data points, for each dose.

View larger version (33K): [in this window] [in a new window]   Fig. 1.   Kinetic data (means ± SE) used for model building and typical y1 response of an optimal model (solid lines in B and C). A: mean plasma Sr concentrations over duration of experiment [administration (AdP) and postadministration (PAdP) periods] for the 4 doses: D1 (), D2 (), D3 (), and D4 (). For convenience, only maximal and minimal daily values are reported during AdP. B and C: expanded time scales cover detailed kinetic data during AdP: the first 10 h for D2 (B) and experimentally monitored part of day 14 for D3 (C). Goodness of fit can be appreciated by comparing behavior of model 12 of compartment 1 with these experimental data.

Twelve-hour urine collections were obtained after the end of treatment, on days 25, 28-32, 35, 46, and 53, and their volumes were estimated. Other data on total and ionized plasma Ca and urinary Ca and creatinine were obtained on days 1, 15, 27, and 54.

	FORMALISM AND NUMERICAL METHODS

TOP ABSTRACT INTRODUCTION PROTOCOL AND DATA FORMALISM AND NUMERICAL METHODS PHYSIOLOGICAL BACKGROUND MODEL BUILDING DISCUSSION APPENDIX A APPENDIX B APPENDIX C REFERENCES

Mathematical Formalism

(1)

(2)

The y_i associated with plasma is y₁. The initial conditions, y_i(0) = y_i(t = 0) and z_i(0) = z_i(t = 0), are calculated with the assumption that systems 1 and 2 are in steady state (dy_i/dt = 0 and dz_i/dt = 0) at time 0 with F(t) = F(t = 0), i.e., before oral Sr has any effect, displacing Sr metabolism from its initial physiological state.

Remark 1. All the experimental plasma data are concentrations. Hence, except for y₁ (compartment 1), any y_i and parameters having a concentration dimension must be used with care. Biological interpretation of system 1 requires an estimate of the apparent distribution volume (V_app, assumed to be constant over the experimental period). This enables concentration (y₁) and apparent concentration (C_app for other y_i) to be transformed to a mass dimension. This conversion was done by using urinary excretion data with their 12-h integral values calculated from the Sr concentration and volume (see Eq. 3). It is also possible that V_app does not apply to system 2. Under these conditions and in the absence of direct experimental data for system 2, variations in z have no meaning as absolute values; only z kinetic behavior is of interest.

Remark 2. When K, H, or input G is dependent on y and/or z variables, their dependence on system 1 and/or system 2 must be specified; e.g., K_ij = f(k,y,z) with l = 1, 2, 3, ... and k having a constant value 0. Otherwise, they are strictly linear, have constant values 0 (e.g., K_ij = k_ij), and are referred to as the input rate (g) and the fractional transfer coefficient (FTC) called k or h. For any FTF in system 1, the term intrinsic nonlinearity (INL) or extrinsic nonlinearity (ENL) has been used to specify its dependence on its own (intrinsic) variables (y) or on the (extrinsic) variables associated with system 2 (z).

Remark 3. The general form (systems 1 and 2) defined above was considered only if the kinetics of some z variables do not satisfy the assumption of a pseudo-steady-state assumption: indeed, if the kinetics of a given z are very rapid, at each time, dz/dt  0, and z approximates its asymptotic value, [e.g., as assumed for formulation of the M-M equation (32)]. Then the dependence of system 1 on z [through some K = f(z)] can be directly formulated as a time-implicit nonlinear function of y, i.e., equivalent to an INL, as opposed to its dependence on the time-explicit differential form of z, with its corresponding ENL. Otherwise, when the kinetics of y and z are not very different, dz/dt  0 must be considered to take into account the effects of z changes on y kinetics, and the general form of both differential systems must be used.

Remark 4. For a given set of parameter values together with the value of the initial Sr plasma concentration, y₁(t = 0)  0, an estimate of the other initial conditions y_i(0) and the initial values of input functions F_i0(t = 0) are required. If Sr metabolism is approximated to be physiologically in a constant steady state, an analytic calculation is possible for some simple situations (a single positive and real solution exists for dy_i/dt = dz_i/dt = 0). When multiple input functions operate within a given model structure, the relations specifying the interdependence between the inputs and the irreversible outputs must be defined a priori.

Remark 5. The compartmental formalism presented here is similar to that successfully used to model endocrine-metabolic systems such as glucose metabolism and its hormonal control (4). It also has certain similarities to the conceptual framework developed for modeling non-steady-state pharmacokinetic/pharmacodynamic data (33).

Numerical Methods

The fit of the model response to the experimental data was carried out by nonlinear parameter estimation with y₁ (compartment 1) compared with plasma Sr kinetics. The method of Levenberg-Marquardt (16, 22) was used to minimize the weighted residual sum of squares (WRSS) as criterion function (8). WRSS was computed using the mean plasma Sr concentration and the reciprocal of the measured error variance at each sampling time and for all the doses studied, except in the preliminary study of stage 1 of model building, during which each dose of Sr was considered separately. Because the optimization method is described as a local one, the initial set of parameter values was chosen manually for its ability to describe the overall behavior fairly satisfactorily. The minimization procedure was repeated using a randomly noised initial set of parameter values to avoid obtaining a local and/or false minimum WRSS. This was done until the values of WRSS and the identified parameter values did not differ between several successive iterations.

To examine the goodness of fit to experimental data, we used, in addition to visual inspection, the run test to check the independence of weighted residuals. We also detected misfitting by observing the values of the partial WRSS computed for each set of data, because the model predictions are usually related to four sets of data (1 per Sr dose). As an indication of validity for increasing complexity in the model structure, the significance of the improvement of fit was checked with an F test.

The accuracy of the identified parameter values (practical or a posteriori identifiability) was estimated from the computed variance-covariance matrix, and precision was expressed as coefficient of variation (CV) in terms of percent fractional standard deviation (7). Theoretical (a priori) identifiability was not analyzed because of the high degree of complexity reached by the nonlinear model structures. Thus the unique estimate of the set of parameter values cannot be completely ensured, despite use of various initial sets of parameter values for optimization. Moreover, when the tested model structure was nonidentifiable, some parameter values diverged during optimization and/or there was a very large (much more than 100%) computed CV. Here, a model is ruled out only if one or several of its parameter values has a CV  100%. A correct accuracy is assumed if the calculated CV is such that the parameter differs significantly from zero for P < 0.05.

	PHYSIOLOGICAL BACKGROUND

TOP ABSTRACT INTRODUCTION PROTOCOL AND DATA FORMALISM AND NUMERICAL METHODS PHYSIOLOGICAL BACKGROUND MODEL BUILDING DISCUSSION APPENDIX A APPENDIX B APPENDIX C REFERENCES

It is generally assumed that Ca and Sr share the same main metabolic pathways, involving the gastrointestinal (GI) tract, kidneys, and bone, that are organized around an IDP (Fig. 2). However, this scheme is an oversimplification. For instance, the central IDP is, a priori, a multicompartmental substructure made up of the extracellular fluids (plasma included), soft tissues, intracellular sectors, and the so-called exchangeable mineral at the bone surface. It will also be useful to take into account a measure of complexity in the gut during model building and interpretation, because intestinal mineral absorption and secretion vary along the GI tract (6). The same is true for bone, because bone mineral metabolism involves, on the one hand, bone accretion and resorption, mainly linked to bone remodeling in adults via local osteoblast-osteoclast activities (26), and, on the other hand, the physicochemical processes by which the bone mineral solid phase is formed (nucleation and crystal growth) from solute ions, matures (changes in size, shape, and chemical composition), and dissolves (12). These processes also depend on the dynamic equilibrium-nonequilibrium (25) at the bone surface that governs its relations with ions in the adjacent extracellular fluid. All these processes can give rise to a complexity not shown in Fig. 2. On the contrary, this scheme may be simplified. For example, renal glomerular filtration and tubular reabsorption are so rapid that only the net urinary mineral excretion may be obtained from our kinetic data.

View larger version (20K): [in this window] [in a new window]   Fig. 2.   Overview of body Ca and/or Sr metabolism. IDP, internal distribution pool; GI, gastrointestinal tract.

Finally, the plasma Sr concentration in the normal adult is 0.5-1 µM, with >95% of the total body Sr (~400 mg) in the bone, as for Ca, and a normal diet provides a daily Sr intake of ~2 mg (10).

	MODEL BUILDING

TOP ABSTRACT INTRODUCTION PROTOCOL AND DATA FORMALISM AND NUMERICAL METHODS PHYSIOLOGICAL BACKGROUND MODEL BUILDING DISCUSSION APPENDIX A APPENDIX B APPENDIX C REFERENCES

Visual inspection of the mean plasma Sr values at each sampling time indicates a dose-dependent rise in plasma Sr concentration during AdP with, as expected, a slowing of the rate with treatment duration. It is followed by a regular decrease during PAdP (Fig. 1). However, there are indications for nonlinear dose dependence(s) of the kinetic behavior during each of these periods. The concentrations reached after each dose during the last part of AdP (mean value computed with data from day 14) show a departure from strict linear dose dependence, and the curve tends to flatten for the highest doses. During PAdP, the influence of the dose on the fall in plasma Sr may be estimated from the concentrations normalized to the mean value during the last part of AdP. The decrease in plasma Sr concentration is significantly slower for the smallest dose than for the other three doses (see Fig. 4B). These changes indicate that nonlinear dose-dependent processes are involved in plasma Sr kinetics. One of our objectives is to develop a model capable of justifying and describing these nonlinearities.

Our model of in vivo Sr metabolism was therefore built in three stages. The two preliminary stages used partial kinetic data to obtain the information required to construct the final model(s) on the basis of all available data. The general idea underlying stages 1 and 2 was to identify fairly qualitative (structural) features, such as the number of compartments, the nature of the involved nonlinearity (or nonlinearities), and its (their) optimal location(s) within the model that helped fit these partial kinetic data, including the dose dependence. We started with a minimal model structure and increased its complexity, inasmuch as each model was required to have a sound physiological meaning.

Stage 1: IDP Compartmental Substructure

Each dose was first considered separately, and the ability of linear or nonlinear compartmental models to fit each of the four sets of PAdP plasma data was tested. At least three compartments were required for strictly linear models, except for one dose, for which only two compartments were needed. Not only the Sr input function [f(t) in Eq. A4], but also all the other identified FTC, varied with the dose under these conditions. When one extrinsic z nonlinearity was introduced within the model, we chose to use system 2 in its generic form (see Assumptions 1 and 2 in APPENDIX A). With inclusion of the time-explicit (Eq. A1, a and b) or time-implicit (Eq. A2) version of z, in no case, did the goodness of fit obtained with this kind of nonlinear model significantly improve compared with that obtained by the best linear model. The parameter values defining system 2 and the corresponding system 1 FTF also showed large uncertainties (CV > 100%). However, the optimal structure for the nonlinear models was a two-compartment model at all the doses examined. These results seem to indicate that an extrinsic z nonlinearity is not detrimental to reproducing, separately, the rates of Sr loss during PAdP, although such nonlinearity is not entirely justified under this condition.

When we used another strategy, we obtained very different results. We examined the capacity of a given model structure plus a single set of parameter values [except for the input function f(t), which obviously differed with the dose] to fit the set of the four dose-related PAdP plasma data. Including one nonlinearity within the model significantly improved the goodness of fit over that obtained using a strictly linear model under these conditions, whatever the form (time implicit or time explicit) used for the generic z. The location of the nonlinearity within the structure was also important, because the criterion differs according to the kind of FTF that z_(n+1) modulates. For instance, the fit was less good when the irreversible (K₀₁) exit depends on z_(n+1), rather than the reversible (K₂₁) exit. Interestingly, we obtained the best fit for both peculiar model structures and when the time-explicit, rather than the time-implicit, form of z_(n+1) was employed.

The main properties and fitting characteristics of these structurally different models (number of compartments and presence, nature, and form of z nonlinearity) are shown in Table 1. These models had the same physiologically relevant feature (basic 3-compartment model; Fig. 3). These include two irreversible exits: one from compartment 1 and the other from compartment 2, which may reflect urinary excretion and internal uptake (probably by bone). The optimal criterion value was obtained with the nonlinear three-compartment structure using system 2 under the logistic form (Eq. A1, c and d) or the generic form (Eq. A1, a and b), provided it was implemented in its time-explicit form. The goodness of fit differed significantly from the other linear or nonlinear versions of the two- or three-compartment model. In other respects, the time-explicit form of z_(n+1) gives an index of cooperativity (p; see Assumption 2 in APPENDIX A) that is much greater than 1, i.e., a positive cooperativity, whereas a negative cooperativity (p < 1) is associated with the time-implicit form.

View larger version (11K): [in this window] [in a new window]   Fig. 3.   Basic 3-compartment nonlinear structure for internal distribution of Sr as identified from PAdP data. Bold arrow, K21 when modulated by z(n+1) as extrinsic nonlinearity (ENL) or intrinsic nonlinearity (INL). For linear model, K21 = k21.

As expected, when linear models with the three-compartment structure (Fig. 3) were used, k₀₁ and k₀₂ values were very inaccurate (CV largely >100%). This is due to a structural indetermination: only the sum of these irreversible exits could be identified. This indetermination seemed to be removed by using a nonlinear FTF. For instance, the set of parameter values defining system 1 was accurate enough (significantly different from 0) when K₂₁ was modulated by the time-explicit form of the logistic z_(n+1). However, this was not so for the accuracy of parameters estimated for system 2.

Different criterion values were produced (Fig. 4A) when optimization was carried out from the basic structure (Fig. 3), with K₂₁ as a logistic time-explicit ENL and various, but fixed, p values. However, these did not differ significantly from the optimal one for sufficiently high p values (p  2). A similar minimal criterion was maintained, despite different p values because of compensatory changes in the values of other parameters defining system 1 and/or system 2. The entire set of parameter values (systems 1 and 2) was well defined for fixed values of p  3. Otherwise, some parameters varied greatly with CV > 100% for system 2. Figure 4B illustrates the goodness of fit to plasma Sr data obtained with this model for p = 3. For each of the four doses, the data and the model response (compartment 1) were normalized to the theoretical y₁ values at 612 h to better appreciate the dose dependence of the kinetic behaviors. The statistical tests used to estimate the goodness of fit were satisfactory.

View larger version (10K): [in this window] [in a new window]   Fig. 4.   Fitting of 3-compartment model (see Fig. 3) with 1 logistic ENL, K21. A: minimal criterion values [square root of weighted residual sum of squares (WRSS)] as a function of different, fixed cooperativity orders (p). Dashed lines, confidence interval (P < 0.05) of minimal WRSS square root obtained after fitting PAdP data (see Table 1). B: fit of normalized model response (lines) to normalized PAdP data (symbols) for each of the 4 groups after withdrawal of oral Sr. Cooperative order p = 3. , Smallest dose (D1 = 170.5 mg/day); , D2-D4; error bars, SE for D1 and D4. See MODEL BUILDING for normalization procedure.

Thus a nonlinear three-compartment model, including an ENL showing positive cooperativity in its dependence on compartment 1, is able to account for the observed PAdP plasma Sr kinetics, including dose dependence, from a single set of parameter values. Obviously, it is too early for a physiological interpretation of the model. Figure 5A illustrates the predicted early kinetic behavior of z_(n+1) (z₄ for the optimal model using the logistic nonlinearity). The four doses studied resulted in very different behaviors. For the three higher doses, the expected z₄ maximum value (1.0) is reached at different times after initiation of oral Sr administration, but all were reached long before the beginning of PAdP. In contrast, the lowest dose resulted in a very late increase, with z₄ < 1 at the end of the AdP. Among the set of parameters, only f(t), i.e., the part of the entry relative to Sr administration, differs according to the dose: the input functions predicted from the nonlinear model are not strictly proportional to the dose, and the dose dependence presents a smooth, regular, S-shaped curve, at variance with the curve predicted from a linear three-compartment model (Fig. 5B). This result indicates that other nonlinear processes, probably related to the GI system, should be considered at a further stage of model building. This is all the more likely, given that, when searching to identify an optimal nonlinear model with f(t) compelled to a strict linear dose dependence, the minimal criterion reached is significantly higher than the optimal value obtained in the absence of this constraint.

View larger version (9K): [in this window] [in a new window]   Fig. 5.   Predictions from optimal logistic model retained from PAdP data. A: early z(n+1) kinetics for oral Sr administration [D1 (solid line), D2 (dashed line), D3 (dashed-dotted line), and D4 (dashed-dotted-dotted line)]. B: dose dependence of constant input (µM/h), as identified during AdP, from nonlinear model (dashed line) or its linear version (solid line). Error bars, SE estimated by variance-covariance matrix.

Finally, at this stage, it is possible, using the urinary data, to transform C_app to mass of Sr. With a model structure (Fig. 3) that 1) dissociates the irreversible exit from compartment 1 (k₀₁), including at least urinary excretion from another (internal) exit (k₀₂), and 2) optimally fits the PAdP plasma kinetic data, it is possible to estimate the relation of V_app to k₀₁. Because no significant variation in urinary creatinine (experimental data not shown) appears during the PAdP, Sr urinary clearance (0.286 l/h or 4.77 ml/min) can be estimated from linear regression between each experimental 12-h urinary Sr PAdP value and the predicted time-corresponding integral of y₁. The following relation is then retained

(3)

which gives rise to V_app > 60 liters with the above model.

Stage 2: GI Compartmental Substructure

The problem of the GI structure is not trivial, because structural indetermination becomes evident once the system becomes even moderately complex, when linear processes alone were assumed. Fortunately, the goodness of fit obtained with a model that includes an irreversible exit from compartment 1, k, in addition to the urinary excretion, k (see Assumption 8 in APPENDIX B) and a one- or two-compartment linear structure for the IDP, is significantly improved when the FTF from the GI compartment to compartment 1 becomes an INL via an M-M equation (Eq. B3) compared with the corresponding strictly linear model (Table 2). The increased number of GI compartments also improves the fit, in contrast to the IDP, for which one compartment is sufficient. Furthermore, the fit was not significantly improved (2.50 vs. 2.76) when this three-compartment minimal model (Fig. 6A) was identified without the f(t) dose-dependence constraint (Assumption 4 in APPENDIX B) and the estimated values of f(t) remained linearly correlated with the Sr intake.

View larger version (11K): [in this window] [in a new window]   Fig. 6.   Optimal nonlinear model obtained from FAdP data. A: structure. K15 (bold arrow) includes saturable process due to Sr intestinal absorption as described in Eq. B3. Compartments 4 and 5 belong to GI tract. At this stage (stage 2) of model building, k51 (dotted arrow) is not different from 0. B: fit of model compartment 1 (lines) to plasma Sr data for 4 doses of Sr. See Fig. 5A legend for explanation of lines.

The parameter values of the strictly linear model are very imprecise, whereas the presence of a saturable transfer function alone (Eq. B3, with k = 0) results in satisfactory accuracy (the parameters significantly differ from 0) for an optimal nonlinear structure with one IDP compartment (compartment 1) and two GI compartments (compartments 4 and 5). Compartment 4 does not exchange directly with compartment 1 (Fig. 6A). Accepting the input F₄₀(t), i.e., physiological and supplementary Sr entry, it serves as a source for the rest of the system, especially for the other GI compartment (compartment 5), from which Sr is nonlinearly absorbed into compartment 1 or excreted into the environment (feces).

The irreversible transfer from the first to the second GI compartment (k₅₄, in Fig. 6A) introduces a lag time [mean residence time ( = 1/k_ij) in compartment 4 of ~1 h] between the moment of oral Sr intake and its absorption, which is linked to the observed time of the plasma concentration peak (~4 h). An illustration of the fit with the FAdP experimental data is given in Fig. 6B. The relation of the GI compartment to compartment 1 is unidirectional (k₅₁ = 0) and restricted to the saturable part of Eq. B3 (k = 0). Indeed, the consideration of an intestinal absorption that is partly linear does not significantly improve the criterion value and makes the value of this parameter very inaccurate. This, together with the poor characterization of the IDP substructure and the lack of significance of the additional irreversible exit, k, is probably due to the short time (10 h) over which plasma kinetics were examined (FAdP). Therefore, we have not confined the use of these findings to the minimal structure (Fig. 6A), but we believe that a more complex relation between compartment 1 and the GI tract could be identified from the overall kinetics data (stage 3).

Stage 3: Overall Compartmental Structure

Figure 7A illustrates the minimal structure used to start this model building. It combines the two substructures from the preceding stages and includes five compartments for system 1: two (compartments 4 and 5) describe the GI tract, and the others (compartments 1-3) refer to IDP. We used only two nonlinear FTF: K₁₅, associated with an INL involving saturable and nonsaturable processes according to Eq. B3, and K₂₁, associated with a logistic ENL (Eq. A1, c and d) illustrated by the (n+1)th compartment.

View larger version (17K): [in this window] [in a new window]   Fig. 7.   Nonlinear multicompartmental models for Sr metabolism. A: initial 5-compartment structure used at the start of stage 3 of model building includes minimal structure characterized for IDP (stage 1, compartments 1-3) and for GI (stage 2, compartments 4 and 5). Compartment (n+1) is associated with z(n+1) extrinsic logistic variable modulating K21. Intrinsic nonlinearity of K15 with its saturable part obeys Michaelis-Menten Eq. B3. Pseudocompartment 6* integrates internal use of Sr from IDP. B: optimal 7-compartment structure obtained at the end of the modeling procedure. In addition to compartments in model described in A, this model is characterized by 1) addition of a linear compartment (compartment 7) and transformation of pseudocompartment 6* to compartment 6, 2) introduction of intrinsic Langmuir-type nonlinearity that acts at K21 and obeys Eq. C2b, and 3) addition of a 2nd extrinsic logistic variable z(n+2) that modulates K15; z(n+1) modulates K21 and K51.

The structural complexity was increased to obtain an adequate fit of the model response to all the data when the four Sr intakes were analyzed simultaneously (Table 3). Only one extra compartment was added to system 1 (or 2 extra compartments, if Eq. C3 was applied with the pseudocompartment 6* changed to compartment 6). This additional compartment was needed only when the overall criterion (23.81) approached the best value obtained (17.98). It is an IDP compartment (compartment 7) in a reversible linear relation to compartment 1, having a higher turnover than other IDP processes. A second logistic ENL [z_(n+2)] was included, with parameters different from those of the initial one, z_(n+1), because a sigmoidal influence on the relation of the GI compartment to compartment 1 appeared to be necessary. Finally, a simplified version of the Langmuir-type INL was introduced (Eq. C2b). Only the inhibition term of the complete Eq. C2 was used, not the saturable process. The Langmuir-type INL operates at the level of K₂₁, on an FTF that is also modulated by z_(n+1).

The final structural arrangement that includes these modifications is shown in Fig. 7B. It was developed using iterative optimization, testing various changes required to give a satisfactory behavior for a reasonable degree of complexity. The framework in which this was done is summarized in Table 3. Each of the nonlinearities (the M-M INL or the logistic ENL) can greatly improve the global criterion (models 2 and 3 in Table 3) over a strictly linear structure (model 1). However, this does not ensure a satisfactory fit to the PAdP data. The logistic ENL has parameter values very different from those identified in stage 1 in this situation: p is close to 1, and g and h_0(n+1) (Eq. A1, c and d) have much higher values than previously observed. The fit tended to be improved when both nonlinearities were simultaneously included in the model. However, elementary combination of the two kinds of nonlinearities (model 4) or simplistic interaction between them (models 5 and 6) did not give consistent results. Indeed, the sets of identified parameter values were erratic, with very large and unrealistic estimates of V_app (V_app > 200 liters), even though the saturable part of K₁₅ was, according to Eq. C1, converted to a logistic ENL through an extrinsic z variable that was identical to [z_(n+1) in model 5] or different from [z_(n+2) in model 6] that modulating K₂₁. The parameter values were also very inaccurately detected for systems 1 and 2 (CV > 100%). The three criterion values decreased substantially when the logistic z_(n+1) influenced K₁₅ and the reverse process K₅₁ [with K₅₁ = k₅₁z_(n+1), Eq. A3], whereas the computed CV of the set of parameter values became <100% for system 1, despite the addition of an extra compartment (model 7). On the contrary, some parameters in system 2 remained very inaccurate, with CV > 100%. Therefore, the Langmuir-type INL (Eq. C2) that regulates the FTF K₂₁ was added. This slightly morphed the overall criterion, but the parameter values were very inaccurate. This was resolved by using the simplified normalized version Eq. C2b, with the integral of the irreversible exit from compartment 2 toward the pseudocompartment 6* as inhibition variable (model 8). The Langmuir-type INL was further kept, because a posteriori examination showed that the system 1 parameter values were accurate, including the inhibition constant k in Eq. C2b. Explicit consideration of two logistic z had interesting consequences for the model-fitting properties. In fact, the way INL and ENL were combined had a great influence on the accuracy of parameter values, whereas the total number of parameters did not vary (models 9-11). The fit was moderately improved, but an optimal criterion associated with parameter values significantly different from 0 for system 1 and with CV < 100% for system 2 was obtained only when the K₁₅ and K₅₁ FTF involved in the relation of the GI compartment to compartment 1 depended on distinct extrinsic variables [z_(n+1) and z_(n+2) and conversely]. The values and accuracy of parameters defining system 1 (IDP and GI) and system 2 of a model in which K₂₁ and K₅₁ depended on z_(n+1), even though K₁₅ was modulated by the other extrinsic variable, z_(n+2), are reported in Table 4. Under such conditions, the system 2 parameters [mainly, g and g] vary greatly, a problem reminiscent of the indetermination observed during stage 1 of model building (Fig. 4A) and linked to the fact that the parameters p_(n+1) and p_(n+2) were not fixed at a given value. Furthermore, comparing model 12 with model 11 (Table 3) reveals that an explicit bone compartment (compartment 6), instead of a simple irreversible exit from compartment 2, does not significantly alter the fit. The calculated supplementary FTC, k₁₆, is very small, and changing the inhibition variable (compartment 6 for model 12 and compartment 3 for model 13) used for the Langmuir-type INL (Eq. C2b) has very little effect on criterion values, provided the inhibitory compartment has a low enough turnover.

The compartmental structure shown in Fig. 7B summarizes the main characteristics of the retained optimal models 12 and 13 (Table 3). These models give essentially identical theoretical responses using global criterion values. However, the values of partial criteria corresponding to FAdP and PAdP only approach the optimal values obtained at stages 1 and 2 of model building. The predicted kinetic behavior of compartment 1 compared with various mean values of plasma Sr throughout the experiment (AdP + PAdP) is shown in Fig. 1A for the four doses studied. Figure 1B shows a comparison of the theoretically predicted and experimental data for two of the doses used and at different times of AdP for the period immediately following the morning oral Sr dose.

Finally, we used another strategy to check the optimality of these structures. We progressively simplified the structure of model 12 (Table 3) by reducing the number of variables in systems 1 and 2 or the transfer processes that are "targets" for INL and ENL. Only one of the feasible simplifications (model 14 in Table 3) had a relatively good value of the global criterion, although it was significantly higher than that obtained with model 12 or 13. This model has the same number of variables and parameters as the optimal model. It differs in that K₅₁ is not modulated by z_(n+1). Under this condition, k₅₁ became zero, making the relation of the GI compartment to compartment 1 irreversible. K₁₅ thus defined the net intestinal absorption of Sr, including a nonsaturable process that is linearly dependent on the GI Sr concentration and a saturable process that is modulated by an extrinsic z variable. Also, the FAdP and PAdP criteria are far from their minimal value obtained at stages 1 and 2, and the system 2 parameter values were very inaccurate, with CV > 100% (model 14 in Table 3).

In summary, the final model (Fig. 7B) that fits all the experimental data well and gives satisfactorily accurate parameter estimates (statistically, at least, for system 1) requires several intricate nonlinear processes. The partial nonlinear structures identified during the preliminary stages 1 and 2 (Fig. 3 and 6A) can be included in a single model, provided the behavior of IDP and the GI compartment are influenced by ENL, i.e., by z functions kinetically expressed and showing high cooperativity in their dependence on plasma Sr concentration.

	DISCUSSION

TOP ABSTRACT INTRODUCTION PROTOCOL AND DATA FORMALISM AND NUMERICAL METHODS PHYSIOLOGICAL BACKGROUND MODEL BUILDING DISCUSSION APPENDIX A APPENDIX B APPENDIX C REFERENCES

In contrast to earlier kinetic studies, most recent compartmental modeling studies of in vivo mineral metabolism have used multiple tracers (stable or radioactive isotopes) given orally and/or intravenously and different sites of experimental measurements (21, 31, 39). These protocols are designed to develop complete models that represent all the metabolic pathways within some high-order, linear compartmental structures. These must be theoretically and practically identifiable (7) while remaining consistent with known physiology. The present modeling study has used quite different experimental data, because they concern the variations in plasma and urinary concentration of the mineral (Sr) itself after various oral doses of the mineral. Nevertheless, the complexity of the identified structure (Fig. 7B) is comparable with that of models for other minerals [Zn (19) and Mo (39)], even though the theoretical identifiability of the resulting model(s) has not been formally stated. The Sr metabolism described by system 1 has 7 compartments and 14 FTF. Thus, insofar as the present modeling procedure is confined to the characterization of a system structure plus the estimate of parameters and variables of physiological interest, this study seems satisfactory. This is mainly due to the qualitative and quantitative aspects of kinetic and dynamic information provided by the plasma and urinary data. We should recall that the data describe the rise in plasma Sr concentration during oral ingestion of Sr and its decrease after termination of oral Sr administration; we had a total of 38 samples per individual covering the experimental period. The short-term effect of oral Sr was also monitored, and data were obtained for four different doses of Sr. Finally, a total of 38 individuals (~10 per Sr dose) were examined. A major feature of our modeling procedure was that we looked for a single model, with just one set of parameter values defining a given compartmental structure, that could describe the entire set of data, including dose dependence. As asserted by Wagner (41), nonlinearities can be recognized in such a situation, but the corollary is that these nonlinearities significantly influence the identified kinetic behavior and the a priori identifiability. The nonlinearities of intrinsic or extrinsic origin become part of the model structure and are thus important for the model kinetic response. Conversely, a compartmental model built using each set of data separately will probably have a structure that is different from the nonlinear one described here. Use of separate data sets disregards one major expression of the nonlinearities, and the linear characteristics of the model should be strengthened with, as a consequence, inherent difficulties in the choice of model structure linked to structural indetermination or a priori unidentifiability. For instance, an unidentifiable linear model may become identifiable once it contains a nonlinearity, as shown by Godfrey (13).

We broke down the overall problem of model building by using two preliminary stages. This allowed us to look for structural and nonlinear characteristics of two subsystems before constructing a single integrated model. These subsystems deal with the multicompartmental substructure that includes the internal distribution of Sr, obtained from the later PAdP data, and the GI substructure, from the earlier FAdP data.

Each of these preliminary stages revealed interesting nonlinear properties of human Sr metabolism in a non-steady state. For instance, accounting for the PAdP data not only requires a peculiar nonlinearity, but we went further into the use of the nonlinear compartmental formalism (Eqs. 1 and 2). Indeed, this nonlinearity must be time explicit, which means that the extrinsic variable responsible for this nonlinearity does not follow the classical rapid-equilibrium (pseudo-steady-state) assumption but requires its time-explicit expression via the differential system 2 (see Remark 3). This is particularly important, because incorrect application of the asymptotic (time implicit) formulation of system 2, with the consecutive substantially improved fit compared with that obtained with linear models (Table 1), could result in an erroneous characterization of this nonlinearity: the time-implicit form shows a negative cooperativity in its dependence on compartment 1 Sr concentration, instead of the highly positive one found with the time-explicit version.

The second preliminary stage involved the kinetic analysis of the short-term effects of Sr immediately after the first oral dose. It revealed a direct link between our study and the present use of Sr as a marker of intestinal Ca absorption (42). Vitamin D increases the intestinal absorption of Ca and Sr (3), providing indirect evidence that Sr and Ca use the same transcellular pathway. Our analysis of the data collected during FAdP for all four Sr doses does ascertain that a saturable M-M-like process (Eq. B3) is involved in the human intestinal absorption of Sr. Our companion article (36) also shows that the parameters defining this process are quantitatively in agreement with Sr directly using some of the transcellular processes involved in the physiological regulation of Ca absorption.

With the overall nonlinear compartmental structure describing all the experimental data (Fig. 7B), it is possible to verify its biological significance, because the parameter values of the optimal models (model 12 or 13 in Table 3) have satisfactory CV values (CV always <30%; Table 4 for system 1 in model 12). The present study, as the first part of a more detailed study (36), checks the model validity without directly addressing the physiological interpretation of the INL and ENL. This may be done by examining the predictions of the model for human Sr metabolism under physiological conditions (steady state) and by simulating the effects of long-term oral Sr.

Initial Steady-State Predictions

Long-Term Predictions

View larger version (15K): [in this window] [in a new window]   Fig. 8.   Long-term predictions of Sr mass in bone (compartment 6) obtained by simulation of model 12 (bold curves) and model 13 (thin curves) when oral Sr is continued for 1 yr. See Fig. 5A legend for explanation of lines.

Thus the quality of the model predictions for Sr metabolism alone makes it possible to examine the physiological interpretation of the model nonlinear features. The physiological nature of the INL and ENL in the model, together with their interactions, must be identified to make our model significant. As for the regulation of many biological systems, these nonlinearities could involve more-or-less direct self-regulatory processes. However, this could be physiologically questionable for Sr for at least three reasons: 1) we know nothing of the physiological role of Sr; 2) there is no evidence for Sr homeostasis; and 3) the daily intakes of oral Sr used to reveal Sr-dependent nonlinearities were up to 1,000 times the normal intake. Thus we need to know whether the processes are Sr specific or directly related to the physiological regulation of Ca metabolism and its homeostasis, inasmuch as Sr and Ca are physicochemically similar.