\section{Introduction}

When a species or pathogen invades a novel environment, determining
when and where it was first introduced is of interest to scientists,
policymakers, and the public. Knowing the location and date of introduction
would further our understanding of the ecology of the process, increase
our capacity to build predictive models, and could inform policies
that aim to prevent new introductions. Although the location and time
of introduction is known in some cases, such as the introduction of
myxomytosis to control rabbits in Australia (\citealt{ratcliffe1952myxomatosis}),
there are many examples where the location and time of introduction
are unknown and must be inferred from indirect data.

The general problem of modeling the spatio-temporal dynamics of invasion
is well-represented in the literature (e.g., \citealt{mollison1986modelling,gibson1996fitting,wikle2003invasion});
however, when the introduction is not directly observed there is a
small number of limited methods that might be used to determine the
location and time of introduction from indirect data. For example,
a method known as geographic profiling has been introduced to estimate
the position of sources of invasive species and diseases using occurrence
records (\citealt{le2011geographic,stevenson2012geographic,verity2014spatial}),
but this approach ignores the temporal dynamics and therefore cannot
estimate the initial source (however see \citealt{mohler2012geographic}).
Similarly, Bayesian models have been developed to estimate the spatio-temporal
dynamics of colonization (e.g., \citealt{cook2007bayesian,catterall2012accounting,broms2016dynamic}),
but these methods are limited. For example, the \citet{cook2007bayesian}
approach requires that the colonization time be known. \citet{catterall2012accounting}
demonstrate how to infer the time of colonization; however, this method
was developed to understand the dynamics after the initial introduction
and not to obtain estimates of the time and location of the first
introduction. Lastly, \citet{caley2015inferring} developed a Bayesian
spatio-temporal model to infer the time of introduction and geographic
distribution of an invasive species from occurrence records; however,
this approach does not provide a framework for inferring the location
of introduction.

If the time and location of introductions were known, then tools commonly
used to map the distribution of species and diseases could be used
to build predictive models and infer the factors that influence the
occurrence of novel introductions. Novel introductions, however, are
rarely directly observed. Instead, the outcome of a successful introduction
is observed by the collection of indirect data such as species or
disease occurrence records. In principal, indirect spatio-temporal
data and models can be used infer historical events such as infection
histories (\citealt{salje2018reconstruction}), changes in climate
(\citealt{tingley2012piecing,tipton2016reconstruction}), and, as
we demonstrate, introductions of pests and pathogens.

Our objective is to show how point process models, a commonly used
statistical approach for modeling the distribution of species and
diseases, can be placed within a hierarchical modeling framework to
obtain statistically principled estimates of the date and location
of introduction. Although we provide a general framework, our work
is motivated by the need to understand when and where $\textit{Pseudogymnoascus destructans}$,
the pathogen that causes white-nose syndrome (WNS), was first introduced
into the United States. White-nose syndrome was first detected in
2006 using photographic evidence of a bat from Howes Cave, near Albany,
New York (Fig. 1; \citealt{blehert2009bat,frick2010emerging}). The
pathogen, $\textit{P. destructans}$, has since spread throughout
the eastern and midwestern United States resulting in high mortality
rates among several species of cave-hibernating bats. Statistical
methods capable of estimating the date and location where $\textit{P. destructans}$
was first introduced into the United States could further our understanding
of the etiology of the disease and provide critical information on
the process of introduction. For example, WNS was recently detected
on a western subspecies of little brown bat (\textit{Myotis lucifugus})
occurring in the state of Washington that is likely the result of
a novel introduction originating from the eastern United States (\citealt{lorch2016first}).
Our modeling approach could eventually be used to test the hypothesis
that the Washington WNS case is the result of a novel introduction
and, if confirmed, identify areas at high risk for future introductions.
\begin{figure}[H]
\begin{centering}
\includegraphics[scale=0.7]{Fig_1}
\par\end{centering}
\noindent \textbf{Fig. 1.} Location and infection status of individual
bat samples from four species tested for white-nose syndrome. The
first detection in 2006 was based on photographic evidence only. The
locations of individual samples were jittered so that cases from the
same location did not overlap. The four species of bats tested include:
little brown bat (\textit{Myotis lucifugus}), big brown bat (\textit{Eptesicus
fuscus}), northern long-eared bat (\textit{Myotis septentrionalis}),
and tri-colored bat (\textit{Perimyotis subflavus}).
\end{figure}

\section{Dynamic Spatio-temporal Statistical Models}

Dynamic spatio-temporal statistical models are used in fields such
as ecology (\citealt{wikle2003invasion,hooten2008doves,williams2017integrated,hefley2017mechanism}),
epidemiology (\citealt{smith2002predicting}), and meteorology (\citealt{wikle1998hierarchical,wikle2001spatiotemporal}).
The defining property of a dynamic statistical model is an explicit
description of how future observations (spatial data) evolve from
past observations (\citealt{wikle2010general,mohler2012geographic,cressie2011statistics,wikle2018stwr}).
For example, statistical implementations of dynamic agent-based models
are constructed using a set of deterministic or stochastic rules governing
the actions of individual agents at each time step (\citealt{hooten2010agent,hooten2010statistical,keith2013agent,heard2015agent,mcdermott2017hierarchical}).

Consider the classic example of a dynamic agent-based model where
an individual performs a random walk. At each time step, the individual
moves a random distance in a random direction. An inherent characteristic
of dynamic models is that the initial conditions must be specified.
For example, when performing a random walk, the location and time
the individual begins walking must be specified. Although a Eulerian
(individual-based) frame of reference can be used to understand dynamic
processes, a Lagrangian (population-level) perspective is also common
(\citealt{hooten2017animal}). In ecology, for example, dynamic models
of populations are expressed as partial differential equations (PDEs),
which are used to describe the global behavior of a large number of
individuals (\citealt{holmes1994partial,hooten2017animal}).

In data-driven modeling applications, hierarchical (conditional) model
specifications can be used to estimate unknown parameters by connecting
data to the underlying dynamic process (\citealt{cressie2011statistics,wikle2018stwr}).
Many applications of dynamic statistical models have focused on estimating
parameters that govern the process (e.g., rate of spread), but rarely
has formal statistical inference about the initial conditions been
of interest (but see \citealt{hefley2017spatial}). Instead, the initial
conditions are modeled phenomenologically and treated as nuisance
parameters (\citealt{wikle2003hierarchical,hooten2008doves,williams2017integrated,hefley2017mechanism}).

We demonstrate how incorporating scientific knowledge about the initial
conditions of a dynamic process provides novel insight into the outbreak
of an infectious disease or invasion of a pest. By justifiably assuming
that the pathogen or pest was introduced at points occurring at unknown
locations and times, we specify a point process to model the initial
conditions. Point process models are commonly used for species distribution
and disease mapping (\citealt{renner2015point}), but can also be
incorporated into hierarchical Bayesian models to learn about latent
(unobserved) processes (e.g., \citealt{BB2L}). In the following sections,
we provide a specific example using a PDE of ecological diffusion;
however, our approach is applicable to commonly used dynamic models
that require the specification of initial conditions such as agent-based
models, contact process models, rule-based approaches (e.g., cellular
automaton), integral equations, and integro-difference equations.

\subsection{Ecological Diffusion}

Diffusion is a dynamic process that can be used to describe the movement
of particles, animals, pathogens, or other objects from a Lagrangian
perspective (\citealt{hooten2017animal}). A specific type of diffusion,
called ecological diffusion (\citealt{hooten2017animal}), has been
used to explain the spread of disease (\citealt{garlick2011homogenization,garlick2014homogenization,hefley2017spatial,hefley2017mechanism}),
the outbreak of pests (\citealt{powell2014phenology,hooten2013computationally,powell2017differential}),
the invasions of species (\citealt{neupane2015invasion}), and the
reintroduction of extirpated species (\citealt{williams2017integrated}).
Ecological diffusion is expressed by the PDE
\begin{equation}
\frac{\partial}{\partial t}u(\mathbf{s},t)=\left(\frac{\partial^{2}}{\partial s_{1}^{2}}+\frac{\partial^{2}}{\partial s_{2}^{2}}\right)[\mu(\mathbf{s},t)u(\mathbf{s},t)]
\end{equation}

\noindent where $u(\mathbf{s},t)$ is the intensity of the dispersing
pathogen (or population), $s_{1}$ and $s_{2}$ are the spatial coordinates
contained in the vector $\mathbf{s}$, and $t$ is the time. The intensity,
$u(\mathbf{s},t)$, can be linked to the expected pathogen abundance
by integrating $u(\mathbf{s},t)$ over a spatial domain of interest
(i.e., $\int_{\mathcal{A}}u(\mathbf{s},t)d\mathbf{s}$). The diffusion
coefficient, $\mu(\mathbf{s},t)$ could depend on location-specific
covariates that control the rate of spread and vary over time. Although
equation (1) only describes ecological diffusion, PDEs can be modified
to capture other important components such a population growth or
long-range dispersal (\citealt{holmes1994partial}). For example,
in our WNS data example we modify equation (1) to include an exponential
growth component to enable the pathogen $\textit{P. destructans}$
to increase in abundance.

Ecological diffusion describes how the spatial distribution of a pathogen
or population evolves over time. To initiate this dynamic process,
the conditions at the time of introduction, $t_{0}$, must be known
or estimated. In situations of disease outbreak neither $t_{0}$ or
$u(\mathbf{s},t_{0})$ are known. A realistic assumption is to presume
that the pathogen was introduced at points, which can be formulated
mathematically as
\begin{equation}
u(\mathbf{s},t_{0})=\begin{cases}
\theta_{j} & \text{if}\;\mathbf{s}=\boldsymbol{\omega}_{j}\\
0 & \text{if}\;\mathbf{s}\neq\boldsymbol{\omega}_{j}
\end{cases}\:,
\end{equation}

\noindent where $\theta_{j}$ is the unknown initial number of pathogen
particles (or animals) and $\boldsymbol{\omega}_{j}$ is the unknown
coordinate of the $j^{\text{th}}$ location of introduction ($j=1,2,\ldots,J$).
In data-driven applications, the time, location, initial number of
particles, and potentially the number of points are unknown, thus
the parameters $t_{0}$, $\boldsymbol{\omega}_{j}$, $\theta_{j}$,
and $J$ must be estimated.

Estimation of the initial conditions can be understood heuristically
as follows. By observing data generated from a diffusion process at
multiple time points (Fig. 2), it is feasible to estimate the rate
at which the process is diffusing (spreading). Once the diffusion
rate is estimated, equation (1) can predict the future (forecast)
or past (backcast) spread of the pathogen. Unlike forecasting, backcasting
ultimately results in a highly constrained state with a functional
form that we have knowledge about. Specifically, if the diffusion
processes started at points, then backcasting has a known endpoint
that constrains the functional form of $u(\mathbf{s},t_{0})$ to equation
(2). This constraint, formulated as a prior or process-level distribution
when using a hierarchical Bayesian statistical model, enables the
posterior distribution of the time, location, and number of points
to be calculated from data. Next we introduce the data set that motivates
our work. Understanding the specifics of the data is required to specify
an appropriate hierarchical Bayesian model that links the ecological
diffusion PDE to observed data and enables us to estimate the date
and location of introduction.
\begin{figure}[H]
\begin{centering}
\includegraphics[scale=0.6]{Fig_2}
\par\end{centering}
\noindent \textbf{Fig. 2.} Illustration showing an initial introduction
of $\theta=100$ particles and how the one-dimensional spatial distribution
of pathogen or population intensity ($u(s,t)$) evolves over time
according to the ecological diffusion partial differential equation.
Greater values of the diffusion rate $\mu(s)$, shown in the inset
plot, result in faster spread whereas spread is inhibited when $\mu(s)=0$.
The diffusion rate depends on the environmental variables at the location
($s$), which results in the spatial variability (roughness) visible
in $u(s,t)$. The pathogen intensity, $u(s,t)$, can be linked to
the expected pathogen abundance by integrating $u(s,t)$ over an interval
of interest (i.e., $\int_{a}^{b}u(s,t)ds$). In this example, $\int_{-0.5}^{0.5}u(s,t)ds=100$
for all time points.
\end{figure}

\section{White-nose Syndrome Data Example}

Surveillance for WNS in the United States began in 2007 using a combination
of passive and active surveillance methods. During 2007\textendash 2012,
samples were obtained from individual bats associated with morbidity
or mortality investigations occurring year-round at underground hibernacula
or on the above-ground landscape, and consisted of bat carcasses,
biopsies of wing skin, or tape lifts of fungal growth on muzzles.
A small number of samples were also obtained from target species (including
$\textit{Myotis}$ spp., $\textit{Perimyotis subflavus}$, $\textit{Eptesicus fuscus}$)
admitted to rehabilitation facilities or state diagnostic laboratories
for rabies testing from approximately December to May. A positive
or negative diagnosis of WNS in individual bats was determined by
observing characteristic histopathologic lesions in skin tissues using
light microscopy (\citealt{meteyer2009histopathologic}). A diagnosis
of suspect WNS was assigned to individuals with clinical signs suggestive
of the disease that either had ambiguous skin histopathology or were
not evaluated in this manner but for which the causative agent, $\textit{P. destructans}$,
was detected by fungal culture or polymerase chain reaction (\citealt{lorch2010rapid}).

We illustrate our modeling approach using a subset of the WNS surveillance
data collected during 2008\textendash 2011 that includes individual
samples of little brown bats (\textit{Myotis lucifugus}), big brown
bats (\textit{Eptesicus fuscus}), northern long-eared bats (\textit{Myotis
septentrionalis}), and tri-colored bats (\textit{Perimyotis subflavus}).
This resulted in $\geq50$ individual samples for each species and
a total of 397 samples with 184 positive or suspected positive cases
of WNS (Fig. 1). We limited our illustration to the four most common
species sampled because the remaining species were uncommon and included
less than 20 individuals and no more than five positive cases for
each species. In our analysis, we chose to exclude the first case
detected in 2006 because it was based on photographic evidence as
well as samples from 2007 because diagnostic case criteria for WNS
were not established until 2008. Finally, we used data collected from
2012 for the same four species to test the ability of our model to
forecast future spread of $\textit{P. destructans}$ and prevalence
of WNS. This resulted in $\geq3$ individuals for each species and
a total of 54 individuals with 37 positive or suspected positive cases.

\subsection{Hierarchical Bayesian Statistical Model}

We take a hierarchical Bayesian approach and specify a statistical
model that facilitates simultaneous parameter estimation of the location
and date of introduction as well as the diffusion and growth rate
(\citealt{hooten2013computationally}). Our model is tailored to the
WNS data (Fig. 1), but the hierarchical Bayesian approach is flexible
and can be adapted to a wide range of studies (e.g., \citealt{Hobbs2015,williams2017integrated}).
For the WNS data, we specify the hierarchical model:\begin{eqnarray}
y_{i} & \sim & \text{Bernoulli}\left(p_{i}\right)\\
p_{i} & = & g^{-1}(u(\mathbf{s}_{i},t_{i})e^{\mathbf{x}^{\prime}_{i}\boldsymbol{\beta}})\\
\frac{\partial}{\partial t}u(\mathbf{s},t) & = & \left(\frac{\partial^{2}}{\partial s_{1}^{2}}+\frac{\partial^{2}}{\partial s_{2}^{2}}\right)[\mu(\mathbf{s})u(\mathbf{s},t)]+\lambda(\mathbf{s})u(\mathbf{s},t)\\
\text{log}\left(\mu(\mathbf{s})\right) & = & \alpha_{0}+\mathbf{z}(\mathbf{s})^{\prime}\mathbf{\boldsymbol{\alpha}}\\
\lambda(\mathbf{s}) & = & \gamma_{0}+\mathbf{w}(\mathbf{s})^{\prime}\mathbf{\boldsymbol{\gamma}},
\end{eqnarray}where $y_{i}$ is equal to $1$ if the $i^{\text{th}}$ tested bat
is WNS-positive or suspected positive and $0$ if negative. The probability
that $y_{i}=1$, $p_{i}$, depends on the pathogen intensity $u(\mathbf{s}_{i},t_{i})$
defined by the PDE in equation (5) and species susceptibility $e^{\mathbf{x}_{i}^{\prime}\boldsymbol{\beta}}$,
where $\mathbf{x}_{i}$ is a vector that contains a binary indicator
variable identifying the species of the $i^{\text{th}}$ tested bat
and $\boldsymbol{\beta}\equiv(\beta_{1},\ldots,\beta_{p_{\beta}})^{'}$
is a vector of species-specific susceptibility coefficients. Differences
among species susceptibility to WNS has been reported although the
mechanism(s) for these differences remain unclear (\citealt{langwig2012sociality,johnson2015antibodies,langwig2015invasion}).
Regardless of the mechanism(s), our model accounts for heterogeneity
among species susceptibility by inclusions of $e^{\mathbf{x}_{i}^{\prime}\boldsymbol{\beta}}$.

For the ecological diffusion PDE, the product of the pathogen intensity
($u(\mathbf{s}_{i},t_{i})$) and species susceptibility ($e^{\mathbf{x}_{i}^{\prime}\boldsymbol{\beta}}$)
in equation (4) is greater than zero and therefore the inverse link
function $g^{-1}(\cdot)$ maps the positive real line $[0,\infty)$
to a probability between 0 and 1. For an inverse link function, we
used the cumulative distribution function of a standard log-normal
distribution (i.e., $g^{-1}(x)=\frac{1}{x\sqrt{2\pi}}\int_{0}^{x}e^{-\frac{1}{2}\text{log}(x)^{2}}dx$).
Although any appropriate link function could be used, our choice enables
efficient implementation because the full-conditional distribution
for $\boldsymbol{\beta}$ are available in closed-form (see Appendix
S1 for details; \citealt{BB2L}). The diffusion rate ($\mu(\mathbf{s}$))
in equation (5) depends on the regression-type equation in (6) that
has an intercept term ($\alpha_{0}$), coefficients ($\boldsymbol{\alpha}$$\equiv(\alpha_{1},\ldots,\alpha_{p_{\alpha}})^{'}$),
and location-specific covariates $\mathbf{z}(\mathbf{s})$. For spreading
diseases, the diffusion rate is always positive, thus $\mu(\mathbf{s})>0$
for all $\mathbf{s}$, which motivates the log link function in (6).
In addition to ecological diffusion, we add an exponential growth
component to the PDE to account for increasing (or decreasing) pathogen
intensity. This addition results in a growth rate ($\lambda(\mathbf{s})$)
that depends on a similar regression-type equation (7) that has an
intercept term ($\gamma_{0}$), coefficients ($\boldsymbol{\gamma}\equiv(\gamma_{1},\ldots,\gamma_{p_{\gamma}})^{\prime}$),
and location-specific covariates $\mathbf{w}(\mathbf{s})$, which
can be the same as or different from $\mathbf{z}(\mathbf{s})$. To
allow the growth rate to be positive or negative, we used the identity
link function so that the support of $\lambda(\mathbf{s})$ is unconstrained
and can take on values from $-\infty$ to $\infty$.

We assume Dirichlet boundary conditions and set $u(\mathbf{s},t)=0$
at the boundary of the United States. Similar to initial conditions,
boundary conditions could be an important component and future studies
may need to estimate the boundary conditions. For computational efficiency,
we prefer Dirichlet boundary conditions because existing techniques
for estimation would require estimating a large number of parameters
(\citealt{wikle2003hierarchical}).

To fully specify a hierarchical Bayesian model, prior distributions
must be carefully chosen for all parameters. We describe and justify
our choice of prior distributions in Appendix S2, with the exception
of the prior distribution for the location of introduction, which
we describe below. We specify the initial conditions as described
in equation (2). The natural choice for the prior of the latent (unobserved)
locations of introduction is a point process, which is a probability
distribution that generates locations in geographic space as random
variables (\citealt{cressie2011statistics}). For the WNS data, we
assume a single point source introduction (i.e., $J=1$), which results
in a conditional point process (i.e., a point process with a known
number of points, in this case just one; \citealt{cressie1993statistics}
p. 651). For an introduction that occurs at a single location, the
intensity function of the point process must be specified as a hyperparameter;
we specify the intensity function of the point process to be constant
across the United States which is equivalent to assuming that $\textit{P. destructans}$
was equally likely to have been introduced at any location within
the United States.

If multiple introductions occur, it may be feasible to estimate the
intensity function of the point process (Appendix S1). For example,
the intensity function may be modeled using a regression-type equation
to estimate how location-specific risk factors influence the likelihood
of novel introductions. This approach could then be used to predict
the geographic areas and identify the environmental conditions where
new introductions are likely to occur. Finally, our approach can be
modified to handle an unknown number of point sources that are introduced
on different dates by using a spatio-temporal Poisson point process
(Appendix S1). This modification may eventually be used to explain
the recent positive cases of WNS in the state of Washington (\citealt{lorch2016first}).

\subsection{Model Implementation}

The covariates $\mathbf{z}(\mathbf{s})$ and $\mathbf{w}(\mathbf{s})$,
which influence the diffusion and growth rate in equations (6) and
(7), respectively, can include any location-specific environmental
factor (e.g., temperature). We used karst terrain and deciduous forest
cover calculated from digitized maps (\citealt{tobin2004digital})
and the 2011 National Land Cover Database (\citealt{NLCD2011}), respectively
(Fig. 3). We expect that the spatial availability of caves and deciduous
forest habitat will influence the diffusion of the pathogen due to
the ability of bats to disperse; however, we acknowledge many other
factors could be included depending on the goals of the study.
\begin{figure}[H]
\begin{centering}
\includegraphics{Fig_3}
\par\end{centering}
\noindent \textbf{Fig. 3.} The Environmental covariates proportion
of deciduous forest (top panel) and karst terrain (bottom panel).
\end{figure}

For model implementation, we first use an analytical technique called
homogenization that results in a computationally efficient approximation
to the solution of the ecological diffusion PDE (\citealt{garlick2011homogenization,hooten2013computationally}).
To solve the homogenized PDE, we use a standard finite-difference
scheme and chose 30 time steps per month, broad-scale grid cells of
100 km by 100 km, and the fine-scale grid cells of 10 km by 10 km
(\citealt{Farlow1993PDE,garlick2011homogenization,hooten2013computationally}).
Finer spatial discretization may be implemented, but the resolution
we chose is a reasonable trade-off between the location accuracy of
our data (mostly reported at the county-level) and what is computationally
feasible using accessible high-performance computing. For implementation
of our Bayesian model, we used a Markov chain Monte Carlo (MCMC) algorithm
with three chains. For each chain, we obtained 320,000 samples from
the posterior distribution and discarded the first 20,000 samples
to ensure that the MCMC algorithm was sampling from the target (stationary)
distribution. All computations were conducted using a program in R
(Appendix S1 and S2) with optimized basic linear algebra subprograms
and subroutines written in C++ and called from R (\citealt{TeamR}).

\subsection{Model and Software Validation}

As noted by \citet{cook2006validation}, implementing non-standard
models such as ours often requires developing the necessary software
\textquotedblleft from scratch,\textquotedblright{} which may increase
the chance of making errors when programming. Furthermore, when using
a complex statistical model it is important to check that the parameters
are identifiable (\citealt{lele2010model}). One way to check our
software for errors and our proposed model for identifiability of
parameters is to use a simulated data set that mimics the WNS surveillance
data, but where the parameter values are known. If our software and
model are able to recover the known values of the parameters from
the simulated data that mimics the WNS surveillance data, then there
is evidence that the software is error free and that the parameters
in our model are identifiable. Appendix S1 contains a simulated data
example to demonstrate that our model is capable of recovering known
values of the parameters (see Figs. S1 \& S2).

In addition to a single simulated data set that mimics the WNS surveillance
data, we conducted a simulation experiment that involved simulating
several thousand data sets under different settings and fitting our
proposed model to each data set. A simulation experiment can be used
to evaluate frequentist properties of the estimated parameters, such
as bias and coverage probability, which ensures statistical inference
from our model is valid (\citealt{little2006calibrated}). We conducted
a simulation experiment employing a simplified version of the model
we used for WNS surveillance data, but included a setting for multiple
introductions (see Table S3 and Figs. S6\textendash S7 for simulation
settings in Appendix S3). For each setting, we evaluated bias of the
estimated time and location as well as the coverage probability of
the 90\% credible intervals (or regions) for each parameter. We expect,
unless something is pathologically wrong with our model (e.g., unidentifiable
parameters; \citealt{gustafson2015bayesian}), that point estimates
will be unbiased and the 90\% credible intervals will cover the true
value of the parameters with a probability of 0.90.

For all settings, our simulation experiment shows that point estimates
of the time and location of introduction are apparently unbiased (Figs
S8\textendash S9 in Appendix S3) and the 90\% credible intervals cover
the true location and time with a probability of $\geq0.86$. Appendix
S3 contains details of our simulation experiment, full results, and
computer code needed to reproduce our experiment.

\subsection{Forecast Comparison}

Accurate forecasts during the early stages of disease outbreak or
invasion of a novel species may be useful to inform management recommendations
and surveillance sampling. In addition to estimating the location
and date of introduction, an auxiliary benefit of Bayesian implementations
of dynamic models is that probabilistic predictions and forecasts
can be obtained as a derived quantity (\citealt{Hobbs2015,hefley2017mechanism}).
To demonstrate the accuracy of forecasts from our model, we used the
samples collected from 2012 that were reserved for this single purpose
(see WNS data description for more detail). To facilitate comparison
of our model with other non-Bayesian approaches, we calculated the
zero-one score using the posterior mean of the probability of infection
(i.e., $\text{E}(p_{i}|\mathbf{y})$) and forecasted that an individual
sample would be WNS ``positive'' if $\text{E}(p_{i}|\mathbf{y})\geq0.5$
and ``negative'' if $\text{E}(p_{i}|\mathbf{y})<0.5$ (\citealt{gneiting2007strictly,gneiting2011making}).
Although any local and proper scoring rule could be used, we chose
the zero-one score because it can be used to calculate the misclassification
rate, which is easy to interpret and communicate to practitioners.
The misclassification rate is calculated using the 2012 data by summing
the total number of samples incorrectly classified and dividing by
the total number of samples.

As a standard of comparison, we followed \citet{hefley2017mechanism}
and compared the forecasting ability of our model to phenomenological
regression-type models developed for classification of binary events
(e.g., positive/negative) from the statistics and machine learning
literature (\citealt{james2013introduction}). Specifically, we compared
forecasts from our model to three approaches, which include: (1) a
generalized linear model (GLM) that included deciduous forest cover
and karst terrain at the sample location and bat species as predictor
variables; (2) a generalized additive model (GAM) that included a
linear effect of the same predictor variables as the GLM, but that
accounted for spatial autocorrelation by using two-dimensional splines
on a sphere (\citealt{wood2017generalized}); and (3) the same GAM
specification as in 2, but with an additional smooth effect of time
to account for increasing prevalence and temporal autocorrelation.
Details associated with the exact model specifications and computations
are reported in Appendix S2.

For each approach, we used the samples collected from 2012 and calculated
the zero-one score using the maximum (or restricted maximum) likelihood
estimate of the predicted probability of infection ($\hat{p}_{i}$)
and forecasted that an individual sample would be WNS ``positive''
if $\hat{p}_{i}\geq0.5$ and ``negative'' if $\hat{p}_{i}<0.5$.
To compare forecasts with our model, we report the misclassification
rate for each approach. Although comparing forecasts among different
models is useful for quantifying forecast accuracy, this comparison
is does not explicitly evaluate the reliability of the inferred time
and location of introduction.

\section{Results}

The year with the highest probability of introduction was 2002 (1992\textendash 2005;
90\% credible interval), which is 4 years before the first photographic
evidence documented the presence of $\textit{P. destructans}$ in
the United States (Fig. 4). When compared to the location of first
detection based on photographic evidence from 2006, the posterior
distribution shows an area covering $123,798$ km$^{2}$ that has
an equal or higher probability of containing the point where $\textit{P. destructans}$
was introduced (Fig. 4). This $123,798$ km$^{2}$ area represents
a 74\% credible region and intersects the location of first detection,
but indicates that location of introduction that are equally or more
probable could be up to 510 km east of the first detection. For comparison,
the 90\% posterior credible region covers an area of $182,493$ km$^{2}$.
This $182,493$ km$^{2}$ regions represents 2.3\% of the area covered
by the 90\% prior credible region, which demonstrates a large reduction
in uncertainty due to the WNS surveillance data.
\begin{figure}[H]
\begin{centering}
\includegraphics[scale=0.65]{Fig_4}
\par\end{centering}
\noindent \textbf{Fig. 4.} Posterior distributions of the location
and date of introduction obtained from fitting our model to the data
shown in Fig. 1. The top panel shows the joint posterior distribution
of the location of introduction with red indicating a higher probability
than yellow. The beige polygon outlines an area that is equally or
more likely to contain the point of introduction when compared to
the location where white-nose syndrome was first detected in 2006
(blue dot). The gray areas are karst formations with caves, which
is likely winter bat habitat. The bottom panel shows the marginal
posterior distribution of the date of introduction. The blue vertical
line represents the earliest known photographic evidence of white-nose
syndrome on bats in the United States. Years with low probability
are not shown (i.e., before 1980).
\end{figure}

In addition to estimating the location and date of introduction, our
approach simultaneously estimates species-specific susceptibility
(Fig. 5), diffusion and growth rates of the pathogen (Fig. 6), and
can be used to forecast the spread of $\textit{P. destructans}$ (Fig.
7 and animation in Appendix S4). As it relates to species-specific
susceptibility, our results show that the northern long-eared bat
was the most susceptible to $\textit{P. destructans}$ whereas the
big brown bat was the least susceptible (Fig. 5). The diffusion and
growth rates are positive in many regions and show spatial variability
that is attributable to deciduous forest cover and karst terrain that
contains caves (cf., Fig. 3 and Fig. 6, Fig. S3). As a result of positive
growth and diffusion rates, forecasts show that $\textit{P. destructans}$
spread rapidly across the contiguous United States (Fig. 7), causing
an overall increase in the probability of infection for all four species
of bats over time (see animations in Appendix S4).

Our forecast comparison revealed a misclassification rate of 26\%
for our hierarchical Bayesian model. For comparisons, the three alternative
approaches we used resulted in misclassification rate of 43\% for
method one (GLM), 28\% for method 2 (GAM that included karst, forest
and species as predictors and that accounted for spatial autocorrelation)
and 30\% for method three (same as method 2, but accounted for temporal
dynamics and autocorrelation). A visual comparison of the forecasts
for all methods is given in Appendix S2 (Fig. S5).
\begin{figure}[H]
\begin{centering}
\includegraphics[scale=0.55]{Fig_5}
\par\end{centering}
\noindent \textbf{Fig. 5.} The estimated species-specific susceptibility
depends on the intensity of $\textit{P. destructans}$ ($u(\mathbf{s},t)$).
Shown is the expected value of the posterior distribution (black lines)
and 95\% credible intervals (gray shading).
\end{figure}
\begin{figure}[H]
\begin{centering}
\includegraphics[scale=0.9]{Fig_6}
\par\end{centering}
\noindent \textbf{Fig. 6.} The estimated diffusion rate (top panel)
determines how fast the pathogen $\textit{P. destructans}$ will spread,
whereas the growth rate (bottom panel) determines how quickly the
intensity of $\textit{P. destructans}$ will increase. Deciduous forest
and karst terrain (Fig. 3) influence the posterior distributions of
the diffusion and growth rate leading to the spatial variability appearing
in the maps.
\end{figure}
\begin{figure}[H]
\begin{centering}
\includegraphics[scale=0.7]{Fig_7}
\par\end{centering}
\noindent \textbf{Fig. 7.} Estimated intensity of $\textit{P. destructans}$
($u(\mathbf{s},t)$) for January 2002\textendash 2012. See Appendix
S4 for animations showing the contiguous United States for all months
over the period 2000\textendash 2012 as well as forecasts of the probability
of infection for all four species of bats. Note that the pathogen
intensity, $u(\mathbf{s},t)$, can be linked to the expected pathogen
abundance by integrating $u(\mathbf{s},t)$ over a spatial domain
of interest (i.e., $\int_{\mathcal{A}}u(\mathbf{s},t)d\mathbf{s}$).
\end{figure}

\section{Discussion}

Our results reveal that it is possible to estimate the date and location
of introduction shortly after the release of a pathogen using indirect
data from 397 samples collected over 4 years. Our results also demonstrate
the need for statistical thinking because the estimated date shows
that the introduction likely occurred several years before the first
detected case. Conversely our results show that the location of first
detection is a plausible location of introduction, but there is a
moderate level of uncertainty surrounding the exact location.

For WNS, the estimated location and date of introduction, along with
the associated uncertainty, constitutes knowledge generated by harnessing
the information within disease surveillance data using statistically
valid techniques. Similar surveillance data is routinely collected
for diseases of wildlife, plants, livestock, and humans, but filling
knowledge gaps using these data requires statistical methods that
capture dynamic processes and quantify uncertainty.

The introduction of $\textit{P. destructans}$ is suspected to have
originated from a single point source due to the clonal nature of
the fungal pathogen in North America (\citealt{drees2017phylogenetics}).
For other pathogens or pests, introductions may occur at multiple
locations. When multiple introductions occur, we explain in Appendix
S1 how our approach can be modified accordingly and potential issues
that may arise when fitting the model to data. When multiple introduction
occur, our approach could identify areas and environmental conditions
where new pathogens or pest introductions are likely to occur\textendash effectively
providing a map that is similar to Fig. 4, but showing areas at high
risk for new introductions rather than a single past introduction.
With this information, new introductions could be anticipated and
perhaps avoided or optimal monitoring and control programs could be
established (\citealt{williams2017monitoring}; \citealt{laber2018optimal}).

Although our study system was epidemiological, dynamic statistical
models have been used to explain the outbreak of pests (\citealt{powell2014phenology,hooten2013computationally,powell2017differential}),
the invasions of a species (\citealt{neupane2015invasion}), the reintroduction
of an extirpated species (\citealt{williams2017integrated}), and
to determine optimal monitoring of ecological processes (\citealt{williams2017monitoring}).
Such a range of applications indicate that dynamic statistical models
can be used for a variety of situations including to determine the
location of remnant populations that enable outbreaks of pests, estimate
the date, location, and propagule size where invasive species were
introduced, or determine the habitat characteristics where an extirpated
species can re-establish.

\section*{\vspace{-0.5cm}
Acknowledgements}

We thank all state, federal and other partners for submitting samples
and the U.S. Geological Survey National Wildlife Health Center for
processing the samples. At the time of publication data were not available
for release because of the sensitivity related to the locations of
the bat species sampled. Data can be requested from the following
agencies that submitted samples. Agencies that submitted samples include
Acadia National Park, Aerospace Testing Alliance Conservation, Alabama
Department of Conservation and Natural Resources, Alabama Division
of Wildlife and Freshwater Fisheries, Alaska Department of Fish and
Game, U.S. Fish and Wildlife Service (USFWS) Alligator River National
Wildlife Refuge, Arizona Game and Fish Department, Arkansas Department
of Health, Arkansas Game and Fish Commission, U.S. Air Force, Avon
Park Air Force Range, Bat Conservation International, Bat World Sanctuary,
U.S. Department of Agriculture (USDA) Forest Service, Black Hills
National Forest, Bureau of Land Management, Boston University, Buffalo
National River (National Park Service), California Department of Fish
and Game, Canadian Cooperative Wildlife Health Centre, Channel Islands
National Park, Clinton Wildlife Management Area, Colorado Department
of Natural Resources, USDA Forest Service Colville National Forest,
Connecticut Department of Environmental Protection, Delaware Department
of Natural Resources and Environmental Control, Department of the
Army, Mississippi Department of Wildlife, Fisheries and Parks, Devils
Tower National Monument, East Stroudsburg University, Ecological Sciences,
Inc., Florida Fish and Wildlife Conservation Commission, Frenso Chaffee
Zoo, Furman University, Georgia Department of Natural Resources, USDA
Forest Service Gifford Pinchot National Forest, Great Smoky Mountains
National Park, USDA Forest Service Hoosier National Forest, Humboldt
State University, Idaho Department of Fish and Game, Illinois Department
of Natural Resources, Illinois Natural History Survey, Indiana State
University, Iowa Department of Natural Resources, U.S. Air Force,
Langley, Kansas Department of Parks, Wildlife, and Tourism, Lake of
the Ozarks State Park, Lake Vermilion-Soudan Underground Mine State
Park, Lava Beds National Monument, USDA Forest Service, Lincoln National
Forest, Louisiana Department Baton Rouge, Maine Department of Inland
Fisheries and Wildlife, Maryland Department of Agriculture, Maryland
Department of Natural Resources, Massachusetts Division of Fish and
Wildlife, Michigan Department of Natural Resources, Ministere des
Ressources naturelles et de la Faune, Minnesota Department of Natural
Resources, Missouri Department of Conservation, Missouri State University,
Montana Department of Fish, Wildlife and Parks, Nature Conservancy,
Nebraska Game and Parks Commission, Nevada Department of Wildlife,
New Hampshire Fish and Game Commission, New York Department of Environmental
Conservation, North Carolina Wildlife Resources Commission, North
Dakota Game and Fish, North Mississippi Refuges Complex, USFWS Noxubee
National Wildlife Refuge, EL Malpais National Monument, Ohio Department
of Natural Resources, Ohio Wildlife Center, Oklahoma State University,
Oklahoma Department of Wildlife Conservation, Olentangy Wildlife Research
Station, Oregon Caves National Monument, Oregon Department of Fish
and Wildlife, Oregon Parks and Recreation Department, Ozark National
Scenic Riverway, Ozark Plateau National Wildlife Refuge, Pennsylvania
Game Commission, Pennsylvania State Animal Diagnostics Laboratory,
Pima County Natural Resources, Parks, and Recreation, Rhode Island
Department of Environmental Management, Rhode Island Division of Fish
and Wildlife, USFWS San Bernardino/Leslie Canyon National Wildlife
Refuge, South Dakota Department of Game, Fish and Parks, Session Woods
Wildlife Management Area, USDA Forest Service, Shawnee National Forest,
Sleeping Bear Dunes National Lakeshore, South Carolina Department
of Natural Resources, South Mountain Wildlife Rehabilitation Center,
Tennessee Wildlife Resources Agency, Texas Parakeet Raisers, Texas
Parks and Wildlife, University of Central Oklahoma, University of
Montreal, University of Nebraska, USFWS Upper Mississippi River National
Wildlife and Fish Refuge, USDA Forest Service National Forests and
Grasslands in Texas, USDA Forest Service Hot Springs Arkansas, USDA
Animal and Plant Health Inspection Service Wildlife Services, USDA
Forest Service North Central Research Station, USDA Forest Service
Pacific Southwest Research Station, USFWS Arkansas Field Office, USFWS
Ecological Services, USFWS Environmental Contaminants, Utah Division
of Wildlife Resources, Vermont Fish and Wildlife, Virginia Department
of Conservation and Recreation, Virginia Division of Natural Heritage,
Virginia Game and Inland Fisheries, Washington Department of Fish
and Wildlife, West Virginia Department of Natural Resources, USFWS
Wheeler National Wildlife Refuge, USDA Forest Service White River
National Forest, Wisconsin Dept of Natural Resources, Wisconsin State
Lab of Hygiene, Wyoming Fish and Game. We thank Cydney Alexis, Nora
Bello, Hannah Birg$\acute{\text{e}}$, Jeffrey Lorch, Bob O'Hara,
Daniel Walsh, and Perry Williams for insightful discussions and comments
that improved this manuscript. The authors acknowledge support for
this research from USGS G16AC00413. Any use of trade, firm, or product
names is for descriptive purposes only and does not imply endorsement
by the U.S. Government.\vspace{-0.5cm}

\section*{Supplementary Material}

\noindent \textbf{Appendix S1}

\noindent Tutorial with R code to simulated data and fit the Bayesian
ecological diffusion model.

\noindent \textbf{Appendix S2}

\noindent Tutorial and R code to reproduce the white-nose syndrome
data analysis.

\noindent \textbf{Appendix S3}

\noindent Simulation experiment description, results, and R code.

\noindent \textbf{Appendix S4}

\noindent Supplementary animations of Fig. 7.

\vspace{-0.5cm}

\begin{singlespace}
\noindent \bibliographystyle{apa}


\section{Introduction}

The outbreak of COVID-19 is disrupting the personal lives of people around the world in a variety of ways and has made profound impacts on the world socially and economically.  As the number of confirmed COVID-19 cases in the U.S. continues to rise in early 2020, many states declared states of emergency and shutdown orders or stay-at-home orders were issued to slow the spread of COVID-19. Many schools, workplaces, and public gathering spaces across the U.S. were closed for an extended time. In June 2020, as states gradually reopen the economy, the number of COVID-19 cases starts to increase again in many states.

In June 2020, Arizona emerged as one of the country's newest COVID-19 hot spots with an alarming resurgence of COVID-19 after the stay-at-home order was lifted in the middle of May. The number of people hospitalized is also climbing. On June 15, Arizona reported 2,497 new COVID-19 cases, a new all-time high record, up from the previous highest number of cases of 2,140 cases in a single day on June 12. Previously, Arizona experienced the first rapid rise in the number of confirmed cases of COVID-19 infection on March 16 with a total of 1,157 confirmed cases including 20 deaths, according to the Arizona Department of Health Services (\cite{ADHS2020}). The spike of the new COVID-19 case is believed to result from the fact that Arizona took a more aggressive approach to reopen the state in the middle of May. In addition, there is also no state-level requirement for everyone to wear face masks. Many Arizonans have different views on practice social distancing which is the only way to prevent infections in the absence of a vaccine. Although under Arizona's reopening plan, businesses are only advised to follow federal guidance on social distancing, the trade-offs between health and the economy have forced many businesses to choose not fully implementing the social distancing requirements during the COVID-19 pandemic. Therefore, quantifying the impact of social distancing on minimizing the epidemic impact can make governments and convince more people to better understand the significance of social distancing and personal precautions. In the present paper, we examine the COVID-19 situation in the state of Arizona in the United States and predict the number of cases in Arizona after the reopening of the economy.

There have been plenty of mathematical models to describe the spread of infectious diseases. To reply to the current COVID-19 pandemic, many modeling studies using ordinary differential equations (ODE) have proposed for predicting the COVID-19 cases (e. g. \cite{He2020MBE,Li2020MBE,Wang2020MBE}). The classical susceptible infectious recovered model (SIR)  (\cite{Huang2020BullWHO}) and susceptible exposed infectious recovered model (SEIR) (\cite{Lai2020Nature, Omori2020MedRxiv, Yang2020JTD}) are the most widely adopted ones for characterizing the outbreak of COVID-19. The extension of the classical SEIR model with age-stratified model (\cite{Prem2020LancetPH}) and meta-population model (\cite{Pujari2020medRxiv}) were also introduced.

In this paper, we present a spatio-temporal model, specifically, a partial differential equation (PDE) model, of COVID-19 based on county-level clusters.  The proposed model describes the combined effects of transboundary spread among county clusters in Arizona and human activities on the transmission of COVID-19. This is, in particular, important for Arizona as well as the other states in the U.S. because new COVID-19 cases in the U.S. are still continuing to head upward after the lockdown ended while most countries had peaked their new cases after lockdowns. One of the factors contributing to the continued spike of the U.S. cases is that personal precautions are not consistent across the U.S. once the lockdown ends. In Arizona, the face mask requirement is implemented at the local government level. A number of counties in Arizona do not require face masks. We divided the counties in Arizona into three regions. The three regions represent different geographical and social characteristics in terms of the spread of COVID-19. These disparities, which reflect differing regional levels of risk, can be better modeled by a spatio-temporal model. As a result, the characteristics of the PDE model make it suitable for COVID-19 prediction. The localized results of the spatio-temporal model provide more useful information for the local governments to closely monitor new COVID-19 clusters and quickly reinstating lockdowns at the local level when epidemic spikes.

Many spatio-temporal models characterize infectious diseases by PDEs, describing the dynamics of susceptible, infected, and recovered populations (e.g. \cite{Brauer2019book, holmes1994ecology, Wang2020mdpi, Zhu2017MBE}). Our previous work (\cite{Haiyan2019book,wang2016regional}, \cite{Wang2020mdpi}) applies PDE models to make a regional level of influenza with geo-tagged data.  The PDE model we develop in the present work focuses on the spread of COVID-19 and incorporating social distancing factors.  While there is a rich literature on the application of PDE to modeling the spatial spread of infectious diseases, to the best of our knowledge, this work is the first attempt to apply PDE models on COVID-19 prediction with real-data validation.

In addition, our model is different from the existing models in that we include the data from the Google Community Mobility Report  (\cite{Google2020Mobility}) which reflects the effects of human activities. To help to combat the spread of COVID-19, Google releases the COVID-19 Community Mobility Reports which provide daily, county-level aggregated data on time spent at different categories of activities, compared with a baseline period before the epidemic.  These activities include retail and recreation, groceries and pharmacies, parks, transit stations, workplaces, and residential areas. Mobility trends obtained from location history are dynamic in time and reflect real-time changes in social behavior, making them a crucial factor in analyzing COVID-19 spread and its countermeasure. There are several ODE and statistical models which utilize the COVID-19 Community Mobility data~(\cite{picchiotti2020arxiv,abouk2020arxiv,voko2020geroscience}). To the best of our knowledge, this is the first PDE model incorporating COVID-19 Community Mobility data to predict the number of COVID-19 cases.

\section{Data Description}
\label{data}
\begin{figure}[!h]
\centering
\includegraphics[width=0.5\textwidth]{AZplot.png}
\caption{Three regions in Arizona
  }
\label{regions}
\end{figure}
There are 15 counties in the U.S. state of Arizona. In general, Arizona is divided into Northern, Central, and Southern Arizona. Although the boundaries of the regions are not well-defined, as we aim to predict COVID-19 cases in Arizona, we group the counties of Arizona into the three regions in Fig. \ref{regions}. The three regions have better balanced COVID-19 cases because Maricopa County, where the city of Phoenix is located, has about 60\% of the population in Arizona. More importantly, they represent different geographical and social characteristics in terms of the spread of COVID-19. Much of the Northern region is National Forest Service land, parkland, and includes part of Navajo reservations. For the per-capita COVID-19 infection rate in the US, the Navajo Nation, which spans parts of Arizona, New Mexico, and Utah, is greater than that of New York City in May 2020. The Southern Arizona includes the city of Tucson and several other small cities close to the border of the U.S. and Mexico.  In the following, we will use NAZ (Region 1), CAZ (Region 2), and SAZ (Region 3) to denote the clusters of the Arizona counties.

\begin{itemize}
\item  NAZ (Region 1)=[Mohave County, Coconino County, Yavapai County, La Paz County, Gila County, Navajo County, Apache County].
\item  CAZ (Region 2)=[Maricopa County].
\item  SAZ (Region 3)=[Yuma County, Pima County, Pinal County, Graham County, Santa Cruz County, Greenlee County, Cochise County].
\end{itemize}

We compute the daily new cases of each region by adding the new COVID-19 cases of all counties belonging to a region. We use the COVID-19 data from the New York Times at the state and the county level over time. The New York Times is compiling the time series data from state and local governments and health departments in an attempt to provide a complete record of the ongoing outbreak. The data can be downloaded from https://github.com/nytimes/covid-19-data and https://www.nytimes.com/article/coronavirus-county-data-us.html.

In addition, we incorporate the Google Community Mobility Report into our model. These Community Mobility Reports provide insights into how people's social behaviors have been changing in response to policies aimed at combating COVID-19. The reports provide the changes in movement trends compared to baselines overtime at the U.S. county level, across different categories of activities. These activities include retail and recreation, groceries and pharmacies, parks, transit stations, workplaces, and residential. The relevant data can be downloaded from https://www.google.com/covid19/mobility/. Then we compute the daily changes of each region by adding the changes of all counties belonging to a region.

We create two time series from the Google Community Mobility Report. We accumulate the data for retail \& recreation, groceries \& pharmacies, parks, transit stations, workplaces, which are believed to increase the number of COVID-19 cases. On the other hand, the data from Google mobility for residential activities describes stay-at-home activities that are believed to prevent COVID-19 epidemics. In some cases, the increase of residential activities may result from events (such as parties) that could promote COVID-19 cases. However, the Google data only indicates residential activities compared to baselines. In our model, we use a Michaelis–Menten function to limit the effect of residential activities.  Because it may take about two weeks for the effect of the activities to take place, we subtract the time $t$ in our PDE model by 14 to reflect the incubation period of COVID-19 in most cases.

\section{PDE Model}
\label{models}
In this section, we introduce a specific PDE model to characterize the spatio-temporal dynamics of COVID-19 cases. To apply a PDE model to the interaction of the dynamics of COVID-19 cases, one must embed the three regions NAZ, CAZ, SAZ into  Euclidean space and arrange them in a meaningful order.  In this paper, the three regions are mapped onto a line where the locations of NAZ, CAZ, SAZ are 1, 2, 3 respectively. The arrangement will make the three regions stay as close as possible for ensuring that the continuous model can capture the spread of COVID-19 cases between the three regions.  Projecting the three regions, geographically located in the north-south direction of Arizona, onto the x-axis of the Cartesian coordinates, as shown in Fig. \ref{embed}. We project the three regions as the order of geographical north-south direction instead of other orders mainly due to the following fact: Maricopa county has much more population than any other county in Arizona. It is more likely that COVID-19 spreads from the metro Phoenix area to both northern and southern counties, at least at the beginning of 2020.

\begin{figure}[!h]
\centering
\includegraphics[page=1,scale=0.7, trim=7cm 7cm 7cm 7cm, clip=true ]{Diffusioncluster.pdf}
\caption{Embedding of three regions into the $x$-axis and two spreading processes}
\label{embed}
\end{figure}

Having embedded the three regions to the Euclidean space, we will develop the spatio-temporal model for the spread of COVID-19 cases according to the balance law: the rate at which a given quantity changes in a given domain must equal the rate at which it flows across its
boundary plus the rate at which it is created, or destroyed, within the domain. The PDE model can be conceptually divided into two processes: an internal process within each region and an external process among regions. Similar derivation for the PDE model has been used in our previous work for PDE models for information diffusion in online social networks in \cite{Haiyan2019book}.

Let $C(x,t)$ represent the COVID-19 cases in Arizona region $x$ at a given time $t$. The changing rate of $C(x,t)$ depends on two processes as in Fig.~\ref{embed}:
\begin{itemize}
\item  Global process: the social interactions of people between regions, such as traveling and commuting between regions, that contribute to the spread of COVID-19;
\item Local process:  in each Arizona region, people become newly infected through social interactions with infected people within a region; and people may take personal precautions to reduce and mitigate COVID-19 spread.
\end{itemize}

Combining the two processes, the dynamics of COVID-19 cases can be captured by Eq.~(\ref{new_model}).
\begin{eqnarray}
\label{new_model}
\left\{\begin{array}{ll}
&\frac{\partial C(x,t)}{\partial t}=\frac{\partial }{\partial x}\big(d(x)\frac{\partial C(x,t)}{\partial x}\big)+r(t)l(x)a(x,t-14)C(x,t)-c\frac{h(x,t-14)C(x,t)}{k+C(x,t)},
\\
&C(x,1)=\psi(x),1<x<3,
\\
&\frac{\partial C}{\partial x}(1,t)=\frac{\partial C}{\partial x}(3,t)=0,t>1,
\end{array} \right.
\end{eqnarray}
where
\begin{itemize}
          \item The term $\frac{\partial}{\partial x} \left(d(x)\frac{\partial C(x,t)}{\partial x}\right)$ denotes the spread of COVID-19 cases between different Arizona regions, where  $d(x)$ measures how fast COVID-19 spreads across different Arizona regions. In epidemiology~(\cite{Brauer2019book} and \cite{natr2002murray}), the term $\frac{\partial}{\partial x} \left(d(x)\frac{\partial C(x,t)}{\partial x}\right)$ has been widely used for describing the spatial spread of infectious diseases. Here we assume $d(x)$ to be constant, i.e., $d(x)\equiv d>0$.

       \item $r(t)l(x)C(x,t)l(x)a(x,t-14)$ represents the new COVID-19 cases from a local Arizona region at location $x$ and time $t$. This type of function is widely used to describe the growth of bacteria, tumors, or social information over time~(\cite{natr2002murray}).
          \begin{itemize}
         \item[$\ast$] The function $r(t)>0$ represents the growth rate of COVID-19 cases at time $t$ for all Arizona regions. For simplicity, we assume that $r(t)$ increases with time $t$ as the COVID-19 cases increase. Therefore, we choose $r(t)=g(b_1+b_2 t)$ and $g(u)=1/(1+exp(-u))$ to describe the pattern with parameters $b_1>0, b_2>0$ to be determined by the collected COVID-19 data.
         \item[$\ast$]The location function $l(x)$ describes the spatial heterogeneity of COVID-19, which depicts different infection rates in the three Arizona regions. The function $l(x)$ is built through a cubic spline interpolation, which satisfies $l(x_i)\equiv l_i, i=1,2,3,$ where $x_i$ represents the location of Arizona region $i$.  $l(x)$ is determined by the collected COVID-19 data.
         \item[$\ast$] $a(x,t-14)$ is accumulated data from Google mobility for retail \& recreation, groceries \& pharmacies, parks, transit stations, workplaces, which are believed to increase the COVID-19 cases.  $t-14$ again reflects the incubation period of COVID-19 in most cases.
        \end{itemize}

    \item The function $c \frac{h(x,t-14)C(x,t)}{k+C(x,t)}$ is the rate of decrease of COVID-19 cases due to human efforts. During the period of the epidemic, Arizona people are advised to take COVID-19 precaution actions, such as wearing masks, staying six feet apart, or staying at home to reduce the contact rate. $c\frac{h(x,t-14)C(x,t)}{k+C(x,t)}$ to describe the potential reduction in COVID-19.
          \begin{itemize}
          \item[$\ast$] In most cases, the data from Google mobility for residential activities are believed to prevent the spread of COVID-19. In some cases, they may result from events (such as home parties) that could promote COVID-19 cases.  Here we use Michaelis–Menten function to limit the effect of home activities in the spread of COVID-19.
           \item[$\ast$]  $c>0$ represents the maximum reduction rate of COVID-19 cases due to government measures and personal precautions to maintain social distancing including wearing masks, staying six feet apart. $k$ is the number of COVID-19 cases at which the reduction rate is $\frac{1}{2} c $. By adjusting the parameters, this model allows us to examine the effect of social distancing to prevent COVID-19 cases.
           \item[$\ast$] $h(x,t-14)$ describes home activities which are believed to prevent COVID-19 epidemics.  $h(x,t-14)$ is accumulated data from Google mobility for home activities.  $t-14$ is because it takes two weeks for the effect of the activity to take place.
          \end{itemize}

        \item $a(x,t-14), h(x,t-14)$ are determined by Google Mobility data with the interpolation function of Matlab.
        \item $d,k,c$ and parameters of $r(t), l(x),$ are determined by the known historical data of the COVID-19 cases with the best fitting function of Matlab.

       \item Neumann boundary condition $\frac{\partial C}{\partial x}(1,t)=\frac{\partial C}{\partial x}(3,t)=0, \quad t>1$ is applied in~ (\cite{natr2002murray}). For simplicity, we count the cases imported from neighbor states as local Arizona cases and assume that no COVID-19 spreads across the boundaries at $ x=1, 3$.

       \item Initial function $C(x,1)=\psi(x)$ describes the initial states of COVID-19 in every Arizona region, which can be constructed from the historical data of COVID-19 cases by cubic spline interpolation.
\end{itemize}

\section{Model Prediction}
\label{prediction}

\begin{figure}[!h]
\centering
\subfigure[]{\includegraphics[width=3.2in, trim=3cm 8cm 3cm 8cm, clip=true]{./Region1.pdf}}\
\subfigure[]{\includegraphics[width=3.2in, trim=3cm 8cm 3cm 8cm, clip=true]{./Region2.pdf}}\
\subfigure[]{\includegraphics[width=3.2in, trim=3cm 8cm 3cm 8cm, clip=true]{./Region3.pdf}}\
\caption{Predictions of COVID-19 cases in Arizona from April 24, 2020, to June 19, 2020, for three regions. The blue lines represent the predicted COVID-19 cases and red lines represent the actual COVID-19 cases.}\label{Predictions}
\end{figure}

\begin{figure}[!h]
\centering
\subfigure[]{\includegraphics[width=3.1in, trim=3cm 8cm 1cm 6cm, clip=true]{./RegionSp1.pdf}}\
\subfigure[]{\includegraphics[width=3.1in, trim=3cm 8cm 1cm 6cm, clip=true]{./RegionSp2.pdf}}\
\subfigure[]{\includegraphics[width=3.1in, trim=3cm 8cm 1cm 6cm, clip=true]{./RegionSp3.pdf}}\
\caption{The relative accuracy in each region: NAZ (Region 1), CAZ (Region 2), and SAZ (Region 3) from April 24, 2020 to June 19, 2020.}\label{RA}
\end{figure}

The basic mathematical properties of the proposed PDE model in Eq.~(\ref{new_model}) such as existence and uniqueness can be established from the standard theorems for parabolic PDEs in~\cite{friedman2008partial}. Below, we evaluate the robustness of our PDE-based predictive model and validate if the model has acceptable prediction performance with the COVID-19 case report in Arizona. The procedure of predictive modeling for the COVID-19 case in Arizona is summarized as follows:
\begin{itemize}

\item \textbf{Prediction Procedure:} In the prediction process of the research time period, the model parameters are to be determined by the historical aggregated data of COVID-19 cases in the three regions under the same structured PDE. In the current experiment, we predict the COVID-19 cases three days ahead. In order to forecast the COVID-19 of a given day, we first train the parameters of the PDE model and then solve the PDE model for prediction. For example, days 1--7, 2--8, \ldots, are used as the training data, and we predict the number of COVID-19 cases for the following days 10, 11, \ldots, respectively.

The process of performing the prediction can be divided into two major processes. We first use an optimization method to fit parameters in the PDE model with historical data of COVID-19 cases. In essence, this is a multi-parameter inverse problem of parabolic equations. We integrate the local and global methods to search for the best fitting parameters.  In the present work for COVID-19 modeling problems, we take a hybrid approach: first, a tensor train global optimization (\cite{oseledets2011tensor}) is used to explore the parameter space thus to locate the starting points and then Nelder-Mead simplex local optimization method (\cite{lagarias1998convergence}) is used to search the local optimization of COVID-19 modeling problem. The Nelder-Mead simplex method is implemented in the fminsearch function in Matlab. Once the model parameters are determined, we use the fourth-order Runge-Kutta method to solve the PDE for one-step forward prediction numerically.

 \item \textbf{Accuracy Measurement:}
We need to quantify the accuracy of model prediction of the COVID-19 cases by comparing the predicted COVID-19 cases with the observed aggregated COVID-19 cases for three regions, which are the ground truth. The mean absolute percentage error $$1-\left|\frac{x_{real}-x_{predict}}{x_{real}}\right|$$ is applied to measure the prediction accuracy, $x_{real}$ is the observed COVID-19 cases at every data collection time point and $x_{predict}$ is the predicted cases.

\end{itemize}

Fig.~\ref{Predictions} illustrates predictions of COVID-19 cases in Arizona from April 24, 2020, to June 19, 2020, for three regions. The period covers the period of the end of stay-home order when the COVID-19 case started to increase quickly. The blue lines represent the predicted COVID-19 cases and red lines represent the actual COVID-19 cases.  Note that we normalize the data to be between $[0,1]$. Parameters in each prediction step are different, here we just provide the parameters in the last prediction for June 19.  The parameters for prediction on June 19 are: $d=0.013556543719$, $k=4.655312010389$, $b_1=0.503402937243188$, $b_2=0.236417906835504$, $c=0.821575737918722$ and $l_1=0.0826411111659247$, $l_2=0.347970638795282$ and $l_3=0.28107732159028$. Fig.~\ref{RA} illustrates the prediction accuracy in NAZ, CAZ, and  SAZ every day from April 24, 2020, to June 19, 2020. The Arizona stay-home order started on March 31, 2020, and ended on May 15, 2020. The average relative accuracy of three regions with three days prediction are 94\%, 97\%, and 95\% respectively. Clearly, the PDE-based prediction model has a good short-term prediction ability during active COVID-19 epidemics.

\section{Effects of Social Distancing and Face Covering}

\begin{figure}[!h]
\centering
\subfigure[]{\includegraphics[width=3.2in, trim=3cm 8cm 3cm 8cm, clip=true]{./RedRegion1.pdf}}\
\subfigure[]{\includegraphics[width=3.2in, trim=3cm 8cm 3cm 8cm, clip=true]{./RedRegion2.pdf}}\
\subfigure[]{\includegraphics[width=3.2in, trim=3cm 8cm 3cm 8cm, clip=true]{./RedRegion3.pdf}}\
\caption{Effects of preventive measures on the growth of COVID-19 cases in Arizona from April 24, 2020, to June 19, 2020. The red lines represent the actual COVID-19 cases.  (a), (b), and (c) capture the predicted COVID-19 cases (blue lines) under additional social distancing and face covering to reduce COVID-19 for Region 1, 2 and 3. }\label{PolicyWhole}
\end{figure}

\begin{figure}[!h]
\centering
\includegraphics[width=0.5\textwidth]{Mask_plot.png}
\caption{Arizona counties with face covering order}
\label{mask}
\end{figure}

To demonstrate the effectiveness of intervention measures such as face masks and social distancing on decreasing COVID-19 cases, we predict the COVID-19 cases with  an increased value of $c$ which measures the effect of the stay-at-home-order and social distancing.  We do not have specific data from Google Community Mobility to incorporate our model for the effectiveness of social distancing and face mask. In general, it is believed that face masks and proper social distancing could significantly reduce COVID-19 cases. Specifically, based on each group of model parameters obtained in the above prediction process, we retain them unchanged in Eq.~(\ref{new_model}) except doubling the value of $c$ (the rate of effectiveness of social distancing).  The doubled value of $c$ suggests that, to reduce COVID-19 cases, people in Arizona take the extra effort of practicing social distancing and wearing face covering even after the lift of the stay-at-home-order.

Fig.~\ref{PolicyWhole} illustrates the COVID-19 cases with possible interventions, where the blue lines represent the predicted COVID-19 cases from April 24, 2020, to June 19, 2020, and the red lines represent the actual COVID-19 cases.   Fig.~\ref{Predictions}  can serve as baseline references. Fig.~\ref{PolicyWhole}  (a), (b), (c) captures the predicted COVID-19 cases (blue lines) under double measurements to reduce COVID-19 for Region 1, 2 and 3 respectively.   Fig.~\ref{PolicyWhole} shows that if stricter precaution measures are taken, the COVID-19 cases will be controlled to some extent. As a result, our analysis confirms the positive short-term effects of precaution measures on the spread of COVID-19.

Unfortunately, there is no mandatory requirement of face-covering at the state level.  About one month after the stay-at-home-order expired on May 15, local governments established their own health and safety measures such as face mask mandates in an effort to slow the spread of COVID-19. As a result, the state keeps breaking a record of daily cases and hospitalizations suspected to be related to COVID-19.  CAZ (Maricopa county) continues to have the highest number of cases in Arizona. On June 22, multiple officials from cities across the state began implementing requirements for face coverings in public. There have been five counties including Maricopa County that have a mask requirement in place as shown in Fig.~\ref{mask} as of June 29, 2020.  Many cities and towns across metro Phoenix already had enacted their own face-covering-mandates in public earlier but Maricopa County (CAZ) took the lead in sending a clear message to every resident in the county to wear a mask. We hope the COVID-19 cases will begin to go downward in the near future.

\section{Discussion}
\label{discussions}
In this paper, the proposed PDE model is built on three regions in Arizona for the prediction of COVID-19 cases. The PDE model is capable of capturing the effects of various geographical and social factors. To ensure high accuracy predictions of COVID-19 cases, daily cases of COVID-19, and the effects of social distancing are included in the prediction model.  In addition, Google Community Mobility data is incorporated to reflect the effects of human activities. As a result, the PDE model can be useful in making localized policies for reducing the spread of COVID-19.

We regroup 15 counties in Arizona into three regions: Northern Arizona, Central Arizona, and Southern Arizona. One may be able to use ordinary differential equations on multiple patches for the prediction of COVID-19 cases in Arizona. Here we choose partial differential equations instead of ordinary differential equations mainly because of the geographical and social characteristics in Arizona in terms of the spread of COVID-19. The three regions are geographically located in the north-south direction of Arizona and can be naturally projected onto the x-axis of the Cartesian coordinates. In addition, Central Arizona, home to the metro Phoenix area, has about 60\% of the population in Arizona.  As a result, the spread of COVID-19 between Central Arizona and other two regions are more likely predominant in terms of transboundary spread. These geographical characteristics, along with the disparities of social distancing and face mask requirements among the Arizona counties, fit squarely into the framework of reaction-diffusion equations for the spread of COVID-19.

While capturing the trend of the spread of COVID-19, we were able to avoid large scale computation by regrouping the counties, which allowed us to use a computational-heavy PDE model.  Thus, the diffusion process of COVID-19 between county-regions as well as the increase of COVID-19 cases within a region are incorporated into the PDE model. The present paper aims to develop, for the first time in the literature, a specific PDE model between county-regions for forecasting daily COVID-19 cases and further quantifying the influences of the mobility of the community on the changing the number of COVID-19 cases in the state of Arizona. In addition,  In the present study, we study the effectiveness of human precautions such as face masks and social distancing on decreasing COVID-19 cases and hope to provide analytical results to help slow the spread of COVID-19 in Arizona.

In this paper, we choose three regions of counties in Arizona for the spatial model because of the availability of both COVID-19 data and Google Community Mobility data. To gain more meaningful predictions with finer granularity, one may use zip codes or cities for the construction of regions or clusters. We could choose clustering algorithms to determine appropriate regions to ensure the connected zip codes stay as much as possible. However, the Google Community Mobility data is not available at the zip code level. Therefore, our model for prediction for regional COVID-19 is based on three regions of counties which is the most cost-effective method for prediction of COVID-19 in Arizona.

Arizona has emerged as a COVID-19 hot spot since the stay-home orders were lifted in mid-May. It is important to study how to slow the virus’ spread. This study first presents a PDE model to describe  COVID-19  cases from community transmission, transboundary transmission between regions, and human activities. We apply this model for the prediction of COVID-19 cases to see the effect of social distancing on COVID-19 cases. The numerical results show that our new method of predicting COVID-19  cases is of acceptable high accuracy and can predict the effects of social distancing. To conclude, the prediction accuracy and the policy simulation suggest that a combining PDE model with real data can provide substantial information on COVID-19 spatio-temporal dynamics.

\section*{Acknowledgments}
The first author was supported by NSF grant \#1737861.
\section*{Conflict of interest}
 The authors declare that they have no conflict of interest.