Analysis of the Incidence of Poxvirus on the Dynamics between Red and Grey Squirrels

Abstract

A model for the interactions of the invasive grey squirrel species as asymptomatic carriers of the poxvirus with the native red squirrel is presented and analyzed. Equilibria of the dynamical system are assessed, and their sensitivity in terms of the ecosystem parameters is investigated through numerical simulations. The findings are in line with both field and theoretical research. The results indicate that mainly the reproduction rate of the alien population should be drastically reduced to repel the invasion, and to achieve disease eradication, actions must be performed to reduce the intraspecific transmission rate; also, the native species mortality plays a role: if grey squirrels are controlled, increasing it may help in the red squirrel preservation, while the invaders vanish; on the contrary, decreasing it in favorable situations, the coexistence of the two species may occur. Preservation or restoration of the native red squirrel requires removal of the grey squirrels or keeping them at low values. Wildlife managers should exert a constant effort to achieve a harsh reduction of the grey squirrel growth rate and to protect the remnant red squirrel population.

1. Introduction

The Eurasian red squirrel (Sciurus vulgaris Linnaeus, 1758) is an arboreal mammal with a large Palearctic range, and in Europe, it is a native and widespread species in most areas, with a few exceptions. It is absent from the southwest of Spain and from the majority of Mediterranean islands, and it occurs only sporadically in the Balkans [1]. The species is known to locally disappear because of the competition with the Eastern grey squirrel (Sciurus carolinensis Gmelin, 1788), which has been introduced from North America into the British Isles, mainly from 1890 onwards [2,3], and in Italy [4], where the species was first reported in 1948. Because of its severe impacts on the red squirrel and on forest biodiversity [5], the grey squirrel was previously reported by the International Union for Conservation of Nature (IUCN) in the list of 100 worst invasive species [6]. More recently, it has been included in the list of invasive alien species of European concern (EU Regulation No. 1143/2014; EU Implementing Regulation No. 1141/2016).

The grey squirrel out-competes the smaller red squirrel in the use of trophic resources. Both species feed on tree seeds, and they largely overlap in their use of habitat [7,8,9]. The work in [10] suggested that interspecific competition for scatterhoarded seeds, with grey squirrels pilfering red squirrels’ food caches, can play a role in the replacement of red by grey squirrels. In sympatric conditions, when the core area of red squirrel home ranges largely overlaps with grey squirrel ranges, these authors observed a reduction in the daily energy intake by red squirrels. The spring body mass of red squirrels was also negatively correlated with the percentage of interspecific core area overlap. In a high-quality mixed deciduous woodland in Italy, [9] estimated a large tree-species niche overlap (about 70%), and taking into account both activity and foraging behavior in time and space, there was no evidence of interspecific niche partitioning. Thus, red squirrels seemed unable to avoid competition with grey squirrels. Furthermore, the invasive species often occurs at higher densities than the native one, and it has higher breeding rates. When red squirrels tend to avoid good quality habitats with high grey densities, their demography may be negatively affected, with a decrease in juvenile recruitment [11]. Available data thus suggest that competition between red and grey squirrels for habitat and food is very likely and that it may affect red squirrel population dynamics, finally playing a role in the observed replacement of the red species by the grey squirrel.

In the United Kingdom, the interspecific competition between squirrels is also mediated by a squirrel poxvirus (SQPV), the Chordopoxviridae genus [12,13]. The virus causes exudative erythematous dermatitis [14], characterized by ulcerated and hemorrhagic scabs affecting the skin around the eyes, nose and lips first, then spreading to the ventral thorax, inguinal area and feet [15,16]. The disease was unrecorded before the grey squirrel introduction [16]; it was then reported in England since 1930 [2] and later confirmed by electron microscopy in several regions [14,16]. Recently, the virus was also detected in Scotland [17]. Serological evidence supported the theory that the virus was introduced to the U.K. with the grey squirrel [13] and [18] proved that the parapoxvirus can cause a debilitating disease in red squirrels, whilst having no apparent impact on greys. According to the experimental tests of these authors, in red squirrels, the disease effects would include ulceration and infection of the skin lesions, along with lethargy, almost certainly causing mortality in the wild within 10–20 days. Only a few red squirrels could survive the infection [19] or at least they could experience a prolonged disease course, so that ultimately the disease may also be transmitted between red squirrels in the wild [18]. Transmission of the poxviruses differs between genera, being spread either by aerosol, contaminated fomites, direct contact or arthropods [13]. In particular, for the SQPV, current data support a direct transmission mechanism. According to [20], in England and Wales, cases of SQPV disease occurred in red squirrels only in geographic areas with seropositive grey squirrels, and a critical community size of squirrels may be required for disease epidemics in red squirrels to occur. This is also suggested by the delay observed between the establishment of invading grey squirrels in new areas and the cases of the disease in the red squirrels. Grey squirrel must probably reach a threshold density or number before the virus is transmitted to reds. Anyway, no evidence was found that the density or number of both squirrel species was controlling the onset of disease outbreaks, as expected when micro-parasites have a reservoir host, and the even spread of cases of disease across months also suggested that a direct rather than vector-borne transmission route was more likely.

Consistent with the hypothesis of an endemic infection of low pathogenicity in most grey squirrel populations, in several United Kingdom populations [19], it was found that 61% of grey squirrels appeared to have been exposed to SQPV, whereas only four (about 3%) of the red squirrels tested had antibodies to parapoxvirus. Out of these four animals, three also had parapoxvirus-associated disease. Overall, these data clearly supported the hypothesis that parapoxvirus spill-over from grey squirrel reservoirs of infection may significantly affect red squirrel populations. As a consequence, in England, the decline of the latter has been attributed not only to direct competition, but it has also been associated with the epidemic outbreaks of the disease [15,21]. Indeed, clinical signs observed in red squirrels by [18] were the same observed on squirrels recovered dead in the wild [16,20,21,22,23], and detailed studies demonstrated that the competitive interactions alone could not account for the observed rate and pattern of replacement observed in some areas of the United Kingdom [24,25].

Since this viral disease had such a significant impact on the decline of the red squirrel in the U.K., it is particularly interesting to focus on the epidemiological dynamics of the red-grey squirrel system. In this respect, deterministic approaches such as the one developed by [26] are particularly useful to provide a clear understanding of the interaction mechanisms. Here, in particular, we develop a mathematical model for a three-population system, including grey squirrels and both healthy and infected red squirrels, in order to highlight the potential effects of the SQPV transmission on the long-term population dynamics of the native species. Even though, so far, the virus has been detected only in the United Kingdom, the development of general models is fundamental, and because grey squirrels have been mainly introduced outside North America as pets, multiple releases in the wild have occurred and are still expected. As a consequence, it is not possible to exclude that the virus will appear in new areas in the future. In this respect, our model provides insights into the dynamics of the two-species system, identifying the parameter space allowing for species coexistence and for disease eradication. Management alternatives to achieve red squirrel conservation are also discussed.

In particular, this paper attempts a different approach than other current research in the field [25,27]. While these models include space in the model formulation and tend to be more specific (e.g., [27] for the Anglesey peninsula in U.K.), here we try to concentrate on a more ample view. The main aim is indeed to obtain more general results and to possibly provide an indication of how to approach the invasion problem in different scenarios. Furthermore, in setting up our model, we consider a frequency-dependent term for disease transmission, in epidemiology also known as standard incidence. This assumption is justified in the discussion on the model formulation, but here, we stress that this approach is rather different from previous model setups that have appeared in the grey squirrel invasion literature. In this way, red squirrels can be wiped out also by the epidemics. This remark prompts the field ecologists to perform new investigations, to ascertain actual disease-related extinction, thus absolving one of the tasks of theoretical models, i.e., suggesting the possible appearance of ecological phenomena. Finally, although recent investigations suggest that squirrelpox prevalence may fluctuate in grey populations, thereby implying also a change in the disease transmission, our basic assumption removes the possibility that some of the grey squirrels are virus-free. This may very well not be the case, but apart from the mathematical simplification that ensues from having a three-dimensional system rather than dealing with four differential equations, we are actually looking at the worst possible invasive scenario, in order to preserve the utmost conservative viewpoint. The analysis of a worst scenario is also justified if we consider that, where SQPV has been present, the decline in red squirrel distribution can be 17–25 times faster than in areas where SQPV is known not to be present [28].

2. Materials and Methods

2.1. The Model

We introduce at first the model using the grey squirrel G, healthy red squirrel R and infected red squirrel I populations. All the populations are counted by numbers. As a consequence, essentially all the parameters denoting rates are frequencies, measured in $t^{-1}$. The main point is that the red squirrels are partitioned into two subpopulations, to investigate possibly the epidemic effects on them. The grey squirrels are instead assumed to be immune from the disease, although asymptomatic carriers of the virus. Therefore, interactions between healthy red squirrels and grey squirrels may lead to new cases of the disease and move the red individual to the class of infected.

Note that all squirrels are local creatures. Indeed, they usually establish home ranges whose extension varies as a function of habitat type, (e.g., from 2 to about 30 hectares for red squirrels; [29,30,31]), and they hardly travel long distances. Dispersal occurs usually once or twice a year, and movements are confined within the range of a few kilometers (e.g., for the grey squirrel, see [32,33]). Basically then, for the interaction terms between populations due to possible competition for resources, we assumed a functional response similar to the classical Holling Type II, thereby modeling the fact that if a population is large, because these are essentially sedentary animals, only the influence of the closest neighbors will most directly determine the dynamics of an individual of the other population. The disease transmission is assumed to be modeled by a standard incidence function, for the same reasons. Both of these assumptions, as in standard epidemic models, lead to a singularity in the cases respectively that both populations or the disease-affected one vanish. Mathematically, this means that the origin and the case $G=N=0$ must be excluded from the domain of the dynamical system (1). Assuming further that all the parameters are nonnegative, the model reads as follows:

$$\begin{array}{ccc}\hfill \frac{dG}{dt}& =& rG\left(1-\frac{G}{E}\right)-\frac{kRG}{G+c(R+I)},\hfill \\ \hfill \frac{dR}{dt}& =& sR\left(1-\frac{R}{K}\right)-\frac{aRG}{G+c(R+I)}-R\frac{\lambda I+\beta G}{R+I+G},\hfill \\ \hfill \frac{dI}{dt}& =& R\frac{\lambda I+\beta G}{R+I+G}-\mu I.\hfill \end{array}$$

(1)

Note that we can take the initial condition $I\left(0\right)=0$ to indicate that initially the disease may very well be absent.

The first equation describes the grey squirrel dynamics, while the remaining two illustrate respectively the healthy and infected red squirrel populations’ evolution. For the grey squirrel, logistic growth is assumed, and then, competition only with the healthy red squirrels is taken into account. Indeed, the point is that the disease is highly virulent for the latter, weakens and kills them in just a short time, so that any of their possible interference for sharing resources with the grey squirrels can be easily ruled out. Note that we take the competition term to be a saturating function both in terms of G and R because the squirrels are mainly residing in a fixed place, for which their interactions with possible neighbors are limited. Thus, even if either one of the populations grows at high values, there is a saturation in the contact rate. The direct competition with infected squirrels instead is disregarded because the latter are too weak, but we take them into account in the denominator, because if there is a large number of them, the interactions between the susceptible red squirrels and the grey ones will diminish, as more infected red squirrels will be the neighbors of the grey, and therefore, the competition damage to the grey population will be reduced.

The healthy red squirrels also follow a logistic growth in which, for the very same reasons just outlined, the intraspecific competition is confined only to their healthy similars. Note that the two carrying capacities of the grey and red squirrels are assumed to be different, respectively E and K, because in suitable habitats, the grey squirrel can live at higher densities (e.g., 1.45–2.99 individuals /ha, up to 7–10 or even more than 10 ind./ha [34,35,36,37]) than the red squirrel, whose densities are usually lower than 1–2 ind./ha [7,27,30,31,38,39,40,41,42]. Furthermore, their net reproductive rates differ, r and s, respectively [43]. The competition has the same form as in the first equation because similar considerations as for the grey squirrel competition term hold for this corresponding red squirrel term in this equation. In addition, the red squirrels suffer from the competition of the grey ones, third term of the equation, at a rate a that once more differs from the corresponding one holding for the grey squirrels, k, because of the latter having a larger body size and being more efficient in the use of trophic resources [10,44]. Further, we need to consider the disease effects. Note that there are two ways that a susceptible red squirrel can contract the disease. This can occur if it comes in contact with an infected consimilar at rate $\lambda$, which is modeled by the fourth term, or also at rate $\beta$ by the interaction with a grey individual, which is always an asymptomatic carrier, as stated earlier, the last term of the second equation. Rates of viral shedding from infected grey and red squirrels may be different [45], but a transmission rate of 3.27 between and within each species was obtained by matching model and field data on the seroprevalence in the U.K. [27,46]. Note that we take the interspecific transmission term to be in a rational form because the squirrels are mainly residing in a fixed place, for which their interactions with possible neighbors are limited. Thus, even if either one of the populations grows at high values, there is a saturation in the contact rate.

Infected red squirrels are recruited via the two infection mechanisms just described for the healthy red ones and appear as the first two terms of the last equation. In addition, they experience natural plus disease-related mortality at rate $\mu$. In view of their weakness, the fast replication process of the virus in their bodies and high lethality of the disease, no other vital dynamics is possible for this class of individuals. Indeed, pregnancy lasts longer than the disease itself: on average, they indeed respectively last 38–40 [47] and 10–20 days [18].

Table 1. A summary of the model parameters and their interpretation.

Name	Interpretation	Unit
r	grey squirrels’ net growth rate	$t^{-1}$
E	grey squirrels’ carrying capacity	pure number
s	red squirrels’ net growth rate	$t^{-1}$
k	grey squirrels’ damage due to interspecific competition	$t^{-1}$
a	red squirrels’ damage due to interspecific competition	$t^{-1}$
K	red squirrels’ carrying capacity	pure number
$\lambda$	intraspecific disease transmission rate for red squirrels	$t^{-1}$
$\beta$	interspecific disease transmission rate for red squirrels	$t^{-1}$
$\mu$	red squirrels’ natural plus disease-related mortality	$t^{-1}$
c	competitivity weight among different squirrel species	$t^{-1}$

summarizes the model parameters.

2.2. Model Reparametrization

We now reformulate Model (1) in terms of the total red squirrel population $N=R+I$, the disease prevalence $i=I{N}^{-1}$ among the latter, while keeping the grey squirrels’ population unchanged. This transformation is rather common in epidemic [48] and ecoepidemic models [49,50], in spite of the fact that it breaks down when $N=0$, because it allows distinguishing the two different vanishing modes of this population, namely when the disease remains endemic in it, affecting the whole vanishing population, or when some of the red squirrels remain still healthy in the process. Two equilibria arise for the situation in which the local population is eradicated by the successful replacement of the invaders, both corresponding to red squirrels’ disappearance in the original Model (1). They are not artificial because they allows us to understand better the role played by the epidemic in the native species disappearance. Thus, they represent a valuable piece of information that, without this transformation, we would not be able to gather. The reformulation leads to the following system:

$$\begin{array}{ccc}\hfill \frac{dG}{dt}& =& G\left[r\left(1-\frac{G}{E}\right)-\frac{kN(1-i)}{G+cN}\right],\hfill \\ \hfill \frac{dN}{dt}& =& N\left[s(1-i)\left(1-\frac{N(1-i)}{K}\right)-\frac{a(1-i)G}{G+cN}-\mu i\right],\hfill \\ \hfill \frac{di}{dt}& =& (1-i)\left[\frac{\lambda iN+\beta G}{G+N}-\mu i-si\left(1-\frac{N(1-i)}{K}\right)+\frac{aiG}{G+cN}\right],\hfill \end{array}$$

(2)

once again initially assuming no disease in the red squirrel population, $i\left(0\right)=0$.

2.3. Study of the Ecosystem Behavior

Taking into account the described model, we analyze the existence and stability of equilibrium solutions, and we also perform numerical simulations to better investigate the current situation. For the latter, we use a set of reference parameter values:

$$\begin{array}{c}\hfill r=1.2,E=1.25,k=0.61,s=1,K=0.5,a=1.65,\lambda =3.27,\beta =3.27.\end{array}$$

(3)

These values are obtained from the literature, mainly from [27,51]. Instead, the natural plus disease-related mortality for the red squirrels and the weight of the two populations are arbitrarily assumed:

$$\begin{array}{c}\hfill \mu =2,c=1.\end{array}$$

(4)

Note that the natural plus disease-related mortality $\mu$ is chosen at a high value in view of the high virulence and lethality of the epidemics, for the red squirrels. In addition, note that the literature [27,43,51] reports that some of them vary in the following suitable ranges:

$$\begin{array}{c}\hfill r\in [1.2,1.3],E\in [0.6,2.5],K\in [0.35,0.65].\end{array}$$

(5)

Thus, we also investigate the sensitivity of the system with respect to some of the ecosystem parameters. In all simulations the initial conditions are chosen as follows, with no initial infectives, as this was the pristine situation before the introduction of the invasive grey population:

$$G\left(0\right)=0.07,G\left(0\right)=2.80,G\left(0\right)=0.$$

3. Results

3.1. Equilibria

As noted earlier, the origin $P_0$ and the point at which the ecosystem collapses, but with the disease remaining endemic among the red squirrels, $P_2=(0,0,1)$, are outside the dynamical system domain. There are only four possible remaining equilibria: the healthy-red-squirrels-only point $P_1=(0,K,0)$, the best possible ecosystem outcome, then an equilibrium in which the red population collapses, but the disease remains endemic, while grey squirrels thrive, for which, however, a pair of points is detected, $P_3=(E,0,1)$ and $P_4=(E,0,i_4)$; the grey-squirrels-free point $P_5=(0,N_5,i_5)$, which is another desirable ecosystem outcome, implying the invasive grey squirrel extinction, although the disease remains endemic among the red ones, and the coexistence of the three populations $P_{*}=(G_{*},N_{*},i_{*})$, which would also be less welcome, but still acceptable, as the local individuals would thrive anyway. Note that the pair of points $P_3$ and $P_4$ denotes two really different situations: although the native squirrel population disappears, in the former case, the disease affects the whole vanishing population, while in the latter, it is still endemic, but part of the red squirrels are still healthy.

3.1.1. Feasibility

Specifically, for the red-squirrels-free equilibrium, observe that there are two cases depending on the status of the disease in the population, while the species itself is disappearing. In both, the grey squirrels attain carrying capacity. In the former, $P_3$, the totality of red squirrels is infected, while the population declines to zero, while for the second, equilibrium $P_4$, for the disease prevalence among the vanishing red squirrels, we find the value:

$$i_4=\frac{\beta }{\mu +s-a}$$

indicating that if $i_4<1$, a fraction of the healthy red squirrels is still present in the vanishing population. This restriction implies:

$$a<\mu +s$$

(6)

for feasibility, while:

$$a+\beta \le \mu +s$$

(7)

to have really a different outcome than $P_3$.

For $P_5$, we find:

$$N_5=\frac{K\lambda (s+\mu -\lambda )}{s\mu },i_5=1-\frac{\mu }{\lambda }.$$

Feasibility requires then:

$$s+\mu \ge \lambda \ge \mu .$$

(8)

The possible coexistence equilibrium is investigated numerically, both for feasibility as well as for stability.

3.1.2. Stability

The Jacobian $J=\left(J_{k\ell }\right)$ has the following entries:

$$\begin{array}{c}\hfill J_{11}=r\left(1-\frac{2G}{E}\right)-\frac{ck{N}^2(1-i)}{{(G+cN)}^2},J_{12}=-\frac{k{G}^2(1-i)}{{(G+cN)}^2},J_{13}=\frac{kGN}{G+cN},J_{21}=-\frac{aN(1-i)}{{(G+cN)}^2},\\ \hfill J_{22}=s(1-i)\left(1-\frac{2N(1-i)}{K}\right)-\frac{a{G}^2(1-i)}{{(G+cN)}^2}-\mu i,J_{23}=-sN\left(1-\frac{2N(1-i)}{K}\right)+\frac{aGN}{G+cN}-\mu N,\\ \hfill J_{31}=(1-i)N\left[\frac{\beta -\lambda i}{{(G+N)}^2}+\frac{aic}{{(G+cN)}^2}\right],J_{32}=(1-i)\left[\frac{(\lambda i-\beta )G}{{(G+N)}^2}-\frac{iacG}{{(G+cN)}^2}+\frac{is}{K}(1-i)\right],\\ \hfill J_{33}=i\left[\mu +s\left(1-\frac{N(1-i)}{K}\right)-\frac{aG}{G+cN}\right]-\frac{\lambda Ni+\beta G}{G+N}+(1-i)\left[\frac{\lambda N}{G+N}-(\mu +s)+\frac{sN}{K}(1-2i)+\frac{aG}{G+cN}\right].\end{array}$$

Again, although the points $P_0$ and $P_2$ are not equilibria, in their neighborhood, the trajectories are repulsed away from them, because the dominant terms are the linear ones, and they depend on the positive reproduction rates $r,s>0$, so that the populations will rebound if reduced to very low terms.

We now evaluate J at each equilibrium, to assess their local stability. $P_1$ is conditionally stable; one eigenvalue is $-s<0$, and the remaining ones provide its stability conditions:

$$\lambda <\mu ,cr<k.$$

(9)

$P_3$ has two negative eigenvalues, $-r$ and $-\mu$, while the third one provides the stability condition:

$$\mu +s<\beta +a.$$

(10)

There is an obvious transcritical bifurcation between $P_3$ and $P_4$; compare (7) and (10).

Furthermore, for $P_4$, the eigenvalues are explicit, $-r$, $s-a-\beta$ and $\beta -(\mu +s-a)$, leading to the following stability conditions:

$$\begin{array}{c}\hfill s<a+\beta <\mu +s.\end{array}$$

(11)

At $P_5$, again, one eigenvalue is explicit, giving:

$$\begin{array}{c}\hfill cr\lambda <k\mu \end{array}$$

(12)

for stability. Other inequalities come from the Routh–Hurwitz conditions on the remaining minor $\widehat{J}\left(P_5\right)$. The trace condition is satisfied, since the diagonal entries are negative in view of the feasibility conditions (8):

$$\begin{array}{c}\hfill J_{22}=(s+\mu -\lambda )\left(\frac{\mu }{\lambda }-2\right)<0,J_{33}=\frac{{\mu }^2}{{\lambda }^2}(\lambda -s-\mu )<0\end{array}$$

while the determinant ultimately gives the stability condition:

$$\begin{array}{c}\hfill 4\lambda {\mu }^2+2s\lambda \mu +2{\lambda }^3<5{\lambda }^2\mu +{\mu }^3+s{\mu }^2+s{\lambda }^2\end{array}$$

(13)

Remark 1.

The equilibria $P_1$, $P_3$ and $P_4$ have real eigenvalues, so that no Hopf bifurcation is possible in their neighborhood.

Proposition 1.

Taking λ as the bifurcation parameter, there exists a Hopf bifurcation at $P_5$, for the critical threshold:

$${\lambda }_{†}=s+\mu .$$

Proof. 

Note that the trace is:

$$\mathrm{tr}\left(\widehat{J}\left(P_5\right)\right)=(s+\mu -\lambda )\varphi \left(\lambda \right),\varphi \left(\lambda \right)=\frac{1}{{\lambda }^2}\left({\lambda }^2-\mu \lambda +{\mu }^2\right)>0,\forall \lambda .$$

Now, for ${\lambda }_{†}$, the trace vanishes, while for the determinant, we find:

$$det\left(\widehat{J}\left(P_5\right)\right)=-{(s+\mu -{\lambda }_{†})}^2\left(\frac{\mu }{\lambda }-2\right)\frac{{\mu }^2}{{\lambda }_{†}^2}-N_5(s+\mu -2{\lambda }_{†})\frac{s{i}_5}{K}{(1-i_5)}^2=N_5{\lambda }_{†}\frac{s{i}_5}{K}{(1-i_5)}^2>0.$$

Hence, the eigenvalues are pure imaginary. ☐

3.2. Simulations

In this section, the ecosystem behavior is explored by means of simulations carried out by standard MATLAB routines for the integration of ordinary differential equations.

The first results are shown in Figure 1. It is clear that the red squirrels disappear with all individuals contracting the epidemics. The role of the mortality appears only in the speed of disappearance; the faster it is, the higher the mortality rate. For this, compare the two frames, where the right one shows the simulations for a value of the mortality, which is one tenth of the one used in the left frame. The results shown in the left frame of Figure 2 indicate nevertheless that even a one thousand-fold reduction in this mortality does not prevent extinction, but as can also be seen by a close look at the prevalence graphs in Figure 1, favor a faster epidemics propagation, because the disease affects the whole population in a much shorter time. A ten-fold increase in this mortality rate prevents the disease from reaching every individual of the red squirrel population, but it clearly favors the vanishing of this population; see the right frame of Figure 2. However, this can hardly be of any help because the native squirrel population cannot be preserved anyway. This situation is known in communities, where the optimum for an (infected) individual, namely its life span, becomes a danger for the community, because more individuals would be infected by its presence. Conversely, the worst for the single individual is in this case optimal from the community point of view: by being removed, it cannot further spread the epidemics.

In all these situations, the grey squirrels are unaffected, as should be expected. Indeed, they might suffer from the red squirrels’ competition, but as the latter disappear, after a short transient, they reach their natural carrying capacity.

We then investigate the sensitivity of the system to the choice of the ecosystem parameters, keeping all the parameters that do not vary as given in (3), (4) and in (5) when not otherwise specifically stated.

For the arbitrary parameters c and $\mu$, we use a range that encompasses the conditions (4). Figure 3 shows that if the red squirrels have much more weight in the interspecific competition (small c), for moderate to large values of the natural plus disease-related mortality, they could outcompete the invading species, at the same time eradicating the disease, as well. However, this occurs for a value of c perhaps too small. The next question is to evaluate how the system responds to changes in the reproductive rate of grey squirrels, so we compare again the range of c with a suitable range of r around the value given in (5). However, the red squirrels vanish altogether with the disease affecting the whole population, over the entire range of the parameter space. The results are not shown. Combining c with the disease intraspecific transmission rate $\lambda$, again, a small stripe of low parameter values for c and moderate for $\lambda$ ensure that the native species thrives with low endemicity and invading species almost vanishing, the results being shown in Figure 4.

The next question is to evaluate how the system responds to changes in the reproductive rate of grey squirrels. The simulations show that no change in the behavior occurs; letting r vary as indicated in (5) and taking $\mu \in [0.1,10.0]$, no substantial changes occur in the ecosystem behavior. To better explore this situation, we explore also what happens if a reduction in the reproductive rates is implemented in the wild. Figure 5 shows the continuation plot when $r\in [0.01,1.01]$, for the same range of $\mu$. If artificially, the growth rate of the grey squirrel is controlled, there might be hope to restore the pristine situation, for a suitably medium value of the natural-plus-disease-related mortality.

Simultaneous perturbations of the two carrying capacities do not lead to any change whatsoever. No changes are observed either in the $K-\lambda$ parameter space; native squirrels are wiped out, all carrying the disease in the process. Essentially, the same occurs in the $K-c$ and in $K-\mu$ parameter space; in the latter, only an extremely narrow stripe, surprisingly for low K, ensures red squirrel survival, but with the disease infecting the whole population; otherwise, the red squirrels disappear anyway even by increasing their own carrying capacity; see Figure 6.

In the $E-r$ parameter space, no changes in the ecosystem behavior are observed; for the other parameters given in (3), the invasion of the grey squirrels is successful.

Letting the grey squirrel carrying capacity E vary as suggested in (5), coupling it once again with $\mu$ leads to almost a direct linear relationship with the grey squirrel population size, this time independent of $\mu$. The results (not shown) indicate that the red squirrel population is wiped out, and the prevalence is heavily affected by $\mu$ and almost independent of changes in E. Furthermore, in the $\beta -r$ parameter space, independent of the grey squirrel reproduction rate r, a small transmission rate induces a reduction in the prevalence, but the red squirrels are wiped out anyway. Combining the interspecific transmission rate $\beta$ with K, the same effect occurs: for low values of $\beta$, the disease prevalence drops, independent of the value of the red squirrels’ carrying capacity, but the latter do not survive anyway.

We now investigate the combined changes of both reproductive rate and carrying capacity of grey squirrels, at first with a low value of the natural plus disease-related mortality, $\mu =0.2$. In this case, the red squirrels are wiped out with an epidemic that in the process affects the whole population. The grey squirrel population instead grows linearly with its carrying capacity, almost independently of the reproduction rate (result not shown).

Taking a larger intermediate value, $\mu =5.0$, instead shows that there is again a narrow stripe in the parameter space in which the red squirrels thrive. It is obtained for low values of the grey squirrel reproduction rate, $r\le 0.2$ approximately, while the grey squirrel population increases linearly with E; see Figure 7. We finally take $\mu =10.0$. The behavior is similar as for $\mu =5.0$, but now for small values of E and moderate values of r, the red squirrels survive, at a similar level as before; compare Figure 7 and Figure 8, the latter corresponding to this simulation. In this case, the disease is also heavily affected, since it is eradicated for very low r, but also when it is endemic, its prevalence is very much reduced; note indeed that the vertical scale for i in Figure 8 is less than half of the corresponding one in Figure 7. The grey squirrel population, when surviving, is hardly influenced by the change of r.

We consider also changes in the respective competitive rates, a and k, but in the present conditions (3) and (4), no significant changes are observed; the red squirrels are wiped out, although the disease does not affect the whole population.

We finally tried to see how the influence of the disease transmission is felt by the ecosystem. By comparing the respective influences of the two transmission rates, trying to lessen the interspecific transmission rate $\beta$ appears to induce more important changes in the ecosystem than reducing the intraspecific one, $\lambda$.

Furthermore, a drastic reduction of the interspecific transmission rate $\beta$ combined with higher natural-plus-disease-related mortality $\mu$ cannot help in preserving the native population. The analysis of the combination of $\lambda$ and $\mu$ does not lead to native squirrels’ survival, as should be expected from the previous above results on the intraspecific disease transmission rate, but their prevalence decreases as the mortality increases.

Decreasing r, E and $\beta$, while increasing $\mu$ from the reference values (3) and (4), namely taking $r=0.3$, $E=0.5$, $\mu =10.0$, $\beta =0.6$, the red squirrel population is preserved with the epidemics disappearing, while the grey squirrels are also eliminated; see the left frame of Figure 9.

Increasing r alone, in these same conditions, however, does not lead to coexistence, as past a value of the invasive species reproduction rate that lies in the range $[0.57,0.58]$, the grey squirrels instead eliminate the native species; see the right frame of Figure 9.

We attempt at last multiple changes of the parameters, taking new values as follows:

$$r=0.5,a=0.4,\mu =1.5,\beta =0.4.$$

(14)

For these values, the coexistence of the three populations occurs with the disease not affecting the whole native population; see the left frame of Figure 10. Therefore, for our in silico experiments, from now on, we use the parameters (14). In this four-dimensional parameter space, we study once more the sensitivity. At first, we use the pair $E\in [0.1,0.8]$ and $\beta \in [0.1,1.0]$ that has not yet been considered; the results are not shown. When both parameters attain their lowest value in these respective ranges, the red squirrel population survives, with low endemicity, while the grey squirrels’ competitors are sensibly kept in check. We also once more compare the coexistence as a function of $E\in [0.1,0.8]$ and the mortality $\mu \in [8,12]$; see Figure 11. Coexistence occurs for lower values of E with a thriving native population at moderate levels, low prevalence and the grey squirrels kept in check at reasonably low values. Furthermore, this result appears to be almost independent of the values of the mortality. The most relevant parameter for obtaining this outcome is the interspecific disease transmission rate $\beta$, because if it is low, the native population thrives also for higher values of the grey squirrel carrying capacity.

Next, we compare the ecosystem outcome as a function once more of the grey squirrel reproduction rate r and the competition a suffered by the native population. The native population survives less for increasing competition, but better if r is very small, just in which case, the invaders suffer greatly; the disease prevalence is confined to reasonably low values for intermediate values of the parameters, in the range explored (results not shown).

As functions of the grey squirrels reproduction rate and interspecific disease transmission rate, the populations survive in an intermediate range of this parameter space, with declining native population as the incidence rate grows. The grey squirrels suffer from a very own low reproduction rate, as should be expected, but in this range, while the red squirrel population is also much reduced or vanishes altogether, persistent prevalence oscillations are observed, indicating that the disease remains endemic while the population disappears; see the right frame of Figure 10.

The ecosystem behavior in terms of $\mu$ and $\beta$ shows that the role of mortality is less important; the red squirrels are preserved independently of it if the transmission rate is low enough. Prevalence drops with higher mortality though, the invasive population picks up instead with higher transmission rates.

Low values of r and moderate values of the mortality ensure the disappearance of the invaders and the survival of the red squirrels. Some oscillations for the prevalence are detected in the same range. A simultaneous reduction in both a and $\beta$ leads to coexistence with a somewhat reduced grey squirrel population, while the native species thrives at good levels, and the disease is sensibly reduced.

4. Discussion

Our analyses revealed that the three-populations system has different equilibria. From a conservation and wildlife management perspective, the most desirable one is the healthy-red-squirrels-only point $P_1=(0,K,0)$, followed by $P_5=(0,N_5,i_5)$. Apart from the origin, $P_0$, these are the only equilibria leading to the grey squirrel extinction, and it is thus particularly interesting to evaluate their conditions. Equilibrium $P_1$ is always feasible, and stability can be attained only if conditions (9) are satisfied, while $P_5$ is feasible and stable if (8), (12) and (13) hold (see

Table 2. Summary of the feasibility and stability conditions for the equilibria with at least one vanishing population.

Equilibrium	Feasibility	Stability
$P_1=(0,N_1,0)$	always	$\lambda <\mu$, $cr<k$
$P_3=(G_3,0,i_3)$	$i_3=1$	$\mu +s<\beta +a$
$P_4=(G_4,0,i_4)$	$a<\mu +s$	$s<a+\beta <\mu +s$
$i_4\le 1$	$a+\beta \le \mu +s$
$P_5=(0,N_5,i_5)$, $i_5<1$	$s+\mu \ge \lambda \ge \mu$	$cr\lambda <k\mu$
		$4\lambda {\mu }^2+2s\lambda \mu +3{\lambda }^3<5{\lambda }^2\mu +{\mu }^3+s{\mu }^2+s{\lambda }^2$

). For both points, the conditions often involve the intraspecific disease transmission rate for red squirrels ($\lambda$) and the overall red squirrel mortality ($\mu$), whose combination did not lead to native squirrel survival according to our simulations. Using the reference parameter values (3), the system is rather likely to evolve from a two-population ecosystem towards a grey-squirrel-only situation, with the disappearance of the red squirrel population in a time span that depends on the value of the red squirrel overall mortality (Figure 1 and Figure 2). Indeed, excluding the unstable $P_2$, all other equilibria imply the persistence of the grey squirrel population, with the native species vanishing while the disease remains endemic in it. This is consistent with the results obtained through other models. For instance in [52], it was also found that, in the absence of control of the invading species, native populations are driven to extinction and that this occurs both in the absence of disease (through competition only; see also [25,43]) and more rapidly when the disease is included [28].

Investigating the sensitivity of the system with respect to the parameter values, we highlighted that the red squirrel could outcompete the alien species and at the same time eradicate the disease if it had much more weight in the interspecific competition (small c, Figure 3 and Figure 4). This is unlikely to happen in most of the habitat types where the grey squirrel successfully settles (e.g., where its carrying capacity is the highest, as in mixed and broadleaved woodlands), but it could locally occur in coniferous forests [51].

The system may also respond to changes in the growth rate (r) of the grey squirrel population. The latter is certainly the parameter that can be manipulated most easily in natural conditions. As shown in Figure 5 and Figure 6, if artificially, the growth rate of the grey squirrel is controlled, the alien species may disappear from the system, with cascade effects on the transmission of the disease to reds. The native species can thus recover and return to its natural carrying capacity. However, to obtain such a result, the reduction in r must be very harsh: the parameter value should be smaller than 0.2. This is required for all values of E (Figure 7), i.e., independent of the habitat type and of its carrying capacity for the grey squirrel. Interestingly, decreasing the growth rate (r), the grey squirrel carrying capacity (E) and the intraspecific transmission rate ($\beta$) at the same time and increasing the red mortality ($\mu$), the red squirrel population can be preserved even for $r=0.3$ (Figure 9). In this case, the grey squirrels are eliminated, and the epidemics disappears. Again, [51,52] also found that grey squirrel control can protect red populations in strongholds, although in their specific case study, periodic outbreaks of the disease could not be prevented because of occasional recolonization of the stronghold by grey squirrels from neighboring areas. The work in [52] also highlighted that there is a threshold level of control (expressed as the proportion of grey squirrels to be removed), above which the invading population can be prevented from establishing and the native species can be protected. Another study [28] also investigated what level of population control of grey squirrels would be necessary to suppress the disease-induced decline in red squirrels. They found that the kill rate required to avoid the spread of the invasive species should be in excess of 60% grey squirrels. For any level of control below this threshold, the grey squirrels were predicted to expand with a concomitant decline in red squirrel populations.

Since achieving $r\le 0.3$ could still be difficult, managers could also aim at a larger 0.5 value. Assuming that the reduction in the growth rate is achieved through the removal of the grey squirrels, the red squirrel damage due to interspecific competition (a) and the intraspecific transmission rate ($\beta$) could also be smaller (e.g., 0.4 as in (14)), and the overall red mortality ($\mu$) should also be reduced. Under these conditions, the coexistence of the three populations is ensured, and the grey squirrels are kept in check at reasonably low values. We also note that the reduction in the interspecific transmission rate ($\beta$) appears more important in the ecosystem than the intraspecific one, $\lambda$, and more important than E (Figure 11). In this respect, [25] also indicated that the progress of the disease within individual populations is largely dependent on the probability of virus transmission between the two species.

Different combinations of the parameter values may thus allow wildlife managers to achieve different objectives. Simulations highlighted that the most desirable outcome could be achieved only in favorable conditions for the red squirrel (very small c) or for very small values of r (e.g., 0.2 or 0.3). How the growth rate of grey squirrels could be reduced, in practice, depends on the specific circumstances, including the stage of the invasion. Actually, two strategies have been proposed to eradicate or at least to control grey squirrels in the wild.

The first strategy involves the removal of squirrels through capture and subsequent euthanasia. Evidence for its effectiveness is provided by data from the U.K. [46], and other mathematical models highlighted its potential to protect red populations and to prevent the spread of SQPV [27,51,52]. According to our results, the removal control approach could ensure the conservation of the native red squirrel population if r is substantially reduced ($r\le 0.5$). In this respect, it should be noted that the removal of the animals is advantageous, since it allows for a reduction in r, and in turn, it could affect other model parameters, such as a and $\beta$, in a favorable way for the conservation of red squirrels. Culling can indeed remove in parallel both the competitive and disease threats posed to red by grey squirrels [46]. Its application is also consistent with the provisions of EU Regulation No. 1143/2014. According to this regulation, rapid eradication should be now implemented when an alien species of European concern is detected in a new area. In the first phases of an invasion, the number of alien individuals should be limited, and their population should be in the lag phase that typically precedes rapid expansion and demographic growth. It should thus still be possible to affect the growth rate of the population by removing a few individuals. On the contrary, if the invasive population is large, reducing r would require a huge capture effort.

The second suggested strategy relies on the reduction of the reproductive rate, as achieved through surgical sterilization or inhibition of fertility, e.g., using immunocontraceptive vaccines. In particular, gonadotrophin-releasing hormone (GnRH) vaccines seem to offer great promise for the control of populations. They interfere with the regulation of reproductive hormones, and they have been successful at reducing fertility in several mammals [53]. Available data on their potential impact on the population dynamics of other squirrel species [54] suggest they could reduce populations, without eradicating them, if at least 70% of the females were treated [54,55,56,57]. At present, this requires the capture of the animals in order to inject the vaccine, and it could be necessary to capture and treat each animal more than once, since the long-term effects of the vaccine have not been demonstrated so far [58,59]. On the basis of a theoretical model, [56] concluded that 90% of females should be treated in order to eradicate a grey squirrel population within 7/8 years, ensuring a constant sterilization effort. The need to ensure a constant effort in reducing the reproductive rate is in agreement with the results of our model, since we found that, if the reproductive rate of the grey squirrel is artificially controlled, there might be hope to restore the pristine situation, for a suitably medium value of the natural-plus-disease-related mortality of the red squirrels. Anyway, sterilization also implies the subsequent release of the individuals in the wild, so that through this approach, we could not expect prompt effects on the model parameters, other than r. In particular, in this case, we would expect no immediate reduction in the occurrence of the disease and, thus, no major impact on the SQPV interspecific transmission rate.

In conclusion, our model confirms that the native red squirrel populations can be preserved or restored, under specific parameter conditions, if the grey squirrel population is removed or kept in check at low values. Independent of the adopted methodology, wildlife managers should aim at a harsh reduction of the grey squirrel growth rate, and a constant effort should be ensured to protect the remnant red squirrel populations, if the complete removal of the alien species is not achieved at the early stages of invasion. As concerns the SQPV, actions should mainly aim at the reduction of the interspecific transmission rate, one of the most important parameters in determining the ecosystem response to the disease, although how to effectively carry out this task in practice may still be an open question. At the moment, programs for the grey squirrel control, or even local eradication, are implemented both in the U.K. and in Italy. The insights provided by our general model could be considered by these programs to improve the planning of management actions. Although modeling outcomes have already been used by other authors to critically analyze the efficiency of removal activities (e.g., [27]), as implemented in specific contexts, future research could concentrate on the analysis of a more general system, including also the control strategy in its formulation, as is commonly used for instance in the case of epidemics [60,61,62]. This approach may provide useful indications to the ecologists that deal with this problem in practice.