Identifying the sources of structural sensitivity in partially specified biological models
Quantifying structural sensitivity in models with uncertain component functions
In general, we consider a system of the form:
$$begin{aligned} dot{{mathbf{x }}}={mathbf{G}} left( h_1left( {mathbf{x}} right) ,h_2left( {mathbf{x}} right) ,ldots ,h_pleft( {mathbf{x}} right) , f_1left( {mathbf{x}} right) ,f_2left( {mathbf{x}} right) ,ldots ,f_{m-p}left( {mathbf{x}} right) right) , end{aligned}$$
(1)
where ({mathbf {x}}in {mathbb {R}}) is the vector of d state variables, (h_i, f_i:{mathbb {R}}^{d_i}rightarrow {mathbb {R}}) are the m different component functions describing the inflows and outflows of biomass, energy or individuals due to certain biological processes, with ({mathbf{G}} :{mathbb {R}}^mrightarrow {mathbb {R}}^d) being a composition function describing the general topology of the system. We consider that the precise mathematical formulation of the functions (f_i) are known (or at least postulated) with the only related uncertainty being the precise choice of their parameters. The functions (h_1,ldots h_p) are considered to have unspecified functional form. Instead, they are represented by bounds on their derivatives matching the qualitative properties we would expect from such a function. For example, the per-capita reproduction rate of a population is generally decreasing, at least at large population numbers, while a feeding term described by a Holling type II functional response of a predator should be increasing and decelerating. The (h_i) may also have quantitative bounds on their values:
$$begin{aligned} h_i^{text {low}}left( {mathbf{x}} right)1)), an important question remains: which of the unknown functions contribute the most to the degree of structural sensitivity in the system? The degree of structural sensitivity does not distinguish between the various sources of uncertainty and therefore cannot quantify the relative contributions of the unknown functions to the uncertainty in the model dynamics.
To determine the contribution of each unknown function (h_i), one can allow the error terms (left( varepsilon _1,varepsilon _2,ldots ,varepsilon _pright)) to vary with the goal of investigating how the degree of sensitivity varies with them. For the purpose of this section, let us denote the initial error terms by (varepsilon _i^0). We might be tempted to use the dependence on the (varepsilon _i) to perform global optimisation under certain constraints to find the best possible reduction of (left( varepsilon _1,varepsilon _2,ldots ,varepsilon _pright)). However, one should bear in mind that this analysis would depend on the base functions (hat{h}_i) considered. While these functions are ideally fitted to experimental data, they are only accurate within the error terms (varepsilon _i^0). Excessively reducing the (varepsilon _i) will force all admissible functions to conform strongly in their shape to these base functions far beyond their demonstrated accuracy of fit.
The dependence of the degree of sensitivity on (varepsilon _i) should therefore only be evaluated locally by calculating the gradient (left( frac{partial Delta }{partial varepsilon _1},ldots ,frac{partial Delta }{partial varepsilon _p}right) |_{left( varepsilon _1^0,ldots varepsilon _p^0right) }) giving the direction for the best local reduction of the errors. To adjust for the fact that the error terms may be of different orders of magnitude, when handling the vectors of error terms we should use the norm
$$begin{aligned} left| varvec{varepsilon }right| = sqrt{sum _{i=1}^{p} left( frac{varepsilon _i }{varepsilon _i^0}right) ^2}. end{aligned}$$
(9)
Working in this norm, the gradient needs to be weighted by the initial error terms to provide the direction for the best local reduction of the error terms, this is described by the following structural sensitivity gradient.
Definition 2
The structural sensitivity gradient in a model with p unknown functions each having an error of magnitude (varepsilon ^0_i) is defined as
$$begin{aligned} left( -varepsilon _1^0cdot frac{partial Delta }{partial varepsilon _1},ldots ,-varepsilon _p^0cdot frac{partial Delta }{partial varepsilon _p}right) |_{left( varepsilon _1=varepsilon _1^0,ldots ,varepsilon _p=varepsilon _p^0right) }, end{aligned}$$
(10)
where (Delta left( varepsilon _1,ldots ,varepsilon _pright)) is the degree of structural sensitivity of the system considered as a function of the error terms (varepsilon _i).
One possible problem with the structural sensitivity gradient is that the degree of structural sensitivity in the system may not be an increasing function of the magnitude of the errors. Consider the case that the exact system is structurally unstable, e.g. at a bifurcation point. Then no matter how small the error terms are, there may still be very high levels of structural sensitivity, while larger error terms may cause the level of uncertainty to decrease5. In this case, the structural sensitivity gradient will indicate that one or more of the functions has a negative contribution to the uncertainty of the system, and cannot be taken as a basis for sensitivity analysis.
An alternative approach to quantifying the individual impact of unknown functions which avoids this issue is the computation of partial degrees of sensitivity with respect to each (h_k). To do this, we fix every unknown function except (h_k), a set which we denote ({mathbf{H}} _{sim k}), by fixing the (x_j^*), (h_ileft( {mathbf{x}} ^*right)), and (frac{partial h_i}{partial x_j} left( {mathbf{x}} ^* right)) that are consequently determined by the isocline equations. Denoting by (V_k) the cross-sections of V where only (h_k) varies, and the cross-sections for ({mathbf{H}} _{sim k}) by (V_{sim k}), the local partial degree of structural sensitivity can be defined as follows.
Definition 3
The local partial degree of structural sensitivity with respect to (h_k), is the degree of structural sensitivity in the model when (h_k) is unspecified and all other functions (h_i in {mathbf{H}} _{sim k}) are fixed:
$$begin{aligned} Delta _k({mathbf{H}} _{sim k}):= 4 cdot int _{V_{k_{text {stable}}}} rho _{mathbf{H} _kvert {mathbf{H}} _{sim k}} , d{mathbf{H}} _k cdot left( 1 – int _{V_{k_{text {stable}}}} rho _{mathbf{H} _kvert {mathbf{H}} _{sim k}} , d{mathbf{H}} _k right) , end{aligned}$$
(11)
where (rho _{mathbf{H} _kvert {mathbf{H}} _{sim k}}) is the conditional probability density function on ({mathbf{H}} _{sim k}):
$$begin{aligned} rho _{mathbf{H} _kvert {mathbf{H}} _{sim k}} = frac{rho }{int _V rho d{mathbf{H}} _{sim k }}, end{aligned}$$
with (rho) the (joint) probability distribution over V.
The local partial sensitivity (Delta _k({mathbf{H}} _{sim k})) is a function of ({mathbf{H}} _{sim k}) in that it depends upon the particular values at which the elements of (V_{sim k}) are fixed. As with the degree of structural sensitivity, it can be interpreted as either the probability that the stability of the given equilibrium will be different for two independent choices of the function (h_k) when the (h_{sim k}) are fixed at the given values, or in terms of variance as (Delta _k({mathbf{H}} _{sim k})=4cdot {text {Var}}_k(Yvert {mathbf{H}} _{sim k})) (Y is the Bernoulli variable for stability). If the joint probability distribution (rho) is uniform in V, then (Delta _k({mathbf{H}} _{sim k})) can be expressed purely in terms of the fraction of the volume of (V_k) which gives a stable equilibrium:
$$begin{aligned} Delta _k({mathbf{H}} _{sim k}) = 4 cdot frac{int _{V_{k_{text {stable}}}} d{mathbf{H}} _k}{int _{V_k} d{mathbf{H}} _k} cdot left( 1 – frac{int _{V_{k_{text {stable}}}} d{mathbf{H}} _k}{int _{V_k} d{mathbf{H}} _k} right) . end{aligned}$$
(12)
To obtain a global measure for the sensitivity of the model to (h_k), we can take the average of (Delta _k) over (V_{sim k}).
Definition 4
The partial degree of structural sensitivity with respect to (h_k) is given by
$$begin{aligned} bar{Delta }_k := int _{V_{sim k}} rho _{mathbf{H} _{sim k}} cdot Delta _k({mathbf{H}} _{sim k}) , d{mathbf{H}} _{sim k} end{aligned}$$
(13)
where (rho _{mathbf{H} _{sim k}}) is the marginal probability density function of ({mathbf{H}} _{sim k}).
Recalling the variance-based interpretation of the degree of sensitivity, we obtain (bar{Delta }_k = 4cdot E_{sim k}({text {Var}}_kleft( Yvert {mathbf{H}} _{sim k}right) )). In other words, (bar{Delta }_k) gives the scaled average variance when all functions except (h_k) are fixed. We can also relate the partial degree of structural sensitivity to indices used in conventional variance-based sensitivity analysis. Dividing (bar{Delta }_k) by the overall degree of structural sensitivity in the model gives us (frac{bar{Delta }_k}{Delta }=frac{E_{sim k}({text {Var}}_kleft( Yvert {mathbf{H}} _{sim k}right) )}{{text {Var}}(Y)}=S_{T_k}), the total effect index23 of (h_k) on the stability of the equilibrium. This is a measure of the total contribution of (h_k) to the sensitivity—both alone and in conjunction with the other functions ({mathbf{H}} _{sim k}). However, since the space of valid functions V is in general not a hypercube, the functions (h_i) are not independent factors, and a total decomposition of variance is not possible. Indeed, even if the joint probability distribution (rho) is uniform, the marginal probability distribution (rho _{mathbf{H} _{sim k}}) will generally not be: instead it will equal the volume of the corresponding cross-section (V_k) for ({mathbf{H}} _{sim k}), divided by the volume of V. An alternative to using the partial degrees of sensitivity would be to consider the first-order sensitivity indices (S_k=frac{{text {Var}}_kleft( E_{mathbf{H} _{sim k}}left( Yvert h_k right) right) }{{text {Var}}(Y)}). However, these do not take into account possible joint effects of the (h_i) on the structural sensitivity of the system, so a small (S_k) does not indicate that (h_k) is not a source of sensitivity, whereas (bar{Delta }_k=0) means that (h_k) does not contribute to the structural sensitivity in the system at all.
Similarly to the gradient of the total degree of sensitivity (Delta) as a function of the respective error tolerances, the vector (left( -bar{Delta }_{h_1},ldots , -bar{Delta }_{h_p}right)) needs to be scaled by the elements of (varepsilon ^0) to give us the optimal direction of decrease in (Delta) if the error terms (varepsilon _i) are subject to a proportional reduction. This is described by (left( -varepsilon _1^0 bar{Delta }_{h_1},ldots , -varepsilon _p^0bar{Delta }_{h_p}right)).
Outline of an iterative framework of experiments for reducing structural sensitivity
When dealing with partially specified models, an important practical task is the reduction of the overall uncertainty in the system by decreasing the uncertainty in the system processes (i.e. the unknown model functions). Here we propose an iterative process of such a reduction based on improving our empirical knowledge of the uncertain functions (h_k).
As a starting point, we assume that experiments have produced data on the unknown functions (h_1,ldots ,h_p), to which we can fit some base functions (hat{h}_1,ldots ,hat{h}_p) with initial errors (varepsilon _1^0,ldots ,varepsilon _p^0). We assume that it is possible to perform additional experiments on all uncertain processes in order to obtain more data such that the (varepsilon _i) can be decreased, but with the natural constraint that the total error can only be reduced by a magnitude of (0 More
