§ ROBOTIC BLOSSOM THINNING SYSTEM OVERVIEW

A selective blossom thinning system needs to be flexible to offer different thinning strategies depending on the flower intensity, fruit set capability, and crop-load bearing capacity of fruit tree canopies such that the thinning intensity could be tailored to achieve desired crop-load. In pome fruit crops such as apples and pears, flowers are closely spaced in clusters and are often located in constrained working spaces caused by canopy parts such as trunks and branches, and trellis structures such as trellis wires and posts. Once flowers are in full bloom, dense clusters with high flower-to-flower occlusions and larger cluster sizes make individual flower delineation and removal challenging.

While it is comparatively easier for humans to identify and perform dual arm motion to remove unopened flowers, detection of individual flowers, maneuvering of the robotic arm, estimating precise position and orientation of individual flowers, and positioning of an end-effector is challenging, slow, and computationally intensive. Hence instead of locating and removing individual flowers, the proposed robotic system is developed to remove a proportion of flowers from a cluster of flowers during full bloom. This approach reduces the operation cycle time since the robotic system will have a single target cluster instead of multiple flowers within a cluster. Furthermore, the proportional thinning approach reduces the chances of the robot being in singularity and reduces the computational complexity for the manipulator path planning and obstacle avoidance algorithms. While mechanical blossom thinning systems have demonstrated potential for variable blossom thinning intensities, the major concern from the growers was the damage mechanical thinners could cause to the canopies due to the end effector size. To address the size issue, a miniature electrically actuated end-effector operable within the flower clusters was designed. We propose that the use of robotic system equipped with miniature mechanical thinners with small operational space could be beneficial to the growers, allowing for removal of proportion of flowers en masee while substantially reducing the damage to the spurs and leaf tissue.

Overview of the proposed robotic blossom thinning system. The system consisted of three major components a machine vision system, a manipulator and motion planning system, and an end-effector system. Mask R-CNN based deep learning algorithm was investigated to identify and delineate cluster boundaries followed by cluster pose (position and orientation) estimation. The cluster visit sequence was determined and motion plan for the robotic manipulator was developed to navigate end effector to desired thinning position and orientation to conduct thinning.

Figure <ref> illustrates the overview of the proposed robotic blossom thinning system consisting of three major subsystems: machine vision system, robotic manipulator and motion planning system, and an end-effector system. The machine vision system included associated sensors and an image processing pipeline for sensing, identifying, and delineating flower clusters. The identified flower clusters were examined for safe reachability using automated and manual (human in loop) approach. Clusters beyond 1m depth, and the clusters with invalid depth measurements due to camera sensor error were automatically excluded. The remaining clusters were further manually evaluated to ensure the safe reachability of the manipulator for any collisions with other orchard components, including the trunk, branches, trellis wires, and training posts.
The flower clusters in orchards were randomly distributed with random position and orientation. Hence, the cluster pose (position and orientation) was estimated such that the end effector could be positioned depending upon the cluster orientation to achieve desire thinning result. In this work the end effector was positioned orthogonal to the cluster surface during the thinning operation. Based on the cluster pose information from the machine vision system, visit to the individual clusters were prioritized by computing shortest path to be followed by the manipulator such that each clusters were visited once.  Once the cluster visit was prioritized, the corresponding motion plan for six degrees of freedom robotic manipulator was computed and executed to navigate the end effector to the desired position and orientation to perform thinning. The end-effector system actuated the end-effector until all the thinning waypoints were visited depending upon the thinning strategy. The thinning waypoints were defined depending on the employed thinning strategies. Two thinning strategies (Boundary Thinning and Center Thinning) were studied to demonstrate the ability to control thinning intensity to perform selective thinning from target clusters. For boundary thinning the end effector was navigated and actuated around the cluster boundaries while for center thinning, the end effector was navigated to the cluster centroid and actuated for fixed 2 seconds. The following sections discuss the integrated system’s experimental environment and various hardware and software subsystems in more detail.

§ EXPERIMENTAL ENVIRONMENT

A commercial apple orchard growing Envy (Scilate cultivar) apple variety in Prosser, WA was selected for field test of the proposed robotic flower thinning system. The apple trees were grown in a modern trellis-trained high-density planar architecture as shown in Figure<ref>. The trees were trained to a V-trellised system leading to an inverted pyramidal structure of the canopies. Traditionally, apple trees are planted in rows with wide spacing and trunks and branches are allowed to grow in all directions, which leads to a three-dimensional canopy structure resulting in 4D working space, making it highly challenging for any robotic operations. The modern orchard systems use trellis wires to support and train fruit trees to create a narrow vertical plane of canopy with flower and fruit growing mostly in the outer surface (also called a fruiting wall architecture). This simplified working space helps improve accessibility, efficiency and productivity of both human workers as well as mechanical and robotic machines.

The experimental orchard (Envy apple variety) used for evaluating the robotic blossom thinning system. The trees in this experimental site were trellis-trained to create a V-shaped fruiting wall orchard.

Attributes of the V-trellised orchard used for field evaluation of the robotic thinning system (from <cit.>)

Variety    Canopy Angle    No. of Trellis    Trellis Wire Spacing    Tree Height    Tree Spacing    Row Width

(Degree)    Wires    (cm)    (cm)    (cm)    (cm)

Envy    75    7    46$\pm$1    366    142$\pm$19    353$\pm$3

§ MACHINE VISION SYSTEM

A machine vision system is critical for any robotic system to operate. The vision system provides a capability to robots to locate objects it would like to manipulate (apple flowers in this case), the objects that need be avoided (e.g., trunk, branch, trellis wire, support posts), and the other components (leaves) in the robot’s workspace. As a first step in developing an autonomous robot for blossom thinning, this work focused on developing deep learning-based models to detect, segment, and localize apple flower clusters in images collected from outdoor orchard environment.

§.§ Eye-in-hand camera system

The machine vision system developed for this work used an Intel RealSense D435i RGB-D camera as eye in hand sensor that was rigidly attached to the tool flange of the robotic manipulator. Such attachment allowed fixed transformation between the camera to the robot coordinate space while providing multi-view imaging of objects when the robot is in motion. The camera system used a single RGB imaging sensor capable of providing a maximum resolution of $1920 \times 1080$ pixels for the given Field of View (FOV) of  $69\degree  \times 42\degree$. Depth was estimated using active stereo IR technology with a maximum resolution of  $1280 \times 720$ with FOV of $87\degree \times 58\degree$. In this work, the RGB and depth image resolution of $1280 \times 720$ pixels was used to speed up the object segmentation pipeline while covering the region within the reachability of the manipulator.

§.§ Flower cluster detection and segmentation

Since the integrated system requires the manipulator to navigate the end-effector that can engage with the individual clusters in a unique way, each segmented cluster requires a unique identifier. This work used the Mask R-CNN model <cit.> for flower cluster detection and pixel-level instance segmentation. The instance segmentation technique was selected over semantic segmentation since instance segmentation separates the multiple instances of the same object class with unique identifiers. Furthermore, the end-effector had to interact within the cluster more precisely (e.g., following the flower cluster boundary); a pixel-level segmentation algorithm was well suited instead of bounding box-based object detection methods such as Faster R-CNN and You Only Look Once (YOLO).

In our previous work, the Mask R-CNN model showed promising accuracy and robustness in detecting and segmenting flower clusters in outdoor orchard images <cit.>.  However, when the same model was applied to the images acquired from the eye-in-hand camera attached to the robotic manipulator, the model struggled to identify clusters. One of the major reasons for the decreased cluster detection was due to the closer imaging distance $(\sim50cm - \sim 75cm)$ for the eye-in-hand camera, which differed from the $\sim 1m$ distance used in the previous work <cit.>. This was verified by manually moving the camera away from the canopy which resulted in improved cluster detection. Hence to train the Mask R-CNN model for close range images, new image dataset was collected from a commercial apple orchard in Prosser, WA, using Intel RealSense D435i RGB-D camera at the resolution of $1280 \times 720$ for both color and depth images. Tree canopies were continuously scanned at a distance of $\sim 70 cm$ (similar to distance of the eye in hand camera) from the trunk in a daylight, and the images frames were extracted later from the recorded data stream. The flower clusters in these images were annotated as polygons using the Labelbox image annotation tool (https://labelbox.com/). The image dataset consisted of 636 images with a total of 2,118 flower clusters, which were randomly divided into training, validation, and test sets with a split ratio of $70\% : 15\% : 15\%$ images.

A Mask R-CNN model was then trained extending the existing Mask R-CNN implementation from Matterport Inc.(Sunnyvale, CA) released under MIT License <cit.>. The system employed TensorFlow 2.0 with Keras as backend engine and ResNet101 as the backbone. The ResNet101 backbone was pre-trained with the
MSCOCO dataset. The network was trained up to 100 epochs using Stochastic Gradient Descent (SGD) search technique with a momentum of 0.9 and a Learning Rate (LR) of 0.00005. Different image augmentation techniques, including random horizontal and vertical flips, random scaling ($0.5$ – $1.5$), and random rotation ($-60 \degree$ – $+60\degree$), were used to increase the dataset variability.

The segmented cluster boundary consisted of a large number of boundary points as the Mask R-CNN model was trying to closely and smoothly locate the cluster boundary. However, to perform robotic thinning, it was not necessary and practical for the end-effector to navigate to all the vertices of the polygon representing the flower cluster boundaries. First, navigating to all the vantage points (vertices) would result in additional computational overhead, which is complex and time-consuming for the motion planner. Second, in most cases, the end-effector was large enough to cover multiple vertices on the boundary at once, making individual visits to vertices redundant. Hence, a Convex Hull algorithm was used to simplify the cluster boundary polygons. The objective of the Convex Hull algorithm was to estimate the smallest convex polygon such that all the cluster boundary points were either on or within the polygon representing the cluster boundary. The implemented Convex Hull algorithm was available in the OpenCV library and had $O(NLogN)$ complexity.

§.§.§ Flower Cluster Reachability

Once the potential flower clusters were segmented, individual clusters were evaluated for potential to reach them safely (without any robot-canopy collision). The flower clusters that did not satisfy the robot reachability requirements were filtered out from the target clusters for thinning. First, clusters beyond 1m depth and having invalid depth due to camera sensor error were first filtered out automatically. The remaining clusters were further evaluated for reachability by the system’s operator as a part of the human-in-the-loop system. The objective of the human in the loop system was to examine the detected clusters from the perspective of proximity to the trunks, branches, and trellis wires, and remove clusters from the cluster target list if navigation to, and thinning of those clusters could cause potential collision between the robotic system and the canopy parts. The effort of the current human-in-loop system could be reduced by detecting and segmenting other canopy parts such as trellis wires, trunks, branches, leaves, and support posts, and making automated thinning decisions based on the operation environment. The current deep learning-based machine vision system developed for cluster segmentation is flexible and can be extended to identify other canopy parts with few modifications.

§.§ Cluster Localization

Flower clusters grow in random/variable locations and orientation within constrained canopy spaces. It is critical to determine the precise position and orientation of the clusters in 3D space such that the thinning end-effector can be positioned and oriented precisely for effective thinning. In this work, the cluster position was estimated by projecting the detected cluster pixels in 3D space and averaging the corresponding 3D cluster position. The cluster orientation was then computed by estimating the segmented cluster’s surface normal vector as illustrated in Algorithm <ref>. Let $P_1 \dots P_N$ be the 3D position of all points in the flower cluster with camera viewpoint location $(v_p)$. To estimate the surface normal of the entire cluster, the cluster was divided into sub-sections defined by radius $(R)$ and number of neighborhood points $(k)$. The surface normal of each cluster subsection was estimated and then averaged to compute the surface normal of the entire cluster (see Algorithm <ref>).

Flower cluster orientation estimation
[1]

$P_1 \dots P_N$, camera viewpoint $(v_p)$
Cluster position ($\textbf{P}$), Cluster normal ($\vec{\textbf{n}}$)
$\textbf{P}=\frac{1}{N}\sum_{i=1}^N {P_i}$

$R$ $\gets$ $0.1$
$k$ $\gets$ $30$
Normal vector list $n_l$ $\gets$ $[\;]$

$i \gets 1$ to $N$
Build a $k-D$ tree
Conduct radius $(R)$ neighbourhood search for $k$ neighbours
Estimate covariance matrix $C$

\begin{equation}
C=\frac{1}{k}\sum_{l=1}^k {(p_l-\overline{p})}\cdot {(p_l-\overline{p})^T},\overline{p}=\frac{1}{k}\sum_{l=1}^k{p_l}, p_l\in P^k
\end{equation}

Estimate eigenvalues of $C$:  $C\cdot \vec {v_j}=\lambda\cdot \vec {v_j}, \textit{j}\in {\{0,1,2\}}$, where, $0\leq\lambda_0\leq\lambda_1\leq\lambda_2$

Surface normal $\vec{n_i}={\{n_x,n_y,n_z\}}$; vector $\Vec{v_0}$

Orient $\vec{n_i} \ni \vec{n_i} \cdot (v_p-p_i)>0$
Append $\vec{n_i}$ to $n_l$
Cluster normal $\vec{\textbf{n}}=\frac{1}{m}\sum_{i=1}^m{\vec{n_i}}, \vec{n_i}\in n_l$

$\textbf{P}, \vec{\textbf{n}}$

§ MECHATRONIC SYSTEM

§.§ Robotic Manipulator and Motion Planning

For the given position and orientation of the flower cluster, the objective of the manipulator was to navigate the end-effector to the thinning location orthogonal to the cluster surface to conduct blossom thinning. An UR5e collaborative robotic manipulator manufactured by Universal Robots (Odense, Denmark) was used in this study. The manipulator has a small footprint, a 5KG payload that was sufficient for orchard operations such as thinning, pruning, and harvesting, and safety functions related to joint speed, torque limit, Tool Center Point (TCP) pose, speed, and force limit, among others. Table <ref> details the robotic manipulator specifications. The experimental orchard used was trained to 2D planar structures making a large proportion of flowers accessible for the robotic thinning. The robot motion planning was achieved using the default ROS MoveIt Open Motion Planning Library (OMPL).

UR5e robot technical specifications <cit.>

DOF    Payload (kg)    Reach (mm)    Footprint $\phi$ (mm)    Weight (kg)    Power Consumption (W)

6     5    850    149    20.5    200

§.§ End-effector design

The major requirement of the blossom thinning end-effector was to exert a high impact on the target flowers to be removed while minimally affecting the flowers to be saved for pollination as well as neighboring leaves and spurs. Flowers are tightly located in clusters, often overlapping each other. Hence, the end-effector should also be small and precise, impacting only at the desired location in flower clusters. As a preliminary experiment, three different actuation mechanisms were considered: water-actuated, air-actuated, and electrically-actuated end-effectors. For both the water and air actuated system, water and air were pressurized in the high-pressure tank and were released as a water jet/air stream using a sub-millimeter size nozzle. Although both water jet and air stream approaches removed some of the flowers, two major problems were observed with these thinning approaches; i) While removing the unnecessary flowers, the flowers that should have been saved were frequently impacted/damaged by the incoming water jet and air stream often dragging them with the water/air flow; ii) The incoming water jet or air stream also induced damage to the flowers in proximal clusters and the clusters in the background. After the preliminary experiments, waterjet and air stream techniques were, therefore, dropped and only the electrically actuated brushing mechanism was further investigated.

Electrically actuated end-effector design with spindle string. Strings used were the pieces of off-the-selves grass trimming lines with grooved edges.

The electrically actuated end-effector system consisted of a miniature spindle string configuration. Commercial-grade string trimmer lines were tied together in an aluminum spindle at $\sim 7mm$ spacing, as shown in Figure <ref>. The strings have grooved edges for faster trimming. This end-effector was coupled with a high-speed motor with a rigid shaft coupler to provide the necessary force to thin flowers quickly and efficiently. During the preliminary experiment, the end-effector was held manually to mimic robotic thinning motion—a clockwise sweeping motion around the cluster boundary with the approach direction orthogonal to the cluster surface. Three different experiments were conducted by varying the motor speed: i) Low speed (3000rpm-4000rpm), ii) Medium speed (5000rpm-6000rpm), and iii) High speed (7000rpm-8000rpm). For each thinning speed, 30 flower clusters were thinned to evaluate thinning efficiency. High-speed motor actuation showed better thinning efficiency compared to others with faster trimming of excess flowers while impacting the neighboring flowers in clusters minimally. Furthermore, it was observed that the trimmer line strings were rigid enough to remove flowers, and minimally damaged branches and spurs in small areas in the event they were present.

§ EXPERIMENTAL SETUP FOR FIELD EVALUATION

Integrated system during the field evaluation in V-trellised fruiting wall architecture commercial apple orchard in Prosser, WA. The UR5e robotic manipulator was attached rigidly at a back of a utility vehicle.

Figure <ref> shows the experimental setup during the field test in a commercial apple orchard in Prosser, WA. All the hardware components were mounted on an electric utility vehicle (John Deere, IL, USA). The mobile platform was driven along the apple tree rows during the test. The experiments were conducted in stop, execute, and go fashion, where the mobile platform was stationary when the manipulator and end-effector were actuated for blossom thinning. The weather condition during the field test was sunny to cloudy, with temperatures ranging from $10-25\degree C$. The wind was variable, ranging from minimal (5mph) to high (18mph).

Figure <ref> shows the hardware components of the robotic system, which primarily included a 6DOF – UR5e manipulator, UR5e controller, UR5e teach pendant, laptop, Intel RealSense 435i RGB-D camera, H-bridge motor driver and Arduino Uno, DC 777 motor, and an end-effector. The manipulator and the laptop were electrically powered by Honda EU2200i (2200 W, 120 V AC) portable generator (Honda Motor Co. Ltd, Tokyo, Japan), while a 12V DC battery powered the end-effector motor. The laptop used was Dell Alienware15R4, with Intel Core i7-8750 CPU, 2.2 GHz processor, 64-bit operating system, 32GB RAM, and 8GB NVIDIA GeForce GTX 1080 Graphics Processing Unit (GPU) running Ubuntu 18.04 LTS operating system. The laptop interacted with the manipulator via Local Area Network (LAN) connection via the UR5e controller box. Communication between the laptop and Intel RealSense camera and Arduino microcontroller was established via a USB 3.0 connection. Since Arduino could only supply up to a 5V signal, the H-bridge motor driver was used as a bridge to control the 12V supply from the battery based on the signal from Arduino.

Hardware components used during field test and communication protocol. Dell Alienware15R4 laptop was used as the central processing system of the integrated system.  The laptop communicated with the Intel RealSense camera and Arduino via USB connection while LAN connection was established between the UR5e controller and the laptop to control the UR5e manipulator.

§.§ Thinning approaches

As discussed in Section <ref>, selective thinning requires the ability to customize the thinning intensity based on different factors such as flower density, fruit set capability, cultivar, and size and distribution of branches and trunks of target fruit trees. To demonstrate the system’s ability to perform variable thinning, two thinning strategies were investigated for flower removal.

* Boundary thinning: In this approach, size of individual clusters was estimated, and the end-effector was actuated around the cluster boundary in a clockwise sweeping motion.

* Center Thinning: In center thinning approach, the end-effecter was moved to the cluster centroid and rotated for a specific duration (2 sec) irrespective of the cluster size and without sweeping motion used in the boundary thinning.

The end effector approach direction for both boundary and center thinning was orthogonal to the cluster surface. In this field evaluation, a total of 148 flower clusters were attempted for thinning (Boundary Thinning:72, Center Thinning: 76). Out of 148 samples, three samples were discarded since the robot started unpredictable or erratic behavior. Since the robotic manipulator has a reachability of 850 mm, the thinning experiments were performed in the clusters located on the branches trained to three trellis wires (second, third, and fourth from the bottom) of the 7-wire apple tree canopies that ensured good reachability of the manipulator  (see Figure <ref>). Figure <ref> shows the images captured during the thinning process covering different sections of the tree canopy.

V-trellised commercial apple orchard in Prosser, WA used for field evaluation of the robotic system. Red lines highlight the section of the canopy where the integrated system was tested.

Examples images depicting the integrated system in action during field evaluation in a commercial apple orchard around three trellis wires: Lower trellis wire (Left), Center trellis wire (Middle), Upper trellis wire (Right)

§ SOFTWARE ARCHITECTURE

Figure <ref> details the software architecture of the integrated system used in this study. All the software components were encapsulated and executed using a Dell Alienware15R4 laptop within ROS Melodic. However, two distinct ROS workspaces and computational environments were set up: Default, Virtual. The Virtual environment was responsible for executing the majority of 2D and 3D machine vision algorithms and operations, such as 2D cluster detection and segmentation and 3D pose estimation. The Default environment was the core of the system, which communicated and managed the peripheral devices such as the IntelRealSense camera, UR5e controller and manipulator, Arduino end-effector controller, and the Virtual environment. Communication among different components in the software layer was managed using the ROS ecosystem.

At the beginning of an experiment, the IntelRealSense D435i camera acquired an RGBD image of the canopy. As discussed in section <ref>, the 2D RGB image was used as the input to the instance segmentation algorithm, which detected and estimated the cluster boundary. The estimated cluster boundary was simplified using the Convex Hull algorithm.

Flowchart of the software architecture used in this study. The diagram illustrates information flow among different components of the robotic thinning system.

Given a set of 3D positions of the feasible/accessible clusters, and the starting position of the manipulator, the next step was to determine the shortest possible route such that the end-effector visits each cluster exactly once and returns to the start position with minimal travel cost (Euclidean distance). This problem was modeled as the classic optimization problem known as the Traveling Salesman Problem (TSP), which is an “NP-hard” problem because of the difficulty in finding solutions for large instances. The TSP was solved by employing Local Search algorithm with OR-Tools, an open-source software to solve optimization problems <cit.>). The Local Search algorithm started with a random end-effector path that visits each cluster and iteratively improves the path by swapping the locations of a pair of clusters and checking if the new path was better until no further improvements could be made. Based on the solution from solving the TSP, the cluster visit sequence for the end effector was prioritized, and the thinning was subsequently performed accordingly.

The thinning process of clusters were modeled as a state machine, as illustrated in Figure <ref>.  During each thinning cycle, the system was in one of the four states: Home, Approach, Thin, or Retract. At the beginning and the end of the experiment, the system was in an idle Home state, waiting for commands. During the initial Home state of the robot, the system carried out various tasks including image acquisition, cluster detection and segmentation, pose estimation, as well as post-processing procedures such as filtering out unreachable clusters, implementing Convex Hull-based boundary simplification, and cluster visit prioritization.

Let $[T_0,T_1,T_2.........T_N]$ be $N$ target clusters selected and prioritized for thinning. During the thinning operation, the end-effector was navigated to three vantage points: Approach Point, Thinning Start Point, and Retract Point.

The system states Approach, Thin, and Retract were defined for motion planning and navigation of the end-effector to these three vantage points. The Approach Point was chosen at a 10 cm offset distance from the cluster centroid along the cluster surface normal. The Approach Point was introduced as an intermediate target such that the end-effector was aligned orthogonally to the cluster surface and ready to approach the cluster. In this study, the canopy was assumed a 2D planar structure. Orienting the end-effector along the cluster surface normal at the Approach Point before physically interacting with the cluster required only a simplified linear motion from the Approach Point to the Thinning Start Point. The Thinning Start Point was the first point within the cluster where the end-effector physically interacted with the flowers in the target cluster. For center thinning, the Thinning Start Point was the cluster centroid, while for Boundary Thinning, the Thinning Start Point was the first vertex in the polygonal boundary.

Simplified illustration of the thinning process as a state machine. This state machine was executed after the target flower clusters were finalized and prioritized.

The thinning process started after the end-effector reached the Thinning Start Point, and the system entered to the Thin state. The system would be in the Thin state for a 2-second for center thinning or until all the boundary points were visited for boundary thinning. The Retract Point was same as the Approach Point, except it was visited after the thinning process was completed. Approach Point and Retract Point ensured the system’s safety by forcing majority of the motion planning tasks to occur outside of the canopy, thereby reducing the risk of collision with canopy parts such as branches, trunks, and trellis wire. Once thinning process was successfully completed or failed for a target cluster, the cluster was removed from the Target List, and the process was repeated until all clusters were visited. The system returned to the Home state at the end of the thinning cycle.


# Modeling the effect of the vaccination campaign on the Covid-19 pandemic

Mattia Angeli<EMAIL_ADDRESS>Georgios Neofotistos  Marios
Mattheakis  Efthimios Kaxiras  John A. Paulson School of Engineering and
Applied Sciences, Harvard University, Cambridge, Massachusetts 02138, USA
Department of Physics, Harvard University, Cambridge, Massachusetts 02138, USA

###### Abstract

Population-wide vaccination is critical for containing the SARS-CoV-2
(Covid-19) pandemic when combined with restrictive and prevention measures. In
this study we introduce SAIVR, a mathematical model able to forecast the
Covid-19 epidemic evolution during the vaccination campaign. SAIVR extends the
widely used Susceptible-Infectious-Removed (SIR) model by considering the
Asymptomatic (A) and Vaccinated (V) compartments. The model contains several
parameters and initial conditions that are estimated by employing a semi-
supervised machine learning procedure. After training an unsupervised neural
network to solve the SAIVR differential equations, a supervised framework then
estimates the optimal conditions and parameters that best fit recent
infectious curves of 27 countries. Instructed by these results, we performed
an extensive study on the temporal evolution of the pandemic under varying
values of roll-out daily rates, vaccine efficacy, and a broad range of
societal vaccine hesitancy/denial levels. The concept of herd immunity is
questioned by studying future scenarios which involve different vaccination
efforts and more infectious Covid-19 variants.

###### keywords:

Covid-19 $|$ Machine Learning $|$ Neural Networks $|$ Vaccines

††journal: Chaos, Solitons & Fractals

##  1 Introduction

The World Health Organization (WHO) declared the SARS-CoV-2 (Covid-19)
outbreak in Wuhan to be a pandemic on March 11, 2020. Since then, Covid-19 has
become a serious global health threat due to its rapid spread, transmission
through asymptomatic infected individuals and complex epidemiological
dynamics. As of May 2021, already more than 3 million lives have been lost due
to the virus. The spread of SARS-CoV-2 has thus far been extremely difficult
to contain.
By the end of 2020, the successful development of effective vaccines and the
onset of their widespread distribution in most of the world’s countries, was
hailed as the decisive mean to contain the pandemic. However, important
questions linger on whether the vaccination effort will succeed in effectively
eradicating the disease. The appearance and wide spread of more contagious
SARS-Cov-2 strains, the onset and scale of the vaccine deployment and high
levels of vaccine hesitancy/denial in the society, are among the key factors
hindering the vaccination effort and the achievement of herd immunity.
Modeling the impact of these key factors on the evolution of the pandemic is
of critical importance for assessing the vaccination effectiveness against it.
In studying past epidemics, scientists have systematically applied ”random
mixing” compartmental models which assume that an infectious individual can
spread the disease to any susceptible member of the population before becoming
recovered or removed, as originally considered by Kermack and McKendrick [1].
These models constrain the total population in compartments by considering
stages of the infection and flows among them.

In the present study we propose a new model named SAIVR, which incorporates
two important characteristics of the Covid-19 epidemic, namely the
considerable transmission of the disease by asymptomatic infected individuals
and the vaccination campaign with World Health Organization (WHO) approved
vaccines. More recent modeling approaches involve agent-based simulations [2],
heterogeneous social networks [3, 4, 5, 6, 7, 8], and Bayesian inference
models [9]. Although a large number of research studies are currently
investigating the Covid-19 epidemiological characteristics [10, 11, 12, 13,
14, 15, 16, 17, 18], we believe that a simple but efficient model, which can
capture the basics of the complex behavior of the pandemic including the
vaccine roll-out, can offer useful guidance for the pandemic’s near-term and
longer-term evolution. By using a recently developed semi-supervised machine
learning approach [19, 20, 21] we systematically reproduced the pandemic
dynamics during the 2021 spring in several different countries. We then used
the model to assess the importance of a rapid vaccination campaign to prevent
future outbreaks driven by more infectious variants.

The work is organized as follows. In Sec. 2 we introduce the SAIVR model and
its parameters. The machine learning approach that we used to reproduce the
infectious curves of 27 selected countries/states is thoroughly described in
Sec. 3. In Sec. 5 we study future scenarios involving more infectious variants
making quantitative arguments on how they might affect herd immunity. Sec. 5
is devoted to concluding remarks.

##  2 The SAIVR model

One of the first attempts to mathematically describe the spread of an
infectious disease is due to Kermack and McKendrick [1]. In 1927 they
introduced the so-called Susceptible-Infectious-Removed (SIR) model. The SIR
model describes the dynamics of a (fixed) population of $N$ individuals split
into three compartments:

  * 1.

$S(t)$ is the Susceptible compartment that counts the number of individuals
susceptible but still not infected by the disease;

  * 2.

$I(t)$ is the Infectious compartment that counts the number of infectious
individuals;

  * 3.

$R(t)$ is the Removed compartment. It represents the number of those who can
no longer be infected either because they recovered and gained long-term
immunity or because they passed away.

Figure 1:  Illustration of the SAIVR model compartments and their inter-
dependencies denoted by incoming and outgoing arrows and relevant flow
parameters.

The model involves two positive parameters, $\beta$ and $\gamma$ which govern
the flow from one compartment to the other:
\- $\beta$ is the transmission rate or effective contact rate of the disease:
an infected individual comes into contact with $\beta$ other individuals per
unit time (the fraction that are susceptible to contracting the disease is
$S/N$);
\- $\gamma$ is the removal rate. $\gamma^{-1}$ is the mean number of days who
is infected spends in the Infectious compartment.
The SIR model obeys the following system of ordinary differential equations
(ODE):

$\displaystyle\frac{dI}{dt}$ $\displaystyle=\beta I\frac{S}{N}-\gamma I$ (1a)
$\displaystyle\frac{dS}{dt}$ $\displaystyle=-\beta I\frac{S}{N}$ (1b)
$\displaystyle\frac{dR}{dt}$ $\displaystyle=\gamma I$ (1c)

Although the SIR model has been adopted to study epidemic outbreaks in many
previous works [22, 23, 24, 25, 26, 27, 2], it lacks few important aspects of
the current ongoing pandemic. First of all, it has been reported [28, 29] that
an important fraction of those who are carrying the virus is asymptomatic.
Since they often avoid contact tracing due to the absence of symptoms, they
can spread the disease while remaining undetected. Furthermore, in December
2020 a global vaccination campaign has started. Vaccinating is a safe way to
transfer people from the Susceptible to the Removed compartment bypassing the
Infectious one thus reducing the likelihood of an outbreak.

The SAIVR model extends the SIR model by incorporating the two aforementioned
additional compartments:

  * 1.

$A(t)$ is the Asymptomatic/Undetected compartment that counts the number of
those individuals that despite being infected are not tested/traced. This
mainly occurs to those who recover from the infection without suffering any
symptoms.

  * 2.

$V(t)$ is the Vaccinated compartment. It takes into account those that have
received a vaccine shot but are still not fully immunized by it.

The SAIVR model ODEs read:

$\displaystyle\frac{dI}{dt}$
$\displaystyle=\beta_{1}I\frac{S}{N}+\alpha_{2}A\frac{S}{N}+\zeta
I\frac{V}{N}-\gamma I$ (2a) $\displaystyle\frac{dA}{dt}$
$\displaystyle=\alpha_{1}A\frac{S}{N}+\beta_{2}I\frac{S}{N}+\eta
A\frac{V}{N}-\gamma A$ (2b) $\displaystyle\frac{dS}{dt}$ $\displaystyle=-\beta
I\frac{S}{N}-\alpha A\frac{S}{N}-\delta\frac{S}{N}+(1-\lambda)\epsilon V$ (2c)
$\displaystyle\frac{dV}{dt}$ $\displaystyle=\delta\frac{S}{N}-\eta
A\frac{V}{N}-\zeta I\frac{V}{N}-\epsilon V$ (2d) $\displaystyle\frac{dR}{dt}$
$\displaystyle=\gamma I+\gamma A+\lambda\epsilon V$ (2e)

The compartment inter-dependencies and flow are presented in Fig. 1. The
parameters of the SAIVR model are the following:
\- $\beta_{1}$ describes the rate at which individuals are exposed to
symptomatic infection. An infected symptomatic individual comes into contact
and infects $\beta_{1}$ susceptible individuals per unit time;
\- $\alpha_{1}$ is the asymptomatic infection rate. An infected asymptomatic
individual comes into contact with $\alpha_{1}$ susceptible individuals per
unit time;
\- $\beta_{2}$ describes the rate at which susceptible individuals become
asymptomatic infected after entering in contact with a symptomatic individual;
\- $\alpha_{2}$ describes the rate at which who’s susceptible becomes
symptomatic after entering in contact with an asymptomatic individual;
\- $\gamma$ retains the same meaning as in the SIR model, representing the
mean removal rate. $\gamma^{-1}$ is the mean amount of time individuals spend
either in the Infectious or Asymptomatic compartments;
\- $\zeta$ is the rate at which a vaccinated (but still not immune) individual
enters in contact with a symptomatic infectious;
\- $\eta$ describes the transmission rate at which who’s asymptomatic comes
into contact and infects vaccinated (but still not immune) individuals;
\- $\delta$ is the first shot vaccination rate;
\- $\lambda$ is the vaccine efficacy;
\- $\epsilon^{-1}$ is the mean amount of time an individual spends in the
Vaccinated compartment before reaching immunity and moving to the Removed
compartment.

Countries and states do not respond to the disease as static entities
passively facing the pandemic. They react by actively imposing (and relaxing)
restrictive measures, learning how to effectively treat the infected,
adjusting social interactions and by launching vaccination campaigns. Finally,
the virus itself evolves in more infectious variants [30].
Country-specific parameters can be obtained by fitting the SAIVR model to a
selected infectious wave occurred in a given country. SAIVR has 14 adjustable
parameters or initial conditions that needs to be estimated; given the
scarcity of data (only the infectious and vaccinated populations are known)
optimizing them presents a challenging problem. To address this, we either
fixed some of them or employed a novel fitting method based on semi-supervised
neural networks, which we present in the following section.

##  3 Solving the SAIVR model with Neural Networks

Figure 2:  Semi-supervised network architecture. During the unsupervised
procedure (blue box), a time sequence $t$, a set of initial conditions
$\mathcal{Z}_{0}$ and parameter bundles $\Theta$ are fed as an input to a 6
layers fully connected network (FCN). The output of the network
$\mathcal{Z}_{NN}$ is multiplied by a function $f(t)$ to become a tentative
parametric solution $\hat{\mathcal{Z}}(t)$ of the system of ODE in Eq. 2. The
quality of $\hat{\mathcal{Z}}(t)$ is probed by the loss function
$\mathcal{L}$. When the network has learned the solutions, the inverse problem
is then solved (red box). An optimization algorithm selects the initial
conditions and parameters in the bundle that best fit a given data-set. The
loss $\mathcal{L}_{inv}$ depends on the infectious population of a given
country/state $I_{Data}$ during the time sequence $t$.

In order to apply the SAIVR model we need a realistic estimate of the
parameters and initial conditions for the system of Eqs. 2. To obtain them we
employed machine learning, a powerful method which has been extensively used
for disease modeling [31, 32, 33, 21] and dynamical system forecasting [34,
35, 19]. Our approach employs a semi-supervised procedure which determines the
optimal set of initial conditions and parameters of the SAIVR model, yielding
solutions that best fit a given data-set. A sketch of this procedure is shown
in Fig. 2.

###  3.1 Unsupervised learning

The unsupervised part (blue box) consists of a data-free Neural Network (NN)
that is trained to discover solutions for an ODE system of the form:

$\frac{d\mathcal{Z}}{dt}=g(\mathcal{Z}),\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mathcal{Z}(t=0)=\mathcal{Z}_{0}$
(3)

where $\mathcal{Z}=(S(t),A(t),I(t),V(t),R(t))$ and $g(\mathcal{Z})$ is given
in Eqs.2. The NN takes as an input a time sequence $t$, a set of initial
conditions $\mathcal{Z}_{0}$, and modeling parameters $\Theta$.
As we’ll see in the following, $t$ is the set of days involved in a given
epidemic wave going from $t_{0}$ to $t_{0}+\Delta t$, $Z_{0}$ are the initial
compartment populations and $\Theta$ some parameters of the SAIVR model. The
initial conditions and parameters are randomly sampled at each iteration $n$
over predefined intervals called bundles [19, 21], so that the network learns
an entire family of solutions. The inputs propagate through the network until
an output vector $\mathcal{Z}_{NN}$ of the same dimensions as the target
solutions $\mathcal{Z}$ is produced. The learned solutions $\hat{\mathcal{Z}}$
satisfy the initial conditions identically by considering parametric solutions
of the form:

$\hat{\mathcal{Z}}=\mathcal{Z}_{0}+f(t)(\mathcal{Z}_{NN}-\mathcal{Z}_{0})$ (4)

where $f(t)=1-e^{-t}$ [20]. The loss function:

$\mathcal{L}=\bigg{\langle}\left(\frac{d\mathcal{Z}}{dt}-g(\mathcal{Z})\right)^{2}\bigg{\rangle}$
(5)

solely depends on the network predictions averaged ($\langle..\rangle$) over
all the iterations $n$, providing an unsupervised learning framework. Time
derivatives are computed using the automatic-differentiation and back-
propagation techniques [36].

###  3.2 Fitting a dataset

Once the NN is trained to provide solutions for the system of Eq. 2, its
weights and biases are fixed and the trained network is used to develop a
supervised pipeline for the estimation of the initial conditions and
parameters, leading to solutions $\tilde{\mathcal{Z}}(t)$ that fit given data.
This procedure is illustrated in the red box in Fig. 2. A solution
$\hat{\mathcal{Z}}$ is generated by the network starting from
$\bar{\mathcal{Z}_{0}}$ and $\bar{\Theta}$ randomly selected in the bundles. A
stochastic gradient descent optimizer then adjusts $\bar{\mathcal{Z}_{0}}$ and
$\bar{\Theta}$ in order to minimize the loss function:

$\mathcal{L}_{inv}=\bigg{\langle}\left(\tilde{I}(t)-I_{Data}(t)\right)^{2}\bigg{\rangle}$
(6)

where $I_{Data}(t)$ is the infectious population of a given country/state and
$\tilde{I}(t)$ is its NN fit.

###  3.3 Applying the method to real data

We used this method to reproduce the most recent Covid-19 waves in 27
countries or states. To test the generality of the model we selected epidemic
waves that occurred in a broad range of geopolitical conditions, restrictive
measures, time periods and vaccination efforts. The bundles and fixed
parameters used during the unsupervised training of the network are listed in
Table. 1.

We found that the model is weakly sensitive on the choice of most parameters
and thus, we kept some of them fixed during the training. The value of the
vaccine efficacy $\lambda$ is based on Refs.[37, 38], where a vaccine efficacy
of $94.8\%$ and $94.1\%$ is reported for the Pfizer-BioNTech and Moderna mRNA
vaccines. The $V\rightarrow I$ and $V\rightarrow A$ rates $\zeta$ and $\eta$
are derived by considering the order-of-magnitude ratio of the infected
individuals in the vaccinated and placebo cohorts of Ref. [37], with
$\beta_{2}$ and $\alpha_{2}$ set to $0.001$ and $0.01$ respectively, which are
in order-of-magnitude agreement with the aforementioned results of the
clinical trials. The $V\rightarrow R$ rate $\epsilon$ is the inverse time an
individual takes to acquire vaccine protection after the first shot. We set it
to $\epsilon^{-1}=21$ to reflect the fact that the second shot of vaccines is
usually administered about three weeks after the first one.

Figure 3:  Infectious population percentage (red dots) of some selected
countries in which a high percentage of their population has received a
vaccine shot. The infectious population is expressed as a function of time
(days). The date at which the wave began is pointed out on the horizontal
axis. The neural network fits and predictions are shown by black solid and
dashed lines respectively.  Table 1: Bundles and fixed parameters

$I_{0}$ | $A_{0}$ | $V_{0}$ | $R_{0}$ | $\beta_{1}$ | $\gamma$ | $\alpha_{1}$ | $\delta$
---|---|---|---|---|---|---|---
$[0.1~{}\%,2~{}\%]$ | $[0.1~{}\%,2~{}\%]$ | $[0~{}\%,60~{}\%]$ | $[0~{}\%,30~{}\%]$ | $[0.1,0.25]$ | $[0.07,0.12]$ | $[0.1,0.25]$ | $[0,0.03]$

$\epsilon$ | $\lambda$ | $\eta$ | $\zeta$ | $\beta_{2}$ | $\alpha_{2}$
---|---|---|---|---|---
1/21 | 0.95 | 1e-2 | 5e-3 | 1e-3 | 1e-2

Figure 4:  Infectious population percentage (red dots) of ten selected US
states. The infectious population is expressed as a function of time (days).
The date at which the wave began is pointed out on the horizontal axis. The
neural network fits and predictions are shown by black solid and dashed lines
respectively.

The remaining parameters $\Theta=(\alpha_{1},\beta_{1},\delta,\gamma)$
strongly depend on what kind of restrictive measures are taken or how fast the
vaccination campaign is, i.e. they are country dependent. We therefore
selected them in bundles so that the network could learn solutions
corresponding to a broad range of parameters and fit multiple countries. The
decay rate $\gamma$ is the inverse of the removal time which is about 1-2
weeks [39]. The main symptomatic infection rate $\beta_{1}$ was sampled in an
interval consistent with previous reports [2]. The earlier estimates that
$80\%$ of infected population is asymptomatic has been considered too high and
have since been revised down [28, 29]; the initial studies estimating this
proportion were limited by heterogeneity in case definitions, incomplete
symptom assessment, and inadequate retrospective and prospective follow-up of
symptoms. We selected the main asymptomatic infection rate $\alpha_{1}$ to be
varying in the same interval of $\beta_{1}$. The first shot vaccination rate
$\delta$ was selected based on known vaccination reports (see A).
Finally, the initial condition bundles
$\mathcal{Z}_{0}=(S_{0},A_{0},I_{0},V_{0},R_{0})$ are defined over broad
intervals able to cover the expected ($S,A,I,V,R$) populations at any given
time for all the cases considered. Although the initial infected population
$I_{0}$ is known, we still included it in the set of quantities to be fit by
the network. We found that by doing so, the network generalizes better
improving the fit of a given epidemic wave.

We then performed the fitting procedure described in 3.2 using the infectious
populations of 27 countries/states. The data (including the number of vaccine
shots administered) is retrieved from the ‘Our World in Data’ GitHub
repository [40]. Strictly speaking, the Infectious population of the SAIVR
model is the amount of people that are actively infected by the virus on a
given day, a number that should not be confused with the daily new cases. As
such, it was computed as the difference between the total number of cases and
of recovered/dead individuals.

We first applied the method to study the most recent Covid-19 wave in some of
the countries with the fastest vaccination campaigns and that managed to
inoculate the first shot in at least $30\%$ of their population. Fig. 3
presents real data (red points), fits (black solid line), and some predictions
(black dashed line) for the infectious populations of Israel, UK, Hungary,
France, Romania and Serbia. We also studied the USA, although due to the large
size of the country we focused only on the largest states or those with the
highest vaccination rates in the first quarter 2021, see Fig. 4. To assess the
generality of the model and the fitting procedure, we applied it to other 12
countries spread throughout the world and which had at the end of spring 2021
a high number of cases. The corresponding fits are shown in Figs. S1-S2 of the
Supplementary Material. As can be seen, the model is able to well reproduce
all these epidemic curves despite missing some abrupt and rapid events that
can be captured by more sophisticated multiple-wave models [2]. All the
parameters determined by the network can be found in the Supplementary
Material; their values are within the bundles of Table. 1. In particular,
$\gamma$ was found in the $~{}10-12$ days range, $\beta_{1}$ oscillating in
$[0.14-0.19]$ while alpha was more volatile.

##  4 Insights on the future: vaccine hesitancy, herd immunity and new
variants

Figure 5: Total infected population as a function of vaccination rate, vaccine
efficacy and vaccine denial population percentage. Results are obtained by
numerically solving the SAIVR model for $I_{0}=10^{-5}N_{pop}$ and
$A_{0}=0.2\times I_{0}$, where $N_{pop}=10^{6}$. The parameters of the model
used are those obtained by applying machine learning on the epidemic curves.
a) Infected population vs. vaccination rate $\delta$ and vaccine efficacy
$\lambda$.
b) Infected population as a function of the percentage of the population that
avoids getting vaccinated and the vaccine rate $\delta$.

In this section we use the results of the analysis performed in the previous
section to study how the vaccination campaign is affecting the pandemic and
its future evolution. Unless otherwise specified, we set $\beta_{1}=0.16$,
$\alpha_{1}=0.2$ and $\gamma=1/12$, the average values retrieved from fitting
real data. We start by pointing out how the vaccine efficacy is a key factor
in halting the spread of the virus and how hesitancy is challenging the
vaccination campaigns. Finally, we discuss the concept of herd immunity and
how it is affected by more infectious Covid-19 variants.

###  4.1 Vaccination efficacy and hesitancy

Figure 5 presents the total infected ($I+A$) population under increasing
values of vaccination onset times ($T_{0}$), vaccination daily rates
($\delta$), vaccine efficacy ($\lambda$) and of vaccine hesitancy/denial
population percentage. In the top panel, total infected population is shown as
a function of the vaccination rate $\delta$ and vaccine efficacy $\lambda$. As
can be seen, even vaccines with a relatively low efficacy can rapidly reduce
the infected population. In Fig. 5 b) we show how the number of those infected
evolves as a function of $\delta$ and the percentage of population that avoids
getting vaccinated. These findings suggest that vaccine hesitancy, which
accounts for a significant proportion of the population might seriously
threaten the reach of herd immunity, especially if the situation is worsened
by the appearance of more infectious Covid-19 strains.

###  4.2 Herd immunity and new Covid-19 variants

Figure 6:  Top row: Time evolution of epidemics following the introduction of
infected individuals in a population that has been already vaccinated at
50$\%$ (left panel), 60$\%$ (middle panel), and 70$\%$ (right panel), with
permanent immunity and no further vaccine roll-out, for four different numbers
of newly infected individuals ($I_{0}=1-0.001\%$ of total population, the
color code is presented in the legends). Results are obtained from numerically
solving the SAIVR model for $\beta_{1}=0.16$, $\alpha_{1}=0.20$, which
represent the average of the countries’ fitted values (see text for details
and for the values of the other parameters). Middle row: Same as in top row
but with $\beta_{1}=0.25$, $\alpha_{1}=0.25$, which represent a more
contagious Covid-19 variant. As is shown, vaccinated coverage of $50\%$ and
$60\%$ cannot prevent the resurgence of outbreaks. Bottom row: Same as in
middle row, but with ongoing vaccination roll-out with rate $\delta=0.001$,
per day, as it is shown, a continuing vaccine roll-out lowers the intensity of
the resurgent waves and prevents the resurgence of subsequent outbreaks.

The achievement of herd immunity has been hailed as the ultimate goal of a
successful vaccination campaign. ‘Herd immunity’, also known as ‘population
immunity’, is the indirect protection from an infectious disease that happens
when a sufficient portion of the population is immune either through
vaccination or immunity developed through previous infection. Once the herd
immunity threshold is met, the spread of the infectious disease is kept under
control, current outbreaks will extinguish and endemic transmission of the
pathogen will be interrupted. Earlier estimates of the threshold found values
of about $60-70\%$ of the population [41, 42, 43]. In reality, highly
transmissible strains tend to increase the threshold value, possibly keeping
this goal out of reach. Furthermore, persistent hesitancy about vaccines makes
vaccinating more than the $60-65\%$ of the population unlikely even in
countries which are at the global forefront of the vaccination effort.

We quantitatively investigate the likelihood of incurring resurgent Covid-19
epidemics after having immunized $50\%$, $60\%$, and $70\%$ of the population,
under different new infection introductions, Covid-19 variants and ongoing
vaccine deployment pace. Herd immunity protection is affected by the initial
value of the removed population ($R_{0}$ at $t=0$), which comprises both
recovered as well as fully vaccinated individuals, assuming permanent immunity
for both cases. In each scenario, we study the epidemic evolution after the
introduction of a cluster of newly infected individuals in the population.
$I_{0}$ represents the newly infected load at $t=0$
($I_{0}=1\%,0.1\%,0.01\%,0.001\%$ of the total population).

In the first part of the study we considered the less infectious variants
spreading during the 2021 spring by using the parameters retrieved in Sec.
3.3. The second scenario involves a more infectious Covid-19 strain such as
the Delta variant, which has been reported to be able to spread the virus more
efficiently [44]. Finally, we explore cases that involve or not further
(continuing) vaccine roll-outs.

Figure 6 presents the results obtained by numerically solving the SAIVR model
for the aforementioned cases. The top row presents the time evolution of an
outbreak in a population where the 50$\%$ (left panel), 60$\%$ (middle panel),
and 70$\%$ (right panel) of the individuals have been immunized, for different
numbers of initially infected individuals $I_{0}$. Results are obtained by
solving the SAIVR model for $\beta_{1}=0.16$, $\alpha_{1}=0.20$, and
$\gamma=1/12$; the average values obtained for the countries considered in
Sec. 3.3 and listed in Supplementary Material. As it can be seen, when the
immune portion of the population is only $50\%$, the outbreaks are contained
but not eradicated as the virus spreads in low intensity waves making the
disease endemic. Given the contagiousness of the less infectious variants, an
immunity threshold larger than $60\%$ is enough to eradicate the disease.
The middle row of Fig. 6 presents the evolution of outbreaks driven by a more
contagious variant; here $\beta_{1}=0.25$, and $\alpha_{1}=0.25$. As it is
shown, if the immunized portion of the population is only $50\%$ or $60\%$,
the resurgence of outbreaks cannot be prevented ($60\%$ immunity protection
makes the disease endemic). Only when the $70\%$ of the population is
immunized the disease is eradicated.
In both the top and the middle rows, the vaccine deployment is not taking
place during the outbreaks. The bottom row instead considers the highly
infectious variant disease evolution but with constant vaccine roll-out (with
rate $\delta=0.001$). As it can be seen, since even after getting only the
first vaccine shot individuals are partially protected, continuing the vaccine
administration rapidly lowers the intensity of the resurgent waves and helps
preventing subsequent outbreaks.

Although recent reports on highly infectious variants claim that the efficacy
of most vaccines is still about $90\%$ in preventing serious illnesses [45],
it is still not clear their performance on halting asymptomatic transmission.
In this study, we assumed that the vaccine efficacy on protecting from more
infectious variants is the same as for the less infectious ones. Despite this
optimistic assumption, the herd immunity threshold is moved to higher values
by simply increasing the infection rates.

##  5 Discussion

Compartmental models are efficient tools to deal with the time evolution of
disease outbreaks. They provide us with useful intuition on the impact of non-
pharmaceutical intervention in decreasing the number of infectious incidence
rates.
In this work, we have augmented the classic SIR model with the ability to
accommodate asymptomatic transmission and vaccinated individuals. The SAIVR
model is a straightforward deterministic model, which does not take into
consideration age, gender or geographic clustering. Despite this, its
simplicity and the insights it offers on how key epidemiological variables
affect individuals are among its main strengths. Its power also lies in the
fact that, as factors such new variants are added to the model, it is easy to
adjust its parameters and provide with best fit curves between the data and
the model predictions.
Since the inclusion of the Asymptomatic and Vaccinated compartments enlarged
the number of parameters and initial conditions of the model, we employed a
novel semi-supervised framework to estimate most of them. An unsupervised
neural network solves the model’s differential equations over a range of
parameters and initial conditions. A supervised approach then incorporates
data and determines the optimal initial conditions and modeling parameters
that best fit the 27 epidemic curves considered. As expected due to the
heterogeneity of the countries sample, the resulting parameters fit are
dissimilar although they follow similar trends.
We used these results to shed light on the impact of the vaccination campaign
on the future of the pandemic. We pointed out how vaccine hesitancy is one of
the most important hurdles of the campaign and further efforts should be done
to support people and give them correct information about vaccines. Because of
this, vaccinating the critical number of people that have to be immune in
order to prevent future outbreaks (i.e. herd immunity), is likely to be out of
reach. Widely circulating coronavirus variants are also a threat as they move
the herd immunity threshold to higher values. This points out the importance
of rapidly reducing the infection rate by any means, such as by imposing
restrictive measures in case highly infective new variants appears before the
herd immunity threshold is reached. These results manifest the need for
continuing the vaccination effort and the drive for achieving high vaccination
coverage in order to contain outbreaks generated by new and possibly more
infectious variants.

Data availability
The code used to perform the fitting is available on GitHub [46]. All study
data are either included in the article and supporting information or
available in Ref. [40].

Author Contribution
E.K. conceived the proposed model and supervised the study. M.A, G.N. and M.M.
designed the code. M.A. and G.N. performed numerical experiments, collected
data and analyzed the results. M.A. wrote the initial draft of the manuscript.
All authors critically revised, improved, and reviewed the manuscript in
various ways, and gave final approval for publication.

Competing interests
The authors declare no competing interest.

##  Appendix A Computational methods

We implemented a fully connected feed forward neural network that consists of
six hidden layers with 48 neurons per layer and Sigmoid activation functions.
The code is written in PyTorch [36] and published on GitHub [46]. In the
following we list the technical details of the learning procedure described in
Sec. 3:

###  A.1 Unsupervised learning

The network parameters (weights and biases) are updated with an adaptive
learning rate using the Adam algorithm [47] until the loss $\mathcal{L}$ of
Eq. 5 becomes smaller than $10^{-8}$. The learning rate ranged from 0.001 to
$10^{-}6$. The (initial) infectious and asymptomatic populations involved in
realistic situations are a small fraction of the total population. The ODEs in
Eqs. 2 are highly non-linear and extremely sensitive on initial conditions.
Subsequently, training the unsupervised network for the initial conditions
given in Tab. 1 is very challenging. To cope with this, we started by training
the network with a different choice of initial condition bundles taken in
‘safer’ intervals ($I_{0},A_{0}=[10\%,30\%]$) where the network was able to
quickly learn the solutions with high-accuracy. Then, we gradually decreased
$I_{0}$ and $A_{0}$ until they meet the values in Tab. 1, while we keep
training the network using the previously trained weights and biases. After
performing this annealing procedure, the network was able to learn the
solutions of the family of ODE defined by Tab. 1.

###  A.2 Supervised learning and fitting procedure

The SAIVR parameters and initial conditions are updated with a stochastic
gradient descent algorithm. Whenever the infection wave started before the
vaccination campaign, we divided the time series in two parts: before
($\delta=0$) and after ($\delta>0$) the vaccination campaign began. The first
shot vaccination rate $\delta$ can abruptly change in time based on logistics,
social or political reasons. Therefore, instead of fitting $\delta$ as we did
with the other parameters in the bundles we estimated it using real data, with
the aim of reproducing the average number of vaccinated people. The Vaccinated
compartment counts those individuals who have received the first shot but have
not been fully immunized yet. We defined the average number of people in the
vaccinated compartment $\tilde{V}$ in a given time interval as
$\tilde{V}=\int_{t_{0}}^{t_{f}}V(t)\epsilon~{}dt$, where $\epsilon$ is the
$V\rightarrow R$ rate. Then, we computed with a similar procedure
$\tilde{V}_{Data}$ using data relative to a given country and where
$V_{Data}(t)$ is the difference between the total number of people who got at
least one shot and those fully vaccinated. $\delta$ and $V_{0}$ were selected
to match $\tilde{V}$ with $\tilde{V}_{Data}$.
Although $I_{0}$ was provided in the data, we did not fix it during the
fitting. We have empirically observed that the network generalizes better in
this case.
The fitting procedure was performed 20 times for any given data-set in order
to select the parameters and initial conditions that fit the ground-truth data
with the lowest loss of Eq. 6.

## References

  * [1] W. O. Kermack, A. G. McKendrick, A contribution to the mathematical theory of epidemics, Proceedings of the Royal Society of London A: mathematical, physical and engineering sciences 115 (772) (1927) 700–721.
  * [2] E. Kaxiras, G. Neofotistos, Multiple epidemic wave model of the COVID-19 pandemic: Modeling study, Journal of Medical Internet Research 22 (7) (2020) e20912.  doi:10.2196/20912.
URL https://doi.org/10.2196/20912

  * [3] M. Barthélemy, A. Barrat, R. Pastor-Satorras, A. Vespignani, Dynamical patterns of epidemic outbreaks in complex heterogeneous networks, Journal of Theoretical Biology 235 (2) (2005) 275–288.  doi:10.1016/j.jtbi.2005.01.011.
URL https://doi.org/10.1016/j.jtbi.2005.01.011

  * [4] M. J. Ferrari, S. Bansal, L. A. Meyers, O. N. Bjørnstad, Network frailty and the geometry of herd immunity, Proceedings of the Royal Society B: Biological Sciences 273 (1602) (2006) 2743–2748.  doi:10.1098/rspb.2006.3636.
URL https://doi.org/10.1098/rspb.2006.3636

  * [5] E. Volz, SIR dynamics in random networks with heterogeneous connectivity, Journal of Mathematical Biology 56 (3) (2007) 293–310.  doi:10.1007/s00285-007-0116-4.
URL https://doi.org/10.1007/s00285-007-0116-4

  * [6] E. Tagliazucchi, P. Balenzuela, M. Travizano, G. Mindlin, P. Mininni, Lessons from being challenged by COVID-19, Chaos, Solitons & Fractals 137 (2020) 109923.  doi:10.1016/j.chaos.2020.109923.
URL https://doi.org/10.1016/j.chaos.2020.109923

  * [7] Q.-H. Liu, M. Ajelli, A. Aleta, S. Merler, Y. Moreno, A. Vespignani, Measurability of the epidemic reproduction number in data-driven contact networks, Proceedings of the National Academy of Sciences 115 (50) (2018) 12680–12685.  arXiv:https://www.pnas.org/content/115/50/12680.full.pdf, doi:10.1073/pnas.1811115115.
URL https://www.pnas.org/content/115/50/12680

  * [8] M. Zhang, A. Verbraeck, R. Meng, B. Chen, X. Qiu, Modeling spatial contacts for epidemic prediction in a large-scale artificial city, Journal of Artificial Societies and Social Simulation 19 (4) (2016) 3.  doi:10.18564/jasss.3148.
URL http://jasss.soc.surrey.ac.uk/19/4/3.html

  * [9] C. GROENDYKE, D. WELCH, D. R. HUNTER, Bayesian inference for contact networks given epidemic data, Scandinavian Journal of Statistics (2010) no–nodoi:10.1111/j.1467-9469.2010.00721.x.
URL https://doi.org/10.1111/j.1467-9469.2010.00721.x

  * [10] S. Sanche, Y. T. Lin, C. Xu, E. Romero-Severson, N. Hengartner, R. Ke, High contagiousness and rapid spread of severe acute respiratory syndrome coronavirus 2, Emerging Infectious Diseases 26 (7) (2020) 1470–1477.  doi:10.3201/eid2607.200282.
URL https://doi.org/10.3201/eid2607.200282

  * [11] Q. Li, X. Guan, P. Wu, X. Wang, L. Zhou, Y. Tong, R. Ren, K. S. Leung, E. H. Lau, J. Y. Wong, X. Xing, N. Xiang, Y. Wu, C. Li, Q. Chen, D. Li, T. Liu, J. Zhao, M. Liu, W. Tu, C. Chen, L. Jin, R. Yang, Q. Wang, S. Zhou, R. Wang, H. Liu, Y. Luo, Y. Liu, G. Shao, H. Li, Z. Tao, Y. Yang, Z. Deng, B. Liu, Z. Ma, Y. Zhang, G. Shi, T. T. Lam, J. T. Wu, G. F. Gao, B. J. Cowling, B. Yang, G. M. Leung, Z. Feng, Early transmission dynamics in wuhan, china, of novel coronavirus–infected pneumonia, New England Journal of Medicine 382 (13) (2020) 1199–1207, pMID: 31995857.  arXiv:https://doi.org/10.1056/NEJMoa2001316, doi:10.1056/NEJMoa2001316.
URL https://doi.org/10.1056/NEJMoa2001316

  * [12] Z. Du, L. Wang, S. Cauchemez, X. Xu, X. Wang, B. J. Cowling, L. A. Meyers, Risk for transportation of coronavirus disease from wuhan to other cities in china, Emerging Infectious Diseases 26 (5) (2020) 1049–1052.  doi:10.3201/eid2605.200146.
URL https://doi.org/10.3201/eid2605.200146

  * [13] C. Rothe, M. Schunk, P. Sothmann, G. Bretzel, G. Froeschl, C. Wallrauch, T. Zimmer, V. Thiel, C. Janke, W. Guggemos, M. Seilmaier, C. Drosten, P. Vollmar, K. Zwirglmaier, S. Zange, R. Wölfel, M. Hoelscher, Transmission of 2019-ncov infection from an asymptomatic contact in germany, New England Journal of Medicine 382 (10) (2020) 970–971.  arXiv:https://doi.org/10.1056/NEJMc2001468, doi:10.1056/NEJMc2001468.
URL https://doi.org/10.1056/NEJMc2001468

  * [14] L. Wynants, B. V. Calster, G. S. Collins, R. D. Riley, G. Heinze, E. Schuit, M. M. J. Bonten, D. L. Dahly, J. A. Damen, T. P. A. Debray, V. M. T. de Jong, M. D. Vos, P. Dhiman, M. C. Haller, M. O. Harhay, L. Henckaerts, P. Heus, M. Kammer, N. Kreuzberger, A. Lohmann, K. Luijken, J. Ma, G. P. Martin, D. J. McLernon, C. L. A. Navarro, J. B. Reitsma, J. C. Sergeant, C. Shi, N. Skoetz, L. J. M. Smits, K. I. E. Snell, M. Sperrin, R. Spijker, E. W. Steyerberg, T. Takada, I. Tzoulaki, S. M. J. van Kuijk, B. C. T. van Bussel, I. C. C. van der Horst, F. S. van Royen, J. Y. Verbakel, C. Wallisch, J. Wilkinson, R. Wolff, L. Hooft, K. G. M. Moons, M. van Smeden, Prediction models for diagnosis and prognosis of covid-19: systematic review and critical appraisal, BMJ (2020) m1328doi:10.1136/bmj.m1328.
URL https://doi.org/10.1136/bmj.m1328

  * [15] J. Koh, S. U. Shah, P. E. Y. Chua, H. Gui, J. Pang, Epidemiological and clinical characteristics of cases during the early phase of covid-19 pandemic: A systematic review and meta-analysis, Frontiers in Medicine 7 (2020) 295.  doi:10.3389/fmed.2020.00295.
URL https://www.frontiersin.org/article/10.3389/fmed.2020.00295

  * [16] F. Riccardo, M. Ajelli, X. D. Andrianou, A. Bella, M. Del Manso, M. Fabiani, S. Bellino, S. Boros, A. M. Urdiales, V. Marziano, M. C. Rota, A. Filia, O. Andrea, F. Trentini, G. Guzzetta, P. Poletti, P. Stefanelli, M. R. Castrucci, A. Ciervo, C. Di Benedetto, M. Tallon, A. Piccioli, S. Brusaferro, G. Rezza, S. Merler, P. Pezzotti, the COVID-19 working group, Epidemiological characteristics of covid-19 cases and estimates of the reproductive numbers 1 month into the epidemic, italy, 28 january to 31 march 2020, Eurosurveillance 25 (49) (2020).  doi:https://doi.org/10.2807/1560-7917.ES.2020.25.49.2000790.
URL
https://www.eurosurveillance.org/content/10.2807/1560-7917.ES.2020.25.49.2000790

  * [17] M. Khalili, M. Karamouzian, N. Nasiri, S. Javadi, A. Mirzazadeh, H. Sharifi, Epidemiological characteristics of covid-19: a systematic review and meta-analysis, Epidemiology and Infection 148 (2020) e130.  doi:10.1017/S0950268820001430.
  * [18] J. Li, D. Q. Huang, B. Zou, H. Yang, W. Z. Hui, F. Rui, N. T. S. Yee, C. Liu, S. N. Nerurkar, J. C. Y. Kai, M. L. P. Teng, X. Li, H. Zeng, J. A. Borghi, L. Henry, R. Cheung, M. H. Nguyen, Epidemiology of covid-19: A systematic review and meta-analysis of clinical characteristics, risk factors, and outcomes, Journal of Medical Virology 93 (3) (2021) 1449–1458.  arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/jmv.26424, doi:https://doi.org/10.1002/jmv.26424.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/jmv.26424

  * [19] C. Flamant, P. Protopapas, D. Sondak, Solving Differential Equations Using Neural Network Solution Bundles, arXiv e-prints (2020) arXiv:2006.14372arXiv:2006.14372.
  * [20] M. Mattheakis, D. Sondak, A. S. Dogra, P. Protopapas, Hamiltonian neural networks for solving differential equations (2020).  arXiv:2001.11107.
  * [21] A. Paticchio, T. Scarlatti, M. Mattheakis, P. Protopapas, M. Brambilla, Semi-supervised neural networks solve an inverse problem for modeling covid-19 spread (2020).  arXiv:2010.05074.
  * [22] M. M. Saito, S. Imoto, R. Yamaguchi, H. Sato, H. Nakada, M. Kami, S. Miyano, T. Higuchi, Extension and verification of the SEIR model on the 2009 influenza a (h1n1) pandemic in japan, Mathematical Biosciences 246 (1) (2013) 47–54.  doi:10.1016/j.mbs.2013.08.009.
URL https://doi.org/10.1016/j.mbs.2013.08.009

  * [23] H. Fang, J. Chen, J. Hu, Modelling the SARS epidemic by a lattice-based monte-carlo simulation (2005).  doi:10.1109/iembs.2005.1616239.
URL https://doi.org/10.1109/iembs.2005.1616239

  * [24] A. Smirnova, L. deCamp, G. Chowell, Forecasting epidemics through nonparametric estimation of time-dependent transmission rates using the SEIR model, Bulletin of Mathematical Biology 81 (11) (2017) 4343–4365.  doi:10.1007/s11538-017-0284-3.
URL https://doi.org/10.1007/s11538-017-0284-3

  * [25] S. A. Alanazi, M. M. Kamruzzaman, M. Alruwaili, N. Alshammari, S. A. Alqahtani, A. Karime, Measuring and preventing COVID-19 using the SIR model and machine learning in smart health care, Journal of Healthcare Engineering 2020 (2020) 1–12.  doi:10.1155/2020/8857346.
URL https://doi.org/10.1155/2020/8857346

  * [26] A. Palladino, V. Nardelli, L. G. Atzeni, N. Cantatore, M. Cataldo, F. Croccolo, N. Estrada, A. Tombolini, Modelling the spread of covid19 in italy using a revised version of the sir model (2020).  arXiv:2005.08724.
  * [27] I. Cooper, A. Mondal, C. G. Antonopoulos, A SIR model assumption for the spread of COVID-19 in different communities, Chaos, Solitons & Fractals 139 (2020) 110057.  doi:10.1016/j.chaos.2020.110057.
URL https://doi.org/10.1016/j.chaos.2020.110057

  * [28] D. Buitrago-Garcia, D. Egli-Gany, M. J. Counotte, S. Hossmann, H. Imeri, A. M. Ipekci, G. Salanti, N. Low, Occurrence and transmission potential of asymptomatic and presymptomatic sars-cov-2 infections: A living systematic review and meta-analysis, PLOS Medicine 17 (9) (2020) 1–25.  doi:10.1371/journal.pmed.1003346.
URL https://doi.org/10.1371/journal.pmed.1003346

  * [29] O. Byambasuren, M. Cardona, K. Bell, J. Clark, M.-L. McLaws, P. Glasziou, Estimating the extent of asymptomatic covid-19 and its potential for community transmission: systematic review and meta-analysis, medRxiv (2020).  arXiv:https://www.medrxiv.org/content/early/2020/09/13/2020.05.10.20097543.full.pdf, doi:10.1101/2020.05.10.20097543.
URL https://www.medrxiv.org/content/early/2020/09/13/2020.05.10.20097543

  * [30] M. K. Hossain, M. Hassanzadeganroudsari, V. Apostolopoulos, The emergence of new strains of sars-cov-2. what does it mean for covid-19 vaccines?, Expert Review of Vaccines 0 (0) (2021) 1–4, pMID: 33896316.  arXiv:https://doi.org/10.1080/14760584.2021.1915140, doi:10.1080/14760584.2021.1915140.
URL https://doi.org/10.1080/14760584.2021.1915140

  * [31] R. Wang, D. Maddix, C. Faloutsos, Y. Wang, R. Yu, Bridging physics-based and data-driven modeling for learning dynamical systems (2021).  arXiv:2011.10616.
  * [32] Z. Yang, Z. Zeng, K. Wang, S.-S. Wong, W. Liang, M. Zanin, P. Liu, X. Cao, Z. Gao, Z. Mai, J. Liang, X. Liu, S. Li, Y. Li, F. Ye, W. Guan, Y. Yang, F. Li, S. Luo, Y. Xie, B. Liu, Z. Wang, S. Zhang, Y. Wang, N. Zhong, J. He, Modified seir and ai prediction of the epidemics trend of covid-19 in china under public health interventions, Journal of Thoracic Disease 12 (3) (2020).
URL https://jtd.amegroups.com/article/view/36385

  * [33] D. Zou, L. Wang, P. Xu, J. Chen, W. Zhang, Q. Gu, Epidemic model guided machine learning for covid-19 forecasts in the united states, medRxiv (2020).  arXiv:https://www.medrxiv.org/content/early/2020/05/25/2020.05.24.20111989.full.pdf, doi:10.1101/2020.05.24.20111989.
URL https://www.medrxiv.org/content/early/2020/05/25/2020.05.24.20111989

  * [34] I. Ayed, E. D. Bézenac, A. Pajot, P. Gallinari, Learning partially observed PDE dynamics with neural networks (2019).
URL https://openreview.net/forum?id=HyefgnCqFm

  * [35] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, D. K. Duvenaud, Neural ordinary differential equations, in: S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, R. Garnett (Eds.), Advances in Neural Information Processing Systems, Vol. 31, Curran Associates, Inc., 2018.
URL
https://proceedings.neurips.cc/paper/2018/file/69386f6bb1dfed68692a24c8686939b9-Paper.pdf

  * [36] A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, A. Lerer, Automatic differentiation in pytorch (2017).
  * [37] F. P. Polack, S. J. Thomas, N. Kitchin, J. Absalon, A. Gurtman, S. Lockhart, J. L. Perez, G. Pérez Marc, E. D. Moreira, S. Zerbini, P. R. Dormitzer, U. Şahin, K. U. Jansen, W. C. Gruber, Safety and efficacy of the bnt162b2 mrna covid-19 vaccine, New England Journal of Medicine 383 (27) (2020) 2603–2615, pMID: 33301246.  arXiv:https://doi.org/10.1056/NEJMoa2034577, doi:10.1056/NEJMoa2034577.
URL https://doi.org/10.1056/NEJMoa2034577

  * [38] L. R. Baden, H. M. El Sahly, B. Essink, K. Kotloff, S. Frey, R. Novak, D. Diemert, S. A. Spector, N. Rouphael, C. B. Creech, J. McGettigan, S. Khetan, N. Segall, J. Solis, A. Brosz, C. Fierro, H. Schwartz, K. Neuzil, L. Corey, P. Gilbert, H. Janes, D. Follmann, M. Marovich, J. Mascola, L. Polakowski, J. Ledgerwood, B. S. Graham, H. Bennett, R. Pajon, C. Knightly, B. Leav, W. Deng, H. Zhou, S. Han, M. Ivarsson, J. Miller, T. Zaks, Efficacy and safety of the mrna-1273 sars-cov-2 vaccine, New England Journal of Medicine 384 (5) (2021) 403–416, pMID: 33378609.  arXiv:https://doi.org/10.1056/NEJMoa2035389, doi:10.1056/NEJMoa2035389.
URL https://doi.org/10.1056/NEJMoa2035389

  * [39] N. George, N. K. Tyagi, J. B. Prasad, Covid-19 pandemic and its average recovery time in indian states, Clinical Epidemiology and Global Health 11 (2021) 100740.  doi:https://doi.org/10.1016/j.cegh.2021.100740.
URL https://www.sciencedirect.com/science/article/pii/S2213398421000440

  * [40] Our world in data github repository: https://github.com/owid/covid-19-data (2021).
URL https://github.com/owid/covid-19-data

  * [41] S. B. Omer, I. Yildirim, H. P. Forman, Herd Immunity and Implications for SARS-CoV-2 Control, JAMA 324 (20) (2020) 2095–2096.  arXiv:https://jamanetwork.com/journals/jama/articlepdf/2772167/jama\\_omer\\_2020\\_it\\_200033\\_1605821558.92096.pdf, doi:10.1001/jama.2020.20892.
URL https://doi.org/10.1001/jama.2020.20892

  * [42] S. Chowdhury, S. Roychowdhury, I. Chaudhuri, Universality and herd immunity threshold: Revisiting the sir model for covid-19, International Journal of Modern Physics C 0 (0) (0) 2150128.  arXiv:https://doi.org/10.1142/S012918312150128X, doi:10.1142/S012918312150128X.
URL https://doi.org/10.1142/S012918312150128X

  * [43] R. Aguas, R. M. Corder, J. G. King, G. Gonçalves, M. U. Ferreira, M. G. M. Gomes, Herd immunity thresholds for sars-cov-2 estimated from unfolding epidemics, medRxiv (2020).  arXiv:https://www.medrxiv.org/content/early/2020/11/16/2020.07.23.20160762.full.pdf, doi:10.1101/2020.07.23.20160762.
URL https://www.medrxiv.org/content/early/2020/11/16/2020.07.23.20160762

  * [44] P. Arora, A. Kempf, I. Nehlmeier, A. Sidarovich, N. Krüger, L. Graichen, A.-S. Moldenhauer, M. S. Winkler, S. Schulz, H.-M. Jäck, M. V. Stankov, G. M. N. Behrens, S. Pöhlmann, M. Hoffmann, Increased lung cell entry of b.1.617.2 and evasion of antibodies induced by infection and bnt162b2 vaccination, bioRxiv (2021).  arXiv:https://www.biorxiv.org/content/early/2021/06/23/2021.06.23.449568.full.pdf, doi:10.1101/2021.06.23.449568.
URL https://www.biorxiv.org/content/early/2021/06/23/2021.06.23.449568

  * [45] J. L. Bernal, N. Andrews, C. Gower, E. Gallagher, R. Simmons, S. Thelwall, J. Stowe, E. Tessier, N. Groves, G. Dabrera, R. Myers, C. Campbell, G. Amirthalingam, M. Edmunds, M. Zambon, K. Brown, S. Hopkins, M. Chand, M. Ramsay, Effectiveness of covid-19 vaccines against the b.1.617.2 variant, medRxiv (2021).  arXiv:https://www.medrxiv.org/content/early/2021/05/24/2021.05.22.21257658.full.pdf, doi:10.1101/2021.05.22.21257658.
URL https://www.medrxiv.org/content/early/2021/05/24/2021.05.22.21257658

  * [46] Mlcd machine learning covid dynamics https://github.com/mattangeli/mlcd-machine-learning-covid-dynamics , year = 2021,.
URL https://github.com/mattangeli/MLCD-Machine-Learning-Covid-Dynamics

  * [47] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization (2017).  arXiv:1412.6980.