Metaecoepidemic models with sound and infected prey migration

Available online at www.sciencedirect.com

ScienceDirect Mathematics and Computers in Simulation 126 (2016) 14–44 www.elsevier.com/locate/matcom

Original articles

Metaecoepidemic models with sound and infected prey migration✩ Ezio Venturino Dipartimento di Matematica “Giuseppe Peano”, Universit`a di Torino, via Carlo Alberto 10, 10123 Torino, Italy Received 21 July 2011; received in revised form 10 February 2016; accepted 23 February 2016 Available online 11 March 2016

Abstract We present a two-patch metaecoepidemic model, in which predators, sound and infected prey occupy the first one. Both kinds of prey can migrate into a second habitat, which constitutes a refuge from the predators. Equilibria of the system and their stability properties are determined analytically and also heavily studied numerically, when too cumbersome. The results indicate that in practical situations it is necessary to carefully assess the environmental state before proceeding to any modification of it, as the consequences might be far from those foreseen. c 2016 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights ⃝ reserved. Keywords: Ecosystems; Fragmentation; Diseases; Predator–prey; Ecoepidemiology

1. Background In a recent investigation, [31], the author has introduced the concept of metaecoepidemic models. These account for interacting populations affected by a disease and living in separate environments, thereby constituting independent ecoepidemic models, see Chapter 7 of [23], that become interdependent since migration among them is allowed. In [31], only sound prey are assumed to migrate between the different environments, while in [5] predators only do. Here we extend the investigation by allowing also the infected prey to migrate, while at the same time modifying the model, by assuming a Holling type I functional response for the predation instead of a Holling type II. Note that the disease is assumed to propagate by contact among the prey, but it does not affect the predators. The rationale behind this new concept lies in the investigations on metapopulations, performed since a number of years to study the effects that landscape fragmentation has on ecological communities, [32,33], since the latter may threaten the species survival. In some cases however, persistence in the whole environment can be obtained, although locally populations may experience extinction [8,13,15,32,34]. Heterogeneous environments may arise naturally due to physical unfavorable events, such as landslides, floods and so on, but also human activity is frequently responsible for substantially altering the ecosystems. Construction of roads and human artifacts is the first example, but most agricultural activities cause even more damage to the natural landscape. As a consequence, populations living in the unperturbed environment find themselves partitioned, with one ✩ Part of this research was presented at the DIEMBM 2010 Conference in Samos, Greece; the paper is dedicated to the organizer of the meeting, Athena Makroglou. E-mail address: [email protected].

http://dx.doi.org/10.1016/j.matcom.2016.02.006 c 2016 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights 0378-4754/⃝ reserved.

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

15

subpopulation being cut off from their similar that lie on the other side of the newly arisen barriers, and therefore start to evolve independently of each other, but may become also less apt to survive in adverse conditions. Metapopulations represent a way to investigate the possible evolution in time of such situations, [19]. They have been successfully applied to study a few species, the spotted owl (Strix occidentalis), the mountain sheep (Ovis canadensis), [14]. In other cases however, the concept of incidence function is used instead, see the case of the butterflies Melitaea cinxia and other species in Finland, [16,26,27]. But it is very difficult to collect real data on these models, as remarked in the literature of metapopulation theory, [8,13], so that the importance of mathematical models in this context becomes clear. Further, ecoepidemic models deal with interacting populations that are affected by a disease. Since the latter occur in nature, they cannot be neglected when studying the dynamics of populations. Furthermore, they must be considered also in metapopulation models. Let us mention in this context the case of mixomatosis, which was humanly inoculated in the wild rabbit, Oryctolagus cuniculus (L.), in Australia to try to control its increasing population, [20,21,25,29]. The rabbits are hunted by the red fox Vulpes vulpes (L.). Further examples are provided by Ovis canadensis, Strix occidentalis and butterflies as mentioned above. The predators of the former are mainly wolf (Canis lupus), coyote (Canis latrans), bear (Ursus), Canada lynx (Lynx canadensis), mountain lion (Puma concolor), golden eagle (Aquila chrysaetos), [9]. Their parasites are mainly nematode lungworms, Protostrongylus stilesi, P. rushi, Cysticercus tenuicollis, Wyominia tetoni, Nematodirus archari, N. davtiani, Trichostrongylus sp., Protostrongylus rushi; Dermacentor albipictus and D. venustus, [4,10]. The principal parasites of Strix occidentalis are instead the helminths, which are assumed by predation on small rodents, while the most important predator is the great horned owl, Bubo virginianus, which hunts the spotted owl to feed on it but also to eliminate it as a possible competitor for resources, [18]. Lepidoptera are hunted by birds, bats, small mammals, reptiles, and their larvae become the main diet of Parus caerulens, P. major, [7]. Butterflies defend themselves by several means: mimetism at the moth stage, diurnal butterflies assume compounds that make them unpalatable for their predators, e.g. Melitaea cinxia feeding on Cotesia melitaearum, the latter getting alkaloids from the plant on which it feeds, [22]. Other defense mechanisms are shapes, colors or sounds which discourage the predators, [28]. Myrmecophily is an association with ants in which caterpillars find refuge near the ant nests, since they are a predator- and parasitoid-free space, [3]. Parasitoids attack Lepidoptera by implanting in them an egg, which later on kills the host, [12,28]. Lepidoptera are also affected by viruses, for instance nuclear polyhedrosis virus (NPV), cytoplasmatic polyhedrosis virus (CPV), granulosis virus (GV), entomopox virus (EPV), and common bacteria in nature (Bacillaceae), such as Bacillus thuringiensis var. kurstaki, which upon ingestion poisons the caterpillars’ digestive system. The previous examples illustrate the necessity of investigating diseases also in the context of metapopulations, extending therefore the concept of ecoepidemics to this wider situation. A first attempt has been made in [31]. With respect to the former study, here we modify the picture by allowing also infected prey to move freely among the patches. Mathematically, this leads to a more complicated system, containing five nonlinear differential equations. The paper is organized as follows. We introduce the model in the next section. In Section 3 we adimensionalize it and investigate its limiting cases, together with the two systems that possibly constitute the independent patches of the model if communications are interrupted. Section 4 contains the analysis of the simplest possible system’s equilibria. The study of the most complicated ones is postponed to the simulations Section 5. The investigation of the effects of parameter changes in the system’s dynamics appears in Section 6 and a final discussion of the results concludes the paper. 2. The model We consider here a two-habitat environment in which two populations coexist, prey and predators. Furthermore, a disease spreads among the prey. In the first patch, both species are present, while in the second one, only the prey, so that the latter can be considered as their refuge, inaccessible by the predators, see Fig. 1. Migration is allowed between the two patches, both for sound and infected prey subpopulations. The following features on the disease are further assumed. It is transmitted by contact, and it is recoverable. Infected individuals are so weak that they are not able to reproduce, nor to compete for resources with the sound ones, and among themselves. Both disease incidence and recovery rates are independent of the patch. The disease incidence is taken in the simple mass action form.

16

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

Fig. 1. The scheme of the system. Thick arrows denote predation, thin arrows instead denote the infection process. The arrows outside the patches denote instead migration.

Sound prey reproduce logistically, with different carrying capacities and reproduction rates in the two patches, while the predators are assumed not to have any other source of food. They however feed indiscriminately upon sound and infected prey, and from the latter they receive no harm. The patch common to all populations is denoted by the index 1 while the index 2 is used for the prey refuge unreacheable by the predators. We distinguish the two subpopulations living in the two patches by suitable subscripts. The new features here with respect to [31] are firstly the fact that both kinds of prey are allowed to migrate, not just the healthy ones, therefore the disease does not constitute an obstacle for movement between patches. The second important change is that here we do not introduce a Michaelis–Menten term to model the predators feeding behavior, but rather assume a Holling type I response function. This function better models situations in which intermingling of the populations allows contacts between all the individuals. This is reasonable if the population sizes are not large, which we assume in this context. Finally, also the intrinsic structure of the system is different, since the disease in the prey is present here in both patches, while in the first model of [31] it was confined in the patch inaccessible to the predators. These aspects therefore make a direct comparison of the systems’ results impossible in this case. However, to discuss the changes due to the predators’ migration, we will study also the corresponding model when only healthy prey migrate and no epidemic is present in the first patch, Section 5.1. We also study the differences in system behavior with the last model of [31], in Section 5.2. Let  S(τ ) denote the sound prey,  I (τ ) the infected ones  ) the predators, where τ denotes time. The model reads then and P(τ    d S1 S 1 − m 21 12 − λ S1  I1 +  ν I1 −  a S1 P S1 + m S2 , (1) = r1 S1 1 − 1 dτ K d I1 − = I1 [ λ S1 −  µ1 −  ν − cP n 21 ] +  n 12  I2 , dτ    dP  =P e S1 +  g I1 −  b , dτ    d S2 S 2 21 12 = r2 S2 1 − − λ S2  I2 +  ν I2 + m S1 − m S2 , 2 dτ K d I2 = I2 [ λ S2 −  µ2 −  ν − n 12 ] +  n 21  I1 . dτ Since we count populations as pure numbers the units of the above original parameters in the model are frequencies. The evolution of the prey population in the first patch is represented by the first equation, accounting for logistic growth, infection and recovery processes, predation and emigration and immigration processes. Note that in the logistic correction term in the sound prey evolution equation, no contribution from the infected is present, since infectives do not contribute to intraspecific competition, so that sound individuals do not feel their presence. The pa1 the carrying capacity in patch 1; rameters have the following interpretation:  r1 is the prey net reproduction rate and K

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

17

the infection process occurs at rate  λ, its recovery at rate  ν; the predation rate is represented by  a , and migration 21 and m 12 denotes the rate in the opposite direction. Similarly the fourth equation defrom patch 1 occurs at rate m scribes the sound prey dynamics in the second patch, with obvious changes in notation and with the exception of the predation. The second equation describes the infected in the first patch. Their number increases only via the immigration process at rate  n 12 , expressed by the last term, or by new “successful” contacts between an infective and a sound individual; infected leave this class either via natural plus disease-related mortality  µ1 or by recovering at rate  ν, they are also hunted by predators, at a different rate  c than for sound individuals, and finally they emigrate at rate  n 21 . Note that these rates differ from those of the sound individuals, since the sickness may prevent the affected individuals to move freely as their sound fellows do. Again, the fifth equation expresses the same dynamics in the second patch, with the exception of the hunting process. The third equation describes the predators dynamics, accounting for the reward obtained by hunting both sound and infected prey, although at different rates,  e and  g respectively, and natural mortality  b. The different benefit predators have from hunting may be due to the fact that infected might be easier to capture, since they are sick, or alternatively also to the fact that they may be less palatable than sound ones. Clearly, the type of chosen response functions changes the final behavior of the system, but as already remarked in [5], even the quadratic nonlinearities may lead to limit cycles, which in standard predator–prey models arise only if the Holling type II response function is assumed, or alternatively other types of nonlinearities, like the very recently introduced square root interaction terms modeling herd behavior, [1]. 3. Preliminary remarks 3.1. Model adimensionalization We rescale the model via the following substitutions, for k = 1, 2 Sk = αk  Sk ,

Ik = γk  Ik ,

 P = β P,

t = ητ,

where now Sk = Sk (t), Ik = Ik (t), P = P(t). By choosing √  a c 1 , β= , η = ν, α1 = α2 = γ1 = γ2 = 1  ν K we can keep the disease incidence and recovery the same in the two patches. Defining the new parameters, for j = 1, 2  1   12  21  rj λK m a m n 12 , m 12 = σ = , h= , m 21 = , n 12 = , rj = ,   ν  ν c  ν  ν  ν 1 1  2     µj n 21 eK gK b K n 21 = , θj = , e= , g= , b= , k= , 1  ν  ν  ν  ν  ν K the simplified system with only the above 14 parameters takes the final form d S1 dt d I1 dt dP dt d S2 dt

= r1 S1 (1 − S1 ) − σ S1 I1 + I1 − h S1 P − m 21 S1 + m 12 S2 , = I1 [σ S1 − θ1 − 1 − h −1 P − n 21 ] + n 12 I2 , = P (eS1 + g I1 − b) ,   S2 = r2 S2 1 − − σ S2 I2 + I2 + m 21 S1 − m 12 S2 , k

d I2 = I2 [σ S2 − θ2 − 1 − n 12 ] + n 21 I1 . dt

(2)

18

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

Its Jacobian is given by  J11 1 − σ S1  σ I1 J22  J = e P g P  m 21 0 0 n 12

−h S1 −h −1 I1 J33 0 0

m 12 0 0 J44 σ I2

 0 n 12   0   1 − σ S2  J55

(3)

with J11 = r1 (1 − 2S1 ) − σ I1 − h P − m 21 , J22 = σ S1 − 1 − h −1 P − n 21 − θ1 ,   2 J33 = eS1 + g I1 − b, J44 = r2 1 − S2 − σ I2 − m 12 , k J55 = σ S2 − 1 − n 12 − θ2 .

(4)

3.2. Limit cases for the populations First of all, under suitable assumptions, the system does not blow up, so it is biologically meaningful. To see this, let us define the total ecosystem population T = S1 + I1 + P + S2 + I2 . The technique is standard, we can refer for instance to [17], but we reemphasize here the main steps for the benefit of the reader. By adding all the equations of (2) for an arbitrary 0 < α we obtain dT + αT = (r1 + α)S1 − r1 S12 + (e − h)S1 P + (g − h −1 )I1 P + (α − b)P dt r2 + (r2 + α)S2 − S22 + (α − θ1 )I1 + (α − θ2 )I2 . k Now taking α < min{b, θ1 , θ2 }, and assuming e < h, g < h −1 , taking the maxima M( j) , j = 1, 2 of the two parabolae in S1 and S2 which appear on the right hand side, namely M( j) =

(r j + α)2 , 4r j

for M∗ = M(1) + M(2) we have the differential inequality dT + αT ≤ M∗ dt from which for some constant C it follows T (t) ≤ M∗ α −1 + C exp(−αt), i.e. T (t) ≤ max{M∗ α −1 , T (0)}. Since the sum of the populations is bounded above and each one of them is nonnegative, it follows that each subpopulation cannot grow to infinity. In particular note that to achieve this result it is sufficient that the (adimensionalized) predation rate on susceptible prey exceeds the return predators obtain from hunting on sound prey, which is a quite natural requirement to ask, as not the whole prey is consumed, and at the same time it is bounded above by the reciprocal of the predation rate on infected prey. The interval [e, g −1 ] is in fact not empty, since presumably the gain from eating sick prey is less than the one obtained from consumption of sound individuals. Indeed g < e < 1 implies ge < 1 from which the result. Secondly, we ask the question whether and under which circumstances the ecosystem may disappear. The system equilibrium corresponding to this result is clearly the origin, which is feasible, in view of the fact that (2) is a homogeneous system. Now, to have ecosystem disappearance, the origin must be a stable equilibrium. Upon evaluation of (3) one negative eigenvalue is immediately found, −b. The remaining minor factorizes, to give the two quadratic equations, λ2 − λ[r1 + r2 − m 12 − m 21 ] + (r1 − m 21 )(r2 − m 12 ) − m 12 m 21 = 0, λ2 + λ[2 + n 21 + θ1 + n 12 + θ2 ] + (1 + n 21 + θ1 )(1 + n 12 + θ2 ) − n 12 n 21 = 0.

19

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

Now the latter has only roots with negative real part, while the Routh–Hurwitz conditions applied to the former give the following necessary and sufficient stability conditions for stability r1 + r2 < m 12 + m 21 ,

r1r2 > r1 m 12 + r2 m 21 .

(5)

The former of the above conditions shows that to maintain this ecosystem alive, it is sufficient that the total prey reproduction rate exceeds the total migration rate, i.e. the sum of the fluxes in the two directions. This can be expressed in other terms by introducing a new parameter M, the total “activity” in the system. When this activity exceeds a threshold, M > 1 then the ecosystem may collapse. Define indeed M=

m 12 + m 21 , r1 + r2

M (1) =

m 21 , r1

M (2) =

m 12 . r2

The latter two quantities represent instead the “activity” in each patch, namely the ratio between newborns coming to the light in the patch and outgoing migration. Remark 1. The second (5) can then be restated as M (1) + M (2) < 1. Again, to keep the ecosystem going, it is then sufficient that in one of the patches the emigration activity exceeds the local reproduction rate, so as to have for some j, M ( j) > 1. Remark 2. From (5) it is also clear that pure imaginary eigenvalues will appear, but they cannot possibly give rise to Hopf bifurcations, because when M crosses the value 1 at the same time M (1) + M (2) < 1 must hold. However, for M = 1 we get r1 = m 12 + m 21 − r2 and substitution into the other equation gives 0 > r22 + m 212 + m 21 m 12 , which cannot be satisfied. 3.3. The independent patches In the model there are two patches which if communications are closed, would develop independently. In patch 1 we would find a system very close to the one studied in [30], while in patch 2 there is an almost standard epidemic model, with the exception that it does exhibit some demographic behavior contrary to the classical models. However, such types of systems have been investigated quite some time ago, [6,11,24]. We briefly summarize the two submodels behavior, since there might be slight differences with the models of the literature. In patch 1, the system is given by the first three equations of (2) in which the migration terms are excluded. Its Jacobian is the top principal minor of order three of (3). This ecoepidemic system allows the following conditionally stable equilibria in the S1 , I1 , P phase space, the origin being always unstable in view of the eigenvalues r1 , −(θ1 +1), −b. The points E 1 = (1, 0, 0), always feasible, the coexistence equilibrium E 4 and      b r1 b θ1 + 1 r1 E2 = , 0, 1− , E3 = , 2 (θ1 + 1) (σ − θ1 − 1) , 0 . e h e σ σ θ1 Let (1)

R0 ≡

σ , 1 + θ1

RP ≡

e . b

(6)

The eigenvalues of E 1 are always real, −r1 , e − b and σ − θ1 − 1, so that no Hopf bifurcations can arise. E 1 is (1) stable for R P < 1 and R0 < 1. E 2 is feasible for R P ≥ 1. Thus there is a transcritical bifurcation for which E 2 emanates from E 1 when R P = 1. The eigenvalues of E 2 are σ S1 − θ1 − h −1 P and the roots of the quadratic λ2 + λr1 S1 + eh P S1 = 0. Since all its coefficients are positive, the roots have negative real parts. Furthermore r1 S1 > 0 unless the reproduction rate vanishes or the susceptible population collapses, so that in this case no Hopf bifurcations can occur. Finally, E 2 is stable when the first eigenvalue is negative, i.e. for b 1 < [h 2 (θ1 + 1) + r1 ]. 2 e σ h + r1 (1)

(1)

E 3 instead is feasible for R0 ≥ 1. Here another transcritical bifurcation occurs, E 3 arises from E 1 when R0 = 1. E 3 has the eigenvalues eS1 + g I1 − b and the roots of λ2 + λ(2r1 S1 + I1 σ − r1 ) + σ I1 (σ S1 − 1) = 0. However,

20

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

note that σ S1 − 1 = θ1 > 0 and 2r1 S1 + I1 σ − r1 = r1 S1 + I1 S1−1 > 0, so again the corresponding eigenvalues have negative real part. The latter inequality is strict for S1 > 0, I1 > 0. Thus Hopf bifurcations never occur here as well. E 3 is thus stable for   (1) (1) 2 gr1 − (eθ1 + gr1 )R0 + bθ1 R0 > 0. For the coexistence equilibrium E 4 the population levels are I1 =

S1 [r1 (1 − S1 ) − h P], σ S1 − 1

P = h(σ S1 − θ1 − 1),

and S1 solution of the quadratic AS12 + B S1 + b = 0,

A = eσ − g(r1 + h 2 σ ), B = gr1 + gh 2 (1 + θ1 ) − e − bσ.

Independently of the sign of B, for A < 0 there is only one positive root, for A > 0 and B < 0 there are two. The remaining feasibility conditions are (1)

R0 S1 ≥ 1,

r1 ≥ r1 S1 + h P.

In patch 2, the system is given by the last two equations of (2) with no migration terms. Its Jacobian is the bottom principal minor of order two of (3). Apart from the origin that is always unstable, since the eigenvalues are −1 − θ2 and r2 , we find the following equilibria,   θ2 + 1 r2 (θ2 + 1) kσ − θ2 − 1 , . Q 1 = (k, 0), Q2 = σ σ θ2 kσ Now Q 1 has eigenvalues −r2 and σ k − θ2 − 1, so that it is stable for (2)

R0 =

kσ < 1. θ2 + 1

(7)

(2)

Q 2 is feasible for R0 ≥ 1, indicating thereby a transcritical bifurcation between Q 1 and Q 2 . The Routh–Hurwitz conditions for the latter are σ I2 (σ S2 − 1) > 0, which is trivially true since it reduces to θ2 + 1 > 1 and 2r2 k −1 S2 + σ I2 − r2 > 0, from which σ k + θ22 > 1, a condition that always holds as a consequence of feasibility, namely σ k + θ22 ≥ kσ ≥ θ2 + 1 ≥ 1. Here the equality would hold whenever θ2 = 0, but in such case the equilibrium does not exist, since except for a particular choice of the parameters, namely kσ = 1, leading to the line of equilibria (σ −1 , I¯2 ) with I¯2 arbitrary. Thus no Hopf bifurcation arises here as well. Further, Q 2 when feasible is always stable. Needless to say, in view of these results and in particular of the transcritical bifurcation, whenever locally asymptotically stable, these two equilibria are also globally asymptotically stable. 4. Possible system configurations 4.1. Coexistence only of the sound prey At this equilibrium, E C = (S1C , 0, 0, S2C , 0), the levels of the two sound prey populations are seen to come from the positive intersections of the two parabolae     S2 1 1 χ (S1 ) = S1 [r1 (S1 − 1) + m 21 ], Π (S2 ) = S2 r2 − 1 + m 12 . m 12 m 21 k They always occur whenever at least one of their zeros is positive, while when they are negative, observing their derivatives at the origin, it is seen that no feasible intersection can exist. The feasibility conditions become, specifically, the three following alternatives M (1) < 1,

M (2) < 1;

M (1) > 1,

M (2) < 1;

M (1) < 1,

M (2) > 1.

Note that the last two conditions are incompatible with the stability of the origin, see indeed the remark after (5).

(8)

21

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

One eigenvalue is explicitly calculated as eS1C − b, the remaining ones are the roots of the two quadratic equations, referring to (4) with the notation emphasizing the quantities evaluated at this point, C C C C λ2 − λ(J11 + J44 ) + J11 J44 − m 12 m 21 = 0,

C C C C λ2 − λ(J22 + J55 ) + J22 J55 − n 12 n 21 = 0.

The Routh–Hurwitz conditions can be used, so that the stability conditions are explicitly S1C < 1>

b , e

M (1) 1 − 2S1C n 21

σ S1C − 1 − θ1

+

+

M (2) 1 − 2k −1 S2C n 12

σ S2C − 1 − θ2

< 1 < 2(ρ1 S1C + ρ2 k −1 S2C ) + M, ,

(9)

σ (S1C + S2C ) < 2 + n 12 + θ1 + n 21 + θ2 ,

where for j = 1, 2 we have set rj ρj = . r1 + r2 In this case it is also seen that in principle Hopf bifurcations could be possible, imposing either one of the conditions 1 = 2(ρ1 S1C + ρ2 S2C ) + M,

(10)

σ (S1C + S2C ) = 2 + n 12 + θ1 + n 21 + θ2 ,

(11)

with the remaining Routh–Hurwitz conditions (9) still holding. 4.2. The disease-free state At this equilibrium both sound prey and predators survive, namely E D = (S1D , 0, P, S2D , 0). We have explicitly     1 b b m 12 e D+ D D S1 = , P = S , r1 1 − − m 21 + e h e bh 2    k 4 D± r2 − m 12 ± (r2 − m 12 )2 + r2 m 21 b . S2 = 2r2 ek Only S2D+ can be feasible, and for its feasibility we need S2D+ ≥

b [e(m 21 − r1 ) + r1 b]. e2 m 12

(12)

The Jacobian at this point factorizes into a cubic and a quadratic. The Routh–Hurwitz stability conditions of the latter become D D D D J22 + J55 < 0, J22 J55 > n 12 n 21 . 3 The cubic i=0 ai λi = 0 has coefficients

(13)

D D a2 = b − eS1D − J44 − J11 , D D D D a1 = eh S1D P D + (eS1D − b)(J44 + J11 ) + J11 J44 − m 12 m 21 , D D D J11 )(eS1D − b) − eh S1D P D J44 , a0 = (m 12 m 21 − J44

so that the Routh–Hurwitz conditions become for it a0 > 0,

a2 > 0,

a2 a1 > a0 .

(14)

In principle, Hopf bifurcations could arise, provided the remaining stability conditions are still satisfied, if one of the following two cases holds D D J22 + J55 = 0,

a2 a1 = a0 .

(15)

22

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

Fig. 2. Damped oscillations in the two uncoupled patches 1 and 2, left and center, and in the system (1), right, are obtained for the parameter values: r1 = 2.5, r2 = 1, K 1 = 100, a = 0.5, c = 0.0, K 2 = 150, m 21 = 3.0, m 12 = 2.0, n 21 = 0, n 12 = 0, λ = 0.5, µ1 = 0.0, µ2 = 0.3, ν = 0.3, g = 0.0, e = 1.3, b = 0.1.

4.3. The predator-free state and the coexistence equilibrium The predator-free state E P = (S1P , I1P , 0, S2P , I2P ) and the ecosystem coexistence E ∗ = (S1∗ , I1∗ , P ∗ , S2∗ , I2∗ ) represent the most complicated situations, in which most, or all, the populations thrive. To assess analytically these equilibria appears to be a rather nontrivial task, as they amount to solve nonlinear systems with four and five unknowns respectively. In the former case, one eigenvalue can however be easily determined, eS1 + g I1 − b, so that if positive, the predator-free equilibrium is certainly not stable. We will investigate these situations numerically. Note that clearly, the choice of the parameters may influence the simulation results. But here, via these numerical experiments, we want essentially to show the different possible outcomes of the system, rather than obtaining precise quantitative information on the model behavior in a specific situation. 5. Simulations We have run experiments to assess the system’s outcomes using demographic values from the literature for the case of the rabbits and foxes mentioned in the introductory part of the paper, [2]. The remaining parameters have been given hypothetical values. 5.1. Comparison with the HTII model In this section we take n i j = 0, for i, j = 1, 2, i ̸= j as well as I1 (0) = 0. In this way, the model (1) is equivalent to the first one presented in [31], except for the type of functional response chosen. In Fig. 2 we show the behavior of system (1) for the same parameter values given in Fig. 3 of [31] that give persistent oscillations in the first patch as well as in the fully coupled model. Evidently, due to the choice of the different functional response, in this case we get damped oscillations that quickly decay to an equilibrium, in the two uncoupled patches and in the full system. This demonstrates that the type of functional response has a role in determining the system’s outcome. 5.2. Comparison with the HTI model In this section we take n i j ̸= 0, for i, j = 1, 2, i ̸= j and I2 (0) = 0. In this way, we can compare the model (1) with the last one presented in [31], and investigate what happens when the infected can migrate. In Fig. 3 we show the behavior of system (1) in the above conditions, when the infected migration rates do not vanish. We take the same parameters used in Figure 8 of [31], and allow infected migration between the patches. The parameters are r1 = 100, K 1 = 100, a = 4.5, c = 0, g = 0, e = 0.8, b = 0.2, r2 = 20, K 2 = 150, m 21 = 3, m 12 = 5, λ = 0.25, ν = 0.1, µ1 = 0.2, µ2 = 0.0. Evidently, the infected migration rates do not influence the susceptible prey equilibrium in patch 1. There is a quick transient for which for low values of n 21 the remaining populations in both patches are quickly reduced, except for the infected in patch 2 that attain instead a plateau. In particular, in such situation the disease is eradicated in patch 1, as well as the healthy prey in patch 2. The predators instead are reduced, but still thrive. The other migration rate n 12 seems to play a less important role.

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

23

Fig. 3. Behavior of the system (1) when infected are allowed to migrate back and forth between the two patches. This is obtained with the parameters of Fig. 8 of [31]. First row: patch 1, left to right S1 , I1 , P; second row: patch 2, left to right S2 , I2 .

Fig. 4. Behavior of the system (1) when infected are allowed to migrate back and forth between the two patches. This is obtained with the parameters of Fig. 9 of [31]. First row: patch 1, left to right S1 , I1 , P; second row: patch 2, left to right S2 , I2 .

Similarly in Fig. 4 we show again the behavior of system (1), when the infected migration rates do not vanish, taking the same parameters used in Figure 9 of [31], and allow infected migration between the patches. The parameters are r1 = 100, K 1 = 10, a = 4.5, c = 0, g = 0, e = 0.2, b = 0.2, r2 = 20, K 2 = 150, m 21 = 3, m 12 = 5, λ = 0.25, ν = 0.7, µ1 = 0.2, µ2 = 0.0. Here again the susceptible prey in patch 1 are unaffected by infected migration. As in the former case there is a sharp drop in all population values for increasing n 21 , but for the infected in patch 2, that show an increase but with respect to n 12 instead decrease at first to attain a plateau. For Figure 10 of [31], with parameters r1 = 100, K 1 = 10, a = 4.5, c = 0, g = 0, e = 0.2, b = 0.2, r2 = 20, K 2 = 150, m 21 = 100, m 12 = 0.01, λ = 0.95, ν = 0.7, µ1 = 0.2, µ2 = 0.0, the results of allowing infected

24

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

Fig. 5. Behavior of the system (1) when infected are allowed to migrate back and forth between the two patches. This is obtained with the parameters of Fig. 10 of [31]. First row: patch 1, left to right S1 , I1 , P; second row: patch 2, left to right S2 , I2 .

migration are shown in Fig. 5. Note that the predators are essentially extinguished, the scale on their vertical axis is 10−6 . Thus here healthy prey in patch 1 attain a plateau, while in patch 2 they quickly vanish even for low values of n 21 . The disease is eradicated in patch 1, infected in patch 2 show an increase with respect to both migration parameters. Another situation is depicted in Fig. 6, for Figure 11 of [31], for the case of Ovis canadensis or Strix occidentalis and the parameter values r1 = 2.5, K 1 = 100, a = 1.5, r2 = 1, K 2 = 150, γ = 0.9, ν = 0.1, µ = 0.2, e = 0.8, b = 0.9 and susceptible migration rates m 21 = 0.03, m 12 = 0.06. Predators, healthy prey in patch 2 and infected in patch 1 are extinguished, infected in patch 2 attain a plateau while healthy prey in patch 1 are at a constant level. In such situation thus the disease could be eradicated, at least in one of the patches, by biological controls, at the expense of losing also the predators’ population. If the latter is a pest, this would also be a way for successfully fighting it. 5.3. The sound-prey-only equilibrium Fig. 7 shows that the sound-prey-only equilibrium exists in both the coupled as well as the uncoupled model, for the parameter values r1 = 100, r2 = 20, k = 2, h = 5, m 21 = 13.6, m 12 = 43.13, n 21 = 0, n 12 = 0, σ = 0.1, θ1 = 0, θ2 = 0, g = 0.15, e = 4.5, b = 5.3. Note that in this case M = 32.4324 > 1 and M (1) + M (2) = 0.0890 < 1, so that the origin is unstable. 5.4. The disease-free equilibrium Fig. 8 shows that for e > b we obtain E 2 in place of E 1 and Q 1 as before in the uncoupled models. In the metaecoepidemic model we get E D , with the same S1 and a larger value of P at the expense of a smaller one for S2 . The parameters are: r1 = 100, r2 = 20, k = 2, g = 0.15, e = 4.5, h = 5, b = 4.3, m 21 = 13.6, m 12 = 43.13, n 21 = 0.0, n 12 = 0.0, σ = 0.3, θ1 = 0.0, θ2 = 0.0. The initial conditions are S1 (0) = 2.0, I1 (0) = 1, P(0) = 1, S2 (0) = 1.5, I2 (0) = 1. 5.5. The predator-free equilibrium Equilibrium E P can be attained by the parameter choice r1 = 100, σ = 5.9, h = 5, θ1 = 3.9, e = 4.5, g = 0.15, b = 14.3, r2 = 20, k = 0.2, θ2 = 0.5, m 21 = 3.6, m 12 = 3.13, n 21 = 0.0, n 12 = 0.0, as it is seen in Fig. 9. Note that the disease is wiped out in the second patch when isolated, but it remains endemic in it when migration is allowed, although only the healthy prey migrate.

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

25

Fig. 6. Behavior of the system (1) when infected are allowed to migrate back and forth between the two patches. This is obtained with the parameters of Fig. 11 of [31]. First row: patch 1, left to right S1 , I1 , P; second row: patch 2, left to right S2 , I2 .

Fig. 7. The sound-prey-only equilibrium E C exists in both the coupled as well as the uncoupled model. Left: patch 1 in the uncoupled model, for the parameter values r1 = 100, σ = 0.1, h = 5, θ1 = 0.0, e = 4.5, g = 0.15, b = 5.3. Center: patch 2 in the uncoupled model, for the parameter values r2 = 20, k = 2.0, θ2 = 0.0 and σ = 0.1 as for patch 1. Right: the fully coupled model, with the same parameter values as in the two separate patches, with migration rates given by m 21 = 13.6, m 12 = 43.13, n 21 = 0.0, n 12 = 0.0.

Fig. 8. The disease-free state E D exists in both the coupled as well as the uncoupled model. Left: patch 1 in the uncoupled model, for the parameter values r1 = 100, σ = 0.3, h = 5, θ1 = 0.0, e = 4.5, g = 0.15, b = 4.3. Center: patch 2 in the uncoupled model, for the parameter values r2 = 20, k = 2.0, θ2 = 0.0 and σ = 0.3 as for patch 1. Right: the fully coupled model, with the same parameter values as in the two separate patches, with migration rates given by m 21 = 13.6, m 12 = 43.13, n 21 = 0.0, n 12 = 0.0.

26

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

Fig. 9. The predator-free state E P exists in both the coupled as well as the uncoupled model. Left: patch 1 in the uncoupled model, for the parameter values r1 = 100, σ = 5.9, h = 5, θ1 = 3.9, e = 4.5, g = 0.15, b = 14.3. Center: patch 2 in the uncoupled model, for the parameter values r2 = 20, k = 0.2, θ2 = 0.5 and σ = 5.9 as for patch 1. Right: the fully coupled model, with the same parameter values as in the two separate patches, with migration rates given by m 21 = 3.6, m 12 = 3.13, n 21 = 0.0, n 12 = 0.0.

Fig. 10. The whole ecosystem thrives, both in the two uncoupled patches as well as in the fully coupled model. Left: patch 1 in the uncoupled model, for the parameter values r1 = 100, σ = 1.5, h = 5, θ1 = 1.0, e = 4.5, g = 0.19, b = 5.3. Center: patch 2 in the uncoupled model, for the parameter values r2 = 20, k = 0.8, θ2 = 0.1 and σ = 1.5 as for patch 1. Right: the fully coupled model, with the same parameter values as in the two separate patches, with migration rates given by m 21 = 1.0, m 12 = 0.1, n 21 = 1.04, n 12 = 1.75.

5.6. Ecosystem’s coexistence The ecosystem coexistence E ∗ = (S1∗ , I1∗ , P ∗ , S2∗ , I2∗ ) is attained, see Fig. 10. Here the parameter values are r1 = 100, h = 5, n 12 = 1.75,

r2 = 20, b = 5.3,

k = 0.8, m 21 = 1.0,

σ = 1.5,

θ1 = 0.1,

g = 0.19, m 12 = 0.1,

e = 4.5,

(16)

n 21 = 1.04,

θ2 = 0.1.

However, note that in other simulations, not reported here, the coexistence occurs in the coupled system although the first patch, when disconnected from the other one, becomes disease-free. 6. Effects of parameter changes In this section we analyze the bifurcations diagrams in the parameter space. Starting from the system’s coexistence equilibrium, with parameter values given by (16), we take pairs of parameters and plot the change in the equilibrium values in this two-dimensional parameter space. In the remainder of this section we illustrate the changes that occur in the system dynamics based on a suitable selection of these pairs of parameters, focusing on migration rates alone, migration rates coupled with demographic or epidemic parameters, combinations of demographic and epidemic parameters and finally of demographic parameters alone and epidemic parameters alone. 6.1. Changes in the migration rates In particular, in Figs. 11–13 we consider just the effects of the changes in some of the migration rates. The situation corresponding to the uncoupled patches would correspond in these figures to the values of the surfaces at the origin. In

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

27

Fig. 11. Behavior in the n 21 –n 12 parameter space. The coupled system, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 .

these cases we do not plot the situation in the uncoupled model, as the population equilibrium values are independent of these parameters, and therefore the corresponding surfaces are horizontal planes. Fig. 11 shows the system behavior corresponding to introducing also the infected individuals migration. Susceptibles in patch 2 are seen to decrease for larger values of n 21 and vice versa to grow with increasing values of n 12 . Thus the emigration rate of infected from patch 1, n 21 , has the effect of contributing to decrease the susceptibles in the arrival patch, while the emigration rate of infected from patch 2, n 12 contributes to have a larger population of susceptibles in that patch. These quite counterintuitive results can however be explained by observing that there will be less infected from which to contract the disease, and therefore the susceptibles increase their numbers. Another situation occurs instead for the susceptibles in patch 1, decreasing sharply for a higher immigration rate of infected from patch 2, while their equilibrium population is less affected by the emigration rate of infected from patch 1, n 21 . The same above motivation could be applied in this situation as well, observing that contacts with infected in this case become more frequent. Further, in patch 1, the infected prey and predators populations attain a plateau almost independently of these migration rates. As expected, instead, the infected in patch 2 decrease sharply when their emigration rate becomes large and increase with a larger immigration rate. We consider next the changes due to migration of susceptibles, Fig. 12. Here infected in both patches, predators and susceptibles in patch 2 behave all similarly, namely increase with higher emigration rates from patch 1 and decrease for larger values of the corresponding immigration rate into patch 1. The result appears to be quite counterintuitive for the populations in patch 2. Instead, susceptibles in patch 1 as expected decrease whenever the emigration rate becomes larger, m 21 , and increase when the immigration rate raises. In Fig. 13 the effects of combined migration rates of both susceptibles and infected are investigated, specifically n 12 and m 21 . Again, in patch 1 predators and infected prey are scarcely affected by changes in these rates: after a sharp raise for low values of n 12 , their populations quickly attain a plateau which remains essentially constant. Susceptibles in patch 2 rise with increasing values of both migration rates, whether be it immigration of susceptibles or infected emigration. With larger m 21 this is expected, but it is remarkable that a higher emigration rate of infected contributes to raise the susceptibles population. Susceptibles in patch 1 are depleted by their emigration rate m 21 , and essentially unaffected by the immigration rate of infected, apart from a quick sharp decline that occurs for low values of this rate. The infected in patch 2 are almost unaffected by the susceptible migration rate, while they become extinguished when their emigration rate becomes larger and larger, as it should be expected. This shows that a possibility of fighting

28

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

Fig. 12. Behavior in the m 12 –m 21 parameter space. The coupled system, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 .

10

2.5

8

2

6

1.5

P1

1.15 1.1 1.05 1 0.95 0.9 0.85 0.8 0.75 10

I1

S1

Coupled system: Patch 1

4 2 0 10

8 6

4 m21

2

0 0

2

0 10 8

8

6

4

1 0.5

6

n12

8 4

2

2

0 0

m21

8

8

6

6

4

4

2 0 0

m21

n12

2

4

6

n12

1.1 1.05 1 0.95 0.9 0.85 0.8 0.75 0.7 10

120 100 80

I2

S2

Patch 2

60 40 20 0 10

8 6

4 m21

2

0 0

2

4 n12

6

8

8 6

4 m21

2

0 0

2

4

6

8

n12

Fig. 13. Behavior in the n 12 –m 21 parameter space. The coupled system, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 .

a disease in one patch by means of biological control for instance would be to allow larger emigration rates of the infected. But this is kind of compensated by a higher endemic value in the other patch.

29

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44 Coupled system: Patch 1

10

0.95

2

0.85

P1

8

I1

S1

2.5

9

0.9

7

1.5

0.8

6

1

0.75 10

5 10

0.5 10

8 6

4 m21

2

0 0

2

8

8

6

4

6

n21

8 4

2

2

0 0

m21

8

8

6

6

4

4

2 0 0

m21

n21

2

4

6

n21

1.1 1.05 1 0.95 0.9 0.85 0.8 0.75 0.7 10

50 40

I2

S2

Patch 2

30 20 10 0 10

8 6

4 m21

2

0 0

2

4

6

8

8

n21

6

4 m21

2

0 0

2

4

6

8

n21

Fig. 14. Behavior in the n 21 –m 21 parameter space. The coupled system, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 .

Studying the combined effect of the migration rates of both infected and susceptibles with a different pair of parameters, namely the emigration rates from patch 1 for both the susceptibles and infected, we discover that susceptibles in patch 1 decrease as expected with their larger emigration rates and are almost unaffected by the infected emigration, Fig. 14. The latter do not seem to affect too much also the infected in patch 1, while higher emigration rates of susceptibles from patch 1 lead to a higher number of infected in patch 1. Again this counterintuitive remark can be explained by observing that in this case there are more contacts with infected per susceptible. The predators population when either migration rate is small will become smaller, but when both migration rates are sizeable, seems to increase with both parameters increasing. In patch 2 susceptibles decline with a higher infected immigration rate n 21 , and as expected increase with a higher susceptible immigration m 21 . Infected instead seem to be less sensitive to the susceptible immigration rate, but increase with a corresponding infected immigration rate, as it can be expected. 6.2. Combinations of migration rates with demographic parameters In Fig. 15 we consider the combination of predators’ mortality and healthy prey emigration rate from patch 2. Comparing the uncoupled patches and the connected system, we see that the healthy and infected prey in patch 1 behave in the very same way. In patch 2 both types of prey are at a constant level in the disconnected system, while instead when the patches communicate the healthy prey decrease for a larger emigration rate m 12 , which is to be expected. The infected however behave like in patch 1, although at lower equilibrium population values. The behavior is independent of the susceptibles migration rate. A transcritical bifurcation occurs at a value of the predators’ mortality rate b∗ ≈ 4, in both the uncoupled and connected systems. In the former case the disease-free equilibrium is replaced by the system coexistence when this mortality exceeds b∗ . Its higher values lead eventually to predators’ extinction. The sharp decline in the predators population is more marked in the coupled system, and therefore the coexistence of all the populations in this case is prevented, as the predators become extinct even for values of their mortality below b∗ . Fig. 16 contains the analysis in the n 12 –r1 parameter space. Since the demographic parameter concerns only patch 1, as in the previous case when the system is disconnected, the populations in patch 2 attain a constant value. In patch 1 instead we find obviously that the populations behavior is independent of the migration rate. There appears to be a transcritical bifurcation for which the predator-free equilibrium is replaced by coexistence of the three populations when

30

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

1

20

20

0.8

15

15

P1

0.6

I1

S1

Uncoupled Patches: Patch 1

10

0.4 5

0.2 0 10

5

0 10 8 6

4

2

m12

2

0 0

0 10

8

8

6

4

10

6

b

8 4

2

4

2

0 0

m12

8

8

6

6

4

2 0 0

m12

b

2

4

6

b

Patch 2

1

15

0.8 10

I2

S2

0.6 0.4

5

0.2 0 10 5 m12

4

2

0 0

0 10

8

6

5 m12

b

0 0

8

6

4

2

b

1

20

0.8

15

P1

0.6

I1

S1

Coupled system: Patch 1

10

0.4 5

0.2 0 10

0 10 8 6

4 m12

2

4

2

0 0

8

8

6

6

b

400 350 300 250 200 150 100 50 0 10

8 4

2 0 0

m12

8

8

6

6

4

2

4

2

m12

b

0 0

2

4

6

b

0.9

12

0.8

10 8

0.7

I2

S2

Patch 2

0.6

0.4 10

6 4

0.5

2 0 10

8 6

4 m12

2

0 0

2

4 b

6

8

8 6

4 m12

2

0 0

2

4

6

8

b

Fig. 15. Behavior in the b–m 12 parameter space. Top frame: the uncoupled patches, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 . Bottom frame: the coupled system, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 .

the prey net reproduction rate r1 exceeds the value r1∗ ≈ 5.2. When the patches are connected, the susceptible prey in patch 1 quickly reach a plateau; they exhibit more or less a similar behavior also in patch 2. Also the infected in patch 1 attain a plateau, but this occurs after a longer transition, in terms of the parameter r1 ; note that for low values of the latter, the disease is eradicated. In patch 2 instead and as intuitively expected there is a sharp decline of the endemicity

31

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44 Uncoupled Patches: Patch 1

0.78

0.2

15

0.15

P1

10

0.76

I1

S1

0.77

0.75

5

0.74 0 10

0.73 10 5 r1

2

0 0

6

4

0.1 0.05 0 10

8 5

n12

2

0 0

r1

6

4

8 5

n12

0 0

r1

2

4

6

8

n12

Patch 2

15

1

10

0.6

I2

S2

0.8

0.4

5

0.2 0 10

0 10 5 2

0 0

r1

4

6

8 5

n12

0 0

r1

2

6

4

8

n12

Coupled System: Patch 1

1

15

0.15 0.1

0.6

P1

10

I1

S1

0.8

0.4

5

0

0.2 0 10 5 r1

2

0 0

6

4

0 10

8

0.05

–0.05 10 5

n12

2

0 0

r1

6

4

8 5

n12

0 0

r1

2

4

6

8

n12

Patch 2

30

0.9

20

I2

S2

1

10

0.8 0.7 10 5 r1

0 0

2

4 n12

6

8

0 10 5 r1

0 0

2

4

6

8

n12

Fig. 16. Behavior in the n 12 –r1 parameter space. Top frame: the uncoupled patches, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 . Bottom frame: the coupled system, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 .

of the disease for higher emigration rates, n 12 . These combined with low values of r1 ensure also in patch 2 the disease eradication, but this leads also to extinction of predators. For them to survive, indeed, higher values of the prey net reproduction rate r1 are needed. This is kind of unexpected, but it can be explained by the fact that the predators are assumed to be specialists, i.e. they feed only on the prey modeled in the system. Note that very low values of the infected immigration rate into patch 1 n 12 lead to predators’ extinction even in the presence of relatively high values of r1 .

32

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

10

2

8

1.5

0.6

6

0.4

P1

1 0.8

I1

S1

Uncoupled Patches: Patch 1

4

0.2

1 0.5

2

0 10

0 10

0 10 5

4

2

0 0

m12

8

6

5

r2

0 0

m12

4

2

8

6

5

r2

0 0

m12

6

4

2

8

r2

Patch 2

5 1

4 3

0.6

I2

S2

0.8

2

0.4 1

0.2

0 10

0 10 5 0 0

m12

4

2

8

8

6

8

6

6

4

4

2

r2

2 0 0

m12

r2

7

0.96

6.5

0.94

6

0.92

1.6 1.4 1.2

P1

0.98

I1

S1

Coupled System: Patch 1

5.5

0.6

0.9

5

0.88 10

4.5 10 8

0.4

m12

2

2

0 0

6

4

4

2

m12

r2

8

6

6

4

4

2

0 0

8

8

6

6

4

0.2 10 8

8

6

1 0.8

r2

4

2

m12

0 0

2

r2

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 10

4.5 4 3.5

I2

S2

Patch 2

3 2.5 2 1.5 10

8

8

6

6

4 m12

4

2 0 0

2

r2

8

8

6

6

4 m12

4

2 0 0

2

r2

Fig. 17. Behavior in the r2 –m 12 parameter space. Top frame: the uncoupled patches, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 . Bottom frame: the coupled system, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 .

Combining the healthy prey net reproduction rate in patch 2, r2 , with their emigration rates m 12 , Fig. 17, shows that patch 1 is totally unaffected in the disconnected case, patch 2 has a constant healthy prey population, while the

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

33

infected population grows linearly with r2 and is clearly not affected by the migration rate because the latter for the uncoupled system vanishes. When the system is coupled, we find that in the range explored there is always coexistence of all the populations. The healthy prey in patch 1 increase with increasing immigration rate and are independent of r2 , their reproduction rate in the other patch. All other populations in the system are also scarcely affected by r2 but show the opposite behavior with respect to the migration rate m 12 , always decreasing for its larger values. For the healthy prey in patch 2 this remark is natural, for the remaining populations in the system this is an interesting observation. 6.3. Combinations of migration rates with epidemic parameters We provide two such examples. In the first one, Fig. 18, the disease recovery rate in patch 2 is matched with the infected emigration rate from patch 1. When the system is disconnected, in patch 1 evidently the populations are unaffected by these parameter changes. In patch 2 again the populations are clearly insensitive to n 21 , but they show a narrow strip of coexistence for very low values of the disease recovery rate. When the latter exceeds θ2∗ ≈ 0.2 the disease is eradicated. When the patches are connected, this behavior is repeated in patch 2, with only a change due to a linear growth of the infected with respect to the immigration rate n 21 as it should be expected. In patch 1 instead, the healthy prey attain again a plateau, but both the infected and predators show a trough; only small values of either one of the two parameters in consideration here allow these populations survival. Hence in this case disease eradication could be achieved by tuning both these parameters to moderately high values, and coupling the patches. This is obtained of course at the price of sacrificing the predators as well. Considering instead the contact rate σ together with the infected emigration rate from patch 2, n 12 , Fig. 19, again shows that in the uncoupled system the latter rate has no effect on the populations. In patch 1 susceptibles decline with higher contact rate, infected instead grow with it, predators initially grow to attain a peak but then decline. In patch 2 this behavior is mimicked by the infected prey, while susceptible prey behave like in patch 1. Connecting the patches does not change the populations’ equilibria in patch 1, nor the healthy prey behavior in patch 2. Only, and importantly, the infected population in patch 2 is depleted, surviving only for moderately large values of the disease transmission rate coupled with low values of their emigration rate. In this case a possible tool for the disease eradication, but only in one of the patches, would be available for the system administrator. 6.4. Combinations of demographic and epidemic parameters We have then considered changes involving pairs of parameters, one of demographic nature and the other one epidemiological. At first we consider the disease recovery rate in patch 1 θ1 combined with two demographic parameters, the healthy prey net reproduction rate in patch 2, r2 , Fig. 20, and the predators’ mortality b, Fig. 21. In the former case patch 1 in the uncoupled system is independent of r2 and from coexistence quickly attains a healthy prey-only equilibrium for moderate values of the disease recovery rate in the same patch, θ1 through a transcritical bifurcation occurring at θ1∗ ≈ 0.8. This behavior does not change significantly if the patches are connected. In patch 2, when disconnected from the other one, susceptible prey attain a constant level, this is a bit unexpected as they seem to be insensitive to their own net reproduction rate r2 . The main changes involve the infected in patch 2: when the system is disconnected, they are obviously not affected by the disease recovery rate in patch 1, but they grow linearly with the susceptibles net reproduction rate r2 , which is kind of interesting. If the migration is allowed, however, they become quickly extinguished for moderate values of θ1 , which is also a remarkable situation, and do not seem instead to be affected by their own susceptibles net reproduction rate in the very same environment, r2 . Thus in this case the disease could be fought by suitably tuning these parameters. When studied together with the predators’ mortality b, the disease recovery rate in patch 1, θ1 shows the coexistence, the disease-free and the healthy prey-only equilibria in patch 1 if not connected with the other one. More precisely, the disease-free equilibrium is attained for values of b below b∗ ≈ 4, independently of the disease recovery rate; this is a remarkable observation. Coexistence occurs in a narrow strip for large b and small θ1 , the healthy prey only equilibrium is attained in the remaining domain, i.e. for b > b∗ and moderate to large values of θ1 . In patch 2 the populations are not affected by these parameters and thus achieve a constant value. When the system is connected, the behavior in patch 1 hardly changes, in patch 2 disease endemicity is found only for low θ1 and high predators’ mortality rate b, which is a fact to be noted. Again this shows a possible tool for disease eradication.

34

10

2

8

1.5

0.6

6

0.4

4

0.2

2

0 10 5 0

0

6

4

2

1 0.5

0 10

8

n21

P1

1 0.8

I1

S1

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

5

θ2

n21

0

0

6

4

2

0 10

8

8 5

θ2

n21

0

0

6

4

2

θ2

250

0.8 0.78

200 150

0.74 0.72 0.7

I2

S2

0.76 100 50

0.68 0.66 10

0 10

8 4

n21

2

2

0.9

8 2 0

θ2

8

8

6

n21

4 2

4 0

2

0

40

0.85

30

0.8

20

0.75

10 0 10 8

8 2

6

4 0

0

θ2

1

0.5

50

n21

2

0

2

0.9

4

4

1.5

θ2

6

2 0

0.95

0.7 10

6

4

2.5

6

4

2

0

S2

n21

4

n21

I2

6

8 7 6 5 4 3 2 1 0 10

8 6

6

P1

I1

S1

0.95

8

4

0 0

1

0.85 10

8

8

6

2

θ2

0 10

8

6

n21

θ2

8

6 4 2

8

6

4

θ2

8

6

n21

0

0

2

4

2

2 0

0

4

6

θ2

Fig. 18. Behavior in the θ2 –n 21 parameter space. Top frame: the uncoupled patches, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 . Bottom frame: the coupled system, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 .

Then we consider the disease transmission rate σ in conjunction with the healthy prey carrying capacity in patch 2, k, and with the healthy prey net reproduction rate in patch 1, r1 . In the former case, Fig. 22, the behaviors of the populations are the same whether the patches are connected or not, with the exception of the predators. The latter when the system is uncoupled are unaffected by k while their equilibrium raises up and then drops with increasing σ . Instead if the system is coupled k influences linearly their behavior, leading to higher values when it gets larger. Further, in patch 1 the prey are insensitive to k too, the healthy

35

25

1

20

0.8

15

P1

1.2

I1

S1

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

0.6

10

0.4

5

0.2 10

0 10

8

8

6 4

2

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 10

4 2

2 0

4

4

2

6 4

4

2

2

0 0

5 4

15

3

10

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 10

6

σ

6

σ

2 1 0 10

8

8

6

6

n12

σ

4

4

2

2

0 0

6

8

8

6

n12

σ

4

2

4 0

0

2

6

σ

250 200

I2

0 0

S2

n12

2

20

8

4

8 4

5

8

2

0 0

6

n12

σ

25

0 10

2

n12

σ

P1

0

8

6

6

8

6

4

2

I1

S1

n12

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 10

40 35 30 25 20 15 10 5 0 10

8

6

4

2

0 0

8 4

8

8

6

n12

σ

0 0

S2

2

4

I2

n12

8

6

4.5 4 3.5 3 2.5 2 1.5 1 0.5 0 10

150 100 50 0 10

8

8

6

n12

4

4

2 0 0

2

6

σ

8 6

n12

8 4

4

2 0 0

2

6

σ

Fig. 19. Behavior in the σ –n 12 parameter space. Top frame: the uncoupled patches, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 . Bottom frame: the coupled system, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 .

ones drop with increasing σ , while instead the infected grow. In patch 2 the healthy prey increase linearly with k, as it should be expected, but drop significantly to attain a low plateau for moderate to large values of the disease contact rate. A similar behavior occurs for the infected, but they attain a higher equilibrium value than the healthy ones. This infected equilibrium is also larger if the patches become connected. In Fig. 23, we find the system behavior as a function of the healthy prey net reproduction rate in patch 1, r1 , together with the disease contact rate. Susceptible prey behave in both patches, whether connected or not, in the same way,

36

7 6 5 4

P1

1.02 1 0.98 0.96 0.94 0.92 0.9 0.88 10

I1

S1

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

3 2 1 0 10

8

r2

8

8

6 4

2

0

0

r2 4

θ1

8

8

6

6

4 2

1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 10

2 0

4

2

0

8

6

6

4

r2

θ1

2 0

0

6

4

2

θ1

5

0.8

4

0.6

3

I2

S2

6 1

2

0.4

1

0.2 0 10 5 2

0

1 0.98 0.96 0.94 0.92 0.9 0.88 0.86 10 6

r2

2

2

0

S2

0

4

8

r2

θ1

8

6

6

r2

θ1

2 0

10 8

r2

8 4 2

0

0

4 2 0

0

2

6

4

θ1

4

2

0

6

4

2

θ1

1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 10 8

8

8 4

2 1.8 1.6 1.4 1.2 1 0.8

6

8

6

8 7 6 5 4 3 2 1 0 10

8 4

8

6

I2

8

4

P1

0

I1

S1

r2

0 10

6

6

r2

θ1

4 2

0 0

2

4

6

θ1

4.5 4 3.5 3 2.5 2 1.5 1 0.5 0 10 8

6

r2

8 4 2

4 0

0

2

6

θ1

Fig. 20. Behavior in the θ1 –r2 parameter space. Top frame: the uncoupled patches, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 . Bottom frame: the coupled system, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 .

namely r1 does not clearly influence them except slightly for patch 1 when migration is allowed. For larger values of the disease transmission rate σ their equilibrium value goes down, as expected. Infected disappear for low values of the contact rate in every situation. When σ grows, their equilibrium raises up and then declines slightly, while it is insensitive to r1 in patch 2 when the system is disconnected, and in both patches when the migration is allowed. Infected in patch 1 instead grow linearly with σ only if the two patches are segregated. Predators are extinguished except for a small wedge corresponding to moderate values of σ and high values of r1 when migration does not occur.

37

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

8

6

b

4

2

0 0

2

15 10 5

8

8

6

4

20

40 35 30 25 20 15 10 5 0 10

P1

1.4 1.2 1 0.8 0.6 0.4 0.2 0 10

I1

S1

Uncoupled Patches: Patch 1

6

b

θ1

4

2

0 0

6

4

2

0 10

8

8

6

4

2

b

θ1

0 0

8

6

4

2

θ1

Patch 2

1

15 10

0.6

I2

S2

0.8 0.4

5

0.2 0 10 5

b

4

2

0 0

0 10

8

6

5

θ1

2

0 0

b

8

6

4

θ1

1

0.6

I1

S1

0.8

0.4 0.2 0 10

8

6

4

b

2

2

0 0

8

6

4

40 35 30 25 20 15 10 5 0 10

P1

Coupled system: Patch 1

8

6

4

b

θ1

2

0 0

6

4

2

8

30 25 20 15 10 5 0 10

8

6

b

θ1

4

2

0 0

2

4

6

8

θ1

0.86 0.84 0.82 0.8 0.78 0.76 0.74 0.72 10

20 15

I2

S2

Patch 2

10 5

8

0 10 6

b

4

2

0 0

2

4

θ1

6

8

8

6

b

4

2

0 0

2

4

6

8

θ1

Fig. 21. Behavior in the θ1 –b parameter space. Top frame: the uncoupled patches, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 . Bottom frame: the coupled system, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 .

They have a similar behavior if the system is connected, but in this case they survive also for smaller values of r1 . Thus in this case there is a way of controlling the predators population, in case it is considered a pest. 6.5. A pair of epidemic parameters Here we observe the system’s behavior as a function of the disease transmission rate and its recovery rate in patch 1, θ1 , Fig. 24. The situation in patch 1 is the same for both the connected and coupled cases. Along a straight line in the parameter space there is a transcritical bifurcation for which for high values of the disease recovery rate and even

38

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

25

1

20

0.8

15

0.6

10

0.4

5

0.2 10

P1

1.2

I1

S1

Uncoupled Patches: Patch 1

0 10 8

6

4

k

2

2

0 0

6

4

8

8

6

4

k

σ

2

0 0

6

4

2

4.5 4 3.5 3 2.5 2 1.5 1 0.5 0 10

8

8

8

6

k

σ

4

2

0 0

6

4

2

σ

Patch 2

10 8

I2

S2

6 4 2 0 10

8

6

4

k

2

0 0

8

8

6

4

2

400 350 300 250 200 150 100 50 0 10

6

4

k

σ

2

0 0

8

6

4

2

σ

25

1

20

0.8

15

0.6

10

0.4

5

0.2 10

P1

1.2

I1

S1

Coupled system: Patch 1

8

0 10 6

4

k

2

4

2

0 0

6

8

8

6

4

k

σ

2

2

0 0

6

4

8

12 10 8 6 4 2 0 10

8

6

k

σ

4

2

0 0

2

4

6

8

σ

Patch 2

10

6

I2

S2

8

4 2 0 10

8

6

k

4

2

0 0

2

4

σ

6

8

35 30 25 20 15 10 5 0 10

8

6

4

k

2

0 0

2

4

6

8

σ

Fig. 22. Behavior in the σ –k parameter space. Top frame: the uncoupled patches, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 . Bottom frame: the coupled system, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 .

low values of the contact rate, the disease is eradicated, healthy prey attaining a plateau and predators are extinguished as well. For high σ and low θ1 , the three populations thrive together. In patch 2, when the system is uncoupled, the populations are clearly insensitive to θ1 , while with increasing σ the infected first raise up to a peak and then slowly decrease, while the susceptibles expectedly show a corresponding decrease. Note that for low values of σ the disease is eradicated. Connecting the system only slightly affects the healthy prey, but sensibly reduces the infected equilibrium

39

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

25 20 15

P1

1.4 1.2 1 0.8 0.6 0.4 0.2 0 10

I1

S1

Uncoupled Patches: Patch 1

10 5 0 10 8

6

4

r1

2

2

0 0

6

4

8

8

6

r1

σ

4

2

2

0 0

6

4

σ

8

0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 10

8

6

4

2

r1

0 0

2

8

6

4

σ

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 10

I2

S2

Patch 2

8

6

r1

4

2

2

0 0

8

8

6

4

40 35 30 25 20 15 10 5 0 10

6

r1

σ

4

2

0 0

8

6

4

2

σ

25 20

0.6

15

0.4

10

0.2

5

0 10

0 10 8

6

r1

4

2

4

2

0 0

6

8

0.6 0.5 0.4 0.3 0.2 0.1 0 –0.1 10

P1

1 0.8

I1

S1

Coupled system: Patch 1

8

6

r1

σ

4

2

0 0

6

4

2

8

8

6

4

2

r1

σ

0 0

2

4

6

8

σ

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 10

I2

S2

Patch 2

8

6

r1 4

2

0 0

2

4

σ

6

8

16 14 12 10 8 6 4 2 0 10

8

6

r1

4

2

0 0

2

4

6

8

σ

Fig. 23. Behavior in the σ –r1 parameter space. Top frame: the uncoupled patches, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 . Bottom frame: the coupled system, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 .

values, decreasing with higher values of the disease transmission and the area in the parameter space where the disease gets eradicated becomes larger. In this situation then there would be the possibility of controlling the disease, as well as the predators population.

40

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

6.6. Effects of one pair of demographic parameters Lastly, we observe the system behavior in the k −r1 parameter space, Fig. 25. In patch 2 the two populations behave in the same way, independently of the fact that the system is connected or not, namely are clearly insensitive to their net reproduction rate r1 in the other patch and instead grow with increasing k, a fact that is certainly expected. Patch 1 when migration is not allowed shows a behavior independent of k, again as we should expect, with a transcritical bifurcation for which for increasing values of r1 we find at first the predator-free equilibrium followed by coexistence of the three populations. When the patches are connected, the behavior unexpectedly becomes essentially insensitive to r1 , except for the susceptible prey quickly attaining a constant value, and increasing with larger values of k for all the populations. 7. Discussion The model in which a refuge is present for both sound and infected prey is prone to disappear in particularly unfavorable circumstances. We have identified threshold quantities M, M (i) , i = 1, 2 which guarantee the survival of the ecosystem. This occurs if the activity in at least one patch exceeds one, where activity means the ratio of the migration rate over the net reproduction rate. Alternatively, the combined net reproduction rates in both patches must be larger than the combined migration rates among patches for sound prey. The system allows several equilibria. The sound prey-only equilibrium is feasible for a few combinations of the activities in each patch, see (8). To ensure stability for this point, additional conditions must be satisfied, (9), among which the size of the prey in the first patch must be low enough, as well as the combined sizes in the two patches. Further restrictions must be imposed on the infected migration rates, and on the activities of sound prey. The simulations reveal that when cutting the communications between patches, the populations settle at different levels. The presence of the predators in the environment cannot be assured only by changes in the migration rates. The disease-free state differs from the former in view of the fact that it contains also the predators. Only one such point can be feasible, although they come in pairs, provided that the sound prey in the refuge exceed a lower bound. Stability is too complicated to study analytically, but the equilibrium has been shown graphically in Fig. 8. This is an important finding, as the disease gets eradicated and moreover this happens also in the separate patches. Thus under these circumstances, there is no big harm in fragmenting the environment. But with a more virulent disease our simulations, not reported here, show that the disease would become endemic in the fragmented habitat, invading also the second patch, which instead is an epidemic-free habitat when communications are cut, and further the predators would be wiped out. Thus the previous conclusion on the harmless operation of secluding environments must be weighted, since it depends on how large the disease contact rate is. Fig. 9 shows that the predator-free equilibrium can be attained. Comparing this model with the ones presented in [31], we have found that the choice of the functional forms has indeed a major influence on the system’s outcome, as it should be expected. Moreover, allowing the infected to migrate may radically change the system’s equilibria. In addition it allows more degrees of freedom so that it becomes possible to use these additional control parameters to drive the system either to disease eradication, if the problem is to fight it, or alternatively to predators’ extinction, if the latter are a pest to be eradicated. We have also examined the effects of changing parameters in pairs, so as to see the outcome on the system coexistence equilibrium. We have discovered that transcritical bifurcations may occur, for which the coexistence equilibrium becomes either the disease-free or the healthy-prey-only equilibrium, when one of the parameters crosses a critical threshold. For instance, predators extinction could be achieved by suitable choices of the disease contact rate and disease recovery in patch 1, Fig. 24, or the infected emigration rate from patch 2 combined with the prey reproduction rate in patch 1, Fig. 16 or finally the disease contact rate with the prey reproduction rate in patch 1, Fig. 23, in addition to some other possible situations. Instead, the infected emigration rate from patch 2 is seen to be critical for disease eradication in several situations. The disease removal in patch 2 can be obtained for example by combinations of this parameter with the corresponding infected immigration rate into patch 2, Fig. 11, the immigration rate of susceptibles into patch 2, Fig. 13, or disease contact rate, Fig. 19. The system becomes disease-free in both patches if n 12 is suitably combined with the prey reproduction rate in patch 1, Fig. 16, or if the infected immigration rate into patch 2, n 21 is suitably coupled with disease recovery rate in patch 2, Fig. 18, but a disease-free environment can be obtained also combining an

41

1.2

25

1

20

0.8

10

0.6 0.4

5

0.2 10

0 10

2

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 10

6

θ1

σ

8

6

4

2

4

2

40 35 30 25 20 15 10 5 0 10

8

6

4

θ1

σ

2

0 0

2

2

0 0

S2

θ1

4

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 10

8

6

4

2

2

0 0

2

8

4

σ

8

6

4

6

4

σ

25

5

20

4

15

3 2

5

1

0 10

0 10

8

6

θ1

σ

4

2

4

2

0 0

6

8

8

6

θ1

σ

4

2

0 0

2

4

6

8

σ

15 10

I2

6

6

θ1

σ

10

8

8

8

6

4

2

0 0

8

6

4

2

0 0

I1

θ1

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 10

8

8

6

4

2

0 0

S2

4

4.5 4 3.5 3 2.5 2 1.5 1 0.5 0 10

P1

6

I2

8

θ1

S1

15

P1

I1

S1

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

5

8

6

θ1

4

2

0 0

2

4

σ

6

8

0 10

8

6

θ1

4

2

0 0

2

4

6

8

σ

Fig. 24. Behavior in the σ –θ1 parameter space. Top frame: the uncoupled patches, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 . Bottom frame: the coupled system, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 .

epidemiological parameter like disease recovery θ1 and demographic parameters, r2 Fig. 20, or b, the predators’ mortality, Fig. 21. Using two epidemiological parameters leads possibly also to the same result, Fig. 24.

42

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44 Uncoupled Patches: Patch 1

12

0.78

10

0.77

I1

S1

0.76

P1

8 6

0.75

4

0.74

2

0.73 10

0 10

8

6

r1

4

2

0 0

8

8

6

4

2

8

6

4

r1

k

2

8

6

6

4

2

0 0

0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 10

r1

k

4

2

2

0 0

8

6

4

k

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 10

I2

S2

Patch 2

8

8

6

140 120 100 80 60 40 20 0 10

8

6

r1 4

2

2

0 0

8

6

r1 4

4

k

2

2

0 0

4

6

k

20 15

0.6 0.4

10 5

0.2 0 10

P1

1 0.8

I1

S1

Coupled System: Patch 1

8

6

r1

0 10

8 4

2

0 0

4

2

8

6

6

r1

k

4

2

2

0 0

6

4

8

4 3.5 3 2.5 2 1.5 1 0.5 0 10

8

6

r1 4

k

2

0 0

2

4

6

8

k

1.4 1.2 1 0.8 0.6 0.4 0.2 0 10

25 20

I2

S2

Patch 2

15 10 5

8

6

r1 4

2

0 0

2

4

k

6

8

0 10

8 6

4

r1

2

0 0

2

6

4

8

k

Fig. 25. Behavior in the k–r1 parameter space. Top frame: the uncoupled patches, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 . Bottom frame: the coupled system, first row patch 1, populations (left to right) S1 , I1 , P1 ; second row patch 2, populations (left to right) S2 , I2 .

These results indicate that fragmenting the environment in such circumstances could be a useful tool for eradicating the epidemic, but also the extinction of predators could arise as a consequence. This could be a warning against possible unwanted outcomes of natural environment modifications. But also the reverse possibility could arise, with a high prevalence disease in the refuge patch being eradicated and the ecosystem surviving, by decoupling the model.

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

43

The numerical simulations indicate that the whole ecosystem would survive, with the disease affecting also the first patch when communications between both environments are possible, while in the uncoupled model the disease would remain endemic only in the refuge patch. In this case there is certainly a benefit, in that the very high prevalence of the disease in the refuge patch is dramatically reduced when coupling the systems. In spite of the theoretical possibility of having Hopf bifurcations in some situations, sustained oscillations have not been discovered in our extensive simulations. Our simulations further show that one equilibrium is attained in every case for all the parameter choices, suggesting the global stability of each equilibrium. In view of all these different possibilities, this study strongly suggests the necessity in practical situations of carefully evaluating the environmental state and relationships among its constituents, before proceeding to any alteration of the landscape, since the outcomes can be very different from the ones one might naively expect. Acknowledgments The author is very much thankful to his biology colleagues Guido Badino, Francesca Bona, Marco Isaia and Simona Bonelli, who provided valuable material for the biological background of this research. The project was partly funded by the grant “Analisi Numerica per le Scienze della Vita” (Numerical analysis for the life sciences) of the Dipartimento di Matematica. References [1] V. Ajraldi, M. Pittavino, E. Venturino, Modelling herd behavior in population systems, Nonlinear Analysis RWA 12 (2011) 2319–2338. [2] G. Amori, L. Contoli, A. Nappi, Fauna d’Italia, (Italian Fauna) Mammalia III Erinaceomorpha, Soricomorpha, Lagomorpha, Rodentia, Calderini, Bologna, Italy, 2008. [3] P.R. Atsatt, Lycaenid butterflies and ants: selection for enemy-free space, Amer. Nat. 118 (1981) 638–654. [4] W.W. Becklund, C.M. Senger, Parasites of ovis canadiensis in montana, with a checklist of the internal and external parasites of the Rocky mountain bighorn sheep in North America, J. Parasitol. 53 (1967) 157–165. [5] F. Bianco, E. Cagliero, M. Gastelurrutia, E. Venturino, Metaecoepidemic models: infected and migrating predators, Int. J. Comp. Math. 89 (13–14) (2012) 1764–1780. [6] S. Busenberg, P. van den Driessche, Analysis of a disease transmission model in a population with varying size, J. Math. Biol. 28 (1990) 257–270. [7] R.J. Cowie, S.A. Hinsley, Feeding ecology of great tits (parus major) and blue tits (parus caeruleus), breeding in suburban gardens, J. Anim. Ecol. 57 (1988) 611–626. [8] J.T. Cronin, Movement and spatial population structure of a prairie planthopper, Ecology 84 (2003) 1179–1188. [9] T. Dewey, L. Ballenger, ovis canadensis (on-line), animal diversity web, accessed december 04, 2009. http://animaldiversity.ummz.umich. edu/site/accounts/information/Ovis canadensis.html. [10] M. Festa-Bianchet, Bighorn sheep, in: D. Wilson, S. Ruff (Eds.), The Smithsonian book of North American Mammals, Smithsonian Institution Press, Washington, 1999, pp. 348–350. [11] L. Gao, H. Hethcote, Disease transmission models with density-dependent demographics, J. Math. Biol. 30 (1992) 717–731. [12] I.D. Gauld, Evolutionary patterns of host utilization by ichneumonoid parasitoids (hymenoptera: Ichneumonidae and braconidae), Biol. J. Linnean Soc. 35 (1988) 351–377. [13] E.J. Gustafson, R.H. Gardner, The effect of landscape heterogeneity on the probability of patch colonization, Ecology 77 (1996) 94–107. [14] R.J. Guti´errez, S. Harrison, Applying metapopulation theory to spotted owl management: a history and critique, in: D.R. McCollough (Ed.), Metapopulations and Wildlife Conservation, Island Press, Washington, 1996, pp. 167–185. [15] I. Hanski, Single-species spatial dynamics may contribute to long-term rarity and commonness, Ecology 66 (1985) 335–343. [16] I. Hanski, A. Moilanen, T. Pakkala, M. Kuussaari, Metapopulation persistence of an endangered butterfly: a test of the quantitative incidence function model, Conserv. Biol. 10 (1996) 578–590. [17] M. Haque, E. Venturino, Effect of parasitic infection in the Leslie-Gower predator–prey model, J. Biol. Systems 16 (2008) 425–444. [18] E.P. Hoberg, G.S. Miller, E. Wallner-Pendleton, O.R. Hedstrom, Helminth parasites of northern spotted owls (strix occidentalis caurina) from oregon, J. Wildl. Dis. 25 (1989) 246–251. [19] P. Kareiva, Population dynamics in spatially complex environments: theory and data, Philos. Trans. R. Soc. Lond. Ser. B 330 (1990) 175–190. [20] P. Kerr, G. McFadden, Immune responses to myxoma virus, Viral Immunol. 15 (2002) 229–246. [21] P.J. Kerr, J.C. Merchant, L. Silvers, G.M. Hood, A.J. Robinson, Monitoring the spread of myxoma virus in rabbit oryctolagus cuniculus populations on the southern tablelands of New South Wales, Australia. II. Selection of a strain of virus for release, Epidemiol. Infect. 130 (2003) 123–133. [22] G. Lei, I. Hanski, Metapopulation structure of cotesia melitaearum, a parasitoid of the butterfly melitaea cinxia, Oikos 78 (1997) 91–100. [23] H. Malchow, S. Petrovskii, E. Venturino, Spatiotemporal Patterns in Ecology and Epidemiology, CRC, Boca Raton, FL, 2008. [24] J. Mena-Lorca, H. Hethcote, Dynamic models of infectious diseases as regulator of population sizes, J. Math. Biol. 30 (1992) 693–717. [25] J.C. Merchant, P.J. Kerr, N.G. Simms, G.M. Hood, R.P. Pech, A.J. Robinson, Monitoring the spread of myxoma virus in rabbit Oryctolagus cuniculus populations on the southern tablelands of New South Wales, Australia. III. Release, persistence and rate of spread of an identifiable strain of myxoma virus, Epidemiol. Infect. 130 (2003) 135–147.

44

E. Venturino / Mathematics and Computers in Simulation 126 (2016) 14–44

[26] A. Moilanen, I. Hanski, Habitat destruction and competitive coexistence in a spatially realistic metapopulation model, J. Anim. Ecol. 64 (1995) 141–144. [27] A. Moilanen, A. Smith, I. Hanski, Long-term dynamics in a metapopulation of the american pika, Amer. Nat. 152 (1998) 530–542. [28] M.J. Scoble, The Lepidoptera: Form, Function, and Diversity, Oxford University Press, Oxford, 1992. [29] R.T. Williams, J.D. Dunsmore, I. Parer, Evidence for the existence of latent myxoma virus in rabbits (oryctolagus cuniculus (l.)), Nature 238 (1972) 99–101. [30] E. Venturino, Epidemics in predator–prey models: disease among the prey, in: O. Arino, D. Axelrod, M. Kimmel, M. Langlais (Eds.), Mathematical Population Dynamics: Analysis of Heterogeneity, Theory of Epidemics, Vol. 1, Wuertz Publishing Ltd, Winnipeg, Canada, 1995, pp. 381–393. [31] E. Venturino, Simple metaecoepidemic models, Bull. Math. Biol. 73 (2011) 917–950. [32] J.A. Wiens, Wildlife in patchy environments: metapopulations, mosaics, and management, in: D.R. McCullough (Ed.), Metapopulations and Wildlife Conservation, Island Press, Washington, 1996, pp. 53–84. [33] J.A. Wiens, Metapopulation dynamics and landscape ecology, in: I.A. Hanski, M.E. Gilpin (Eds.), Metapopulation Biology: Ecology, Genetics, and Evolution, Academic Press, San Diego, 1997, pp. 43–62. [34] J. Wu, Modeling dynamics of patchy landscapes: linking metapopulation theory, landscape ecology and conservation biology, in: R. Wang (Ed.), Yearbook in Systems Ecology, English edition, Chinese Academy of Sciences, Beijing, 1994, pp. 97–116.