On the Computation of Output Bounds for Compartmental In-Series Models under Parametric Uncertainty

In this work, the problem of obtaining tight output bounds for compartmental in-series models under parametric uncertainty is addressed. It is well-known that current methods used to compute a solution envelope may produce a significant overestimation. However, monotonicity analysis enables us to estimate a tight solution envelope. Our main aim is to get an equivalent model to the initial one, which is usually non-monotone, by means of a suitable combination of equations. In this new model the system monotonicity with respect to the uncertain parameters depends on the elimination rate values of the original model. If the equivalent model is monotone, no overestimation occurs in the computation of the output bounds. Keywords-Uncertainty; Parametric Uncertainty; Compartmental systems; Interval simulation; Monotonicity


I. INTRODUCTION
Mathematical models have appeared in many different real situations emerging from biology, economics, engineering, medicine, human sciences and many other research fields.The most common mathematical models used to mimic real processes are compartmental systems, in which each compartment represents a state of the system.
However, as a mathematical model is usually a simplified version of a real process, a mismatch between the behaviour of the model and the reality is produced.This mismatch yields non-modelled dynamics.Moreover, this kind of processes is also characterized by its variability, leading to parametric uncertainty.Therefore, the exact values of the model parameters are unknown, but they can be bounded by intervals.While there is only one possible behaviour for a model with constant parameters, parametric uncertainty produces a large set of different possible solutions.
Traditionally, Monte Carlo methods have been used to deal with uncertainty [1], owing to the fact that a large number of solutions can be easily computed.However, independently of the number of simulations executed, the output bounds obtained cannot ensure the inclusion of all the possible solutions [2].This inclusion guarantee is needed for error-bounded parametric identification and constraint-satisfaction problems.In the former, the range (or a tight enclosure) of the output trajectory must be computed and compared with measurements to estimate intervals for the model parameters guaranteeing data consistency.In the latter, the computed range must be compared with the constraints to be satisfied so as to obtain an inner and outer approximation of the output set for the decision variables.Furthermore, the computational cost of Monte Carlo methods increases proportionally to the number of simulations performed to cover the uncertain input space sufficiently.For these reasons, other methods have been considered to compute output bounds, such as region-based and trajectory-based approaches [3] mainly founded on interval analysis [4] and monotone systems theory [5], [6].
The aim of this work is to compute output bounds for compartmental in-series models with parametric uncertainty.That is, a tight solution envelope must be computed to ensure the inclusion of all the possible solutions for the model as well as to minimise the overestimation.Otherwise, if the overestimation is high, it could not be useful from a practical point of view, for instance, in an insulin therapy for diabetes patients [7].
This work has been organised as follows: In Section 2, a monotonicity analysis approach is introduced.In Section 3, compartmental in-series models are presented.In Section 4, a new method is proposed for the analysis of the system monotonicity with respect to the parameters.In Section 5, the proposed method is applied to compute the output bounds of a linear glucose model.Finally, Section 6 outlines the conclusions of this study.

II. UNCERTAIN SYSTEMS
Continuous-time systems under parametric uncertainty are described by an initial-value problem (IVP): where f is the vector function with components f i , x is the state vector, p is the parameter vector, and n p is the number of parameters.The solution of (1) is denoted by x(t; t 0 , x 0 , p).
We consider that the parameters and the initial conditions are unknown, but they can be bounded by intervals.Representing intervals in bold, interval vectors p and x 0 include all the possible values for the parameters p and for the initial conditions x 0 of the model, respectively.The set of possible solutions considering parametric uncertainty is denoted by x(t; t 0 , x 0 , p): The computation of solution envelopes plays a key role in the simulation of systems under parametric uncertainty.Such a computation can be performed by using one-step-ahead iteration based on previous approximations of a set of point-wise trajectories generated by the selection of particular values of the parameters p ∈ p and initial conditions x 0 ∈ x 0 by using heuristics such as a monotonicity analysis of the system [8].
Monotone systems have very robust dynamical characteristics, since they respond to perturbations in a predictable way.The interconnection of monotone systems may be studied in an analytical way [9], by considering a flow x(t) = φ(x 0 , t).A system is monotone if x 0 y 0 ⇒ φ(x 0 , t) φ(y 0 , t) for all t ≥ 0, where is a given order relation.Cooperative systems form a class of monotone dynamical systems [5] in which ∂f i ∂x j ≥ 0, for all i = j, t ≥ 0.
In order to calculate a solution envelope, an upper bounding model and a lower bounding model are computed.In an upper bounding model, the cooperative states with respect to the output take their upper bound value, while the monotone but non-cooperative states, known as competitive states, take the value of their lower bound.On the other hand, a lower bounding model is obtained taking account of the lower bound of the cooperative states, and the upper bound of the competitive states.In both cases, the non-monotone states are still computed as intervals that produce a significant overestimation.
The model parameters are considered as invariant states to carry out the monotonicity analysis, where

III. COMPARTMENTAL IN-SERIES MODELS
It is well known that a compartmental system consists of a finite number of interconnected subsystems called compartments.The interactions among compartments are transfers of material according to the law of conservation of mass.These are natural models useful for many application areas subject to that law, which appear in physiology, chemistry, medicine, epidemiology, ecology, pharmacokinetics and economy [10].The state variables of these systems represent the amount of material contained in each compartment and then they are restricted to be non-negative over time; that is, they belong to the broader class of positive systems.
A general in-series model composed of n compartments is represented in Figure 1.
e n u u n Fig. 1.Diagram of a compartmental in-series model.
An in-series model is named bidirectional if the fluxes between the compartments go forward and backward.However, if the fluxes just go forward, the in-series system is called unidirectional.Bidirectional in-series models are given by the equations: for i ∈ {2, ..., n − 1}, where the states of the model Q j (t), j ∈ {1, ..., n}, are the in-series compartments, and Q n (t) is the output of the model.Furthermore, u(t) and u n (t) represent the inputs and the parameters e j , j ∈ {1, ..., n}, are the elimination rates for each compartment, while k i,i+1 (•) and k i+1,i (•), i ∈ {1, ..., n−1}, are non-negative scalar functions that represent the flux from the compartment i to the compartment j and they may depend on the states of the model, i.e., k i,j (

IV. ANALYSIS OF THE SYSTEM MONOTONICITY
In this section, we analyse compartmental in-series models by focusing on the monotonicity of the dynamical system with respect to the states and the parameters.
The general system described by equations ( 2) can be non-monotone with respect to the states since it is not possible to determine the exact sign of the partial derivatives ∂ Qi(t) ∂Qj(t) , i, j ∈ {1, ..., n}, i = j.Notice that, for instance, the sign of the following partial equation cannot be determined: Therefore the monotonicity analysis cannot be accomplished with respect to the states and the parameters of the model.Nevertheless, as we are focused on the output of the model, this fact can be avoided by the transformation of this system into an equivalent system having the same output, given by: (3) for i ∈ {2, ..., n − 1}, where S i = n j=i Q j (t), ∀i ∈ {1, ..., n}.It is worth mentioning that all the fluxes k i,j in this new system may depend on the new states S i , such that k i,j (•) = k i,j (S 2 (t) − S 1 (t), . . ., S i+1 (t) − S i (t), . . ., S n (t)) Note that this equivalent system has been obtained by the combination of equations ( 2) and the output compartment is not modified, because S n (t) = Q n (t), which enable us to compute the output bounds of original system (2) by means of the conclusions obtained by the monotonicity analysis of new equivalent system (3) for i ∈ {2, ..., n − 1}: Biomath 1 (2012), 1210043, http://dx.doi.org/10.11145/j.biomath.2012.10.043 Under the conditions ∂ki,i+1(•) ∂Sj ≥ 0 and ∂ki+1,i(•) ∂Sj ≤ 0, ∀i, j : i ∈ {1, ..., n − 1}, j ∈ {1, ..., n}, we remark that the compartments of model ( 3) are cooperative among each other if e j ≥ e j+1 , ∀j ∈ {1, ..., n − 1}, since the partial derivative equations ∂ Ṡi(t) ∂Sj(t) , i, j ∈ {1, ..., n}, i = j are always non-negative.Furthermore, in this same case, the inputs u(t) and u n (t), and the functions k j,j+1 (•), j ∈ {1, ..., n − 1} are also cooperative, while the elimination rates e j , j ∈ {1, ..., n} and the functions k j+1,j (•), j ∈ {1, ..., n − 1} are competitive.

V. LINEAR GLUCOSE MODEL
In the sequel, we illustrate the result presented in the previous section through a linear glucose example.Namely, we turn the non-monotone system describing the Cobelli et al. model [11] into an equivalent monotone system, in which output bounds are easily computed without overestimation.
Plasma glucose concentration in blood is maintained within a narrow range with the help of the insulin hormone.Insulin is secreted by the pancreas with the role of reducing glucose concentration in blood.Under normal circumstances, a decrease in plasma glucose concentration is followed by a decrease in insulin secretion.On the other hand, insulin secretion increases when plasma glucose concentration increases, for instance after an ingestion.

. Diagram of the linear glucose model developed by Cobelli et al.
The analysis of glucose kinetics is essential to analyse the insulin secretion by the pancreas.Cobelli et al. developed a physiological model to study the insulin system and the control exerted by glucose on insulin secretion.This compartmental model describes the nonaccessible portion of the system, and it is composed of three compartments, as seen in Figure 2. The output compartment is the concentration of the accessible pool, displayed in the central position.The mass balance and measurement equations are given by: where Q 1 (t) is the accessible pool of the plasma glucose, Q 2 (t) and Q 3 (t) illustrate peripheral compartments, respectively, in rapid and slow equilibrium with the accessible pool, and the output of the model is given by the plasma glucose concentration G(t), which depends on the central compartment Q 1 (t).The parameter V I is the volume of plasma in the accessible compartment, the parameter EGP denotes the input, the constant parameters k 1,2 , k 1,3 , k 2,1 and k 3,1 are the fluxes between the compartments, while the parameters k 2,0 and k 3,0 stand for the elimination rates of the peripheral compartments.In this model there is no elimination rate in the accessible pool.The parameters values have been obtained from [12].
Performing a monotonicity analysis, it is deduced that the system is cooperative with respect to the compartments, as the partial derivatives among the model compartments are all non-negative.Furthermore, the input EGP is also a cooperative parameter, while V I , and the elimination rates k 2,0 and k 3,0 are competitive parameters.The monotonicity evaluation with respect to the fluxes k 1,2 , k 1,3 , k 2,1 or k 3,1 is not possible.
Cobelli et al. model ( 4) can be analysed as two compartmental in-series models interconnected, where the central compartment is the output of both in-series models.This system is equivalent to where As k i,i+1 (•) and k i+1,i (•), i ∈ {1, ..., n − 1}, are constant parameters then ∂ki,i+1 ∂Sj = 0 and ∂ki+1,i ∂Sj = 0.Moreover, as there is no loss in the output compartment and k 2,0 ≥ 0 and k 3,0 ≥ 0, then the Lemma IV.1 conditions hold.Thus, equivalent system ( 5) is cooperative with respect to the parameters k 2,1 and k 3,1 , the initial conditions and the input EGP .Furthermore, the equivalent system is competitive with respect to the parameters k 1,2 and k 1,3 , the elimination rates k 2,0 and k 3,0 , and the volume V I .The black dashed lines in Figure 3 display the computed output bounds, while the light grey lines represent several numerical simulations executed by the variation of the parameters and initial conditions values.First of all, the computation of output bounds is performed following the traditional monotonicity approach, where the parameters k 1,2 , k 1,3 , k 2,1 and k 3,1 have to be considered as intervals.The solution envelope in Figure 3a illustrates the overestimation produced in this case.On the other hand, when lemma IV.1 is applied the system is monotone with respect to all the states and parameters, thus none of them have to be considered as intervals.Therefore, it is possible to compute the output bounds without overestimation, as shown in Figure 3b.
Biomath 1 (2012), 1210043, http://dx.doi.org/10.11145/j.biomath.2012.10.043VI.CONCLUSION Different approaches in the literature have tackled the problem of parametric uncertainty for ordinary differential equations.In this work, a new method for the computation of output bounds on the compartmental in-series models is proposed.This method has been compared with previous approaches in a linear glucose model.
The most common method in the literature to compute a tight solution envelope is to perform a monotonicity analysis of the system for a trajectory-based approach.After this method is applied, only non-monotone compartments and parameters produce an overestimation in the computation of output bounds.This happens in compartmental in-series models, as they have non-monotone compartments and parameters.
Our proposal consists in obtaining an equivalent model by the combination of the differential equations of the original model, but without altering the output compartment.Then, by this way, a monotonicity analysis of the equivalent model is performed, obtaining that the new model is monotone with respect to all the compartments and parameters (cooperative or competitive) if the lemma IV.1 conditions meet.Thus, this approach allows us to compute a solution envelope adjusted to numerical simulations, in which no overestimation is produced, and computing just two simulations, one for the upper bound and another one for the lower bound.
Our proposed method outperforms previous approaches for the computation of output bounds on compartmental in-series models, as it computes the solution envelope without overestimation.