“Does a Respiratory Virus Have an Ecological Niche, and If So, Can It Be Mapped?” Yes and Yes

Abstract

Although the utility of Ecological Niche Models (ENM) and Species Distribution Models (SDM) has been demonstrated in many ecological applications, their suitability for modelling epidemics or pandemics, such as SARS-Cov-2, has been questioned. In this paper, contrary to this viewpoint, we show that ENMs and SDMs can be created that can describe the evolution of pandemics, both in space and time. As an illustrative use case, we create models for predicting confirmed cases of COVID-19, viewed as our target “species”, in Mexico through 2020 and 2021, showing that the models are predictive in both space and time. In order to achieve this, we extend a recently developed Bayesian framework for niche modelling, to include: (i) dynamic, non-equilibrium “species” distributions; (ii) a wider set of habitat variables, including behavioural, socio-economic and socio-demographic variables, as well as standard climatic variables; (iii) distinct models and associated niches for different species characteristics, showing how the niche, as deduced through presence-absence data, can differ from that deduced from abundance data. We show that the niche associated with those places with the highest abundance of cases has been highly conserved throughout the pandemic, while the inferred niche associated with presence of cases has been changing. Finally, we show how causal chains can be inferred and confounding identified by showing that behavioural and social factors are much more predictive than climate and that, further, the latter is confounded by the former.

1. Introduction

Recently there has been debate [1,2,3,4,5,6] as to whether Species Distribution Models (SDM) are appropriate tools in the study of the COVID-19 pandemic. Of course, such a debate raises the important question of when and under what circumstances an SDM, or an Ecological Niche Model (ENM), are likely to be valid and/or useful in the study of disease in general? Just which diseases, or aspects of diseases, can be usefully studied using SDM/ENMs? Which pathogens have ecological niches and which do not? In addition, if some do not, why do they not? Taking the broad view: If ecology is the study of the relations between organisms, both among themselves and with their environment, and an ecological niche is the full set of biotic and abiotic factors that favour the presence of an organism [7,8], then clearly all disease pathogens must have an ecological niche. The question, rather, is: how is that niche to be characterised and quantified? What are the appropriate niche dimensions and can a meaningful and useful ENM and, consequently an SDM, be constructed?

Although there is controversy as to the application of ENMs and SDMs to the case of COVID-19, there have been varying degrees of success in applying them to other human diseases. In particular, to zoonoses [9,10], where non-human species are involved in the transmission cycle. The principal applications have been to the creation of ENMs and SDMs for important agents in the transmission cycle of the zoonosis, such as known vectors [11] and hosts [12], considered as risk factors, as well as the pathogen itself [13]. In all these cases the corresponding niche variables were abiotic—climate and environmental layers—in spite of the fact that biotic factors clearly play an important role as niche variables for the disease itself. Biotic factors were successfully incorporated into ENMs and SDMs for zoonoses in [14,15,16,17,18,19,20] Importantly, the capacity to include biotic factors into ENMs and SDMs led not only to the creation of more predictive and explainable models but also, for example, to the prediction and confirmation of multiple new hosts of several zoonoses [16]. However, all these models, independently of the inclusion of biotic factors, were based on presence–absence data and on the implicit assumption that the system was in “equilibrium” [21].

In all of these cases of human disease there has been a component that is considered “ecological”, either in the targets of the ENM/SDM, such as vectors or hosts of a zoonosis, as well as in the niche variables themselves, such as climatic and other environmental variables, as well as potential hosts and vectors. However, in the case of a disease such as COVID-19, which is based on human-to-human transmission, there is a question as to the relevance of such factors. Additionally, abundance, rather than just presence–absence is a fundamental concern. Finally, it is clear that for epidemics or pandemics of this nature in no way is an equilibrium assumption valid. It is these criticisms, among others, that have led to the questioning of an ecological ENM/SDM approach and to the conclusion that an epidemiological rather than ecological approach is to be preferred, which leads us to a further question: when is an ecological versus an epidemiological approach more appropriate? The difference between the two is stated quite clearly in the epidemiological literature [22,23]. Put simply, ecological models are naturally associated with the “where” of a disease, while epidemiological models are more concerned with the “who”. This difference is aptly captured in [24] with the question: “Is it better to have a heart attack in the United States or Canada?”, the emphasis being that this type of question is of an ecological nature, whereas the question: How many people out of a cohort of individuals diagnosed with COVID-19 will subsequently die? is an epidemiological question. Indeed, this ecological “where” perspective has a long history in the social sciences as a whole [25,26]. From a modelling viewpoint, the difference between the two approaches lies in what statistical ensemble will be used to draw inferences. Epidemiological models use an ensemble of “individuals”, while ecological models use an ensemble of “places”, where each place is associated, either explicitly or implicitly, with a population. The “places” can vary depending on context: from countries and other political units, to eco-regions, any fixed-area spatial cells, and to pixels on a raster.

Additionally, there is a subsequent question of “when”, that is relevant for both the ecological and epidemiological points of view. Standard SIR-type modelling [27], for instance, is associated with the “who” and the “when”, with the “who” being associated with a fixed-number of states—susceptible, infected, etc. that are associated with individuals [28]. It is the “when” and, relatedly, “how many” that has dominated the modelling of the COVID-19 pandemic, with standard SIR models, and variations thereof [29,30,31], playing an important role [32,33]. However, more sophisticated techniques, such as deep learning, have also been used [34]. All these methods however, are based on modelling the time series of the data to be predicted, such as cases or deaths. In no way can they account for the high degree of multi-factoriality involved in the evolution in space and time of the pandemic. In other words, they cannot account for the “whys” that accrue from the direct and indirect causes of this evolution. Thus, although such dynamical models can be extended to consider “place” by constructing a SIR-type model for a particular place [35,36]; such models, however, do not account for the distinguishing features of that particular “place”, as does an SDM/ENM. On the other hand, currently, SDM/ENMs do not account for the dynamics of a disease and are also unnecessarily restrictive in the set of niche variables considered. Finally, there is, behind the “who”, “where”, and “when”, the “why” that explains them. This, we argue, is the most important role for the use of ENMs in the study of human disease.

In this paper, we show that an ecological approach, using SDM/ENMs, can be usefully applied to any human or non-human disease, transmissible or not, taking as a specific example the spatio-temporal distribution of COVID-19 in Mexico to answer in the affirmative that a respiratory virus does, indeed, have an ecological niche, and that it can be mapped. We show that a wide variety of habitat variables, both environmental, behavioural, and social, and of different types and spatio-temporal resolutions, can be included, to yield a deeper understanding of the factors that drive the COVID-19 pandemic. Moreover, we show how to generalise SDM/ENMs to incorporate and predict the dynamics of the disease. Finally, we show how the formalism can be used to disentangle the complex causal chains that are a fundamental part of a complex adaptive system.

2. Materials and Methods

Several of the methodological elements used in this paper have been used previously for creating ENM/SDMs in multiple contexts [15,16,17,18]. Further details can be found in these papers and in the Supplementary Material.

2.1. Defining a Spatial Grid

All SDM/ENMs are based on the notion of co-occurrence between a target, C, and, one or more, predictors/habitat variables, $\mathbf{X}=(X_1,X_2,\dots ,X_N)$. Normally, co-occurrence is considered purely in spatial terms, although the concept can be extended to co-occurrences in time. In either case, to specify whether there is a co-occurrence or not requires a specification of a grid that divides a spatio-temporal region into cells. A spatial grid can consist of cells of arbitrary shape, as long as they form a partition; i.e., each spatial point is a member of one and only one cell. A partition may be uniform, such as formed by rectangular cells of a given area, or irregular, as is the case for political/administrative units, such as municipalities, counties, states, etc. In standard SDM/ENM modelling [37] this partition is implicit, corresponding, for example, to the pixels of environmental rasters. Such a partition, however, creates a barrier when wishing to include biotic factors, such as point collection data, that cannot be naturally represented as a raster. In short, we may wish to ask: what is the relative importance of average annual temperature, as taken from WorldClim, versus the presence of a prey species in an SDM/ENM for a vagile carnivore?

In Stephens et al. [38], a methodology has been developed for incorporating spatial data layers of arbitrary type and spatial resolution. A spatial grid is overlaid on a chosen spatial region and co-occurrences defined with respect to a given cell. Thus, if there is a co-occurrence between C and a particular variable, $X_i$, in a particular cell, then $N\left(C{X}_i\right)=1$. In the case of a continuous variable, such as species abundance, temperature, or precipitation, the variable, $X_i$, is coarse grained into a set of n discrete bins, leading to n discrete “presence/absence” variables. Thus, $X_i^m\left(\alpha \right)=0,1$ represents the presence/absence in a cell, $\alpha$, of values of $X_i$ in the range defined by the mth bin. Any categorical variable can be left as is, or also coarse grained into a smaller number of categories, if necessary. The criterion for fixing the bin distribution of a variable are that it allows for the best discrimination—dependence of C on $X_i$—while maintaining an appropriate degree of statistical significance—number of cells associated with a given $X_i^m$. Furthermore, the target C can also be discretised if necessary in the same way, with bins $C_i^m$. In summary, all variables become categorical, with each category being associated with a binary variable. An advantage of this categorisation is that no relation is assumed between one bin and another, as would be the case in a regression-based approach for example. Using a non-uniform grid however, can introduce some bias, as some municipalities are bigger than others. We have, however, used spatially uniform grids without substantial changes in our results.

2.2. Constructing SDMs and ENMs

Having transformed all variables to this binomial form, counts, $N\left(C{X}_i^m\right)$, can be made over the region of interest, corresponding to the number of cells that contain a presence of the target $C_i^m$ and the variable $X_i^m$. The number of cells associated with a given $C_i^m$ and/or $X_i^m$ may be fixed or varying depending on the target class. For instance, if we define the class using a relative measure, such as the 10% of cells with the highest abundance for a given species, or the highest number of confirmed cases or deaths in the case of COVID-19, then the class will always be associated with 10% of the cells, independently of time. However, if the class is based on an absolute measure, such as presence/absence of the species, or confirmed cases of a disease, the number of cells in the class may vary in time as, for example, in the case of an invasive species, or an emerging epidemic. Similar considerations hold for $X_i^m\left(t\right)$. If it represents a relative measure, such as the 10% of cells with highest average annual temperature, then it will always cover 10% of cells. However, if it represents the presence/absence of an invasive species, the number of associated cells will change. We will consider these cases in more detail below.

From the co-occurrence counts, probability distributions may be calculated, such as $P\left(C{X}_i^m\right)$, $P(C\mid X_i^m)$, and $P(X_i^m\mid C)$, which are related through Bayes theorem. These distributions can be compared to a null hypothesis and a binomial test used, for instance, to determine the statistical significance of the deviation from this hypothesis. For example, for the posterior, a natural null hypothesis is $P\left(C\right)$ (This is equivalent to the null hypothesis of type SIM2 in the classification of Gotelli [39]. This null hypothesis is one that leads to lower rates of Type I errors and corresponds, in the framework of presence–absence matrices, to keeping the number of observations fixed but randomising their location.) and a statistical diagnostic for determining if the habitat variable $X_i^m$ is significant or not is

$$\epsilon (C\mid X_i^m)=\frac{N_{X_i^m}(P(C\mid X_i^m)-P\left(C\right))}{\sqrt{(}{N}_{X_i^m}P\left(C\right)(1-P\left(C\right)))}$$

(1)

In the case that the binomial distribution can be approximated by a normal distribution, $\mid \epsilon \left(C\right|X_i^m)\mid >1.96$ is equivalent to the 95% confidence interval. For a multivariate niche the corresponding distributions of interest are: $P\left(C\mathbf{X}\right)$, $P\left(C\right|\mathbf{X})$, and $P\left(\mathbf{X}\right|C)$. Although these exist formally from a frequentist perspective, i.e., $P\left(C\right|\mathbf{X})=N(C\mathbf{X})/N(\mathbf{X})$, in the case of a high-dimensional habitat, both $N\left(C\mathbf{X}\right)$ and $N\left(\mathbf{X}\right)=0,1$, which means that a direct statistical estimation is impossible. To overcome this problem, in [40], the likelihood $P\left(\mathbf{X}\right|C)$ has been estimated by assuming a factorisation of the form

$$P(\mathbf{X}\mid C)=\prod _{i=1}^{N_{\xi }}P(\xi \mid C)$$

(2)

where $\xi$ is a combination of a small number of variables. Thus, the abiotic and biotic habitat variables are partitioned into a set of $N_{\xi }$ non-overlapping combinations. In the case that $N_{\xi }=N$, this corresponds to the well known Naive Bayes approximation [41] (Although a complete factorisation of the likelihood may seem a strong assumption, the Naive Bayes approximation has been shown to be a robust performer even in cases where there are strong correlations between variables. An explanation for its surprisingly good performance can be found in [41]). Usually, in order to calculate $P\left(C\right|\mathbf{X})$, as the evidence function $P\left(\mathbf{X}\right)$ is independent of C, to omit it, the following “score” function is often used

$$\begin{array}{ccc}\hfill S(C\mid \mathbf{X})& =& ln(\frac{P(C\mid \mathbf{X})}{P(\overline{C}\mid \mathbf{X})})=ln(\frac{P(\mathbf{X}\mid C)}{P(\mathbf{X}\mid \overline{C})})+ln(\frac{P\left(C\right)}{P\left(\overline{C}\right)})\hfill \\ & =& {\displaystyle \sum _{i=1}^m}s\left(X_i^m\right)+ln(\frac{P\left(C\right)}{P\left(\overline{C}\right)})\hfill \end{array}$$

(3)

where $\overline{C}$ is the set complement of C, with $P\left(\overline{C}\right)=1-P\left(C\right)$, and $s\left(X_i^m\right)=ln(\frac{P(X_i^m\mid C)}{P\left(X_i^m\right|\overline{C})})$ is the contribution (“score”) to the overall $S\left(C\right|\mathbf{X})$ from the variable $X_i^m$. If $s\left(X_i^m\right)>0,<0$ then the factor $X_i^m$ contributes positively/negatively to the occurrence of C. Everything now is based on simple cell counts: $P\left(X_i^m\right|C)=N\left(C{X}_i^m\right)/N\left(C\right)$, $P\left(C\right)=N\left(C\right)/N_s$, $P\left(X_i^m\right)=N\left(X_i^m\right)/N_s$, where $N\left(C\right)$ is the number of cells with a presence of the target C, $N\left(X_i^m\right)$ is the number of cells with a presence of the variable $X_i^m$ and $N_s$ is the total number of cells in the spatial grid. $P\left(C\right|\mathbf{X})$ can be determined directly from $S\left(C\right|\mathbf{X})$ by deriving the relation $P\left(C\right|S\left(C\right|\mathbf{X}\left)\right)=N\left(CS\right)/N\left(S\right)$, discretizing the score range into bins and then counting how many cells in a given score range also have a presence of the target. $\Delta (C,\mathbf{X})=\left(P\right(C\left|\mathbf{X}\right)-P\left(C\right))$ is a measure of how niche-like the conditions specified by $\mathbf{X}$ are. The most niche-like conditions, ${\mathbf{X}}_n$, are those where $P\left(C\right|\mathbf{X})$ reaches its maximum value, while we can term those conditions, ${\mathbf{X}}_{an}$, where $P\left(C\right|\mathbf{X})$ reaches its minimum value, as the most “anti-niche” like. $P\left(C\right|\mathbf{X})$ represents a height function for a given point in an N-dimensional Hutchinsonian ecological niche space. The corresponding “Niche Landscape” is our ENM. As every spatial cell $\alpha$ can be assigned a corresponding niche profile, $\mathbf{X}\left(\alpha \right)$, $P\left(C\right|\mathbf{X}\left(\alpha \right))$ for different spatial cells also now yields an SDM over our spatial region of interest.

A significant advantage of assuming a factorisation of the likelihood is that the subsequent model is completely transparent, with each factor contributing separately, so that each factor can be compared and contrasted with the rest. Moreover, the same applies for groups of factors, so we can compare the relative predictive and discriminative value of climate factors versus socio-economic factors, or air contamination versus mobility.

Of fundamental importance in the construction of $P\left(C\right|\mathbf{X})$ is the statistical ensemble from which the counts will be made. Although we consider a spatial ensemble, the data assigned to each cell invariably have a temporal dimension. In the case of standard SDM/ENM, with a target class defined through point collection data, each data point, $d(\mathbf{x},t)$, is associated with a spatial specification, usually latitude and longitude, and a collection date, t. Similarly, abiotic factors are also time-stamped. A model for $P\left(C\right|\mathbf{X})$ assumes that the distributions in time of both C and $\mathbf{X}$ are both statistically constant and representative.

In the case that the target variable is metric, such as number of cases, $N_C$, the classifier $S\left(C\right|\mathbf{X})$ can also be used to predict $N_C(\alpha ,t)$ for a given spatial cell and over a particular time period. In the case of a spatial model, the cells may be divided up into training and test sets. The points $(N_C\left(\alpha \right),S(C\mid \mathbf{X}\left(\alpha \right)))$ for each cell, $\alpha$, of the training set can then be plotted, and a regression performed to determine the relationship $N_C=F\left(S\left(C\right|\mathbf{X})\right)$. This relation can now be used to predict $N_C$ for any cell in the test set given we have $S(C\mid \mathbf{X})$ there.

2.3. Dynamical ENMs and SDMs

A key aspect of the COVID-19 pandemic is its highly dynamic nature, which can manifest itself in several different ways. First of all, the habitat variables themselves may be time dependent, $X_i\equiv X_i\left(t\right)$. Secondly, the target itself may be dynamical, $C\equiv C\left(t\right)$. In the latter case, this may represent the fact that a disease or a species/disease is present at a given spatial point at one point in time but not another. Thus, the case of an invasive species would fit into this category. In the former, the disappearance of food resources or climate change would naturally fit. In general, we wish to model $P\left(C\left(t\right)\right|\mathbf{X}\left(t^{\prime }\right))$. It is important to point out that there is a difference, however, between having a ENM that is dynamical versus just considering a different configuration of the habitat variables substituted into a given model. In the case of climate change, for instance, this is usually performed [42] by determining a static ENM, equivalent to $P\left(C\right|\mathbf{X}\left(t\right))$, which can be applied to any spatial point, $\alpha$, then using a climate change model to determine $\mathbf{X}(\alpha ,t)\to X(\alpha ,t^{\prime })$ for some future time $t^{\prime }$. In other words, we assume the ENM does not change, only the spatial distribution of the habitat variables, which is then used to determine the new spatial distribution of the target with the same original model. In other words, we determine an SDM at $t^{\prime }$ using an ENM derived at time t.

To move beyond this limited context, we must return to the question of spatio-temporal cells as opposed to just a spatial grid. Considering a timespan T, we divide T into $N_T$ intervals. We now have a total of $N_{Ts}=N_s\times N_T$ spatio-temporal cells. The simplest division is into uniform intervals, such as a year into 12 months. In a static setting, the interpretation of $P\left(C\right|\mathbf{X})$ is that it represents an equilibrium relation between C and $\mathbf{X}$, even though the data used to calculate it, more often than not, spans substantial and, effectively, non-commensurate time periods. For example, using WorldClim data from a given year as habitat variables for a niche model of a species represented by collection data taken from a 150-year time period. Without assuming equilibrium we should calculate $P\left(C\right(t\left)\right|\mathbf{X}\left(t\right))$. However, we may also consider calculating $P\left(C\left(t\right)\right|\mathbf{X}\left(t^{\prime }\right))$, i.e., to predict the effect of the habitat variables $\mathbf{X}\left(t^{\prime }\right)$ at time $t^{\prime }<t$ on our target $C\left(t\right)$. Of course, an important question at the heart of this is: how are changes in $\mathbf{X}\left(t^{\prime }\right)$ reflected in $C\left(t\right)$? This requires longitudinal observations and/or an understanding of the underlying causal relationships between the $X_i$ and C.

There are four distinct paths to incorporating time into the SDM/ENM: (i) assume equilibrium, and thereby ignore time dependence; (ii) construct a model $P\left(C\right(t\left)\right|\mathbf{X}\left(t\right))$ using a time slice/history to predict a spatial distribution $C(\alpha ,t)$ on that time slice, as a function of the habitat variables on the same time slice; (iii) predict in time, assuming niche conservatism, i.e., construct $P\left(C\right(t\left)\right|\mathbf{X}\left(t\right))$ and use it to predict the spatial distribution, $C(\alpha ,t^{\prime })$ at some later time using the habitat $\mathbf{X}(\alpha ,t^{\prime })$; (iv) predict in time, without assuming niche conservatism, by constructing $P\left(C\left(t\right)\right|\mathbf{X}\left(t^{\prime }\right))$, where $t^{\prime }<t$. Models of type (i) and (ii), we term spatial prediction models. They are analogous to standard SDM/ENMs, in that they predict a spatial distribution. However, in distinction, in case (ii) they do so using a time slice of data that permits us to compare and contrast $C\left(t\right)$, $\mathbf{X}\left(t\right)$ and $P\left(C\right(t\left)\right|\mathbf{X}\left(t\right))$ over time. For instance, we may use data only from May 2021 and compare with data from May 2020, or we may consider data from all of 2021, etc. Models of type (iii) and (iv) we term time prediction models. In this case they are an ENM/SDM equivalent to a SIR-type model, predicting the evolution of “where” in time, as opposed to “who”. The case with niche conservatism would be equivalent to a SIR type model, where the parameters of the model do not change, whereas the non-conservative case would correspond to the case where they do change and are fitted to dynamic data.

The notion of a time prediction model also naturally leads us to consider ENM/SDMs that are associated with changes in the distribution of a species/disease. For instance, we may consider the set of cells, $C(t-1)$, that have a presence in the time interval $t-1$ and the set of cells, $C\left(t\right)$, that have a presence in the time interval t. ${\Delta }_C(t,t-1)$ may then represent those cells that had a presence/absence at $t-1$ and, in contrast, an absence/presence at t. This could, for example, model the range expansion of a species. In this case, an ENM $P\left({\Delta }_C(t,t-1)\right|\mathbf{X}(t-1))$ represents a model for predicting changes in the distribution due to the habitat variables. This model can then be applied to produce an SDM $P\left({\Delta }_C(\alpha ,t+1,t)\right|\mathbf{X}(\alpha ,t))$ using the habitat variables at t and predicting the change in distribution between t and $t+1$.

We give more details about how to construct the different spatio-temporal models in the Supplementary Material. As discussed above, spatially, we used a grid corresponding to the municipalities of Mexico. For the temporal partition, we chose a month, though any other timescale of interest could have been used.

2.4. Testing Model Performance

The way in which our models are tested depends on the type of statistical ensemble that is used to train and then test the model. In the case of an ENM on a given time slice, $P\left(C\right(t\left)\right|\mathbf{X}\left(t\right)$, the model will be created on a training set that corresponds to a randomly chosen fraction, $f_{train}$, of spatial cells with data associated with a time period t. The model is then tested on the remaining fraction of cells $f_{test}=1-f_{train}$ from the same time period. The test and train sets can be chosen in multiple ways. Here, we consider a 70%/30% split. Standard performance statistics can then be determined, such as from a confusion matrix or from the area under the ROC, among others [43,44]. However, we can also use the ENM $P\left(C\right(t\left)\right|\mathbf{X}\left(t\right))$ trained on 100% of spatial cells and apply it to a future time period $t^{\prime }$. In this case, the entire set of spatial cells at $t^{\prime }$ is the test set. The luxury of a problem such as COVID-19 is that we have relatively accurate data about the spatial distribution of the disease as a function of time. In this case, the ENM $P\left(C\right(t\left)\right|\mathbf{X}\left(t\right))$ is used by substituting at $t^{\prime }$, $\mathbf{X}(\alpha ,t^{\prime })$ for each spatial cell, $\alpha$, to obtain for that cell $P\left(C(\alpha ,t^{\prime })\right|\mathbf{X}(\alpha ,t))$. We can now test the quality of the prediction using any of the above standard metrics, given that $P\left(C(\alpha ,t)\right|\mathbf{X}(\alpha ,t^{\prime }))$ is a classifier, and we can compare with the actual distribution $C(\alpha ,t^{\prime })$. We will expect the time t ENM to yield good performance at $t^{\prime }$ only if the niche is conserved over this time period.

Aside from the pure assumption of niche conservatism, where we apply the ENM $P\left(C\right(t\left)\right|\mathbf{X}\left(t\right))$, we can compare and contrast ENMs from different time slices to determine the degree to which they are changing. This would correspond to determining if the niche is actually changing over time. This can be performed for each habitat variable by comparing $s\left(X_i^m\left(t\right)\right)$ with $s\left(X_i^m\left(t^{\prime }\right)\right)$. If the niche is changing then we must consider just how fast it is changing. This can be deduced by comparing the $s\left(X_i^m\left(t\right)\right)$ across time. If $s\left(X_i^m\left(t\right)\right)\sim s\left(X_i^m\left(t^{\prime }\right)\right)$ for two time periods t and $t^{\prime }>t$, then the niche is conserved over the interval $(t^{\prime }-t)$. If they differ substantially, then we may reduce $(t^{\prime }-t)$ and determine how much change there has been over this shorter interval.

For time prediction models that do not assume niche conservatism, the training set consists of choosing two time periods $t-1$ and t. An ENM $P\left(C\right(t\left)\right|\mathbf{X}(t-1))$ is created on all spatial cells. This model is then applied to all spatial cells for the time periods t and $t+1$ as test set. Thus, we use the ENM $P\left(C\right(t\left)\right|\mathbf{X}(t-1))$, substituting $\mathbf{X}(\alpha ,t)$ to create the SDM $P\left(C\right(\alpha ,t+1\left)\right|\mathbf{X}(\alpha ,t))$ for prediction of $C(\alpha ,t+1)$. We can then apply any of the above standard metrics to evaluate the performance of this ENM/SDM.

2.5. Predicting Abundances

ENMs that are based on presence/no presence are, naturally, binary classification models. Although a metric variable, such as abundance, can be treated as such by considering multiple classes, it is also possible to use a binary prediction model, such as the top 10% of municipalities with highest confirmed cases, to construct a relation, $N\left(S\right(C\left|\mathbf{X}\right))$ between the score, $S\left(C\right|\mathbf{X})$, and the number of cases on a training set. This model can then be applied to a test set, predicting the expected number of cases for a given cell, $\alpha$, using its habitat profile $\mathbf{X}\left(\alpha \right)$. This procedure can be used for both spatial and time prediction models. In the former case, the abundance predictions are associated with a given time slice t using a split of the cells into training and test sets, while for the latter we predict from a training set that consists of all cells on a time slice t to predict abundances on a test set of all cells on a time slice $t^{\prime }$.

2.6. Inferring Causality in ENMs

Another criticism of the application of standard SDM/ENMs to the pandemia has been the lack of plausibility of the relation between, say, infection rates and habitat variables, such as climate or contamination. As pointed out in [40], this is a problem with correlative approaches in general. This criticism however, can be applied to a phenomenological study of any complex adaptive system, where cause–effect relations are multi-layered and, therefore, often indirect. Epidemiology and medicine in general are rife with problems of causal inference and there are two basic approaches to it: a “classic” approach [45] and a “modern” approach [46,47,48]. An important criterion from our perspective is that of “strength of association” [45], where, although a small association does not mean that there is no causal effect, the larger the association, the more likely that it is causal. This is key to our use and understanding of both $\epsilon \left(C\right|X_i^m)$ and $s\left(X_i^m\right)$ and their multi-factorial counterparts $\epsilon \left(C\right|\mathbf{X})$ and $s\left(\mathbf{X}\right)$. Thus, for two niche factors, $X_{\alpha }$ and $X_{\beta }$, if $s\left(X_i^m\right)>s\left(X_j^n\right)$, then we will judge $X_i^m$ to be causally closer to C than $X_j^n$. An illustrative example of this would be a food chain: carnivore → herbivore → plant food → climate, as considered in [40].

In [40,49], we have proposed a formalism for examining causality that is particularly natural for application to ENMs. If we have two niche factors $X_i^m$ and $X_j^n$, we can better understand their potential causal relations with a target C by considering the following relations

$$\epsilon \left(C\right|X_i^m{X}_j^n;P_N)=\frac{N_{X_i^m{X}_j^n}(P\left(C\right|X_i^m{X}_j^n)-P_N)}{\sqrt{(}{N}_{X_i^m{X}_j^n}P\left(C\right|X_i^m)(1-P\left(C\right|X_i^m)))}$$

(4)

where $P_N$ represents the null hypothesis with respect to which we will determine the predictability of the combination $X_i^m{X}_j^n$. If $P_N=P\left(C\right)$ we are gauging the consistency of $P\left(C\right|X_i^m{X}_j^n)$ with the null hypothesis that the distribution of the target species is independent of the variable combination $X_i^m{X}_j^n$. However, we may also choose as null hypotheses $P_N=P\left(C\right|X_i^m)$ or $P_N=P\left(C\right|X_j^n)$, in which case the null hypothesis is that $P\left(C\right|X_i^m{X}_j^n)$ is independent of $X_j^n$ and $X_i^m$, respectively. With this approach, we can determine the degree to which $X_j^n$ confounds $X_i^m$ or vice versa. For example, if $P_N=P\left(C\right|X_i^m)$ and $\epsilon \left(C\right|X_i^m{X}_j^n;X_i^m)>1.96$, we can conclude that within a 95% confidence interval that $P\left(C\right|X_i^m{X}_j^n)$ is not consistent with the null hypothesis and therefore the habitat variable $X_j^n$ is predictive of the distribution of C beyond what is explained by the habitat variable $X_i^m$. As $X_i^m$ may be a biotic variable and $X_j^n$ an abiotic one, we have used this to show that, very often, biotic factors are confounders for abiotic factors and not vice versa [40]. Although, here, we are concentrating on causal relations and confounding with respect to pairs of habitat variables the formalism naturally extends to larger numbers of variables.

In any spatial cell $\alpha$ we may determine if there is a presence or absence of a given coarse grained bin for either variable, $X_i^m\left(\alpha \right)=0,1$ and $X_j^n\left(\alpha \right)=0,1$. For example, $X_8^4\left(\alpha \right)X_{10}^2\left(\alpha \right)=11$ would represent a presence of bin 4 of variable 8 and a presence of bin 2 of variable 10 in the cell $\alpha$, while $X_8^4\left(\alpha \right)X_{10}^2\left(\alpha \right)=01$ would represent an absence of bin 4 of variable 8 and a presence of bin 2 of variable 10 in the cell. Thus, there are 4 possible combinations for any pair of habitat variables: presence–presence = 11, presence–absence = 10, absence–presence = 01, absence–absence = 00. By comparing and contrasting the different combinations, we may determine, for example, if presence of one habitat variable is more predictive than the other. Note $X_j^n\left(\alpha \right)=0$ does not imply absence of the variable itself in the cell, just the range denoted by the nth coarse grained bin.

2.7. Data and Habitat Variables

In this paper, considering the evolution of the pandemia in Mexico, we used an irregular spatial partition consisting of the 2458 municipalities of Mexico. The advantage of this partition is that it is most aligned with publicly available socio-demographic and socio-economic factors that can serve as corresponding niche variables. Its disadvantage is that there is a potential bias in the spatial distribution of different municipalities according to their area. For the epidemiological data, used to define the target classes, we used data from the General Directorate of Epidemiology, Secretariat of Health in Mexico [50]. For the socio-demographic and socio-economic data, we used 124 variables taken from the 2010 Mexican Census [51] at the municipality level. For mobility data, we used 12 variables provided by the Institute of Geography of UNAM, that represent the average, daily labour flows between a pair of municipalities. For the air contamination factors, we used three atmospheric pollutants (formaldehyde (HCHO) nitrogen dioxide (NO${}_2$), sulphur dioxide (SO${}_2$)) [52,53,54], while for climactic data, we used 19 bioclimatic variables from the WorldClim database (www.worldclim.org (accessed on March 2022)) with a spatial resolution of 30 arc-seconds (≈1 km) [55], which includes 11 temperature and eight precipitation indices that express annual trends (e.g., annual mean temperature and precipitation), seasonality (e.g., annual temperature and precipitation ranges), and environmental extremes (e.g., highest and lowest values of temperature for the warmest and coolest months).

As all the socio-demographic, socio-economic, and mobility data were metric and continuous, for each variable, we ranked the 2458 municipalities from highest to lowest value then divided the ranked list into deciles, with 10% of all municipalities in each decile. By doing this, as opposed to dividing the metric interval for the variable into equal parts, we assure equal statistical weight to each coarse grained category. For the air pollution variables, we divided each raster layer into 20 ranges, which were chosen to have roughly equal number of pixels in each one. Climate variables were also divided into 20 ranges. Thus, in this way we mapped both target class and habitat variables in their entirety into a set of binomial presence/absence variables for each cell on our grid. We also included in as a habitat variable the decile of confirmed cases associated with a given municipality, but at period $t-1$, as opposed to the target variable which was associated with period t. In this way, we could show, in principle, how the history of the habitat can also be included as a predictor. Indeed, this is the first step at showing how the present methodology may be developed to include SIRs-type modelling properties.

3. Results

An important aspect of a complex, adaptive phenomenon, such as the COVID-19 pandemic, is that there are many questions of “where”, “when”, and “why” that require answers, with each question in turn requiring its own SDM/ENM. Here, we will present representative results for the following: (i) Where are the highest number of confirmed cases in a given time period to be found? (ii) Where are confirmed cases in a given time period to be found? (iii) Where will most confirmed cases be found in a future time period? (iv) Where will cases be found in a future time period where previously there were none? In a more ecological language and taking confirmed cases of COVID-19 as the “species” of interest: (i) Where is the highest abundance of the species of interest to be found? (ii) Where is the species to be found in a given time period independently of its abundance? (iii) Where will the species be found in highest abundance in a future time period? (iv) Where will the species be found where it was not present previously? (In the online system, EpI-PUMA, publicly available in a Platform-as-a-Service environment ((http://covid19.c3.unam.mx), accessed on 5 April 2021, the above and many more different SDM/ENMs are available that use the methodology described in this paper). These five representative models illustrate important differences about how ENM/SDMs can be used to answer where, when, and why questions. Firstly, we may use a model developed on a region $R_{train}\left(t\right)$, using data from a time period t, to predict on another region $R_{test}\left(t\right)$, as is standardly performed, and where the complete spatial region under consideration is $R\left(t\right)=R_{train}\left(t\right)\cup R_{test}\left(t\right)$. Secondly, we may use a model developed on a region $R\left(t\right)$ to predict on $R\left(t^{\prime }\right)$, the same region at some later time.

3.1. Dynamic Biotic and Abiotic Factors Can Be Included in an SDM/ENM for COVID-19

To answer each of the above questions, we must extend standard ENM/SDM modelling to include both dynamical target classes, $C\left(t\right)$, and associated dynamical habitat variables, $\mathbf{X}\left(t^{\prime }\right)$, and develop SDM/ENMs, $P(C\left(t\right)\mid \mathbf{X}\left(t^{\prime }\right))$, that relate them. With question (i) the target class will be the 10% of Mexican municipalities that have the highest number of cases in month t. In that context we are characterising niche by the relative abundance of the target species, with those locations where it is highest being characterised as most niche-like. This target class is time dependent, as it may be that those places with the highest abundance at one time are not the same as those at another. With question (ii) the target class will be the presence of confirmed cases in month t. Questions (iii) and (iv) refer to these target classes but in a future month $t^{\prime }$. In particular, for the presence of confirmed cases, we will be interested in predicting those municipalities that had no previous cases in time period t but do in period $t^{\prime }$.

Thus, the targets for our ENM/SDMs are all, in principle, time-dependent. For the habitat variables, we use those presented in the Methods section below, with an illustrative dynamic habitat variable being the decile of municipalities corresponding to a given range of number of confirmed cases in the period $t-1$.

3.2. Predictive SDM/ENMs Can Be Created for COVID-19

SDM/ENMs were created for the target classes (i)–(v). In Figure 1 (left) we see the results for a purely spatial ENM/SDM model with target class, corresponding to question (i), being the 10% of municipalities with the highest number of confirmed cases, considered for three different, representative months of the pandemic: March 2020, June 2020, and January 2021. In this case, performance is based on a 5-fold 70%/30% train/test split of the 2458 municipalities. Figure 1 (left) shows the corresponding ROC curves and the corresponding AUC for five separate SDM/ENMs for each test set, with: socio-demographic/economic variables only, mobility only, climate only, air contamination only, and all factors together, for each of the three considered months. The relative performance of each sub-model is clear, with the mobility and socio-demographic models being the best performers, followed by the air contamination model, then the climate model, with the total model equivalent to the mobility and socio-demographic models.

Table 1. Top 18 most niche-like and bottom 18 most anti-niche-like habitat factors for the ENM for top 10% of municipalities with highest number of cases ranked by score in March 2020 and showing the corresponding rank in June 2020 and January 2021.

Percentile	Type	Variable	Rank (March)	Rank (June)	Rank (January)
Top	Mobility	10-Internal-labour-flow 25719:664418	1	1	1
Top	Mobility	10-Labour-flow-inward 3279:508685	2	2	2
Top	Mobility	10-Inward-centrality 0.026:0.347	3	3	5
Top	Mobility	10-Inwards-municipalities-connected 64:849	4	4	6
Top	Socio-Demographics	10-Inhabited private homes that have Internet	5	7	3
Top	Socio-Demographics	10-Inhabited private homes that have a computer, tablet or laptop	6	8	4
Top	Mobility	10-Labour-flow-out 4954:287854	7	9	7
Top	Socio-Demographics	10-Total-centrality 0.053:0.416	8	5	10
Top	Socio-Demographics	10-Total municipalities-connected 130:1018	9	6	11
Top	Socio-Demographics	10-Male population aged 18 and over with post-basic education	10	11	9
Top	Socio-Demographics	10-Population aged 18 and over with post-basic education	11	10	8
Top	Socio-Demographics	10-Inhabited private dwellings that have a cell phone	12	12	12
Top	Socio-Demographics	10-Female population aged 18 and over with post-basic education	13	13	13
Top	Socio-Demographics	01-Illiterate population aged 15 and over	14	15	16
Top	Socio-Demographics	01-Illiterate male population aged 15 and over	15	17	19
Top	Socio-Demographics	01-Population aged 15 and over without schooling	16	23	25
Top	Socio-Demographics	10-Population from 15 to 64 years old	17	14	22
Top	Socio-Demographics	10-Inhabited private dwellings that have a television	18	18	14
Bottom	Mobility	01-Total municipalities-connected 2:14	1773	1702	1699
Bottom	Mobility	02-Labour-flow-inward 22:51	1774	1725	1726
Bottom	Socio-Demographics	10-Population with limitation to walk, go up or down	1775	1726	1727
Bottom	Socio-Demographics	10-Population born in the entity	1776	1727	1728
Bottom	Socio-Demographics	02-Inhabited private homes that have a washing machine	1777	1728	1729
Bottom	Socio-Demographics	01-Inhabited private dwellings that have drainage	1778	1754	1757
Bottom	Socio-Demographics	10-Inhabited private dwellings that do not have drainage	1779	1755	1758
Bottom	Socio-Demographics	01-Inhabited private dwellings that have a refrigerator	1780	1756	1759
Bottom	Socio-Demographics	01-Population from 6 to 11 years old that does not attend school	1781	1610	1587
Bottom	Socio-Demographics	07-Population with limitation	1782	1773	1774
Bottom	Mobility	03-Internal-labour-flow 698:1179	1783	1783	1783
Bottom	Mobility	01-Internal-labour-flow 9:329	1784	1784	1784
Bottom	Mobility	01-Labour-flow-out 1:20	1785	1785	1785
Bottom	Mobility	01-Labour-flow-inward 1:21	1786	1786	1786
Bottom	Mobility	01-Outwards-centrality 0:0.003	1787	1789	1787
Bottom	Mobility	01-Outwards-municipalities-connected 1:7	1788	1790	1788
Bottom	Mobility	01-Inward-centrality 0:0.002	1789	1667	1650
Bottom	Mobility	01-Inwards-municipalities-connected 1:5	1790	1668	1651

shows the most niche-like and most anti-niche-like habitat factors, as ranked by the statistical measure of co-occurrence $\epsilon$ (see Section 2) (see the Methods section) for each month. As can be seen, the top 18 and bottom 18 ranked factors are almost identical for the three different periods.

In Figure 2 (bottom), we show the values of the score $s\left(X_i^m\left(t\right)\right)$ (see Section 2) for the 10 deciles, $m=[1,10]$, of the habitat variable $X_i=$ Internal labour flow of the municipality, observing that the score contributions are highly conserved across the three time periods. This niche conservation effect is shown in an even more striking fashion in Figure 3, where we see the correlation between the $\epsilon \left(X_i^m\left(t\right)\right)$ values (left) and the score values, $s\left(X_i^m\left(t\right)\right)$, (right) for the ENM models for each individual habitat variable for the three different periods. For instance, we observe that the $\epsilon \left(X_i^m\left(t\right)\right)$ values from March 2020 explain 96% of the variance in the $\epsilon \left(X_i^m\left(t\right)\right)$ values for January 2021. As a further test of niche conservation, we used an ENM trained on all of the cells in one month to predict a future month. In Figure 1 (right), we see the ROC and AUC for the five models using different categories of habitat variable (mobility, socio-demographic, air contamination, climate, and all) trained on data from all spatial cells for the period March 2020 and tested on the periods June 2020 and January 2021. Similarly, we show the result of analogous models trained on the period June 2020 and tested on the period January 2021.

In Figure 4 (left), we see the results for a purely spatial ENM/SDM model with target class, corresponding to question ii), being those municipalities with the presence of confirmed cases of COVID-19 for the months: March 2020, June 2020 and January 2021. Performance is again based on a 70%/30% train/test split of the 2458 municipalities, with Figure 4 (left) showing the corresponding ROC curves and the corresponding AUC for five separate SDM/ENMs for each test set, with: socio-demographic/economic variables only, mobility only, climate only, air contamination only, and all factors together, for each of the three considered months. Again, the relative performance of each sub-model is clear, with the same relative performance as for the previous models. However, we note that for the overall performance of the mobility and socio-demographics of all models is less than their counterparts in the case of top 10% of municipalities with most cases as class, while the performance of the climate and air contamination models is similar.

In Figure 5 (left), in analogy with Figure 3, we see the correlations between $\epsilon \left(X_i^m\left(t\right)\right)$ values for all habitat variables for the model for the presence of COVID-19 cases for two different months, while in Figure 5 (right), we show the correlations between $s\left(X_i^m\left(t\right)\right)$ values for all habitat variables for the same model. Note that the correlation between the $\epsilon \left(X_i^m\left(t\right)\right)$ and $s\left(X_i^m\left(t\right)\right)$ distributions in March 2020 and June 2020 and March 2020 and January 2021 are now significantly smaller than their counterparts for the top 10% model.

Furthermore, in Figure 2 (Top) we show the values of the score $s\left(X_i^m\left(t\right)\right)$ for the 10 deciles, $m=[1,10]$, of the habitat variable $X_i=$ Internal labour flow of the municipality to compare and contrast the score contributions of this key niche variable as a function of time. We observe that the relative contributions from each decile differ significantly when compared to the results for The top 10% of municipalities with highest number of cases, as shown in Figure 2 (Bottom).

Turning now to the above ENMs at time t as predictors of the species distribution at $t^{\prime }$, in Figure 4 (Right), we see the results of a model trained on data from those municipalities with confirmed cases in March 2020 and June 2020 to predict which municipalities would have cases in June 2020 and January 2021 and January 2021, respectively. In this case, we observe that the ENM for month t is less predictive for month $t^{\prime }$ than the corresponding model for identifying the top 10% of municipalities with the highest number of cases.

3.3. ENMs Can Be Used to Predict Numbers of Cases

As well as using $S\left(C\right|\mathbf{X})$ as a pure classifier (see Section 2), to identify a particular class, such as presence/absence, top 10% highest cases, etc., we can also use it to predict the actual number of cases. To perform this, we regress the total score $S\left(C\right|\mathbf{X})$ against the number of cases on our training set and then apply it to the test set. An example is shown in Figure 6, for the classification model with $C=$ top 10% of municipalities with the highest number of cases, where the training set is a randomly chosen 70% of spatial cells and the test set the remaining 30%. All habitat variables were used. Figure 6 (left) shows the relation between score and number of confirmed cases for our three example months—March 2020, June 2020, and January 2021. First, municipalities were ranked according to their score, then grouped into 10 deciles. An exponential function fit was used. Figure 6 (right) shows the results of applying the model to the 30% out of sample test set, showing the relation between predicted and actual abundances.

3.4. Causal Relationships Can Be Deduced

As an illustration of the formalism for disentangling causal relationships (see Section 2), we consider the relative contributions of climatic factors, as represented by the WorldClim variable Average annual temperature, and human mobility factors, as represented by Internal labour flow of the municipality, to prediction of the top 10% of municipalities with the greatest number of cases. The temperature variable is divided into 20 coarse-grained bins, while the mobility variable is divided into 10 bins. In any spatial cell, $\alpha$, we may then determine if there is a presence or absence of a given coarse grained bin for either variable. Thus, in Equation (4), we have $X_i^m$ representing the 10 different bins of the mobility variable and $X_j^n$ representing the 20 bins of the temperature variable. As, for a given spatial cell, there are four possibles states, corresponding to presence/absence of the two habitat variables, $X_i^m{X}_j^n$ can be represented by four $20\times 10$ matrices.

For each cell of each matrix, we calculate $P\left(C\right|X_i{X}_j)$ and also $\epsilon \left(C\right|X_i^m{X}_j^n)$, using as null hypothesis $P\left(C\right)$, as shown in Figure 7, with data from January 2021. Note that the absence of a variable bin, $X_i^m$, corresponds to those cells where there is no presence of $X_i^m$. However, there will be a presence of at least one other $X_i^{m^{\prime }}$, for $m^{\prime }\ne m$. For example, in the matrix $Clim0\_Mob0$, the cell $Tmp.p04=0,08\_Inter-Mob=0$ corresponds to those cells that are not in the average annual temperature range corresponding to $Tmp.p04$ and also not in the Internal labour flow range corresponding to $08\_Inter-Mob$. However, that set of cells could have presences of temperature ranges other than $Tmp.p04$ and mobility ranges other than $08\_Inter-Mob$.

The fact that the climate variable is confounded by the mobility variable is manifest in the matrix $Clim1\_Mob0$, where there is a presence of a particular climate variable bin and an absence of the corresponding mobility variable bin. For instance, the cell $Tmp.p04=1,10\_Inter-Mob=0$ in that matrix corresponds to those municipalities that are not in the highest decile of mobility, but are in the fourth decile of average annual temperature. $\epsilon \left(C\right|X_i^4{X}_j^{10})=-8.05$ and $P\left(C\right|X_i^4{X}_j^{10})\sim 0$ there, indicating that the probability of being in the top 10% of highest number of cases is very low. Indeed, we can see that for any temperature range, for $10\_Inter-Mob=0$, there is very little chance of the corresponding municipalities being in the top 10%. On the contrary, for $Clim0\_Mob1$, we see that $P\left(C\right|X_i^m{X}_j^{10})\sim 0.8$ for any m, thus indicating that high mobility is very niche-like, independently of the value of the temperature variable. Although for $Clim1\_Mob0$ the range $Tmp.p04$ can give rise to statistically significant $\epsilon \left(C\right|X_i^4{X}_j^n)$ values for $m\ne 10$, by examining the case $Clim1\_Mob1$, we see that the $\epsilon \left(C\right|X_i^4{X}_j^n)$ values for $n\ne 9,10$ are all negative and statistically significant. Thus, climate may appear to be niche-like but, in fact, is confounded by mobility, and now we can intuit why a climate model can create predictability even though climate is not predictive in and of itself.

A model based on climate only is equivalent to the combination $X_i^m=1$, $X_j^n=\ast$, where ∗ denotes that we have marginalized over this value. For instance, $P\left(C\right|X_i^m)={\sum }_{X_j^n=0,1}P\left(C\right|X_i^m{X}_j^n)=P\left(C\right|X_i^m{X}_j^n=0)+P(C\mid X_i^m{X}_j^n=1)$. Although, $P\left(C\right|X_i^m{X}_j^n=0)\sim 0$ for $n\ne 9,10$, $P\left(C\right|X_i^m{X}_j^{10}=1)\sim 0.6-0.8$ for any m. Thus, the primary reason why the climate model has some degree of predictability is that those places of highest mobility have a climate “profile”, i.e., they are not equally distributed across all average annual temperature ranges. It is not, however, the climate that is causal. Another way to see this is shown in Figure 8, where we show the marginal probabilities $P\left(C\right|X_i^m)$ and $P\left(C\right|X_j^n)$, which would correspond to the results associated with an ENM based only on the climate variable Average annual temperature or on the human mobility variable Internal labour flow of the municipality, respectively.

We can note that, although the degree of predictability in the climate-only model is weak, there is a variation, with $\epsilon \left(C\right|X_i^m=1)$ ranging from $0.9$ to $-2.9$. However, the origin of the most niche-like value of $0.9$ is associated with the residual predictability from decile 10 of the mobility variable.

4. Discussion

Our goal in this article was to demonstrate that the questions: “Does a respiratory virus have an ecological niche, and if so, can it be mapped?” can be answered in the affirmative. We have explicitly created several ENMs and SDMs for COVID-19 that are both predictive and contain habitat factors that are more causally plausible than climate, for instance. In order to achieve this, we introduced several innovations compared to standard niche and species distribution modelling. Firstly, we showed how to extend niche and species distribution modelling to “non-equilibrium” situations, where both target and niche variables are potentially time varying, as well as the relation between them. Secondly, we created models with habitat variables that were represented by quite different data types and associated spatial resolutions. Finally, we showed how causal relations and confounding can be better understood by introducing a hierarchy of conditional probabilities and the associated intuition that a more causally direct factor should have a bigger effect than an indirect one.

We constructed ENMs as Niche Landscapes, $P\left(C\right(t\left)\right|\mathbf{X}\left(t\right))$ - Bayesian posteriors which serve as a height function on a Hutchinsonian ecological space with 2749 dimensions, spanning air pollution, climate, mobility, socio-economic, and socio-demographic data. Usually in ENM/SDMs, the target class, C, is a binary variable, such as presence/no presence of a species. In Coro [5], the corresponding variable was “high” infection rate, defined in a binary fashion with respect to a reference infection rate. We too, have used binary class variables, by choosing a particular subgroup of spatial cells, corresponding to presence/absence of confirmed cases, or if a spatial cell was in the top 10% of highest total infections (In the EpI-PUMA system, publicly available in a Platform-as-a-Service environment (http://covid19.c3.unam.mx, accessed on 5 April 2021), 72 different SDM/ENMs are available that use the methodology described in this paper.). However, the binary nature of the target variable is a choice rather than a restriction. For instance, taking infection rate as a continuous variable, this can be divided into as many quantiles as we please, say $n_c$. The target variable now consists of $n_c$ “presence/absence” variables, each with its own ENM/SDM, $P\left(C_i\right|\mathbf{X})$, $i\in [1,n_c]$. Even in the case of a binary decomposition of a continuous variable; however, the metric nature of the variable leaves an imprint on $P\left(C\right|S\left(C\right|\mathbf{X}\left)\right)$, such that higher score, as a continuous metric variable, corresponds to higher infection rate or number of cases, as we have demonstrated, leading to a model that can predict abundances.

In criticising ENM/SDMs as a useful tool for the COVID-19 pandemic, or any other, it is important to distinguish between applicability of ENM/SDMs in general versus a particular instance of an ENM/SDM. It is appropriate to criticise a model that includes climate, and which has been used to infer a corresponding causal effect on the pandemic, without any analysis of possible confounders. This does not mean, however, that with appropriate habitat variables and a methodology for disentangling confounding, that useful ENM/SDMs cannot be constructed, as we have shown here. Our results clearly show that models for predicting where the highest number of cases will be, that are built on mobility, socio-demographic and socio-economic data, are much more predictive than models built on climate and/or air contamination, as can be seen in the results of Figure 1, and that this has been true throughout the pandemic. As seen in Table 1, the most important niche/anti-niche factors for this target are all associated with the highest/lowest levels of mobility, as proxied by inter- and intra-municipal labour flows, and a particular socio-economic profile. The fact that a climate model can exhibit some degree of predictability does not mean that climate is the direct driver. On the contrary, we can ask if any apparent predictability due to climate factors is confounded by human-based factors. As we have shown in Figure 7 and Figure 8, this is indeed the case. Just the distribution of the habitat variable scores themselves tells us that if there is confounding then it is the human factors that confound climate and not vice versa, as a confounder should be causally closer to its effect than the confounded variable and therefore have a higher score, as is implicit in the Bradford–Hill criteria.

With respect to the socio-demographic and socio-economic factors, we see that they have a similar predictive power to the mobility factors and reflect the socio-demographic and socio-economic conditions where COVID-19 cases are highest. For instance, the habitat variables % of households that have a computer and % of households that have internet access are both significant niche factors, as are other factors that correspond to a more educated population with higher economic status. Of course, the interpretation of these factors is not as intuitive as mobility and we certainly do not wish to attribute direct causality to them. However, there are a variety of relevant factors for COVID-19 that can be related to internet access and computer usage for example, such as age, educational status, population density [56], as urban areas have better infrastructure, and mobility itself [57,58]. As with any epidemiological or ecological model, the interpretation of the predictors as representing direct versus indirect interactions is highly non-trivial. The formalism we have introduced for disentangling confounding can help in this regard.

The places where COVID-19 is in highest abundance represent one particular characterisation of the spatio-temporal distribution of the COVID-19 “species” and its associated ecological niche. Where it is present and where it will be in the future compared to where it is now are two others. Hutchinson [7,59] defined niche as that region of ecological space where the net growth rate of the species is $r\ge 0$ at low density. For a pathogen that is not capable of free movement and is dependent on a host, there are two natural growth rates: one associated with the pathogenic load within a host and another that is measured by the number of infected hosts. Obviously, in this paper we are concerned with the latter, where the growth rate is characterised by the basic reproduction rate, $R_0$, or the effective reproduction rate, $R_t$ [60]. If we defined niche for COVID-19 through an analogue of $r\ge 0$, such as $R_t>1$, we would clearly be in a situation where we were passing from “niche” to “anti-niche” ($R_t>1\to R_t<1$) and vice versa, continuously in space and time due to a multitude of factors, including public health interventions, such as lockdowns and vaccinations, as well as resistance to vaccinations, new variants, and a host of others. That these factors alter the ecological conditions in a certain place at a certain time is undeniable. However, to keep track of such changes and how they impinge on how niche-like conditions are in a certain spatio-temporal cell, $(\alpha ,t)$, requires the corresponding data. Here, we have preferred to use characteristics of the pathogen distribution that are easier to measure—number or presence of cases—but with which we can characterise the concept of niche, and a set of habitat variables that go beyond those previously considered. We take highest case abundance to distinguish those conditions in ecological space, and in geographical space and time, that favour higher abundance of the pathogen, as proxied by abundance of cases. This is a relative measure, as it may be that a municipality, $\alpha$, has $R_t\left(\alpha \right)<1$ but such that it is higher than any other.

Holt [61] has suggested that besides an “Establishment Niche”, that corresponds to the original Hutchinsonian niche, a “Population Persistence Niche”, associated with the range of niche space in which populations above some threshold density, $N>0$, can persist, is a complementary notion. We take where COVID-19 is in highest abundance to be closer to this Population Persistence Niche, whereas where it is present is closer to the original Hutchinsonian characterisation and, especially, as it is portrayed in the majority of ENM/SDMs [37], where presence/absence is used to characterise both. On the other hand, where it is now versus in the past corresponds to neither, but is taken to reflect the potential range expansion/contraction of the species. We believe that these examples, and more, show that there are multiple characteristics of a species distribution that can and should be modelled and can be used to characterise complementary notions of niche.

The methods we have exhibited and the corresponding examples also allow us to understand to what degree a niche is conserved. In the case of species abundance, we saw that the niche associated with those places where the relative number of COVID-19 cases was highest was highly conserved, with very little difference across the entire pandemic. Moreover, we showed how this conservatism could be quantified using our statistical diagnostic $\epsilon (C\mid X_i^m)$ and the score functions, $s\left(X_i^m\left(t\right)\right)$, associated with the habitat variables $X_i^m$, with the time dependence of the associated score function reflecting the relationship between target and habitat variables. For example, if $s\left(X_i^m\left(t\right)\right)$ is strongly positive at one time, t, and not another, $t^{\prime }$, then we may say that the variable $X_i^m$ is niche-like with respect to C at time t, but not at $t^{\prime }$. Niche conservatism with respect to $X_i^m$ is then quantified by $s\left(X_i^m\left(t\right)\right)\sim s\left(X_i^m\left(t^{\prime }\right)\right)$. This niche conservatism is also manifest in that an ENM created at time t provides a SDM that is just as predictive at a later time $t^{\prime }$ as at t. Thus, from a Hutchinsonian perspective, the relation $P\left(C\right(t)\mid \mathbf{X})$ may be conserved in ecological space, even though the spatial distribution of $C(\alpha ,t)$ and/or $\mathbf{X}(\alpha ,t)$ could change in time.

In the case of presence/absence as target class, the corresponding realised niche is not conserved, in that $s\left(X_i^m\left(t\right)\right)\ne s\left(X_i^m\left(t^{\prime }\right)\right)$. This is seen in Figure 5, when compared to Figure 3, where the regression slopes for the March 2020–June 2020 and March 2020–January 2021 comparisons have lower $R^2$ values and also slopes $<1$, indicating that the scores in January or June are only about 60% of their values in March. The $R^2$ values and slopes for the $\epsilon$ comparisons in the same figures show the same effect. These differences just reflect the range expansion of COVID-19, where it has been argued that: “there is no unsuitable habitat” for COVID-19 [2]. This is linked, however, to the notion of presence, not to abundance. The fact that the score and $\epsilon$ contributions are diminishing in the presence model from March 2020 to June 2020 and January 2021 is due to the fact that the habitat variables are less able to discriminate between those cells where cases are present versus absent. The differences between June 2020 and January 2021 are much less as, by this time, a large majority of municipalities now had confirmed cases. We can also see this niche non-conservation in Figure 3 (top), using as an example the variable Internal labour flow of the municipality. We see that in March 2020 only the first three deciles of municipalities ranked by that score are associated with a positive score—niche-like—whereas in June 2020 and January 2021 60% have positive scores. Similarly, we see that in the most anti-niche-like deciles, $Decile\_02$ and $Decile\_01$, the scores are becoming less negative, indicating that those municipalities with the lowest internal labour flows are becoming less anti-niche-like. If the pandemic had a presence in every municipality, then every $s\left(X_i^m\left(t\right)\right)\to \infty$, corresponding to the fact that there can be no discrimination between where the species is present and where it is absent. Everywhere is niche-like. This would not be the case, however, for abundance, as is seen in Figure 3 (bottom), where the change is due to the fact that our relative abundance measure is associated with the top 10% of municipalities with the highest number of cases. The range expansion of COVID-19 presence also has an impact on the performance of the corresponding SDMs, as seen in Figure 4, where we see that for the all, socio-demographic, and mobility models, the corresponding AUC for the presence/absence spatial models is much less than their abundance counterparts.

We have also shown how a classification-based ENM can be used to predict abundances, with an example being the number of confirmed cases of COVID-19 for a given month. We see that the score based on all habitat variables explains approximately 90% of the variance. We believe that the ENM/SDM formalism we have developed here has the capacity to be truly epi-ecological/eco-epidemiological given the right habitat variables. The static habitat variables we have chosen cannot account for the dynamic expansion and contraction that is characteristic of epidemics and which naturally emerges from a mechanistic SIR-type modelling formalism. What is required are habitat variables that are the analogue of what enters into a differential equations type modelling environment: changes in abundances from one time period to another, for example, or even changes in those changes, as it it these variables that account for the underlying dynamics. Indeed, an ENM built on a given time slice that does not explicitly account for such variables will either underestimate or overestimate abundances, depending on whether the epidemic is in an expansion or contraction phase. We will return to this in a future publication.

We have also shown how it is possible to distinguish confounding and how this can be used as a tool to disentangle causality. As a test case we considered combining a behavioural/socio-economic/socio-demographic variable—Internal labour flow of the municipality—and a climatic variable, showing how climate as an apparently predictive habitat variable is confounded by the more predictive socio-economic and mobility factors.

An apparent limitation of our work is that we have worked at a quite coarse-grained level, that of municipalities. This is a limitation of the data used, however, not the methodology. In principle, a much finer spatial resolution, at the level of “census blocks”, for instance, could be used if data were available at that resolution for the target class and habitat. In the same vein, dynamical data that represented both changes in the local environment, such as lockdowns and hospital occupation rates, or changes in behaviour, such as mask wearing compliance, cell phone data as a proxy for short-term mobility, or vaccination rates, or a host of other factors could also be included. It would be interesting for example to determine to what extent adaptations in the pathogen were potentially reflected as changes in its niche. There is also the question of the validity of the target class data as maintaining accurate counts of confirmed cases is difficult. However, this would have little to no impact on the overall conclusions of our ENM/SDMs.

In summary, the five principal advantages of our methodology relative to standard ENM/SDMs are: (i) First and foremost, our methodology can account for the highly multifactorial nature of COVID-19 as a complex adaptive system, where there are many directly or indirectly causal factors that affect both its spread and its magnitude; (ii) It is Bayesian in nature and therefore has the advantage of being adaptive, where the incorporation of new information can be achieved naturally using Bayes theorem, which, in addition, also allows for the incorporation of human expertise by means of Bayesian priors, as well as the addition of data-based evidence; (iii) By considering the incorporation of new information in time, any time variation in the relation between a target class C and its corresponding niche variables, $\mathbf{X}$ can be tracked to determine how conserved the niche is; (iv) It permits a more profound analysis of causality and confounding through the consideration of a hierarchy of conditional probabilities; (v) It naturally permits the construction of niche models associated with different notions of niche – presence/absence, abundance, etc.—as these are simply representable as distinct classes of interest, C. The models we have presented can predict ”absolutes”, such as the number of cases in a given municipality, transversally, i.e., on the same time slice. Additionally, they can predict “relatives”, such as the relative number of cases in a given municipality longitudinally. We have shown that the corresponding models are accurate, with a high degree of correlation between predictions and actual numbers. Furthermore, that accuracy is intimately related to the incorporation of relevant niche variables such as mobility, socio-economic, and socio-demographic factors. Our goal here was not to offer a gold-standard model for prediction of the pandemic. As mentioned, there are more predictive and directly causal factors than we have included here. However, the niche of COVID-19 is immensely complex and adaptive and the incorporation of the vast array of relevant factors that determine it is a huge challenge in data collection and integration. What we have shown is that if that data can be represented in space and time, $\mathbf{X}(\alpha ,t)$, then it may be incorporated into an ENM/SDM using the methodology we have shown here. Moreover, the innovations we have shown are independent of the specific use case of COVID-19 considered here. They represent general extensions of current niche and species distribution modelling, and can be applied to any ecological system, where they are necessary. In particular, they can be applied to situations where the target and/or niche variables are changing in time. Aside from disease dynamics, invasive species, habitat, and biodiversity loss are other prime areas where time dependence is crucial and where our approach can be used. Moreover, we believe that taken in its fullest sense, where a niche/anti-niche represents the complete set of both biotic and abiotic drivers that favour where a “species” is or is not, or at least should be, a universally applicable concept, with the SDM determining the “where” and the ENM the corresponding “why”. Indeed, its applicability is only restricted in an ecological setting by just what we mean by “ecological”. If we take ecology to cover any interaction between biota and the environment then we should accept factors as mask-wearing compliance as a potential niche factor for example. Furthermore, when thinking of a “Species” distribution model, we may encouraged to be liberal in our notion of “species”, where it may represent, for instance, cases of non-transmissable diseases, such as diabetes or heart disease.

In conclusion, there is a difference between stating that ENM/SDMs generally are inappropriate vehicles for modelling a dynamic phenomenon such as the COVID-19 pandemic versus stating that a particular ENM/SDM is inappropriate. We have shown that ENM/SDMs can be generated which overcome such criticisms.