# Autonomously and Simultaneously Refining Deep Neural Network Parameters by a
Bi-Generative Adversarial Network Aided Genetic Algorithm

Yantao Lu, Burak Kakillioglu, Senem Velipasalar
Department of Electrical Engineering and Computer Science
Syracuse University
Syracuse NY, 13244-1240
e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>This work has been
funded in part by National Science Foundation (NSF) under CAREER grant
CNS-1206291 and NSF Grant CNS-1302559.

###### Abstract

The choice of parameters, and the design of the network architecture are
important factors affecting the performance of deep neural networks. Genetic
Algorithms (GA) have been used before to determine parameters of a network.
Yet, GAs perform a finite search over a discrete set of pre-defined
candidates, and cannot, in general, generate unseen configurations. In this
paper, to move from exploration to exploitation, we propose a novel and
systematic method that autonomously and simultaneously optimizes multiple
parameters of any deep neural network by using a GA aided by a bi-generative
adversarial network (Bi-GAN). The proposed Bi-GAN allows the autonomous
exploitation and choice of the number of neurons, for fully-connected layers,
and number of filters, for convolutional layers, from a large range of values.
Our proposed Bi-GAN involves two generators, and two different models compete
and improve each other progressively with a GAN-based strategy to optimize the
networks during GA evolution. Our proposed approach can be used to
autonomously refine the number of convolutional layers and dense layers,
number and size of kernels, and the number of neurons for the dense layers;
choose the type of the activation function; and decide whether to use dropout
and batch normalization or not, to improve the accuracy of different deep
neural network architectures. Without loss of generality, the proposed method
has been tested with the ModelNet database, and compared with the 3D Shapenets
and two GA-only methods. The results show that the presented approach can
simultaneously and successfully optimize multiple neural network parameters,
and achieve higher accuracy even with shallower networks.

##  I Introduction

Deep learning-based techniques have found widespread use in machine learning.
Even before the convolutional approaches becoming popular, simple multi-layer
perceptron neural networks (MLPNN) were widely used for classification tasks
for several reasons. First, they are very easy to construct and run since each
layer is represented and operated as a single matrix multiplication. Second, a
neural network can become a complex non-linear mapping between input and
output by the introduction of non-linear activation functions after each
layer. Regardless of how large the input size (i.e. the number of features) or
the output size (i.e. the number of classes) are, a neural network can
discover relations between them when the network is sufficiently large; and
enough samples, which cover the problem domain as much as possible, are
provided during training.

Although a MLPNN is successful to form a complex non-linear relationship
between the feature space and the output class space, it lacks the ability of
discovering features by itself. Until recent years, the classical approach was
to provide either the input data directly or some high-level descriptors,
which are extracted by applying some algorithm on the input data, as the
source of features. Using the raw input as features does not guarantee to
yield any satisfactory mappings and the latter approach would require the
investigation of multiple hand-crafted descriptor extraction algorithms for
different applications. The introduction of convolutional layers in neural
networks removed the necessity of having prior feature extractors, which are
not easy to craft for different applications. Convolutional layers are
designed to extract features directly from the input. Since they have been
proven to be successful feature extractors and thanks to much faster
computation of the operations of convolutional neural networks (CNNs) on
specialized processors such as GPUs, the use of CNNs exploded recently. After
Krizhevksy et al. [1] achieved a significant increase in the classification
accuracy on the ImageNet Large Scale Visual Recognition Challenge (ILSVRC) [2]
in 2012, many others followed their approach, by creating different
architectures and applying them to numerous applications in many different
domains.

It is well-known that the training of deep learning methods requires large
amounts of data, and they usually perform better when training data size is
increased. However, for some applications, it is not always possible to obtain
more data when the dataset at hand is not large enough. In many cases, even
though the raw data can be collected easily, the labeling or annotation of the
data is difficult, expensive and time consuming. Successors of [1] yielded
better accuracies with less number of parameters on the same benchmark with
some architectural modifications using the same building blocks. This shows
that the choice of parameters, and the design of the architecture are
important factors affecting the performance. In fact, the design of a CNN
model is very important to achieve better results, and many researchers have
been working hard to find better CNN architectures [3, 4, 5, 6, 7, 8, 9] to
achieve higher accuracy.

However, there has not been much work on developing an established and
systematic way of building the structure of a neural network, and this task
heavily depends on trial and error, empirical results, and the designer’s
experience. Considering that there are many design and parameter choices, such
as the number of layers, number of neurons in each layer, number of filters at
each layer, the type activation function, the choice of using drop out or not
and so on, it is not possible to cover every possibility, and it is very hard
to find the optimal structure. In fact, often times some common settings are
used without even trying different ones. Moreover, the hyper-parameters in
training phase also play important role on how well the model will perform.
Likewise, these parameters are also tuned manually in an empirical way most of
the time.

In this work, we focus on optimizing the network architecture and its
different parameters for any neural network model. We propose a novel and
systematic way, which employs a revised generative adversarial networks,
referred to as Bi-GAN, together with a Genetic Algorithm (GA). The proposed
method can autonomously refine the number of convolutional layers, the number
and size of filters, number of dense layers, and number of neurons; decide
whether to use batch normalization and max pooling; choose the type of the
activation function; and decide whether to use dropout or not.

###  I-A Related Work

There have been works focusing on optimizing neural network architectures.
Most of the proposed approaches are based on the GAs [10], or evolutionary
algorithms, which are heuristic search algorithms. Benardos and Vosniakos [11]
proposed a methodology for determining the best neural network architecture
based on the use of a GA and a criterion that quantifies the performance and
the complexity of a network. In their work, they focus on optimizing four
architecture decisions, which are the number of layers, the number of neurons
in each layer, the activation function in each layer, and the optimization
function. Magnier and Haghighat [12] presented an optimization methodology
based on a combination of a neural network and a multi-objective evolutionary
algorithm. The methodology was used for the optimization of thermal comfort
and energy consumption in a residential house. Leung et al. [13] presented the
tuning of the structure and parameters of a neural network using an improved
GA. Ritchie et al. [14] proposed a method to automate neural network
architecture design process for a given dataset by using genetic programming.
Islam et al. [15] employed a genetic algorithm for finding the optimal number
of neurons in the input and hidden layers. They apply their approach to power
load prediction task and report better results than a manually designed neural
network. Yet, their approach is used to optimize only the number of neurons
for input and hidden layers, and optimization of other important design
decisions such as the number of layers or activation function type are not
discussed. Stanley and Miikkulainen [16] presented the NEAT algorithm for
optimizing neural networks by evolving topologies and weights of relatively
small recurrent networks. In a recent work, Miikkulainen et al. [17] proposed
CoDeepNEAT algorithm for optimizing deep learning architectures through
evolution by extending existing neuroevolution methods to topology, components
and hyperparameters.

The genetic algorithm-based optimization uses a given set of blueprints and
models, i.e. it performs a finite search over a discrete set of candidates.
Thus, GAs, in general, cannot generate unseen configurations, and they can
only make a combination of preset parameters. GAs are good at searching better
solutions from limited possibilities, such as type of layers and activation
functions. However, it cannot search for a solution which is not defined
before. In addition, the complexity of GAs increases significantly when the
number of choices increases to large scale. Rylander [18] has shown that the
generations needed for convergence increases exponentially with the node size.

Apart from the genetic algorithms, Bergstra and Bengio [19] have proposed
random search for hyper-parameter optimization, and stated that randomly
chosen trials are more efficient for hyperparameter optimization than trials
on a grid. Yan and Zhang [20] optimized architectures’ width and height with
growing running time budget through submodularity and suparmodularity.

Figure 1: Proposed Bi-GAN aided GA network for refining deep neural network
parameters.

Generative Adversarial Networks (GANs) [21] are one of the important
milestones in deep learning research. In contrast to CNNs, which extract rich
and dense representations of the source domain, and may eventually map source
instances into some classes, GANs generate instances of the source domain from
small noise. They employ deconvolution operators, or transposed convolutions,
to generate N-D instances from 1-D noise. GAN’s power comes from the
competition with the discriminator, which decides whether the generated
instance belongs to the source domain. Discriminator acts like the police who
is trying to intercept counterfeit money, where in this case the generator is
the counterfeiter. Generator and discriminator are trained together until
discriminator cannot distinguish the generated instances from the instances in
the source domain. GANs have been adapted in many applications [22, 23, 24,
25, 26].

In this paper, to move from exploration to exploitation, we propose a novel
and systematic method that autonomously and simultaneously optimizes multiple
parameters of any deep neural network by using a GA aided by our proposed bi-
generative adversarial network (Bi-GAN). In contrast to traditional GANs, our
proposed Bi-GAN involves two generators, and two different models compete and
improve each other progressively with a GAN-based strategy to optimize the
networks during GA evolutions. The Bi-GAN allows the autonomous exploitation
and choice of the number of filters and number of neurons from a large range
of values. Our proposed approach can be used to autonomously refine the number
of convolutional layers, the number and size of filters, number of dense
layers, and number of neurons; decide whether to use batch normalization and
max pooling; choose the type of the activation function; and decide whether to
use dropout or not, to improve the accuracy of different deep neural network
architectures.

For this work, without loss of generality, we tested the performance of our
approach by using the ModelNet database, and compared it with the 3D Shapenets
and two GA-only methods. The results show that the presented approach can
simultaneously and successfully optimize multiple neural network parameters,
and achieve increased accuracy even with shallower networks. The rest of this
paper is organized as follows: The proposed method is described in Sec. II.
The experimental results are presented in Sec. III, and the paper is concluded
in Sec. IV.

##  II Proposed Method

The overall structure of the proposed method is shown in Fig. 1. It involves a
GA aided by a Bi-GAN. While Bi-GAN is used to set/optimize the parameters from
a large set covering a large range, GA is applied to make discrete decision
from a small set of choices. The decisions the GA makes include the choice of
the activation function; whether or not to use batch normalization, dropout
and max pooling; number of convolutional layers and dense layers; and the
kernel size of convolutional layers. Table III shows the set of parameters,
and the discrete set of values that they can take.

TABLE I: Parameter Choices Activation function | {“relu”, “leaky relu”, “sigm.”, “tanh”}
---|---
Batch norm. | {True, False}
Dropout | {True, False}
Max pooling | {True, False}
Num. of conv layers | {1,2,3}
Num. of dense layers | {1,2,3}
Kernel size | {3,5}

The parameter set of the network $i$ for the GA has the following form:

$\begin{split}p_{i}^{GA}=&[prm_{i}^{conv^{1}},prm_{i}^{conv^{2}}...,prm_{i}^{conv^{C}},\\\
&prm_{i}^{dense^{1}},prm_{i}^{dense^{2}}...,prm_{i}^{dense^{D}}],\\\
&i\in\\{1,2,...,n_{m}\\}\end{split}$ (1)

where $C$ and $D$ are the maximum number of possible convolutional layers and
dense layers, respectively, $n_{m}$ is the number of network models in the
population and

$\begin{split}prm_{i}^{conv^{j}}=&[1/0\text{ (conv. layer exists or
not)},\text{kernel size},\\\ &\text{activation func.},1/0\text{ (for batch
norm.)},\\\ &1/0\text{ (for max pooling)}],\\\ &i\in\\{1,2,...,n_{m}\\}\\\
&j\in\\{1,2,...,C\\}\end{split}$ (2)
$\begin{split}prm_{i}^{dense^{k}}=&[1/0\text{ (dense layer exists or not)},\\\
&\text{activation func.},1/0\text{ (for batch norm.)},\\\ &1/0\text{(for
dropout)}],\\\ &i\in\\{1,2,...,n_{m}\\}\\\ &k\in\\{1,2,...,D\\}.\end{split}$
(3)

For the discrete set of parameters given in Table III, the values of $C$ and
$D$ are $3$.

First, $n_{m}$-many models are randomly initialized for the first network
population $N^{1}$. This is done by choosing the values of convolutional and
dense parameters, in (2) and (3), randomly, from the possible choices. Then,
the number of filters for the convolutional layers, and the number of neurons
for the fully connected layers are determined by the Bi-GAN as will be
described in Section II-A.

Then, the $n_{m}$ models are trained, and evaluated to obtain their accuracy
scores. Based on the accuracy scores, the GA is applied as detailed in Section
II-B.

###  II-A Bi-GAN network

We propose a novel and modified generative adversarial network (GAN), referred
to as Bi-GAN, to find the optimal network parameters, that have a large range
of values. The proposed Bi-GAN network for refining different neural network
parameters is shown in Fig. 2. It is composed of a generative part, an
evaluation part and a discriminator. In contrast to traditional GANs, there
are two generators ($G_{1}$ and $G_{2}$), two evaluators ($E_{1}$ and
$E_{2}$), and one discriminator ($D$). The input to the two generators is
Gaussian noise $z\sim p_{noise}(z)$. On the other hand, the input to the
evaluators is the training data $x\sim p_{data}(x)$.

Figure 2: Proposed Bi-GAN incorporating two generators and one discriminator.

As will be discussed in more detail below, the generators $G_{1}$ and $G_{2}$
have the same network structure. From input noise $p_{noise}(z)$, $G_{1}$ and
$G_{2}$ generate the input network parameters $G_{1}(z)$ and $G_{2}(z)$ to be
used and evaluated by $E_{1}$ and $E_{2}$, respectively. $E_{1}$ and $E_{2}$
have the structure of the neural network whose parameters are being optimized
or refined. They calculate the classification accuracy on the training data
$x$. $E_{i}(x,G_{i}(z)),i\epsilon\\{1,2\\}$ represents the classification
accuracy obtained by the evaluator $E_{i}$ when the parameters $G_{i}(z)$ are
used. The generator resulting in higher accuracy is marked as more accurate
generator $G_{a}$, and the other generator is marked as $G_{b}$, where
$a\in{1,2},b=!a$.

We define the discriminator $D$ as a network, which is used for binary
classification between better and worse generator. $G_{1}(z)$ and $G_{2}(z)$
are fed into the discriminator $D$, and the ground truth label about which is
the better generator comes from the evaluators. The discriminator $D$ provides
the gradients to train the worse performing generator.

####  II-A1 Generative part

The two generators $G_{1}$ and $G_{2}$ have the same neural network structure
shown in Fig. 3. Their input is a Gaussian noise vector $z$, and their outputs
are $G_{1}(z)$ and $G_{2}(z)$. As seen in Fig. 3, generators are composed of
fully connected layers with leaky relu activations. At the output layer,
$tanh$ is employed so that $G_{i}^{j}(z)\in(-1,1)$, where
$j\in\\{1,2,...,length(G_{i}(z))\\}$ and $i\in\\{1,2\\}$. Then, the range of
$G_{i}(z)$ is changed from $(-1,1)$ to $(pm_{min}^{j},pm_{max}^{j})$ by using

$G_{i}(z)^{\prime}=[G_{i}(z)\times\frac{pm_{max}-pm_{min}}{2}+\frac{pm_{max}+pm_{min}}{2}].$
(4)

In (4), $pm_{max}$ and $pm_{min}$ are preset maxima and minima values, which
are defined empirically based on values that a certain parameter can take, so
that the value of the refined parameters can only change between $pm_{max}$
and $pm_{min}$. For instance, in the case of the number of neurons for a fully
connected layer, $pm_{max}$ and $pm_{min}$ are 4000 and 10, respectively. The
re-scaled values $G_{1}(z)^{\prime}$ and $G_{2}(z)^{\prime}$ are then used as
the parameters of evaluator networks. The length of $G_{i}(z)$ is determined
by the number of network parameters that are refined, and is set at the
generator network’s last fully connected layer.

Figure 3: Generator network

Generators are trained/improved by the discriminator, which is a binary
classifier used to differentiate the results from generator outputs $G_{1}(z)$
and $G_{2}(z)$. Labels “a” and “b” represent the generators with higher
accuracy and lower accuracy results, respectively. The generator, which has
the worse performance and is labeled by “b”, is trained by stochastic gradient
descent (SGD) from the discriminator to minimize $log(1-D(G_{b}(z))$ by using

$\bigtriangledown_{G_{b}}\frac{1}{m}\displaystyle\sum_{j=1}^{m}log(1-D(G_{b}(z^{(j)}))),$
(5)

where $m$ is the number of epochs.

When $G_{a}^{\prime}(z)$ becomes equal to $G_{b}^{\prime}(z)$ for two
consecutive iterations, the weights of $G_{b}$ will be re-initialized to
default random values. The purpose of this step is to prevent the optimization
stopping at a local maxima and also prevent the vanishing tanh gradient
problem.

####  II-A2 Evaluation part

As mentioned above, one of the strengths of the proposed approach is that it
can be used to refine/optimize parameters of different deep neural network
structures. In other words, the evaluator networks have the same structure as
the neural network whose parameters are being optimized or refined.

Evaluator networks are built by using the parameters $G_{1}(z)^{\prime}$ and
$G_{2}(z)^{\prime}$ provided by the generators. The training data $x\sim
p_{data}(x)$ is used to evaluate these network models. We employ an early
stopping criteria. More specifically, if no improvement is observed in $c$
epoches, the training is stopped.

We then obtain the accuracies $acc_{i}=\mathbb{E}_{x\sim
p_{data}(x)}E_{i}(x)$, $i=\\{1,2\\}$. Let $a$ be the value of $i$ resulting in
higher accuracy, and $b=!a$. Then “a” is used as the ground truth label for
the discriminator, which marks the generator with better parameters, and
trains the worse generator $G_{b}$.

####  II-A3 Discriminator

We define the discriminator $D$ as a network (seen in Fig. 4, whose output is
a scalar softmax output, which is used for binary classification between
better generator and worse generator. $G_{1}(z)$ and $G_{2}(z)$ are fed into
the discriminator $D$, and the ground truth label about which is better
generator comes from the evaluators. Let $D(G(z))$ represent the probability
that $G(z)$ came from the more accurate generator $G_{a}$ rather than $G_{b}$.
We train $D$ to maximize the probability of assigning the correct label to the
outputs $G_{1}(z)$ and $G_{2}(z)$ of both generators. Moreover, we
simultaneously train the worse generator $G_{b}$ to minimize
$log(1-D(G_{b}(z))$. The whole process can be expressed by:

$min_{G_{a}}max_{D}\mathbb{E}_{z\sim
p_{z}(z)}(log(D(G_{a}(z)))+log(1-D(G_{b}(z)))),$ (6)

where, $a=argmax_{i=\\{1,2\\}}(\mathbb{E}_{x\sim p_{data}(x)}E_{i}(x))$,
$b=!a$.

Figure 4: Discriminator network

The pseudo code for the proposed Bi-GAN is provided in Algorithm 1.

while  _in the iterations_ do

Generate $m\times 2$ noise samples
$\\{Z_{1}^{(1)},Z_{1}^{(2)},...,Z_{1}^{(m)}\\}$ and
$\\{Z_{2}^{(1)},Z_{2}^{(2)},...,Z_{2}^{(m)}\\}$ from Gaussian white noise

while  _j in range(m)_ do

Build evaluators $E_{1}^{(j)}$ and $E_{2}^{(j)}$ based on parameters from
$G_{1}(Z_{1}^{(j)})$ and $G_{2}(Z_{2}^{(j)})$

Calculate $acc_{1}^{j}$ and $acc_{2}^{j}$ from $\mathbb{E}_{x\sim
p_{data}(x)}E_{i}(x)$

End if no acc. impr. after $c$ epoches

end while

Calculate mean value $acc_{i}=(1/m)\sum_{j=1}^{m}acc_{i}^{j}$, $i=\\{1,2\\}$

Find $G_{a}$ as $G_{argmax(acc_{1},acc_{2})}$ and $G_{b}$ as the other one.

Update Discriminator by SGD:
$\bigtriangledown_{D}\frac{1}{m}\displaystyle\sum_{j=1}^{m}(log(D(G_{a}(z^{(j)})))+log(1-D(G_{b}(z^{(j)}))))$

Update Generator $G_{b}$ by SGD:
$\bigtriangledown_{G_{b}}\frac{1}{m}\displaystyle\sum_{j=1}^{m}log(1-D(G_{b}(z^{(j)})))$

end while

Algorithm 1 Bi-GAN Algorithm.

###  II-B Genetic Algorithm

As mentioned above, we use a GA to make discrete decisions from a set of
choices shown in Table III. Within each GA evolution, our proposed Bi-GAN is
used to set/optimize the values of the number of filters for the convolutional
layers and the number of neurons for the dense layers. Then, the network
models are trained and evaluated to obtain their accuracy scores. Based on the
accuracy scores, the GA is applied.

####  II-B1 Initial Population

The first generation of the networks, $N^{1}$, is generated randomly such that
$N^{1}=\\{N_{1},N_{2},...,N_{n_{m}}\\}$, where $n_{m}$ is the number of
models. This is done by choosing the values of convolutional and dense
parameters, in (2) and (3), randomly, from the possible choices.

####  II-B2 Bi-GAN optimization

Our proposed Bi-GAN is used, as described in Sec. II-A, to update the number
of neurons for the fully connected layers, and the number of filters for the
convolutional layers, of the $n_{m}$ network models.

####  II-B3 Evaluation

After the number of neurons for the fully connected layers, and the number of
filters for convolutional layers are determined by Bi-GAN, each generated
network model $N_{i}$ ($i\in\\{1,...,n_{m}\\})$ will be evaluated by the
fitness function $fitness=f(N_{i})$, which is a measure of the accuracy of
each model. Models with better performance will have higher values. Thus,
$E=\\{E_{1},E_{2},...,E_{n_{m}}\\}$, will hold the fitness scores
$E_{i}=f(N_{i})$, where $i\in\\{1,2,...,n_{m}\\}$.

####  II-B4 Selection

In the selection part, $t$-many top ranked models are selected from the
$sorted(E)$ and $r$-many models are selected randomly from the rest of the
network models. Then, $d$-many models are dropped in order to prevent over-
fitting and getting stuck at a local optimum. The remaining selected models
are the parent models ($P$), which will be used to create new models for the
next generation.

####  II-B5 Crossover and Mutation

Crossover is applied to generate $n_{m}$-many child network models from the
parents. The choice of parents is performed as follows: Instead of always
choosing two parents randomly from the parent pool, we associate a counter
$C_{P}$ with each parent $P$, and initialize it to zero. This counter is
incremented by one each time a parent is used for crossover. First, two
parents are selected randomly from the $t+r-d$ many parents. A new ‘child’
network is generated from the parents via crossover, and the counters of the
parents are incremented by one. Then, two parents, whose counter is still
zero, are selected randomly from the parent pool. Another network is generated
from them via crossover, and the counters of the parents are incremented. If
there is only one network model left with counter equal to zero, and the
number of children is still less than $n_{m}$, then this model is chosen as
one of the parents, and the other parent is chosen randomly from the rest of
the models who have a counter value of one. If there are no more parents left
with counter equal to zero, and the number of children is still less than
$n_{m}$, then two parents, whose counter is one, are picked randomly, and
their counter is incremented to two after crossover. This process is repeated
until the number of children models reaches $n_{m}$.

The crossover between parent models $a$ and $b$ is performed as illustrated in
Fig. 5. First two integers (ID1 and ID2) are picked randomly between $1\text{
and }C$ and $1\text{ and }D$, respectively. Then, the parameters of the child
network is set so that

$\begin{split}p_{child}^{GA}=&[prm_{a}^{conv^{1}},prm_{a}^{conv^{2}},...,prm_{a}^{conv^{ID_{1}}},\\\
&prm_{b}^{conv^{ID_{1}+1}},...,prm_{b}^{conv^{C}},\\\
&prm_{a}^{dense^{1}},prm_{a}^{dense^{2}},...,prm_{a}^{dense^{ID_{2}}},\\\
&prm_{b}^{dense^{ID_{2}+1}},...,prm_{b}^{dense^{D}}].\end{split}$ (7)

After all the $n_{m}$-many child networks are obtained via crossover, 20% of
the population is chosen randomly to perform mutation. As seen in (2) and (3),
there are five different convolutional layer parameters, and four different
dense layer parameters. Thus, there are $5*C+4*D$-many possible parameters
that can be mutated. An integer is picked randomly between $1$ and $5*C+4*D$,
and the corresponding parameter type is chosen randomly from the possible
choices in Table III. For instance, if the random number corresponds to the
filter size parameter, then its value is chosen randomly from $\\{3,5\\}$.

Then, the entire process is repeated by using this new population, updating
the number of neurons and the number of filters for each network model in the
population by using our propose Bi-GAN, and so on. The pseudo code for the
entire process is provided in Algorithm 2.

Figure 5: Crossover.

Randomly initialize $n_{m}$ models for population $N^{1}$.

while  _$i^{th}$ iteration_ do

Update hyper parameters by Bi-GAN (Alg. 1).

Train and evaluate $N^{i}_{1},N^{i}_{2},...,N^{i}_{n_{m}}$ by fitness function
$f(N^{i}_{j})$ and obtain scores $E$.

Select $t$ top scored networks $N_{top}=N^{i}(argmax(E))$

Randomly choose $r$ networks $N_{rand}$ from the rest of population $N^{i}$

Merge $N_{top}$ and $N_{rand}$ and then drop $d$ networks ($N_{drop}$)

Form $N_{parent}={(N_{top}\bigcup N_{rand})-N_{drop}}$

Choose parents from $N_{parent}$ for crossover and generate $n_{m}$ new
networks and add them to $N^{i+1}$

Choose $20\%$ of the networks in $N^{i+1}$, and perform mutation on them.

end while

Algorithm 2 GA-BiGAN algorithm.

##  III Experimental Results

Without loss of generality, we have applied the proposed approach on a 3D
convolutional neural network by using the voxelized version of ModelNet40
dataset, which contains 3D CAD models of 40 object classes. Wu et al. [27]
voxelized each object from the ModelNet at 12 different orientations (around
gravity axis) for data augmentation. Voxel grids are $30\times 30\times 30$,
and every object is fitted into this range. We used this pre-voxelized version
of the dataset in our experiments. The dataset contains 40 subfolders for
different objects. Each of these 40 subfolders contains 2 subfolders for
training and testing. Train:test ratio differs for each object, but overall
train:test ratio is around 3:1. Some example voxelized objects from the
ModelNet40 dataset are shown in Fig. 6.

Figure 6: Sample voxelized objects from ModelNet40 dataset.

We have compared our proposed approach with two other approaches, which are
based only on GAs and referred to as small-set GA, and large-set GA. The
small-set GA is a basic genetic algorithm with limited number of choices.
Large-set GA is given a larger set of choices for the number of neurons, and
the number of filters. As for the activation function, batch normalization,
dropout and max pooling decisions, the number of convolutional and dense
layers, and the kernel size, the parameter choices are the same as in Table
III for both small- and large-set GA. The difference between small- and large-
set GA is the parameter choices for the number of neurons and the number of
filters. For the small-set GA these choices are as follows:

Num. of neurons: | {16,32,64,128,256,512,1024,2048,4096}
---|---
Num. of kernels: | {1,4,16,64,256}

For the large-set GA, the number of neurons can be any integer between $16$
and $4096$, and the number of filters can be any integer between $1$ and
$256$. In our proposed approach, the number of neurons and the number of
filters are determined by our proposed Bi-GAN method from this continuous
range.

These three approaches were run by using the same data, for the same amount of
time to compare their performances. The parameters used in Algorithm 1 and 2
are as follows: For the Bi-GAN part, $m=100$, and $c=5$. For the GA part, the
parameters used are $n_{m}=25$, $t=4$, $r=2$ and $d=1$. The results are
summarized in Table II. Same population size was used for all the GAs. As can
be seen, our proposed approach provides the highest accuracy, and performs
better compared to only GA-based approaches. Figure 7 shows the accuracy of
each method over time. The proposed method determines the number of neurons
and the number of filters without requiring a discrete set of choices.

TABLE II: Accuracy values obtained with different networks | Accuracy
---|---
ShapeNet model | 0.8417
Small-set GA | 0.8294
Large-set GA | 0.3641
Proposed Method | 0.8520

Figure 7: Accuracy of different network refinement approaches over time.

Another important point is the performance comparison with respect to
Shapenets [27], which is a hand-crafted network with four convolutional
layers, and two dense layers. Table II also includes the accuracy obtained
when the Shapenet [27] model is used. As can be seen, by autonomously refining
the network parameters, our proposed approach provides higher accuracy than
the manually designed model by using a shallower network. The parameter values
that each of the three approaches (proposed method, small-set GA and large-set
GA) ends up choosing/using are provided in Table IV (the first entry is the
filter size). As can be seen, the proposed method achieves this higher
accuracy of 85.2% by using only three convolutional layers as opposed to the
four-convolutional-layer Shapenet model.

Iterations | 5-pops | 10-pops | 20-pops
---|---|---|---
| Small-set GA | Proposed | Small-set GA | Proposed | Small-set GA | Proposed
10 | 0.5813 | 0.6077 | 0.5945 | 0.6563 | 0.6409 | 0.7945
20 | 0.6969 | 0.7601 | 0.6183 | 0.7798 | 0.7689 | 0.7984
50 | 0.7626 | 0.8124 | 0.8069 | 0.8197 | 0.8113 | 0.8429
TABLE III: Comparison of the proposed method with the Small-set GA for different population sizes. TABLE IV:  Approach | Final Parameter Values
---|---
Sm-set GA | [C1: 5,‘Relu’, NO Batch norm,Maxpool,64],
| [C2: 5,‘Relu’,No Batch,Max pool,64],
| [C3: 3,‘Relu’,0,1,256],
| [D1: Relu, No Batch, Dropout, 256],
| [D2: Leaky relu, no batch, no dropout, 64]
Lrg-set GA | [C1: 5,’Leaky Relu’,No Batch n., max pool,30],
| [C2: 3,‘Relu’,No Batch, Max pool, 31],
| [C3: 5,‘Relu’,No batch, No max pool, 9],
| [D1: Leaky Relu, Batch norm, Dropout, 31]
Prop.Meth. | [C1: 5,’Relu’,No Batch n,max pool,80],
| [C2: 3,‘Relu’,No Batch,Max pool, 105],
| [C3: 3,‘Leaky Relu’,No batch n.,Max pool,202],
| [D1: relu, no batch, no dropout, 601],
| [D2:leaky relu, no batch, dropout, 240]

Since small-set GA performs better than the large-set GA, in the remainder of
the experiments, we compared our proposed method with the small-set GA. We
have tried three more population sizes while keeping the other parameter
choices same as before. The results are summarized in Table III. As can be
seen, the proposed method provides the higher accuracy rates for all different
population sizes. Figures 8 and 9 show the change in accuracy and loss of
these approaches, respectively, with each evolution (when population size is
20). As can be seen, the proposed method performs better during its
evolutions.

Figure 8: Accuracy with every evolution of the proposed method and the Small-
set GA.

Figure 9: Loss with every evolution of the proposed method and the Small-set
GA.

We also performed one more experiment, where we let the maximum number of
convolutional layers to be 5. In this experiment, the proposed method achieved
an accuracy of 86.62%.

##  IV Conclusion

In this paper, we have presented a novel and systematic method that
autonomously and simultaneously optimizes multiple parameters of any given
deep neural network by using a genetic algorithm (GA) aided by a novel Bi-
Generative Adversarial Network (GAN) with two generators, which is referred to
as Bi-GAN. The proposed Bi-GAN allows the autonomous exploitation and choice
of the number of neurons, for the fully-connected layers, and number of
filters for the convolutional layers, from a large range of values. Our
proposed approach can be used to autonomously refine the number of
convolutional layers and dense layers, number and size of kernels, and the
number of neurons; choose the type of the activation function; and decide
whether to use dropout and batch normalization or not, to improve the accuracy
of different deep neural network architectures. Without loss of generality,
the proposed method has been tested with the ModelNet database, and compared
with the 3D Shapenets and two GA-only methods. The results show that the
presented approach can simultaneously and successfully optimize multiple
neural network parameters, and achieve increased accuracy even with shallower
networks.

## References

  * [1] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton,  “Imagenet classification with deep convolutional neural networks,”  in Advances in neural information processing systems, 2012, pp. 1097–1105.

  * [2] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei,  “Imagenet: A large-scale hierarchical image database,”  in Computer Vision and Pattern Recognition, 2009. CVPR 2009. IEEE Conference on. IEEE, 2009, pp. 248–255.
  * [3] Karen Simonyan and Andrew Zisserman,  “Very deep convolutional networks for large-scale image recognition,”  arXiv preprint arXiv:1409.1556, 2014.
  * [4] Matthew D Zeiler and Rob Fergus,  “Visualizing and understanding convolutional networks,”  in European conference on computer vision. Springer, 2014, pp. 818–833.
  * [5] Min Lin, Qiang Chen, and Shuicheng Yan,  “Network in network,”  arXiv preprint arXiv:1312.4400, 2013.
  * [6] Christian Szegedy, Wei Liu, Yangqing Jia, Pierre Sermanet, Scott Reed, Dragomir Anguelov, Dumitru Erhan, Vincent Vanhoucke, Andrew Rabinovich, et al.,  “Going deeper with convolutions,”  Cvpr, 2015.
  * [7] Christian Szegedy, Vincent Vanhoucke, Sergey Ioffe, Jon Shlens, and Zbigniew Wojna,  “Rethinking the inception architecture for computer vision,”  in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2016, pp. 2818–2826.
  * [8] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun,  “Deep residual learning for image recognition,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 770–778.
  * [9] Shaoqing Ren, Kaiming He, Ross Girshick, and Jian Sun,  “Faster r-cnn: Towards real-time object detection with region proposal networks,”  in Advances in neural information processing systems, 2015, pp. 91–99.
  * [10] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan,  “A fast and elitist multiobjective genetic algorithm: Nsga-ii,”  IEEE Transactions on Evolutionary Computation, vol. 6, no. 2, pp. 182–197, April 2002.
  * [11] PG Benardos and G-C Vosniakos,  “Optimizing feedforward artificial neural network architecture,”  Engineering Applications of Artificial Intelligence, vol. 20, no. 3, pp. 365–382, 2007.
  * [12] Laurent Magnier and Fariborz Haghighat,  “Multiobjective optimization of building design using trnsys simulations, genetic algorithm, and artificial neural network,”  vol. 45, pp. 739–746, 03 2010.
  * [13] F. H. F. Leung, H. K. Lam, S. H. Ling, and P. K. S. Tam,  “Tuning of the structure and parameters of a neural network using an improved genetic algorithm,”  IEEE Transactions on Neural Networks, vol. 14, no. 1, pp. 79–88, Jan 2003.
  * [14] Marylyn D Ritchie, Bill C White, Joel S Parker, Lance W Hahn, and Jason H Moore,  “Optimizationof neural network architecture using genetic programming improvesdetection and modeling of gene-gene interactions in studies of humandiseases,”  BMC bioinformatics, vol. 4, no. 1, pp. 28, 2003.
  * [15] Badar Ul Islam, Zuhairi Baharudin, Muhammad Qamar Raza, and Perumal Nallagownden,  “Optimization of neural network architecture using genetic algorithm for load forecasting,”  in Intelligent and Advanced Systems (ICIAS), 2014 5th International Conference on. IEEE, 2014, pp. 1–6.
  * [16] Kenneth O Stanley and Risto Miikkulainen,  “Evolving neural networks through augmenting topologies,”  Evolutionary computation, vol. 10, no. 2, pp. 99–127, 2002.
  * [17] Risto Miikkulainen, Jason Liang, Elliot Meyerson, Aditya Rawal, Dan Fink, Olivier Francon, Bala Raju, Arshak Navruzyan, Nigel Duffy, and Babak Hodjat,  “Evolving deep neural networks,”  arXiv preprint arXiv:1703.00548, 2017.
  * [18] Bart Ian Rylander,  Computational Complexity and the Genetic Algorithm,  Ph.D. thesis, Moscow, ID, USA, 2001,  AAI3022336.
  * [19] James Bergstra and Yoshua Bengio,  “Random search for hyper-parameter optimization,”  Journal of Machine Learning Research, vol. 13, no. Feb, pp. 281–305, 2012.
  * [20] Junqi Jin, Ziang Yan, Kun Fu, Nan Jiang, and Changshui Zhang,  “Neural network architecture optimization through submodularity and supermodularity,”  arXiv preprint arXiv:1609.00074, 2016.
  * [21] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio,  “Generative adversarial nets,”  in Advances in neural information processing systems, 2014, pp. 2672–2680.
  * [22] Leon A Gatys, Alexander S Ecker, and Matthias Bethge,  “A neural algorithm of artistic style,”  arXiv preprint arXiv:1508.06576, 2015.
  * [23] Alec Radford, Luke Metz, and Soumith Chintala,  “Unsupervised representation learning with deep convolutional generative adversarial networks,”  arXiv preprint arXiv:1511.06434, 2015.
  * [24] Justin Johnson, Alexandre Alahi, and Li Fei-Fei,  “Perceptual losses for real-time style transfer and super-resolution,”  in European Conference on Computer Vision, 2016.
  * [25] Phillip Isola, Jun-Yan Zhu, Tinghui Zhou, and Alexei A Efros,  “Image-to-image translation with conditional adversarial networks,”  arXiv preprint, 2017.
  * [26] Jun-Yan Zhu, Taesung Park, Phillip Isola, and Alexei A Efros,  “Unpaired image-to-image translation using cycle-consistent adversarial networks,”  arXiv preprint arXiv:1703.10593, 2017.
  * [27] Zhirong Wu, Shuran Song, Aditya Khosla, Fisher Yu, Linguang Zhang, Xiaoou Tang, and Jianxiong Xiao,  “3d shapenets: A deep representation for volumetric shapes,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2015, pp. 1912–1920.


# Structural sensitivity in the functional responses of predator-prey models

Sarah K. Wyse${}^{\text{a},\star}$, Maria M. Martignoni${}^{\text{b}}$, May
Anne Mata${}^{\text{c}}$,
Eric Foxall${}^{\text{a}}$, and Rebecca C. Tyson${}^{\text{a}}$

###### Abstract

In mathematical modeling, several different functional forms can often be used
to fit a data set equally well, especially if the data is sparse. In such
cases, these mathematically different but similar looking functional forms are
typically considered interchangeable. Recent work, however, shows that similar
functional responses may nonetheless result in significantly different
bifurcation points for the Rosenzweig-MacArthur predator-prey system. Since
the bifurcation behaviours include destabilising oscillations, predicting the
occurrence of such behaviours is clearly important. Ecologically, different
bifurcation behaviours mean that different predictions may be obtained from
the models. These predictions can range from stable coexistence to the
extinction of both species, so obtaining more accurate predictions is also
clearly important for conservationists. Mathematically, this difference in
bifurcation structure given similar functional responses is called structural
sensitivity. We extend the existing work to find that the Leslie-Gower-May
predator-prey system is also structurally sensitive to the functional
response. Using the Rosenzweig-MacArthur and Leslie-Gower-May models, we then
aim to determine if there is some way to obtain a functional description of
data so that different functional responses yield the same bifurcation
structure, i.e., we aim to describe data such that our model is not
structurally sensitive. We first add stochasticity to the functional responses
and find that better similarity of the resulting bifurcation structures is
achieved. Then, we analyze the functional responses using two different
methods to determine which part of each function contributes most to the
observed bifurcation behaviour. We find that prey densities around the
coexistence steady state are most important in defining the functional
response. Lastly, we propose a procedure for ecologists and mathematical
modelers to increase the accuracy of model predictions in predator-prey
systems.

${}^{\text{a}}$ Department of Computer Science, Mathematics, Physics and
Statistics, Irving K. Barber Faculty of Science, University of British
Columbia Okanagan, SCI354, 1177 Research Road, Kelowna V1V 1V7, BC, Canada
${}^{\text{b}}$ Department of Mathematics and Statistics, Memorial University
of Newfoundland, 230 Elizabeth Ave, St. John’s, A1C 5S7, NL, Canada
${}^{\text{c}}$ Department of Math, Physics, and Computer Science, University
of the Philippines Mindanao, Mintal, Davao City, Philippines
⋆ Corresponding author.
Email addresses<EMAIL_ADDRESS>(S.K. Wyse<EMAIL_ADDRESS>(M.M. Martignoni),
<EMAIL_ADDRESS>(M.A. Mata<EMAIL_ADDRESS>(E. Foxall),
<EMAIL_ADDRESS>(R.C. Tyson).
ORCID IDs: 0000-0002-5500-2711 (S.K. Wyse), 0000-0002-0192-8228 (M.M.
Martignoni), 0000-0002-2967-344X (M.A. Mata), 0000-0003-1610-6662 (E. Foxall),
0000-0002-7373-4473 (R.C. Tyson).

Keywords – Structural sensitivity, bifurcation structure, functional response,
coexistence steady state, stochasticity, prey density domain

Funding: This work was supported by The University of British Columbia
Okanagan [Undergraduate Research Award] and Natural Sciences and Engineering
Research Council of Canada [grant numbers RGPIN-2016-0577, RGPIN-2018-04480].
The funding sources were not involved in the conduction of research or
preparation of this article.

##  1 Introduction

Mathematical modeling of biological phenomena necessitates making assumptions
about the biological system so that the underlying mechanisms can be expressed
mathematically. Indeed, the modeler must decide which features of a system to
include and how to realistically formulate them (Krkošek et al., 2007;
Barclay, 2001; Molnár et al., 2014; Cook et al., 1995; Baumgaertner et al.,
2016; Majmudar et al., 2020). This process becomes more complicated when there
are many apparently equivalent mathematical functions that may be used
(Fussmann and Blasius, 2005). In this paper, we are interested in systems
where the data is typically sparse and/or noisy, making it impossible to
select one best-fitting functional curve (O’Donoghue et al., 2001). We focus
here on the case of predator-prey systems and, in particular, on the predation
interaction. The mathematical expression of this interaction is often referred
to as the functional response.

A common assumption in modeling is that, if mathematical functions possess
similar features and are appropriately parametrized to resemble one another,
then the models utilizing these functional responses will produce essentially
the same behaviour. Hence, many mathematicians typically choose a convenient
functional response that fits the predation data. It was found, however, that
similar looking functional responses may produce very different results
(Fussmann and Blasius, 2005). Indeed, Fussmann and Blasius (2005) found that,
for a given prey density, their model can produce high amplitude oscillations,
bistability, or stable coexistence at low densities, depending on the form of
the functional response, even if these responses are visually
indistinguishable. This difference in model behaviour means that the
assumption that similar curves can be used interchangeably does not always
hold. This behaviour is called structural sensitivity. A structurally
sensitive model has the property that a small perturbation in model functions,
here the functional response, can lead to a different bifurcation structure.
Ecologically, predictions made using models with different bifurcation
structures could range from stable coexistence to extinction of both species.

The predator-prey system in Fussmann and Blasius (2005) is a Rosenzweig-
MacArthur model that makes use of three similar functional response curves
that are visually indistinguishable with respect to data fitting. The three
functional responses are: Holling Type II, which is derived by considering how
the predator’s activity is split between searching for and handling prey;
Ivlev, which is based on the maximum digestion rate of the predator; and a
hyperbolic tangent (trigonometric), which is a phenomenological curve that has
no biological basis but is used in some population models (Aldebert et al.,
2018). Although the theory behind each of these functional responses is
different, they are considered interchangeable because it is difficult to
determine if a predator is limited by its handling time or its digestion time.
Jeschke et al. (2002) have developed a model that includes both of these
limitations on predation together but we do not consider it in this
manuscript. These functional responses have five specific characteristics in
common: zero predation at zero prey density, monotonically increasing
predation with increasing prey density, a negative second derivative,
saturating behaviour at high prey density, and smoothness (Fussmann and
Blasius, 2005; Seo and Wolkowicz, 2018).

Fussmann and Blasius (2005) found that the Rosenzweig-MacArthur model with
these three functional response curves is structurally sensitive. That is,
using the carrying capacity as the bifurcation parameter, this model produces
different bifurcation diagrams for all three very similar looking functional
responses. The model has a Hopf bifurcation of the coexistence steady state
for all three functional response curves, but the critical value of the
carrying capacity at which the Hopf bifurcation occurs varies by as much as
2000% from the smallest to largest critical transition value. For large values
of the carrying capacity, the model that uses the Holling functional response
produces large oscillations with limit cycles that approach very close to the
prey-only and extinction steady states. On the other hand, the model using the
trigonometric functional response initially produces smaller oscillations with
limit cycles further from the prey-only and extinction steady states. The
model using the Ivlev functional response produces oscillations of
intermediate size. Additionally, the trigonometric model is the only model
that has bistability. The difference in these results is problematic because
the model disagreements make it impossible to determine, for example, if a
real system is nearing a critical transition, such as the rapid emergence of
large oscillations or extinction (Scheffer et al., 2009). This dilemma leads
us to ask if there is a better way to mathematically describe data such that
there is better similarity of the bifurcation structures obtained from
different functional responses with respect to the bifurcation points, the
coexistence steady state, and the branches of the limit cycles.

In recent literature, Adamson and Morozov (2014, 2013) have shown the
existence of structural sensitivity in a range of mathematical models and have
developed a test to determine whether or not a model is structurally
sensitive. This test, however, does not address the question of how to alter
the functional responses in a model to reduce or remove structural
sensitivity, it only suggests that care must be taken in interpreting model
predictions.

In this paper, we investigate three different approaches for translating a
description of the data into a functional response curve. Our goal is to
determine if any of our approaches can decrease or even remove structural
sensitivity. First, we repeat the work of Fussmann and Blasius (2005) on
another predator-prey model, the Leslie-Gower-May model, and find that it is
also structurally sensitive to its functional response. We then base our study
of structural sensitivity on the behaviour of these two different predator-
prey models. Next, we study the effect of stochasticity of various types and
amplitudes on the system and observe the resulting changes in structural
sensitivity. Lastly, we investigate what part(s) of the functional response
might be related to structural sensitivity.

In Section 2, we introduce the models and functional responses. In Section 3,
we explain the methods we use to investigate structural sensitivity in the
functional response. In Section 4, we state our results. In Section 5, we
discuss and relate our findings to previous work and suggest a method to
improve predictions obtained from structurally sensitive models.

##  2 Models

In this paper, we consider two classic predator-prey models: the Rosenzweig-
MacArthur model and the Leslie-Gower-May model.

In the Rosenzweig-MacArthur model, prey density, $x$, grows via logistic
growth with a growth rate $r$ and carrying capacity $K$. Predators and prey
interact via a predation term, $f_{i}(x)$, which is discussed in more detail
below. Lastly, predator density, $y$, increases through growth related to
predation and decreases through natural death with mortality rate $m$. The
Rosenzweig-MacArthur model is written (Turchin, 2013):

$\displaystyle\frac{dx}{dt}$
$\displaystyle=rx\left(1-\frac{x}{K}\right)-f_{i}(x)y,$ (1a)
$\displaystyle\frac{dy}{dt}$ $\displaystyle=f_{i}(x)y-my.$ (1b)

The parameter values given in Fussmann and Blasius (2005) are used throughout
this paper (for the Rosenzweig-MacArthur model) unless specified otherwise.

The prey density equation in the Leslie-Gower-May model is the same as that in
the Rosenzweig-MacArthur model. The predator density equation, on the other
hand, has logistic growth of the predator with a growth rate $s$ and carrying
capacity $x/q$ that decreases as prey density decreases. Here, the parameter
$q$ is the equilibrium density ratio (Turhcin and Hanski, 1997). The Leslie-
Gower-May model is written (Turchin, 2013):

$\displaystyle\frac{dx}{dt}$
$\displaystyle=rx\left(1-\frac{x}{K}\right)-f_{i}(x)y,$ (2a)
$\displaystyle\frac{dy}{dt}$ $\displaystyle=sy\left(1-\frac{qy}{x}\right).$
(2b)

The parameter values given in Vitense et al. (2016) are used throughout this
paper (for the Leslie-Gower-May model) unless specified otherwise.

The predation functional response, $f_{i}(x)$, is our main focus. Following
Fussmann and Blasius (2005), we consider the following three functional
responses: Holling, Ivlev, and trigonometric. These functional responses are
described, respectively, by the following three functions:

$\displaystyle f_{H}(x)$ $\displaystyle=\frac{a_{H}x}{1+b_{H}x},$ (3a)
$\displaystyle f_{I}(x)$ $\displaystyle=a_{I}(1-\text{exp}(b_{I}x)),$ (3b)
$\displaystyle f_{T}(x)$ $\displaystyle=a_{T}\text{tanh}(b_{T}x),$ (3c)

where $a_{i}$ and $b_{i}$ are the two parameters through which the curves can
be fit to data, or to each other. The Rosenzweig-MacArthur functional
responses have units of prey/predator/day and the Leslie-Gower-May functional
responses have units of prey/predator/year. Note that in application, the
three functional responses fit data equally well. Further, since each
functional response has two parameters, it is not possible to distinguish the
models using criteria such as the Akaike Information Criteria or
Schwarz/Bayesian Information Criteria as these criteria rely greatly on the
number of parameters used in each functional response (Akaike, 1974).

##  3 Methods

Structural sensitivity has not previously been established in the Leslie-
Gower-May model. In this paper, we first show that model (2) is structurally
sensitive. Then we study the structural sensitivity of the Rosenzweig-
MacArthur and Leslie-Gower-May models by altering the functional responses (3)
in the following three ways:

  1. 1.

Adding various types of stochasticity (Section 3.2).

  2. 2.

Varying the parameters $a_{i}$ and $b_{i}$ to fit the functional responses
over a specified subset of the prey density domain (Section 3.3).

  3. 3.

Using continuous curves defined piecewise from two functional responses,
switching from one to the other at a point where the functional response
curves intersect (Section 3.4).

Method 1 uses stochasticity to create more overlap in the functional responses
and hypothesizes that if the functional responses are more similar, then the
resulting model behaviours will also become more similar. Methods 2 and 3
consist of focusing on different intervals of prey densities by applying a
different weight to each part of the functional response. Our goal with these
methods is to determine if there is a portion of the functional response that
dictates the bifurcation structure.

Before investigating these methods to determine if we can decrease or remove
structural sensitivity, we first fit our functional responses. To do this
fitting, one of the functional responses is fixed and the other two are fit to
the first response via nonlinear least squares to solve for the unknown
$a_{i}$ and $b_{i}$ values. We use the Nelder-Mead algorithm (Nelder and Mead,
1965) to apply nonlinear least squares and chose the interval of prey
densities to be [0,4] to match the previous work by Fussmann and Blasius
(2005).

In the case of the Rosenzweig-MacArthur model, the fixed functional response
is the Ivlev response and we fit the Holling and trigonometric functional
responses to the fixed Ivlev functional response. As expected, we recover the
same parameter values reported in Fussmann and Blasius (2005). In the case of
the Leslie-Gower-May model, we use parameter values from Vitense et al. (2016)
to define and fix the Holling functional response. Then, we fit the Ivlev and
trigonometric functional responses to the fixed Holling functional response.
The resulting Rosenzweig-MacArthur and Leslie-Gower-May functional responses
appear in Figures S1a and 1(a).

(a) Original Functional Responses

(b) Stochastic Functional Responses

(c) Fitted Functional Responses

(d) Piecewise Identical Functional Responses

Figure 1: Functional responses used in the Leslie-Gower-May models: a) the original model; b) stochastic functional responses with $\sigma=50$ where the shaded region represents one standard deviation; c) the functional responses with nonlinear least squares applied over prey densities in [0.3,1] shown by the cyan fitting region; and d) sample piecewise identical functional responses. Parameter values are given in Table 1 and intersection values for the piecewise identical responses are given in Table 3. The Rosenzweig-MacArthur functional responses are similar with a difference in time scale and are given in the Supplementary Material (Figure 1).  | Rosenzweig-MacArthur | Leslie-Gower-May
---|---|---
| Holling | Ivlev | Trigonometric | Holling | Ivlev | Trigonometric
$a_{i}$ | 3.05 | 1 | 0.99 | 1683.333 | 451.447 | 446.182
$b_{i}$ | 2.68 | 2 | 1.48 | 3.333 | 2.313 | 1.743
Table 1: Values of the parameters $a_{i}$ and $b_{i}$ used to fit the original
functional responses. Note that the Rosenzweig-MacArthur model uses a
functional response with units prey/predator/day and the Leslie-Gower-May
model uses a functional response with units prey/predator/year.

###  3.1 Structural Sensitivity

To determine if the Leslie-Gower-May model is structurally sensitive, we plot
a bifurcation diagram (shown in Figure 3(a)) using the carrying capacity , K,
as the bifurcation parameter for each of the sub-models produced using the
three functional responses (3) in model (2), and note the differences and
similarities between them.

We find that the Leslie-Gower-May model is structurally sensitive: the Hopf
bifurcation occurs at a very different carrying capacity for each sub-model.
Figure 3(a) shows that the coexistence steady state in the model using the
Holling functional response destabilizes at $K=5.55$, in the Ivlev model at
$K=12.09$, and in the trigonometric model at $K=26.36$. Additionally, the
Ivlev and trigonometric models have saddle node bifurcations of limit cycles
at $K=10.51$ and $K=12.94$, respectively, and this bifurcation is not present
in the Holling model. Following the formation of these limit cycles,
bistability is observed in the Ivlev and trigonometric models. The Ivlev model
is bistable over prey density values $K\in[10.51,12.09]$ and the trigonometric
model is bistable over prey density values $K\in[12.09,26.36]$.

###  3.2 Stochastic Model Methods

The addition of stochasticity to the functional responses makes them more like
real data and increases their overlap. We add stochasticity to models (1) and
(2) by rewriting the model equations in the following form:

$\displaystyle d\xi$ $\displaystyle=-\xi(t)dt+\sigma dW,$ (4a)
$\displaystyle\frac{dx}{dt}$ $\displaystyle=F(x,y,\xi(t)),$ (4b)
$\displaystyle\frac{dy}{dt}$ $\displaystyle=G(x,y,\xi(t)),$ (4c)

i.e., $\xi(t)$ is an Ornstein-Uhlenbeck process and $(x,y)$ solves a system of
ODEs driven by $\xi$. We use $\sigma$ to alter the amount of noise included in
each time step. $F$ is the right-hand-side of equations 1a and 2a (depending
on the model in question) with $f_{i}(x)+\xi(t)$ in place of $f_{i}(x)$.
Similarly, $G$ is the right-hand-side of equation 1b with $f_{i}(x)+\xi(t)$ in
place of $f_{i}(x)$. Note that $f_{i}(x)$ does not appear in equation 2b, so
it is unchanged. Hence, model 4 is a system of one stochastic differential
equation (SDE) and two ordinary differential equations (ODEs). At each time
step, we sample a Gaussian random variable, $dW=\sqrt{dt}Z$ where $Z\sim
N(0,1)$, and solve the SDE (equation 4a) using the Euler-Maruyama scheme to
obtain $\xi(t)$. Then, we plug $\xi(t)$ into the system of ODEs and
numerically solve the system using Euler’s method with a step size small
enough to achieve stable solutions. Note that if the solutions result in prey
or predator densities less than zero due to the stochastic element, then the
prey or predator density is set to zero. For the Rosenzweig-MacArthur model,
we use $\sigma=0.01$ and for the Leslie-Gower-May model, we use $\sigma=50$
unless specified otherwise. Sample stochastic functional responses are given
in Figures 1(b) and S1b. The corresponding results are given in Section 4.1.

###  3.3 Fitted Functional Response Methods

The behaviour of both the Rosenzweig-MacArthur and the Leslie-Gower-May models
is stable coexistence for values of the carrying capacity, $K$, less than the
destabilizing bifurcation. We hypothesize that if the model behaviour is more
similar leading up to the destabilizing bifurcation, then perhaps the
destabilizing bifurcation will occur for more similar values of the carrying
capacity. If our hypothesis is correct, it suggests that a key portion of the
functional response is in the vicinity of the coexistence steady state.

We choose an interval containing the prey density values of the coexistence
steady state to perform the fit. We call this interval the fitting region and
choose it such that it is smaller than the original $[0,4]$ fitting region.
Note that we choose this region regardless of whether the coexistence steady
state is stable or unstable. In models that do not use data, we choose a
functional response to be fixed. We will refer to the sub-model where Holling
is chosen as the fixed functional response as the fixed Holling case and the
other sub-models follow the same naming scheme. Note that in the applied case
of this method, we do not have to choose a fixed functional response since we
can fit the functional responses to data. We then apply nonlinear least
squares to the functional responses fit to either the fixed functional
response or to data over the fitting region. The average prey density values
of the coexistence steady states in the Rosenzweig-MacArthur and Leslie-Gower-
May models are $x=0.05$ and $x=0.65$, respectively. So, the corresponding
fitting regions are taken to be $x\in[0,0.1]$ and $x\in[0.3,1]$. The parameter
values obtained from this fit to the fixed Holling functional response are
given in Table 2. The fixed Ivlev and fixed trigonometric functional responses
are given in Figures S1c and S4. The corresponding parameter values are given
in Table S1. The Holling fitted functional response results are given in
Section 4.2.

| Rosenzweig-MacArthur | Leslie-Gower-May
---|---|---
| Holling | Ivlev | Trigonometric | Holling | Ivlev | Trigonometric
$a_{i}$ | 3.05 | 0.6282 | 0.3707 | 1683.333 | 399.8009 | 384.1465
$b_{i}$ | 2.68 | 4.8204 | 7.6281 | 3.333 | 3.1656 | 2.4028
Table 2: Values of the parameters obtained from the fitted functional response
where the Holling functional response is fixed and the other functional
responses are fit to the Holling response via the parameters $a_{i}$ and
$b_{i}$.

###  3.4 Piecewise Identical Functional Response Methods

Here, rather than fitting the functional responses to each other, we modify
the functional responses so that they are exactly the same over portions of
the prey density domain $[0,4]$. Each functional response has three
intersection points with each of the other two functional responses. All three
functional responses intersect at (0,0), since there is zero predation at zero
prey density. The values for the non-trivial intersections of the Rosenzweig-
MacArthur and the Leslie-Gower-May functional response are given in Table 3.
As a result, we can consider piecewise functional responses made up of
contributions from two functional responses. These piecewise functional
responses give us another way to ask if there is a key part of the functional
response that determines the model behaviour. That is, we can examine how each
portion of the functional response, (i.e., low, medium, and high prey density)
affects the model’s bifurcation structure.

We consider three cases for each functional response for each model (18 cases
total). To form each case we choose two functional responses and interchange
the functional response we use for each of the three regions. The
intersections noted earlier become the transition points between the regions.
Sample piecewise functional responses are given in Figures 1(d) and S1d. In
these figures, we introduce notation that is used from here forward to
describe the piecewise identical sub-models. This notation is a three letter
code that corresponds to each of the three prey density regions we have chosen
to investigate, with Holling, Ivlev, and trigonometric functions each
represented by their first letter. For example, HII represents the functional
response that uses the Holling functional response for low prey densities, and
uses the Ivlev functional response for medium and high prey densities. The
corresponding results are given in Section 4.3.

Rosenzweig-MacArthur | Holling/Ivlev | Ivlev/Trig | Trig/Holling
---|---|---|---
Extinction | 0 | 0 | 0
Low to medium prey density | 0.6191 | 0.5651 | 0.5942
Medium to high prey density | 2.5799 | 2.1597 | 2.4695
Leslie-Gower-May | Holling/Ivlev | Ivlev/Trig | Trig/Holling
Extinction | 0 | 0 | 0
Low to medium prey density | 0.5726 | 0.4458 | 0.5152
Medium to high prey density | 2.4486 | 2.2611 | 1.8077
Table 3: Prey density values for the intersections in the functional responses
used to form the piecewise identical functional responses.

##  4 Results

###  4.1 Stochastic Model Results

In the case of the Rosenzweig-MacArthur model, we added Gaussian stochasticity
with $\sigma=0.01$. We choose this amplitude because stochasticity of smaller
amplitude had virtually no effect on the bifurcation structure. $\sigma=0.01$
stochasticity adds some variance to the time series but does not significantly
alter the qualitative behaviour for carrying capacities less than the Hopf
bifurcation in the deterministic model. That is, given initial conditions that
lead to the stable coexistence steady state in the deterministic case, the
stochastic model also results in trajectories that converge to the coexistence
steady state. For values of the carrying capacity greater than that of the
bifurcation in the deterministic model, the lower side of the limit cycle
reaches very low prey densities of approximately $10^{-60}$ at times
($10^{-58}$ times smaller than the stochasticity). As a result, the
stochasticity causes extinction within 100 time steps in most of the time
series of the stochastic model for carrying capacities greater than that of
the deterministic model’s destabilizing bifurcation. These results mean that
we are unable to clearly observe the stochastic behaviour and we cannot
conclude whether or not stochasticity reduces structural sensitivity in the
Rosenzweig-MacArthur model.

In the case of the Leslie-Gower-May model, we added Gaussian stochasticity
with $\sigma=50$. As in the case of the Rosenzweig-MacArthur model, we chose
this amplitude of stochasticity because it is large enough to affect the
bifurcation structure, but not so large that all of trajectories rapidly
terminate in extinction of the predator or both populations. Results obtained
from both smaller and larger amplitudes of stochasticity are provided in the
Supplementary Material (Figure S3).

The bifurcation diagrams (quasi-steady state behaviour) for the stochastic
Leslie-Gower-May model are given in Figure 3(b). This figure shows better
similarity in the bifurcation structures with the addition of stochasticity.
The bifurcation structures are obtained by plotting the maximum and minimum of
the stochastic time series (after removal of transients). Thus, these windows
show the smallest value of the carrying capacity for which the coexistence
steady states may destabilize. Limit cycles in our stochastic model appear to
have greater amplitude than those of the deterministic model. We find that the
stochastic bifurcation structures show better similarity than the
deterministic model’s bifurcation structures. This similarity suggests that
there is a decrease in structural sensitivity of the model.

Plots of time series (Figure S2) indicate that bistability remains a property
of the stochastic Ivlev and trigonometric models. In the bistable regions, the
solution trajectories of the Ivlev and trigonometric models alternate randomly
between the coexistence steady state and the stable limit cycle, and also
spend some time at the unstable limit cycle. This ability to alternate between
stable and unstable steady states is unique to the stochastic model, and has
been observed in other work (Sharma et al., 2015; Abbott and Nolting, 2017).
These observations occur for small to intermediate stochastic amplitudes,
larger amplitudes lead to rapid extinction, and more generally, to the
stochastic behaviour washing out the deterministic structure.

(a) Original Functional Response

(b) Fitted Functional Response

(c) Piecewise Identical Functional Response

Figure 2: Bifurcation diagrams for the Rosenzweig-MacArthur model (1): a) the
original model; b) sub-models using least squares applied to the functional
responses over the prey density interval $[0,0.1]$ where the Holling
functional response is fixed; and c/d) sub-models using piecewise identical
functional responses. The corresponding functional responses are given in
Figure S1 and the bifurcation values are given in Table 4. As described in
Section 4.1, the addition of stochasticity causes extinction in the
Rosenzweig-MacArthur model and thus, no bifurcation diagrams could be
produced.

(a) Original Functional Response

(b) Stochastic Functional Response

(c) Fitted Functional Response

(d) Piecewise Identical Functional Response

Figure 3: Bifurcation diagrams for the Leslie-Gower-May model (2): a) the
original model; b) our stochastic Leslie-Gower-May model with $\sigma=50$
where both the deterministic and stochastic models were obtained by running
time series and discarding the transient solutions, so these diagrams only
show the quasi-steady state behaviour; c) sub-models using least squares
applied to the functional responses over the prey density interval $[0.3,1]$
where the Holling functional response is fixed; and d) sub-models using
piecewise identical functional responses. The corresponding functional
responses are given in Figure 1 and the bifurcation values are given in Table
4.

###  4.2 Fitted Functional Response Results

We only consider the fixed Holling case here as the fixed Ivlev and fixed
trigonometric cases are similar and given in the Supplementary Material
(Figures S4 and S5).

The Rosenzweig-MacArthur fitted functional response sub-models yield the
bifurcation plots shown in Figure 2(b). The bifurcation values for these sub-
models are given in Table 4. These fits result in better similarity of the
bifurcation structures and less structural sensitivity as compared to the
original model (Figure 2(a)). For example, when the Ivlev functional response
is fit to the fixed Holling functional response, the Ivlev fit to Holling sub-
model shows behaviour similar to the Holling only sub-model as there is a
leftward shift in the destabilizing Hopf bifurcation. This result holds across
each of the fits. Hence, considering the prey densities in the chosen fitting
region appears to be a better metric to describe functional responses in the
Rosenzweig-MacArthur model; in the sense that model behaviour is less
sensitive to small changes in the functional response using this method.

The Leslie-Gower-May fitted functional response sub-models yield the
bifurcation diagrams given in Figure 3(c) and the value of the bifurcation
parameter for each bifurcation is given in Table 4. The resulting sub-models
show better similarity in bifurcation structure when compared to the fixed
functional response. For example, Ivlev fit to Holling shows a leftward shift
of the Hopf bifurcation as compared to the Ivlev only bifurcation structure
and the disappearance of the saddle bifurcation of limit cycles. So, the Ivlev
fit to Holling sub-model shows behaviour similar to the Holling only sub-
model. The same result holds for each of the other cases, and there is better
similarity in model behaviour when compared to the fixed functional response.
Therefore, as in the Rosenzweig-MacArthur model, fitting the functional
responses over the chosen fitting region appears to decrease structural
sensitivity in the Leslie-Gower-May model.

We would like to note that these theoretical results are subjective to which
functional response is chosen to be fixed. That is, the model behaviour will
always become more similar to the fixed model. In the case of the results
shown in this section, the model behaviour always becomes more similar to the
Holling sub-model behaviour. Our method when applied to a model that uses data
to fit the functional responses, however, is not subjective. We will discuss
this application of our method to data further in Section 5.2.

| Rosenzweig-MacArthur | Leslie-Gower-May
---|---|---
Original Model (2(a), 3(a)) | Holling | Ivlev | Trig | Holling | Ivlev | Trig
Hopf | 0.4452 | 1.071 | 10.12 | 5.554 | 12.09 | 26.36
Limit Cycle Saddle | - | - | 2.644 | - | 10.51 | 12.94
Holling Fixed (2(b), 3(c)) | Holling | Ivlev | Trig | Holling | Ivlev | Trig
Hopf | 0.4452 | 0.4632 | 0.7729 | 5.554 | 5.857 | 6.052
Limit Cycle Saddle | - | - | 0.5057 | - | - | 5.964
Pw-I Holling (2(c), 3(d)) | HTT | THT | HHT | HTT | THT | HHT
Hopf | 0.4451 | 0.4451 | 0.4451 | 26.36 | 26.36 | 5.554
Limit Cycle Saddle | - | - | - | 7.95 | 7.95 | -
Pw-I Trigonometric (2(c), S8) | THH | HTH | TTH | THH | HTH | TTH
Hopf | 10.12 | 10.12 | 10.12 | 5.554 | 5.554 | 26.36
Limit Cycle Saddle | 2.431 | 2.431 | 2.644 | - | - | 13.13
Table 4: Values of the bifurcation parameter, the carrying capacity, for each
bifurcation shown in the respective figures. The values in parentheses give
the figure number of the corresponding bifurcation diagrams. Pw-I stands for
piecewise identical, H stands for Holling, and T stands for trigonometric.

###  4.3 Piecewise Identical Functional Response Results

Suppose that medium prey densities determine the bifurcation structure of a
model. Then, according to our hypothesis, we expect YXY and ZXZ to show model
behaviours more similar to the X model behaviour (be that H, I, or T). Based
on the results in Section 4.2, we might actually expect low prey densities to
be the most important, in which case we expect XYY and XZZ to have better
similarity to the X model behaviour. However, our results show that none of
the three regions we consider separately determine the model behaviour. In
this section, we describe only the piecewise identical functional responses
made up of the Holling and trigonometric functional responses. The bifurcation
diagrams for these cases of the Rosenzweig-MacArthur and Leslie-Gower-May sub-
models are given in Figures 2(c) and 3(d). The rest of the cases are similar
and are given in the Supplementary Material (Figures S6-S8).

In the case of the Rosenzweig-MacArthur model, there is not better similarity
of the Holling to trigonometric sub-models (i.e., sub-models of the form HTT,
THT, or HHT) behaviour towards the trigonometric behaviour. This result
suggests that both the first and second regions determine the resulting
bifurcation behaviour of the Rosenzweig-MacArthur model and that these regions
are the most structurally sensitive.

However, the Rosenzweig-MacArthur trigonometric to Holling piecewise identical
sub-models (i.e., sub-models of the form THH, HTH, or TTH) are slightly
different. In the THH plot, there is a small leftward shift of the saddle
bifurcation of limit cycles from the original value (Table 4) at $K=2.644$ to
$K=2.431$. There is, however, no change in the location of the Hopf
bifurcation. In the HTH plot, the same leftward shift occurs. Lastly, the TTH
plot has a resulting bifurcation structure that matches that of the
trigonometric only functional response model behaviour. Overall, the
trigonometric functional response has a strong effect on the bifurcation
structure. Further, when the Holling functional response takes up two of the
three regions (in the THH and HTH cases), the Holling functional response has
an effect, though limited, on the bifurcation structure. This effect causes
the aforementioned shift in the saddle bifurcation of limit cycles. The
piecewise identical functional responses seem to affect the resulting model
behaviour without following any clear rules. Hence, we are unable to make any
claims on which prey density regions determine the Rosenzweig-MacArthur model
behaviour.

In the case of the Leslie-Gower-May model, the HTT and THT sub-models show
better similarity when compared to the trigonometric bifurcation structure
while the resulting bifurcation structure for the HHT sub-model matches that
of the Holling only functional response. We found that in the Leslie-Gower-May
model, the functional response that is used across more prey density regions
is the one that determines where the Hopf bifurcation occurs and whether there
exists a saddle bifurcation of limit cycles. However, the functional response
used across only one prey density region causes a shift in the saddle
bifurcation of limit cycles (if it exists) for every case, except surprisingly
the IIH sub-model where there is no shift from the Ivlev behaviour (Figure
S8). This result means that the sub-models with the Holling functional
response describing two prey density regions result in only Holling model
behaviour, but the sub-models with either the Ivlev or trigonometric
functional response describing two prey density regions result in mixed
dynamics.

##  5 Discussion

###  5.1 Reducing Structural Sensitivity

Structural sensitivity is an intrinsic feature of a wide range of ecological
models (Adamson and Morozov, 2013, 2014). In our work, we focus on predator-
prey models and investigate the structural sensitivity that occurs with
respect to the predation function. Our work has attempted to answer how to fit
data in such a way to gain confidence in the results, regardless of which
curve is used. We found that only the sub-models that use the functional
responses fit around the coexistence steady state consistently decrease
structural sensitivity across both the Rosenzweig-MacArthur and Leslie-Gower-
May models. Using this result, we also give suggestions as to which intervals
of prey density are most relevant when fitting a functional response to data
and which intervals of prey density are most relevant for ecologists to
collect data.

####  5.1.1 Stochasticity

One of the differences between the functional response curves in Fussmann and
Blasius (2005) and real data is that the latter include stochasticity. A
model’s steady states and deterministic structure only provide a partial
understanding of the stochastic structure and so the stochastic system may
behave in unanticipated ways (Sharma et al., 2015). Hence, it is important to
consider a stochastic model when analyzing biological systems. We examine the
effect of adding stochasticity to the functional responses of the Rosenzweig-
MacArthur and Leslie-Gower-May models and, specifically, whether or not this
addition can remove or decrease structural sensitivity.

In our work, we observe that the destabilizing of the coexistence steady state
can occur earlier in the stochastic system than the deterministic system. Most
notably, there is greater similarity in the stochastic bifurcation structures
from the three functional responses implying decreased structural sensitivity.
We find that the maximum and minimum values of the deterministic system lie
within the interval determined by the maximum and minimum values of the
stochastic system. This result supports the observation that trajectories in
the stochastic system may approach the steady states of the deterministic
system (Aguirre et al., 2013; Predescu et al., 2007). How stochasticity alters
a model’s behaviour from the deterministic system depends on the amplitude and
distribution of the stochastic random variable. For example, white noise
weakens bistable behaviour while red noise can amplify the bistability
observed in a system (Sharma et al., 2015). Using a different distribution can
change the shape and depth of potential wells, and thus produce different
stability properties (Sharma et al., 2015). We do not, however, study other
stochastic distributions here and leave this investigation for future work.

White noise is widely used in studies of stochastic models (van der Bolt et
al., 2018). It is characterized by having the same variance at all
frequencies. However, in recent years, there has been interest in varying the
colour, or auto-correlation, of applied noise. Red noise, which has more auto-
correlation and is characterized by low frequencies, has been found to better
describe the stochasticity found in ecological systems (Vasseur and Yodzis,
2004; Greenman and Benton, 2003; Grove et al., 2020; Jonsson and Wennergren,
2019). In this paper, we only consider white noise. An interesting future step
is to consider the effect of red noise on the bifurcation structure of the
stochastic system.

####  5.1.2 Fitted Functional Responses

It is sometimes possible to improve mathematical model predictions through
describing data in a particular way. In order to help ecologists make accurate
decisions regarding conservation efforts, it is important to obtain the best
model predictions possible. We discuss the use of data in these model
predictions further in Section 5.2 below.

In our models, the slopes of the functional response curves are greatest near
the coexistence steady state, meaning that this region is also where the
predator response to a change in prey density is strongest (Aldebert et al.,
2016). This property suggests that having a higher density of data points in
the region of the functional response domain near the coexistence steady state
prey density value, or giving these data points more weight, is important.
Further, since the initial behaviour of all three models is similar and all
are stable at the coexistence steady state before any destabilizing occurs, we
hypothesized that if the corresponding part of the functional response (region
around the stable or unstable coexistence steady state prey density value)
held more weight in the nonlinear least squares fitting, then the bifurcation
structures would be more similar. That is, if the fitted functional responses
method (Section 3.3) is a better metric to decrease structural sensitivity, we
expect to observe better similarity of the resulting model’s bifurcation
structure. We found that this hypothesis is supported in both the Rosenzweig-
MacArthur and the Leslie-Gower-May models.

####  5.1.3 Piecewise Identical Functional Responses

We found that more similar bifurcation structures are obtained for functional
responses that fit very closely at low and intermediate prey density values,
but also disagree wildly at large prey density values near saturation. So, the
piecewise identical functional response results suggest that low and medium
prey densities are the most important regions to fit to decrease structural
sensitivity.

Our results, however, appear to be in contrast from previous findings,
indicating that it is most important to have data at saturation or high prey
densities (Aldebert and Stouffer, 2018). However, this contradiction is only
apparent. In our work, we attempt to decrease structural sensitivity and
describe the functional responses such that any of the three can be used to
provide equally accurate predictions. On the other hand, Aldebert and Stouffer
(2018) aim to choose one functional response that fits the model best. To do
so, they suggest collecting data around saturation such that the functional
responses result in models with vastly different bifurcation structures. Since
structural sensitivity in the resulting models is increased, it becomes easier
to determine which functional response and model fit the data best. While the
authors are able to choose a best functional response to fit data, the model
remains structurally sensitive.

###  5.2 Using Data in Models

The main application of our work is in models that use data to fit the
functional response. As we noted in Section 4.2, the choice of the fitted
functional response is subjective if our method is used merely theoretically.
However, in the case of applied models, our method is no longer subjective
since we fit the functional responses to data instead of a subjectively chosen
fixed functional response. Once all three functional responses are fit using
our method to make their bifurcations structures more similar, the resulting
bifurcation structure should show the behaviour of the data rather than a
single chosen functional response.

Ecological data, however, is typically sparse as it is often difficult to
collect (Ovaskainen and Soininen, 2011). Furthermore, the data is also
typically very noisy. As a result, it is important that ecologists know what
data is most important to collect to build accurate models that can be used
with confidence to predict future population dynamics. In Section 4.2, we
found that model predictions are most consistent if the functional responses
match over the coexistence steady state prey density region. The implication
is that better predictions can be obtained if there is additional data near
the coexistence steady state. We believe that this region of prey densities is
most important in fitting the functional responses because it corresponds to
the steady state behaviour for values of the carrying capacity less than
destabilizing bifurcation. The idea being that if the model behaviour is more
similar before the destabilizing bifurcation, then the destabilizing
bifurcation may occur for more similar values of the carrying capacity. This
result opposes previous work that states saturation is the most important part
to match of the functional responses (Aldebert and Stouffer, 2018).

To apply our result to ecological data and to make predictions in predator-
prey models, we propose the following procedure:

  1. 1.

Collect data as per normal methods,

  2. 2.

Fit the three functional responses to this data and determine the parameter
values for each $a_{i}$ and $b_{i}$,

  3. 3.

Use the functional responses to choose a fitting region around the stable or
unstable coexistence steady state prey density values across all three models,

  4. 4.

Do one or both of the following:

    1. (a)

Collect more data for prey densities in the coexistence steady state prey
density region,

    2. (b)

Give more weight to the data in the coexistence steady state prey density
region,

  5. 5.

Fit the three functional responses again and use these functional responses to
make predictions regarding the system.

This algorithm takes into account our result that collecting data for prey
densities around the coexistence steady state appears to be most important
while also considering the difficulty of collecting data (step 4b). Note that
we have only considered structural sensitivity with respect to the functional
response in this work. Models that are structurally sensitive with respect to
other terms, such as prey growth, may require a different algorithm.

##  6 Conclusion

Structural sensitivity in the functional response of a model can cause similar
functional responses to produce very different bifurcation structures. In this
paper, we investigated three potential methods to achieve better similarity of
resulting bifurcation structures. We found that 1) stochasticity results in
better similarity of the bifurcation structures, 2) functional responses fit
over the coexistence steady state prey density region seem to decrease
structural sensitivity, and 3) piecewise identical functional responses do not
consistently decrease structural sensitivity. Overall, adding stochasticity
and fitting functional responses over the coexistence steady state prey
density region appear to reduce structural sensitivity. Further investigation
into the fitted functional response method is necessary to determine if there
is an optimal interval size that decreases structural sensitivity the most.

Investigating other types of stochasticity, specifically red noise or noise
with stronger auto-correlation, remains an interesting question for further
investigation. Working with ecological data to test our proposed procedure and
make predictions regarding the population dynamics of a predator-prey system
is a very useful application of this work. It is also of interest to
investigate more predator-prey models and more functional response curves to
determine if our results generalize. It is also possible to have structural
sensitivity occur with respect to the prey growth function (Adamson and
Morozov, 2013). This type of structural sensitivity remains an area for future
investigation.

As part of this work, we also investigated the derivatives of the functional
responses. That is, we fit the parameters $a_{i}$ and $b_{i}$ such that the
derivatives of the functional responses were as similar as possible. This
method of fitting did not show promise and requires further investigation to
determine if it can decrease structural sensitivity.

## References

  * Abbott and Nolting (2017) Abbott, K.C., Nolting, B.C.: Alternative (un)stable states in a stochastic predator–prey model.  Ecol. Complex. 32, 181–195 (2017).  DOI 10.1016/j.ecocom.2016.11.004
  * Adamson and Morozov (2013) Adamson, M.W., Morozov, A.Y.: When can we trust our model predictions? unearthing structural sensitivity in biological systems.  Proc. Math. Phys. Eng. Sci. 469(2149) (2013).  DOI 10.1098/rspa.2012.0500
  * Adamson and Morozov (2014) Adamson, M.W., Morozov, A.Y.: Defining and detecting structural sensitivity in biological models: developing a new framework.  J. Math. Biol. 69, 1815–1848 (2014).  DOI 10.1007/s00285-014-0753-3
  * Aguirre et al. (2013) Aguirre, P., González-Olivares, E., Torres, S.: Stochastic predator–prey model with allee effect on prey.  Nonlinear Anal. Real World Appl. 14(1), 768–779 (2013).  DOI 10.1016/j.nonrwa.2012.07.032
  * Akaike (1974) Akaike, H.: A new look at the statistical model identification.  IEEE Trans. Automat. Contr. 19(6), 716–723 (1974).  DOI 10.1109/tac.1974.1100705
  * Aldebert et al. (2018) Aldebert, C., Kooi, B.W., Nerini, D., Poggiale, J.C.: Is structural sensitivity a problem of oversimplified biological models? insights from nested dynamic energy budget models.  J. Theor. Biol. 448, 1–8 (2018).  DOI 10.1016/j.jtbi.2018.03.019
  * Aldebert et al. (2016) Aldebert, C., Nerini, D., Guaduchon, M., Poggiale, J.C.: Structural sensitivity and resilience in a predator-prey model with density-dependent mortality.  Ecol. Complex. 28, 163–173 (2016).  DOI 10.1016/j.ecocom.2016.05.004
  * Aldebert and Stouffer (2018) Aldebert, C., Stouffer, D.B.: Community dynamics and sensitivity to model structure: towards a probabilistic view of process-based model predictions.  J. R. Soc. Interface 15(149), 20180,741 (2018).  DOI 10.1098/rsif.2018.0741
  * Barclay (2001) Barclay, H.J.: Modeling incomplete sterility in a sterile release program: interactions with other factors.  Popul. Ecol. 43, 197–206 (2001).  DOI 10.1007/s10144-001-8183-7
  * Baumgaertner et al. (2016) Baumgaertner, B.O., Tyson, R.C., Krone, S.M.: Opinion strength influences the spatial dynamics of opinion formation.  J. Math. Sociol. 40(4), 207–218 (2016).  DOI 10.1080/0022250X.2016.1205049
  * Cook et al. (1995) Cook, J., Tyson, R., White, J., Rushe, R., Gottman, J., Murray, J.: Mathematics of marital conflict: Qualitative dynamic mathematical modeling of marital interaction.  J. Fam. Psychol. 9(2), 110–130 (1995).  DOI 10.1037/0893-3200.9.2.110
  * Fussmann and Blasius (2005) Fussmann, G., Blasius, B.: Community response to enrichment is highly sensitive to model structure.  Biol. Lett. 1(1), 9–12 (2005).  DOI 10.1098/rsbl.2004.0246
  * Greenman and Benton (2003) Greenman, J.V., Benton, T.G.: The amplification of environmental noise in population models: Causes and consequences.  Am. Nat. 161(2), 225–239 (2003).  DOI 10.1086/345784
  * Grove et al. (2020) Grove, M., Borg, J.M., Polack, F.: Coloured noise time series as appropriate models for environmental variation in artificial evolutionary systems (2020).  DOI 2006.16204
  * Jeschke et al. (2002) Jeschke, J.M., Kopp, M., Tollrian, R.: Predator functional responses: Discriminating between handling and digesting prey.  Ecol. Monogr. 72(1), 95–112 (2002).  DOI 10.2307/3100087
  * Jonsson and Wennergren (2019) Jonsson, A., Wennergren, U.: Approximations of population growth in a noisy environment: on the dichotomy of non-age and age structure.  Theor. Ecol. 12, 99–110 (2019).  DOI 10.1007/s12080-018-0391-2
  * Krkošek et al. (2007) Krkošek, M., Ford, J.S., Morton, A., Lele, S., Myers, R.A., Lewis, M.A.: Declining wild salmon populations in relation to parasites from farm salmon.  Science 318(5857), 1772–1775 (2007).  DOI 10.1126/science.1148744
  * Majmudar et al. (2020) Majmudar, J.R., Krone, S.M., Baumgaertner, B.O., Tyson, R.C.: Voter models and external influence.  J. Math. Sociol. 44(1), 1–11 (2020).  DOI 10.1080/0022250X.2019.1625349
  * Molnár et al. (2014) Molnár, P.K., Lewis, M.A., Derocher, A.E.: Estimating allee dynamics before they can be observed: Polar bears as a case study.  PLOS ONE 9(1), 1–12 (2014).  DOI 10.1371/journal.pone.0085410
  * Nelder and Mead (1965) Nelder, J.A., Mead, R.: A simplex method for function minimization.  Comput. J. 7(4), 308–313 (1965).  DOI 10.1093/comjnl/7.4.308
  * O’Donoghue et al. (2001) O’Donoghue, M., Boutin, S., Murray, D.L., Krebs, C.J., Hofer, E.J., Breitenmoser, U., Christine Breitenmoser-Würsten, G.Z., Doyle, C., Nams, V.O.: Mammalian predators: Coyotes and lynx.  In: Krebs, C.J., Boutin, S., Boonstra, R. (eds.) Ecosystem Dynamics of the Boreal Forest: The Kluane Project, 1 edn., chap. 13, pp. 275–323. OUP (2001)
  * Ovaskainen and Soininen (2011) Ovaskainen, O., Soininen, J.: Making more out of sparse data: Hierarchical modeling of species communities.  Ecology 92(2), 289–95 (2011).  DOI 10.1890/10-1251.1
  * Predescu et al. (2007) Predescu, M., Sirbu, G., Levins, R., Awerbuch-Friedlander, T.: On the dynamics of a deterministic and stochastic model for mosquito control.  Appl. Math. Lett. 20(8), 919–925 (2007).  DOI 10.1016/j.aml.2006.12.001
  * Scheffer et al. (2009) Scheffer, M., Bascompte, J., Brock, W.A., Brovkin, V., Carpenter, S.R., Dakos, V., Held, H., van Nes Egbert, H., Rietkerk, M., Sugihara, G.: Early-warning signals for critical transitions.  Nature 461(7260), 53–59 (2009).  DOI 10.1038/nature08227
  * Seo and Wolkowicz (2018) Seo, G., Wolkowicz, G.S.K.: Sensitivity of the dynamics of the general rosenzweig–macarthur model to the mathematical form of the functional response: a bifurcation theory approach.  J. Math. Biol. 76, 1873–1906 (2018).  DOI 10.1007/s00285-017-1201-y
  * Sharma et al. (2015) Sharma, Y., Abbott, K.C., Dutta, P.S., Gupta, A.K.: Stochasticity and bistability in insect outbreak dynamics.  Theor. Ecol. 8, 163–174 (2015).  DOI 10.1007/s12080-014-0241-9
  * Turchin (2013) Turchin, P.: Trophic interaction.  In: Complex Population Dynamics, Monogr. Popul. Biol., chap. 4, pp. 95–100. PUP (2013)
  * Turhcin and Hanski (1997) Turhcin, P., Hanski, I.: An empirically based model for latitudinal gradient in vole population dynamics.  Am. Nat. 149, 842–874 (1997).  DOI 10.1086/286027
  * van der Bolt et al. (2018) van der Bolt, B., van Nes, E.H., Bathiany, S., Vollebregt, M.E., Scheffer, M.: Climate reddening increases the chance of critical transitions.  Nat. Clim. Change. 8, 478–484 (2018).  DOI 10.1038/s41558-018-0160-7
  * Vasseur and Yodzis (2004) Vasseur, D., Yodzis, P.: The color of environmental noise.  Ecology 85(4), 1146–1152 (2004).  DOI 10.1890/02-3122
  * Vitense et al. (2016) Vitense, K., Wirsing, A.J., Tyson, R.C., Anderson, J.J.: Theoretical impacts of habitat loss and generalist predation on predator–prey cycles.  Ecol. Modell. 327, 85–94 (2016).  DOI 10.1016/j.ecolmodel.2016.02.002