Controlling infection in predator-prey systems with transmission dynamics.

We propose in this paper a prophylactic treatment strategy for a predator-prey system. The objective is to fight against the propagation of an infectious disease within two populations, one of which preys on the other. This propagation is modeled by means of an SIS (susceptible-infectious-susceptible) epidemic model with vital dynamics and infection propagation in both species through contact and predation, including mortality rates in both populations due directly to the disease. Treatment strategies are represented by new parameters modeling the uptake rates in the populations. We analyze the effect of various treatment strategy scenarios (prey only, predator only, or both) via their uptake rates and possible cost structures, on the size of the infected populations. We illustrate if and when applying such preventive treatments lead to a disease prevalence drop in both populations. We conduct our study using an optimal control model seeking to minimize the treatment cost(s), subject to the transmission dynamics and predator-prey dynamics.

Introduction

The decrease in Earth’s natural resources and the depletion in the ecology of wildlife increase rapidly because of the massive expansion of human urbanization areas, the water/air pollution, and deforestation. Consequently, there is a noticeable increase in the rate of infectious diseases that are affecting the wild animal populations (for example see (Leighton, 2011)). In this regard, numerous studies have discussed the dynamics of the population of wildlife animals with respect to infectious diseases. For instance, a significant study in (Wiethoelter, Beltrán-Alcrudo, Kock, & Mor, 2015) sheds light on infectious diseases in wildlife and livestock. In contrast, the pathogen outbreak and long-term prevalence in populations as in (Becker & Hall, 2014) is developed in an explicit model to discuss the influence of the interaction of three provisions: host demography, contact behavior, and the immune defense.

There is a considerable and significant effort to study the natural animal population dynamics by using disease modeling, as well as the impact on the human population. Two noticeable mathematical analyses are presented in (Han, Ma, & Hethcote, 2001), based on predator-prey models with stable equilibria. The predator’s population could acquire an infection from the prey population through the predation process. Along these lines, the epidemiological model in (Chauhan, Misra, & Dhar, 2014) looked at the stability analysis of an SIR model (susceptible-infected-recovered), with and without a preventive treatment, depending on the basic reproduction number . An SIS epidemiological model that is discussed in (Zhou & Liu, 2003), presents the dynamical behavior of an impulsive system to characterize the pulse preventive treatment process.

Our research is inspired by several existing diseases and pathogens. The epidemics caused by pathogens extensively affect the fauna of many areas around the world, including North America. The work in (Sapp et al., 2016) studies a roundworm (Baylisascaris procyonis) which infects many species like raccoons and birds. This type of roundworm can produce a fatal neural larva in the host’s blood and is known to be transmissible to humans as well (Bauer, 2013). This pathogen can be lethal for transient hosts (like mice, humans, etc.) and, to a lesser extent, for definitive hosts (raccoons, bears, panda bears). Chronic wasting disease in deer, and occasional cows, is fatal to most deer. It is unclear if it is transmitted to predators or humans (Belay et al., 2004; Oraby, Vasilyeva, Krewski, & Lutscher, 2014). In a “worst-case” scenario the infection could move to apex predators and/or humans. Finally, toxoplasmosa gondii (T. gondii) (Bevins et al., 2012; Flegr, Prandota, Sovičková, & Israili, 2014) is a parasite common for instance in populations of deer and cougars, and is transmitted by both predation and contact. It is also common in domestic felines from which it can be transmitted to humans - especially during pregnancy. However, T. gondii is suspected to be fatal to sea otters (see (Miller et al., 2002), a study based on data and correlation discovery between T. gondii presence in brains of dead sea otters). In (Miller et al., 2002) the authors infer that the otters may acquire the parasite from coastal freshwater runoffs carrying the parasite, but there is no specific information on inter-species transmission by predation: it is known that otters eat rodents and birds, and T. gondii is highly prevalent in rodents, thus they could acquire the parasite by predation.

The preventive treatments are introduced in our epidemic model as controls in the form of treatment uptake rates. The history of optimal control problems goes back to the 50s (Leitmann, 1962), when Pontryagin et al., studied the calculus of variations depending on Euler-Lagrange equations. Since then, extensive studies have been conducted to enrich the concept of the classical optimal control problem with necessary conditions such as the Pontryagin Maximum Principle (Barbu, 2012). In most cases, when the differential system of the control problem is nonlinear, the analytical solution is too difficult to be reached in closed form. Hence, a numerical approach is needed. The numerical method we use is the steepest descent method, which belongs to the class of indirect method (see (Kirk, 2012; Wang, 2009)).

The rest of this paper is organized as follows. In Section 2, we discuss the predator-prey model considered in this paper. In Section 3, the initial model is modified to introduce treated compartments in both populations and infection specific death rates. In Section 4, we introduce an optimal control formulation which aims at highlighting scenarios of minimizing the cost of treatment while reducing the infection. In Section 5, we present numerical results for the optimal control formulation. Finally, we conclude in Section 6 with some discussions on the results and future work ideas.

Description of the predator-prey model with infection

The predator-prey model with infection

We analyze a mathematical model that describes an infectious disease in predator-prey populations by building on the model proposed in (Han et al., 2001), with some important modifications. We consider the prey’s population to be of size and to consist of the susceptible prey group and the infected prey group . Correspondingly, the predator’s population is of size and consists of the susceptible predator group and the infected predator group . The differential equations that reflect the SIS predator-prey model with standard incidence, vital dynamics and predation are given by:where and . The parameters used to describe the model are sum up in Table 1 below.

Table 1

The parameters of the SIS predator-prey model.

Parameter	Description
b_1	The prey’s natural birth rate
d_i, i=1,2	The natural death rate in each population
r_1=b_1−d_1	The net growth rate in prey
K_1	The carrying capacity of the prey
a	Predation rate
k	The efficiency of predation
β_i, i=1,2	The transmission coefficient
γ_i,i=1,2	The recovery rate in each population
α	The transmission coefficient from prey to predator
μ_i,i=1,2	The death rates in population i due directly to infection

Equivalently, we can write the system as:

In the sequel, we will focus on (1) as it simplifies the presentation. The solutions of the system (1) can be found in the nonnegative orthant:

Equilibrium states of the system

We highlight all possible equilibrium states of this system, as they arise following a closed-form analysis using straightforward substitution methods. We will denote a generic equilibrium point of (1) by: . The immediate equilibria to deduce are , and . Outside of these two evident guesses, let us use equations (1), (3)) from system (1) to isolate:and

We see that if both and in equilibrium, then we are left with

Looking at the possibility of equilibria in which either: and , or and leads to contradictions, which means that the only other possible equilibria are the endemic ones of the type

To illustrate the emergence of numerically, we present in Fig. 1 the trajectories of the system (1) with the parameters given in Table 2 and with as an initial point.

Fig. 1

The solution trajectory of system (1) reaching with values from Table 2 starting from and time days. In this case, the conditions above are giving the following equilibrium point .

Table 2

Specific parameters of the SIS predator-prey model with mass action incidence.

b_1	d_1	d_2	K_1	a	k	β_1
0.0367	0.03	0.0023	1000	0.0016	0.0056	0.1
		γ_2	α	μ_1	μ_2
0.0334	0.02	0.01	0.025	0.00334	0.00034

From this figure, we can see that the system reaches an endemic equilibrium point where the disease in both populations persists. Therefore, we aim to introduce an optimal preventive treatment plan to this model by treating the prey population to minimize the total infected populations in both prey and predator.

Introducing a preventive treatment into the SIS predator-prey model

We propose that the preventive treatment is offered to the prey and the predator. The treated subgroups have to now be modeled in the system (1) via new compartments. Thus, (1) can be rewritten as follows:where now and . The parameters and represent the treatment uptake rates and the variables and are the “vaccinated” compartments (for details see (Zhou & Liu, 2003)).

By considering fixed parameter values in the equilibrium conditions above we obtain additional expressions for the endemic equilibria as follows:

In Fig. 2, we show the evolution of the system for fixed values of and in three different cases. First, we consider the model with one treated compartment only for prey (so ), second we consider a model with one treated compartment only for predators () and lastly, we consider the model with treated compartments for both prey and predator.

Fig. 2

The solution trajectory of system (2) with values from Table 2 starting from and time days for several values of (upper panel), (middle panel), (lower panel). Final values are given in Table 3.

Table 3

Size of each population from the system (2) with parameter values from Table 2, starting from after days for several values of .

	N_1(T)	I_1(T)	V_1(T)	N_2(T)	I_2(T)	V_2
u_1=u_2=0	281.87	83.23	–	2.28	2.25	–
u_1=0.2,u_2=0	258.14	0.00	219.75	2.94	1.08	–
u_1=0,u_2=0.2	265.02	77.96	–	2.36	1.03	1.32
u_1=0.2,u_2=0.2	243.20	0.00	206.94	3.01	0.00	2.97

This preliminary experiment is done in an ideal case in the sense that we do not consider the cost of treatment. In the following section, we consider an optimal control model to include a cost function forms for the treatments to be implemented. We will specifically look at treatment costs for the prey population first, then at treatment costs for the prey populations second, presenting our specific considerations for these two scenarios.

Optimal control of infectious populations

Optimal control formulation

We can now reformulate system (2) as part of an optimal control problem as in (Barbu, 2012). The general idea in the control theory is to minimize some cost function governed by a system of ordinary differential equations. A general form of this problem is as follows:where is the Lagrangian functional, is a final time cost, and .

In our case, , the ordinary differential equation is the state system (2) while our performance measure with respect to the infection is given by:

Furthermore, we assume that the cost of treatment is given by a quadratic function of the form

The quadratic cost, with a U-shape, is a classical example of an average cost function. Close to zero the cost is relatively low to encourage the investment while it begins rising more rapidly as the investment becomes larger. Hence, the cost function becomes

We presume that the initial condition of the system (2) and the initial time are clearly identified with . On the other hand, we suppose that the control vector is defined as:

Hence, we can rewrite the problem (3) with respect to our SISV model in the following way:

In the following section, we compute the numerical solution of the optimal control (4).

Numerical solution of the optimal control problem

A classical approach to compute the numerical solution of the optimal control (4) consists of solving the first-order optimality conditions of the problem derived via the Pontryagin’s Maximum Principle. This is known as an indirect approach in the difference with the direct approach which consists of a discretization of the control problem. The numerical determination of the optimal control and optimal trajectory of a nonlinear control problem with necessary conditions can be found for instance in (Kirk, 2012) in detail. A numerical method to solve this problem is the steepest descent (SD) method whose algorithm is presented in (Kirk, 2012) and (Wang, 2009). We used this method in similar work to solve an optimal control of a social norm game in (Jaber & Cojocaru, 2018) and to solve the optimal control of a preventive treatment game in (Cojocaru & Jaber, 2018) to maximize the number of players who choose to vaccinate in an epidemiological model.

The Hamiltonian which is crucial converting (4) into a two-point boundary-value problem is defined as follows:where is the costate. Consequently, we can define the necessary conditions for the problem (4) as follows:

Combining the optimality conditions and the transversality condition, we can formulate the following two-point boundary-value problem:

For a given u and ρ, the first two conditions constitute a classical ordinary differential equation with an initial condition, while for u and x fixed, the third and fourth conditions constitute a backward differential equation. Both can be solved using a classical ordinary differential equation solver. Following (Wang, 2009), we use MATLAB© to solve the differential equations. The constraint on the control is handled using a Lagrangian penalization.

We now consider the resolution of the optimal control problem introduced above with a required precision of . For the numerical experiment we use a scaled version of the system where the variables are .

Numerical solution of the optimal control of the predator-prey SISV

We now run the numerical approach described previously on the optimal control formulation of the 2-species infectious disease model. In particular, we observe how the endemic equilibrium reacts to our treatment strategies. In all of the following experiments, we consider the equilibrium above as an initial point at for the control problem (4) and we consider a time horizon of days.

Constant uptake rates of treatment

We first study controlling the infection assuming enough resources, to sustain an average constant uptake rate of treatments: over a 1 year time horizon.

For an easier read of our results, we present first our results in tabulated form (Table 4), where we highlight the new, controlled, endemic equilibria arising from each of the three cases highlighted above and we note that the simulated results were already presented in Fig. 2.

Table 4

Comparison of the infection and the costs for the optimal control vaccine formulation starting from after days.

Cases	N_1	I_1	V_1	N_2	I_2	V_2
no treatment	281.87	83.23	–	2.28	2.25	–
u_1=0.2,u_2=0	333.51	0.00	284.51	2.70	1.00	–
u_1=0,u_2=0.2	277.24	81.79	–	2.36	1.34	0.99
u_1=u_2=0.2	327.99	0.00	279.76	2.78	0.04	2.69
model	Ψ_u	Ψ_I	Ψ
no treatment	–	467.97	467.97
u_1=0.2,u_2=0	7.30	229.11	236.41
u_1=0,u_2=0.2	7.30	375.67	382.97
u_1=u_2=0.2	14.6	116.44	131.04

From these results, we can draw a few conclusions. Overall, treatment reduces the infection in both the preys and the predators’ populations. Treatment at a fixed cost allows us to eliminate the infection in the preys. However, the infection is persistent in the predators, especially when treatment is only applied to the prey population (Case 1 above). As expected, the scenario in which both populations are treated achieves the lowest infected population fractions, but the vaccination cost is higher than the other scenarios.

Optimally controlled uptake rates of treatment

Following our analysis so far, we see that it confirms the need to search for an optimal uptake strategy (i.e. finding an optimal and ). As before, the base scenario is the endemic equilibrium state which corresponds to the case without any treatment strategy. The next three cases apply the optimal control techniques described in the previous section to find and ), costs and , as well as the states in three different cases:

The case of (OC of );

The case of (OC of );

The case of (OC of )

We tabulated our results below (see Table 5), whereas in Fig. 3, we show the evolution of the populations using the optimal strategy for the vaccine rates in the three scenarios described above.

Table 5

Comparison of the infection and the costs for the optimal control vaccine formulation starting from after days.

Cases	N_1	I_1	V_1	N_2	I_2	V_2
no treatment	281.87	83.23	–	2.28	2.25	–
OC of u_1	332.64	0.00	256.78	2.58	1.02	–
OC of u_2	280.68	82.87	–	2.29	2.23	0.03
OC of u_1,u_2	330.90	0.00	255.30	2.62	0.20	2.08
model	Ψ_u	Ψ_I	Ψ
no treatment	0	467.97	467.97
OC of u_1	4.91	286.19	291.10
OC of u_2	0.0029	464.79	464.80
OC of u_1,u_2	5.08	232.71	237.79

Fig. 3

The solution trajectory of system (2) with values from Table 2 starting from after days for several scenarios from the top to the bottom: with only , only , and both controls. Final values are given in Table 5.

From these results, we can see that the optimized control achieves a tradeoff between the reduction of the infection and reducing the costs. As in the previous experiment, the vaccination of the prey seems the most efficient to reduce the infection in both populations.

Conclusion

In this paper we studied a well-known model of predator-prey populations with an infectious disease that can affect both species, and where transmission happens both by contact and by predation. We highlight two frameworks in which possible treatments could be analyzed, once they are introduced to either populations or to both. In the first framework, we chose to ignore costs of treatment implications, assuming that there are enough resources and treatments available to both species. The expected conclusion here was that treating both species leads to lower levels of infection, at the highest cost.

In the second framework, we used an optimal control formulation which allows to minimize the treatment costs and the infectious levels in both populations. We find that there is a tradeoff between reducing the infection in both species versus controlling the treatment cost.

In both frameworks, with fixed rates and with optimal rates, vaccinating the preys only is an efficient way of also lowering the infection among the predators, while keeping the cost away from maximal levels.

This model is immediately amenable to further investigations. First off, to inform decision making, the model needs to be applied to real populations. Our inspiration for the numerical values were populations of white-tale deers and cougars in British Columbia, Canada (Negri & Hornocker, 2010; Robinson, Wielgus, & Gwilliam, 2002). The tradeoff between treating only the prey rather than both populations would depend on the specific populations, as well as on the availability and ease of implementation of the treatments in question. For this purpose, our model can also be generalized, for instance, including a time delay for the reproduction of predators after predation. Although this type of extension modified the equilibria and the stability of the system (Xiao & Chen, 2001), it is compatible with our vaccination strategy and alters only the strategies involving a vaccination in the predators. The differential equation in the optimal control problem might become harder to solve. Second, there are immediate multiyear, budget-type scenarios that can be simulated. For instance, the treatment could be seasonally applied, depending on the pathogen, which may necessitate the inclusion of resource allocation scenarios over several years. Last but not least, cost functions need not be as regular as assumed here (note we assumed a quadratic function, which is convex and twice differentiable, which makes the numerical method for solving the optimal uptake rate somewhat simpler). Another possible and widely used cost function is of logistic type, which generally is neither convex nor concave, thus making solving for the optimal uptake rate more difficult.

Declaration of competing interest

None.