## Abstract

Mathematical models in biology are highly simplified representations of a complex underlying reality and there is always a high degree of uncertainty with regards to model function specification. This uncertainty becomes critical for models in which the use of different functions fitting the same dataset can yield substantially different predictions—a property known as structural sensitivity. Thus, even if the model is purely deterministic, then the uncertainty in the model functions carries through into uncertainty in model predictions, and new frameworks are required to tackle this fundamental problem. Here, we consider a framework that uses partially specified models in which some functions are not represented by a specific form. The main idea is to project infinite dimensional function space into a low-dimensional space taking into account biological constraints. The key question of how to carry out this projection has so far remained a serious mathematical challenge and hindered the use of partially specified models. Here, we propose and demonstrate a potentially powerful technique to perform such a projection by using optimal control theory to construct functions with the specified global properties. This approach opens up the prospect of a flexible and easy to use method to fulfil uncertainty analysis of biological models.

## 1. Introduction

Mathematical models of various ecological systems based on differential equations often have the troublesome property of structural sensitivity in which the use of functional forms that are quantitatively close and qualitatively similar yield contradictory dynamical behaviour [1–4]. Ecological models are particularly vulnerable to structural sensitivity for two main reasons. First, biological processes have a high level of complexity, and so the precise form of any function chosen to represent them cannot be justified, because it is necessarily a simplification of the true relation. Second, biological data often have substantial error terms, so significantly different functions can fit the same dataset equally well. Generally, although structural sensitivity is fairly well acknowledged, the conventional approach in mathematical biology is to stick with a particular functional form, with parameter variation being the full extent of any attempt to deal with uncertainty [5–7]. Simply varying parameters, however, is insufficient to check for and deal with structural sensitivity, because models have been shown to be highly sensitive to the formulation of model functions while remaining robust with respect to parameter perturbations [2,3]. The use of a particular parametrization cannot even be fully trusted if it is derived mechanistically, because such mechanistic derivations always involve many simplifications, so we should not expect it to remain valid if we take into account heterogeneity of a population, larger time or space scales and fluctuating environmental factors (see the Introduction to [8] for a more detailed explanation).

A better approach is to explicitly include the uncertainty in model functions by considering partially specified models [4,9]. In such models, we leave unknown functions unspecified apart from requiring that they satisfy some qualitative criteria inherited from the biological problem being modelled—we may require a function to be increasing, for instance, or to pass through the origin, etc. Note that this approach is based on similar ideas to the seminal works of Gause and Kolmogorov as early as the 1930s [10,11]. Many properties of models are locally determined, such as the number and linear stability of equilibrium points, and it is quite easy to deal with such properties in partially specified models. Indeed, investigation of the stability of coexistence stationary states is a central part of ecological modelling [12–14]. Near a (hyperbolic) equilibrium, the behaviour of the system is equivalent to that of its linearization, which is completely determined by the value of the equilibrium density and the local values of the unknown functions and their derivatives at this point [15]. We can then treat these values as independent parameters and construct a generalized bifurcation diagram in this new parameter space. This is the basis of the approach known as ‘generalized modelling’ [16–18], which links generalized bifurcation diagrams to real-world data by considering a transformation of model functions in which the Jacobian depends not only on the derivatives of unspecified functions at the equilibrium, but also on the ‘elasticities’, which can theoretically be measured from data. However, because all generalized parameters still need to be measured at equilibrium states, obtaining these measurements can be impossible in the case of an unstable equilibrium; so while we can still try to check the range of the generalized parameters that we consider to be realistic, this will be a somewhat subjective choice.

A more appealing approach would be to consider the entire possible span of the shape of the unspecified biological function—without a need to take measurements from the whole system in which the process is embedded. Based on this concept, a new framework was recently developed to perform bifurcation analysis while taking into account the entire span of the shape of the function [4,8]. The crux of the method is to take the entire infinite dimensional space of functions admitted by the data and respecting the global constraints of biological realism, i.e. constraints on the function over its whole domain, and to project it into the generalized bifurcation space. This gives us a closed region in the bifurcation space that consists only of those generalized parameters that can be taken by functions fitting the data and the constraints of biological realism. Using this principle, one can quantify the uncertainty in the system and even carry out a probabilistic bifurcation analysis of the model, for example defining the probability of having oscillations in the case where the exact shape of the function is not specified [4,19].

This new method of structural sensitivity analysis can rigorously cover all possible model functions, instead of requiring model functions to be restricted to certain (often arbitrary) equations, which makes it particularly useful when modelling biological systems with a high degree of uncertainty. However, its widespread practical application will hardly be possible without resolving the key mathematical question of whether or not functions exist satisfying certain local and global restrictions—that is, they and their derivatives take certain values at certain points, and stay within given bounds across the whole domain. In previous papers on this framework for structural sensitivity analysis, all proofs and derivations (using geometrical arguments) were limited to the particular case where the global bounds themselves satisfied the properties required from the model function—for instance, for a functional response of Holling type II shape, the quantitative bounds must themselves be increasing, saturating and concave-down functions. Another drawback is that for each set of qualitative and quantitative function properties, we need to carry out analytical calculations which might well be enough to put off more practically orientated biological modellers, and this inflexibility obstructs the development of software to perform the analysis automatically. Finally, there remains the challenge of constructing appropriate weightings on the neighbourhood of relevant function values in the generalized parameter space: some of these values should be treated as marginal because they correspond to functions that fit the data less well than others, and there is no straightforward way to construct such a weighting.

In this paper, we make a first attempt to develop a general approach to resolve the fundamental issues raised and to point the way towards the development of automated software to make the structural sensitivity framework more accessible for both modellers and biologists. The new approach consists of *approximately* determining whether a valid function satisfying a set of local and global properties exists by constructing a functional that ‘rewards’ functions for adhering to our criteria of validity (both quantitative and qualitative) and penalizes them for straying from the criteria, finding the function that maximizes this functional, and checking whether this ‘highest scoring’ function satisfies the constraints. To find the function maximizing the functional, we use optimal control theory—the area of mathematics concerned with finding control parameters of a differential equation that maximize functionals of the resulting solution—in particular, the Pontryagin's maximum principle. To demonstrate the efficiency of this method, we consider a Rosenzweig–MacArthur predator–prey model [12,20] and use our approach to reveal structural sensitivity in this system when we consider functions within complex bounds. We show that our approach outperforms the usual tactic of detecting sensitivity by simply varying parameters of fixed functional forms. Finally, we show how optimal control theory can be used to estimate a relevant ‘weighting’ of each set of functions used in modelling.

## 2. General framework

### (a) Defining a partially specified biological model

Partially specified models represent an attempt to include uncertainty in models by considering one or more model functions as unspecified functions, which are not given by a particular equation [9]. In this paper, we shall consider partially specified ordinary differential equation models with only one function left unspecified, *x* is a single real-valued variable (note that the method can be easily generalized to an arbitrary number of unspecified functions)
*g*_{1}, … , *g*_{s} are fully specified functions, with only their parameters being unknown; the function *G* determines how the subfunctions are assembled together, and is assumed to be known. Our aim is to check whether the model is sensitive to the functional form of *f*, and thereby determine to what extent we can trust the predictions of a corresponding model with a precise equation chosen for *f*. In order to link the partially specified model to data and theory, we consider that *f* is required to remain between two boundary functions
*f* can be specified in different ways, depending on the situation: they may be taken from experimental data, in which case they will be the corresponding error bounds representing, say, 95% confidence intervals, or they could be used as a way of representing a purely hypothetical dataset in our analysis, so that we can check how the model reacts to changes in the whole data profile. Alternatively,

We may also want to restrict *f* and its derivatives at a certain set of values
*I*_{i} are indexing sets that determine which derivatives of *f* we restrict at each point *x*_{i} is the inflection point. Because we usually cannot know the exact value at which we have an inflection point, we should treat *x*_{i} as a parameter of the investigation in this case, to be varied over a range.

Naturally, it is more complicated to investigate partially specified models than to investigate models in which all functions are fixed, and innovative ways are required to deal with this difficulty. However, when we have high uncertainty in our model functions, the use of only fully specified models may not be sufficient. Consider that most generic function spaces are infinite dimensional. If we have a partially specified model and choose a particular function, determined by two parameters, say, to represent an unspecified function, then this only encompasses a two-dimensional subset of the infinite dimensional space of valid functions. So unless we have a very good reason for choosing this parametrization, its analysis is likely to be highly misleading [2,3]. If we can work directly with the partially specified model, then our analysis will be general, as any such model necessarily includes every model with a particular valid parametrization as a special case.

### (b) Structural sensitivity analysis of partially specified models

The framework first introduced in [8] for the investigation of partially specified models centres on the fact that for many system properties we do not need to know the entire shape of the constituent functions. For instance, local to a hyperbolic equilibrium, system behaviour is determined solely by the values taken by the functions and their derivatives at this equilibrium [15]. In the notation of §2a, we need to know only the values *x** and *f*(*x**) to determine if *x** is an equilibrium density. Given these values and *f* itself. In this way, we can construct bifurcation diagrams in a ‘generalized parameter space’, in which *x**, *f*(*x**) and the relevant derivatives of *f* at *x** are treated as if they were model parameters. This idea is the basis of both the generalized modelling framework [16,18] and earlier classical works [10,11].

Given the generalized parameter values *x**, *f*(*x**), … , *x** = 0.1 and *f*(*x**) = 0.9, then, in this case, clearly no function can take these values and lie within the boundaries. A generalized bifurcation plot should therefore look like figure 3—with the dark blue region representing impossible generalized parameter values, surrounding the finite region of possible ones. Finding this finite region is then the main challenge—afterwards, we can determine the range of possible system behaviour, and compute the relative sizes of the domains of this region that correspond to different dynamical regimes. In the case where different domains have comparable size (or measure), we have significant structural sensitivity in the model. Thus, the main idea of this framework for structural sensitivity analysis is finding a *projection:* from the infinite dimensional space of functions which fit the data and theoretical assumptions, onto the low-dimensional generalized bifurcation space.

To find such a projection, we need a way to determine whether, given a set of the values *x**, *f*(*x**), … , *x** as one of the points in (2.4), at which we locally restrict *f*, *f*′, … , *f*^{(p)}, and check whether *any* function exists satisfying these new constraints. Once we have a criterion for determining the existence of such a function, we can scan the generalized bifurcation space consisting of the values *x**, *f*(*x**), … ,

The main mathematical challenge for the framework is to prove the existence of a function *f* that respects a set of global and local constraints (equations (2.2)–(2.4), respectively). In [2,4], a geometric approach was used to obtain such a projection up to the first derivative for a particular class of qualitative constraints, in the case where the boundaries

### (c) Determining the existence or not of a function satisfying global/local constraints

One potential approach to address the issues mentioned above and make the framework of [4,8] more flexible is to use methods from optimal control theory to find the projection from the space of valid functions into the generalized bifurcation space *approximately*. Here, we shall propose an algorithm for doing this. In this section, we consider the equilibrium *x** along with *f* satisfying a set of global and local constraints in the form of equations (2.2)–(2.4).

In order to address this problem approximately, we apply optimal control theory. Broadly speaking, optimal control theory aims to find the parameter input of a system of differential equations such that the solution maximizes an objective functional, which scores solutions of the system based on some desirable criterion (alternatively we can aim to minimize a cost functional penalizing undesirable criteria). The basic idea to check for the existence of a function *f* satisfying (2.2)–(2.4) is to create an objective functional that rewards functions for adhering to the global bounds (2.2) and (2.3) and build a set of differential equations that yield a candidate for the function *f* as a solution (in such a way that all solutions must satisfy the local restriction provided by (2.4)). The currently considered values

We specify our objective functional as taking the following form:
*F*_{i} are given by
*γ* describes the width of the transition layer at the boundary around *γ* yielding narrower layers and a steeper transition. The inclusion of the (*p* + 1)th derivative is discussed later, along with the form of *R*. Such a function *F*_{i} is shown in figure 1: it imparts high values to functions that satisfy the constraints, without discriminating too much between functions that do (as would be the case if *F _{j}* had a pronounced peak, for instance). The functional

*I*will yield a higher value for a function

*f*if it satisfies constraints (2.2) and (2.3) for more

*x*-values. Note that the construction of this functional involves a softening of the hard bounds in (2.2) and (2.3): we no longer demand our function to satisfy these inequalities strictly. Therefore, we could dispense with the hard bounds in (2.2) and (2.3) altogether, and replace (2.6) with a normal distribution or a similarly peaked function, without changing the framework. In this case, instead of trying to determine the set of all values in the generalized bifurcation space that could be taken by functions fitting the hard bounds, we would be directly assigning these points a score based on how well their optimal functions fit a probability distribution.

Boundary layers that are too narrow can result in ill-posed problems related to non-uniqueness of the optimal solution. They can also result in computational difficulties, because in the vicinity of the optimal solution close functions *f* can provide very close values of *I*—as the gradients of the objective functions are too small (i.e. the top of the function in figure 1 is too flat). To transform an ill-posed problem into a well-posed problem, we use the standard framework based on Tikhonov regularization (see Tikhonov & Arsenin [21]) for a general introduction into Tikhonov's regularization method). In particular, we insert into the objective functional a function R of the highest unrestricted derivative, which is multiplied by a small parameter *I* too strongly. In our computation, we choose the following equation for *R*:
*C*/2 (indeed, this is how Tikhonov regularization works), so it should give the mid-point in the range of the variation of *f*^{(p + 1)}. In the case where we do not have information on this derivative, it is natural to consider *C* to be zero, as we assume in this paper. However, we should stress that variation of this parameter within a broad range does not strongly affect the results when *η* is sufficiently small. The parameter *σ* determines the width of the range of *f*^{(p + 1)} which is allowed. In this paper, we set *σ* = 1 for the sake of simplicity.

In order to represent the problem of finding the function maximizing (2.5) as an optimal control problem, we consider the trivial differential equation linking *f* and each of its derivatives in turn. We treat the highest derivative considered, *u*(*x*), after being integrated *p* + 1 times, yields the function that maximizes *I*? Details are provided in appendix A in the electronic supplementary material. Taking the notation *f*, and applying Pontryagin's maximum principle [22], we find that the solutions to the optimal control problem must satisfy the following boundary value problem
*x*_{min}, *x*_{max}] by the points *x*_{i}, *x*_{i}_{+1}], the derivatives in (2.4) which are fixed at *x*_{i} and *x _{i }*

_{+ 1}must serve as boundary conditions

*j*th derivative is therefore left unfixed. In this way, we obtain 2

*p*boundary conditions over each interval, and each of our boundary value problems is therefore well-posed.

### (d) Demonstration model

In order to demonstrate the method, we shall use it to conduct an investigation into structural sensitivity of the Rosenzweig–MacArthur predator–prey model [12,20]. This is a classical model in mathematical ecology, due both to its elegance and to the complex dynamics it can admit. The equations for the model are
*x* is the prey density and *z* is the predator density; *g* is the growth term (including natural mortality) for the prey species, whereas *f* is the functional response of the predator—the rate of prey consumption as a function of prey density. For our particular choice of growth term *g*, we choose the logistic function, *r* is the initial prey growth rate; *K* is the carrying capacity of the prey; *k* is the trophic efficiency coefficient; and *d* is the natural predator mortality rate.

We assume that the functional response *f* is unspecified, but we require it to satisfy the following restrictions:

Condition (2.12) signifies that negative feeding is not possible; (2.13) and (2.14) come from the assumption that *f* is a functional response with similar shape to a Holling type II function; therefore, it should be increasing—more prey should always result in a greater consumption rate—and should be saturating [23]—the predator gets diminishing returns for large prey numbers owing to the handling time; (2.15) is a natural requirement. To demonstrate that our method works for complex error bounds *f*′(*x*) and *x*) are, respectively,

Our aim is to check the sensitivity of the stability of the interior equilibrium in this partially specified system to the choice of functional response term. The model is well studied in the literature [14] with the key results here being that we have a single interior equilibrium only, this equilibrium is linearly stable if
*x** and *f*′(*x**), because (2.11) implies that *f*(*x**) = *d*/*k*. We need to project our space of valid functions onto this space by scanning the *x** and *f*′(*x**) values (we set *f′*(*x**) = *F*, *G* and *Q* are given by equation (2.6) (with *F* = *F*_{0}, *G* = *F*_{1}, *Q* = *F*_{2}), with *γ* = 50. *R* is defined as (2.7), with *C* = 0. As our local restrictions are at *x* = 0 (from (2.15)) and *x* = *x** (we fix *x**, *f*(*x**) and *f′*(*x**), because we are checking if they correspond to a valid function), we split the domain into [0, *x**] and [*x**, *x*_{max}], and apply the Pontryagin's maximum principle to obtain a boundary value problem over each of these domains, as per §2c. More details of how we apply the framework outlined in §2c here, along with the equations and conditions for the boundary value problem, are contained in appendix B in the electronic supplementary material.

The resultant nonlinear boundary value problem is not easily investigable analytically, but several numerical techniques exist for the solution of such problems. The method we apply here is the nonlinear shooting method, a standard technique that is covered in many textbooks on numerical analysis, for instance [24]. In this way, we can compute the function *I*—so

## 3. Implementation of the method

### (a) Stability plots

For the Rosenzweig–MacArthur model with the function bounds shown in figure 2, and parameter values *d* = 0.1, *k* = 0.3 and *K* = 0.58, the region of generalized bifurcation space corresponding to valid functions is shown in figure 3. In figure 3, dark blue indicates that the optimal function which takes these *x** and *a* and *b* as far as the function stays within the constraints. Note that the Monod function is very popular in ecological literature [23].

From figure 3, it is immediately clear that there is a high degree of structural sensitivity in the system—the presence of both green and red regions in our projected domain indicates that the system can admit both a stable and an unstable interior equilibrium when different valid functional responses are chosen. It is also clear that our method outperforms the approach of varying parameters of the fixed Monod function, because our projected domain covers the entirety of the light blue region, and more besides. The reason for this is that varying the parameters of the Monod function, or of any other fixed function, can only perturb the function in an extremely constrained way, which artificially restricts the range of values *x** and *f*′(*x**) that can be taken by such a function while staying between the upper and lower bounds. Fixing the function to a particular equation means that local perturbations cannot be made without inducing global perturbations elsewhere. For example, with the Monod function, the derivative at some *x** could be increased by decreasing the half saturation constant, *b*. However, because we have fixed our function to a precise equation, this decrease in *b* will necessarily increase the slope of the function at the origin, and might take the function above the upper bounds. Similarly, increasing the maximal feeding rate *a* may take the function above the upper bound at high values of *x*. In the case of figure 3, because the light blue region is entirely contained in the stable region, all Monod functions fitting our data range will yield a stable interior equilibrium and varying its parameters will not detect the structural sensitivity in the system.

### (b) Probabilistic bifurcation analysis

Based on the stability plot shown in figure 3 that demonstrates the existence of structural sensitivity in the system, one can quantify this sensitivity by introducing measures of the probability of having different model behaviour, i.e. a stable and an unstable equilibrium. Introducing a probability density function on our generalized bifurcation space also allows us to follow changes in the probability of having different dynamics with variation of other models parameters. In the simplest case, we can consider a uniform probability distribution of functions in that case the probability of having a certain type of model behaviour will be given by the area/volume of the corresponding proportion of the domain in generalized bifurcation space taken by valid functions. For instance, in our investigation, the probability of having a stable equilibrium would be computed from figure 3 by taking the ratio between the area of the green region and the area of the red and green regions combined.

One of the salient features of the Rosenzweig–MacArthur model is a loss of stability of the interior equilibrium through a Hopf bifurcation as the carrying capacity, *K*, of the prey species is increased—a phenomenon known as the ‘paradox of enrichment’ [2,13]. Therefore, it is worthwhile to consider how the probability of having a stable equilibrium in our partially specified model varies with *K* (fixing all the other parameters for simplicity). This is given by the blue curve in figure 4. We see that for low carrying capacities, the probability of stability is 1, as in this case, the entire projected domain in generalized bifurcation space yields stable dynamics. As *K* is increased, we then see the possibility arises of having an unstable equilibrium, as the Hopf bifurcation curve enters our projected domain as is the case in figure 3. A further increase in *K* causes the red region of the projected domain to advance until it fills the entire domain, at which point all valid functional responses will yield an unstable interior equilibrium. Therefore, we can still observe the paradox of enrichment in the partially specified Rosenzweig–MacArthur model, but as a monotone decrease in the probability of observing a stable equilibrium, rather than a concrete stability loss, as would be the case in the standard Rosenzweig–MacArthur model, with all model functions fixed. The bifurcation taking place via a shift in probability reflects the fact that we retain the uncertainty in the formulation of the functional response throughout our analysis.

### (c) The degree of structural sensitivity

To quantify the structural sensitivity in the model, one can evaluate the ‘degree of structural sensitivity’ which is defined as twice the probability of two randomly chosen functions (according to the probability distribution we are considering in the generalized bifurcation space) yielding conflicting model behaviour [8]. In the case that we are interested in equilibrium stability, and considering a uniform distribution over the points in the projected domain for simplicity of notation, this will be given by
*V*_{stable} is the volume (or area) of the stable region of the projected domain, and *V*_{total} is the volume of the total projected domain. In §3d, we shall consider how to improve the degree of structural sensitivity by weighting points according to how well the corresponding model functions can fit the constraints, instead of assuming a uniform distribution. Note that the degree of sensitivity (i) will equal 0 if the stable region takes up either all or none of the projected domain; (ii) has a maximum of 1, which is attained whenever the stable region takes up exactly half of the projected domain; (iii) essentially only depends on *V*_{stable}/*V*_{total}, the probability of having a stable equilibrium; and (iv) will be unaltered if we replace the volume of the stable region with the volume of the unstable region in the calculation, because we necessarily have
*K* as the red curve in figure 4. Note that the degree of sensitivity will also depend on the other model parameters that are fixed in figure 4.

### (d) Weighting functions in the generalized bifurcation space

Quantification of the degree of sensitivity of a partially specified model may depend strongly on the choice of the probability density function assumed in the generalized bifurcation space. The assumption of a uniform probability distribution in (3.2) is somewhat hard to justify, so we should introduce a probability distribution *ρ* that *weights* points *x*-values) fraction of the possible values that has some valid function taking the values *x**.

We can also use optimal control theory to efficiently compute the functional density of points. We first partition the domain by points *ν*_{l}, *l* = 1, … , *p*, at which we check the range of values *f*(*ν*_{l}) can take. For each point *ν*_{l}, and for a given value *f*_{l}, we can check whether *f*_{l} can be taken by a valid function satisfying *ν*_{l} in the partitioning set along with 0, *x** and *x*_{max}, and specify *ν*_{l}, *f*_{l}) that best fits our quantitative and qualitative constraints. If *f*_{l}, we can work out the entire range of points such a function can pass through at *ν*_{l} and take the ratio of this range to the total height of the bound at *ν*_{l}. After we do this for all *ν*_{l}, we can then take the average of these values to approximate the functional density.

Figure 5 shows the area that can be covered for two sets of values, *a*, these values are *x** = 0.15, *b*, the values are *x** = 0.16 and *ρ* in (3.4).

## 4. Discussion and conclusion

In this paper, we have presented a method for detecting and quantifying structural sensitivity in biological models by formulating them as partially specified models and applying optimal control theory. This method has an advantage over standard parameter-based approaches [7], because it allows us to cover all function relations and not stick to any particular mathematical formulation. Interestingly, even more advanced non-parametric methods in statistics, which take into account constraints on functions, might not be reliable, because using two different non-parametric methods can result in rather different predictions [25]. Partially specified models work by leaving uncertain functions unspecified apart from some local and global qualitative constraints and some error bounds that the functions must pass between. In order to analyse such models in terms of the number of equilibria and their stability, for instance, we can then simply find the isocline equations and the Jacobian matrix as usual: any values that we need for this which are unknown owing to the function being unspecified can be considered as parameters in a generalized bifurcation space.

In partially specified models, the constraints on the considered functions are an integral part of the model, and not something to be explored as a supplementary investigation, after the analysis is done. The main approach is to project the set of functions fitting the constraints into the generalized bifurcation space. This is in contrast to the well-known framework of generalized modelling and the analogous structural kinetic modelling [16,18,26]. In these frameworks, unknown model functions are also left unspecified and local bifurcation analysis carried out by way of incorporating unknown values into a generalized bifurcation space. However, instead of incorporating the constraints on the functions into the model, and anchoring our generalized bifurcation space to these constraints via a projection, the initial model is transformed, so that the generalized bifurcation analysis comes out in terms of general parameters which are more biologically interpretable and measurable than equilibrium densities, etc. However, there is no escaping the fact that all of these generalized parameters will necessarily be values that need to be measured at the equilibrium density itself. With the approach proposed here, considering the whole possible span of functions, this dependence on measurements at a single density vanishes. It is interesting to briefly compare our sensitivity analysis results to those using the generalized modelling approach (see appendix C for details in the electronic supplementary material). We can see that the generalized modelling approach considers a huge range of values of the generalized parameters *x** and

The key innovation introduced in this paper concerns how we project the set of valid functions—those satisfying qualitative constraints and staying within certain bounds—into the generalized bifurcation space. To find such a projection in general, we need to determine whether or not a differentiable function *f* exists satisfying some global constraints on it and its derivatives (2.2) and (2.3), as well as being fixed at some points (2.4). The precise solution of this problem in general is an extremely challenging mathematical endeavour. Here, we construct an approximation to the initial mathematical problem by constructing a functional that penalizes functions for breaking the constraints or leaving the area between the bounds, finding the function that maximizes it and checking whether this function satisfies the constraints. This is an optimization problem that can then be solved by standard numerical methods.

In [8], the projection of the set of valid functions into the generalized bifurcation space was done exactly, using a geometrical method. This approach is computationally very quick, but requires new calculations by hand for each set of qualitative constraints on functions, and is only valid for error bounds that satisfy these constraints themselves. The optimal control approach presented here is more flexible with regards to the shape of error bounds and qualitative constraints it can handle. In particular, changing the error bounds requires only a small change in the optimal control problem, because the actual shape of the boundary only determines the precise functional we aim to optimize. Because the optimal control approach distils all of the functionals inputted into a boundary value problem to be solved, it opens up the possibility for software to be developed that performs structural sensitivity analysis without a need for extra analytical work by the user.

We should stress here that when building the objective functional (2.5) we ‘soften’ the bounds *h ^{i}*(

*x*) in (2.2) and (2.3): functions and their derivatives no longer need to lie between two curves to be valid. However, the use of smoothed step functions in the objective functional mean that we approximate the hard bounds. The framework presented here could be altered by dispensing with bounds completely and replacing the smoothed step functions in the objective functional with Gaussian-type functions. The generalized bifurcation parameters could then be weighted by how the corresponding optimal functions score in the objective functional. This score is obtained automatically through the optimal control approach, and even in the case of approximately hard bounds, it may serve as an alternative weighting to the functional density introduced in §3d without any extra computation required. However, this involves giving full voice to the best-fitting function, whereas the functional density considers the range of possible functions taking the given generalized bifurcation values.

A potential misgiving of the approach to sensitivity analysis presented here is that we have demonstrated it on a very simple model: the two-component Rosenzweig–MacArthur predator–prey model, where only the functional response is unknown. However, we should stress that we have presented this model for the sake of simplicity. The same method can certainly be applied when multiple functions are unspecified, although we need to consider higher-dimensional domains instead of considering two-dimensional plots such as figure 3. Such domains would be harder to visualize, but we can still readily compute the probability of having a given type of dynamics by computing the (higher dimensional) volumes of the corresponding regions. Moreover, even if the unknown functions are embedded in complex models containing dozens or hundreds of equations, the analysis presented here stays the same.

There are several existing issues with regards to smooth numerical execution of the method presented here. These issues are related to the difficulty of solving high-dimensional boundary value problems numerically. In particular, more advanced numerical methods would need to be developed for efficiently solving these problems: here we have used the simple shooting method [24], which requires reasonable initial guesses as inputs, and shows slow convergence near the margins of the domain of valid functions in the generalized bifurcation space. Creation of more efficient methods would pay off remarkably, because it would enable the development of software to project unknown functions into generalized bifurcation spaces automatically, without any analytical work needed from the user. This would pave the way towards rigorous, flexible and fully automated structural sensitivity analysis being made available to the entire modelling community.

## Ethics

This work did not involve any active collection of any biological data.

## Data accessibility

This work does not have any experimental data.

## Authors' contributions

A.M, M.A. and O.K. conceived the sensitivity analysis method, interpreted the computational results. A.M. and M.A. wrote the paper. M.A. implemented and performed most of simulations and optimization calculations in consultation with A.M. and O.K. All authors gave final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

This research was supported by the Russian Foundation for Basic Research (RFBR grant no. 13-01-12452 ofi_m2).

## Acknowledgements

We highly appreciate the efforts of three anonymous reviewers for their extensive comments that helped us to substantially improve the manuscript.

- Received September 5, 2015.
- Accepted August 18, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.