COVASIM was developed under leadership of the Institute for Disease Modeling and released in May 2020 by Kerr and colleagues [35]. It is one of several open-source ABMs, together with OpenABM-Covid19 [36] or COMOKIT [37]. The model captures the transition from susceptible to infected followed by a split between asymptomatic individuals and various degrees of symptoms, resulting either in recovery or death (Figure 1). The model was created to support interventions offered at the time, which did not include vaccination. We thus modified the model to account for our current understanding of viral dynamics and the use of vaccines over two doses (Figure 1). When instantiating the model to the US population, we used a resolution of 1:500 (ie, each simulated agent accounts for 500 US inhabitants). Given our resolution and target population size, our application exceeded half a million agents and can thus be described as a “large-scale COVID-19 simulation” [38]. Our simulations started on January 1, 2020, using CDC data for the number of infected, recovered, and immunized individuals to date (see subsection Initializing the Model). We then simulate for 6 months, that is, 180 time ticks based on a temporal resolution of 1 day per simulation step (ie, tick). To cope with the computational challenges created by a large-scale stochastic model, a philanthropic grant supports us in performing cloud-based simulations via the Microsoft Azure (Microsoft Corporation) platform.

Apart from being open source, there are two reasons that we selected COVASIM. First, it captures heterogeneity within individuals (eg, assigns an age and uses age-specific disease outcomes) and transmission patterns by placing agents within synthetic networks corresponding to a multiplicity of contexts: work (based on employment rates), school (based on enrollment), home (based on household size), and the general community. However, these high-resolution age-specific contact patterns are not unique to COVASIM. For example, the OpenABM-Covid19 [36] also embeds agents in age-stratified occupation networks (encompassing work and school), household networks, and a general random network. COMOKIT [37] similarly uses the Gen* toolkit from the same team to redistribute populations from census units down to exact buildings such as the nearest school. Thus, the second rationale for choosing this platform is that it has been used in the most peer-reviewed modeling studies to date [39,40], hence providing an additional layer of scrutiny and confidence in the correctness of the model (ie, validation) and its implementation (ie, verification). As detailed in our recent study [34], changes in the evidence base have required alteration in the model to keep it valid. Consequently, we modified three COVASIM parameters to account for the current biological and epidemiological evidence on COVID-19 (Table 1).

In addition to support for heterogeneity, COVASIM implements several nonpharmaceutical interventions. Although our focus is on vaccines, such interventions may be continuing in parallel with the vaccination campaign; hence, we have to take them into account when forecasting case counts. Interventions can be organized into three broad categories: preventative care (eg, social distancing and face masks), lockdown (eg, stay-at-home orders such as remote work or school closures), or testing-related (eg, testing itself, then quarantining and contact tracing) [14,41,42]. In line with our previous work on nonpharmaceutical interventions, we considered all 6 specific interventions. Although all 6 are natively supported by the COVASIM platform, we changed testing delays from their default value (constant) to a distribution (based on a survey across all 50 US states) [43], thus accounting for the variability observed in practice.

Since our focus is on vaccines, our search space is primarily devoted to quantifying the effect of vaccine-related variables (ie, efficacy, compliance, and capacity). As every nonpharmaceutical intervention could lead to several variables (eg, compliance with face masks or efficacy of face masks), considering all variables for every such intervention in addition to vaccine-related variables would lead to an impractical search space. We thus leveraged the systematic assessment of our previous study [34], which simulated all combinations of nonpharmaceutical interventions at two different levels of strength (ie, a binary factorial design of experiments). We analyzed results from this broad search to select 5 scenarios (Table 2) that resulted in five different levels of infection count after 6 months, in the absence of any vaccine (Figure 2). In other words, to circumvent the unwieldy notion of simulating all aspects of vaccines and nonpharmaceutical interventions, we selected 5 scenarios that produce linear to logistic growths in cumulative infections, thereby conducting a parameter sweep across possible growth behaviors. We supplemented these 5 scenarios with an extreme no intervention scenario, which provides an upper bound on the number of cases.

Given that we made minor changes to the biology (incubation and proportion of symptomatic cases) and consider several ongoing intervention scenarios, it is necessary to confirm the validity of the model established using earlier data in previously published studies. Consequently, we ran the modified COVASIM model based on data observed until September 3, 2020, and compared the simulated results with observations until the end of year. Similar trends and orders of magnitude were observed (Figure 3), thus providing qualitative validation. Note that the 5 scenarios chosen (Table 2) bound the growth of COVID-19 in the United States such that we are comprehensively examining possible trends going forward instead of limiting ourselves to the single trend that fit best on previous data.

As detailed in our discussion, there is substantial uncertainty and frequent changes regarding the number of vaccines that may be administered monthly. We thus considered two vaccine availability scenarios, both proposed by federal governments. The first scenario from the former Trump administration, named Operation Warp Speed, stated that vaccines will be available in tiered amounts (20 million in December, 30 million in January, and 50 million every month thereafter). The second scenario from the Biden administration, known as the 100-day goal, proposes that there will be 1 million vaccines every day [44], thus covering 50 million Americans. Although there are other scenarios, they vary from state to state (eg, the governor of New Jersey aspires to vaccinate 70% of the adult population within 6 months [45]) and are subject to frequent revisions. Given the countrywide nature of our simulation, we relied on federal plans while detailing challenges (see also the Discussion section).

In setting the monthly capacity, we noticed the necessity to adjust the schedule of the Operation Warp Speed plan, since the initial aim of 20 million people immunized by the end of December 2020 only resulted in 3 million doses administered. In other words, it would be incorrect to model the monthly capacity of Operation Warp Speed as announced since there is evidence that its initial objective was unmet, due to a variety of logistical challenges. Consequently, we shifted the expectations of the Operation Warp Speed plan by 1 month, such that the capacity for January now corresponds to the initial expectations for December (20 million) and so forth.

At the same time as either vaccination schedule is active, we also have the 6 scenarios listed in the previous sections. As these scenarios include a no-intervention case, we are able to study the interaction between nonpharmaceutical interventions and vaccines. In total, this gives 12 distinct situations. In addition, we also varied two essential parameters regarding vaccines: the percentage of the population that seeks vaccination (which we refer to as vaccine compliance from hereon) and the efficacy of the vaccine. Varying these two parameters across 12 situations in a large-scale ABM results in substantial computing needs. These are challenging to parallelize, as the run time of each experiment is not the same. Therefore, we took advantage of the massive parallelism enabled by the cloud computing platform Azure to accelerate computation. Using this platform, we varied vaccine compliance and vaccine efficacy between the bounds listed in Table 3.

Regarding our approach to vaccine efficacy, we noted that individuals can be infected after their first dose, as has been documented on thousands of cases [46]. We thus used the probability of 52% (observed in clinical trials [47]) to obtain early protection by the vaccine, and otherwise, an individual may still be infected in the waiting period leading to the second dose. After the second dose is applied, we needed to ensure that the agent meets the vaccine efficacy set by our parameters. That is, the probability of obtaining immunity after the second dose was set such that the probability of immunity from the two doses matches the vaccine efficacy.

Although we did not track which of the two approved mRNA COVID-19 vaccines (Pfizer-BioNTech or Moderna) were administered, we varied vaccine efficacy to account for a margin of uncertainty regarding their respective performances. Since the vaccine capacity is either planned to increase (Operation Warp Speed) or be at a high constant rate, a simulated agent given one dose will always be able to come back to get the second dose on time. Should an agent be contaminated or die before the second dose, it is then released for administration to another agent.

We also varied the percentage of the population who seeks vaccination. As noted in a recent study, this percentage has varied among studies: 10.8% did not intend to be vaccinated when asked in April 2020, but this number jumped to 31.1% by May, and an August poll found that only a minority would want to be vaccinated [48]. In addition to changes in the sociopolitical climate and public discourse surrounding vaccination, there will also be changes since “many receptive participants preferred to wait until others have taken the vaccine” [49]. Seeing positive vaccination outcomes in others may in part address the fear of serious side effects, which is a recurring concern for individuals who may not intend to participate in vaccination [50]. Given past variations and changes in the future, we handled uncertainty through a parameter sweep in vaccine compliance.

A simulation model is composed of an initialization (setting characteristics of agents for t=0) and rules governing its update, thereby producing data for analysis. The previous subsections covered the rationale for the inclusion of agents’ characteristics and the design of the rules, while the next subsection focuses on the analysis. This subsection thus briefly covers our approach to initialization such that our results could be independently replicated by other modeling teams.

Our initial time tick t=0 corresponds to January 1, 2020. We thus needed to set the number of agents who have been infected, recovered, or immunized (due to the rollout of vaccines in December) by that time. A COVID-19 case remains infectious within a time window of 2 weeks, after which there is either recovery or complications. From December 18-31, there was a total of 3,311,345 active cases. To appropriately initialize our simulation, we needed to further track when an individual was infected. Incorrectly setting them to be all infected on December 18 would result in nobody being infected when the simulation starts on January 1. At the other extreme, assuming that they were all infected on December 31 would lead to an overestimate of disease spread into 2021. We thus seeded the timing of each infection by using the daily distribution from CDC data between December 18-31 (Table 4). All numbers were divided by 500 since our agent resolution is 1 agent for 500 real-world US inhabitants (1:500). The number of individuals who acquired immunity via recovery was set to the total case count observed by December 17. Individuals who died from COVID-19 are grouped together with recovered ones (ie, we did not subtract them from the count) since our simulations track the number of new infections; dead individuals do not alter these results as they can neither be infected nor infect others. The total number of individuals immunized from vaccination was set to 2 million (ie, 4000 agents).

To quantify the spread of the disease, we fitted the progression of cumulative infection to a logistic growth model, which is a simple yet effective model describing resource-limited growths in natural processes and has been used on several occasions for COVID-19 [51-53]. Let the cumulative infection be P = P(t), then the logistic model stipulates that P is the solution of the differential equation:

where is the time derivative of P, r is the growth rate (proportional to the maximum value attained by ), and K is the carrying capacity. As our simulations produce the complete time series for P, we can estimate using finite differences, thereby extracting parameters r and K through a linear regression as equation 1 suggests. In the regression, the independent and dependent variables are P and / P, respectively. In addition, we measured the goodness of fit as that of the linear regression. Since the simulation is stochastic, multiple replications are needed for each configuration to obtain an average behavior. We used the CI method [54] to perform enough replications so that for every time step t, the 95% CI of P at time t falls within 5% of the average. Therefore, we performed the fitting for each individual run and computed the average r and K across all runs.

Although we report the carrying capacity K in Multimedia Appendices 1 and 2, the interpretation of this variable can be difficult for a broader audience. The growth rate r is proportional to the maximum fraction of the carrying capacity K that is infected on the worst day. In other words, it is an indication of how fast the disease spreads at its peak, based on another variable. For ease of interpretation, we focused on the adjusted growth rate whose unit is directly in number of individuals. The adjusted growth rate reported in this paper is obtained as:

For instance, an adjusted value of 200,000 means that at most 200,000 individuals will be infected on the worst day.

As the early steps of the simulation witness a shift from a vaccine-naïve population to one that gradually builds vaccine-based immunity, early trends differ from the longer ones that are the focus of this study. This is a typical situation in modeling, whereby estimating the long run performance measures requires to first run the model for a certain amount of time (known as the warm-up period) [55]. We empirically determined that a warm-up period of 20 days was sufficient to start the curve fitting; that is, we created the time series for P starting from t≥20. As evidenced by Figure 4, this warm-up period results in very good fit for the logistic model under both federal plans. This approach also generalizes better, since the reported r and K can accurately characterize the spread of the disease for most time periods instead of being skewed by the first few days.

An essential aspect of a return to normalcy is about the conditions under which that is achieved. If the disease is left uncontrolled, and simplifying the matter of variants, we would still return to normalcy within 6 months because a large share of the population would already have been infected and either recovered or died (Figure 5). The goal is thus not only to eventually achieve stability in the number of cases but to achieve it at a minimal level (Figure 5; bottom blue curve).

To efficiently orchestrate simulations over the Microsoft Azure cloud computing platform, we used a distributed scheme shown in Figure 6. The setup starts by creating a manager, which uses queues to organize the two types of work that need to be performed.

Available workers contact the manager, who will assign work (Figure 6a) by prioritizing the job queue and then the timing queue. For example, if a worker notifies the manager that it is available and there is a simulation run to perform in the job queue, then the manager will hand that one run to the worker (Figure 6b). If a worker is available and all queued simulations have been performed, then the manager will task the worker with identifying how many simulations are necessary for the next configuration (Figure 6c), which will refill the job queue.

The incoming CDC director predicted half a million deaths by mid-February 2021 [56], thus stressing the urgency of vaccination. However, vaccination is an unprecedented and complex endeavor whose success depends on many other variables such as vaccine compliance, vaccine efficacy, and the ongoing presence of nonpharmaceutical interventions. In line with expectations, our large-scale agent-based simulations showed that vaccination can reduce the total number of infections across all possible scenarios. The capacity pledged under the new Biden plan (one million doses a day) would have a greater impact than the plan of the previous administration (Operation Warp Speed) when accounting for its initial delays.

Two key findings of our study are as follows. First, we demonstrated the necessity to maintain nonpharmaceutical interventions over the next 6 months. As interventions are relaxed (from scenario 1 offering the most control to scenario 6 offering no control), there is an increase in case count such that a return to normalcy is not achieved through vaccination but rather through a very high number of infected individuals. Second, there is an unexpected interplay between vaccination strategies, nonpharmaceutical interventions, and vaccination availabilities. As nonpharmaceutical interventions lose momentum (scenarios 3 and above), an increase in vaccine compliance leads to an unexpected increase in infections due in part on the low availability of vaccines and the priority on vaccinating older adults. More so than the observation that tighter nonpharmaceutical interventions result in the slower spread of infections, this result further delineates the necessity of preparing the population to continuing nonpharmaceutical interventions even as the vaccination progresses.

There are three main limitations to our current understanding of the COVID-19 pandemic and the vaccination campaign that affect how our simulations account for (1) the number of vaccines that can be administered each month, (2) biological aspects, and (3) healthy or asymptomatic carriers.

First, an unprecedented vaccine campaign comes with logistical challenges and uncertainty given the complex array of factors involved. As a result, fewer than the expected number of doses may be administered: federal officials aimed at giving the first dose to 20 million people during December 2020, but various delays resulted in fewer than 3 million people receiving a first dose [57]. It was recently reported that “federal officials say they do not fully understand the cause of the delays” [57] and that the administration “pledged to immediately distribute millions of COVID-19 vaccine doses from a stockpile that the U.S. health secretary has since acknowledged does not exist” [58]. This situation has resulted in views that “much of the narrative earlier this year regarding Warp Speed’s preparation appears to be a sham” [59], reinforced by reports that the Biden administration found no vaccine distribution plan upon taking over from their predecessors [60]. Some of the factors causing a delay are known: there can be shipping delays or delays in administering doses due to a lack of hospital staff members, as they are already caring for individuals infected with COVID-19. Other factors may be more surprising, such as the intentional destruction of vaccine doses by hospital staff [61]. As any simulation model is necessarily a simplification, we did not include factors whose value would be entirely unknown (eg, what will be the shipment delay?) or whose existence is anecdotal given the total number of doses (eg, intentional destruction or storage errors). We were limited in our ability to use real-world numbers on how many individuals received the vaccine, as this data is captured at the state level, and several states’ reporting systems have experienced errors [62]. Although there are efforts at centralizing data (eg, national news outlets aggregate data across states [63]), the level and nature of errors differ across states, which is a challenge to estimate overall model uncertainty.

We have thus followed the federal plan for the number of individuals who can get vaccinated each month. Out of all the doses that are planned, fewer may be distributed and an even lower number may ultimately be administered. Our simulations are thus likely representing an upper bound on the number of vaccines administered, leading to more optimistic results than in reality. The gap is particularly pronounced in December 2020 and may remain significant in January 2021, but early logistical issues and delays should be gradually addressed, such that the gap between federal expectations and actual implementation narrows over time.

Second, all biological aspects of the virus are based on the strains that dominated throughout 2020. Epidemiological studies from these strains have informed parameters such as transmissibility, incubation period, the proportion of asymptomatic carriers, the severity of symptoms and hence the course of the disease, and the efficacy of treatments or vaccines. The existence of different strains is well established, as phylogenies have shown seven distinct lineages [64,65], but there has not yet been a documented need to ascribe different parameter values (ie, different viral behaviors) to each strain. There are two possible reasons. First, there are relatively few mutations and thus a limited chance of a drastically different outcome naturally occurring: the virus is “considered a slowly-evolving virus as it possesses an inherent proofreading mechanism to repair the mismatches during its replication” [65]. Second, there has been little selective pressure on the virus, as it was spreading through a population that had never been exposed to an antigen (ie, immunologically naïve). Both arguments are now changing.

A new strain from the lineage B.1.1.7, named Variant of Concern 202012/01 (denoted VOC-202012/01), emerged with an unusually large number of 23 changes in its genomes (including mutations and deletions) [66]. Some of the biological changes make it easier for the virus to attach to its targets and enter cells, which is captured through epidemiological indicators as increased transmissibility [67,68]. This is relevant for our study, as this more contagious COVID-19 strain has been spreading in the United States and may dominate by March 2020 [69]. To date, there is no peer-reviewed evidence of an impact on disease severity or vaccine efficacy over a large population sample, but the function for some of the mutated parts remains unknown (hence the possibility of an impact on severity), and early studies over 20 volunteers suggest that antibodies from vaccines are only one-third as effective on some variants [70]. In addition, vaccination means that the virus is no longer spreading through an immunologically naïve population, thus creating selective pressure for functional mutations that can help the virus adapt. Our simulation results are thus optimistic as they use a lower transmissibility than provided by the new strain, and we did not worsen any of the other parameters to account for possible selective pressure.

Third, our model considers that individuals who were successfully immunized can act as a buffer in the spread of the epidemic. Reality may be more nuanced, as viral transmission from a vaccinated host to an unvaccinated one may be possible. At the time of writing (March 2021), we do not yet have conclusive findings about this possibility. As trials continue, we may find that immunized individuals should be treated in a model as healthy carriers for a period. We also noted that the immunity conferred by the vaccine appears to have a different response than the immunity acquired by recovering from a natural infection. That is, a vaccine promotes the production of antibodies in the blood, but a natural immunity may lead to developing antibodies in the mucosal regions [71], which are the first site of infection (in the nose and mouth). From a modeling viewpoint, the two immunities may thus have to be treated differently in the future.

Finally, we note that our model is built specifically for the United States. It would not be accurate when transposed to another country with minimal changes (eg, only reducing the population size). For example, stark differences in vaccine rollout strategies exist between the United Kingdom and the United States, which would affect our simulations. In the United States, two doses of the same vaccine are normally administered, as the CDC stated that “mRNA COVID-19 vaccines are not interchangeable with each other or with other COVID-19 vaccine products” [72]. However, new guidance from the United Kingdom allows a mix-and-match vaccine regimen in which the second dose may be from a different vaccine in exceptional circumstances (eg, if the vaccine from the first dose is not available upon the patient’s return), even though clinical trials for mixed regimens are only due to be conducted at a later, unspecified time [73]. Another difference is that the United Kingdom front-loads the vaccine by delivering as many first doses as possible, which thus no longer guarantees that a patient can receive the corresponding second dose upon return (hence raising the need for a mix-and-match) and potentially delays the delay before a second dose up to 12 weeks [73]. In contrast, the United States is against delaying the second dose [74], thus our model operates on the assumption that a patient can complete treatment on time.

Our simulation of half a million agents qualifies as large-scale in the context of COVID-19 ABMs. In another context, the scale may be different as the computational costs of the simulation or historical practices in a research community may differ. For example, in HIV research, simulations have used half a million cells for about 20 years on personal computers, so a large-scale may be a more appropriate qualifier for simulations with billion cells [75,76]. As noted by Gumel and colleagues [77] in their extensive discussion on modeling methods for COVID-19, ABMs “are computationally-intensive”; thus, we may expect a smaller simulated population than in compartmental models or meta-population models, given the same hardware and simulation time.

Many ABMs for COVID-19 are in the scale of several hundred agents [78-83] to tens of thousands of agents [37,84,85]. Fewer studies have over 100,000 agents [86], and only a paucity of studies has a number of agents that is about equal (eg, the model of Hoertel and colleagues [87] used 500,000 agents) or greater than (eg, one million agents in a February 2021 simulation of Bogota) in this study [38,87,88]. Due to this distribution of agent population across studies, the qualifier of large is applied as we get to the scale of 500,000 or more agents [38]. It should not be interpreted to suggest that this is the largest population size achieved to date. Indeed, a few high-profile studies have modeled their target populations with such a fine resolution that the simulation may qualify as a digital twin. For example, Chang et al [89] used over 24 million agents by adding a COVID-19 component (AMTraC-19) to an existing model and running it over 4264 compute cores.

Although a justification for the scale is a recommended best practice in ABM for artificial societies [90], such a justification is not always present in published studies. The studies that have justified their choice of scale have often done it based on the size of the target population (eg, single city or campus) or implicitly invoked the notion of a computational burden when downscaling. Explicit mentions of computational costs have been made by the developers of frameworks, such as Comokit, who stated that 10-20,000 agents could be simulated on one laptop within 10 minutes [37].

A desirable return to normalcy would be achieved via immunization rather than through a very high number of infected cases and their natural immunity. Our extended ABM shows that vaccines are not sufficient to return to normalcy while avoiding a high number of cases. Nonpharmaceutical interventions are necessary and require a high level of compliance to ensure that immunity from vaccination outpaces the immunity from infections. Although our findings account for different vaccination capabilities, compliance levels, and vaccine efficacy, they are nonetheless based on a simulation model, which is necessarily a simplification of reality. Simplifications here include the logistics of vaccine dissemination, variants, and the presence of healthy carriers (vaccinated) and asymptomatic cases (not vaccinated).