[trim=0 6 0 0,clip, width=0.95]fig/alg/overview.png

Figure : Overview. The proposed framework is composed of two key modules: an
offline StyleGAN-based multi-view generation and an image-conditioned neural
implicit network. During training, the neural implicit network learns unknown
geometry, texture, and camera poses for the objective of approximating the
multi-view pseudo images. At inference time, the learned neural implicit
function performs $3$D reconstruction for the object from a single image.

##  1 Proposed Method

We start with an overview of the proposed framework (Fig. ) and then present
the individual modules in detail. We first introduce an offline and effective
multi-view generator based on the pre-trained StyleGAN models, which produce
plausible multi-view images with full control over viewpoints. Then, we detail
the proposed image-conditioned neural implicit field learning framework,
including a neural implicit network, differentiable rendering procedure, and
an uncertainty prediction module. These three modules work jointly for the
objective of exploiting the pseudo images to learn generic object shape priors
and perform $3$D reconstruction from a single input image.

###  StyleGAN based Multi-view Generation

We briefly review the embedding space of the StyleGAN [karras2019style,
karras2020analyzing]. Typically, a generator $G(\cdot)$ samples a latent code
$\mathbf{z}$ from a pre-defined distribution $\mathcal{Z}$ such as the normal
distribution, and produces an output image $\mathbf{I}$. The code $\mathbf{z}$
is first mapped to an intermediate latent space $\mathcal{W}$ via a Multi-
layer Perceptron (MLP), and then $\mathcal{W}$ is transformed to
$\mathcal{W}^{+}$ space by $16$ learned affine transformations. The generator
$G(\cdot)$ projects $\mathbf{W}$ to the final image:
$G(\mathbf{W})\\!=\\!\mathbf{I}$. Such latent codes have been shown to learn
various disentangled semantics [shen2021closed, harkonen2020ganspace]. For
instance, StyleGANRender [zhang2020image] finds that the latent codes
$\mathbf{W}_{v}\\!:=\\!(\mathbf{w}_{1},\mathbf{w}_{2},\mathbf{w}_{3},\mathbf{w}_{4})\\!\in\\!\mathbb{R}^{4\times
512}$ in the first $4$ layers control camera viewpoints. That is, given a
source and reference generated image pair $(\mathbf{I}^{S},\mathbf{I}^{R})$
with their latent codes $(\mathbf{W}^{S},\mathbf{W}^{R})$, we can generate an
image of the source object with the reference viewpoint by swapping
$(\mathbf{W}_{v}^{S},\mathbf{W}_{v}^{R})$ and keeping the rest dimensions of
$\mathbf{W}^{S}$. We denote this multi-view generation strategy as _Baseline_.
However, while $\mathbf{W}_{v}$ indeed alters the object viewpoint, it still
perceives the shape of the reference object, as shown in Fig. 1. It is
difficult to develop a multi-view stereo system from noisy multi-view images
with inconsistent shapes. To tackle this issue, inspired by SeFa
[shen2021closed], we propose a novel offline multi-view generation mechanism
that generates plausible images with enhanced cross-view consistency. SeFa
suggests that the weight parameters in early affine transformations contain
essential knowledge of image variations. One can obtain interpretable
directions in the latent space by computing eigenvectors of their weight
matrices and selecting eigenvectors with the largest eigenvalues. With this
observation, we propose to filter the viewpoint-irrelevant features in
$\mathbf{W}_{v}^{R}$ guided by the eigenvectors of the $k$ largest
eigenvalues, which are computed from the weight parameters in the first $4$
transformations. It can be formulated as:

$\argmin_{\alpha}||\hat{\mathbf{W}}_{v}-\mathbf{W}_{v}^{R}||^{2},\quad\hat{\mathbf{W}}_{v}=\alpha\mathbf{V}+\mathbf{W}_{v}^{S},$
(1)

where $\hat{\mathbf{W}}_{v}$ is the enhanced viewpoint latent code.
$\mathbf{\alpha}$ is a $k$-dim viewpoint coefficient.
$\mathbf{V}\in\mathbb{R}^{k\times 512}$ denotes the eigenvectors of
$\mathbf{A}^{T}\mathbf{A}$ associated with the $k$ largest eigenvalues.
$\mathbf{A}\in\mathbb{R}^{m\times 512}$ are the weights of transformations.
$\alpha\mathbf{V}$ is duplicated into four rows in Eqn. 1 and $\alpha$ can be
solved by gradient descent. Finally, we generate a new image by combining
$\hat{\mathbf{W}}_{v}$ and remaining the dimensions of $\mathbf{W}^{S}$.

[b]0.48 [width=]fig/alg/multiview_examples.jpg [b]0.48
[width=]fig/alg/implict_function.png

Figure :  Figure :  Figure :  (a) Multi-view generator comparisons of
_Baseline_ [zhang2020image] and our approach (_Ours_). As can be observed,
besides the viewpoint, the baseline perceives shape cues from the references
(yellowyellow circle, best view in zoom in). While our generated images show
more consistency in object shape across views, which benefits shape learning.
Please refer to Supp for more categories. (b) The architecture of our neural
implicit fields and differentiable renderer.

Training Pseudo Images Given a pre-trained StyleGAN model, we first synthesize
training images and filter out images that have more than one instance or an
unrealistic instance, resulting in $N$ training samples
$\\{\mathbf{I}^{j}\\}_{j=1}^{N}$. Then, we manually select $n$ reference view
samples, which roughly cover the common object viewpoints ranging from
$0-360^{\circ}$ in azimuth. Finally, for each training sample
$\mathbf{I}^{j}$, we produce $n$ multi-view images of the same object with
fixed camera poses as the pseudo images
$\\{\mathbf{I}^{*j}_{i}\\}_{i=1}^{n}\in\mathbb{R}^{W\times H\times 3}$. We
further apply the work [qin2020u2] to obtain corresponding instance
segmentation $\\{\mathbf{M}^{j}_{i}\\}_{i=1}^{n}$ of pseudo images. Pose
Sampling Space We assume a pinhole camera model
$\mathcal{C}=(\mathbf{t},\mathbf{K})$, where $\mathbf{K}\in\mathbb{R}^{3\times
3}$ is the intrinsic parameter. We represent the camera pose/extrinsic
parameters $\mathbf{t}=(\mathbf{c}_{0},\mathbf{q})$ based on its $3$D position
$\mathbf{c}_{0}\in\mathbb{R}^{3}$ and its rotation from a canonical view.
$\mathbf{q}\in\mathbb{R}^{4}$ is the quaternion vector representing the camera
rotation. We assume the observed object is approximately inside the unit
sphere. We further assume a known intrinsic and the principal point at the
image center, as commonly assumed in stereo systems [kar2017learning]. For the
pose initialization, we manually annotate the camera poses
$\\{\mathbf{\tilde{t}}_{i}\\}_{i=1}^{n}$ for the $n$ reference images, which
initializes the pseudo images’ poses,
$\mathbf{t}_{i}^{j}=\mathbf{\tilde{t}}_{i},j\\!\in\\![1,N]$. _Since we only
need to annotate $n$ reference samples per object category, it is far more
practical than prior works that require camera pose label per training sample_
(see Tab. LABEL:tab:3D_modeling_review).

###  Image-conditioned Neural Implicit Field

$3$D Geometry Representation As illustrated in Fig. 1, the object geometry to
be reconstructed is represented by the function [yariv2020multiview]:
$\mathcal{F}:\gamma(\mathbf{x})\rightarrow(s,\mathbf{f})$ that maps a point
$\mathbf{x}\in\mathbb{R}^{3}$ to its signed distance value $s$ to the object
surface and a local geometry feature $\mathbf{f}\in\mathbb{R}^{d_{f}}$.
$\gamma(\cdot)$ denotes a positional encoding operator on $\mathbf{x}$ with
$6$ exponentially increasing frequencies introduced in NeRF
[mildenhall2020nerf]. The surface $\mathcal{S}$ is represented as the zero
level set of MLP $\mathcal{F}$ with learnable parameters $\theta$:

$\mathcal{S}=\\{\mathbf{x}\in\mathbb{R}^{3}|\mathcal{F}^{(s)}_{\theta}(\gamma(\mathbf{x}))=0\\}.$
(2)

. Neural Renderer Given a pixel $p$ of a masked input image, we march a ray
$\mathbf{r}=\\{\mathbf{c}_{0}+t\mathbf{v}|t\geqslant 0\\}$, where
$\mathbf{c}_{0}$ is the camera position and $\mathbf{v}$ the viewing
direction. $\mathbf{\hat{x}}$ denotes the first intersection between the ray
$\mathbf{r}$ and the surface, which can be efficiently detected via the sphere
tracing [hart1996sphere] and implemented in a differentiable manner. As in
Fig. 1, the rendered color of the pixel $p$ is encoded as [yariv2020multiview]

$\mathcal{G}:(\mathbf{\hat{x}},\mathbf{\hat{n}},\mathbf{f},\mathbf{v})\rightarrow\mathbf{c},$
(3)

which is a function of the surface properties at $\mathbf{\hat{x}}_{p}$,
including the surface point $\mathbf{\hat{x}}_{p}$, the surface normal
$\mathbf{\hat{n}}_{p}\in\mathbb{R}^{3}$, the local geometry feature
$\mathbf{f}_{p}$, and a viewing direction $\mathbf{v}_{p}\in\mathbb{R}^{3}$.
Similarly, the function $\mathcal{G}$ is implemented as an MLP with learnable
parameters $\phi$. The surface normal $\mathbf{\hat{n}}_{p}$ can be computed
by the spatial derivative
$\frac{\delta\mathcal{F}^{(s)}}{\delta\mathbf{x}_{p}}$ via back-propagation
through the network $\mathcal{F}$. Incorporating the surface normal and view
direction enables $\mathcal{G}$ to represent the light reflected from a
surface point $\mathbf{x}$ across different viewpoints [mildenhall2020nerf,
yariv2020multiview]. Also, according to [yariv2020multiview], introducing the
local geometry feature vector $\mathbf{f}$ allows the renderer to handle more
complex appearances, which might appear in GAN-generated images. The neural
implicit field can be learned through back-propagation without any $3$D
supervision by comparing the rendered images with the multi-view pseudo
images. Minimizing this error encourages multi-view photo consistency because
only when the point is on the actual surface, the MLPs
$(\mathcal{F},\mathcal{G})$ will predict accurate color with fewer multi-view
variations. Image Encoder We propose to utilize architecture to condition on a
single image, such that the learned neural implicit field can generalize to a
new object instance without re-training. Specifically, as shown in Fig. , we
use an image encoder $\mathcal{E}$ as a hyper network [ha2016hypernetworks,
lin2020sdf] to predict $\theta$ (the parameters of $\mathcal{F}$), written as
$\theta=\mathcal{E}_{\Phi}(\mathbf{I})$, where $\Phi$ is the neural network
weights.

###  Model Learning

Given $N$ sets of masked pseudo images
$\\{\mathbf{I}^{j},\\{\mathbf{I}^{*j}_{i}\\}_{i=1}^{n}\\}_{j=1}^{N}$, along
with their initial camera poses
$\\{\\{\mathbf{t}^{j}_{i}\\}_{i=1}^{n}\\}_{j=1}^{N}$, we optimize the encoder
parameters $\Phi$, texture MLP parameters $\phi$ and camera poses
$\\{\\{\mathbf{t}^{j}_{i}\\}_{i=1}^{n}\\}_{j=1}^{N}$ by minimizing the loss:

$\mathcal{L}=\sum_{j=1}^{N}(\mathcal{L}_{RGB}+\lambda_{mask}\mathcal{L}_{Mask}+\lambda_{eik}\mathcal{L}_{Eik}),$
(4)

where $\mathcal{L}_{RGB}$ is photometric loss, $\mathcal{L}_{Mask}$ is
silhouette loss, $\mathcal{L}_{Eik}$ is Eikonal regularization, and
$\lambda_{*}$ are loss weights. Photometric Loss Let $\mathbf{I}^{*}_{p}$,
$\mathbf{M}_{p}\\!\in\\!\\{0,1\\}$ be the RGB and silhouette values of pixel
$p$ in an image sample $\mathbf{I}^{*}_{i}$ taken at the view direction
$\mathbf{v}_{p}$ associated with camera $\mathcal{C}_{i}$. $p\in P$ indexes
all pixels in the input image set $\\{\mathbf{I}^{*}_{i}\\}_{i=1}^{n}$. The
photometric loss is defined on mini-batches of pixels in $P$:

$\mathcal{L}_{RGB}=\frac{1}{|P|}\sum_{p\in
P^{in}}|\mathbf{I}^{*}_{p}-\mathcal{G}(\mathbf{\hat{x}}_{p},\mathbf{\hat{n}}_{p},\mathbf{f}_{p},\mathbf{v}_{p})|,$
(5)

where $\mathbf{f}_{p}$, $\mathbf{\hat{x}}_{p}$, $\mathbf{\hat{n}}_{p}$ are
defined in Eqn. 3. $P^{in}\subset P$ represents the subset of pixels $P$ where
intersection has been found and $\mathbf{M}_{p}=1$. $|\cdot|$ denotes the
$L_{1}$ loss. Silhouette Loss We define the silhouette loss as

$\mathcal{L}_{Mask}=\frac{1}{|P|}\sum_{p\in
P^{out}}CE(\mathbf{M}_{p},\mathbf{\hat{M}}_{p}),$ (6)

where $\mathbf{\hat{M}}$ is the masked rendering. $P^{out}=P-P^{in}$
represents the indices in the mini-batch for which there is no ray-geometry
intersection or $\mathbf{M}_{p}=0$. $CE(\cdot,\cdot)$ denotes the cross-
entropy loss. Conventionally, given the ray
$\mathbf{r}_{p}=\\{\mathbf{c}_{0}+t\mathbf{v}_{p}|t\geqslant 0\\}$ of pixel
$p$, $\mathbf{\hat{M}}_{p}$ is defined as:

$\mathbf{\hat{M}}_{p}=\left\\{\matrix{1}&\mathbf{r}_{p}\cap\mathcal{S}\\\
0\textup{otherwise}.\right.$ (7)

To make this differentiable, we follow [yariv2020multiview] and compute

$\mathbf{\hat{M}}_{p}=\textup{sigmoid}(-\beta\,\underset{t\geq
0}{\min}\mathcal{F}^{(s)}(\mathbf{c}_{0}+t\mathbf{v}_{p})).$ (8)

When $\beta\rightarrow\\!\infty$, $\mathcal{F}^{(s)}\\!<\\!0$ means inside the
surface and $\mathcal{F}^{(s)}\\!>\\!0$ outside. Eikonal Regularization A
special property of signed distance functions is their differentiability with
a gradient of unit norm, satisfying the Eikonal equation
$||\nabla\mathcal{F}||_{2}=1$ [osher2004level, gropp2020implicit]. We thus
encourage our implicit geometry representation to satisfy the Eikonal
property:

$\mathcal{L}_{Eik}=\sum_{\mathbf{\tilde{x}}}\Vmatrix||\nabla_{\mathbf{\tilde{x}}}\mathcal{F}^{(s)}(\mathbf{\tilde{x}})||_{2}-1_{2}^{2},$
(9)

where $\mathbf{\tilde{x}}$ is uniformly sampled at the $3$D region of
interest.

###  Uncertainty Prediction Module

The pseudo images are inherently aleatoric, i.e., there might be areas with
either notable artifacts in one image or with inconsistent shape/texture
across images. To adaptively treat these problematic areas in learning, we
propose to use Bayesian learning [kendall2017uncertainties, xu2021digging] to
model this aleatoric uncertainty. Specifically, we exploit the feature space
of the image encoder $\mathcal{E}$ and train a shallow decoder to estimate an
uncertainty map $\mathbf{U}$, which has the same size as pseudo images.
Formally, we model the observed color $\mathbf{I}^{*}_{p}$ at pixel $p$ with a
likelihood function $\mathcal{P}(\mathbf{I}^{*}_{p})$ that follows the
Laplacian distribution with ray-dependent variance $\sigma_{p}$:

$\mathcal{P}(\mathbf{I}^{*}_{p})=\frac{1}{2\sigma_{p}}\textup{exp}\left(-\frac{|\mathbf{I}^{*}_{p}-\mathcal{G}(\mathbf{\hat{x}}_{p},\mathbf{\hat{n}}_{p},\mathbf{f}_{p},\mathbf{v}_{p})|}{\sigma_{p}}\right),$
(10)

where $\sigma_{p}$ denotes the uncertainty. Since $L_{1}$ distance is less
sensitive to outliers, which is more suitable for optimizing the rendered
appearance against the pseudo RGB values. Thus, we adopt Laplacian likelihood
to model the inconsistent uncertainty distribution. To find the parameters
best explaining the model, we maximize the likelihood function, i.e.,
minimizing the negative log-likelihood:

$-\\!\textup{log}(\mathcal{P}(\mathbf{I}^{*}_{p}))\\!=\\!\frac{|\mathbf{I}^{*}_{p}\\!-\\!\mathcal{G}(\mathbf{\hat{x}}_{p},\mathbf{\hat{n}}_{p},\mathbf{f}_{p},\mathbf{v}_{p})|}{\sigma_{p}}\\!+\textup{log}\sigma_{p}\\!+\textup{log}2.$
(11)

Therefore, we update the photometric loss as:

$\mathcal{L}_{RGB}\\!=\frac{1}{|P|}\\!\sum_{p\in
P_{in}}\\!\left(e^{-\mathbf{U}_{p}}|\mathbf{I}^{*}_{p}\\!-\mathcal{G}(\mathbf{\hat{x}}_{p},\mathbf{\hat{n}}_{p},\mathbf{f}_{p},\mathbf{v}_{p})|+\mathbf{U}_{p}\right).$
(12)

We train this loss on mini-batches of pixels in $P$. Here $P$ indexes all
pixels in the input multi-view image set
$\\{\mathbf{I}_{i=1}^{*}\\}_{i=1}^{n}$, which contribute to the same object’s
depth, texture and uncertainty learning. During training, for the surface
point $\hat{\mathbf{x}_{p}}$, the first term $e^{-\mathbf{U}_{p}}$ can be seen
as a weighted distance which assigns larger weights to less uncertain pixels.
The second term $\mathbf{U}_{p}$ is a penalty term. $\mathcal{L}_{RGB}$
encourages the neural implicit representation to bring together the texture
information of all cross-view corresponding pixels in all images.
Consequently, the pixel-wise uncertainty value is able to mine the multi-view
inconsistencies among those pixels which contribute to the same $3$D surface
point. In practice, we train a $2$-layer convolutional network to predict the
log variance $\mathbf{U}_{p}:=\textup{log}\sigma_{p}$ (please refer to Supp
for more details).

###  Implementation Details

The encoder $\mathcal{E}$ is implemented as a ResNet-$18$ [he2016deep]
followed by fully-connected layers. Both $\mathcal{F}$ and $\mathcal{G}$ have
$4$ fully-connected layers. For the main experiment, we set $N{=}2{,}000$,
$n{=}40$, $k{=}5$, $W{=}H{=}256$, $d_{f}{=}256$, $\beta=50$,
$\lambda_{mask}{=}0.01$, $\lambda_{eik}{=}0.1$. We implement our model in
Pytorch and use the Adam optimizer with a learning rate of $1e{-}4$ for both
network and camera pose parameters.


# The Symmetry-Preserving Mean Field Condition for Electrostatic Correlations
in Bulk

Zhonghan Hu<EMAIL_ADDRESS>Qingdao Institute for Theoretical and
Computational Sciences (QiTCS), Shandong University, Qingdao, 266237, P. R.
China  Institute of Frontier and Interdisciplinary Science, Shandong
University, Qingdao, 266237, P. R. China

###### Abstract

Accurate simulations of a condensed system of ions or polar molecules are
concerned with proper handlings of the involved electrostatics. For such a
Coulomb system at a charged planar interface, the Coulomb interaction averaged
over the lateral directions with preserved symmetry serves as a necessary
constraint in building any accurate handling that reconciles a simulated
singlet charge density with the corresponding macroscopic charge/dielectric
response. At present, this symmetry-preserving mean-field (SPMF) condition
represented in the reciprocal space, is conjectured to be necessary for a
simulated bulk system to reproduce correctly the charge structure factor of
the macroscopic bulk, as well. In this work, we further examine analytically
the asymptotic behavior of the charge structure factor at small wavenumbers
for an arbitrary charge-charge interaction. In light of our theoretical
predictions, simulations with a length of nearly $0.1$ micron are carried out
to demonstrate that, typical efficient methods violating the SPMF condition
indeed fail to capture the exact charge correlations at small wavenumbers for
both ionic and polar systems. However, for both types of systems, these
existing methods can be simply amended to match the SPMF condition and
subsequently to probe precisely the electrostatic correlations at all length
scales.

There have been much advances regarding efficient and accurate treatments of
the long-ranged electrostatics for molecular dynamics and Monte-Carlo
simulations in recent yearsAllen and Tildesley (2017); Yi  _et al._ (2015);
Lowe  _et al._ (2018); Maggs (2002); Fukuda (2013); Wang  _et al._ (2016a, b);
Liang  _et al._ (2015); Girotto  _et al._ (2017); Bakhshandeh  _et al._
(2018); Zhang  _et al._ (2018); Wang  _et al._ (2018); Urano  _et al._ (2020);
Yuan  _et al._ (2021); Jin  _et al._ (2021). These techniques differ from each
other in various means: for bulkMaggs (2002); Fukuda (2013); Wang  _et al._
(2016a, b); Zhang  _et al._ (2018); Wang  _et al._ (2018); Urano  _et al._
(2020); Jin  _et al._ (2021) versus for interfacesLiang  _et al._ (2015);
Girotto  _et al._ (2017); Bakhshandeh  _et al._ (2018); Yuan  _et al._ (2021),
efficient algorithms based on existing formulationsWang  _et al._ (2016b);
Urano  _et al._ (2020); Yuan  _et al._ (2021) versus newly designed functional
forms of the simulated electrostaticsFukuda (2013); Liang  _et al._ (2015);
Girotto  _et al._ (2017); Bakhshandeh  _et al._ (2018), for classical non-
polarizable systemsMaggs (2002); Fukuda (2013); Wang  _et al._ (2016a, b);
Urano  _et al._ (2020); Jin  _et al._ (2021) versus for systems with
dielectric contrastLiang  _et al._ (2015); Girotto  _et al._ (2017);
Bakhshandeh  _et al._ (2018); Yuan  _et al._ (2021), and methods well
characterized by a simulation Hamiltonian in the framework of classical
statistical mechanicsFukuda (2013); Liang  _et al._ (2015); Wang  _et al._
(2016a, b); Girotto  _et al._ (2017); Bakhshandeh  _et al._ (2018); Urano  _et
al._ (2020); Yuan  _et al._ (2021) versus others that might not be easily
characterized in the same wayMaggs (2002); Zhang  _et al._ (2018); Wang  _et
al._ (2018); Jin  _et al._ (2021). The development and implementation of these
techniques have not only enhanced greatly our knowledge of complex Coulomb
systems but also persistently called for a conceptual understanding of the
intrinsic connections among the many techniques themselves from a transparent
theoretical viewpoint.

The recent symmetry-preserving mean-field theoryHu (2014a); Yi  _et al._
(2017a); Pan  _et al._ (2017, 2019) might shed some light on the connections.
The core of this theory applied to the interfacial electrostatics is the so-
called SPMF condition — the lateral average of the pairwise charge-charge
interaction, $\nu({\mathbf{r}})$ involved in an accurate simulation must equal
that of the Coulomb interactionPan  _et al._ (2019)

$\left<\nu({\mathbf{r}})\right>_{\rm sp}\equiv\frac{1}{A}\iint
dxdy\,\nu({\mathbf{r}})=\frac{-\lvert z\rvert}{2A\varepsilon_{0}},$ (1)

where $A=L_{x}L_{y}$ is the cross-sectional area of the simulation box and
$\varepsilon_{0}$ is the vacuum permittivity. $\left<\,\right>_{\rm sp}$
stands for the symmetry-preserving (sp) average over the lateral directionsHu
(2014a). Whenever it operates on a two-dimensional (2D) Fourier series,
integrations of trigonometric functions all vanish and only the first Fourier
coefficient is left overPan  _et al._ (2019). $\nu({\mathbf{r}})$ can always
be defined by rewriting the electrostatic energy in the simulation asHu
(2014b); Yi  _et al._ (2017b); Yuan  _et al._ (2021)

${\cal U}_{\rm elec}=\sum_{i<j}q_{i}q_{j}\nu({\mathbf{r}}_{ij}),$ (2)

where $q_{i}$ is the charge of the $i$-th particle and
${\mathbf{r}}_{ij}={\mathbf{r}}_{i}-{\mathbf{r}}_{j}$ denotes the relative
vector between the $i$\- and $j$-th particles.

Figure 1: The electric potentials (top) and the corresponding electric fields
(bottom) produced by a pair of unit charges (red solid lines for $q=+1$ and
blue solid lines for $q=-1$) separated at $d$ interacting through the
laterally averaged Coulomb interaction (left), the quadratic infinite boundary
term (middle) and the periodic series (right) of Eq. (4), respectively. Any
potential (black) produced by the pair is the sum of the two individual
potentials (colored). The dash lines indicate the periodic images of the pair.

The SPMF condition is necessary for a finite simulation consist of a few
thousands to billions of particles to reproduce the charge/dielectric response
of the corresponding macroscopic Coulomb system ($\sim 10^{23}$ particles)Pan
_et al._ (2019). To illustrate this statement precisely, let us consider an
example of mobile ions confined between two macroscopic plate electrodes with
fixed surface charges, $\pm\sigma$. In this model system, the mobile ions
driven by the electric field must form a nontrival charge distribution along
the normal ($z$) direction in the interfacial region around each electrode. At
thermal equilibrium under normal conditions, the fixed surface charge of each
electrode ($\pm\sigma$) should be screened completely by the compensating
mobile surface charge ($\mp\sigma$) resulting from the integration of the
equilibrium singlet charge density over the corresponding interfacial region,
such that, mobile ions in the bulk region feel no net electrostatic force on
average. In a numerical simulation employing $\nu({\mathbf{r}})$ for the
charge-charge interaction under the usual 2D periodic boundary condition
(PBC), the SPMF condition becomes necessary and sufficient for producing the
exact compensating surface charge, $\mp\sigma$. On the other hand, the
condition is necessary but insufficient for producing the detailed variation
of the equilibrium singlet charge density along the $z$ direction. However, as
long as $\nu({\mathbf{r}})$ differs from $\nu^{\rm e2d}({\mathbf{r}})$
involved in the formally exact Ewald2D sum methodArnold  _et al._ (2002);
Mináry  _et al._ (2002); Hu (2014b); Pan and Hu (2014) by a slowly varying
long-ranged component, the SPMF theory has suggested that interfacial
structure, thermodynamics as well as dynamics in the $z$ direction can often
be obtained accurately from the approximated but much more efficient methodHu
(2014a). In terms of $\nu({\mathbf{r}})$, all existing methods for interfacial
electrostatics can be classified into three categories: the formally exact
Ewald2D method and its efficient algorithms, approximated SPMF methods (e.g.
Bengtsson (1999); Zhang  _et al._ (2020) and more in ref. Pan  _et al._
(2019)), and other methods which violate the SPMF condition but could be
amended to work well for interfacesPan  _et al._ (2019).

While the guiding constraint, Eq. (1), and the associated SPMF theory provide
a transparent analysis of the interfacial electrostatics, the question arises
as to whether or not the electrostatics in the bulk phase could be simply
analyzed as well. Since the equilibrium singlet charge density in bulk
vanishes absolutely, the existence of any possible constraint must be related
to the charge correlations in general. In practice, various methods with
diverse $\nu({\mathbf{r}})$ have all proved successful for treating the local
correlations of Coulomb systemsAllen and Tildesley (2017); D. Wolf and
Eggebrecht (1999); Wang  _et al._ (2016a); Chen  _et al._ (2004) and at first
glance, it may appear very difficult to built any direct relation between
$\nu({\mathbf{r}})$ and the charge correlations. In this work, as opposed to
validating the usual correlations in the real space, we focus on the charge
structure factor defined as the ensemble average in the reciprocal spaceHansen
and McDonald (2006),

$S({\mathbf{k}})=\frac{1}{V}\sum_{ij}q_{i}q_{j}\left<e^{i{\mathbf{k}}\cdot{\mathbf{r}}_{ij}}\right>,$
(3)

where $V=L_{x}L_{y}L_{z}$ is the volume of the cubic simulation cell.
$S({\mathbf{k}})$ at ${\mathbf{k}}=2\pi(k_{x}/L_{x},k_{y}/L_{y},k_{z}/L_{z})$
with $k_{x}$, $k_{y}$ and $k_{z}$ all integers, probes the charge correlation
at the discrete wave vector that is compatible with the reciprocal lattice of
the finite simulation cell. Correspondingly, the averaged Coulomb interaction
represented in the reciprocal space is known to be a combination of a non-
periodic quadratic term and a periodic Fourier seriesPan  _et al._ (2017,
2019),

$\frac{-\lvert
z\rvert}{2\varepsilon_{0}A}=\frac{-z^{2}}{2\varepsilon_{0}V}+\frac{1}{V}\sum_{k_{z}\neq
0}\frac{e^{ikz}}{\varepsilon_{0}k^{2}},$ (4)

where $k=\lvert{\mathbf{k}}\rvert$ with ${\mathbf{k}}$ specified to be
$(0,0,2\pi k_{z}/L_{z})$ normal to the $xy$ plane. Eq. (4) is valid up to a
constant for any $z\in[-L_{z},L_{z}]$. The quadratic term, previously called
the infinite boundary term, identifies with the excluded $k\to 0$ term of the
periodic series, provided that the divergence caused by $1/k$ and $1/k^{2}$
are all removedHu (2014b).

The electric field produced by this quadratic term for a pair of unit charges
($q=\pm 1$) is a constant proportional to the distance between the two charges
in the prescribed direction (see Fig. 1). Each quadratic term in an
instantaneous configuration thus adds up to contribute a constant electric
field proportional to the non-periodic (itinerant) total dipole moment. Being
the reality of a non-periodic quantity, the itinerant dipole moment possesses
no translational invariance and becomes ill-defined upon the PBC transform of
any non-zero charge. As such, the quadratic infinite boundary term, as its
name says, must be merely responsible for the effect of the macroscopic
boundary, which is deemed irrelevant to any well-defined bulk property such as
$S({\mathbf{k}})$ of Eq. (3).

When focusing on the translationally invariant electrostatic correlations
under the full PBC, one might exclude the non-periodic term and subsequently
conjecture that the pairwise interaction accounting for an accurate
$S({\mathbf{k}})$ with ${\mathbf{k}}=(0,0,2\pi k_{z}/L_{z})$ satisfies
necessarily,

$\left<\nu({\mathbf{r}})\right>_{\rm sp}=\frac{1}{V}\sum_{k_{z}\neq
0}\frac{e^{ikz}}{\varepsilon_{0}k^{2}},$ (5)

which defines the SPMF condition for the electrostatics in the bulk phase.
Similar constraints apply to symmetry-preserving averages in other directions
for the purpose of determining precisely $S({\mathbf{k}})$ at the
corresponding wave vectors. Noting that $1/(\varepsilon_{0}k^{2})$ is the 3D
Fourier transform of the Coulomb interaction and $\left<\,\right>_{\rm sp}$
removes all the 3D Fourier coefficients associated with
$e^{i{\mathbf{k}}\cdot{\mathbf{r}}}$ subject to $k_{x}\neq 0$ or $k_{y}\neq
0$, this SPMF condition suggests essentially that, any chosen lateral average
of $\nu({\mathbf{r}})$ represented in the reciprocal space must equal that of
the Coulomb interaction such that the finite simulation under the 3D PBC is
able to reproduce correctly the macroscopic charge correlations in the
corresponding normal direction of the reciprocal space.

It turns out to be informative to examine the existing methods for bulk
electrostatics in terms of Eq. (5). The pairwise interaction of the well known
Ewald3D sum method with the tinfoil boundary condition (e3dtf) readsHu
(2014b); Yi  _et al._ (2017b); Pan  _et al._ (2019)

$\nu^{\rm e3dtf}({\mathbf{r}})=\tau^{\rm
3D}+\frac{1}{V}\sum_{{\mathbf{k}}\neq{\mathbf{0}}}\frac{e^{i{\mathbf{k}}\cdot{\mathbf{r}}}}{\varepsilon_{0}k^{2}},$
(6)

which indeed satisfies Eq. (5) (up to the constant $\tau^{\rm 3D}$) in any
direction. The e3dtf method is therefore considered to be formally exact for
probing the electrostatic correlations in the reciprocal space. This shift of
perspective resolves the long-time concern regarding the artificial anisotropy
introduced by the e3dtf methodCaillol (1992); Allen and Tildesley (2017); Yi
_et al._ (2015): charge correlations in the reciprocal space under the 3D PBC
are as isotropic as they are on the 3D hypersurface of a 4D sphereCaillol
(1992). In another word, any anisotropy introduced by $\nu^{\rm
e3dtf}({\mathbf{r}})$ in the real space won’t be problematic any more once the
quantity of interest is computed properly by the inverse Fourier transform of
the corresponding ${\mathbf{k}}$-dependent ensemble average. When the Ewald3D
sum is associated with other boundary conditions that depend on the itinerant
dipole momentde Leeuw  _et al._ (1980); Caillol (1994); Hu (2014a), Caillol
has argued that the itinerant dipole moment is a decoupled collective variable
that behaves as an independent harmonic oscillator for an electrolyteCaillol
(1994), which well supports our exclusion of the redundant non-periodic
infinite boundary term.

Distinct pairwise interactions involved in other useful methods in principle
violate the SPMF condition. As an example, the recent zero-multipole (zm)
method introduces deformed Coulomb interactions, $\nu({\mathbf{r}}|l,\alpha)$
parameterized by the order of the multipole moment $l$ and the damping factor
$\alpha$Fukuda (2013); Wang  _et al._ (2016a). Its 3D Fourier transform, which
is the Fourier coefficient of the corresponding pairwise interaction, differs
(d) from $1/(\varepsilon_{0}k^{2})$ by

$\hat{\nu}_{\rm
d}({\mathbf{k}}|l,\alpha)=1/(\varepsilon_{0}k^{2})-\hat{\nu}({\mathbf{k}}|l,\alpha).$
(7)

For three typical sets of parameters: $l=0$ (zm0), $l=2$ and $\alpha=0$ (zm2),
and $l=3$ and $\alpha=0$ (zm3)Fukuda (2013); Wang  _et al._ (2016a), the
remaining Fourier coefficients read explicitly

$\hat{\nu}_{\rm
d}({\mathbf{k}}|0,\alpha)=e^{-k^{2}/(4\alpha^{2})}/(\varepsilon_{0}k^{2}),$
(8) $\hat{\nu}_{\rm d}({\mathbf{k}}|2,0)=\frac{15\left[(3-k_{c}^{2})\sin
k_{c}-3k_{c}\cos k_{c}\right]}{\varepsilon_{0}k^{2}k_{c}^{5}},$ (9)

where $k_{c}=kr_{c}$ with $r_{c}$ the cutoff distance, and

$\hat{\nu}_{\rm d}({\mathbf{k}}|3,0)=\frac{105\left[(15-6k_{c}^{2})\sin
k_{c}-(15-k_{c}^{2})k_{c}\cos k_{c}\right]}{\varepsilon_{0}k^{2}k_{c}^{7}},$
(10)

respectively. $\hat{\nu}_{\rm d}({\mathbf{k}}|l,\alpha)$ is nonzero unless
both $\alpha=0$ and $r_{c}\to\infty$. To match the SPMF condition in a
prescribed (e.g. $z$) direction for the general case of finite $\alpha$ and
$r_{c}$, the amended zero-multipole (azm) method must include an extra
electrostatic term

${\cal U}_{\rm extra}=\sum_{i<j}q_{i}q_{j}\nu_{\rm sp}(z_{ij}|l,\alpha),$ (11)

where $\nu_{\rm sp}(z|l,\alpha)$ sums over the remaining Fourier coefficients
in the direction with preserved symmetry as in Eq. (5)

$\nu_{\rm sp}(z|l,\alpha)=\frac{1}{V}\sum_{k_{z}\neq 0}\hat{\nu}_{\rm
d}({\mathbf{k}}|l,\alpha)e^{ikz}.$ (12)

Besides, both electrostatic terms corresponding to $\nu_{\rm sp}(x|l,\alpha)$
and $\nu_{\rm sp}(y|l,\alpha)$ can be added to retain isotropy and
subsequently match the SPMF condition in all three directions.

It is possible to examine analytically the influence of the different pairwise
interactions on the charge structure factor, $S({\mathbf{k}})$ for both ionic
and polar systems. Because various methods have been successful in probing the
short-ranged correlationsFukuda (2013); Wang  _et al._ (2016a); Zhang  _et
al._ (2018); Wang  _et al._ (2018); Jin  _et al._ (2021), which is related to
$S({\mathbf{k}})$ at relatively large wavenumbers, deficiencies of the
pairwise interaction violating the SPMF condition must be most evident when
one computes $S({\mathbf{k}})$ at small wavenumbers that reflect otherwise the
long-ranged correlations. For a conducting ionic fluid interacting through the
Coulomb force, the seminal work of Stillinger and Lovett argued that $S(k)$
has the universal form at small wavenumbersStillinger and Lovett (1968)

$S({\mathbf{k}})=0+k_{b}T\varepsilon_{0}k^{2}+{\cal O}(k^{4}),$ (13)

where $k_{b}$ is the Boltzmann constant and $T$ is the temperature. This exact
asymptotic form of course agrees with the well known Debye-Huckel limitHansen
and McDonald (2006)

$S({\mathbf{k}})=\frac{k_{b}T\varepsilon_{0}}{\lambda^{2}+\hat{\nu}({\mathbf{k}})\varepsilon_{0}},$
(14)

where $\hat{\nu}({\mathbf{k}})=1/(\varepsilon_{0}k^{2})$ and $\lambda$ is the
Debye length defined through $\lambda^{2}=k_{b}T\varepsilon_{0}/I$ with the
ionic strength $I$ given by the sum of squares of charges divided by the
volume. Both forms in Eqs. (13) and (14) characterize the charge correlations
in the conducting fluid — any fixed charge is completely screened by the
surrounding mobile charges as a consequence of the divergence of
$\hat{\nu}({\mathbf{k}})$ at ${\mathbf{k}}\to 0$Hansen and McDonald (2006).
This complete screening effect is irrespective of any details of the short-
ranged non-electrostatic interaction.

In a simulation employing otherwise the deformed Coulomb interaction, the
asymptotic behavior of $S({\mathbf{k}})$ can be readily obtained by replacing
$\hat{\nu}({\mathbf{k}})$ in Eq. (14) with $\hat{\nu}({\mathbf{k}}|l,\alpha)$,
which however remains finite at ${\mathbf{k}}\to 0$ for finite $\alpha$ and
$r_{c}$. Therefore, any short-ranged pairwise interaction must lead to an
incomplete screening which gives the divergence of the screening function,
$S({\mathbf{k}})/k^{2}$ at ${\mathbf{k}}\to 0$:

$\frac{S({\mathbf{k}})}{k^{2}}=\frac{k_{b}T\varepsilon_{0}}{\lambda^{2}k^{2}+\hat{\nu}({\mathbf{k}}|l,\alpha)\varepsilon_{0}k^{2}}+{\cal
O}(k^{2}).$ (15)

When the simulation box is large enough to evaluate $S({\mathbf{k}})$ of Eq.
(3) at a sufficient small $k$, $S({\mathbf{k}})/k^{2}$ for the simulated
conducting fluid must approach the infinite.

Figure 2: The charge structure factors and screening functions (inset) reduced
by the ionic strength and $\varepsilon_{0}/\beta$ respectively for the ionic
fluid simulated with the three zero-multiple methods. $\beta=1/(k_{b}T)$. The
dash line and the solid lines (inset) indicate the analytic asymptotes at
large and small wavenumbers respectively. $\lambda=0.526\,\mbox{\AA}$

Figure 3: Same as Fig. 2 but for the SPC/E model of water. The screening
functions (inset) are additionally reduced by $\epsilon_{r}-1$ with
$\epsilon_{r}\simeq 70$.

Similar things happen to a molecular fluid with a dielectric constant
(relative permittivity) of $\epsilon_{r}$ ($\epsilon_{r}\geqslant 1$).
$\lambda^{2}k^{2}$ in Eq. (15) is substituted with $1/(\epsilon_{r}-1)$

$\frac{S({\mathbf{k}})}{k^{2}}=\frac{k_{b}T\varepsilon_{0}}{1/(\epsilon_{r}-1)+\hat{\nu}({\mathbf{k}}|l,\alpha)\varepsilon_{0}k^{2}}+{\cal
O}(k^{2}),$ (16)

such that $S({\mathbf{k}})$ approaches the correct limit of the dielectric
responseChandler (1977); Rodgers and Weeks (2009)

$S({\mathbf{k}})=k_{b}T\varepsilon_{0}(\epsilon_{r}-1)k^{2}/\epsilon_{r}+{\cal
O}(k^{4}),$ (17)

when $\hat{\nu}({\mathbf{k}}|l,\alpha)$ approaches $1/(\varepsilon_{0}k^{2})$
at $\alpha\to 0$ and $r_{c}\to\infty$. The screening function evaluated in a
simulation using the short-ranged $\nu({\mathbf{r}}|l,\alpha)$ now approaches
the limit $k_{b}T\varepsilon_{0}(\epsilon_{r}-1)$. This limit, away from
divergence for any finite $\epsilon_{r}$, is expected because there is no
unscreened free charge any more in the molecular fluid.

The above analysis emphasizes that the SPMF condition imposed on the pairwise
interaction is a direct consequence of the screening effect of the Coulomb
force. Violation of the SPMF condition leads to improper behaviors of charge
correlations characterized by $S({\mathbf{k}})/k^{2}$ at small wavenumbers. In
light of these theoretical predications, we carry out simulations of the model
ionic fluidHu and Weeks (2010) at $T=5000\,K$ and the SPC/E waterBerendsen
_et al._ (1987) at $T=298.15\,K$ with one length of the simulation box
$930\,\mbox{\AA}$ and $939.6\,\mbox{\AA}$, respectively. In the zm2/zm3
methods, $r_{c}=30\,\mbox{\AA}$ and $r_{c}=12\,\mbox{\AA}$ for the ionic fluid
and the water respectively. In the zm0 method, $\alpha=1/4.5\,\mbox{\AA}^{-1}$
as in other simulations of waterHu and Weeks (2010); Remsing  _et al._ (2016);
Gao  _et al._ (2020). The ionic strengths $I=0.00859\,e^{2}\mbox{\AA}^{-3}$
and $I=0.0359\,e^{2}\mbox{\AA}^{-3}$ for the ionic fluid and the water
respectively. The Debye length of the ionic fluid is
$\lambda=0.526\,\mbox{\AA}$. Complete details of the simulations are available
to the publicnot .

Figs. 2 and 3 display the computed $S({\mathbf{k}})$ and
$S({\mathbf{k}})/k^{2}$ at discrete wavenumbers for both systems. At extremely
large $k$, $S({\mathbf{k}})$ approaches the corresponding ionic strength
indicating the self-correlation. In addition, the damped oscillation of
$S({\mathbf{k}})$ for the SPC/E water approaches the intra-molecular
correlation

$\frac{S({\mathbf{k}})}{I}\sim 1+\frac{4q_{\rm O}q_{\rm H}}{q_{\rm
sum}}\frac{\sin(kd_{\rm OH})}{kd_{\rm OH}}+\frac{2q_{\rm H}^{2}}{q_{\rm
sum}}\frac{\sin(kd_{\rm HH})}{kd_{\rm HH}},$ (18)

where $q_{\rm O}$, $q_{\rm H}$, $d_{\rm OH}$ and $d_{\rm HH}$ are the
parameters of charges and bond lengths in the SPC/E modelBerendsen  _et al._
(1987); not . $q_{\rm sum}=q_{\rm O}^{2}+2q_{\rm H}^{2}$, is the sum of the
squares of the charges in one molecule. $S({\mathbf{k}})$ from all methods are
indistinguishable at the wavenumbers corresponding to distances ($2\pi/k$)
smaller than the local correlation lengths. However, distinct features between
the SPMF condition satisfied and violated methods are found at small $k$ for
both the ionic fluid and the water. These characteristics are in excellent
agreement with the analytical predictions from Eqs. (15) and (16) with
$\hat{\nu}({\mathbf{k}}|l,\alpha)$ determined by Eqs. (7) to (10) and the
specified parameters of the systems.

In conclusion, the SPMF condition provides a simple concept to connect the
electrostatic correlations among Coulomb and simulated Coulomb systems.
Exceptions are simulations employing techniques that are not well represented
by a pairwise interaction. These useful methods include the local maxwell
approachMaggs (2002), the random batch methodJin  _et al._ (2021); Liang  _et
al._ (2021) and the fast growing machine-learning techniquesZhang  _et al._
(2018); Wang  _et al._ (2018). While analytic understandings of the mentioned
methods have yet to be achieved, the present work suggests that it is
convenient to validate them by calculating correlations in the reciprocal
space such as $S({\mathbf{k}})$, the usual structure factor and intermediate
scattering functionsHansen and McDonald (2006). We note that numerical
evidence does show that the local maxwell approach behaves correctly at small
wavenumbersMaggs (2002). Further investigations at all wavelengths for both
time-averaged and time-displaced correlations are expected. We finally hope
that the SPMF condition for interfacial and bulk electrostatics helps to
provide a simple starting point when developing methods to simulate a complex
condensed phase.

This work was supported by the NSFC (Grant No. 21873037)

## References

  * Allen and Tildesley (2017) M. P. Allen and D. J. Tildesley,  _Computer Simulation of Liquids_ , 2nd ed. (Oxford University Press, Oxford, UK, 2017).
  * Yi  _et al._ (2015) Shasha Yi, Cong Pan, and Zhonghan Hu, “Accurate treatments of electrostatics for computer simulations of biological systems: A brief survey of developments and existing problems,” Chin. Phys. B 24, 120201 (2015).
  * Lowe  _et al._ (2018) Benjamin M. Lowe, Chris-Kriton Skylaris, Nicolas G. Green, Yasushi Shibuta, and Toshiya Sakata, “Calculation of surface potentials at the silica–water interface using molecular dynamics: Challenges and opportunities,” Japanese J. App. Phys. 57, 04FM02 (2018).
  * Maggs (2002) A. C. Maggs, “Dynamics of a local algorithm for simulating coulomb interactions,” J. Chem. Phys. 117, 1975–1981 (2002).
  * Fukuda (2013) Ikuo Fukuda, “Zero-multipole summation method for efficiently estimating electrostatic interactions in molecular system,” J. Chem. Phys. 139, 174107 (2013).
  * Wang  _et al._ (2016a) Han Wang, Haruki Nakamura, and Ikuo Fukuda, “A critical appraisal of the zero-multipole method: Structural, thermodynamic, dielectric, and dynamical properties of a water system,” J. Chem. Phys. 144, 114503 (2016a).
  * Wang  _et al._ (2016b) Han Wang, Xingyu Gao, and Jun Fang, “Multiple staggered mesh ewald: Boosting the accuracy of the smooth particle mesh ewald method,” J Chem. Theory Comput. 12, 5596–5608 (2016b).
  * Liang  _et al._ (2015) Yihao Liang, Zhenli Xu, and Xiangjun Xing, “A multi-scale monte carlo method for electrolytes,” New J. Phys. 17, 083062 (2015).
  * Girotto  _et al._ (2017) Matheus Girotto, Alexandre P. dos Santos, and Yan Levin, “Simulations of ionic liquids confined by metal electrodes using periodic green functions,” J. Chem. Phys. 147, 074109 (2017).
  * Bakhshandeh  _et al._ (2018) Amin Bakhshandeh, Alexandre P. dos Santos, and Yan Levin, “Efficient simulation method for nano-patterned charged surfaces in an electrolyte solution,” Soft Matter 14, 4081–4086 (2018).
  * Zhang  _et al._ (2018) Linfeng Zhang, Jiequn Han, Han Wang, Roberto Car, and Weinan E, “Deep potential molecular dynamics: A scalable model with the accuracy of quantum mechanics,” Phys. Rev. Lett. 120, 143001 (2018).
  * Wang  _et al._ (2018) Han Wang, Linfeng Zhang, Jiequn Han, and Weinan E, “Deepmd-kit: A deep learning package for many-body potential energy representation and molecular dynamics,” Computer Phys. Commun. 228, 178–184 (2018).
  * Urano  _et al._ (2020) Ryo Urano, Wataru Shinoda, Noriyuki Yoshii, and Susumu Okazaki, “Exact long-range coulombic energy calculation for net charged systems neutralized by uniformly distributed background charge using fast multipole method and its application to efficient free energy calculation,” J. Chem. Phys. 152, 244115 (2020).
  * Yuan  _et al._ (2021) Jiaxing Yuan, Hanne S. Antila, and Erik Luijten, “Particle–particle particle–mesh algorithm for electrolytes between charged dielectric interfaces,” J. Chem. Phys. 154, 094115 (2021).
  * Jin  _et al._ (2021) Shi Jin, Lei Li, Zhenli Xu, and Yue Zhao, “A random batch ewald method for particle systems with coulomb interactions,” SIAM J. Sci. Comput. 43, B937–B960 (2021).
  * Hu (2014a) Zhonghan Hu, “Symmetry-preserving mean field theory for electrostatics at interfaces,” Chem. Commun. 50, 14397–14400 (2014a).
  * Yi  _et al._ (2017a) Shasha Yi, Cong Pan, Liming Hu, and Zhonghan Hu, “On the connections and differences among three mean-field approximations: a stringent test,” Phys. Chem. Chem. Phys. 19, 18514–18518 (2017a).
  * Pan  _et al._ (2017) Cong Pan, Shasha Yi, and Zhonghan Hu, “The effect of electrostatic boundaries in molecular simulations: symmetry matters,” Phys. Chem. Chem. Phys. 19, 4861 (2017).
  * Pan  _et al._ (2019) Cong Pan, Shasha Yi, and Zhonghan Hu, “Analytic theory of finite-size effects in supercell modelling of charged interfaces,” Phys. Chem. Chem. Phys. 21, 14858 (2019).
  * Hu (2014b) Zhonghan Hu, “Infinite boundary terms of ewald sums and pairwise interactions for electrostatics in bulk and at interfaces,” J. Chem. Theory Comput. 10, 5254–5264 (2014b).
  * Yi  _et al._ (2017b) Shasha Yi, Cong Pan, and Zhonghan Hu, “Note: A pairwise form of the ewald sum for non-neutral systems,” J. Chem. Phys. 147, 126101 (2017b).
  * Arnold  _et al._ (2002) A Arnold, Jason de Joannis, and Christian Holm, “Electrostatics in periodic slab geometries. i,” J. Chem. Phys. 117, 2496–2502 (2002).
  * Mináry  _et al._ (2002) Peter Mináry, Mark E. Tuckerman, Katianna A. Pihakari, and Glenn J. Martyna, “A new reciprocal space based treatment of long range interactions on surfaces,” J. Chem. Phys. 116, 5351–5362 (2002).
  * Pan and Hu (2014) Cong Pan and Zhonghan Hu, “Rigorous error bounds for ewald summation of electrostatics at planar interfaces,” J. Chem. Theory Comput. 10, 534–542 (2014).
  * Bengtsson (1999) Lennart Bengtsson, “Dipole correction for surface supercell calculations,” Phys. Rev. B 59, 12301–12304 (1999).
  * Zhang  _et al._ (2020) Chao Zhang, Thomas Sayer, Jürg Hutter, and Michiel Sprik, “Modelling electrochemical systems with finite field molecular dynamics,” J. Phys.: Energy 2, 032005 (2020).
  * D. Wolf and Eggebrecht (1999) S.R. Phillpot D. Wolf, P. Keblinski and J. Eggebrecht, J. Chem. Phys. 110, 8254 (1999).
  * Chen  _et al._ (2004) Yng-gwei Chen, Charanbir Kaur, and John D. Weeks, “Connecting systems with short and long ranged interactions: Local molecular field theory for ionic fluids,” J. Phy. Chem. B 108, 19874–19884 (2004).
  * Hansen and McDonald (2006) J. P. Hansen and I. R. McDonald,  _Theory of simple liquids_ , 3rd ed. (Academic Press, Inc., Amsterdam, 2006).
  * Caillol (1992) J. M. Caillol, “Asymptotic behavior of the pair‐correlation function of a polar liquid,” J. Chem. Phys. 96, 7039–7053 (1992).
  * de Leeuw  _et al._ (1980) S. W. de Leeuw, J. W. Perram, and E. R. Smith, “Simulation of electrostatic systems in periodic boundary conditions. i. lattice sums and dielectric constants,” Proc. R. Soc. London, Ser. A Math. Phys. Sci. 373, 27–56 (1980).
  * Caillol (1994) Jeanichel Caillol, “Comments on the numerical simulations of electrolytes in periodic boundary conditions,” J. Chem. Phys. 101, 6080–6090 (1994).
  * Stillinger and Lovett (1968) Frank H. Stillinger and Ronald Lovett, “General restriction on the distribution of ions in electrolytes,” J. Chem. Phys. 49, 1991–1994 (1968).
  * Chandler (1977) David Chandler, “The dielectric constant and related equilibrium properties of molecular fluids: Interaction site cluster theory analysis,” J. Chem. Phys. 67, 1113–1124 (1977).
  * Rodgers and Weeks (2009) Jocelyn M. Rodgers and John D. Weeks, “Accurate thermodynamics for short-ranged truncations of coulomb interactions in site-site molecular models,” J. Chem. Phys. 131, 244108 (2009).
  * Hu and Weeks (2010) Zhonghan Hu and John D. Weeks, “Efficient solutions of self-consistent mean field equations for dewetting and electrostatics in nonuniform liquids,” Phys. Rev. Lett. 105, 140602 (2010).
  * Berendsen  _et al._ (1987) H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, “The missing term in effective pair potentials,” J. Phys. Chem. 91, 6269–6271 (1987).
  * Remsing  _et al._ (2016) R. C. Remsing, S. Liu, and J. D. Weeks, “Long-ranged contributions to solvation free energies from theory and short-ranged models,” Proc. Natl. Acad. Sci. USA 117, 2819–2826 (2016).
  * Gao  _et al._ (2020) Ang Gao, Richard C. Remsing, and John D. Weeks, “Short solvent model for ion correlations and hydrophobic association,” Proc. Natl. Acad. Sci. 117, 1293–1302 (2020).
  * (40) The complete simulation code has been uploaded to github.com/zhonghanhu1981/mdcode/zm.ionic.water.tar.gz.
  * Liang  _et al._ (2021) Jiuyang Liang, Zhenli Xu, and Yue Zhao, “Random-batch list algorithm for short-range molecular dynamics simulations,” J. Chem. Phys. 155, 044108 (2021).


Efficacy versus abundancy: comparing vaccination schemes

Omar El Deeb* 1,2  and Maya Jalloul 3

1 Department of Natural Sciences, Lebanese American University, Beirut,
Lebanon

2 Faculty of Technology, Lebanese University, Aabey, Lebanon

3 Department of Economics, Lebanese American University, Beirut, Lebanon

These authors contributed equally to this work.

*<EMAIL_ADDRESS>
## Abstract

###### Abstract

We introduce a novel compartmental model accounting for the effects of vaccine
efficacy, deployment rates and timing of initiation of deployment. We simulate
different scenarios and initial conditions, and we find that higher abundancy
and rate of deployment of low efficacy vaccines lowers the cumulative number
of deaths in comparison to slower deployment of high efficacy vaccines. We
also forecast that, at the same daily deployment rate, the earlier
introduction of vaccination schemes with lower efficacy would also lower the
number of deaths with respect to a delayed introduction of high efficacy
vaccines, which can however, still achieve lower numbers of infections and
better herd immunity.

keywords: Compartmental model; Vaccine; COVID-19; Epidemiology.

##  1 Introduction

Humanity has been struggling with viral pandemics and infectious diseases that
caused great catastrophes throughout the entire history of mankind. Since the
Athenian, Antonine and Justinian Plagues to the Black Death, Spanish Flu,
Cholera, and Smallpox up until the HIV pandemic, SARS, Swine flu, Ebola
outbreak and, most recently, the Coronavirus infection, pandemics have been a
major source of disease, death, economic crises and political turmoils [1, 2,
3]. Nevertheless, those outbreaks also lead to major discoveries and
advancements in sciences and public health, especially in medicine,
pharmaceuticals, vaccines and development of public policies [4, 5, 6].

The Coronavirus disease (COVID-19) associated to the severe acute respiratory
syndrome coronavirus 2 (SARS-CoV-2) continues as a main cause of
hospitalization and death and as a main public health risk since the first
case was registered in Wuhan, China in December 2019 [7]. It was declared as a
global pandemic on March 11, 2020 [8] after spreading from China into several
Asian and European countries and consequently into tens of countries across
the world. By December 22, 2021, the world has suffered about 278 million
registered infections and 5.4 million deaths associated to COVID-19 in 224
different countries, territories and entities [9]. Health systems were put
under enormous pressure, and the consequent mitigation measures and associated
closures and lockdowns taken across the world for lengthy periods of time
contributed to a global economic slowdown and recession in several countries
[10, 11], in addition to a heavy negative impact on education, employment,
tourism and social activities [12, 13, 14] exacerbating existing political
discontent and stirring political unrest [15].

Similarly to the other respiratory pathogens, airborne transmission of SARS-
CoV-2 occurs by inhaling droplets loaded with the virus emitted by infectious
people. Infection can also occur through mucous membranes like the eyes, nose
or mouth [16].

Many variants of the virus emerged with an increased risk to global public
health. The World Health Organization (WHO) characterized variants of interest
(VOI) and variants of concern (VOC) to track and monitor [17]. The main
variants of concern have stronger capabilities of transmissibility and health
risks, and they are labeled by the WHO as the Alpha, Beta, Gamma, Delta and
Omicron variants which were first documented in the United Kingdom, South
Africa, Brazil and India, respectively [18]. The recent emergence of variants
creates a major cause of concern since they can lead to an epidemic rebound
especially with the possibility of vaccine resisting and the emergence of
deadlier and/or more transmissible future variants. Increased viral
transmission means higher probability of emergence of variants as well [19].

Several types of vaccines were developed and approved by global (WHO) or by
national health agencies, using different techniques like the viral vector
vaccines (Oxford/AstraZeneca, Sputnik V/Gamaleya, Janssen), genetic vaccines
(Pfizer/BioNTech, Moderna), inactivated vaccines (Sinovac, Sinopharm, Bharat)
and protein vaccines (Novavax, Sanofi) [20]. By December 22, 2021, a total of
about 8.85 billion doses of those vaccines were administered with main
concentration in China, USA, India and Europe [21]. Other vaccines, inhaled
aerosol dry vaccines and antiviral drugs are under extensive research and
development [22, 23, 24]. Antibody therapies are under development as well and
they offer an effective treatment for the most seriously ill patients, but
they remain expensive and in short supply [25].

A global roll-out of vaccines is needed to guarantee a swift elimination of
the infection, and for this end vaccines need to be available, affordable, and
accessible at a global scale. In high-income countries, large stocks were
produced and purchased, combined with active logistic and public health
resources, but in low-income countries, vaccination is still slow mainly due
to insufficiency of vaccines. The WHO has called for more equity and stronger
support for its COVAX initiative in order to supply wider vaccine access for
poor countries, but statistics show that the gap remains deep as many rich
countries have already bought multiple times of doses needed per person [26].
A similar gap also arises on the national level as health illiteracy,
religious beliefs and sometimes political partisanship are slowing down
vaccinations and creating irregular immunity rates among different regions or
groups [27]. As such, the global efforts to stop the spread remains far from
successful, and the priority should be focused to optimize the use, deployment
and timing of available vaccines, whether those recognized by the WHO or by
national health agencies across the globe in order to prevent further deaths
as well as emergence of new variants.

The spread patterns of infectious diseases have been and are being actively
studied and simulated using several mathematical and computational models in
the interdisciplinary fields of health economics, physical biology and medical
physics, applied mathematics and epidemiology . A wide variety of techniques
were employed ranging from compartmental models, Agent-Based models (ABM),
spatio-temporal analysis to data-driven analysis and artificial intelligence
[28, 29]. The first compartmental model was the famous Kermack–McKendrick
Susceptible-Infectious-Removed (SIR) model which divides the total population
under study into three compartments where agents would move from one
compartment to another upon infection or removal through recovery or death
[30]. Numerous variations of this model were introduced to account for
characteristics and dynamics of different diseases as well as in simulating
similar interactions and spreads in the social and behavioral sciences [31,
32, 33]. The research on COVID-19 spread extensively employed SIR models that
were improved to account for other effects like exposure, travel, quarantine,
vaccination, public measures and many other alterations to best describe the
underlying features of the contagion [34, 35, 36, 37, 38], as well as general
and country specific ABMs, spatio-temporal and data driven studies [39, 40,
41, 42, 43, 44, 45, 46, 47].

However, with the mass introduction of various globally and nationally
recognized vaccines with varying efficacies, and with the varying degrees of
availability and logistic capabilities across the globe, an important question
arises for decision makers: which vaccine to choose to deploy in a certain
country from the available options in that country, given their efficacies,
their availability timeline and expected deployment rate of each. This is a
major concern for public health officials aiming at minimizing the number of
deaths and/or the number of active infections, that was not thoroughly
discussed and analyzed in the literature. This paper studies the trade-off
between vaccine efficacy and abundance and then between efficacy and time of
availability, and the corresponding expected outcomes for deaths and
infections.

The paper is organized as follows: section 2 presents the theoretical model,
section 3 puts forward the results and discussions about the efficacy,
deployment and timing rates while section 4 concludes the paper.

Fig 1: A schematic diagram of the $SS_{V}II_{V}RD$ model showing the six
compartments of the model (Susceptible, Vaccinated susceptible, Infectious,
Vaccinated Infectious, Recovered and Dead) together with their transfer
dynamics and the coefficients relating transfers among different compartments.

##  2 Theoretical Model

We introduce a novel compartmental model to account for the effect of
vaccines, taking into account different vaccination deployment rates ($v$),
efficacies ($e$) under several scenarios of infection spread rates represented
by the reproductive number $(R_{t})$, starting with different initial
conditions in relation to different numbers of infections, deaths and immune
populations in different countries. The $SS_{V}II_{V}RD$ model consists of 6
compartments: the susceptible unvaccinated population $S$, the susceptible
vaccinated population $S_{V}$ who might become infected during the assumed
effective immunity period depending on the efficacy of the vaccine, the
infectious unvaccinated population $I$, the infectious vaccinated population
$I_{V}$, the recovered population $R$ who are either fully protected by the
vaccine during its effective protection time range or who have already
recovered from the infection and still possess immunity, and finally the dead
$D$. Each compartment is normalized with respect to the total population so
that $N=S+S_{V}+I+I_{V}+R+D=1$.

The $6$ compartments are interlinked (Figure 1), and their dynamics can be
modeled by a system of coupled linear ordinary differential equations
described below.

$\begin{cases}\frac{dS}{dt}&=-\beta\left(I+I_{V}\right)S-evS+\frac{evS}{T_{VI}}+\frac{\gamma_{R}I+\gamma_{RV}I_{V}}{T_{II}}\\\
\frac{dS_{V}}{dt}&=\left(1-e\right)vS-\beta^{\prime}\left(I+I_{V}\right)S_{V}\\\
\frac{dI}{dt}&=\beta\left(I+I_{V}\right)S-\left(\gamma_{D}+\gamma_{R}\right)I\\\
\frac{dI_{V}}{dt}&=\beta^{\prime}\left(I+I_{V}\right)S_{V}-\left(\gamma_{DV}+\gamma_{RV}\right)I_{V}\\\
\frac{dR}{dt}&=evS+\gamma_{R}I+\gamma_{RV}I_{V}-\frac{evS}{T_{VI}}-\frac{\gamma_{R}I+\gamma_{RV}I_{V}}{T_{II}}\\\
\frac{dD}{dt}&=\gamma_{D}I+\gamma_{DV}I_{V}\end{cases}$ (1)

The parameter $\beta$ is the incidence rate between susceptible $(S)$ and all
infectious people $(I+I_{V})$ and it can be determined using the reproductive
number of the infectious spread $R_{t}$ and the average rate of recovery of
unvaccinated people$\gamma_{R}$. In this sense, $\beta=\gamma_{R}R_{t}$. The
rate of recovery of vaccinated people is given by $\gamma_{RV}$ and the
incidence rate between vaccinated susceptible people and infectious people is
given by $\beta^{\prime}=\gamma_{RV}R_{t}$. The rate of death of unvaccinated
people is $\gamma_{D}$ and that of vaccinated ones is $\gamma_{DV}$. Finally,
$T_{VI}$ represents the average time duration of acquired vaccine immunity and
$T_{II}$ represents the average duration of acquired post infection immunity.
The change in the total population is $\frac{dN}{dt}=0$, assuming that $N$
remains constant during that period, neglecting other changes due to natural
growth, immigration, etc…

The transition dynamics between different compartments can be summarized
according to the following:

  * •

The susceptible $S$ may become infected $I$ upon contact with the infectious
(vaccinated or non-vaccinated) at rate $\beta$. They may also become recovered
$R$ upon receiving a vaccine with efficacy $e$ given at a daily rate $v$. In
addition, recovered people $R$ would become susceptible again some time after
recovery. Mainly, recovered people due to infection would become susceptible
again in an infection immunity period $T_{II}$ and recovered people due to
vaccination would become susceptible again in a vaccine immunity period
$T_{VI}$.

  * •

The vaccinated susceptible population $S_{V}$ are the susceptibles $S$ who got
vaccinated but are still susceptible to infection, and may become infectious
$I_{V}$ upon a meeting rate $\beta^{\prime}$ with other infectious people.

  * •

The infectious compartment $I$ is populated by unvaccinated susceptibles $S$
encountering other infectious agents at a rate $\beta$, while exit from this
compartment is attributed to deaths at rate $\gamma_{D}$ into compartment $D$
and recoveries at rate $\gamma_{R}$ into compartment $R$.

  * •

The vaccinated infectious compartment $I_{V}$ includes vaccinated susceptibles
$S_{V}$ catching the disease at a rate $\beta^{\prime}$ and diminished by
people dying at rate $\gamma_{DV}$ into $D$ and people recovering a rate
$\gamma_{RV}$ into $R$.

  * •

The recovered population $R$ is formed by effectively vaccinated susceptible
people coming from $S$ and people surviving the infection at rates
$\gamma_{R}$ and $\gamma_{RV}$ from infectious and vaccinated infectious
populations $I$ and $I_{V}$. Simultaneously, recovered people would eventually
lose their acquired immunity on average periods of $T_{II}$ after infection
and $T_{VI}$ after vaccination, thus would exit the recovered compartment $R$
back to the susceptibles $S$.

  * •

Finally, the dead $D$ increase at death rates $\gamma_{D}$ and $\gamma_{DV}$
among infectious and vaccinated infectious populations $I$ and $I_{V}$. In
this model, we only consider the deaths caused by COVID-19, disregarding net
population changes (natural deaths as well as natural births). Consequently,
the dead compartment gets populated from the infected and vaccinated infected
populations.

The vaccination deployment rate $v$ depends on the available supply of the
vaccine as well as the logistical capability and the popular demand at a given
time [48]. In this paper it is assumed to be equivalent to the daily
vaccination percentage of the susceptible population, with different scenarios
representing slow, moderate and fast deployment rates. The vaccine efficacy
varies among different employed vaccines as well as in relation to new
emerging variants in addition to the possibility of supplying a single dose of
a double dosed vaccine under short supply. The model accounts for a range of
scenarios with low, moderate and high efficacy corresponding to the former
cases. We also inspect those scenarios under different reproductive rates
$R_{t}$ which depend on different values of mitigation measures related to
social distancing, protective masks, sanitation, and other factors that alter
the rate of infection spread. We assume different scenarios of low, moderate,
high and alternating (where $R_{t}$ varies in a periodic pattern between high
and low extrema) reproductive rates.

The numerical values of these parameters are taken in relation to available
data and research. The efficacy values vary between $0.5\leq e\leq 0.95$
depending on the available vaccines [49, 50, 51, 52, 53, 54, 55]. The rate of
recovery for vaccinated people $\gamma_{R}=\frac{1}{14}$ in relation of an
average of $14$ days needed for recovery [56], while the death rate amounts to
around $2\%$ [9] of the infected which leads to
$\gamma_{D}=\frac{\gamma_{R}}{50}$. We also assume that those who are infected
after vaccination would need a similar time of recovery, despite the fact that
their symptoms would be much reduced [57, 58], thus $\gamma_{RV}=\gamma_{R}$
but their death rate, as studies reveal, would be considerably lower by
$70-85\%$ [58, 59], which is modeled through
$\gamma_{DV}=\frac{\gamma_{D}}{5}$. We take $T_{VI}=90$ days and $T_{II}=360$
days to represent expected periods of acquired immunity after recovery and
after vaccination respectively [60]. The deployment rate $v$ is assigned
hypothetical values varying between $0.1\%$ of the susceptible population per
day for slow vaccination rollout and $1.5\%$ for the highest pace of
vaccination rollout.

The reproduction rate $R_{t}$ is assigned values of $0.7,1.1$ and $1.5$
corresponding to low, medium and high reproduction rates, in addition to a
fourth scenario where $R_{t}$ alternates sinusoidally between $0.4$ and $2$
according to the relation $R_{t}=0.8\cos(\frac{t}{15})+1.2$ where $t$
represents the time in days, with a period of around $94$ days to simulate the
effect of the consecutive waves of the spread of the infection, as observed
empirically [61, 62]. The periodicity of the function $\cos(at+b)$ is
mathematically defined to be $\frac{2\pi}{a}$. The parameter $\frac{1}{15}$ in
our model simulates a periodicity of $30\pi$ which approximately corresponds
to a $94$ day period. Its dimension is time given in the unit of days. Thus,
this forecast simulates the recurrence of waves of high then low infection
rates on an average period of 94 days. The growth rate of the active cases of
COVID-19 infections is studied in the United Kingdom, South Africa and Brazil
[63, 64]. These articles present the progression of the registered values of
$R_{t}$ since the start of the spread until recently. In both countries,
$R_{t}$ assumed values that varied between $0.5$ up to $3$ depending on the
mitigation measures and closures. The registered data constitutes a motivation
to simulate the possible scenarios of infection spreads in presence of
vaccines for representative values of $R_{t}$ corresponding to low, middle and
high rates of spread.

We simulate this model under different combinations of efficacy $e$,
deployment rate $v$ and reproductive number $R_{t}$ to analyze the
corresponding cumulative numbers of infected, recovered and dead populations.

Fig 2: The relative numbers of active cases, infected vaccinated cases, immune
population and total cumulative deaths are shown for a low reproductive rate
$R=0.7$ under initial conditions of: $73.1\%$ vaccinated population $2.0\%$
currently infected and $0.23\%$ dead (upper row, United Kingdom) and $29.2\%$
vaccinated population $0.02\%$ currently infected and $0.16\%$ dead (lower
row, South Africa) for nine different vaccination scenarios. The forecast
corresponds to vaccine efficacy $(e)$ and daily deployment rate $(v)$ of:
$e=92\%,v=0.5\%$ (in brown), $e=92\%,v=0.3\%$ (in black), $e=92\%,v=0.1\%$ (in
blue), $e=72\%,v=0.7\%$ (in gray), $e=72\%,v=0.5\%$ (in red), $e=72\%,v=0.3\%$
(in green), $e=55\%,v=1.5\%$ (in orange), $e=55\%,v=1\%$ (in purple) and
$e=55\%,v=0.7\%$ (in yellow).

##  3 Results and discussions

###  3.1 Efficacy vs abundance

We simulate the theoretical model introduced in (1) to determine relative
numbers of active cases, infected vaccinated cases, immune population and
total cumulative deaths under nine different combinations of efficacies $(e)$
of available vaccines and their deployment rates $(v)$ given by:
$e=92\%,v=0.5\%$, $e=92\%,v=0.3\%$, $e=92\%,v=0.1\%$, $e=72\%,v=0.7\%$,
$e=72\%,v=0.5\%$, $e=72\%,v=0.3\%$, $e=55\%,v=1.5\%$, $e=55\%,v=1\%$ and
$e=55\%,v=0.7\%$ in order to compare the levels of infection, death and
immunity between higher efficacy vaccines at lower abundance and lower
efficacy vaccines with more abundance or deployment rate, together with middle
values between them.

#### Case Study

We simulate in this paper the actual data of infections, deaths and
vaccination in two countries with high vaccination rate (United Kingdom) and
low vaccination rate (South Africa), according to data publically available in
[9]. The initial conditions corresponding to the United Kingdom are: $73.1\%$
vaccinated population $2.0\%$ currently infected and $0.23\%$ dead, while
those of South Africa are:$29.2\%$ vaccinated population $0.02\%$ currently
infected and $0.16\%$ dead. In this sense, the results of this forecast are
not country specific but are of global significance. They can be applied to
any other country with similar conditions. We repeat this simulation for
different levels of infection spread modeled through low, medium, high and
alternating reproductive rates $R_{t}$ defined before. The results are
displayed in Figures 2,3,4 and 5, respectively.

Figure 2 shows that when the reproductive rate is low ($R_{t}=0.7$), the
number of active cases would fall down quickly under all vaccination schemes,
for various efficacies and deployment rate, both in countries that are in the
middle or early stages of the spread (upper and lower figures respectively).
The number of infected vaccinated people is the highest for the low efficacy
vaccine. Herd immunity would be achieved optimally under fast deployment of
the high efficacy vaccine and the middle efficacy vaccine, while the lowest
immunity is achieved under the slowest deployment of the high efficacy
vaccine. This shows that the deployment rate is an essential factor in
combination with efficacy for achieving immunity. The figure also shows that
the number of deaths is most reduced under the two fastest deployment schemes
of the low efficacy vaccine, followed by the fastest deployment rate of that
of middle efficacy. The latter result shows that to reduce deaths, the most
essential measure is to ensure the the availability and the fast deployment of
the vaccine, while reaching herd immunity and decreasing infections (thus
lowering the pressure on the health sector and decreasing the possibility of
appearance of new mutations) depends relatively more on its efficacy.

Fig 3: The figure shows the relative numbers of active cases, infected
vaccinated cases, immune population and total cumulative deaths for a medium
reproductive rate $R=1.1$ with the same two initial conditions and nine
vaccination scenarios of Figure 2.

In the case of a medium reproductive rate ($R_{t}=1.1$) corresponding to a
slow increase in the number of infections in absence of vaccine, Figure 3
shows that the number of active cases would decrease in countries that are in
their middle stages of infection, while they would keep on increasing for a
limited interval of time in countries that are in their early stages of
spread. The duration needed for the active cases to start decreasing is
shortest for the vaccine with fastest deployment rate (despite its low
efficacy), while that with the slowest deployment rate (with high efficacy)
needs more time to decrease the number of active cases. The highest number of
infected vaccinated people would be obtained through the rapidly deployed low
efficacy for countries in their middle stages of spread, while for early stage
countries, each scenario of deployment of the low efficacy vaccine would lead
to maximal infections under different times. In a similar manner to what
happens for low reproductive rate, immunity will be maximally obtained under
the deployment of high efficacy vaccines at high rates, while the least deaths
are achieved under the fastest deployment rates of vaccines.

Fig 4: The figure shows the relative numbers of active cases, infected
vaccinated cases, immune population and total cumulative deaths for a high
reproductive rate $R=1.5$ with the same two initial conditions and nine
vaccination scenarios of Figure 2.

For a country under a fast spread scenario simulated by $R_{t}=1.5$, Figure 4
shows that the number of active cases would fall in a middle spread stage
country, though it takes more time than that with lower $R_{t}$. However, in
early spread countries, the number of active cases will rise and peak after a
period of time before falling down after a relatively long time (around
$100-200$ days) of vaccination. The slowest vaccination schemes of the low and
middle efficacy vaccines would reach the highest peaks and need more time for
controlling the spread while the fastest vaccination schemes would achieve it
the quickest. The number of infected vaccinated people would be the largest in
the case of the low efficacy vaccine deployed rapidly in a middle spread stage
country, while for a country at early spread stages, the two slowest
deployments of the low efficacy vaccine lead to the highest infections among
the vaccinated population. Immunity will be maximally obtained under the
deployment of high efficacy vaccines at high rates, while the least deaths are
achieved under the fastest deployment rates of vaccines, like what was shown
before for low and medium spread rates, but with a much higher magnitude of
deaths in all scenarios due to a higher reproductive rate.

Fig 5: The figure shows the relative numbers of active cases, infected
vaccinated cases, immune population and total cumulative deaths for an
alternating reproductive rate, varying sinusoidally between $R=0.4$ and $R=2$,
with the same two initial conditions and nine vaccination scenarios of Figure
2.

In the more realistic case of an alternating reproductive number varying
between $0.4\leq R_{t}\leq 2$ and corresponding to consecutive patterns of low
and high spread waves, we realize that in a country at a middle spread stage,
all vaccination strategies would lead to bringing down the number of active
cases with a small peak arising after the return of the next wave, while in a
country at early stages of spread, the slow deployment of vaccines, whether
with high or low efficacies, would cause the re-emergence of several peaks of
infections and high number of active cases. Only the fast deployment of
vaccines of various efficacies would lead to curb the number of infections
just after the first peak of spread. The number of infected vaccinated people
would be the highest for the case of fast deployment of low efficacy vaccines
in countries with middle stages of spread while stronger peaks of infections
among the vaccinated people would occur under slower vaccination scenarios of
low efficacy vaccine during the second wave in countries with early spread
stages. The highest herd immunity levels would be attained under the scenarios
of fastest deployments of high and medium efficacy vaccines in both categories
of countries under early or middle stages of spread after $3-4$ months,
despite the early lead and the sharp rise in the immune population during the
first few weeks of fast deployment of low efficacy vaccines. The cumulative
number of deaths would rise in jumps corresponding to successive waves of
spread especially in early stage countries, but in both cases, the fastest
vaccination schemes of low efficacy vaccines would ultimately lead to the
minimal number of deaths while the slowest pattern of high efficacy vaccine
deployment would result in the highest cumulative rate of deaths. Those
simulations are displayed in Figure 5.

Fig 6: The figure shows the expected relative numbers of active cases,
infected vaccinated cases, immune population and cumulative deaths at a low
reproductive rate $R=0.7$ in relation to four cases simulating the effect of
efficacy and the time of initiation of the vaccination process: a vaccine of
efficacy $e=52\%$ to start deployment immediately (red), a vaccine of efficacy
$e=72\%$ to start deployment immediately (black), a vaccine of efficacy
$e=72\%$ to start deployment after 30 days (blue) and a vaccine of efficacy
$e=92\%$ to start deployment after 30 days (brown), all being deployed at an
equal rate of $0.7\%$ of the susceptible population per day. The upper row
corresponds to initial conditions of $40.6\%$ vaccinated population $1.8\%$
currently infected and $0.18\%$ dead while the lower row corresponds to
$3.1\%$ vaccinated population $0.46\%$ currently infected and $0.11\%$ dead.

###  3.2 Efficacy vs time of availability

An important aspect of vaccination strategies that we consider in this paper
is the question of vaccine efficacy versus time of availability or start of
deployment. In this simulation, we forecast the relative number of active
cases, infected vaccinated people, the immune population and the cumulative
number of deaths under low, medium and high reproductive rates $R=0.7,1.1$ and
$1.5$ respectively in Figures 6, 7 and 8. Four cases simulating the effect of
efficacy and the time of initiation of the vaccination process were
considered: a vaccine of low efficacy $e=52\%$ to start deployment immediately
(red), a vaccine of medium efficacy $e=72\%$ to start deployment immediately
(black), a medium efficacy vaccine with $e=72\%$ to start deployment after 30
days (blue) and a vaccine of high efficacy $e=92\%$ to start deployment after
30 days (brown), all being deployed at an equal rate of the susceptible
population per day. As in previous simulations, the upper row corresponds to
countries at middle stages of spread with initial conditions of $40.6\%$
vaccinated population $1.8\%$ currently infected and $0.18\%$ dead while the
lower row corresponds to countries at early stages with $3.1\%$ vaccinated
population $0.46\%$ currently infected and $0.11\%$ dead.

Under the circumstances of low reproductive rate $R_{t}=0.7$ depicted in
Figure 6, it is clear that the number of active cases will fall significantly
under the four efficacy-timing schemes in both early and middle spread stage
categories. The number of infected vaccinated people in both categories would
occur under the scenario of immediate deployment of low efficacy vaccines.
Herd immunity would be maximally attained through the adoption of high
efficacy vaccines being deployed with a one month period of delay while it
would be the lowest using low efficacy vaccines deployed immediately. On the
level of cumulative deaths, the least number of deaths is realized under the
scenario of immediate deployment of medium then low efficacy vaccines, while
the delay would raise the number of deaths even while using high efficacy
vaccines.

Fig 7: The figure shows the expected relative numbers of active cases,
infected vaccinated cases, immune population and cumulative deaths
corresponding to a medium reproductive rate $R=1.1$ with same initial
conditions, efficacies and time of initiation as those detailed in Figure 6.

Figure 7 simulates same efficacy-timing scenarios under medium reproductive
rate $R_{t}=1.1$. We notice that in countries in their middle stage of
infection spread, the number of active cases would fall quickly under all
scenarios, while in early stage countries, the late deployment scenarios would
cause a rise and a peak in infections before falling down significantly,
whereas the quick deployment of low or medium efficacy vaccines would lower
the active cases faster during the first few months. We also notice that in
about $100$ days, the number of active cases due to a delayed high efficacy
vaccine will catch up and fall below the expected active cases under an
immediately deployed low efficacy vaccine. The number of infected vaccinated
people is the highest for the lowest efficacy vaccine and the lowest for the
highest efficacy vaccine in both country categories. Similarly, for both
categories, herd immunity is maximally attained using the delayed high
efficacy vaccine rather than the immediate low efficacy one which provides the
lowest percentile of immune population. However, regarding deaths, the lowest
numbers of deaths are attributed to the immediate deployment of medium then
low efficacy vaccines, while a delayed deployment will cause more deaths even
while using high and medium efficacy vaccines.

Fig 8: The figure shows the expected relative numbers of active cases,
infected vaccinated cases, immune population and cumulative deaths
corresponding to a high reproductive rate $R=1.5$ with same initial
conditions, efficacies and time of initiation as those detailed in Figure 6.

Under the scenario of high reproductive number simulated by $R_{t}=1.5$, the
number of active cases would fall slowly in countries that are in their middle
stages of vaccination and spread, while it would rise and peak under all
vaccination scenarios in early stage countries. The highest peaks are
attributed to delayed vaccinations even though infections would fall rapidly
once vaccination starts, and eventually the corresponding number of active
cases would fall below that of low efficacy vaccine being deployed
immediately. Immediate deployment also helps to flatten the curve on
infections, thus reducing the expected peak number of cases. The maximum
number of infected vaccinated people would correspond low efficacy vaccine
deployed immediately while the lowest corresponds to the high efficacy vaccine
deployed late, for both country categories. Similarly, herd immunity is
maximally achieved by the high efficacy vaccine despite being introduced late
while the lowest level of immunity is caused by the low efficacy vaccine
despite early introduction. The immediate deployment of medium efficacy
vaccine minimizes the number of deaths, followed by the immediate low efficacy
vaccine, while late deployment raises the number of deaths for both country
categories, but with higher relative differences among the death outcomes for
middle stage countries.

###  3.3 Limitations

This study takes into account all parameters related to vaccine efficacy and
deployment rates under several infectious rates and initial conditions.
However, there is a limiting factor related to the maximal proportion of the
population who are willing or at least who would eventually take the vaccine.
Here, we took this into account indirectly by using a daily deployment rate
$v$ which is proportional to the susceptibles $S$, and not to the total
population $N$. In this sense, in the initial phases of vaccination, the
number of people taking the vaccine would be the highest, but as time
progresses, the number of susceptibles decreases, hence the daily number of
vaccinated people decreases. It is natural to assume that, as when a country
reaches a high level of vaccination, less people will be willing to get the
vaccine. If vaccination rates were only connected to abundance or logistic
infrastructure, they would have been linked to the total population $N$. It is
also important to mention that the reproduction number $R_{t}$ can differ for
different population groups. Its value is dominated by the group in which most
transmission occurs. In this sense, it is not a homogeneous parameter that
applies on a nationwide or social scale but depends on demographic, socio-
economic and other factors.

Due to various reasons ranging from religious and political beliefs, into non-
scientific and anti-vaxxer fears, there might be a sizable sector of the
society who would refuse to get vaccinated [27]. Vaccine hesitancy is not
directly simulated in the model, but it is indirectly represented through
relating the daily deployment rate to the number of susceptibles hence it
decreases as the number of vaccinated people increases. In our scenario,
immunity of this portion of the population would still be achieved through
infection rather than vaccination.

##  4 Conclusion

In this paper we introduced a general novel compartmental model accounting for
the vaccinated population, infected vaccinated population, active infections,
and deaths with various vaccine efficacies and vaccination deployment rates.

We simulated different scenarios and initial conditions, and we showed that
abundance and higher rate of deployment of low efficacy vaccines would lower
the cumulative number of deaths in comparison to slower deployment of high
efficacy vaccines. However, the high efficacy vaccines can better lower the
number of active cases and achieve faster and higher herd immunity.

We also discovered that at the same daily deployment rate, the earlier
introduction of vaccines with lower efficacy would also lower the number of
deaths with respect to a delayed introduction of high efficacy vaccines, which
can, however, lower the number of infections and attain higher levels of herd
immunity.

## References

  * 1. Huremovic, D., Brief History of Pandemics (Pandemics Throughout History). Psychiatry of Pandemics, 7-35 (2019), Springer, https://doi.org/10.1007/978-3-030-15346-5_2.
  * 2. Piret, J. and Boivin, G., Pandemics Throughout History. Front. Microbiol., 11:631736 (2021), doi: 10.3389/fmicb.2020.631736.
  * 3. Hays, J., Epidemics and pandemics: their impacts on human history (2005), ABC-CLIO, ISBN 978-1-85109-658-9.
  * 4. Porter, R., The Cambridge illustrated history of medicine (2001), Cambridge University Press, ISBN: 9780521002523.
  * 5. Stern, A.M. and Markel, H., The History Of Vaccines And Immunization: Familiar Patterns, New Challenges. Health Aff. 24-3 (2005), https://doi.org/10.1377/hlthaff.24.3.611.
  * 6. Plotkin, S.A. and Plotkin S.L., The development of vaccines: how the past led to the future, Nat Rev Microbiol. 9, 889–893 (2011). https://doi.org/10.1038/nrmicro2668.
  * 7. Report of the WHO-China Joint Mission on Coronavirus Disease 2019 (COVID-19), World Health Organization (2019), https://www.who.int/docs/default-source/coronaviruse/who-china-joint-mission-on-covid-19-final-report.pdf, (Accessed on August 20, 2021).
  * 8. WHO Director-General’s opening remarks at the media briefing on COVID-19 - 11 March 2020, World Health Organization (2020), https://www.who.int/director-general/speeches/detail/who-director-general-s-opening-remarks-at-the-media-briefing-on-covid-19—11-march-2020, (Accessed on June 7, 2021).
  * 9. Worldometer, COVID-19 coronavirus pandemic, https://www.worldometers.info/coronavirus/, (Accessed on August 20, 2021)
  * 10. Jackson, J., Weiss, M., Schwarzenberg, A. and Nelsson, R., Global Economic Effects of COVID-19, Homeland Security Digital Library (2020), https://search.bvsalud.org/global-literature-on-novel-coronavirus-2019-ncov/resource/en/grc-740379, (Accessed on July 19, 2021).
  * 11. Donthu, N. and Gustafsson, A., Effects of COVID-19 on business and research, J Bus Res. 117: 284–289 (2020), https://doi.org/10.1016/j.jbusres.2020.06.008.
  * 12. Editorial Board, Research and higher education in the time of COVID-19, Lancet 396, 583 (2020), https://doi.org/10.1016/S0140-6736(20)31818-3.
  * 13. Pepita, B., Schalkwyk, M., Mckee, M. and Labonte R., COVID-19 and the collapse of global trade: building an effective public health response, Lancet Planet. Health 5, 102-107 (2021), https://doi.org/10.1016/S2542-5196(20)30291-6.
  * 14. Tarricone, R. and Rognoni, C., What can health systems learn from COVID-19?, Eur Heart J Suppl. 22, 4-8 (2020), https://doi.org/10.1093/eurheartj/suaa185
  * 15. The pandemic has exacerbated existing political discontent, The Economist, https://www.economist.com/international/2021/07/31/the-pandemic-has-exacerbated-existing-political-discontent, (Accessed on July 19, 2021).
  * 16. Van Doremalen, N., Morris, D., Holbrook, M., Gamble, A., Williamson, B., et al., Aerosol and surface stability of SARS-CoV-2 as compared with SARS-CoV-1, N Engl J Med. 382, 1564-1567 (2020), doi.org/10.1056/NEJMc2004973.
  * 17. Tracking SARS-CoV-2 variants, World Health Organization (2021), https://www.who.int/en/activities/tracking-SARS-CoV-2-variants, (Accessed on July 16, 2021).
  * 18. The delta variant is the most dangerous SARS-CoV-2 mutation yet, The Economist (2021), https://www.economist.com/graphic-detail/2021/06/16/the-delta-variant-is-the-most-dangerous-sars-cov-2-mutation-yet, (Accessed on July 20, 2021).
  * 19. Fontanet, A., Autran, B., Lina, B., Kieny, M.P., Abdool Karim, S., and Sridhar, D., SARS-CoV-2 variants and ending the COVID-19 pandemic, Lancet 397, 592-594 (2021), https://doi.org/10.1016/S0140-6736(21)00370-6
  * 20. British Society for immunology, Types of vaccines for COVID-19 (2021), https://www.immunology.org/coronavirus/connect-coronavirus-public-engagement-resources/types-vaccines-for-covid-19, (Accessed on Jun 29, 2021).
  * 21. Mathieu, E., Ritchie, H., Ortiz-Ospina, E., Roser, M.,Hasell, J., Appel, C., et al. A global database of COVID-19 vaccinations. Nat Hum Behav. 5, 947–953 (2021). https://doi.org/10.1038/s41562-021-01122-8.
  * 22. Heida, R., Hinrichs, W., and Frijlink, H., Inhaled vaccine delivery in the combat against respiratory viruses: a 2021 overview of recent developments and implications for COVID-19 (2021), Expert Rev. Vaccines, https://doi.org/10.1080/14760584.2021.1903878.
  * 23. Kim, J.H., Marks, F. & Clemens, J.D., Looking beyond COVID-19 vaccine phase 3 trials. Nat Med 27, 205–211 (2021). https://doi.org/10.1038/s41591-021-01230-y.
  * 24. Tang, Z., Zhang, X., Shu, Y., Guo, M., Zhang, H., Tao, W., et al., Insights from nanotechnology in COVID-19 treatment, nanotoday 36, 101019 (2021), https://doi.org/10.1016/j.nantod.2020.101019.
  * 25. A new weapon in the war against SARS-CoV-2 has been found, The Economist (2021), https://www.economist.com/science-and-technology/2021/06/16/a-new-weapon-in-the-war-against-sars-cov-2-has-been-found. (Accessed on July 21, 2021).
  * 26. Mullard A., How COVID vaccines are being divided up around the world, Nature (2020), https://doi.org/10.1038/d41586-020-03370-6.
  * 27. America’s vaccination woes cannot be blamed only on politics, The Economist (2021), https://www.economist.com/united-states/2021/07/27/americas-vaccination-woes-cannot-be-blamed-only-on-politics. (Accessed on August 11, 2021).
  * 28. Brauer, F., Castillo-Chavez, C., & Feng, Z., Mathematical Models in Epidemiology, Texts in Applied Math., Springer (2019), New York, https://doi.org/10.1007/978-1-4939-9828-9.
  * 29. Brauer, F., Mathematical epidemiology: Past, present, and future, Infect Dis Model. 2-2, 113-127 (2017), https://doi.org/10.1016/j.idm.2017.02.001.
  * 30. Kermack, W. & Meckendrik A., A contribution to the mathematical theory of epidemics, Proc R Soc Lond. A115, 700–721 (1927), https://doi.org/10.1098/rspa.1927.0118.
  * 31. Brauer, F., The Kermack–McKendrick epidemic model revisited, Math Biosci. 198-2, 119-131 (2005), https://doi.org/10.1016/j.mbs.2005.07.006.
  * 32. Weiss, H., The SIR model and the Foundations of Public Health, Materials Matematics, 17 (2013), https://ddd.uab.cat/record/108432.
  * 33. Altmann, M., Susceptible-infected-removed epidemic models with dynamic partnerships, J. Math. Biology 33, 661–675 (1995). https://doi.org/10.1007/BF00298647.
  * 34. Fanelli, D. & Piazza, F., Analysis and forecast of COVID-19 spreading in China, Italy and France, Chaos Soliton Fract. 134, 109761 (2020), https://doi.org/10.1016/j.chaos.2020.109761.
  * 35. Baba, I.A., Yusuf, A., Nisar, K., Abdel-Aty, A. and Nofal, T., Mathematical model to assess the imposition of lockdown during COVID-19 pandemic, Results Phys. 20, 103716 (2021), https://doi.org/10.1016/j.rinp.2020.103716.
  * 36. Viguerie, A., Lorenzo, G., Aurichhio, F., Baroli, D., Hughes, T., Patton, A., et al., Simulating the spread of COVID-19 via a spatially-resolved susceptible–exposed–infected–recovered–deceased (SEIRD) model with heterogeneous diffusion, Appl Math Lett. 111, 106617 (2021), https://doi.org/10.1016/j.aml.2020.106617.
  * 37. El Deeb, O. & Jalloul M., The dynamics of COVID-19 spread: evidence from Lebanon, Math Biosci Eng. 17-5, 5618-5632 (2020), https://doi.org/10.3934/mbe.2020302.
  * 38. Yang, C. & Yang J., A mathematical model for the novel coronavirus epidemic in Wuhan, China, Math Biosci Eng. 17-3, 2708-2724 (2020), https://doi.org/10.3934/mbe.2020148.
  * 39. Silva, P., Batista, P., Lima, H., Alves, M., Guimaraes, F. and Silva, R., COVID-ABS: An agent-based model of COVID-19 epidemic to simulate health and economic effects of social distancing interventions, Chaos Soliton Fract. 139, 110088 (2020), https://doi.org/10.1016/j.chaos.2020.110088.
  * 40. Cuevas, E., An agent-based model to evaluate the COVID-19 transmission risks in facilities, Comput Biol Med. 121, 103827 (2020), https://doi.org/10.1016/j.compbiomed.2020.103827.
  * 41. Spassiani, I., Sebastiani, G. and Palu, G., Spatiotemporal Analysis of COVID-19 Incidence Data, Viruses 13-3, 463 (2021), https://doi.org/10.3390/v13030463.
  * 42. El Deeb, O., Spatial Autocorrelation and the Dynamics of the Mean Center of COVID-19 Infections in Lebanon, Front Appl Math Stat. 6, 620064 (2021), https://doi.org/10.3389/fams.2020.620064.
  * 43. Gross, B., Zheng, Z., Liu, S., Chen, X., Sela, A., Li, J., et al., Spatio-temporal propagation of COVID-19 pandemics, EPL 131, 58003 (2020), https://doi.org/10.1209/0295-5075/131/58003.
  * 44. Trigger, S. & Czerniawski E., Equation for epidemic spread with the quarantine measures: application to COVID-19, Phys Scr. 95, 105001 (2020), https://doi.org/10.1088/1402-4896/abb2e2.
  * 45. Khan, T., Ullah, R., Zaman G., and El Khatib, Y., Modeling the dynamics of the SARS-CoV-2 virus in a population with asymptomatic and symptomatic infected individuals and vaccination, Phys Scr. 96, 104009 (2021), https://doi.org/10.1088/1402-4896/ac0e00
  * 46. Chaaban, S., Doraiswamy S., Chaabna, K., Mamtani S., and Cheema S., The Impact of COVID-19 School Closure on Child and Adolescent Health: A Rapid Systematic Review, children 8(5), 415 (2021), https://doi.org/10.3390/children8050415.
  * 47. Anastassopoulou, C., Russo L., Tsakris, A., and Seittos, C., Data-based analysis, modelling and forecasting of the COVID-19 outbreak, Plos One 15(3), e0230405 (2020), https://doi.org/10.1371/journal.pone.0230405.
  * 48. Paltiel D., Zheng, A., Schwartz, J., Speed Versus Efficacy: Quantifying Potential Tradeoffs in COVID-19 Vaccine Deployment, Ann Intern Med. 174-3, 568-570 (2021), https://doi.org/10.7326/M20-7866.
  * 49. Chagla, Z., The BNT162b2 (BioNTech/Pfizer) vaccine had 95% efficacy against COVID-19 $\geq$7 days after the 2nd dose, Ann Intern Med. 174-2, JC15 (2021), https://doi.org/10.7326/ACPJ202102160-015.
  * 50. Knoll, M.D. and Wonodi, C., Oxford–AstraZeneca COVID-19 vaccine efficacy, Lancet 397, 72-74 (2021), https://doi.org/10.1016/S0140-6736(20)32623-4.
  * 51. Baraniuk, C., What do we know about China’s covid-19 vaccines?, BMJ 373, n912 (2021), https://doi.org/10.1136/bmj.n912.
  * 52. Covid: Brazil approves and rolls out AstraZeneca and Sinovac vaccines, BBC (2021), https://www.bbc.com/news/world-latin-america-55699535, (Accessed on August 1, 2021).
  * 53. Van Tulleken, C., Covid-19: Sputnik vaccine rockets, thanks to Lancet boost, BMJ 373, n1108 (2021), https://doi.org/10.1136/bmj.n1108.
  * 54. Logunov D.Y., Dolzhikova, I., Zubkova O., Tukhvatulin, A., Shcheblyakov, D., Dzharullaeva, A., et al., Safety and immunogenicity of an rAd26 and rAd5 vector-based heterologous prime-boost COVID-19 vaccine in two formulations: two open, non-randomised phase 1/2 studies from Russia, Lancet 396, 887-897 (2020), https://doi.org/10.1016/S0140-6736(20)31866-3.
  * 55. Chagla, Z., In high-risk adults, the Moderna vaccine had 94% efficacy against COVID-19 $\geq$14 d after the 2nd dose, Ann Intern Med. 174-3, JC28 (2021), https://doi.org/10.7326/ACPJ202103160-028.
  * 56. Center for Disease Control and Prevention (CDC), Interim Guidance on Ending Isolation and Precautions for Adults with COVID-19 (2021), https://www.cdc.gov/coronavirus/2019-ncov/hcp/duration-isolation.html. (Accessed on June 29, 2021).
  * 57. Lacobucci, G., Covid-19: Single dose of Pfizer and Oxford vaccines cuts risk of hospital admission by 80% in over 80s, data suggest, BMJ 372, n612 (2021), https://doi.org/10.1136/bmj.n612.
  * 58. Bernal, J.L., Andrews, N., Gower, C., Robertson, C, Stowe, J., Tessier, E., et al., Effectiveness of the Pfizer-BioNTech and Oxford-AstraZeneca vaccines on covid-19 related symptoms, hospital admissions, and mortality in older adults in England: test negative case-control study, BMJ 373, n1088 (2021), https://doi.org/10.1136/bmj.n1088.
  * 59. Moghadas, S., Vilches, T., Zhang, K., Wells, C., Shoukat, A., Singer, B., et al., The Impact of Vaccination on Coronavirus Disease 2019 (COVID-19) Outbreaks in the United States, Clin Infect Dis. ciab079 (2021), https://doi.org/10.1093/cid/ciab079.
  * 60. Sette, A. & Crotty S., Adaptive immunity to SARS-CoV-2 and COVID-19, Cell 184-4, 861-880 (2021), https://doi.org/10.1016/j.cell.2021.01.007.
  * 61. Vinceti, M., Filippini, T., Rothman, K., Di Federico, S., and Orsini, N., SARS-CoV-2 infection incidence during the first and second COVID-19 waves in Italy, Environ Res. 197, 11097 (2021), https://doi.org/10.1016/j.envres.2021.111097.
  * 62. Mazzaferro, S., Rocca, A., and Bagordo, D., Waves of infection and waves of communication: the importance of sharing in the era of Covid-19, J Nephrol. 34, 633–636 (2021), https://doi.org/10.1007/s40620-021-00974-7.
  * 63. Brooks-Pollock, E., Read, J., McLean, K., Keeling, M., and Danon, L.,, Mapping social distancing measures to the reproduction number for COVID-19, Phil. Trans. R. Soc. B. 376, 20200276 (2021), https://doi.org/10.1098/rstb.2020.0276.
  * 64. Chaturvedi, D. and Chakravarty, U., Predictive analysis of COVID-19 eradication with vaccination in India, Brazil, and U.S.A, Infect. Genet. Evol. 92, 104834 (2021), https://doi.org/10.1016/j.meegid.2021.104834. .