Learning complex spatial dynamics of wildlife diseases with machine learning-guided partial differential equations

Abstract

Emerging wildlife pathogens often display geographic variability due to landscape heterogeneity. Modeling approaches capable of learning complex, non-linear spatial dynamics of diseases are needed to rigorously assess and mitigate the effects of pathogens on wildlife health and biodiversity. We propose a novel machine learning (ML)-guided approach that leverages prior physical knowledge of ecological systems, using partial differential equations. We present our approach, taking advantage of the universal function approximation property of neural networks for flexible representation of the underlying dynamics of the geographic spread and growth of wildlife diseases. We demonstrate the benefits of our approach by comparing its forecasting power with commonly used methods and highlighting the obtained insights on disease dynamics. Additionally, we show the theoretical guarantees for the approximation error of our model. We illustrate the implementation of our ML-guided approach using data from white-nose syndrome (WNS) outbreaks in bat populations across the US. WNS is an infectious fungal disease responsible for significant declines in bat populations. Our results on WNS are useful for disease surveillance and bat conservation efforts. Our methods can be broadly used to assess the effects of environmental and anthropogenic drivers impacting wildlife health and biodiversity.

Impact Statement

This article presents a modern solution to a critical problem in wildlife ecology using scientific machine learning (ML). We developed an ML-guided partial differential equation (PDE) method to learn pathogen dynamics. Theoretical guarantees are derived from PDEs and ML approximation theory. This method is applied to white-nose syndrome, an infectious fungal disease caused by Pseudogymnoascus destructans, which is responsible for significant declines in several bat species in North America. The presented results are relevant for disease surveillance and conservation efforts.

1. Introduction

Different anthropogenic activities are linked to an increased frequency of the emergence of wildlife infectious diseases (Daszak et al., 2000). These diseases can spread rapidly and usually affect hosts that lack adequate response mechanisms, contributing to the decline and extinction of important species (Woodroffe, 1999; Daszak et al., 2000; McCallum, 2012; García-March et al., 2020; Plewnia et al., 2023). Moreover, emerging pathogens exhibit a variety of spatial spread and growth patterns due to landscape heterogeneity and geographic host distribution (Meentemeyer et al., 2011; Lilley et al., 2018). There is a continued need to forecast prevalence and unravel complex spatial disease dynamics (Hefley et al., 2017b). Mathematical ecological models can mechanistically describe disease spread and growth mediated by landscape characteristics (Cosner, 2008). Machine learning (ML) models provide a flexible approach to approximating functions learning from data (Hastie et al., 2009). The integration of scientific concepts from mathematical ecology with the flexibility of ML has the potential to accurately forecast disease prevalence while learning the complex spatial dynamics of disease.

The utilization of scientific knowledge for enhancing the generalization properties of ML has gained relevance in the last years. This paradigm has been presented with different names, such as theory-guided data science (Karpatne et al., 2017) or scientific ML (Dandekar et al., 2020). Hereafter, we adopt the term scientific machine learning (SciML). Recent research in SciML has shown the benefits of integrating non-linear ML algorithms with mechanistic models framed as differential equations (Baker et al., 2019; Dandekar et al., 2020; Rackauckas et al., 2020). The incorporation of prior knowledge of the underlying problem’s structure, through differential equations, yields an interpretable model with robust prediction and forecasting (Dandekar et al., 2020; Rackauckas et al., 2020). However, when using external covariates to inform predictions, the functional forms of covariate impact are usually unknown. The objective of this article is to present a novel method for the flexible incorporation of external covariates into differential equations using ML algorithms. Our approach has the potential to support ongoing research in wildlife disease modeling. For illustration, we apply our method to model white-nose syndrome (WNS) in bats.

Bats have a vital role in ecosystems as pollinators, pest controllers, and nutrient recyclers (Ramírez-Fráncel et al., 2022). Since its emergence in 2006, WNS, an infectious fungal disease caused by Pseudogymnoascus destructans (Pd), has been responsible for a significant reduction of populations of several bat species across North America (Blehert et al., 2009; Cheng et al., 2021). Recent WNS modeling approaches use mechanistic models, based on reaction–diffusion partial differential equations (PDEs), known as ecological diffusion equations (EDEs), to incorporate prior knowledge in statistical models, which helps to understand the role of habitat variation in WNS spread (Oh et al., 2023). EDEs describe pathogen spread mediated by heterogeneous landscapes, allowing spatial covariates to affect PDE coefficients (Garlick et al., 2011, 2014; Hefley et al., 2017a,b). These EDEs mechanistically characterize spatio-temporal processes, which have benefits in forecasting compared to more classical statistical approaches for modeling spatio-temporal processes (Hefley et al., 2017a,b). Nevertheless, all attempts to include spatial variation in EDEs have been limited to additive linear effects, which may not capture the real dynamics of these ecological processes.

In this work, we propose a novel method that leverages prior physical knowledge of ecological systems through PDEs, with flexible representations of underlying dynamics using ML algorithms. By incorporating ML algorithms, we augment EDEs. To demonstrate our approach, we use neural networks (NNs) as universal function approximators, to learn underlying unknown dynamic patterns of spreading wildlife diseases, specifically applying it to Pd in the contiguous U.S. We remark that, while in this article we use NNs, our developed approach easily allows the use of other ML algorithms, such as k-nearest neighbors, regression trees, or random forests. Finally, we derive theoretical results on the approximation error of our method that apply to any supervised ML model embedded in our EDE system.

The rest of this article is organized as follows. In Section 2, we provide some background on EDEs and define notation essential to the presentation of our ML-guided diffusion model. Section 3 describes the case study that motivated this article and presents the modeling details of the proposed method. In Section 4 we summarize and interpret the results of the WNS application. Finally, in Section 5 we discuss the results and present a conclusion. The theoretical results derived in this work can be found in Appendix B. Appendix C provides pseudocode for the derived equations in the background section. Appendix D presents pseudocode for the main method of this article.

2. Background: ecological diffusion equations

This section will serve as background material for understanding the EDEs, which have been used in multiple previous studies for wildlife diseases (Garlick et al., 2011; Hefley et al., 2017b; Oh et al., 2023, 2024). We use such an equation to characterize the spread and growth of wildlife diseases over space and time, influenced by local conditions in varying landscapes. EDEs can merge two important epidemiological processes. First, we have a term that describes the spatial spread (diffusion) of the modeled pathogen density, which contains a spatially varying coefficient that is inversely related to aggregation (high diffusion means little pathogen aggregation, while low diffusion means high aggregation). Second, we have a growth component that characterizes population dynamics with a modulation coefficient (growth rate) that depends on local spatial information. Understanding the growth and spread of pathogens, particularly potential drivers of these processes, is critical for developing effective control strategies.

Let $ \mathcal{S}=\left\{\left(x,y\right)\in {\mathrm{\mathbb{R}}}^2|0\le x\le {L}_1,0\le y\le {L}_2\right\} $ be a spatial domain of interest, and $ \mathcal{T}=\left(0,T\right]\subset {\mathrm{\mathbb{R}}}^{+} $ be a temporal domain. We are interested in the reaction–diffusion PDE system:

(2.1) $$ \frac{\partial }{\partial t}u\left(x,y,t\right)={\nabla}^2\left[\mu \left(x,y\right)u\left(x,y,t\right)\right]+\lambda \left(x,y\right)u\left(x,y,t\right), $$

(2.2) $$ \hskip-11em u\left(x,y,0\right)=\left\{\begin{array}{cc}\eta, & \mathrm{if}\left(x,y\right)=\omega \\ {}0,& \mathrm{if}\left(x,y\right)\ne \omega, \end{array}\right. $$

(2.3) $$ \hskip-21.4em u\left(\left(x,y\right)\in \partial \mathcal{S},t\right),=0\hskip1em \forall t>0, $$

where $ u\left(x,y,t\right) $ represents the density of a dispersing pathogen at location $ \left(x,y\right)\in \mathcal{S} $ and time $ t\in \mathcal{T} $ , $ \partial \mathcal{S} $ represents the boundary of the spatial domain, and $ \omega $ represents the location of the initial introduction of the pathogen with initial density $ \eta . $ Spatial spread is captured by the Laplacian operator $ {\nabla}^2=\frac{\partial^2}{\partial {x}^2}+\frac{\partial^2}{\partial {y}^2} $ . The function $ \mu \left(x,y\right) $ is the spatially dependent diffusion coefficient such that $ 0<\mu <\infty $ , and the function $ \lambda \left(x,y\right) $ is a spatially varying growth rate. In many cases, the spatial variation of both $ \mu $ and $ \lambda $ is at a much finer scale than the scale of variation for the population, $ u $ . Particularly, in the case of wildlife diseases, the population scale is hundreds of kilometers, and the scale of habitat variation, $ \mu $ and $ \lambda $ , is meters. This mismatch creates a computational difficulty resolved by homogenization.

2.1. Homogenization and analytic approximation

In this subsection, we describe the process of homogenization and its application to the EDEs in equations (2.1)–(2.3). Then, we provide an approximate analytical solution that offers an advantage in computation time and memory requirements compared to numerical solvers (Oh et al., 2024).

Homogenization is an asymptotic method (Holmes, 2012) that leverages the two different scales of interest to simplify the EDE system, (2.1)–(2.3). This multi-scale method frames the PDE on the large scale by locally averaging the fine-scale variation, producing locally constant coefficients in a PDE system with known solutions. After solving this homogenized PDE system, the solution is adjusted by the fine-scale variation.

Following known homogenization procedures (Garlick et al., 2011, 2014), we derive a homogenized EDE system:

(2.4) $$ \frac{\partial }{\partial t}c\left(x,y,t\right)=\overline{\mu}\left(x,y\right){\nabla}^2c\left(x,y,t\right)+\overline{\lambda}\left(x,y\right)c\left(x,y,t\right), $$

(2.5) $$ c\left(x,y,0\right)=\left\{\begin{array}{cc}\eta \bar{\mu}\left(x,y\right),& \mathrm{if}\left(x,y\right)=\omega, \\ {}\hskip-3.5em 0,& \mathrm{if}\left(x,y\right)\ne \omega, \end{array}\right. $$

(2.6) $$ \hskip-10em c\left(\left(x,y\right)\in \partial \mathcal{S},t\right)=0,\hskip1em \forall t>0, $$

where $ c\left(x,y,t\right) $ is the homogenized approximated solution, related to the original PDE solution, $ u\left(x,y,t\right) $ , by

$$ u\left(x,y,t\right)\approx \frac{c\left(x,y,t\right)}{\mu \left(x,y\right)}. $$

The new homogenized coefficients are locally constant averages taken over the large scale. The homogenized diffusion coefficient, $ \overline{\mu}\left(x,y\right) $ , is defined as

(2.7) $$ \overline{\mu}\left(x,y\right)=\frac{\int_{\mathcal{A}}1 dxdy}{\int_{\mathcal{A}}\frac{1}{\mu \left(x,y\right)} dxdy}, $$

and the homogenized growth rate, $ \overline{\lambda}\left(x,y\right) $ , is

(2.8) $$ \overline{\lambda}\left(x,y\right)=\frac{\overline{\mu}\left(x,y\right)}{\mid \mathcal{A}\mid }{\int}_{\mathcal{A}}\frac{\lambda \left(x,y\right)}{\mu \left(x,y\right)} dxdy. $$

The homogenization region $ \mathcal{A}\subset {\mathrm{\mathbb{R}}}^2 $ , with area $ \mid \mathcal{A}\mid $ , over which we calculate the homogenized parameters, is defined at an intermediate scale between the large and the small scale. Using separation of variables and the Fourier series expansion, Oh et al. (2024) derive the following approximate solution to the system (2.4)–(2.6). Then, for sufficiently large positive integers $ M $ and $ N $ , a truncated series solution to the homogenized system is

(2.9) $$ {\tilde{c}}_{M,N}\left(x,y,t\right)=\sum \limits_{m=1}^M\sum \limits_{n=1}^N{c}_{m,n}\left(x,y,t\right), $$

where

(2.10) $$ {c}_{m,n}\left(x,y,t\right)=\frac{4\eta \overline{\mu}\left(\omega \right)}{L_1{L}_2}\sin \left(\tilde{m}{\omega}_1\right)\sin \left(\tilde{n}{\omega}_2\right)\sin \left(\tilde{m}x\right)\sin \left(\tilde{n}y\right)\exp \left(\overline{\lambda}\left(x,y\right)t-\overline{\mu}\left(x,y\right)t\left({\tilde{m}}^2+{\tilde{n}}^2\right)\right), $$

with $ \tilde{m}= m\pi /{L}_1 $ , $ \tilde{n}= n\pi /{L}_2 $ , and $ \omega =\left({\omega}_1,{\omega}_2\right) $ . Finally, the analytical approximated solution to $ u\left(x,y,t\right) $ , the solution of the EDE system (2.1)–(2.3), is

(2.11) $$ {\tilde{u}}_{M,N}\left(x,y,t\right)=\frac{{\tilde{c}}_{M,N}\left(x,y,t\right)}{\mu \left(x,y\right)}. $$

Refer to Oh et al. (2024) for a comparison of the analytical approximate solution with the commonly used Forward Time Centered Space solver method (finite differences), showing the computational benefits provided by the analytical approximate solution. Appendix C presents pseudocode for this section.

3. Methodology

In this section of this article, we will present the WNS surveillance data. Then, we introduce the ML-guided PDE method for ecological diffusion, which is the main contribution of the presented work. Finally, we specify the model and describe its fitting process.

3.1. White-nose syndrome data

The dataset used in this article has geo-referenced samples collected across the contiguous U.S. by state and federal agencies and tested for the presence of Pd at the U.S. Geological Survey National Wildlife Health Center (Madison, WI) between April 2007 and April 2020 (Ballmann et al., 2021). Samples were taken from bat carcasses, bat wing punch biopsies, forearm or muzzle skin, wing swabs, guano, roost substrate swabs, cave (wall/ceiling) swabs, and soil, all of them tested for Pd infection using a quantitative PCR (Muller et al., 2013). We use bat samples and environmental samples, as Pd is known to persist in the environment (Lorch et al., 2013), to model the presence or absence of Pd rather than the disease (WNS) that this pathogen causes.

Two outbreaks of the disease have been reported: one in New York in 2006 and another in Washington in 2016. However, we only model the data from the New York outbreak by excluding all positive samples in Washington state. We have the spatial coordinates (latitude and longitude) of 18,515 negative and 1,557 positive observations (Figure 1). Additionally, we consider five spatial covariates that are believed to influence Pd dynamics (Oh et al., 2023): percentage of tree canopy cover (canopy), linear hydrography including streams/rivers, braided streams, canals, ditches, artificial paths, and aqueducts (waterways), topographic ruggedness index which is related to the magnitude of elevation (TRI), number of coal mines per $ 10 $ km $ \times 10 $ km grid cell (mines), and percentage of karst geomorphology (karst). See Additional Figure A1 in Appendix A. We transformed these covariates to lie in the interval $ \left[0,1\right] $ . The covariates were selected following recommendations from bat biologists on the Strategic Pd Surveillance Advisory Team and have been used in previous modeling efforts (Oh et al., 2023).

Figure 1. White-nose syndrome (WNS) is an infectious fungal disease in bats caused by Pseudogymnoascus destructans. Presented here are the geographic locations where WNS samples were collected by the USGS WNS surveillance team between 2006 and 2016. There were 1,557 positive tests, represented by the colored dots, with lighter colors corresponding to more recently taken samples and darker colors corresponding to earlier samples. There were 18,515 negative samples, represented with small black dots.

3.2. Machine learning-guided partial differential equations

Using the notation from Section 2, the analytical approximated solution $ {\tilde{u}}_{M,N}\left(x,y,t\right) $ in (2.11) to the PDE system (2.1)–(2.3) depends on the two functions $ \mu \left(x,y\right) $ and $ \lambda \left(x,y\right) $ , which are the unknown spatial coefficients. Recent research has shown that NNs can be used as function approximators to learn components of differential equation systems (Dandekar et al., 2020; Rackauckas et al., 2020). Following this principle, we propose the use of ML algorithms to learn unknown dynamic patterns related to the diffusion ( $ \mu $ ) and growth ( $ \lambda $ ) in our system. Let $ {NN}_{\theta_1}\left(z\left(x,y\right)\right) $ and $ {NN}_{\theta_2}\left(z\left(x,y\right)\right) $ be two supervised ML algorithms with trainable parameters $ {\theta}_1 $ and $ {\theta}_2 $ respectively, taking an input vector $ z\left(x,y\right)\in {\mathrm{\mathbb{R}}}^p $ of spatial covariates. We consider $ {NN}_{\theta_1} $ to approximate $ \log \mu $ (instead of $ \mu $ directly to ensure positivity), and $ {NN}_{\theta_2} $ to approximate $ \lambda $ . Additionally, considering that we have a set of observed values $ {u}_{\mathrm{data}}={\left\{{x}_i,{y}_i,{t}_i,{u}_i\right\}}_i^n $ that will be modeled using (2.1)–(2.3), let $ \mathrm{\ell}\left({\tilde{u}}_{M,N},{u}_{\mathrm{data}}|\theta \right) $ be a loss function depending on $ \theta =\left({\theta}_1,{\theta}_2\right) $ through $ {\tilde{u}}_{M,N} $ , that measures how well the approximated solution matches WNS surveillance data. We are interested in finding the set of parameters that minimize our loss function, i.e.,

$$ \hat{\theta}=\arg \hskip0.3em \underset{\theta }{\min}\mathrm{\ell}\left({\tilde{u}}_{M,N},{u}_{\mathrm{data}}|\theta \right). $$

Hereafter, we assume that $ {NN}_{\theta_1} $ and $ {NN}_{\theta_2} $ are two feedforward artificial NNs with trainable weights and biases $ {\theta}_1 $ and $ {\theta}_2 $ , respectively. However, we remark that any two ML algorithms can be used, and they do not need to be the same. See Appendix B for the theoretical guarantees for the approximation error for $ u\left(x,y,t\right) $ under the proposed procedure. We showed that the approximation error will be small as long as: (1) homogenization assumptions are satisfied; (2) the truncation bounds $ M $ and $ N $ on the approximate solution are sufficiently large; and (3) the ML algorithms are properly trained. Finally, note that the function $ \mathrm{\ell} $ can be the mean squared error or the negative log-likelihood associated with an assumed distribution in our data (e.g., binomial distribution for binary observations), in which case we would need to perform a post-processing transformation on $ {\tilde{u}}_{M,N} $ using a link function.

3.3. Model specification and fitting

Herein, we describe our statistical model for wildlife diseases, which will be applied to the WNS data. Note that, previously, we assumed data were direct observations of $ u\left(x,y,t\right) $ however, in the WNS data, we observe a function of $ u\left(x,y,t\right) $ , namely, the binary observations of tests for the presence of a pathogen. Thus, we connect the observed binary response to the solution of our PDE by the logit link function defined as $ {g}^{-1}(a)=\exp (a)/\left(1+\exp (a)\right) $ . Let $ {v}_i $ represent the binary response variable (1 for infected, 0 otherwise) sampled at location $ \left({x}_i,{y}_i\right) $ and time $ {t}_i $ for $ i=1,\dots, n $ . We denote $ {u}_{\mathrm{data}}={\left\{{x}_i,{y}_i,{t}_i,{v}_i\right\}}_i^n $ and consider the model:

(3.1) $$ \hskip-8.43em {v}_i\sim \mathrm{Bernoulli}\left({p}_i\right), $$

(3.2) $$ \hskip-7.5em {p}_i={g}^{-1}\left(u\left({x}_i,{y}_i,{t}_i\right)\right), $$

(3.3) $$ \hskip0.43em \frac{\partial }{\partial t}u\left(x,y,t\right)={\nabla}^2\left[\mu \left(x,y\right)u\left(x,y,t\right)\right]+\lambda \left(x,y\right)u\left(x,y,t\right), $$

(3.4) $$ \hskip-13.2em u\left(x,y,0\right)=\left\{\begin{array}{cc}\eta, & \mathrm{if}\left(x,y\right)=\omega, \\ {}0,& \mathrm{if}\left(x,y\right)\ne \omega, \end{array}\right. $$

(3.5) $$ \hskip-21.43em u\left(\left(x,y\right)\in \partial \mathcal{S},t\right),=0\hskip1em \forall t>0, $$

(3.6) $$ \hskip-14em \mathrm{\log}\mu \left(x,y\right)=N{N}_{\theta_1}\left(z\left(x,y\right)\right), $$

(3.7) $$ \hskip-11.3em \lambda \left(x,y\right)={NN}_{\theta_2}\left(z\left(x,y\right)\right), $$

where $ {p}_i=P\left({v}_i=1\right) $ is the probability of observing a positive outcome, $ g $ is the logit link function, and $ u $ is the solution to the PDE system (3.3)–(3.5) with $ \mu \left(x,y\right) $ representing the spatially varying diffusion coefficient, and $ \lambda \left(x,y\right) $ the spatially varying growth rate. $ {NN}_{\theta_1} $ and $ {NN}_{\theta_2} $ represent two ML algorithms with trainable parameters $ {\theta}_1 $ and $ {\theta}_2 $ , respectively, and $ z\left(x,y\right) $ is the 5-dimensional input vector of spatial covariates described at the beginning of this section. The location of the introduction point, $ \omega $ , and the initial density, $ \eta $ , in equations (2.1)–(2.3) were assumed to be known, using estimates from a recent study (Oh et al., 2023). Assuming $ {\tilde{u}}_{M,N} $ approximates $ u $ and denoting $ {\hat{p}}_i={g}^{-1}\left({\tilde{u}}_{M,N}\left({x}_i,{y}_i,{t}_i\right)\right) $ , we minimize the negative log-likelihood

$$ \mathrm{\ell}\left({\tilde{u}}_{M,N},{u}_{\mathrm{data}}|\theta \right)=-\frac{1}{n}\sum \limits_{i=1}^n{v}_i\log \left({\hat{p}}_i\right)+\left(1-{v}_i\right)\log \left(1-{\hat{p}}_i\right). $$

In Figure 2, we illustrate how NN predictions inform our PDE model during model fitting. The two NNs $ {NN}_{\theta_1} $ and $ {NN}_{\theta_2} $ have the same architecture with $ 2 $ hidden layers, $ 512 $ neurons in each hidden layer, and ReLU as the activation functions. The NNs were defined using the Tensorflow (TF) environment in Python. Regarding optimization, we use Adam optimizer with a scheduled learning rate starting at $ 0.0005 $ , with a decay rate of $ 0.96 $ per epoch. We observed a challenging optimization problem when the weights and biases are randomly initialized, specifically for the NN $ {NN}_{\theta_2} $ dedicated to predicting $ \log \left(\mu \right) $ . Therefore, we suggest a careful initialization of the bias of the output layer following the estimation of the wavefront speed heuristic presented in Shigesada and Kawasaki (1997). To do this, we let $ \mathcal{P} $ represent the collection of the $ P $ indexes corresponding to the positive samples in the training dataset, then the bias initialization in the output layer of $ {NN}_{\theta_2} $ is

$$ \log \left(\frac{1}{P}\sum \limits_{i\in \mathcal{P}}\frac{d\left(\left({x}_i,{y}_i\right),\left({\omega}_1,{\omega}_2\right)\right)}{\pi {t}_i}\right)\approx 22.2, $$

where $ \left({\omega}_1,{\omega}_2\right) $ is the coordinate pair for the introduction point, $ \left({x}_i,{y}_i\right) $ is the coordinate pair for the $ i- $ th positive sample, $ {t}_i $ is the time the $ i- $ th positive sample was taken, and $ d\left(\cdot, \cdot \right) $ represents the Euclidean distance. See Appendix D for implementation details of model fitting.

Figure 2. Diagram of neural networks (NNs) informing the partial differential equation (PDE) to model pathogen spread and growth. The covariates, tree canopy cover (canopy), linear hydrography (waterways), topographic ruggedness index (TRI), number of coal mines (mines), and karst geomorphology (karst) are used as input. The NNs (red and blue boxes) take these covariates and predict the unknown spatially varying log-diffusion $ \left(\log \mu \right) $ and growth $ \left(\lambda \right) $ coefficients which characterize the PDE. The solution of the PDE, given the predicted coefficients, is post-processed and contrasted with the observed data using a loss function. The loss guides the training process to refine the NN’s predictions.

3.4. Validation and comparison with PDE with spatial additive linear effects

For model evaluation in a forecasting task, we left the last 12 months of samples as test data. Thus, we evaluated our model’s ability to detect the Pd pathogen by measuring the area under the Receiver Operating Characteristics curve (ROC-AUC) and the area under the Precision/Recall curve (PR-AUC) metrics in both training and testing datasets. We chose these metrics according to Saito and Rehmsmeier (2015) recommendations for imbalanced classification problems. The baseline for ROC-AUC is $ 0.5 $ , meaning that a random classifier would obtain $ 0.5 $ in this metric. For PR-AUC, the baseline depends on the positive/negative ratio: baseline = (number of positives)/(number of positives + number of negatives). Therefore, for the training dataset, the baseline was $ 0.084 $ , while for the testing dataset, the baseline was $ 0.044 $ .

To contrast our proposed method with the existing efforts on EDE for WNS modeling, we considered an alternative version of the model that uses the PDE with spatial additive linear effects. We assumed linear functional forms for the diffusion coefficient

$$ \log \mu \left(x,y\right)={\alpha}_0+{\alpha}^{\top }z\left(x,y\right), $$

and growth rate

$$ \lambda \left(x,y\right)={\gamma}_0+{\gamma}^{\top }z\left(x,y\right). $$

where $ {\alpha}_0 $ and $ {\gamma}_0 $ are intercepts, and $ \alpha $ and $ \gamma $ are vector coefficients are associated with the spatial covariate vector $ z\left(x,y\right) $ . This approach is similar to the ones presented in Oh et al. (2023) and Oh et al. (2024). Hereafter, this model will be called the EDE with spatial additive linear effects model.

4. Results

After fitting the model, we obtained a characterization of the prevalence of Pd, i.e., we can obtain $ {\hat{p}}_i={g}^{-1}\left({\tilde{u}}_{M,N}\left(x,y,t\right)\right) $ at any time $ t $ and space point $ \left(x,y\right) $ . Additionally, we learned two functions $ \log \hat{\mu}\left(x,y\right)={NN}_{{\hat{\theta}}_1}\left(z\left(x,y\right)\right) $ and $ \hat{\lambda}\left(x,y\right)={NN}_{{\hat{\theta}}_2}\left(z\left(x,y\right)\right) $ that can help agencies target management actions by allowing them to understand where the disease can spread faster ( $ \mu $ ) and where it is expected to grow more readily ( $ \lambda $ ). Figure 3 shows the probability of infection maps for April 2010 and April 2020. We observe a heterogeneous diffusive behavior, which is important for differentiating the impact of Pd reaching different zones. Additionally, the upper maps of Figure 3 depict a process with a defined boundary between high and low probabilities of infection, which may enable management agencies to allocate testing resources based on the likelihood of detecting the pathogen (Oh et al., 2023).

Figure 3. Machine learning-guided partial differential equation approximations of the probability of presence of Pseudogymnoascus destructans (Pd) (upper) vs. observed data (lower) for April 2010 (left) and April 2020 (right). The observed data color represents time (earlier times have a darker color). We observe a heterogeneous diffusive behavior, which is important for differentiating the impact of Pd reaching different zones.

The training dataset (samples up to April 2019, Pd imbalanced ratio $ =0.084 $ ) obtained a PR-AUC $ =0.359\pm 0.088 $ and ROC-AUC $ =0.849\pm 0.070 $ , for testing data (samples after April 2019, Pd imbalanced ratio $ =0.044 $ ), we obtained PR-AUC $ =0.085\pm 0.010 $ and ROC-AUC $ =0.729\pm 0.024 $ . The ROC-AUC value represents the probability of correctly predicting a random positive and negative example. Thus, the expected ROC-AUC of a random classifer is 0.5. The higher the ROC-AUC, the better the model’s ability to distinguish between positive and negative samples. The Pd imbalanced ratio is the PR-AUC baseline, which can be interpreted as the expected performance of a random classifier for each dataset. Therefore, considering the average PR-AUC, we perform about $ 4 $ times better than a random classifier for the training dataset (PR-AUC $ =4\times \mathrm{Pd}\;\mathrm{imbalanced}\ \mathrm{ratio} $ in training data), and $ 2 $ times better than a random classifier for the testing dataset (PR-AUC $ =2\times \mathrm{Pd}\;\mathrm{imbalanced}\ \mathrm{ratio} $ in testing data). Table 1 summarizes the validation metric results of 30 runs of our model with median, mean, and standard deviation compared with the linear model. In all metrics, our ML-guided PDE method outperformed the PDE with a spatial additive linear effects model.

Table 1. Validation metrics on Pseudogymnoascus destructans (Pd) detection for our machine learning-guided partial differential equation (PDE) method (machine leaning [ML]-guided) and the PDE with spatial additive linear effects model (linear) for the train and test datasets

^a Area under the receiver operating characteristics curve.

^b Area under the precision/recall curve.

Although improving computation time was not an objective of this development, we report the computation time of our experiments. The PDE model with spatial additive linear effects was completed in an average $ 19.76\pm 6.63 $ minutes, whereas the ML-guided PDE method completed $ 14.16\pm 4.33 $ minutes. Both approaches were trained using the same computational characteristics (Google Colab, RAM 12.7 GB, Disk 78.2 GB, and 2 CPU cores). Additionally, we note that the average loss function value (negative log-likelihood) in the test dataset for our ML-guided PDE method was $ 0.705\pm 0.015 $ , and for PDE with spatial additive linear effects model was $ 0.720\pm 0.034 $ (Figure 4).

Figure 4. Left: loss function value for test dataset in our machine learning-guided partial differential equation (PDE) method (machine learning [ML]-guided) vs. the PDE with spatial additive linear effects model (linear). Right: training computation time (in seconds) for ML-guided vs. linear model.

Furthermore, we did a literature review to externally validate our outcomes with the reported spread and growth of WNS in previous studies. Our NN predictions (upper left, Figure 5) presented a high diffusion and low growth in the karst landscapes which aligns with reported findings about the spread of Pd in the karst regions of Appalachian Mountains (Maher et al., 2012; Hoyt et al., 2021) (arrow 1 in zoomed states in Figure A2 in Appendix A). Similar behavior can be observed with the linear prediction (lower left, Figure 5), but with higher values. In the NN predictions (upper panels in Figure 5), we observed a color pattern characterizing an area with relatively high diffusion and low growth in northern Florida, across Georgia and South Carolina states (arrow 2 in zoomed states in Figure A2 in Appendix A). Such a pattern may indicate a predicted natural barrier because the combination of high diffusion (low pathogen residence time) and low growth rates would result in no aggregation of pathogens in this area. This finding supports the hypothesis of this area being the southern limit of Pd distribution (Hoyt et al., 2021). The linear predictor (lower panels in Figure 5) does not match with this hypothesis, as it shows a low growth area with a relatively low diffusion. Conversely, Pd is known to grow best at low temperatures (12.5–15.8 $ {}^{\circ } $ C) (Verant et al., 2012), and our approximated growth coefficient (in both models) presents high values in the north-central region, which matches zones of low monthly average temperatures.

Figure 5. Maps with approximated log-diffusion coefficient $ \log \mu $ and growth coefficient $ \lambda $ from neural networks (upper) and linear function (lower). Color scales are different to allow easier visualization of details in each figure. Approximated diffusion and growth values are only interpretable within the boundaries of the continental USA. Predicted values are not interpretable for large water bodies (e.g., the ocean and the Great Lakes).

Figure 6 was created to explore the non-linear relationships present in the diffusion and growth coefficients learned by the NNs and compare them with the linear models. We plotted log-diffusion and growth as a function of each covariate, keeping the remaining covariates fixed. We let one covariate vary at a time and fixed the other covariates at the 0.5 level. In general, we observed that the functional forms for all covariates in diffusion and growth were close to a second-degree polynomial, and in a few cases, a third-degree polynomial for the NN. The growth coefficient showed a more complex behavior compared with the diffusion one. Nevertheless, we remark that these relationships should not be over-interpreted.

Figure 6. Functional relationship between each variable versus the log-diffusion coefficient $ \log \mu $ (upper) and growth coefficient $ \lambda $ (lower) from neural networks (NNs) (right) and a linear model (left). Each covariate is varied while the remaining covariates are fixed at 0.5. Note the different values on the y-axis between linear and NN models.

Finally, to provide interpretation from our results, we leveraged the SHAP (SHapley Additive exPlanations) values technique (Lundberg and Lee, 2017). SHAP values are used to explain the covariate contributions to the predictions of each observation. In particular, we used the Deep SHAP approach, which gives a high-speed approximation for SHAP values in deep learning models. Figure 7 displays bee swarm plots of SHAP values for $ \mathrm{5,000} $ randomly sampled locations for log-diffusion and growth NN models. In these figures, covariates are ranked from top to bottom by their mean absolute SHAP value. For each covariate, each of the randomly sampled locations has a point distributed along the horizontal axis by its SHAP value. SHAP values with high density are represented by stacking the points vertically. For each point, the color represents its raw value. Therefore, in Figure 7, we observe that canopy cover is the most influential covariate for log-diffusion and growth. Higher values of canopy are associated with lower values of diffusion, while lower values of canopy are associated with higher values of diffusion. Similarly, large values of karst have a negative association with diffusion, and lower values of karst are associated with higher diffusion rates. For Pd growth, higher values of karst have a positive association, while lower values have a negative association. Lastly, large values of canopy are associated with low growth rates, and lower canopy values have a positive association with growth.

Figure 7. Bee swarm plots of SHapley Additive exPlanations (SHAP) values from $ \mathrm{5,000} $ randomly sampled locations for log-diffusion (upper) and growth (lower) neural network models. SHAP values explain the covariate contributions to the predictions of each observation. Covariates are ranked from top to bottom by their mean absolute SHAP value (shown in parentheses beside the name). For each covariate, each location has a point distributed along the horizontal axis by its SHAP value. SHAP values with high density are represented by stacking the points vertically. Color represents covariate raw values.

5. Conclusions and discussion

We developed and implemented a modern solution to a critical problem in wildlife ecology, using ML-guided PDEs. Our approach exploits prior physical knowledge of ecological systems, encoded through PDEs, using ML algorithms for a flexible representation of interpretable coefficients. Moreover, we provided the theoretical guarantees derived from the analytical approximation for our PDE and the universal approximation theorem B.4. Additionally, we demonstrated the benefits of our method by comparing its forecasting power with the commonly used PDE with a spatial additive linear effects model. However, while our validation metrics outperformed existing methods, there is room for improvement regarding the pathogen detection power of the presented model.

Some of our practical findings are compatible and support hypotheses from existing WNS literature, which gives us confidence that our methods have the potential to be broadly applied to other wildlife diseases or other ecological processes that exhibit growth and spread. Some of the areas that can benefit from our methods are plant and animal population growth and dispersal, modeling the expansion of invasive species, and studying wildlife migration. Furthermore, the proposed method is a general approach that can be adapted to different PDE system structures, including other initial and boundary conditions, with different diffusion and growth mechanisms.

An important upside of the presented method is the modeling flexibility for the diffusion and growth coefficients associated with the EDEs. Unlike using additive linear effects or other hard parametric assumed relationships, our approach does not require users to know or specify a functional form of the impact of covariates on the processes of interest. Thus, the incorporation of the effect of several different covariates is much easier to do with the embedded ML algorithms. Our method showcased the value of SciML in wildlife epidemiology, providing a powerful framework for modeling complex ecological processes. This article may inspire forthcoming research focused on modeling wildlife diseases considering a wide range of covariates related to important environmental processes, such as climate change, and anthropogenic activities, such as pollutant emissions. Finally, we believe that exploring and comparing different ML algorithms, including more complex NN architectures, is a promising future path of investigation with the potential of capturing intricate non-linear relationships that simpler models might miss.