## Abstract

It is well recognized that models in the life sciences can be sensitive to small variations in their model functions, a phenomenon known as ‘structural sensitivity’. Conventionally, modellers test for sensitivity by varying parameters for a specific formulation of the model functions, but models can show structural sensitivity to the choice of functional representations used: a particularly concerning problem when system processes are too complex, or insufficiently understood, to theoretically justify specific parameterizations. Here we propose a rigorous test for the detection of structural sensitivity in a system with respect to the local stability of equilibria, the main idea being to project infinite dimensional function space onto a finite dimensional space by considering the local properties of the model functions. As an illustrative example, we use our test to demonstrate structural sensitivity in the seminal Rosenzweig–MacArthur predator–prey model, and show that the conventional parameter-based approach can fail to do so. We also consider some implications that structural sensitivity has for ecological modelling: we argue that when the model exhibits structural sensitivity but experimental results remain consistent it may indicate that there is a problem with the model construction, and that in some cases trying to find an ‘optimal’ parameterization of a model function may simply be impossible when the model exhibits structural sensitivity. Finally, we suggest that the phenomenon of structural sensitivity in biological models may help explain the irregular oscillations often observed in real ecosystems.

## 1. Introduction

At present, one of the major limitations to the predictive power of biological models is their potential sensitivity to variation of their constituent functions. The conventional way of assessing the possibility of errors in model predictions consists of checking the sensitivity of the results to variation in the model parameters [1–3], but recently it has been shown that such an approach can be somewhat limited—in particular, a large variation in parameters can result in small deviations in model outcomes, while a small variation in the mathematical formulation of functions can sometimes result in a vastly different model behaviour [4–8]. This phenomenon is referred to as structural sensitivity, and it can potentially cause a serious limitation to the predictive power of biological models [5,8].

More precisely, structural sensitivity occurs in a model *M* when a close model *M*_{ε} shows either qualitatively different behaviour (for example, a change in the stability of an equilibrium) or a significant quantitative difference in model predictions, even though the distance between the functions constituting the equations for both models (as well as their first derivatives) does not exceed a certain small *ε* [8], and the qualitative properties of these functions are conserved (the sign of derivatives and the curvature, vanishing at zero, etc.). Structural sensitivity is different from an absence of the well-known property of structural stability, since structural stability holds when all the systems within a *sufficiently small*(in the *C*^{1}-sense) distance of the initial model are qualitatively consistent [9], whereas in the case of structural sensitivity we consider that perturbations of the model functions are in a small but finite range. The concept of structural sensitivity is of vital importance in biological modelling since in most experiments and observations we are unable to distinguish between close functions that differ from one another by less than 5–10% [10–16].

Surprisingly, verification of the existence of structural sensitivity can be necessary even if the use of a particular parameterization of functional dependence has a solid theoretical background. Consider, for instance, the famous Monod equation used to describe the feeding rate of a consumer/predator, which has a concrete theoretical basis in ecology [17–19]. Since we often deal with the dynamics of a population, and not of an individual, we need to use the averaged response over the entire population of consumers, so even in the case where the consumption of each individual is described exactly by a Monod formulation, averaging this nonlinear function over the whole population (since the handling time and the maximum consumption rate will differ between individuals) can result in a different mathematical expression. Also, when considering models on larger space scales we often cannot directly implement the functional relations that are observed on small spatial scales [16,20–22]. Thus we somehow need to scale up the microscale dynamics, and the expression for the functional response may be changed [22]. Finally, to describe species interactions in a real ecosystem, we need to take into account other important factors such as feeding history, complex feeding behaviour, adaptation to food and even short-term evolution [23,24]. Inclusion of each of these factors can seriously alter the well-known Monod formulation, and clearly this holds true for any other functional relation such as a closure term and growth rate.

The major problem with verifying structural sensitivity in models is that for a given mathematical formulation of a functional dependence *f* we need to somehow check for a difference in the model outcomes for ‘all’ close functions. Such a procedure may seem impossible in principle: the functional space to which *f* belongs usually has an infinitely large number of dimensions, and one cannot check them all. For this reason, investigation of structural sensitivity is usually done by choosing a few concrete parameterizations of *f* and comparing the resultant model outcomes in each case [6], but clearly such a sensitivity analysis is rather subjective, since it strongly depends on the choice of the forms of *f* that are compared.

In this paper, we provide a rigorous non-parametric method for the investigation of structural sensitivity in biological models. In particular, we consider how the local stability of stationary states will be altered by perturbing the functions of the model. Our method is based on a semianalytical approach and we quantify structural sensitivity by introducing a concept of the ‘degree’ of sensitivity. We demonstrate our test by using it to reveal structural sensitivity in the seminal Rosenzweig–MacArthur predator–prey model, with respect to variation in both the functional response of the predator/consumer and the per capita growth rate of the prey, and show why the conventional methods of structural sensitivity analysis based solely on variation of the parameters can be misleading in certain cases. In addition to this extensive example, we discuss the implementation of our test for a number of ecological models published previously in the literature, for which it is able to demonstrate the existence of structural sensitivity. We argue that, in the case of structural sensitivity, finding a model function that can accurately mimic experimental/observational results consistently may be impossible (even if we can find a function that will reproduce a given set of results). Finally, we put forward the hypothesis that structural sensitivity in mathematical models may partially account for the irregular oscillations of species densities often observed in natural ecosystems and mesocosms.

## 2. General framework

We consider a biological model based on a system of autonomous ODEs given by
2.1where the vector function *F*=(*F*_{1},*F*_{2},…,*F*_{n}) is taken to be sufficiently smooth. Each *F*_{i} is the summation of terms describing the inflow and outflow of biomass, energy or individuals on a different trophic level, and depends on a set of functions *f*_{ij}, representing growth rates, numerical and functional responses, closure terms, etc.,
2.2We assume that the mathematical formulation of some of these *f*_{ij} is well known (or postulated), and we denote these functions by *g*_{ij}. The only potential uncertainty regarding these functions consists in the correct choice of their parameters. We denote by *h*_{ij} the functions whose exact shape/formulation is unknown, and for which we have only information on general properties such as the sign of their derivatives, whether they vanish at zero and the existence of any thresholds. Such systems are known in the literature as ‘partially specified models’ [25], and, ideally, the number of unknown functions should be much smaller than the number of known ones.

To obtain a concrete biological system from (2.1), we need to specify all functions in the model by choosing parameters in the known functions *g*_{ij}, and by specifying the parameterization of the *h*_{ij}, as well as the corresponding parameters. The resultant model will have a set of attractors—that can be stationary states, periodic, quasi-periodic or chaotic—that is determined by the choice of the right-hand-side functions. In this paper, we shall focus on the influence of uncertainty in the formulation of the unknown functions *h*_{ij} upon the stability of any stationary states in the model.

Let us assume for now that there are *p* functions in the model with an unknown mathematical formulation, which we denote by *h*_{1},…,*h*_{p} (dropping the *h*_{ij} notation here for simplicity). To complete the system, it is usual to choose particular realizations of the functional shapes of *h*_{1},…,*h*_{p} using empirical observation or theoretical reasoning. The question is, how sensitive are the model outcomes to the choice of the formulation of *h*_{1},…,*h*_{p}? In the case of a given equilibrium, we know from the implicit function theorem that, provided the system is structurally stable, if we replace the function set *h*:={*h*_{1},…,*h*_{p}} in the system by some other sufficiently close function set the new system will have a corresponding equilibrium *x** close to and we can ask whether the stability of the new equilibrium *x** is altered from the stability of To answer this question, we need to introduce some suitable definitions of an *ε*-neighbourhood of *h* in functional space that take into account the qualitative properties that the must satisfy:

### Definition 2.1

A function set belongs to the -neighbourhood of if the following conditions are satisfied:
2.3where *Q* denotes a class of functions having certain specific properties (e.g. decreasing functions, convex functions and functions passing through the origin) and *Ω* is a compact (closed and bounded) domain in for some *m*≤*n*—since the function set may depend only upon *m* of the *n* dimensions of the state space. Provided that we require *Q* to contain only functions with second derivatives contained within some interval (*A*_{1},*A*_{2}), the compactness of *Ω* ensures that any in the -neighbourhood of *h* must have restricted first derivatives. Therefore in this case, the metric defined by definition 2.1 is implicitly a *C*_{1}-metric in the sense that functions that are close in the sense of definition 2.1 must also have close derivatives. For this reason, we suggest that *Q* be always defined such that the second derivatives of its functions are within a given range.

In other words, the definition requires that, for any *x*=(*x*_{1},…,*x*_{n}), the corresponding point must be contained in the *p*-dimensional neighbourhood *B*_{ε}(*h*_{1}(*x*),…,*h*_{p}(*x*)), as well as the each satisfying certain desirable properties. This definition considers the absolute distance between *h* and , but we are often interested in the relative error in experimental measurements etc. so we should also propose an alternative definition of the *ε*_{Q}-neighbourhood of *h* that uses the relative distance between the functions.

### Definition 2.2

A function set belongs to the -neighbourhood of *h*:=(*h*_{1},…,*h*_{p}) if the following conditions are satisfied:
2.4where *Q* and *Ω* have the same meaning as before. Hereon, a neighbourhood obtained from a specific choice of definition 2.1 or 2.2 shall be denoted an -neighbourhood/-neighbourhood, while ‘*ε*_{Q}-neighbourhood’ shall be used as a general term which can refer to either.

Finally, we may also consider a mixed definition where for some subdomains of *Ω* we require condition (2.4) and for the others condition (2.3). We may want to consider the relative distance when the functions take small values, for example, but also want to consider the absolute distance when they take very large values.

We base our test on the rigorous definition of structural sensitivity proposed in Cordoleani *et al.* [8]. According to this definition, a model exhibits structural sensitivity if there is another model in its *ε*_{Q}-neighbourhood with an asymptotic behaviour qualitatively different or significantly quantitatively different from that of the original (more specifically, two models exhibit ‘qualitatively different asymptotic’ behaviour if their orbits are not homeomorphic, and ‘significantly quantitatively different asymptotic behaviour’ if the Hausdorff distance between their attractors exceeds a given threshold *σ*). In this paper, we shall focus on the qualitative changes aspect of structural sensitivity, when the stability of an equilibrium point *x** is not consistent across all choices of in the *ε*_{Q}-neighbourhood of *h*, where *ε* is fixed beforehand. The absence of structural sensitivity in this sense is different from structural stability, which requires that the qualitative system dynamics remain consistent under *infinitely* small perturbations of model functions. Absence of structural sensitivity becomes structural stability if we take the limit *ε*→0, so a system that is not structurally stable will necessarily exhibit structural sensitivity for *ε*>0. For this reason, we should assume that the system we are testing is structurally stable. It is well known that, in a non-degenerate case, the stability of an equilibrium is given by the characteristic equation of the corresponding linearized system,
2.5where the are functions that can be analytically computed on the basis of the Jacobian matrix for a specified function set . *x** is stable if and only if the roots of (2.5) have negative real parts [29], and this gives us the condition of stability in terms of *x**, and the first-order partial derivatives
where *m* is the dimension of the union of the domains of the unknown functions. Although it is generally impossible to solve (2.5) analytically, it is easy to do it numerically for fixed *x**, and In particular, if *n*=2, we only need to check the sign of the determinant and the trace of the Jacobian.

In order to check for structural sensitivity in an unspecified system, we need to make an initial choice of function set, *h* (in practice, this choice should be made such that it fits experimental data somehow), and then test each function set in the *ε*_{Q}-neighbourhood of *h*. The space of such function sets usually has infinite dimension—even if they consist of only one function—but we see from (2.5) that the only values given by the function set that determine the stability of a given fixed point *x** are the (*n*+*p*+*m*⋅*p*) values *x**, and (recall that *n* is the dimension of the state space, *p* is the number of unknown functions and *m* is the dimension of the domain of the unknown function set . Additionally, since from (2.1) *F*(*x**)=0 must hold if *x** is to be an equilibrium, this gives us *n* of these values—either analytically or numerically—provided that the others are specified. Instead of considering the functions themselves, we can consider the remaining (*p*+*m*⋅*p*) unknown values and proceed as follows. We consider the range of possible *x** in the vicinity of that are the stationary states of the system with some from the *ε*_{Q}-neighbourhood of *h*, i.e.
2.6where is the system specified by choosing the function set . Similarly, for each possible *x**∈*X*, we can then consider the range of values that the functions can take that allow to stay in the *ε*_{Q}-neighbourhood of *h*, and for any permissible combination of *x** and we can consider the range of values that the first-order partial derivatives, can take with still remaining in the *ε*_{Q}-neighbourhood of *h*. This range will be finite, since we choose *Q* such that its functions have second derivatives that do not exceed some range (*A*_{1}, *A*_{2}).

This approach gives us the projection of the *ε*_{Q}-neighbourhood of *h* into the ‘(*X*-*H*-*H*′) space’—that is, the space with elements —and since we can use characteristic equation (2.5) to verify the stability of *x**, we can then check the entire projected neighbourhood for values that alter the stability of the equilibrium from that of in model (2.1) with the base function set *h*. If the stability is altered for certain values of *x**, and in the projection of the *ε*_{Q}-neighbourhood of *h*, we can conclude that the model exhibits structural sensitivity.

The key advantage of this approach is that the *ε*_{Q}-neighbourhood of *h* in (*X*-*H*-*H*′) parameter space has at most *p*+*m*⋅*p* dimensions, and so may be covered numerically, whereas the corresponding neighbourhood in the function space may have infinite dimensionality and therefore be impossible to check numerically. The only potential difficulty with the approach is determining what the correct *ε*_{Q}-neighbourhood in the (*X*-*H*-*H*′) space actually is, since the definition of the *ε*_{Q}-neighbourhood is non-local: we require and there may be other global constraints in *Q*.

Another issue faced by this approach is how to quantify the *degree* of structural sensitivity displayed by a system. Ideally, we would like the degree of structural sensitivity in a system to reflect the probability that any two function sets—independently chosen from the *ε*_{Q}-neighbourhood at random—yield different predictions for the stability of the given equilibrium. In order to achieve this, we first need to choose a suitable probability distribution over the *ε*_{Q}-neighbourhood in the (*X*-*H*-*H*′) space in order to determine the probability of choosing a point in the stable or unstable region. One of the simplest approaches we can take is to assume that the probability distribution is uniform over the *ε*_{Q}-neighbourhood, in which case the probability of selecting a function projecting to a point in a given region of the *ε*_{Q}-neighbourhood is simply given by the proportion of its area/volume to the total area/volume of the *ε*_{Q}-neighbourhood.

The next logical step in finding an improved probability distribution would be to, where possible, weight each point in the (*X*-*H*-*H*^{′}) space by some measure of the density of functions that we are projecting onto this point. The problem here, however, is that we are dealing with a projection from an infinite dimensional function space into an at-most (*p*+*m*⋅*p*)-dimensional parameter space. As a pragmatic solution to this difficulty, we propose a measure of the density of functions *ρ*, projected to a point in the *ε*_{Q}-neighbourhood of *h*, that relies only on the space : we consider the area/volume of points (*x*,*y*_{1},…,*y*_{p}) in for which the graph of at least one function in the *ε*_{Q}-neighbourhood of *h* can pass for which there exists some function, *f*, in the *ε*_{Q}-neighbourhood of *h* such that (*x*,*y*_{1},…,*y*_{p})=(*x*,*f*_{1}(*x*),…*f*_{p}(*x*))—*while also satisfying* and

Thus, the probability of selecting a function that yields stable dynamics can be defined as
2.7where *ρ* is the functional density, *V* _{1} is the domain corresponding to stability and *V* _{ε} is the total *ε*_{Q}-neighbourhood of *h* after their projection into the finite-dimensional (*X*-*H*-*H*^{′}) space. The case where the probability distribution is taken to be uniform over the *ε*_{Q}-neighbourhood can be obtained from (2.7) by taking *ρ*≡1. Using this probability distribution, the probability of two randomly chosen functions yielding different predictions will be given by 2*V* ⋅(1−*V*), since the choice is considered to be independent. This value will range from 0 (the *ε*_{Q}-neighbourhood will consist of only stable or only unstable regions) to 0.5 (there is an equal probability that a function from the *ε*_{Q}-neighbourhood will predict a stable or unstable equilibrium), so we scale this value by two to obtain the degree of structural sensitivity,
2.8In the next section, we shall show a few examples of how our investigation can be carried out in simple cases in which there is a single unknown function.

## 3. Implementation of the structural sensitivity test in models

Here, we demonstrate our method by using it to reveal structural sensitivity in the classical Rosenzweig–MacArthur predator–prey model [26] given by
3.1and
3.2where *P* and *Z* are the densities of prey and predator, respectively. is the per capita growth rate of prey; is the predator functional response, *k* is the trophic efficiency coefficient and *m* is the predator mortality.

### (a) Investigation of structural sensitivity with respect to the functional response of the predator

Following the work of Fussmann & Blasius [6], which previously uncovered structural sensitivity in such a system, we first consider that the unknown function is the functional response of the predator. In this case we consider that the prey growth rate is fixed as the logistic function, i.e. and that the mortality term is constant. We obtain the function class *Q* by imposing the following constraints on :
3.3which define a Holling type II (i.e. monotonically increasing and saturating) functional response over *Ω*=[0,*P*_{max}] [4,27], with the second derivative bounded below by the parameter *A*<0. To test structural sensitivity of the model we consider the functions from the *ε*_{Q}-neighbourhood of a certain base function *h*, where the class of functions *Q* consists of all functions satisfying (3.3). Figure 1*a*,*b* shows a standard base functional response, together with the upper and lower bounds of any functions within its *ε*_{Q}-neighbourhood—using definitions 2.1 and 2.2 from §2, respectively. The upper and lower limits of in these cases are constructed by plotting the boundaries and respectively. (As with the two types of *ε*_{Q}-neighbourhood, we will generally use *h*_{ε+} and *h*_{ε−} to refer to both kinds of boundary unless the specific notation is necessary.) The possible range of values taken by the stationary prey density *P** is defined by the intersection between the horizontal line *m*/*k* and the curves *h*_{ε+} and *h*_{ε−}, which follows from considering the non-trivial predator isocline from equation (3.2). Clearly, for all , the value is uniquely defined by , so we do not need to consider variation of this value.

For each stationary *P**∈(*P*_{1},*P*_{2}) there exists a range of possible derivatives, but not all of them are feasible. It can be shown (see electronic supplementary material, appendix A) that there exists a function satisfying (3.3) such that its derivative at *P** is equal to *DP* if and only if the three following conditions are met:
3.4By scanning all feasible values of *P** and *DP* we can check the stability of the interior stationary state (*P**,*Z**) for all functions in the *ε*_{Q}-neighbourhood of *h*. The stability condition for this equilibrium (based on a simple analysis of the Jacobian) is given by [26]
3.5The generic behaviour of model (3.1)–(3.2) with a predator functional response satisfying (3.3) is well known. For a fixed the stability of the stationary state is determined by the carrying capacity *K*: for a small *K*, the interior stationary state is stable; an increase in *K* will eventually result in the system's destabilization via a Hopf bifurcation—a phenomenon known as the ‘paradox of enrichment’ [28,29].

We shall implement our test for structural sensitivity with respect to the functional response using two different base functions: the Monod parameterization (3.6) and the hyperbolic tangent function (3.7), two of the most commonly used Holling type II functional response terms in the literature [6,18,30],
3.6and
3.7where *a*_{i} and *b*_{i} are parameters with standard meanings. It is easy to see that both *h*_{1} and *h*_{2} satisfy conditions (3.3). The parameters in the expressions are chosen in such a way that functions *h*_{1} and *h*_{2} are close to each other in terms of absolute difference, which does not exceed *ε*=0.12 (figure 1*c*).

To perform our investigation, we used the conditions (3.4) to numerically compute the region in the (*P**−*DP*) space that corresponds to the *ε*_{Q}-neighbourhood of *h*, and then used criterion (3.5) to determine the regions of stability and instability of the fixed point (*P**,*Z**) in this space. We then calculated the degree of structural sensitivity *Δ*, assuming *ρ*≡1 in (3.8) for simplicity. Figure 2*a* shows the dependence of *Δ* on the carrying capacity, *K*, for different base functions and values of *ε*, using definition 2.1 of the distance between two functions (absolute closeness of functions). We see that, in all cases, *Δ* is low for small values of *K*—this is because, initially, the *ε*_{Q}-neighbourhood is dominated by the region of stability—but as *K* is increased, *Δ* increases to 1 as the region of stability in the *ε*_{Q}-neighbourhood shrinks, before dropping off again as the region of instability starts to dominate. This follows the paradox of enrichment, since, for all choices of functional response, increasing *K* will destabilize the equilibrium *eventually* (*Δ* reaches zero if we consider values of *K* higher than those shown in figure 2*a*), but, significantly, *Δ* is far from zero across a large range of *K*, which indicates that there are significant regions of both stability and instability in these cases (figure 3), and hence we have sensitivity of the model. Therefore, while in this model the paradox of enrichment does take place, the exact carrying capacity, *K*, at which the Hopf bifurcation happens is strongly dependent upon the formulation of the functional response, and can vary across a huge range. Moreover, the degree of structural sensitivity observed is similar regardless of whether the Monod or hyperbolic tangent base function is used, indicating that for close base functions the structural sensitivity test does not depend upon the choice between them.

Figure 2*b* shows the dependence of the degree of structural sensitivity on *K* when definition (2.2) of the closeness of functions is used (i.e. relative closeness). One can see that for small *ε* the results depend on the base function chosen—the Monod function in this case exhibits structural sensitivity for a much smaller range of *K*. This discrepancy can be explained by noting that, while the absolute error between the two base functions is small, the relative error is *ε*≈0.5 or *ε*≈1, depending on which of the two functions is considered the base function, so they are quite far removed from each other's *ε*_{Q}-neighbourhood in terms of relative distance. Therefore this result is not an inconsistency on the part of the method. When *ε* is large enough (*ε*=0.7) to compare with the relative distance between the two base functions, the Monod and hyperbolic tangent base functions again exhibit similar degrees of structural sensitivity. It is of interest that using different concepts of the closeness of functions to define the *ε*_{Q}-neighbourhood of a given base function *h* can result in different predictions regarding structural sensitivity. This can be seen from comparison of the Monod curves in figure 2*a,b*: when relative closeness is considered, the system with a Monod base function exhibits only structural sensitivity for values of *K* between 0.3 and 3—a much smaller range of *K* than for the same system when absolute distance is used.

In figure 3*a*,*b* we have plotted the stability portrait in the (*P**−*DP*) space for the Monod and hyperbolic tangent base functions, respectively, where the absolute distance is used with *ε*=0.1, and *K*=1.2 (figure 2*a*). In the green region, (*P**, *Z**) is a stable equilibrium; the red region indicates an unstable equilibrium; and the dark blue region corresponds to the values of *P** and *DP* that are not found for any function in the *ε*_{Q}-neighbourhood of *h*. It should be noted that all regions are actually projections of the corresponding regions in infinite dimensional function space.

The (*P**−*DP*) representation can shed some light on the limitations of conventional sensitivity analysis—which is based only on a variation of the model parameters for a fixed functional form. We can vary the parameters *a*_{i} and *b*_{i} in equations (3.6) and (3.7) and consider all possible combinations that belong to the corresponding *ε*_{Q}-neighbourhoods: these regions in the (*P**–*DP*) space are denoted by the azure domains in figure 3*a*,*b*. In the cases shown, varying the model parameters does not indicate any structural sensitivity, but it is clear from the more complete analysis that there is extensive structural sensitivity which variation of parameters misses completely. This is because variation of the parameters for a fixed functional form allows us only to vary the values *P** and *DP* in a thin strip directed in a direction roughly parallel to the Hopf bifurcation line. This ‘blinkered’ approach therefore gives a misleading representation of the degree of structural sensitivity in the system, and so we can rely only on a conventional parameter-based analysis of sensitivity when the carrying capacity *K* is close to the Hopf bifurcation value for a given *h*, since in this case the azure domain will intersect the boundary between the stability and instability domains.

### (b) Investigation of structural sensitivity with respect to the growth term of the prey

We now consider the structural sensitivity of the system with respect to variation in the formulation of the growth term of the prey, in the cases where the functional response is fixed either as the Monod or as a hyperbolic tangent functional response with the same parameters *a*_{i} and *b*_{i} as used for the base functions in §3*a*. We impose the following constraints on :
3.8which define a logistic-type growth term, with a second derivative bounded above and below by the parameters *A*_{2} and *A*_{1}, respectively. We take our base function (which we denote by *r*(*P*)) to be the classic logistic growth function
3.9which clearly satisfies conditions (3.8). We require that be within an absolute distance *ε* of this base function, that is, it must satisfy where *r*_{ε−}(*P*)=*r*(*P*)+*ε* and *r*_{ε+}(*P*)=*r*(*P*)−*ε*.

When varying the growth term, since the functional response *h* is fixed, the equilibrium prey density *P** is given by the constant *h*^{−1}(*m*/*k*), and the stability condition in the case of an arbitrary growth function is given by
3.10We consider the values and and note that there exists a valid function in the *ε*_{Q}-neighbourhood of *r* satisfying and if and only if the following conditions hold (see electronic supplementary material, appendix B):
3.11Using these conditions, we can check the stability of (*P**, *Z**) for all the functions in the *ε*_{Q}-neighbourhood of *r*, and plot the corresponding regions in the (*rP–DP*) space.

Figure 4 shows the degree of structural sensitivity *Δ* (constructed for *ρ*≡1) plotted against the carrying capacity, *K* (cf. figure 2). Several values of *ε* are considered, and both the Monod and hyperbolic tangent functional responses from §3*a* were used. In this case, the degree of structural sensitivity shown by the system with respect to variation of the growth term is completely different depending on which functional response is taken—with the Monod functional response, varying the growth term will alter only the stability of the equilibrium for *K* between 0.5 and 2, while with the hyperbolic tangent functional response there is a high degree of structural sensitivity even when *K*=16. This substantial difference is in spite of the fact that the two functional responses used are very close together, and it suggests that, if we have two unknown functions, performing analysis of structural sensitivity with respect to each in turn will not be enough to determine the extent of structural sensitivity with respect to variation of *both* of them together.

In figure 5*a*,*b*,*d*,*e*, we present the *ε*_{Q}-neighbourhood of *r* (with *ε*=0.1) along with its regions of stability and instability in the (*rP–DP*) space for two values of *K* with the Monod and hyperbolic tangent functional responses, respectively. This figure shows us that an increase in *K* does not affect the Hopf bifurcation line in this space, but instead shifts the *ε*_{Q}-neighbourhood towards the *DP*=0 axis. The difference in structural sensitivity between the Monod and hyperbolic tangent functional responses is entirely due to the resultant change in the value of *P**.

For the growth function, it is quite simple to compute the ‘functional density’ *ρ* of each point in (*rP–DP*) space, as discussed at the end of §2: as the domain of the function is one dimensional, we can bound the functions passing through a point by the parabolas of maximum positive and negative curvature, and use this approach to decide which points have functions passing through them which also pass through (*P**,*rP*) with derivative *DP* (see electronic supplementary material, appendix D, for the specific criteria and details of the method in this case). Figure 6 illustrates these points for two given values of (*rP*,*DP*). Figure 5*c*,*f* shows the functional density of the *ε*_{Q}-neighbourhood in the case of figure 5*b*,*d*, respectively. The functional density is shown to be greater towards the centre of the domain, but does not vary a great deal (the range is 0.13–0.59 in figure 5*c*, but only falls below 0.5 in non-red regions, and the range is 0.6–0.82 in figure 5*f*) so including this information will not change the character of the graphs in figure 4 a great deal in this particular model.

## 4. Discussion and conclusions

Although it is well recognized that biological models can be quite sensitive to the mathematical formulation of the model functions [4–8], the conventional sensitivity analysis is mostly undertaken by varying the model parameters for fixed functional forms ([1–3] and many other references). This is mainly because it has been held as ‘evident’ that we cannot check all possible mathematical formulations of a given functional dependence *f*, since the functional space to which *f* belongs has an infinite number of dimensions. In this paper, we have challenged this widespread opinion and suggested a simple but rigorous test of structural sensitivity that reveals the effect of the choice of the mathematical formulation of model functions on the local stability of equilibria. The main idea is to project the infinite dimensional functional space onto a finite low-dimensional space consisting of the values of the functions and their derivatives at the equilibrium. Our method is exhaustive, in the sense that no function which alters the stability conditions will be missed, provided the resolution of the mesh is high enough. Importantly, we consider only functions with the same qualitative global properties as the initial model function (e.g. monotonically increasing/decreasing, sign of curvature, etc.), so the observed violations of stability are more surprising than in some early reported cases, where the perturbations were small in magnitude but could result in a qualitative change of function properties—changing the sign of the curvature, for instance [31].

Implementation of our sensitivity test reveals the main reason for the limitations of the conventional parameter-based methods of sensitivity. As is shown, for instance, in figures 3 and 5, varying only the model parameters for a fixed mathematical formulation often results in a displacement in the functional space mainly confined to a limited number of directions, and if the main direction of displacement in the functional space is roughly parallel to the bifurcation hypersurface, then variation of the model parameters will not result in a crossing of this bifurcation hypersurface (figure 3*a*,*b* and figure 5*b*,*d*,*f*), unless the base function lies very close to it in the first place (figure 5*a*). Therefore, this limited approach can reach a misleading conclusion that the stability of the given equilibrium will be consistent across a wide range of parameters, while considering perturbations of the base function in the *ε*_{Q}-neighbourhood in the direction traversal to the bifurcation hypersurface will reveal a stability change (figure 3*a*,*b* and figure 5*b*,*d*,*f*). For example, in the electronic supplementary material, appendix C, we analytically construct the *ε*_{Q}-neighbourhood corresponding to variation of parameters for a fixed mathematical formulation *f*and derive the conditions of structural sensitivity for small *ε*_{Q} when the number of parameters is equal to two.

Projecting the infinite dimensional function space onto a low-dimensional subspace containing only the equilibrium points and the values of the functions and their derivative(s) at the equilibrium has a drawback in that information about measures of regions in the function space may be lost—in particular, the measure of the neighbourhood considered, and its sub-domains of stability and instability. Formally, we can solve this problem by introducing certain weights of the points in the considered subspace to represent their impact. In this paper, we introduce the functional density *ρ*, which characterizes the relative amount of functions that can be potentially constructed to take the given values and derivatives at the equilibrium (e.g. figure 6). Another important issue arises from the fact that some parameterizations can have more biological significance than others, and so when computing the degree of sensitivity *Δ* we would ideally give more consideration to those functions that are more biologically relevant and disregard less meaningful parameterizations. However, weighting the possible parameterizations according to their biological significance requires detailed information about the nature of the underlying processes which is often not available.

Our sensitivity test requires the choice of a certain base function, which defines the *ε*_{Q}-neighbourhood in which we will consider all possible functions. Therefore it is pertinent to ask whether the choice of base function influences the outcome of our structural sensitivity test. The answer to this fundamental question depends on how the magnitude of the perturbation of the base function compares with *ε*. In any case, very small perturbations of the base function should not influence the qualitative behaviour, since the model is assumed to be structurally stable, but additionally our experience—based on considering a number of well-known models (see below)—leads us to conclude that, provided the perturbation of the base function is not too large relative to *ε*, the result of the structural sensitivity test will not change significantly (figure 2). In practice, both the tolerance, *ε*, and the base function should be determined from the same experimental dataset, and so our sensitivity test can be assumed to be independent of the base function.

One tentative extension of the method is to consider all possible perturbations of the initial model *M*_{0} by relaxing the restrictions on the properties of functions constituting this model. In other words, we can potentially consider all close models *M* such that their right-hand-side functions as well as the first derivatives are close to those of *M*_{0} and test how such perturbations will affect the stability properties: the corresponding sensitivity analysis can be easily done in a way similar to the one seen earlier. The issue with this approach, however, is that such an analysis can be misleading, since we risk focusing on biologically meaningless ‘exotic’ perturbations: for instance, allowing an increase in predator density to result in an increase in the overall rate of change of prey density. This approach may also lead to consideration of models that simply contradict important experimentally obtained knowledge on a given species (e.g. assuming Holling type III functional response for a species that always exhibits Holling type II response). Thus when considering perturbations of models we should always take into account any vital biological information that we have *a priori*.

Along with the simple Rosenzweig–MacArthur predator–prey model we have also tested several other models for structural sensitivity, with varying degrees of complexity. These consisted of several plankton models of the ocean [32,33], chemostat-like models [8,34,35] and one competition model [36]. In particular, we have included our analysis of the model described in Fussmann & Blasius [34] as appendix E electronic supplementary material to demonstrate that our model is still applicable to more complex models (the model in question is four dimensional, and conditions for linear stability cannot be obtained analytically). In each case we varied only one functional relation (either the functional response of consumer/predator or the per capita growth rate), other functions being kept constant. Surprisingly, we found that structural sensitivity can be observed in all the aforementioned models, which leads us to believe that structural sensitivity is a fairly widespread phenomenon in ecological modelling. However, the presence of structural sensitivity in the mentioned models does depend on other parameters—within some parameter ranges it is not observed—as well as the magnitude of *ε*, i.e. the maximum possible perturbation from the base function. Finally, implementation of our test has shown that the bounds of the magnitude of the functional relation's second derivative can play an important role, while very often these are unknown, and can at best be estimated from experimental data.

The key question is what we need to do in the case that we find that a given model exhibits structural sensitivity. Certainly, we should be rather careful regarding the model's quantitative predictions—which may be quite inaccurate given the usual uncertainty of model functions and large scattering of points in laboratory experiments [10,12–14]. One possible course of action could be to use data on the experimental population dynamics to reconstruct the unknown underlying model functions [7,25,37,38]. One can try to reconstruct as closely as possible the ‘true’ model functions that should be used in modelling the given experimental mesocosm or ecosystem. However, in the case that we have structural sensitivity, an attempt to reveal such ‘true’ parameterizations for a given set of equations may simply be in vain. Functional relations that we use in models are often oversimplifications of a large number of factors—for instance, the functional response *h* of a consumer is not only a function of food and/or the population density of the consumer itself, i.e. *h*=*h*(*P*,*Z*)—and even if the environmental conditions (temperature, light intensity, etc.) are kept fixed, other factors such as adaptation and evolution of the prey can influence the functional response [23,39,40]. The influence of spatial scale can also be a crucial factor: implementing parameterizations obtained on a small (i.e. laboratory) spatial scale in modelling the dynamics over larger temporal and spatial scales can be erroneous for a number of reasons [20–22,41]. As a result, the ‘true’ functional relation may simply not exist as a function of the given state variables of the model, and can be defined only up to a certain accuracy *ε*. In this case, uncertainty in the choice of a parameterization may arise not because of some experimental errors, but as a result of internal drift of the functional relations in real ecosystems owing to factors not included into the simplified model. For this reason, any case where an optimal parameterization of model functions mimics the experimental data well, but a close parameterization results in a pronounced deviation, should be considered to be rather suspicious: if the data had been obtained another time, an entirely different optimal parameterization may have been found. On the other hand, if the dynamics observed in a biological system appear to be relatively consistent, structural sensitivity in a corresponding model can be an indicator that something is wrong with the model construction, since we should expect some small variation in the biological functions that, according to a structurally sensitive model, would result in significant variation in the population dynamics. In this case, we probably need to stop searching for a function providing us with a fantastic fitting and make necessary changes to the model structure, by including adaptation and evolutionary factors, for instance [23,40].

One particularly interesting consequence of structural sensitivity in models is that it may help explain the observed variability of population sizes in natural ecosystems. The widespread opinion is that irregular oscillations of species densities are a consequence of either internal chaotic dynamics [42,43], the influence of environmental noise [44,45] or the interplay of both factors [46]. The phenomenon of structural sensitivity allows us to propose another scenario of such irregular species oscillations, since in reality functional relations between system components are not fixed but slowly change in time, because of both environmental fluctuations and evolutionary factors. In corresponding models, we can describe such a drift of functions, for instance, as a random walk in function space within a certain domain of size *ε*_{Q}. Owing to structural sensitivity, small variations in functions may result in transitions between the stability and instability subdomains of the functional space (see figures 2 and 5), and, moreover, within the instability domain the amplitude as well as the period of resulting oscillations can prominently change through small variations in the model functions [6–8]. In this way, small fluctuations in model functions may be amplified and result in large-amplitude irregular oscillations of species densities, which may present themselves in the original biological system. Such a mechanism has been proposed in previous works [47], but fluctuations in model functions were considered to be only due to variation in model parameters.

We strongly encourage researchers who construct and implement various models in the life sciences to implement our structural sensitivity test on their models, along with the conventional parameter-based tests. They may find that some results that seem to be robust to variation of model parameters turn out to be sensitive to a slight variation of model functions, and while this may be an unpleasant surprise at first, identification of structural sensitivity can help modellers better understand the true underlying mechanisms in the model they are studying, and inform potential future refinements. The proposed test can be easily extended to consider other features of Hopf bifurcations, for instance, to check how sensitive the type of Hopf bifurcation is (i.e. subcritical or supercritical) with respect to variation of model functions: in this case we will need information regarding the second derivatives of the model functions [9]. The test can also be easily modified in order to check for structural sensitivity in discrete-time population models.

Finally, the proposed test does have its limitations as well, and more work is required to improve it. The test is designed to detect qualitative changes to system dynamics, but structural sensitivity can also manifest itself in terms of large quantitative changes in a model's predictions, with potentially calamitous consequences to the predictive power of a model [5,8]—developing methods to detect when this is the case would be in the interest of predictive modelling in a wide range of disciplines. In particular, one case in which the approach used here is not applicable is when we need to reveal the sensitivity of oscillatory dynamics (regular or chaotic) to the choice of the model functions: for instance, the amplitude and period of any limit cycles may be sensitive to the choice of parameterization [5,6,8]. In this case, local values of the functions do not give us enough information: we need to know the entire shape of the functions between the maximal and minimal values of the species density oscillations. The same problem emerges when we are interested in the influence of functional variation on transient (i.e. non-asymptotic) dynamics. We are currently working on extending the test presented here, and are planning to address these problems in our future works.

- Received August 25, 2012.
- Accepted September 27, 2012.

- © 2012 The Author(s) Published by the Royal Society. All rights reserved.