Modeling the Ebola zoonotic dynamics: Interplay between enviroclimatic factors and bat ecology

Understanding Ebola necessarily requires the characterization of the ecology of its main enzootic reservoir, i.e. bats, and its interplay with seasonal and enviroclimatic factors. Here we present a SIR compartmental model where we implement a bidirectional coupling between the available resources and the dynamics of the bat population in order to understand their migration patterns. Our compartmental modeling approach and simulations include transport terms to account for bats mobility and spatiotemporal climate variability. We hypothesize that environmental pressure is the main driving force for bats’ migration and our results reveal the appearance of sustained migratory waves of Ebola virus infected bats coupled to resources availability. Ultimately, our study can be relevant to predict hot spots of Ebola outbreaks in space and time and suggest conservation policies to mitigate the risk of spillovers.


Introduction
Zoonoses constitute 75% of emerging infectious diseases and pose a significant threat to public health [1,2]. In particular, the 2014 Ebola epidemic in West Africa has been the largest registered ever: as of December 2016 around 28,000 probable human cases with *75% mortality rates in laboratory confirmed patients [3]. In addition, Ebola virus (EV) decimates the great ape population, thus posing a conservation hazard, it represents a major threat worldwide through the importation of infections and its possible misuse as biological weapon [4], and, altogether, has dramatic economic [5], and humanitarian [6] consequences. Therefore, despite the promising advances to find a vaccine [7,8], understanding the factors and mechanisms underlying Ebola outbreaks and spillovers and developing predictive tools to prevent them are of major interest.
The 2014 EV strain in West Africa has been identified as Zaire's [9]. Notably, this strain originates thousand of miles away, in Central Africa. The source of the outbreak cannot be ascribed to the mobility of infected humans from Central Africa [10] but to a single zoonotic transmission event in Guinea: a human-bat contact [11]. a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 Evidence supports the idea that bats are the main EV reservoir regardless of the large diversity of EV hosts [12]. First, EV has been detected in tree-roosting bats [13,14]. Second, inoculation of the virus has revealed that bats can survive EV infection [15][16][17]. Third, novel filovirus have been only found in bats [18]. Fourth, there is a demonstrated connection of Ebola outbreaks with direct exposure to bats [19]. Finally, as for the expansion of the EV zoonotic niche, satellite telemetry has shown that bats are able to migrate thousands of kilometers annually [20]. Consequently, understanding the Ebola problem requires, on one hand, the characterization of the ecology of its main enzootic reservoir, i.e. bats, and, on the other hand, its interplay with seasonal and enviroclimatic factors [21][22][23][24][25][26] that drive and shape the bat migration patterns.
During the past years the Ebola modeling efforts have been intense [27][28][29][30][31][32][33]. However, these have traditionally focused on human epidemiology (the effect). Here, by focusing on the ecology and its interplay with the environmental conditions (the cause), we aim at shifting the current paradigm. In particular, we explore the feedback between environmental pressure and the dynamics of the EV enzootic niche. Our results show that a bidirectional interplay-i.e., how resources condition the population dynamics and, simultaneously, how the existing population conditions the available resources-may contribute to the expansion of the Ebola zoonotic niche.
The paper is organized as follows. In the Methods section we introduce and characterize theoretically an Ebola zoonotic compartmental model that accounts for the dynamics of the bat population. In the Results section we show by means of numerical simulations how resources variability and seasonality drives bat migration and EV spreading. Finally, the implications and the main conclusions of our research are detailed in the Discussion section.

Ebola zoonotic compartmental model
Here we propose a SIR (Susceptible-Infected-Recovered) compartmental zoonotic model for bats. Fig 1 summarizes our modeling approach. Notice that on top of susceptible, B S , infected, B I , states we also consider a recovered (from infection) state, B R . Our hypothesis is based on data about the dynamics of filovirus infection on bats: most of the infected individuals are older juveniles [23] and bats can survive infection [15][16][17]. A reasonable conjecture to explain the demography of infection is by invoking recovery. During the recovery phase we assume that bats cannot get infected and/or transmit the disease (see Discussion).
In our model we consider the following processes to account for the dynamics of the bat population and the EV infection: birth (rate b K ), death (rate c), competition for resources (K), EV transmission (rate e), recovery from infection (rate d), and retrieval to the susceptible state from a recovery state (rate f). We notice that the obtained analytical results are valid in case d = f = 0, that is, when the recovery state is suppressed. As detailed below, we also explore the mobility of bats by means of diffusive terms, D K . In addition, we consider an equation for the dynamics of resources, K (see Fig 1). By overlooking the mobility terms for the time being, our model for the time evolution of the density of bats in different states reads, where B = B S + B I + B R stands for the total population of bats. Bats mobility is considered in our modeling approach by including in the r.h.s. of Eqs (1)-(3) terms of the form +D K r 2 B Z , where r 2 ¼ @ 2 @x 2 þ @ 2 @y 2 and Z stands for S, I, or R. That is, we consider the simplest form of nondirected mobility, i.e. diffusion. Other forms of transport are certainly possible, e.g. convective terms, but those require additional hypotheses about the causes driving bats migration.
According to the experimental data, EV infection does not modify bats' physiology. Thus, we assume that the birth and death rates (or the mobility coefficient) do not depend on the state of infection. Yet, as detailed below, we assume that the birth rate and the mobility coefficient depend on the available resources, K. Such functional dependence (see details below) effectively summarizes the effect of the environmental pressure in our model. The parameter λ 2 [0, 1] reflects the lack of information about the possibility that bats are born either EV free, λ = 1, or infected, λ = 0. The demographic data about Marburg filovirus infection [23] indicates that most of the infected individuals are older juveniles and not newborns. Consequently, the former possibility, λ = 1, seems to be the more plausible. Note that the absolute value, |b K − c|, in the equations is required to make the Verhulstian quadratic competition term, *B 2 , always negative regardless the sign of the growth rate, a K = b K − c. As for the equation for the resources (carrying capacity), γ stands for the rate of depletion of resources by the bats, and r accounts for the rate at which the resources naturally return to their "bare" value, K 0 . The latter corresponds to the value of the carrying capacity in the absence of bats.

Simulation scheme and parameter values
Our numerical simulations use the FTCS scheme [34] in an hexagonal lattice. Each hexagonal tile represents a surface of 25 km 2 (unit length: l c = 5 km) and the temporal time step is t c = 1 day. Our simulations span a surface of 2.5 Á 10 4 km 2 and 10 years such that a) the simulated spatial domain is large enough to justify the spatial variation of resources, and b) several (annual) seasonality periods are included to avoid transient effects due to initial conditions. The value of the parameters are summarized in Table 1. We try to implement values as realistic as possible. Yet, we must point out that most parameters either have not been measured and/or their characterization has been made for different bat species under different conditions. In any case, the ultimate objective of this study is not to fit factual data about EV infection but to propose a mechanism to explain the interplay between environmental conditions, EV infection and the bat migratory dynamics.

The interplay between resources and the zoonotic population buffers EV infection
In order to analyze the role played by the interplay between resources and the bat population we analyze first the difussionless case, D K = 0. By adding Eqs (1)-(3) the equations for the total bat population, B, and the resources reduce to, That is, a logistic growth equation combined with a relaxational dynamics for the resources.

Parameter Units Value
Birth rate [35] Infection rate [16,23,36,37]  We estimate the onset of infection using two different, complementary, approaches. On one hand, we compute the basic reproduction number, R 0 using the next generation method [41][42][43][44][45]. Eqs (1)-(4) has a single, non-trivial (a K > 0), infection-free stationary solution: The first term on the l.h.s. of Eq (7) collects the transmissions terms and the second one the transitions terms such that [41], The above expression is valid if The condition that sets the boundary that separates these behavior is given by  (5) and (6) has two possible stationary states (attractors): {K st , B st }: {μK 0 , μK 0 } and {K 0 , 0}. By estimating the characteristic polynomial of the Jacobian and implementing the Routh-Hurwitz stability conditions (negative eigenvalues), we conclude that the former attractor is stable is a K > 0 whereas {K 0 , 0} is stable if a K < 0. Thus, if the growth rate is positive, b K = b 1 > c, the carrying capacity relaxes to a value μK 0 < K 0 that determines the maximum density of bats that can be maintained by the resources when taking into account consumption. On the other hand, if the growth rate is negative, b K = b 2 < c, the bat population is extinguished and the resources relax to their bare value, K 0 .
The aforementioned attractors can be split into a higher dimensional phase space in terms of the different states. Thus, Eqs (1) Thus, the stability analysis reveals that the EV infection develops as long as there is a positive growth rate and the infection rate is larger than a critical value, e c . This is the same condition we obtained in terms of the basic reproduction number using the next generation method. As expected the larger the the recovery rate, d, diminishes the possibility of an Ebola outbreak. Also, if bats are born infected, λ = 0, or the resources, K 0 , decrease it favors Ebola infection. As mentioned above, the case λ = 0 seems unplausible according to the demography of infected cases in bats and we assume hereafter that λ = 1, i.e. bats are EV born-free. As for the role played by the interplay between resources and bats, notice that 0 μ 1; if either bats do not condition the amount of available resources, γ = 0, or the resources are quickly replenished, γ ( r, then μ ' 1. Therefore, notably, the consumption of resources by the bats effectively buffers EV infection. Note also that, if either bats consume rapidly the resources, γ ) r or the replenishment rate vanishes, r = 0, then μ ' 0 and the bat population, as well as the resources become extinct.

The environmental pressure modulates EV infection
We explore the effects of the environmental pressure in the dynamics of the bats population assuming that if the available resources are below a certain value, K Ã , then the survival of the population is challenged. Thus, we assume that if . According to the discussion above, in the former case the attractor corresponds to {μK 0 , μK 0 } while in the latter the attractor corresponds to {K 0 , 0}. Since the (stationary) amount of resources is bounded in the interval (μK 0 , K 0 ), we assume that K Ã belongs to that interval. As shown in Fig 2, this mechanism compels bats and resources to chase attractors with an evanescent stability. This in turn induces an oscillatory behavior and a new stable attractor develops. We notice that if K Ã = 2 (μK 0 , K 0 ) then the environmental pressure does not play any role since the induced attractor is unstable and either {μK 0 , μK 0 } or {K 0 , 0} are stable depending on the sign of the growth rate a K .
As for the location of the new attractor, we notice that By using this functional form of b K instead of Heaviside step functions we avoid the difficulties derived by its non-continuous behavior. Hence, by solving Eqs (5) and (6) using this definition we found that, Thus, in the limit n ! 1, we obtain fK st ; B st g ¼ fK Ã ; mðK 0 À K Ã Þ 1À m g. A stability analysis (Jacobian method) reveals that the chasing dynamics of the evanescent attractors is relaxational (the real part of the eigenvalues is negative). Consequently, the induced oscillations are not sustained but damped, i.e. a focus-like behavior, and the system eventually relaxes to this fixed point that is located in the pathway connecting the evanescent attractors. The same procedure can be implemented in the case of Eqs (1) and (4) Fig 3 shows, as a function of the dimensionless parameters κ, μ, and K 0 /K Ã , the regions for which the environmental pressure plays different roles.
We can further analyze the infection dynamics by calculating the basic reproduction number. In that case, the transmission and transition terms read, Hence, In our simulations, see Table 1, R 0 % 1:23. As for the onset for infection, i.e. the condition  R 0 > 1, it reduces to, e > e c + |b 1 + b 2 − 2c|/(2K Ã ). Thus the threshold for infection progagation, as estimated by the basic reproductive number, is more restrictive than the condition found using the stability analysis (see Discussion).

Seasonality and climate variability: Sustained migration waves
As for the spatiotemporal dynamics of our model, if the initial condition is space-dependent then mobility (diffusion) leads to transient migratory waves. Yet, if no further considerations are taken into account, every location reaches the same stationary state and mobility, eventually, does not play any role in the dynamics of infection. However, temperature, precipitation and other environmental indicators display periodicity and, in addition, environmental factors are also affected by the geographical location. As shown below these features, in combination with mobility terms, produced sustained migration waves. In our modeling approach we  Table 1): κ = 17/28, μ = 1/2, and K 0 /K* = 3/2. https://doi.org/10.1371/journal.pone.0179559.g003 implement the seasonality and climate variability by considering the following, simple, spatiotemporal dependence of the bare carrying capacity, K 0 , where hK 0 i stands for an average carrying capacity, 0 < Γ < 1 the fractional variation that we assume constant, ω the seasonality temporal frequency, and F(r) is a phase that accounts for the spatial variability: different regions are subjected to the same periodic behavior of the resources but asynchronously. Here, we will consider the case of two well differentiated regions such that F(r) = πδ r,r 0 , where r 0 indicates all points that belong to one of the regions (see Fig 4). Under these conditions, if D K 6 ¼ 0, traveling waves of migrating bats develop. As for D K , we hypothesize that bats migration must be coupled to environmental pressure: we argue that bats migrate if the lack of resources threatens their survival. In our model we implement this response to the environmental pressure by means of an non-homogeneous diffusion coefficient: D K = D 6 ¼ 0 if K < K Ã and D K = 0 otherwise (see Fig 2). Fig 4 shows the oscillatory behavior of bats population as a function of space and time and reveals their migratory dynamics chasing resources. With the parameters used in our simulations the percentage of EV infected bats oscillates between *16% and *19% depending on the season; the season when the amount of resources is minimal/maximal having the largest fraction of EV infected/healthy bats.

Discussion
From the point of view of the required framework to study the outbreaks and spillovers associated to zoonotic diseases, one possible approach is to use models able to deal with the complex behavior arising from the interactions of individual units (human, animal,. . .). Moreover, in the particular case of Ebola, understanding the various levels of complexity of the problem ultimately requires anthropological and ecological considerations. Thus, the high mobility of humans in Africa through porous borders, cultural practices, socioeconomic conditions, or the high diversity of Ebola zoonotic carriers have been recognized as important elements [47]. Rule-Agent based models can provide such detailed level of description and establish predictive correlations among relevant factors. However, to identify fundamental mechanisms is difficult within such framework. Another option is to use a reductionist modeling approach, less precise from the viewpoint of the predictive character, but able to pinpoint precisely underlying mechanisms. Herein we have chosen this second alternative by using a compartmental model to address a specific question: how do the interplay between seasonality, climate variability, and environmental pressure drive the migration patterns of bats and the Ebola infection dynamics in their population?
Our model first assumes a bidirectional coupling between the available resources, effectively described by means of a carrying capacity variable, and different states of the bats population such that they reshape their respective equilibrium values. Environmental pressure is taken into account by hypothesizing a threshold below which the survival of the population is challenged. We have shown that environmental pressure induces new equilibrium states and an oscillatory behavior. Moreover, our results reveal that the environmental pressure modulates the population dynamics and, given an infection rate, can either promote or hinder the infection levels depending on the balance of three factors: the ratio between the rate of resources consumption by bats and the replenishing rate, the value of the carrying capacity threshold setting environmental pressure with respect to the "bare" resources, and a dimensionless parameter that accounts for the growth rate of the bat population and their ability to recover from infection. Finally, when seasonality and variability are considered as factors regulating the available resources and we assume conditional migration as a response mechanism to the population survival, we have shown that sustained waves of migratory bats develop. The migratory waves are correlated with the dynamics of the resources and our model suggests that infection is more probable when resources are low. The reason for this behavior is an increased competition for resources during those periods.
As a matter of discussion, here we have considered a minimal approach in terms of the possible states of infection: susceptible, infected and recovered. We have included a recovered state based on serological data of bats infected with Marburg filovirus [23]and physiological studies about Ebola infected bats [15]- [17]. We argue that if the death rate of bats is no changed because of the infection and adults are not the population sector with the largest number of infection cases then a recovery state can explain these observations. We acknowledge that this hypothesis could be eventually rebutted or that other serological states are possible (e.g. inmune, exposed). We notice that by presenting a more complex model with additional infection states we would rely on un-tested hypotheses and also on the estimation of additional parameters that are unknown at this point. Also, it is worth noting that some of the main conclusions of our study do not depend on this: the dynamics that set the basic interplay between the bat population and the environment, as given by Eqs (5) and (6), is independent of having a recovery state. In addition, the expressions that determine the concentration of different states or the phase diagram plot (Fig 3) are valid even in the case that the recovery rate is zero.
Our results also reveal an interesting effect in regards of the reproduction number due to the environmental pressure. In the absence of the environmental pressure, the condition found for the onset of infection is the same regardless whether it is obtained either using the Jacobian approach (stability analysis of the stationary infection-free solution) or through the next generation method (basic reproduction number). However, when the birth rate is a function of the available resources (environmental pressure) then the threshold obtained depends on the approach. This issue has been pointed out in other studies [41,48]: while a linear stability analysis (Jacobian method) reveals the stability condition of the disease-free state it does not necessarily provides a value of R 0 that is biologically meaningful. Likewise, the threshold condition found using the basic reproduction number does not necessarily imply an asymptotic stability condition of the disease-free state.
A number of studies have shown that in the context of population dynamics and infectious diseases the diffusion may play a critical role and leads to the emergence of patterns [49-51]. The underlying mechanism is based on the Turing instability that relies on the existence of distinct diffusion scales for different species [52]. In the case of the Ebola infection in bats, data does not suggest that their physiology is altered due to the infection and, consequently, we have not consider different migratory capabilities depending on their infection status. As a result, the emergence of spatial patterns of infection can be discarded in our model and diffusion plays merely an homogenizing role that spreads EV infection when the environmental pressure forces migration. However, we notice that possible extensions of our model might consider a multiple species approach given that the zoonotic niche of Ebola is larger than originally expected [12]. In that case, the migratory capabilities would be certainly different depending on the host and this would open the possibility of finding a pattern emergent behavior.
Our study can be relevant to predict hot spots of Ebola outbreaks in space and time and may also suggest conservation policies to mitigate the risk of spillovers. This depends, ultimately, on the possibility of calibrating our model with the right, realistic, parameters and points out the importance of quantitative ecology and climatology approaches. As for possible extensions of this study, our results suggest the importance that the role played by stochasticity in climate and seasonality since large fluctuation events, while having low probability, may lead to a significant increase in the infected population of bats and therefore of an Ebola outbreak. Work to explore these factor is in progress.