#  System inference via field inversion for the spatio-temporal progression of
infectious diseases: Studies of COVID-19 in Michigan and Mexico

Z. Wang, M. Carrasco-Teja, X. Zhang, G.H. Teichert & K. Garikipati
Mechanical Engineering, Mathematics and the Michigan Institute for
Computational
Discovery & Engineering, University of Michigan Corresponding author,
<EMAIL_ADDRESS>
###### Abstract

We present an approach to studying and predicting the spatio-temporal
progression of infectious diseases. We treat the problem by adopting a partial
differential equation (PDE) version of the Susceptible, Infected, Recovered,
Deceased (SIRD) compartmental model of epidemiology, which is achieved by
replacing compartmental populations by their densities. Building on our recent
work (_Computational Mechanics_ , 66, 1177, 2020), we replace our earlier use
of global polynomial basis functions with those having local support, as
epitomized in the finite element method, for the spatial representation of the
SIRD parameters. The time dependence is treated by inferring constant
parameters over time intervals that coincide with the time step in semi-
discrete numerical implementations. In combination, this amounts to a scheme
of field inversion of the SIRD parameters over each time step. Applied to data
over ten months of 2020 for the pandemic in the US state of Michigan and to
all of Mexico, our system inference via field inversion infers spatio-
temporally varying PDE SIRD parameters that replicate the progression of the
pandemic with high accuracy. It also produces accurate predictions, when
compared against data, for a three week period into 2021. Of note is the
insight that is suggested on the spatio-temporal variation of infection,
recovery and death rates, as well as patterns of the population’s mobility
revealed by diffusivities of the compartments.

##  1 Introduction

Classical epidemiological models, such as the Susceptible-Infected-Recovered
(SIRD) model [1], are ordinary differential equations (ODEs) defined by
specifying the compartmental sub-population numbers over some geographical
region. Spatial effects have typically been introduced by resolving smaller
regions and treating them individually[2, 3, 4, 5, 6]. During the long-lasting
and widespread epidemics, such as the COVID-19 Pandemic, the effects on the
infection rate of imposing–and then lifting–mobility restrictions and social
distancing mandates revolve on the question of the time and spatially varying
mobility of the population. At the finest resolution, this must be approached
via agent-based models[7], using individuals’ mobility data. However, this
data is not available for the entire population, and contact tracing campaigns
face challenges of recruiting workers, access, technology, as well as socio-
political resistance. Against these difficulties, an intriguing question to
explore is whether simple reaction-diffusion models can detect the evidence of
mobility in these data. Such an approach must start with a partial
differential equation (PDE) version of the epidemiological models, which is
easily defined by converting compartmental sub-populations to densities over
sub regions by normalizing with the corresponding areas. To address the
mobility of the population, diffusion terms are introduced to the SIRD model,
which is transformed to a set of reaction-diffusion PDEs in two spatial
dimensions [8, 9, 10, 11].

The widespread availability of data in the public domain [12, 13, 14, 15, 16,
17, 18, 19] has spurred a widespread interest among computational and data
scientists, who have sought to test and refine their methods against these
repositories. This has opened up the possibility that advances in
computational and data science may contribute to the existing and rapidly
expanding body of work in epidemiology, in inferring the dynamics of COVID-19
and making projections. We have similarly sought to build off our recent work
in data-driven and machine learning approaches [20, 21, 22, 23, 24, 25, 26,
27, 28, 29] and presented a class of system identification techniques for
inference of ODE and PDE forms of the SIRD model, as well as Bayesian neural
networks for representation and uncertainty quantification-guided prediction
[11]. That work focused on the US state of Michigan. In this communication, we
revise our approach for inference of the PDE SIRD model with _temporal_ and
_spatial_ evolving parameters and diffusivities. Importantly, instead of
global polynomial representations of PDE SIRD parameters over the spatial and
temporal domains, we adopt field inversion over time intervals that coincide
with the time steps of our underlying numerical implementation. This affords
much greater accuracy over the global polynomial ansatz. Adjoint-based
gradient optimization for field inversion of parameters at each time step
replaces the use of stepwise regression-based system identification in our
previous work. We find that the improved accuracy with respect to the data
over the time interval of inference, as well as of the predictions, is worth
the increased expense. We have brought abundant high-quality, public domain,
data [12, 13, 14, 15, 16, 17, 18, 19] on the evolution of COVID-19 in both,
the State of Michigan, with a population of 9.98 million, distributed in 83
counties, over 250,493 km2, and the country of Mexico, with a population of
126 million, distributed in 32 geographical entities (31 states, plus Mexico
City), on 1,972,550 km2. The temporal resolution by days and spatial
resolution by counties/states have allowed us to study the mobility in these
data using our methods of system inference.

In Section 2 we review our previous work of system inference for the spatio-
temporal SIRD model first, and then extend it by incorporating temporal and
spatial parameters and diffusivities using a finite element representation.
The PDE SIRD model-constrained inference are presented in Section 3. Section 4
is on data preparation. The results for inference of classical SIRD parameters
as well as the diffusivities, and forward prediction are presented in Section
5. Our conclusions appear in Section 6.

##  2 Compartmental differential equations models of infectious disease
dynamics

We begin with the conventional SIRD compartmental epidemiology model. The
population, taken to remain constant at $N$, is divided into four disjoint
compartments with time-dependent sub-populations: $S(t)$ for susceptible,
$I(t)$ for infected, $R(t)$ for recovered and $D(t)$ for deceased individuals.
The governing system of ordinary differential equations (ODEs) is:

$\displaystyle\frac{\text{d}S}{\text{d}t}$
$\displaystyle=-\frac{\beta}{N}SI+\gamma R$ (1)
$\displaystyle\frac{\text{d}I}{\text{d}t}$
$\displaystyle=\frac{\beta}{N}SI-\mu I-\alpha I$ (2)
$\displaystyle\frac{\text{d}R}{\text{d}t}$ $\displaystyle=\mu I-\gamma R$ (3)
$\displaystyle\frac{\text{d}D}{\text{d}t}$ $\displaystyle=\alpha I$ (4)
$\displaystyle N$ $\displaystyle=S(t)+I(t)+R(t)+D(t).$ (5)

This is the canonical form of the model where the sub-populations are assumed
to be well-mixed so that spatial variations can be ignored over the domain of
interest. Here $\beta$ is the infection rate, $\mu$ is the recovery rate,
$\gamma$ is the rate of immunity loss, and $\alpha$ is the death rate.

We have extended the SIRD model to a system of partial differential equations
(PDEs) in two spatial dimensions using the same compartments [11]. However,
the population variables are now replaced with spatio-temporally varying
densities,
$\widehat{S}(\boldsymbol{x},t),\widehat{I}(\boldsymbol{x},t),\widehat{R}(\boldsymbol{x},t),\widehat{D}(\boldsymbol{x},t)$
defined as numbers per unit area.

$\displaystyle\frac{\partial\widehat{S}}{\partial t}$
$\displaystyle=\mathcal{D}_{\text{S}}\nabla^{2}\widehat{S}-\frac{\beta}{\widehat{N}}\widehat{S}\widehat{I}+\gamma\widehat{R}$
(6) $\displaystyle\frac{\partial\widehat{I}}{\partial t}$
$\displaystyle=\mathcal{D}_{\text{I}}\nabla^{2}\widehat{I}+\frac{\beta}{\widehat{N}}\widehat{S}\widehat{I}-\mu\widehat{I}-\alpha\widehat{I}$
(7) $\displaystyle\frac{\partial\widehat{R}}{\partial t}$
$\displaystyle=\mathcal{D}_{\text{R}}\nabla^{2}\widehat{R}+\mu\widehat{I}-\gamma\widehat{R}$
(8) $\displaystyle\frac{\partial\widehat{D}}{\partial t}$
$\displaystyle=\alpha\widehat{I}$ (9)

where $\mathcal{D}_{\text{S}},\mathcal{D}_{\text{I}},\mathcal{D}_{\text{R}}$
are diffusivities of the corresponding compartments, and represent the
mobility of the sub-population via random walks. We define
$\widehat{(\bullet)}=(\bullet)/\int_{\Omega}\text{d}A$ where $\Omega$ is the
domain of study: either the lower peninsula of the State of Michigan, or the
territory of the country of Mexico. Furthermore the population constraint
holds:
$\int_{\Omega}\widehat{N}\text{d}A=\int_{\Omega}\widehat{S}(t)\text{d}A+\int_{\Omega}\widehat{I}(t)\text{d}A+\int_{\Omega}\widehat{R}(t)\text{d}A+\int_{\Omega}\widehat{D}(t)\text{d}A$.
In what follows of this communication, we only consider the PDEs SIRD model,
and, for the sake of readability, we dispense with the hats on the
compartments.

We adopt the weak form, and specifically, the finite element framework for the
above system of PDEs. For a generic, finite-dimensional field $u^{h}$, the
problem is stated as follows: find
$u^{h}\in\mathscr{S}^{h}\subset\mathscr{S}$, where
$\mathscr{S}^{h}=\\{u^{h}\in\mathscr{H}^{1}(\Omega)~{}|~{}u^{h}=~{}\bar{u}\;\mathrm{on}\;\Gamma^{u}\\}$,
such that $\forall~{}w^{h}\in\mathscr{V}^{h}\subset\mathscr{V}$, where
$\mathscr{V}^{h}=\\{w^{h}\in\mathscr{H}^{1}(\Omega)~{}|~{}w^{h}=~{}0\;\mathrm{on}\;\Gamma^{u}\\}$,
the finite-dimensional (Galerkin) weak form of the problem is satisfied. The
variations $w^{h}$ and trial solutions $u^{h}$ are defined component-wise
using a finite number of basis functions,

$w^{h}=\sum_{a=1}^{n_{\mathrm{b}}}c^{a}N^{a},\quad\qquad
u^{h}=\sum_{a=1}^{n_{\mathrm{b}}}d^{a}N^{a},$ (10)

where $n_{\mathrm{b}}$ is the dimensionality of the function spaces
$\mathscr{S}^{h}$ and $\mathscr{V}^{h}$, and $N^{a}$ represents the basis
functions. To obtain the Galerkin weak forms, we multiply each equation in
strong form in (6-9) by a weighting function
$w_{\text{S}}^{h},w_{\text{I}}^{h},w_{\text{R}}^{h},w_{\text{D}}^{h}$,
respectively, integrate by parts, apply boundary conditions appropriately, and
use the Backward Euler method for time-discretization with $(\bullet)_{n}$
denoting a discretized quantity at time $t_{n}$ and $\Delta t$ being the time
step. See Ref. [11] for details. This leads to:

$\displaystyle\int_{\Omega}w^{h}_{\text{S}}\frac{S^{h}_{n}-S^{h}_{n-1}}{\Delta
t}\text{d}s$ $\displaystyle=-\int_{\Omega}\mathcal{D}_{\text{S}}\nabla
w^{h}_{\text{S}}\cdot\nabla
S^{h}_{n}\text{d}s-\int_{\Omega}w^{h}_{\text{S}}\left(\frac{\beta}{N}S^{h}_{n}I^{h}_{n}+\gamma
R^{h}_{n}\right)\text{d}s$ (11)
$\displaystyle\int_{\Omega}w^{h}_{\text{I}}\frac{I^{h}_{n}-I^{h}_{n-1}}{\Delta
t}\text{d}s$ $\displaystyle=-\int_{\Omega}\mathcal{D}_{\text{I}}\nabla
w^{h}_{\text{I}}\cdot\nabla
I^{h}_{n}\text{d}s+\int_{\Omega}w^{h}_{\text{I}}\left(\frac{\beta}{N}S^{h}_{n}I^{h}_{n}-\mu
I^{h}_{n}-\alpha I^{h}_{n}\right)\text{d}s$ (12)
$\displaystyle\int_{\Omega}w^{h}_{\text{R}}\frac{R^{h}_{n}-R^{h}_{n-1}}{\Delta
t}\text{d}s$ $\displaystyle=-\int_{\Omega}\mathcal{D}_{\text{R}}\nabla
w^{h}_{\text{R}}\cdot\nabla
R^{h}_{n}\text{d}s+\int_{\Omega}w^{h}_{\text{R}}\left(\mu I^{h}_{n}-\gamma
R^{h}_{n}\right)\text{d}s$ (13)
$\displaystyle\int_{\Omega}w^{h}_{\text{D}}\frac{D_{n}-D_{n-1}}{\Delta
t}\text{d}s$ $\displaystyle=\int_{\Omega}w^{h}_{\text{D}}\alpha
I^{h}_{n}\text{d}s$ (14)

where the boundary terms vanish because we assume that the sub-populations do
not leave the region–zero flux boundary conditions.

In our previous work [11], we have characterized the coefficients to vary via
a global-in-time polynomial basis. While the inferred model reproduced the
trends, there was a notable error over time of the statewide sub-populations
estimates
$S(t,\mbox{\boldmath$x$}),I(t,\mbox{\boldmath$x$}),R(t,\mbox{\boldmath$x$}),D(t,\mbox{\boldmath$x$})$
obtained by forward simulation with inferred quantities (See Figure 14 and 15
in [11]).

(a) Michigan State

(b) Mexico

Figure 1: Reconstructed maps of Michigan State and Mexico with finite element
mesh.

Additionally, the highly complex geometry of the State of Michigan, and of
Mexico (See maps in Figure 1), and potentially highly nonuniform distributions
of the coefficients in space makes it challenging to characterize the
coefficients with simple basis functions. Global polynomials in space could
not sufficiently resolve the emergence and disappearance of “hot spots” and
“cold spots” [11]. In this communication, we allow the coefficients
$\beta,\gamma,\mu,\alpha,\mathcal{D}_{\text{S}},\dots,\mathcal{D}_{\text{R}}$
of the PDE SIRD model to vary over space via finite-dimensional, locally
supported representations as we do for the primary variables
$S(t,\mbox{\boldmath$x$}),I(t,\mbox{\boldmath$x$}),R(t,\mbox{\boldmath$x$}),D(t,\mbox{\boldmath$x$})$.
Further more, we allow the coefficients to vary daily, leading to:

$\displaystyle\beta^{h}_{n}=\sum_{a=1}^{n_{p}}\beta^{a}_{n}N^{a},\qquad\gamma^{h}_{n}=\sum_{a=1}^{n_{p}}\gamma^{a}_{n}N^{a},\qquad\mu^{h}_{n}=\sum_{a=1}^{n_{p}}\mu^{a}_{n}N^{a}\qquad\alpha^{h}_{n}=\sum_{a=1}^{n_{p}}\alpha^{a}_{n}N^{a}$
(15)
$\displaystyle\mathcal{D}_{\text{S}_{n}}^{h}=\sum_{a=1}^{n_{p}}\mathcal{D}_{\text{S}_{n}}^{a}N^{a},\qquad\mathcal{D}_{\text{I}_{n}}^{h}=\sum_{a=1}^{n_{p}}\mathcal{D}_{\text{I}_{n}}^{a}N^{a},\qquad\mathcal{D}_{\text{R}_{n}}^{h}=\sum_{a=1}^{n_{p}}\mathcal{D}_{\text{R}_{n}}^{a}N^{a}$
(16)

where, as for the primary variables, the subscripts $(\bullet)n$ denote the
coefficients on day $n$. With this, the PDE SIRD equations become:

$\displaystyle\int_{\Omega}w^{h}_{\text{S}}\frac{S^{h}_{n}-S^{h}_{n-1}}{\Delta
t}\text{d}s$ $\displaystyle=-\int_{\Omega}\mathcal{D}^{h}_{\text{S}_{n}}\nabla
w^{h}_{\text{S}}\cdot\nabla
S^{h}_{n}\text{d}s-\int_{\Omega}w^{h}_{\text{S}}\left(\frac{\beta^{h}_{n}}{N}S^{h}_{n}I^{h}_{n}+\gamma^{h}_{n}R^{h}_{n}\right)\text{d}s$
(17)
$\displaystyle\int_{\Omega}w^{h}_{\text{I}}\frac{I^{h}_{n}-I^{h}_{n-1}}{\Delta
t}\text{d}s$ $\displaystyle=-\int_{\Omega}\mathcal{D}^{h}_{\text{I}_{n}}\nabla
w^{h}_{\text{I}}\cdot\nabla
I^{h}_{n}\text{d}s+\int_{\Omega}w^{h}_{\text{I}}\left(\frac{\beta^{h}_{n}}{N}S^{h}_{n}I^{h}_{n}-\mu^{h}_{n}I^{h}_{n}-\alpha^{h}_{n}I^{h}_{n}\right)\text{d}s$
(18)
$\displaystyle\int_{\Omega}w^{h}_{\text{R}}\frac{R^{h}_{n}-R^{h}_{n-1}}{\Delta
t}\text{d}s$ $\displaystyle=-\int_{\Omega}\mathcal{D}^{h}_{\text{R}_{n}}\nabla
w^{h}_{\text{R}}\cdot\nabla
R^{h}_{n}\text{d}s+\int_{\Omega}w^{h}_{\text{R}}\left(\mu^{h}_{n}I^{h}_{n}-\gamma^{h}_{n}R^{h}_{n}\right)\text{d}s$
(19) $\displaystyle\int_{\Omega}w^{h}_{\text{D}}\frac{D_{n}-D_{n-1}}{\Delta
t}\text{d}s$
$\displaystyle=\int_{\Omega}w^{h}_{\text{D}}\alpha^{h}_{n}I^{h}_{n}\text{d}s$
(20)

where the PDE SIRD parameters are interpolated from nodal variables as defined
in Equations (15) and (16).

##  3 System inference by field inversion using adjoint-based gradient
optimization

The system inference problem is to invert for the quantities
$\beta^{a}_{n},\gamma^{a}_{n},\mu^{a}_{n},\alpha^{a}_{n},\mathcal{D}^{a}_{\text{S}_{n}},\dots,\mathcal{D}^{a}_{\text{R}_{n}}$.
Since these quantities are interpolated via Equations (15) and (16) to be
expressed as the corresponding fields
$\beta^{h}_{n},\gamma^{h}_{n},\mu^{h}_{n},\alpha^{h}_{n},\mathcal{D}^{h}_{\text{S}_{n}},\dots,\mathcal{D}^{h}_{\text{R}_{n}}$
in (17-18), the system inference problems is one of field inversion. It is
stated in Equations (21-22) as:

$\displaystyle\text{Given (\ref{eq:SIRDparamsinterp}-\ref{eq:SIRDdiffinterp}),
at
each}\;t_{n}:\quad\left(\beta^{a}_{n},\dots,\mathcal{D}^{a}_{\text{R}_{n}}\right)_{a=1}^{n_{p}}$
$\displaystyle=\text{arg
}\underset{(\beta^{a}_{n},\dots,\mathcal{D}^{a}_{\text{R}_{n}}){a=1}^{n_{p}}}{\min}\text{
}\ell_{i},\quad\text{such that (\ref{eq:S_2D_weak_inv}-\ref{eq:D_2D_weak_inv})
hold}$ (21)

and $\ell_{i}$ is the loss function defined:

$\displaystyle\ell_{i}=$
$\displaystyle\int_{\Omega}W_{\text{S}}\left(S^{h}_{n}-S_{n}^{\text{d}}\right)^{2}+W_{\text{I}}\left(I^{h}_{n}-I_{i}^{\text{d}}\right)^{2}+W_{\text{R}}\left(R^{h}_{n}-R_{n}^{\text{d}}\right)^{2}+W_{\text{D}}\left(D^{h}_{n}-D_{n}^{\text{d}}\right)^{2}\text{d}v$
(22)

where $(\bullet)^{\text{d}}$ denotes data for the corresponding quantity. Due
to the large differences in the magnitudes of different sub-populations, we
choose the weights $W_{\text{S}},\cdots,W_{\text{D}}$ to be:

$\displaystyle
W_{\text{S}}=\frac{I_{n}^{\text{d}}}{\text{mean}\left(S_{n}^{\text{d}}\right)},\quad
W_{\text{I}}=\frac{I_{n}^{\text{d}}}{\text{mean}\left(I_{n}^{\text{d}}\right)},\quad
W_{\text{R}}=\frac{I_{n}^{\text{d}}}{\text{mean}\left(R_{n}^{\text{d}}\right)},W_{\text{D}}=\frac{I_{n}^{\text{d}}}{\text{mean}\left(D_{n}^{\text{d}}\right)}.$
(23)

The weights normalize the sub-populations and prioritize regions with higher
infected populations. These regions are of greater interest for studying the
progression of the disease as they tend to have a higher population density
and, therefore, infected populations.

This PDE-constrained optimization problem is solved iteratively, and requires
the gradient of the PDE constraint, Equations (17-20), with respect to
parameters. We adopt classical adjoint-based gradient optimization. This
approach involves a single linear solution of the adjoint equation of the
original PDE constraint at each iteration, followed by solution of the fields
to be inverted:
$(\beta^{h}_{n},\gamma^{h}_{n},\mu^{h}_{n},\alpha^{h}_{n},\mathcal{D}^{h}_{\text{S}_{n}},\dots,\mathcal{D}^{h}_{\text{R}_{n}})$
and the updated forward solution $S^{h}_{n},I^{h}_{n},R^{h}_{n},D^{h}_{n}$. In
this work we use the L-BFGS-B optimization algorithm from SciPy[30] and the
dolfin-adjoint software library [31] to compute the gradient.

##  4 Data preparation on maps of Michigan and Mexico

First, we constructed two-dimensional meshes for Michigan and Mexico that
fully resolve the counties/states as shown in Figure 1. The data are available
as cumulative sub-population numbers
$I^{\text{d}}_{n},R^{\text{d}}_{n},D^{\text{d}}_{n}$ at the county/state
level. We used a uniform density of each sub-population to compute
$I^{\text{d}}_{n},R^{\text{d}}_{n},D^{\text{d}}_{n}$ within the county/state,
and applied Gaussian filtering to smooth the discontinuities at the
county/state boundaries. Note that the discrete Gaussian filter can not be
applied in a straightforward manner to unstructured meshes. Starting with a
field $u$ that represents any of the four sub-population densities, and
$G(\mbox{\boldmath$x$}_{0},\mbox{\boldmath$x$})=\frac{1}{2\pi\sigma^{2}}e^{-\frac{||\mbox{\boldmath$x$}||^{2}}{2\sigma^{2}}}$
as the two dimensional Gaussian distribution function centered at
$\boldsymbol{0}$ with standard deviation $\sigma$, which is related to the
kernel size in the discrete Gaussian filter, we scale the filtered solution
denoted by $u(\boldsymbol{x}_{0})$ at each finite element node:

$\displaystyle
u(\mbox{\boldmath$x$}_{0})=\frac{1}{{\int_{\Omega}G(\mbox{\boldmath$x$}_{0},\mbox{\boldmath$x$})\text{d}v}}{\int_{\Omega}G(\mbox{\boldmath$x$}_{0},\mbox{\boldmath$x$})u_{\text{raw}}(\mbox{\boldmath$x$})\text{d}v}$
(24)

The spatio-temporal evolution of these fields was used in the system inference
problem as described in Section 3.

##  5 Results

Figure 2 shows the sub-populations
$S(\mbox{\boldmath$x$},t),I(\mbox{\boldmath$x$},t),R(\mbox{\boldmath$x$},t),D(\mbox{\boldmath$x$},t)$
in both Michigan and Mexico obtained by forward simulation with inferred
quantities compared with data on December 29, 2020 ($t=281$ days), where $t=0$
is March 23, 2020, the start of the lockdown in Michigan. The inferred model
for Michigan accurately replicates the initial burst of disease and the
following multiple waves around Detroit (please see the SI movie:
michigan_prediction.mp4). It also captures the second burst in the southwest
of Michigan around the city of Grand Rapids. The high burden of the disease in
these, the largest and second largest cities, respectively, in Michigan,
reflects well-known socio-economic challenges related to Detroit in
particular, and more generally reflected in other urban centers. Similarly,
Mexico City, with highest population density in Mexico ($6,200/\text{km}^{2}$
[18]), was the worst affected area in that country and dominated the evolution
of the disease (See SI movie: mexico_prediction.mp4).

(a) Michigan

(b) Mexico

Figure 2: Comparison of the data on distributions of the susceptible (S),
infected (I), recovered (R) and deceased (D) sub-populations against forward
PDE SIRD simulations with inferred quantities respectively. Data and
simulation results are shown for Day 281 starting from March 23, 2020. The
inset adjacent to the map of Mexico is a magnified view of the region
surrounding Mexico City. The time-resolved dynamics may be seen as movies in
Supplementary Information: michigan_prediction.mp4 and mexico_prediction.mp4.

The low error between the simulation and data leads to greater confidence in
the inferred parameters. Figure 3 shows the inferred infection rate, death
rate and the recovery rate in Michigan every 70 days starting from $t=0$ (the
time-resolved dynamics are shown in SI movie: michigan_parameter.mp4). The
evolution of these inferred parameters reveals that the population’s infection
rate, $\beta(t)$, declined from the initially higher values in highly infected
areas (such as Detroit), and spread to the western parts of Michigan. The
death rate was mostly stable after May 2020 ($t>69$), and remained low in the
more highly infected areas. This can be attributed to the ramp up of the
public health campaign, hospitalizations and emergency response of the medical
system, and prioritization to the more highly infected areas. The recovery
rate around Detroit city evolved in multiple stages: increasing$\to$
decreasing $\to$ increasing, which was consistent with the multiple waves
reflected in the data on the recovered population in this region (SI movie:
michigan_prediction.mp4).

Figure 3: Inferred infection rate ($\beta$), death rate($\alpha$) and recovery
rate ($\mu$) over the lower peninsula of Michigan starting from the lockdown
on March 23, 2020. The time-resolved dynamics may be seen as movies in
Supplementary Information: michigan_parameter.mp4.

At the finest resolution, the mobility of the population during disease
evolution may be approached via agent-based models refined to resolve
individuals. However, given the difficulties encountered in effective contact
tracing, and its acceptance by the population [9, 32, 33], an intriguing
question to explore is whether simple reaction-diffusion models can detect the
evidence of mobility in these data. Figure 4 shows the inferred diffusivities
of the susceptible, infected, and recovered sub-populations. Note that for
field inversion, the population density data for each compartment was taken to
be uniform within each county/state, since no finer grained information was
available, and then subject to Gaussian smoothing before inference. Thus the
density gradients, which drive the inference of diffusivities, arise at the
counties/states scale more than they do at the intra-county/intra-state.
Accordingly, the inferred diffusivities are meaningful on this scale. The
lower Peninsula of Michigan is about 446 km long from north to south and 314
km wide from east to west–scales that can help place the diffusivities in
Figure 4 in perspective. The mobility of the infected population was always
high around the highly infected areas. In Michigan, this infected population
gradually shifted to the southwestern part of the state from the initial burst
around Detroit. This finding is consistent with the second burst around Grands
during the evolution of the pandemic. The recovered population demonstrated a
similar pattern of mobility, and was more active in the southern part of
Michigan around the more highly infected regions. On the other hand,
susceptible population closely tracks the total population. Since the
population at large has low mobility, the susceptible population’s mobility is
low in high population density areas. See SI movie michigan_prediction.mp4 for
these dynamics.

Figure 4: Inferred diffusivities of the susceptible, infected, and recovered
sub-populations in Michigan starting from lockdown, March 23, 2020. The lower
Peninsula of Michigan is about 446 km long from north to south and 314 km wide
from east to west–scales that can help place the diffusivities in perspective.
Time-resolved dynamics maybe seen as movies in Supplementary Information:
mexico_parameter.mp4.

Figure 5 show the inferred infection rate, death rate and the recovery rate of
the inferred model for Mexico. We can clearly see the spreading of the disease
from Mexico City. Similar to the case of Michigan, infection rates, and to a
lesser extent, death rates, were relatively lower in Mexico City, which is the
most densely populated region of the country, than that in the surrounding
cities. The recovery rate was high in Mexico city due to the relatively
greater resources of the medical system there. The infection and death rates
tended to be stable for five months following March 23, 2020, and the recovery
rate gradually increased in more areas. Notably, far from the Mexico city,
Baja California also displayed a high inferred rate of infection. We suspect
this to be because it borders California, USA, and the international border
restrictions did not contain the spread of the virus between the two regions.
Unlike Mexico City, the death rate remained high, and the recovery rate did
not increase to levels comparable to the capital, perhaps because of the
looser restrictions in this popular tourist destination.

Figure 5: Inferred infection rate ($\beta$), death rate($\alpha$), recovery
rate ($\mu$) in Mexico starting from lockdown, March 23, 2020. Time-resolved
dynamics may be seen as movies in Supplementary Information:
michigan_parameter.mp4.

The diffusivities of the corresponding sub-populations of the inferred model
for Mexico are shown in Figure 6. Similar to Michigan, the mobility of
infected and recovered sub-populations are higher around the highly infected
Mexico City. Mexico is about 3000 km long from north to south and 1900 km wide
from east to west–scales that can help place the diffusivities in Figure 6 in
perspective. Unlike the case in Michigan where there were multiple bursts in
different cities, the mobilities of all sub-populations became stable after
about 5 months from March 23, 2020. This may reflect differences in the
proclivity toward domestic/local mobility of the populations of Michigan and
Mexico–two regions with strongly contrasting social, economic and cultural
characteristics.

Figure 6: Inferred diffusivities of the susceptible, infected, and recovered
sub-populations at day 281 at Mexico starting from lockdown, March 23, 2020.
Mexico is about 3000 km long from north to south and 1900 km wide from east to
west–scales that can help place the diffusivities in perspective. The time-
resolved dynamics may be seen as movies in Supplementary Information:
mexico_parameter.mp4.

Finally, taking the inferred parameters on the last day used for inference
(Day 281), we predicted the evolution of sub-populations for three weeks (Days
282 to 303) using the inferred model. Figures 7 and 8 show the predicted
spatio-temporal evolution of the infected population against the raw data for
both Michigan and Mexico. The inferred models captured closely the evolution
of the infected-sub-populations, indicating that the dynamics of the disease
tended to be steady in January 2021. The prediction of recovered and deceased
sub-populations are shown in Figures 9 and 10 under Supplementary Information.

Figure 7: Prediction of infected population for 3 weeks against the data for
Michigan. Figure 8: Prediction of infected population for 3 weeks against the
data for Mexico. The inset adjacent to the map of Mexico is a magnified view
of the region around Mexico City.

##  6 Conclusion

This communication builds upon our previous work [11] on system inference and
machine learning from data to study the progression of COVID-19 across the
state of Michigan. We extended the PDE SIRD model by allowing the infection
rate, death rate and the recovery rate, as well as the diffusivities of the
susceptible, infected, and recovered sub-populations to vary over space and
time. Using field inversion to infer the parameters as finite-dimensional
fields on time scales of a single day, we obtained models to predict the
evolution of disease with high accuracy. This provides us with the ability to
analyze the dynamics of the disease through the inferred parameters, and make
accurate predictions within a reasonable time frame. Particularly, we can
detect the evidence of time and spatially varying mobility of the population
through the simple diffusion-reaction models instead of the relying on the
agent-based models which require individual’s mobility data. The latter can
prove challenging, technically as well as politically, to obtain.

As discussed in Section 5, our inferred models capture the geographical spread
of infection, the number of deaths and the size of the recovered population
starting from one highly infected area to its surrounding cities and
eventually spreading to further areas. Particularly, the higher infection and
death rates in areas with low infection at later times suggests that more
attention is needed in such locations. This may be due to a lack of medical
services, or a lack of compliance with mitigation strategies. Our inferred
models also reveals higher mobility surrounding the highly infected areas
suggesting the importance of quarantine and social distancing.

Finite-dimensional representation allows the parameters to accurately capture
the spatial dependence, however the non-parametric representation makes the
projection of these parameters beyond the data range extremely challenging. A
prediction cannot be made with confidence if the dynamics of the disease
reflected by these parameters are not stable. Of course, extrapolation is
challenging in almost all data-driven methods. One possible alternate is to
develope surrogate models of these parameters via time dependent neural
networks under the constraints of the SIRD model to learn the spatial
variation in time, and thus to make reasonable prediction of the dynamics in
the evolution the disease, such as we have demonstrated previously [11].
Nevertheless, without including factors such as mobility restrictions or other
mandates, only short time predictions may be accurate.

## Acknowledgements

We acknowledge the support of Defense Advanced Research Projects Agency
(DARPA) under Agreement No. HR0011199002, “Artificial Intelligence guided
multi-scale multi-physics framework for discovering complex emergent materials
phenomena”

##  Appendix A Appendix: Additional results

Figure 9: Prediction of the recovered and deceased sub-populations for 3 weeks
against with the data at Michigan.

Figure 10: Prediction of the recovered and deceased sub-populations for 3
weeks against with the data at Mexico. The inset in the Mexico plot is a
magnified view of Mexico City area.

## References

  * Kermack and McKendrick [1927] W.O. Kermack and A.G. McKendrick.  A contribution to the mathematical theory of epidemics.  _Proceedings of the Royal Society of London, Series A_ , 115:700–721, 1927.
  * Eisenberg et al. [2015] M.C. Eisenberg, J.N.S. Eisenberg, J.P. D’Silva, E.V. Wells, S. Cherng, Yu-H. Kao, and R. Meza.  Forecasting and uncertainty in modeling the 2014-2015 ebola epidemic in west africa.  2015\.
  * Eisenberg et al. [2013] Marisa Eisenberg, Gregory Kujbida, Ashleigh R. Tuite, David N. Fisman, and Joseph H. Tien.  Examining rainfall and cholera dynamics in haiti using statistical anddynamic modeling approaches.  _Epidemics_ , 5:197–207, 2013.  doi: 10.1016/j.epidem.2013.09.004.
  * Wesolowski et al. [2012] Amy Wesolowski, Nathan Eagle, Andrew J. Tatem, David L. Smith, Abdisalan M. Noor, Robert W. Snow, and Caroline O. Buckee.  Quantifying the impact of human mobility on malaria.  _Science_ , 338:267–270, 2012.  doi: 10.1126/science.1223467.
  * Colizza et al. [2007] Vittoria Colizza, Alain Barrat, Marc Barthelemy, Alain-Jacques Valleron, and Alessandro Vespignani.  Modeling the worldwide spread of pandemic influenza: Baseline case and containment interventions.  _PLoS Medicine_ , 4:e13, 2007.  doi: 10.1371/journal.pmed.0040013.
  * Hethcote [2000] Herbert W Hethcote.  The mathematics of infectious diseases.  _SIAM review_ , 42(4):599–653, 2000.
  * Hunter et al. [2017] Elizabeth Hunter, Brian Mac Namee, and John Kelleher.  A taxonomy for agent-based models in human infectious disease epidemiology.  _Journal of Artificial Societies and Social Simulation_ , 20(3):2, 2017.  doi: 10.18564/jasss.3414.  URL http://jasss.soc.surrey.ac.uk/20/3/2.html.
  * Viguerie et al. [2020] Viguerie, Alex, Veneziani, Alessandro, Lorenzo, Guillermo, Baroli, Davide, Aretz-Nellesen, Nicole, Patton, Alessia, Yankeelov, Thomas E., Reali, Alessandro, Hughes, Thomas J. R., Auricchio, and Ferdinando.  Diffusion–reaction compartmental models formulated in a continuum mechanics framework: application to covid-19, mathematical analysis, and numerical study.  _Computational Mechanics_ , 66:1131–1152, 2020.  doi: 10.1007/s00466-020-01888-0.  URL https://doi.org/10.1007/s00466-020-01888-0.
  * Zohdi [2020] T. I. Zohdi.  An agent-based computational framework for simulation of global pandemic and social response on planet x.  _Computational Mechanics_ , 66:1195–1209, 2020.  doi: 10.1007/s00466-020-01886-2.  URL https://doi.org/10.1007/s00466-020-01886-2.
  * Chang et al. [2020] Lili Chang, Moran Duan, Guiquan Sun, and Zhen Jin.  Cross-diffusion-induced patterns in an sir epidemic model on complex networks.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 30(1):013147, 2020.  doi: 10.1063/1.5135069.  URL https://doi.org/10.1063/1.5135069.
  * Wang et al. [2020a] Zhenlin Wang, Xiaoxuan Zhang, GH Teichert, Mariana Carrasco-Teja, and Krishna Garikipati.  System inference for the spatio-temporal evolution of infectious diseases: Michigan in the time of COVID-19.  _Computational Mechanics_ , 66(5):1153–1176, 2020a.
  * [12] 1Point3Acres.com.  URL https://coronavirus.1point3acres.com/en.
  * Yang et al. [2020] Tong Yang, Kai Shen, Sixuan He, Enyu Li, Peter Sun, Pingying Chen, Lin Zuo, Jiayue Hu, Yiwen Mo, Weiwei Zhang, Haonan Zhang, Jingxue Chen, and Yu Guo.  Covidnet: To bring data transparency in the era of covid-19, 2020.
  * [14] Johns Hopkins University of Medicine.  COVID-19 dashboard by the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University (JHU).  URL https://coronavirus.jhu.edu/map.html.
  * [15] Michigan State Coronavirus Data.  URL https://www.michigan.gov/coronavirus/.
  * [16] The New York Times.  Coronavirus in the U.S.: Latest Map and Case Count - The New York Times.  URL https://www.nytimes.com/interactive/2020/us/coronavirus-us-cases.html.
  * [17] The Institute for Health Metrics and Evaluation.  COVID-19 Projections.  URL https://covid19.healthdata.org/united-states-of-america.
  * [18] Inegi: Censo de población y vivienda.  URL https://www.inegi.org.mx/programas/ccpv/2020/#Tabulados.
  * [19] Conacyt: Covid-19 méxico.  URL https://datos.covid-19.conacyt.mx/#DownZCSV.
  * Wang et al. [2019] Z. Wang, X. Huan, and K. Garikipati.  Variational system identification of the partial differential equations governing the physics of pattern-formation: Inference under varying fidelity and noise.  _Computer Methods in Applied Mechanics and Engineering_ , 356:44 – 74, 2019.  ISSN 0045-7825.  doi: https://doi.org/10.1016/j.cma.2019.07.007.
  * Wang et al. [2020b] Zhenlin Wang, Bowei Wu, Krishna Garikipati, and Xun Huan.  A perspective on regression and bayesian approaches for system identification of pattern formation dynamics.  _Theoretical and Applied Mechanics Letters_ , 10(3):188–194, 2020b.
  * Wang et al. [2021a] Z. Wang, X. Huan, and K. Garikipati.  Variational system identification of the partial differential equations governing microstructure evolution in materials: Inference over sparse and spatially unrelated data.  _Computer Methods in Applied Mechanics and Engineering_ , 377:113706, 2021a.  ISSN 0045-7825.  doi: https://doi.org/10.1016/j.cma.2021.113706.  URL https://www.sciencedirect.com/science/article/pii/S0045782521000426.
  * Wang et al. [2021b] Z. Wang, B. Martin, J. Weickenmeier, and K. Garikipati.  An inverse modelling study on the local volume changes during early morphoelastic growth of the fetal human brain.  _Brain Multiphysics_ , 2:100023, 2021b.  ISSN 2666-5220.  doi: https://doi.org/10.1016/j.brain.2021.100023.  URL https://www.sciencedirect.com/science/article/pii/S2666522021000034.
  * Wang et al. [2020c] Z. Wang, J.B. Estrada, E.M. Arruda, and K. Garikipati.  Discovery of deformation mechanisms and constitutive response of soft material surrogates of biological tissue by full-field characterization and data-driven variational system identification.  _to appear in Journal of the Mechanics and Physics of Solids_ , 2020c.  doi: 10.1101/2020.10.13.337964.  bioRxiv.
  * Teichert and Garikipati [2019] Gregory H Teichert and Krishna Garikipati.  Machine learning materials physics: Surrogate optimization and multi-fidelity algorithms predict precipitate morphology in an alternative to phase field dynamics.  _Computer Methods in Applied Mechanics and Engineering_ , 344:666–693, 2019.
  * Teichert et al. [2019] Gregory H Teichert, AR Natarajan, A Van der Ven, and Krishna Garikipati.  Machine learning materials physics: Integrable deep neural networks enable scale bridging by learning free energy functions.  _Computer Methods in Applied Mechanics and Engineering_ , 353:201 – 216, 2019.  ISSN 0045-7825.  doi: https://doi.org/10.1016/j.cma.2019.05.019.  URL http://www.sciencedirect.com/science/article/pii/S0045782519302889.
  * Teichert et al. [2020] GH Teichert, AR Natarajan, Anton Van der Ven, and Krishna Garikipati.  Scale bridging materials physics: Active learning workflows and integrable deep neural networks for free energy function representations in alloys.  _Computer Methods in Applied Mechanics and Engineering_ , 371:113281, 2020.
  * Zhang and Garikipati [2020] Xiaoxuan Zhang and Krishna Garikipati.  Machine learning materials physics: Multi-resolution neural networks learn the free energy and nonlinear elastic response of evolving microstructures.  _Computer Methods in Applied Mechanics and Engineering_ , 372:113362, 2020.  doi: https://doi.org/10.1016/j.cma.2020.113362.
  * Zhang and Garikipati [2021] Xiaoxuan Zhang and Krishna Garikipati.  Bayesian neural networks for weak solution of pdes with uncertainty quantification.  2021\.
  * Virtanen et al. [2020] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, CJ Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake Vand erPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1. 0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  _Nature Methods_ , 17:261–272, 2020.  doi: https://doi.org/10.1038/s41592-019-0686-2.
  * Mitusch et al. [2019] Sebastian K. Mitusch, Simon W. Funke, and Jørgen S. Dokken.  dolfin-adjoint 2018.1: automated adjoints for fenics and firedrake.  _Journal of Open Source Software_ , 4(38):1292, 2019.  doi: 10.21105/joss.01292.  URL https://doi.org/10.21105/joss.01292.
  * Chan and Saqib [2021] Eugene Y. Chan and Najam U. Saqib.  Privacy concerns can explain unwillingness to download and use contact tracing apps when covid-19 concerns are high.  _Computers in Human Behavior_ , 119:106718, 2021.  ISSN 0747-5632.  doi: https://doi.org/10.1016/j.chb.2021.106718.  URL https://www.sciencedirect.com/science/article/pii/S0747563221000406.
  * Kretzschmar et al. [2020] Mirjam E Kretzschmar, Ganna Rozhnova, Martin C J Bootsma, Michiel van Boven, Janneke H H M van de Wijgert, and Marc J M Bonten.  Impact of delays on effectiveness of contact tracing strategies for covid-19: a modelling study.  _The Lancet Public Health_ , 5(8):e452–e459, 2020\.  ISSN 2468-2667.  doi: https://doi.org/10.1016/S2468-2667(20)30157-2.  URL https://www.sciencedirect.com/science/article/pii/S2468266720301572.


# Stopping time property of thresholds of Storey-type FDR procedures

Xiongzhi Chen111Corresponding author. Center for Statistics and Machine
Learning, and Lewis-Sigler Institute for Integrative Genomics, Princeton
University, Princeton, NJ 08544. Email<EMAIL_ADDRESS>and Rebecca W.
Doerge222Departments of Statistics and Agronomy, Purdue University, West
Lafayette, IN 49707. Email<EMAIL_ADDRESS>
###### Abstract

For multiple testing, we introduce Storey-type FDR procedures and the concept
of “regular estimator of the proportion of true nulls”. We show that the
rejection threshold of a Storey-type FDR procedure is a stopping time with
respect to the backward filtration generated by the p-values and that a
Storey-type FDR estimator at this rejection threshold equals the pre-specified
FDR level, when the estimator of the proportion of true nulls is regular.
These results hold regardless of the dependence among or the types of
distributions of the p-values. They directly imply that a Storey-type FDR
procedure is conservative when the null p-values are independent and uniformly
distributed.
Keywords: False discovery rate, rejection threshold, stopping time property,
Storey-type FDR procedure.
MSC 2010 subject classifications: Primary 62C25, 62M99.

##  1 Introduction

To control false discovery rate (FDR, Benjamini and Hochberg, 1995) in
multiple testing, it is crucial to show the conservativeness of an FDR
procedure, i.e., that the FDR of an FDR procedure is no larger than a pre-
specified level. To show the conservativeness of their FDR procedures, Storey
et al. (2004) used the stopping time property (STP) of the rejection threshold
that was claimed in their Lemma 4, and Liang and Nettleton (2012) quoted this
lemma as their Lemma 5 and used it to prove their Theorem 7. They assumed that
the p-values corresponding to the true null hypotheses (i.e., null p-values)
are independent and uniformly distributed. However, none of them provided a
formal proof of the STP of the rejection thresholds, and without the STP the
optional stopping theorem (see, e.g., Karatzas and Shreve, 1991) can not be
applied and the martingale based approach in Storey et al. (2004) and Liang
and Nettleton (2012) fails.

In this article, we introduce Storey-type FDR procedure in Section 2 as an
extension of Storey’s procedure in Storey et al. (2004), introduce the concept
of regular estimator of the proportion of true nulls, and show in Section 3
that the STP of the rejection threshold of a Storey-type FDR procedure is
generic when the estimator of the proportion of true nulls is regular. This
provides a formal justification of the use of STP of the rejection thresholds
in Storey et al. (2004) and, if needed, in Chen and Doerge (2014). However, it
also implies that the rejection threshold of an adaptive FDR procedure may not
be a stopping time when the estimator of proportion true nulls is not regular.
Further, we show that a Storey-type FDR estimator at the threshold of its
corresponding Storey-type FDR procedure equals to the targeted FDR level. Some
consequences of these findings are given in Section 4, which include that a
Storey-type FDR procedure is always conservative when the null p-values are
independent and uniformly distributed. We end the article with a short
discussion in Section 5.

##  2 Storey-type FDR procedures

Let there be $m$ null hypotheses $H_{i}$ with associated p-values $p_{i}$ for
$i=1,...,m$, such that only $m_{0}$ among them are true nulls and the rest
false nulls. However, the true status of each $H_{i}$ is unknown, and so is
the proportion of true nulls $\pi_{0}=m_{0}/m$. A one-step multiple testing
procedure (MTP) claims that $H_{i}$ is a false null if and only if $p_{i}\leq
t$ using a rejection threshold $t\in\left[0,1\right]$. It gives
$R\left(t\right)$ as the total number of null hypotheses claimed to be false
and $V\left(t\right)$ as the number of true nulls claimed to be false. The FDR
(Benjamini and Hochberg, 1995) of such an MTP is defined as

$FDR\left(t\right)=\mathbb{E}\left[V\left(t\right)/\max\left\\{R\left(t\right),1\right\\}\right]\text{,}$

where $\mathbb{E}$ is the expectation.

Let $\tilde{\pi}_{0}\left(\lambda\right)$ with a tuning parameter
$\lambda\in[0,1)$ be an estimator of $\pi_{0}$ such that
$\tilde{\pi}_{0}\left(\lambda\right)$ ranges in $\left[0,1\right]$. We define
a “Storey-type FDR estimator” as

$\widetilde{FDR}_{\lambda}\left(t\right)=\min\left\\{1,\dfrac{\tilde{\pi}_{0}\left(\lambda\right)t}{m^{-1}\max\left\\{R\left(t\right),1\right\\}}\right\\}\text{\quad
for \quad}t\in\left[0,1\right].$ (1)

For $\alpha\in\left[0,1\right]$, let

$t_{\alpha}\left(\widetilde{FDR}_{\lambda}\right)=\sup\left\\{t\in\left[0,1\right]:\widetilde{FDR}_{\lambda}\left(t\right)\leq\alpha\right\\}$
(2)

and the decision rule based on $\widetilde{FDR}_{\lambda}$

$\text{claim }H_{i}\text{ as a false null }\Longleftrightarrow p_{i}\leq
t_{\alpha}\left(\widetilde{FDR}_{\lambda}\right).$ (3)

The procedure defined by (1) and (3) is called a “Storey-type FDR procedure”.

A Storey-type FDR procedure allows the use of various estimators of $\pi_{0}$
and is very versatile. For example, (i) it is Storey’s procedure when
$\tilde{\pi}_{0}\left(\lambda\right)$ in (1) is the estimator

$\hat{\pi}_{0}\left(\lambda\right)=\min\left\\{1,\left(1-\lambda\right)^{-1}m^{-1}\sum\nolimits_{i=1}^{m}\mathbf{1}_{\left\\{p_{i}>\lambda\right\\}}\right\\}$
(4)

in Storey et al. (2004); (ii) it is the generalized FDR procedure in Chen and
Doerge (2014), designed for multiple testing based on discrete p-values, when
$\tilde{\pi}_{0}\left(\lambda\right)$ is the estimator of $\pi_{0}$ proposed
there and further studied in Chen and Doerge (2015); (iii) it is the BH
procedure in Benjamini and Hochberg (1995) when $\lambda=0$ is set in (4) in
Storey’s procedure or $\tilde{\pi}_{0}\left(\lambda\right)\equiv 1$ is set;
(iv) it is a “dynamic adaptive FDR procedure” studied in Liang and Nettleton
(2012) when $\lambda$ in (4) is determined from the p-values and can be a
random variable.

##  3 Stopping time property of the rejection threshold

For notational simplicity, we will from now on write
$t_{\alpha}\left(\widetilde{FDR}_{\lambda}\right)$ as
$\tilde{t}_{\alpha}\left(\lambda\right)$. Define the backward filtration

$\mathcal{F}_{t}=\sigma\left(\mathbf{1}_{\left\\{p_{i}\leq s\right\\}},t\leq
s\leq 1,i=1,...,m\right)$ (5)

for each $t\in[0,1)$ and the “stopped backward” filtration
$\mathcal{G}=\left\\{\mathcal{F}_{t\wedge\lambda}:1\geq t\geq 0\right\\}$.
Further, introduce the following

######  Definition 1.

An estimator $\tilde{\pi}_{0}\left(\lambda\right)$ of the proportion $\pi_{0}$
of true nulls with tuning parameter $\lambda\in[0,1)$ is said to be regular if
it is measurable with respect to (wrt) $\mathcal{F}_{\lambda}$.

It is crucial to note that, if $\lambda$ is a functional of the p-values such
that the information contained in $\mathcal{F}_{\lambda}$ is not sufficient to
determine the value of $\tilde{\pi}_{0}\left(\lambda\right)$, then
$\tilde{\pi}_{0}\left(\lambda\right)$ is not measurable wrt
$\mathcal{F}_{\lambda}$. This can happen for dynamic Storey procedures
considered in Liang and Nettleton (2012), where $\lambda$ is adaptively
determined from the data. When this happens,
$\tilde{t}_{\alpha}\left(\lambda\right)$ does not have to be a stopping time
wrt to $\mathcal{G}$, the optional stopping theorem can not be applied, and
the martingale arguments to prove the conservativeness of (dynamic) Storey FDR
procedure in Storey et al. (2004) and Liang and Nettleton (2012) are invalid.
Here is out main result:

######  Theorem 1.

If $\tilde{\pi}_{0}\left(\lambda\right)$ is regular, then
$\tilde{t}_{\alpha}\left(\lambda\right)$ is a stopping time with respect to
$\mathcal{G}$. Further, if $\tilde{t}_{\alpha}\left(\lambda\right)<1$, then

$\widetilde{FDR}_{\lambda}\left(\tilde{t}_{\alpha}\left(\lambda\right)\right)=\alpha.$
(6)

The proof of Theorem 1 is given in Appendix A. Note that
$\tilde{t}_{\alpha}\left(\lambda\right)<1$ excludes the meaningless case
$\alpha=1$ as the targeted FDR level. Therefore,
$\tilde{t}_{\alpha}\left(\lambda\right)<1$ in Theorem 1 essentially is not an
assumption. Theorem 1 shows that, regardless of whether the p-values are
independent, whether they are continuously distributed, or the stochastic
orders of their distributions wrt the standard uniform random variable, the
rejection threshold of a Storey-type FDR procedure is a stopping time wrt to
the stopped backward filtration, when the estimator of the proportion of true
nulls is regular. It also shows that a Storey-type FDR estimator equals the
pre-specified FDR level at the threshold of its corresponding FDR procedure,
when the estimator of the proportion of true nulls is regular and the
corresponding decision rule is statistically meaningfully implemented. Equally
importantly, Theorem 1 shows that the STP may not hold for the rejection
threshold of a Storey-type procedure when
$\tilde{\pi}_{0}\left(\lambda\right)$ is not regular, i.e., when adaptivity of
the estimator $\tilde{\pi}_{0}\left(\lambda\right)$ to $\mathcal{F}_{\lambda}$
is not ensured.

We make four remarks on Theorem 1 and its proof: (i) the measurability of
$\tilde{\pi}_{0}\left(\lambda\right)$ wrt to $\mathcal{F}_{\lambda}$ is
critical to the STP of $\tilde{t}_{\alpha}\left(\lambda\right)$; (ii) the fact
that (6), i.e.,
$\widetilde{FDR}_{\lambda}\left(\tilde{t}_{\alpha}\left(\lambda\right)\right)=\alpha$
for $\tilde{t}_{\alpha}\left(\lambda\right)<1$, holds for general p-values is
a new result, whose consequences will be discussed in Section 4; (iii) the
validity of the STP of $\tilde{t}_{\alpha}\left(\lambda\right)$, regardless of
the joint distribution of the p-values, depends crucially on the intrinsic
structure of the backward filtration $\mathcal{G}$; (iv) the STP of the
threshold (wrt a similar backward filtration) defined by equation (7.2) in
Pena et al. (2011) does hold trivially due to the special structure of the
filtration there.

##  4 Conservativeness of Storey-type FDR procedures

We give a few implications of Theorem 1. For concise statements, we set
assumption A1) as “the null p-values are independent and uniformly
distributed”. Further, for each $t\in[0,1)$ let
$\mathcal{H}_{t}=\sigma\left(\mathbf{1}_{\left\\{p_{i}\leq s\right\\}},0\leq
s\leq t,i=1,...,m\right)$.

######  Corollary 1.

Assume A1), the following hold:

  1. 1.

Lemma 4 of Storey et al. (2004) is valid when
$\hat{\pi}_{0}\left(\lambda\right)$ in (4) is regular.

  2. 2.

Storey’s procedure in Storey et al. (2004) is conservative when
$\hat{\pi}_{0}\left(\lambda\right)$ in (4) is regular and
$\tilde{t}_{\alpha}\left(\lambda\right)<1$.

  3. 3.

If $\lambda$ is a stopping time wrt to
$\left\\{\mathcal{H}_{t}:t\in[0,1)\right\\}$ and lies in a compact subset
$\left[\kappa,\tau\right]$ in $\left(0,1\right)$ for some $\kappa<\tau$ and
$\tilde{\pi}_{0}\left(\lambda\right)$ is regular, then a Storey-type FDR
procedure is conservative when $\tilde{t}_{\alpha}\leq\kappa$.

###### Proof.

The first claim is a direct consequence of Theorem 1. We show the second
claim. By Corollary 1 in Liang and Nettleton (2012),
$\mathbb{E}\left[\widetilde{FDR}_{\lambda}\left(\tilde{t}_{\alpha}\right)\right]\geq
FDR\left(\tilde{t}_{\alpha}\right)$. However, Theorem 1 implies
$\mathbb{E}\left[\widetilde{FDR}_{\lambda}\left(\tilde{t}_{\alpha}\right)\right]=\alpha$.
So, $FDR\left(\tilde{t}_{\alpha}\right)\leq\alpha$. Finally, we show the
third. By Theorem 4 in Liang and Nettleton (2012),
$\mathbb{E}\left[\widetilde{FDR}_{\lambda}\left(\tilde{t}_{\alpha}\right)\right]\geq
FDR\left(\tilde{t}_{\alpha}\right)$. Therefore, from Theorem 1 we see
$FDR\left(\tilde{t}_{\alpha}\right)\leq\alpha$. This completes the proof. ∎

Corollary 1 illustrates the connection between the regularity of the estimator
of the proportion of true nulls, the STP of the rejection threshold, and the
conservativeness of a Storey-type FDR procedure. The third conclusion in
Corollary 1 allows to choose $\lambda$ in
$\tilde{\pi}_{0}\left(\lambda\right)$ from the data in a Storey-type FDR
procedure without potentially sacrificing its conservativeness. Further, it
covers the dynamic Storey’s FDR procedure in Liang and Nettleton (2012) where
$\hat{\pi}_{0}\left(\lambda\right)$ in (4) is used and $\lambda$ is determined
from the data.

Compared to the proofs of the conservativeness of (dynamic) Storey’s procedure
in Storey et al. (2004) and Liang and Nettleton (2012), the proof of Corollary
1 we provide based on Theorem 1 is much shorter and the second conditioning
step there is no longer needed. This illustrates the usefulness of our finding
(6) for general p-values in the martingale arguments to prove the
conservativeness of a Storey-type FDR procedure.

##  5 Discussion

We have shown the STP of the rejection threshold of a Storey-type FDR
procedure for general p-values, that a Storey-type FDR estimator equals the
targeted FDR level at its rejection threshold, and that a Storey-type FDR
procedure is usually conservative, when the estimator of the proportion of
true nulls is regular. This implies that the STP may not hold for the
rejection threshold of an adaptive FDR procedure when an estimator of the
proportion of true nulls is not adaptive to the backward filtration.
Therefore, caution should be taken before applying martingale arguments to
prove the conservativeness of an adaptive FDR procedure. In view of the
increasing use of stochastic processes (such as decision processes) and
suprema related to these processes (such as rejection thresholds) in multiple
testing, to better understand the behavior of these suprema remains an
important task.

## Acknowledgements

This research was funded by a National Science Foundation Plant Genome
Research Program grant (No. IOS-1025976) to RWD. Its major part was completed
when X. Chen was a PhD student at the Department of Statistics, Purdue
University.

##  Appendix A Proofs

In order to prove the STP of the rejection threshold, we first need to
understand the the scaled inverse rejection process
$L\left(t\right)=t\left(\max\left\\{R\left(t\right),1\right\\}\right)^{-1}$
with $t\in\left[0,1\right]$, with which the regularity of
$\tilde{\pi}_{0}\left(\lambda\right)$, the property of the stopped backward
filtration $\mathcal{G}$ and a contrapositive argument will yield Theorem 1.
In what follows, no assumption will be made about the independence between the
p-values, the continuity of the p-value distributions, or their stochastic
orders wrt the standard uniform distribution.

###  A.1 Downward jumps of the scaled inverse rejection process

Order the p-values into
$p_{\left(1\right)}<p_{\left(2\right)}<...<p_{\left(n\right)}$ distinctly,
where the multiplicity of $p_{\left(i\right)}$ is $n_{i}$ for $i=1,...,n$. Let
$p_{\left(n+1\right)}=\max\left\\{p_{\left(n\right)},1\right\\}$ and
$p_{\left(0\right)}=0.$ Define $T_{j}=\sum\limits_{l=1}^{j}n_{l}$ for
$j=1,...,n$.

######  Lemma A.1.

The process $\left\\{L\left(t\right),t\in[0,1]\right\\}$ is such that

$L\left(t\right)=\left\\{\begin{array}[]{lcl}t&\text{if}&t\in[0,p_{\left(1\right)})\text{,}\\\
tT_{j}^{-1}&\text{if}&t\in[p_{\left(j\right)},p_{\left(j+1\right)})\text{ for
}j=1,..,n-1\text{,}\\\
tm^{-1}&\text{if}&t\in[p_{\left(n\right)},p_{\left(n+1\right)}]\text{.}\end{array}\right.$
(A.1)

Moreover, it can only be discontinuous at $p_{\left(i\right)}$, $1\leq i\leq
n$, where it can only have a downward jump with size

$L\left(p_{\left(i\right)}-\right)-L\left(p_{\left(i\right)}\right)=\dfrac{p_{\left(i\right)}n_{i}}{R\left(p_{\left(i\right)}\right)\left[R\left(p_{\left(i\right)}\right)-n_{i}\right]}>0\text{.}$

###### Proof.

Clearly

$R\left(t\right)=\left\\{\begin{array}[]{ccl}0&\text{if}&0\leq
t<p_{\left(1\right)}\text{,}\\\ T_{j}&\text{if}&p_{\left(j\right)}\leq
t<p_{\left(j+1\right)},j=1,..,n-1\text{,}\\\
m&\text{if}&p_{\left(n\right)}\leq t\leq
p_{\left(n+1\right)}\text{,}\end{array}\right.$

and (A.1) holds. Therefore, the points of discontinuities of
$L\left(\cdot\right)$ are the original distinct p-values. This justifies the
first part of the assertion.

Now we show that $L\left(\cdot\right)$ can only have downward jumps at points
of discontinuity. Define

$\varphi\left(t,\eta\right)=\dfrac{t+\eta}{R\left(t+\eta\right)}-\dfrac{t}{R\left(t\right)}=\dfrac{\eta
R\left(t\right)+t\left[R\left(t\right)-R\left(t+\eta\right)\right]}{R\left(t+\eta\right)R\left(t\right)}\text{.}$

From the fact
$R\left(p_{\left(j\right)}\right)-R\left(p_{\left(j\right)}-\right)=n_{j}>0$
but $R\left(p_{\left(j\right)}+\right)-R\left(p_{\left(j\right)}\right)=0$ for
each $1\leq j\leq n,$ it follows that
$\varphi\left(p_{\left(j\right)},0+\right)=\lim_{\eta\downarrow
0}\varphi\left(t,\eta\right)=0$ but

$\varphi\left(p_{\left(j\right)},0-\right)=\lim_{\eta\uparrow
0}\varphi\left(p_{\left(j\right)},\eta\right)=\dfrac{p_{\left(j\right)}n_{j}}{R\left(p_{\left(j\right)}\right)\left[R\left(p_{\left(j\right)}\right)-n_{j}\right]}>0\text{.}$

Thus
$L\left(p_{\left(j\right)}-\right)-L\left(p_{\left(j\right)}\right)=\varphi\left(p_{\left(j\right)},0-\right)>0$
and the proof is completed. ∎

Lemma A.1 shows that the process $L\left(t\right)$ is piecewise linear and can
only have downward jumps as $t$ increases. The conclusion of Lemma A.1 is
right the contrary to the claim in the proof of Theorem 2 in Storey et al.
(2004) that “the process $mt/R\left(t\right)$ has only upward jumps and has a
final value of $1$”, since it says “the process $mt/R\left(t\right)$ has only
downward jumps”. We construct a counterexample to their claim as follows.

For a small increase $c$ in $t$ which results an increase $a_{c}$ in
$R\left(t\right)$, we see that

$L\left(t+c\right)-L\left(t\right)=\dfrac{t+c}{R\left(t\right)+a_{c}}-\dfrac{t}{R\left(t\right)}=\dfrac{cR\left(t\right)-ta_{c}}{\left[R\left(t\right)+a_{c}\right]R\left(t\right)}<0$

if and only if $\dfrac{c}{a_{c}}<\dfrac{t}{R\left(t\right)}$. Construct $m$
p-values with $n\geq 4$ such that there exists some $1\leq j_{0}<n-2$ with
$n_{j_{0}+1}>T_{j_{0}}$ but $p_{\left(j_{0}+1\right)}<1$. Choose $c_{1}$ and
$c_{2}$ such that
$0<c_{1}<\left(p_{\left(j_{0}+1\right)}-p_{\left(j_{0}\right)}\right)/2$ and
$0<c_{2}<p_{\left(j_{0}+1\right)}-2c_{1}$. Let
$t_{0}=p_{\left(j_{0}+1\right)}-c_{1}$ and $c=c_{1}+c_{2}$. Then
$p_{\left(j_{0}\right)}<t_{0}<p_{\left(j_{0}+1\right)}$,
$R\left(t_{0}\right)=T_{j_{0}}$ and
$R\left(t_{0}+c\right)=T_{j_{0}}+n_{j_{0}+1}$. Further, $0<c<t_{0}$ and
$a_{c}=n_{j_{0}+1}$. So $\dfrac{c}{n_{j_{0}+1}}<\dfrac{t_{0}}{T_{j_{0}}}$ and
$L\left(t_{0}+c\right)-L\left(t_{0}\right)<0$. Letting $c\rightarrow 0$ gives
$p_{\left(j_{0}+1\right)}$ as a point of downward jump for $L\left(t\right)$.

###  A.2 Proof of Theorem 1

When $\alpha=0$, $\tilde{t}_{\alpha}\left(\lambda\right)=0$ if
$\tilde{\pi}_{0}\left(\lambda\right)>0$ but
$\tilde{t}_{\alpha}\left(\lambda\right)=1$ if
$\tilde{\pi}_{0}\left(\lambda\right)=0$. In this case,
$\tilde{t}_{\alpha}\left(\lambda\right)$ is already a stopping time and
$\widetilde{FDR}_{\lambda}\left(\tilde{t}_{\alpha}\left(\lambda\right)\right)=\alpha=0$.
On the other hand, $\tilde{t}_{\alpha}\left(\lambda\right)=1$ when $\alpha=1$.
In this case, $\tilde{t}_{\alpha}\left(\lambda\right)$ is a stopping time and
$\widetilde{FDR}_{\lambda}\left(\tilde{t}_{\alpha}\left(\lambda\right)\right)=\tilde{\pi}_{0}\left(\lambda\right)$.
Note that $\tilde{\pi}_{0}\left(\lambda\right)=0$ can not happen when
$0<\tilde{t}_{\alpha}\left(\lambda\right)<1$.

First, we show the claim on the stopping time property. By the previous
discussion, we only need to consider the case where
$0<\tilde{t}_{\alpha}\left(\lambda\right)<1$ and
$\tilde{\pi}_{0}\left(\lambda\right)>0$. Define

$\mathbf{X}_{t}^{\left(m\right)}\left(\omega\right)=\left(\mathbf{1}_{\left\\{p_{1}\leq
t\right\\}}\left(\omega\right),...,\mathbf{1}_{\left\\{p_{m}\leq
t\right\\}}\left(\omega\right)\right)\text{, }\omega\in\Omega\text{,}$

where $\left(\Omega,\mathcal{A},\mathbb{P}\right)$ is the probability space
with $\Omega$ the sample space, $\mathcal{A}$ the sigma-algebra, and
$\mathbb{P}$ the probability measure. Then
$\mathcal{F}_{t}=\sigma\left(\mathbf{X}_{s}^{\left(m\right)}\left(\omega\right),1\geq
s\geq t\right)$, $t\in\left[0,1\right]$ and $\left\\{\mathcal{F}_{t}:0\leq
t\leq 1\right\\}$ is a non-increasing sequence of sub-sigma-algebras of
$\mathcal{A}$. Write
$\widetilde{FDR}_{\lambda}\left(t\right)=\widetilde{FDR}_{\lambda}\left(t,\omega\right)$.

By the definition of $\tilde{t}_{\alpha}\left(\lambda\right)$,

$\left\\{\omega\in\Omega:\tilde{t}_{\alpha}\left(\lambda\right)\leq
s\right\\}=\bigcap\limits_{\left\\{t:s<t\leq
1\right\\}}A_{t}=\tilde{A}_{s}\text{,}$ (A.2)

where
$A_{t}=\left\\{\omega\in\Omega:\widetilde{FDR}_{\lambda}\left(t,\omega\right)>\alpha\right\\}$,
$1\geq t>s$. Thus, we only need to show
$\tilde{A}_{s}\in\mathcal{F}_{s\wedge\lambda}$. Since either
$\mathcal{F}_{s\wedge\lambda}=\mathcal{F}_{s}\supseteq\mathcal{F}_{\lambda}$
when $s\leq\lambda$ or
$\mathcal{F}_{s\wedge\lambda}=\mathcal{F}_{\lambda}\supseteq\mathcal{F}_{s}$
when $s\geq\lambda$, the stopping time property holds once we prove
$\tilde{A}_{s}\in\mathcal{F}_{s}$. Let $\mathbb{Q}$ be the set of all rational
numbers. Since $0<\alpha<1$ and $\tilde{\pi}_{0}\left(\lambda\right)>0$, the
following decompositions are valid:
$A_{t}=\bigcup\nolimits_{r\in\mathbb{Q}}\left(A_{t,r}\cap B_{r}\right)$ and

$\tilde{A}_{s}=\bigcap\limits_{\left\\{t:s<t\leq
1\right\\}}\bigcup\limits_{r\in\mathbb{Q}}\left(A_{t,r}\cap
B_{r}\right)=\bigcup\limits_{r\in\mathbb{Q}}\bigcap\limits_{\left\\{t:s<t\leq
1\right\\}}\left(A_{t,r}\cap B_{r}\right)\text{,}$

where
$A_{t,r}=\left\\{\omega\in\Omega:\dfrac{t}{m^{-1}\max\left\\{R\left(t\right),1\right\\}}\geq
r\right\\}$ and
$B_{r}=\left\\{\omega\in\Omega:\dfrac{\alpha}{\tilde{\pi}_{0}\left(\lambda\right)}<r\right\\}$.
Thus, it suffices to show

$\bigcap\limits_{\left\\{t:s<t\leq 1\right\\}}\left(A_{t,r}\cap
B_{r}\right)=\bigcap\limits_{\left\\{t:s<t\leq 1\right\\}}A_{t,r}\cap
B_{r}\in\mathcal{F}_{s}\text{,}$ (A.3)

We now move to show

$A_{s,r}^{\ast}=\bigcap\limits_{\left\\{t:s<t\leq
1\right\\}}A_{t,r}\in\mathcal{F}_{s}\text{.}$ (A.4)

Define $I_{i}=[p_{\left(i\right)},p_{\left(i+1\right)})$, $i=1,...,n-1$. We
will add $I_{0}=[p_{\left(0\right)},p_{\left(1\right)})$ and
$I_{n}=[p_{\left(n\right)},p_{\left(n+1\right)}]$ when $p_{\left(n\right)}<1$.
When $p_{\left(n\right)}=1$, we take $I_{n-1}$ to be
$\left[p_{\left(n-1\right)},p_{\left(n\right)}\right]$. Obviously there must
be a unique $j^{\ast}$ with $0\leq j^{\ast}\leq n$ such that $s\in
I_{j^{\ast}}$. Given $R\left(1\right)=m$ and
$\tilde{\pi}_{0}\left(\lambda\right)\in\left[0,1\right]$, the properties of
$L\left(\cdot\right)$ in Lemma A.1 imply

$A_{s,r}^{\ast}=A_{s,r}\bigcap\left(\bigcap\nolimits_{j=j^{\ast}+1}^{n+1}A_{p_{\left(j\right)},r}\right)\text{.}$

Consequently, $A_{s,r}^{\ast}\in\mathcal{F}_{s}$, i.e., (A.4) holds.

Recall $\tilde{A}_{s}$ defined in (A.2). If $B_{r}\in\mathcal{F}_{\lambda}$,
then (A.4) implies (A.3), which implies $\tilde{A}_{s}\in\mathcal{F}_{s}$,
i.e., $\tilde{t}_{\alpha}\left(\lambda\right)$ is a STP wrt to $\mathcal{G}$
if $\tilde{\pi}_{0}\left(\lambda\right)$ is measurable wrt to
$\mathcal{F}_{\lambda}$. However, this holds since
$\tilde{\pi}_{0}\left(\lambda\right)$ is regular.

Finally, we show
$\widetilde{FDR}_{\lambda}\left(\tilde{t}_{\alpha}\right)=\alpha$, where we
have written $\tilde{t}_{\alpha}\left(\lambda\right)$ as $\tilde{t}_{\alpha}$.
Clearly, $\widetilde{FDR}_{\lambda}\left(\tilde{t}_{\alpha}\right)\leq\alpha$
by the definition of $\tilde{t}_{\alpha}$. Our arguments next proceed by
contrapositive reasoning, i.e., that
$\widetilde{FDR}_{\lambda}\left(\tilde{t}_{\alpha}\right)<\alpha$ gives a
contradiction. Obviously, there must be a unique $0\leq j^{\prime}\leq n$ such
that $\tilde{t}_{\alpha}\in I_{j^{\prime}}$. Since $\tilde{t}_{\alpha}<1$, we
have $\alpha<1$, $\tilde{\pi}_{0}\left(\lambda\right)>0$,
$p_{\left(j^{\prime}+1\right)}>0$ and
$I^{\ast}=\left[\tilde{t}_{\alpha},p_{\left(j^{\prime}+1\right)}\right)\subseteq$
$I_{j^{\prime}}$. Let $\rho_{m}^{\ast}\left(t\right)=m^{-1}L\left(t\right)$.
Then $\rho_{m}^{\ast}\left(t\right)$ is continuous and strictly increasing on
$I_{j^{\prime}}$ by Lemma A.1. If
$\widetilde{FDR}_{\lambda}\left(\tilde{t}_{\alpha}\right)<\alpha$ and
$\alpha<1$, then there must be some $d^{\prime}>0$ such that
$\tilde{I}^{\ast}=\left[\tilde{t}_{\alpha},\tilde{t}_{\alpha}+d^{\prime}\right]\subseteq
I^{\ast}$ and that

$\widetilde{FDR}_{\lambda}\left(t\right)=\tilde{\pi}_{0}\left(\lambda\right)\rho_{m}^{\ast}\left(t\right)\text{\quad
for all \quad}t\in\tilde{I}^{\ast}.$ (A.5)

However, (A.5) implies that there exists $\hat{t}_{\alpha}\in\tilde{I}^{\ast}$
such that $\hat{t}_{\alpha}>\tilde{t}_{\alpha}$ but

$\widetilde{FDR}_{\lambda}\left(\hat{t}_{\alpha}\right)=\tilde{\pi}_{0}\left(\lambda\right)\rho_{m}^{\ast}\left(\hat{t}_{\alpha}\right)<\alpha.$

This contradicts the definition of $\tilde{t}_{\alpha}.$ Hence
$\widetilde{FDR}_{\lambda}\left(\tilde{t}_{\alpha}\right)=\alpha$ must hold,
which completes the whole proof.

## References

  * Benjamini and Hochberg (1995) Benjamini, Y. and Y. Hochberg (1995).  Controlling the false discovery rate: a practical and powerful approach to multiple testing.  J. R. Statist. Soc. Ser. B 57(1), 289–300.
  * Chen and Doerge (2014) Chen, X. and R. Doerge (2014).  Generalized estimators formultiple testing: proportion of true nulls and false discovery rate.  http://arxiv.org/abs/1410.4274.
  * Chen and Doerge (2015) Chen, X. and R. Doerge (2015).  A weighted fdr procedure under discrete and heterogeneous null distributions.  http://arxiv.org/abs/1502.00973v2.
  * Karatzas and Shreve (1991) Karatzas, I. and S. E. Shreve (1991).  Brownian Motion and Stochastic Calculus (2 ed.).  Graduate Texts in Mathematics. Springer Verlag.
  * Liang and Nettleton (2012) Liang, K. and D. Nettleton (2012).  Adaptive and dynamic adaptive procedures for false discovery rate control and estimation.  J. R. Statist. Soc. Ser. B 74(1), 163–182.
  * Pena et al. (2011) Pena, E. A., J. D. Habiger, and W. Wu (2011).  Power-enhanced multiple decision functions controlling family-wise error and false discovery rates.  Ann. Statist. 39(1), 556–583.
  * Storey et al. (2004) Storey, J. D., J. E. Taylor, and D. Siegmund (2004).  Strong control, conservative point estimation in simultaneous conservative consistency of false discover rates: a unified approach.  J. R. Statist. Soc. Ser. B 66(1), 187–205.


# The van der Waals idea of pseudo associations and the critical
compressibility factor

O. Bakai E-mail<EMAIL_ADDRESS>
(Received October 2, 2019, in final form October 31, 2019)

###### Abstract

The dimensionless value of critical compressibility factor in the van der
Waals theory of gas-liquid critical point is a universal constant,
$Z_{\text{c}}=0.375$. Experimentally measured values of this quantity for
simple fluids are considerably smaller than the theory prediction. Van der
Waals once assumed that this discrepancy can be removed by taking account of
the impact of the molecular pseudo-associations on the fluid criticality but
he did not complete a proper modification of his theory following up on this
idea. The communication is devoted to the filling of this gap.

Key words: van der Waals, critical compressibility factor, heterophase
fluctuations

###### Abstract

Áåçðîçìðíèé êîåôöíò êðèòèíî¿ ñòèñëèâîñò â òåîð¿ Âàí äåð Âààëüñà êðèòèíî¿ òîêè
ãàç-ðäèíà óíâåðñàëüíîþ êîíñòàíòîþ, $Z_{\text{c}}=0.375$. Åêñïåðèìåíòàëüíî
âèìðÿíå çíàåííÿ ö¿ âåëèèíè äëÿ ïðîñòèõ ïëèíâ çíàíî ìåíøèì, íæ òåîðåòèíå
ïåðåäáàåííÿ. Âàí äåð Âààëüñ çðîáèâ ïåðåäáàåííÿ, ùî öÿ ðîçáæíñòü ìîæå áóòè
óñóíóòà, ÿêùî âðàõóâàòè âïëèâ ìîëåêóëÿðíî¿ ïñåâäîàñîöàö¿ íà êðèòèíñòü ïëèíó.
Ïðîòå, âí íå çàâåðøèâ ìîäèôêàö¿ ñâî¿ òåîð¿ âäïîâäíî äî ö¿ äå¿. Öÿ ðîáîòà
ïðèñâÿåíà çàïîâíåííþ ö¿ ïðîãàëèíè.

Ключов слова: Âàí äåð Âààëüñ, êîåôöíò êðèòèíî¿ ñòèñëèâîñò, ãåòåðîôàçí
ôëóêòóàö¿

##  1 Introduction

The van der Waals equation of state represented in the form of the law of
corresponding states contains a universal dimensionless ratio, the critical
compressibility factor,

$Z_{\text{c}}=\frac{P_{\text{c}}v_{\text{c}}}{T_{\text{c}}}=0.375,$ (1)

where $P$, $T$ and $v$ is the pressure, temperature and specific volume. Here,
index “c” denotes quantities evaluated at the critical point.

In his famous Nobel prize lecture, discussing the discrepancy between the
predicted universal value of the critical compressibility factor (1) and its
experimentally measured value for CO2 (which is nearly 30% less), van der
Waals assumed that it is necessary to take account of the formation of the
transient molecular complexes — “pseudo associations” — within the dense fluid
in order to improve the equation of state and the critical compressibility
factor. The pseudo-associations were described as follows [1].

“Let the number of molecules that have combined into a complex be so large
that it is possible to speak of a molecule at the centre surrounded by a
single layer containing almost as many other molecules as is possible
simultaneously. Then, for the surrounding molecules the attraction directed
towards the interior acts only to maintain the complex; and this part of its
attraction is lost for the surface pressure. Only the forces acting outwards
from these molecules can contribute to the formation of the internal pressure.
But of course, for pseudo association as for true association the number of
formed complexes increases with decreasing temperature and volume. At the
critical point, so I was compelled to conclude, only a very small part of the
weight is present as complexes. If pseudo association exists in a substance,
there are at least two types of molecules, namely simple and complex.”

The formulated problem was not properly solved by van der Waals. It is
revisited in the present communication.

##  2 The van der Waals equation of state taking account of the heterophase
fluctuations

In fact, the problem formulated by van der Waals was partially solved by
Frenkel in the theory of heterophase fluctuations (the liquid-like droplets
within the gas and gaseous bubbles within the liquid) which he formulated
considering the pre-transition anomalies of the coexisting gas and liquid far
below the critical point [2]. Identifying the liquid-like droplets within the
gas as the pseudo-associates, one obtains an improved equation of state.
However, Frenkel’s theory is applicable only in the case of weak heterophase
fluctuations when the droplets concentration is so small that they do not
interact with each other (the droplet ideal gas). Frenkel has described the
impact of the heterophase fluctuations on the thermodynamic quantities far
from the critical point.

Contrary to the van der Waals expectation (see above), near the critical point
the concentration of droplets is large, and their interaction cannot be
ignored (see [3] and references therein). In this regime, the coexisting
nominal gas and nominal liquid phases are essentially the states of a
heterophase fluid with comparable values of the liquid-like and gas-like
fractions (l- and g-fractions). Elementary constituents of the l- and
g-fractions are transient droplets and bubbles which have a mesoscopic scale
comparable with the range of short-range order in the fluid. These entities
can be considered as statistically independent species of the fluid. They are
called l- and g-fluctuons, respectively.

The nominal gaseous and liquid phases near the critical point are heterophase
mixtures of mesoscopic bubbles and droplets. Interactions between heterophase
fluctuations lead to a modification of the van der Waals equation of state.

The specific volume of the heterophase fluid,

$v=c_{\text{g}}v_{\text{g}}+c_{\text{l}}v_{\text{l}}\,,$ (2)

is determined by the concentrations of molecules belonging to the bubbles and
droplets, $c_{\text{g}}$ and $c_{\text{l}}$ respectively. Here, $v_{\text{g}}$
and $v_{\text{l}}$ is the specific volume in the bubble and droplet,
respectively.

Let us assume that the specific volume (2) obeys the van der Waals equation of
state while the l- and g-fluctuons are elementary disturbances.

The free energy per molecule of a homogeneous one-component simple fluid in
the van der Waals approximation has a standard form,

$\frac{F_{\text{VdW}}(v,T)}{N}=-T\ln\left(\frac{v-b}{\Lambda^{3}}\right)-\frac{a}{v}\,,$
(3)

where $\Lambda$ is the thermal de Broglie wavelength, N is the number of
atoms, $b$ is the excluded volume, and $a$ is the virial coefficient.

To include the contribution of the heterophase fluctuations to the free energy
at the coexistence curve in a self-consistent manner we have to require that
the chemical potentials of molecules within the l- and g-fluctuons are equal
to the chemical potential of the van der Waals fluid. With this condition, the
free energy of the fluid gets contributions due to the interaction energy,
$\varepsilon_{\text{fl,int}}$, and the mixing entropy, $s_{\text{fl,mix}}$, of
the l- and g-fluctuons,

$\frac{F_{\text{fl}}(v,T,\sigma_{\text{l}})}{N}=f_{\text{fl}}(v,T,\sigma_{\text{l}})=\frac{N_{\text{fl}}}{N}\left(\varepsilon_{\text{fl,{int}}}-Ts_{\text{fl,{mix}}}\right)=\frac{1}{\bar{k}_{\text{fl}}}\left[g_{2}\sigma_{\text{l}}\sigma_{\text{g}}+T\left(\sigma_{\text{l}}\ln\sigma_{\text{l}}+\sigma_{\text{g}}\ln\sigma_{\text{g}}\right)\right],$
(4)

where index “fl” denotes the quantities describing the fluctuons,
$N_{\text{fl}}$ is the number of fluctuons, $\bar{k}_{\text{fl}}$ is the mean
number of atoms per fluctuon, $\sigma_{\text{l}}$ and $\sigma_{\text{g}}$ are
fractions of the l- and g-fluctuons ($\sigma_{\text{l}}+\sigma_{\text{g}}=1$),
and parameter $g_{2}>0$ quantifies the fluctuonic interaction.

The free energy accounting for the fluctuonic contribution is given by

$\displaystyle\frac{F(v,T,\sigma_{\text{l}})}{N}=-T\ln\left(\frac{v-b}{\Lambda^{3}}\right)-\frac{a}{v}+\frac{1}{\bar{k}_{\text{fl}}}\left[g_{2}\sigma_{\text{l}}\sigma_{\text{g}}+T\left(\sigma_{\text{l}}\ln\sigma_{\text{l}}+\sigma_{\text{g}}\ln\sigma_{\text{g}}\right)\right].$
(5)

The number of molecules per l- and g-fluctuon, $k_{\text{l}}$ and
$k_{\text{g}}$, determine the mean number of molecules per fluctuon,
$\bar{k}_{\text{fl}}$,

$\bar{k}_{\text{fl}}=\sigma_{\text{l}}k_{\text{l}}+\sigma_{\text{g}}k_{\text{g}}.$
(6)

The specific volumes of the l- and g-fraction are related to the mean specific
volume $v$ and $\sigma_{\text{l}}$ and $\sigma_{\text{g}}$ as follows:

$v=\frac{\sigma_{\text{l}}k_{\text{l}}v_{\text{l}}+\sigma_{\text{g}}k_{\text{g}}v_{\text{g}}}{\bar{k}_{\text{fl}}}.$
(7)

This equation is equivalent to equation (2). Combining equations (5) and (7)
one has

$\displaystyle\frac{F(v,T,\sigma_{\text{l}})}{N}$ $\displaystyle\equiv
f(v,T,\sigma_{\text{l}})=-T\ln\left(\frac{v-b}{\Lambda^{3}}\right)-\frac{a}{v}$
$\displaystyle+\frac{v}{\sigma_{\text{l}}k_{\text{l}}v_{\text{l}}+\sigma_{\text{g}}k_{\text{g}}v_{\text{g}}}\left[g_{2}\sigma_{\text{l}}\sigma_{\text{g}}+T\left(\sigma_{\text{l}}\ln\sigma_{\text{l}}+\sigma_{\text{g}}\ln\sigma_{\text{g}}\right)\right].$
(8)

As seen from equation (8), the heterophase fluctuations affect the pressure,

$\displaystyle P$ $\displaystyle=\frac{\partial
f(v,T,\sigma_{\text{l}})}{\partial
v}=\frac{T}{v-b}-\frac{a}{v^{2}}+\frac{1}{\sigma_{\text{l}}k_{\text{l}}v_{\text{l}}+\sigma_{\text{g}}k_{\text{g}}v_{\text{g}}}\left[g_{2}\sigma_{\text{l}}\sigma_{\text{g}}+T\left(\sigma_{\text{l}}\ln\sigma_{\text{l}}+\sigma_{\text{g}}\ln\sigma_{\text{g}}\right)\right]$
$\displaystyle=P_{\text{VdW}}(v,T)+P_{\text{fl}}(T,\sigma_{\text{l}}).$ (9)

Here, $P_{\text{VdW}}(v,T_{\text{l}})$ is the van der Waals pressure and

$P_{\text{fl}}(T,\sigma_{\text{l}})=\frac{1}{\sigma_{\text{l}}k_{\text{l}}v_{\text{l}}+\sigma_{\text{g}}k_{\text{g}}v_{\text{g}}}\left[g_{2}\sigma_{\text{l}}\sigma_{\text{g}}+T\left(\sigma_{\text{l}}\ln\sigma_{\text{l}}+\sigma_{\text{g}}\ln\sigma_{\text{g}}\right)\right]$
(10)

is the fluctuonic pressure.

The equilibrium value of $\sigma_{\text{l}}$ is found by minimizing
$f\left(v,T,\sigma_{\text{l}}\right)$ with respect to this parameter. One
should find solutions of the equation

$\frac{\partial
f\left(v,T,\sigma_{\text{l}}\right)}{\partial\sigma_{\text{l}}}=0,$ (11)

that simultaneously satisfy the condition

$\frac{\partial^{2}f\left(v,T,\sigma_{\text{l}}\right)}{\partial\sigma_{\text{l}}^{2}}>0.$
(12)

To this end, the Ising-type effective Hamiltonian,
$g_{2}\sigma_{\text{l}}\sigma_{\text{g}}+T\left(\sigma_{\text{l}}\ln\sigma_{\text{l}}+\sigma_{\text{g}}\ln\sigma_{\text{g}}\right)$,
should be minimized. As known, it possesses the critical (bifurcation)
temperature $T_{\text{c}}=g_{2}/2$ at which
$\sigma_{\text{l}}(T_{\text{c}})=\sigma_{\text{g}}(T_{\text{c}})=1/2$. For
consistency with the van der Waals theory we have to set

$g_{2}=2T_{\text{c,VdW}}=\frac{16}{27}\frac{a}{b}$ (13)

and, as it follows from equations (6), (7),

$\bar{k}_{\text{fl,c}}=\frac{k_{\text{l,{c}}}+k_{\text{g,{c}}}}{2}\,,\quad
v_{\text{c}}=\left({\rm
1}-\varsigma\varsigma_{v}\right)\frac{v_{\text{l,{c}}}+v_{\text{g,{c}}}}{2}.$
(14)

Here,

$\varsigma=\frac{k_{\text{l}}-k_{\text{g}}}{k_{\text{l}}+k_{\text{g}}}\,,\quad\varsigma_{v}=\frac{v_{\text{g}}-v_{\text{l}}}{v_{\text{l}}+v_{\text{g}}}.$
(15)

Let us note that $f_{\text{fl}}(v_{\text{c}},T_{\text{c}})<0$. It means that
the homogeneous fluid is unstable. It transforms into a stable heterophase
fluid consisting of the l- and g-fluctuons. At the critical point, the
fluctuonic contribution to pressure is equal to

$P_{\text{fl,c}}=\frac{2}{k_{\text{l,c}}v_{\text{l,c}}+k_{\text{g,c}}v_{\text{g}}}\left(\frac{1}{4}g_{2}-T_{\text{c}}\ln
2\right)=-\frac{T_{\text{c}}\left(\ln
2-0.5\right)}{\left(1-\varsigma_{\text{c}}\varsigma_{v,\text{c}}\right)}\frac{2}{v_{\text{l,c}}+v_{\text{g,c}}}=-0.193\frac{T_{\text{c}}}{\bar{k}_{\text{fl,c}}v_{\text{c}}}.$
(16)

Thus, the total critical pressure is equal to

$P_{\text{c}}=\frac{a}{27b^{2}}-0.193\frac{T_{\text{c}}}{\bar{k}_{\text{fl,c}}v_{\text{c}}}\,,$
(17)

and the fluctuonic pressure impacts the critical compressibility factor

$Z_{\text{c}}=\frac{P_{\text{c}}v_{\text{c}}}{T_{\text{c}}}=Z_{\text{c,VdW}}-\frac{0.193}{\bar{k}_{\text{fl,c}}}=0.375-\frac{0.193}{\bar{k}_{\text{fl,c}}}.$
(18)

##  3 Recovering $\bar{k}_{\text{fl,c}}$ from experimental data

It is seen that the fluctuonic correction in equation (18) depends on just one
mesoscopic parameter $\bar{k}_{\text{fl,c}}$ and reduces the critical
compressibility factor. Assuming that the difference between
$Z_{\text{c,VdW}}$ and the measured value $Z_{\text{c},\exp}$, is due to the
heterophase fluctuations, one can estimate $\bar{k}_{\text{fl,c}}$ using
equation (18)

$\bar{k}_{\text{fl,c}}=\frac{0.193}{Z_{\text{c,VdW}}-Z_{\text{c},\exp}}.$ (19)

Results for some fluids with London dispersion intermolecular forces are
represented in table 1.

Table 1:  Values of the parameter $\bar{k}_{\text{fl,c}}$ for fluids with dispersion intermolecular forces. Fluid |  $T_{\text{c}}$, K | $Z_{\text{c},\exp}=\frac{P_{\text{c}}v_{\text{c}}}{T_{\text{c}}}$ | $\bar{k}_{\text{fl,c}}$
---|---|---|---
C2H4 | 282 | 0.270 | 1.83
CO2 | 304.19 | 0.275 | 1.93
Xe | 289.8 | 0.288 | 2.2
Kr | 209.48 | 0.288 | 2.2
CH4 | 190.8 | 0.290 | 2.27
N2 | 126.2 | 0.291 | 2.27
Ar | 150.8 | 0.291 | 2.27
O2 | 154.6 | 0.292 | 2.32
Ne | 44.5 | 0.298 | 2.47
H2 | 33.2 | 0.304 | 2.7
He4 | 5.20 | 0.308 | 2.88
He3 | 5.19 | 0.320 | 3.5

The obtained values of $\bar{k}_{\text{fl,c}}$ can be used for making
specification of the parameters of mesoscopic models of fluid [3]. Let us note
that the Lee-Yang lattice-gas model (LY-model) [4] is a partial case of the
lattice-fluctuon model (LFM) [3] and that the last one reproduces the results
of the LY-model at $k_{\text{l}}=1$, $k_{\text{g}}=0$ and
$\varsigma=\varsigma_{v}=1$. In the LY-model
$\bar{k}_{\text{fl,c}}=k_{\text{l}}/2=1/2$. In the LFM
$\varsigma,\varsigma_{v}<1$ and $k_{\text{l}}<2\bar{k}_{\text{fl,c}}$.
Therefore, $k_{\text{l}}<3.66$ for ethylene and $k_{\text{l}}<7$ for He3.

A growth of $\bar{k}_{\text{fl,c}}$ for the quantum liquids H2, He4 and He3
correlates with an increase of the de Boer parameter. Tabulated values of $B$
can be found, e.g., in [5]. For H2, He4 and He3, the parameter $B$ is equal to
$0.78$, $0.99$ and $1.76$, respectively. Correlation between
$\bar{k}_{\text{fl,c}}$ and $B$ can be interpreted as follows. Since the
l-fluctuon radius is nearly equal to the correlation radius of the direct
correlation function, it increases with an increase of the parameter $B$ that
characterizes the molecule delocalization. In the Kac model [6], the van der
Waals equation is obtained for the system of hard spheres with the range of
the attractive force tending to infinity while its strength becomes
proportionally weaker. The parameter $B$ increases with the growth of the
attractive force range and tends to infinity. Therefore, the compressibility
parameter expectedly approaches $Z_{\text{c,VdW}}$ and the accuracy of the van
der Waals equation improves with an increase of the parameter $B$.

##  4 Discussion and summary

The success of the law of the corresponding states urged van der Waals to look
for a general physical reason of the critical compressibility factor non-
universality. This problem is a subject of permanent interest. A number of
recent papers are devoted to this problem (see [5, 7, 8] and references
quoted). Considerations reported here show that the van der Waals conjecture
on the role of pseudo associates is reasonable. Moreover, it allows one to get
an estimate of the number of molecule within the l-fluctuon which is useful
while applying the mesoscopic models of gas-liquid transformations [3].

To conclude, reformulation of the van der Waals equation of state taking
account of the heterophase fluctuations in the vicinity of the critical point
allows one to associate the critical compressibility factor with the
mesoscopic parameters of the l- and g-fluctuons. The formulated equation
revives the idea of the pseudo-associates impact on the critical
compressibility factor.

## Acknowledgement

The author is deeply grateful to Professor J.M.H. Levelt Sengers for providing
the text of her excellent historical overview of the problem of fluid
heterogeneities near the critical point [9].

## References

  * [1] Van der Waals J.D., 1910 Nobel Lecture, In: Nobel Lectures, Physics 1901–1921, Elsevier Publishing Company, Amsterdam, 1967.
  * [2] Frenkel J., J. Chem. Phys., 1939, 7, 538, doi:10.1063/1.1750484.
  * [3] Bakai O., J. Mol. Liq., 2017, 235, 135–148, doi:10.1016/j.molliq.2016.12.002.
  * [4] Lee T.D., Yang C.N., Phys. Rev., 1952, 87, 410, doi:10.1103/PhysRev.87.410.
  * [5] Apfelbaum E.M., Vorob’ev V.S., J. Phys. Chem. B, 2009, 113, 3521–3526, doi:10.1021/jp808817p.
  * [6] Kac M., Uhlenbeck G.E., Hemmer P.C., J. Math. Phys., 1963, 4, 216, doi:10.1063/1.1703946.
  * [7] Apfelbaum E.M., Vorob’ev V.S., J. Chem. Phys., 2013, 139, 046101, doi:10.1063/1.4816526.
  * [8] Kulinskii V.L., J. Chem. Phys., 2013, 139, 184119, doi:10.1063/1.4829837.
  * [9] Levelt Sengers J.M.H., Physica A, 1979, 98, 363–402, doi:10.1016/0378-4371(79)90145-6.

Ukrainian 0

äåÿ ïñåâäîàñîöàö¿ Âàí äåð Âààëüñà êîåôöíò êðèòèíî¿ ñòèñëèâîñò O. Áàêàé

Íàöîíàëüíèé íàóêîâèé öåíòð ‘‘Õàðêâñüêèé ôçèêî-òåõííèé íñòèòóò’’,
âóë. Àêàäåìíà, 1, 61108 Õàðêâ, Óêðà¿íà