key: cord-0944770-5oj6bqta
authors: Getz, W. M.; Salter, R.; Vissat, L. L.; Koopman, J. S.; Simon, C. P.
title: Adaptive vaccination may be needed to extirpate COVID-19: Results from a runtime-alterable strain-drift and waning-immunity model
date: 2021-06-11
journal: nan
DOI: 10.1101/2021.06.07.21258504
sha: 537fd9fab77122726092763cd562bdba0f46d8e8
doc_id: 944770
cord_uid: 5oj6bqta

We developed an elaborated susceptible-infected-recovered (SIR) individual-based model (IBM) with pathogen strain drift, waning and cross immunity, implemented as a novel Java Runtime-Alterable-Model Platform (J-RAMP). This platform allows parameter values, process formulations, and scriptable runtime drivers to be easily added at the start of simulation. It includes facility for integration into the R statistical and other data analysis platforms. We selected a set of parameter values and process descriptions relevant to the current COVID-19 pandemic. These include pathogen-specific shedding, environmental persistence, host transmission and mortality, within-host pathogen mutation and replication, adaptive social distancing, and time dependent vaccine rate and strain valency specifications. Our simulations illustrate that if waning immunity outpaces vaccination rates, then vaccination rollouts may fail to contain the most transmissible strains. Our study highlights the need for adaptive vaccination rollouts, which depend on reliable real-time monitoring and surveillance of strain proliferation and reinfection data needed to ensure that vaccines target emerging strains and constrain escape mutations. Together with such data, our platform has the potential to inform the design of vaccination programs that extirpate rather than exacerbate local outbreaks. Finally, our RAMP concept promotes the development of highly flexible models that can be easily shared among researchers and policymakers not only addressing healthcare crises, but other types of environmental crises as well.

The amount of structure and data needed in complex biological systems' mod-284 els, depends on the questions that these models have been formulated to ad-285 dress [26, 37] . In this paper, we steered away from making specific predictions-286 because universal solutions are not always locally applicable. Rather, we focused 287 on questions relating to the potential efficacy of different vaccination rollouts 288 (both vaccination rates and valencies of vaccines) in the context of strain emer-289 gence and the potential for vaccination programs to go awry. To make specific 290 predictions requires localized data, particularly as it may relate to social dis-291 tancing behavior or other factors that affect contact rates over time [30, 38, 39] . 292 Additionally, we cannot expect models to make strain-specific informative pre-293 dictions, unless they have been designed to do so and are supported by location 294 specific data regarding the relative transmissibility of strains and other strain 295 specific data. In the absence of such data, models become a tool for anticipat- 296 ing pitfalls and avoiding unintended consequences of well-intentioned actions. 297 Equally important, however, is the implementation of our model as a RAMP 298 (runtime alterable model platform), because this greatly facilitates the use of 299 our model by ourselves and others in testing out different hypothesis about the 300 process generating observed population level strain transmission dynamics. 301 The illustrative simulations we present demonstrate that if vaccination roll-302 out rates and vaccine valencies are applied in blanket manner-that is, without 303 monitoring strain emergence and obtaining some assessment of the epidemio-304 logical characteristics of those strains (e.g., relative transmissibility, virulence, 305 waning and cross immunity characteristics)-then the possibility exists for a 306 vaccination program to do more harm than good. This of course, is not to say 307 that vaccination programs should ever be avoided; only that once a vaccination 308 program is embarked upon, it should be implemented at a sufficiently vigorous 309 rate to ensure that it is not outrun by relatively fast waning immunity rates. 310 Additionally, when escape mutations emerge, generally associated with reinfec-311 tion cases [8] , then the valency of vaccines should be switched rapidly enough to 312 8 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint thwart escape mutations, rather than aiding such mutations by thwarting less 313 transmissible competing strains instead. 314 At this time, the primary value of our SEIVD-IBM J-RAMP is in testing var- outbreaks is obtaining strain-specific cross-immunity data (for characterization 327 of the elements c j of the cross immunity matrix C), waning rates (which we 328 have not made strain specific, but our model could be generalized to include The process of making our model both location and strain specific could be 338 undertaken using methods designed to enhance the relevancy of models, such 339 as Appropriate Complexity Modeling (ACM: [26, 37] ). Further, in some cases it 340 may be useful to add spatial or age-structure information to our SEIVD model 341 IBMs or include a contact network [42] , which itself may contain spatial or re-342 fined class category (e.g. age or work category) information. In addition, our 343 current implementation represents strain differences in terms of J loci with two 344 alleles (denoted by 0 and 1 respectively) at each locus. A more realistic repre-345 sentation of the genetic basis of strain differences may involve a more involved 346 genetic representation in which several alleles are possible at each locus.

Although cross-immunity and immune waning are entangled in our immu-348 nity modifier functions (i.e., φ ij ; see Eq. A.11), cross-neutralization data can 349 be used to estimate the cross and waning immunity parameters using appro-350 priate methods [43, 44] . Such data are becoming more widely available through All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. In the coming year, as we obtain more information on the nature of immunity 366 to SARS-CoV-2, it will become more apparent to us whether or not COVID-19 367 will settle into global endemicity [19, 47] . If this is the case, then constant vigi-368 lance and a well designed vaccination program with respect to vaccinating the 369 young and implementing booster vaccinations with appropriate strain valency 370 will become the order of the day. The J-RAMP presented here, with appropri-371 ate elaborations that will become evident through its future application, such as 372 being able to compute the best time to administer booster shots of the same or 

Our SEIVD-IBM in a nutshell 392 We constructed an individual-based model (IBM) of a susceptible-exposed-393 infectious-recovered (i.e., an SEIVD model, where removed R are split into 394 V=immune/vaccinated, and D=dead) epidemiological process [11, 12] in a ho-395 mogeneous population with a random encounter contact rate parameter κ 0 > 0. 396 Our formulation includes a host immunological waning process [6, 19] . It also 397 includes the emergence of pathogen strains due to a mutational process that (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The detailed formulation of our model and its algorithmic implementation is provided in Appendix A. In a nutshell we: 407 1. defined a set of 2 J pathogen strains (user selected value for J ranging from 408 0 to 8; pathogen index j = 0, ..., 2 J −1) with a genetic-relatedness topology 409 of a J-dimensional unit cube-i.e., each pathogen has J-loci that can take eter κ 0 that is Poisson distributed on [t, t + 1), t = 0, 1, · · · , multiplied 425 by an adaptive response function that reduces the contact rate with in-426 creasing disease prevalence, such that the κ(t) is reduced to κ 0 /2 when All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. implementation of the SEIVD-IBM described in the previous subsection.

12 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint

The characteristics we envision for a model platform to be a RAMP are:

1. RAMPs include a set of model parameters (constants) whose values can 484 be selected or specified (sometimes within a predefined range of values) 485 at simulation runtime using a switch, nob, slider, or text-entry window 486 accessed via a platform graphical interface or dashboard (e.g., see Table 1 which apply to our SEIVD-IBM J-RAMP). Table 2 ). platform. An R routine can be formulated to both manage the simulation 504 run and to retrieve and process the resulting data. The RAMP simulation 505 becomes a "virtual package" to the controlling R logic. See Appendix B. 506 We implemented our RAMP using Java and made ample use of all of the fea- The first variable that needs to be determined is the unit of time we use for 517 our simulations because all process rate parameters are scaled by its selection.

Since the time resolution of empirical COVID-19 incidence and mortality data 519 is daily, we selected are unit of time t to be days. Additionally, based on various 520 sources including a metapopulation study of COVID-19 parameter estimates 521 [53], we set the latent and infectious periods to be 4 and 7 days respectively. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint (95% confidence interval [2.69, 3.59]) though these authors suggested that the most recent studies they looked at suggested that R 0 was more likely between 530 2 and 3.

Assuming no seasonal effects (i.e., δ = 0 in Eq. 3), in our single strain runs 532 we kept R 0 ≈ 3 by setting κ = 4 (effective contacts per day, which nominally 533 involves being with 6 feet of one another at least 15 mins) thereby requiring 534 β = 0.3, because in the notation of the continuous time model introduced in 535 Eq. A.29, we have an inverse latent period value of γ ≡ 1/σ E = 1/4 and an 536 inverse infectious period value of ρ ≡ 1/σ I = 1/7 days ( In reality, the contact rate κ is not constant, but is time varying. It is influ- 

To keep things simple in applying our model to a multistrain setting, we 560 assumed that the immune effects on infector shedding and mutation and repli-561 cation within the infectee through a cross-immunity processes, can be tweaked 562 through a single constant c ∈ [0, 1) (making the elements c j of the matrix 563 C strain-dependent will obviously require considerable supporting data). One 564 scenario is to assume cross immunity only applies to nearest neighbors, i.e.:

Nearest-neighbor C: ding. For simplicity's sake, however, we assume that infectee with major strain 570 j will shed minor strains in the immediate neighborhood of j at comparative 571 14 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint rate ζ ∈ [0, 1) and be major strain-independent: i.e., we assume

Third, again in the absence of evidence to the contrary, we assume that 573 seasonal fluctuations are strain independent. In this case, for strain-independent 574 constants δ ∈ [0, 1) and θ ∈ [−π, π] that set seasonal fluctuations related to 575 environmental persistence, Eq. A.13 in the Appendix A simplifies for a single 576 constantη to: (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Figure 1 : An overview of the processes included in our SEIVD-IBM model. The probability π inf ih,jl of A h being infected primarily with pathogen in terms of receiving an effective dose from agent Ai, (η) is computed in terms of a concatenation of shedding rates (ζ i ), environmental persistence rates (ε ), and host transmission (β h ) processes (Eq. A.15) and includes both waning and cross immunity factors. The probability π inv h that the dominant strain emerging in host A h is strain given initial infection with strain is computed in terms of within-host mutation and within-host replication process (Eq. A.16) and also includes both waning and cross immunity factors. These two probabilities are then used to compute the overall probability π ih,j (Eq. A.17) that infector i, infected with major strain j, infects infectee h with major strain . The quantity R eff (t ) is the expected number of individuals each infectious agent is expected to infect around time t ∈ [t + σE, t + σE + σI], where R0 = R eff (0) is estimated for our model using Eq. A.29. Equation references for the other mathematical functions provided in this figure are given in Table 2. 21 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Table 1 ). The top left window of this dashboard contains information on the final state of the population (in this case S = 2135 and D = 165 in a population of N0 = 10, 000), the bottom left bar graph of dashboard panel is the final values of E, I, V and D at epidemic cessation at time t = 166 (days) or the simulation run time, whichever comes first. Dashboard also shows a graph of incidence and cumulative deaths (purple and black: selected using colored buttons below the graph). The bottom ribbon of the dash board has a series of radio buttons that respectively open a Log, a JavaScript (JS), and a Scripting (S) window, Line and Bar graph windows (for multistrain runs), as well as windows for controlling vaccination strategies (V), listing realtime agent information (A), pathogen parameter values (P), monitoring probability computations (Intern), coding and controlling runtime alternative operations (Op), and three runtime buttons (Reset, Step, Run). B. Graphs of prevalence (cut out from main panel when only the red button is on) and C. daily deaths (crimson button) are pasted below the dashboard.

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Table 1 with N = 10, 000. The waning immunity half-life values are t half = 90, 120 or 180, as labeled. In the top two panels the contact rate is constant, while in the bottom two panels it is adaptive with p half I = 0.05 or 0.01 as labeled (Eq. A.10). Note that the adaptive case approaches the constant case as p half I → ∞ and that the vertical axis in the four panels each have different scale, though the horizontal time axes is the same in all four case (365 days). Table 1 with N = 100, 000, p half I = 0.002, and t half v = 365) for the cases where vaccination rates v(t) (indicated by blue lines) are applied during the second and third years only to individuals not previously vaccinated selected at random (A. v(t) = 0.001, B. v(t) = 0.002) and to individuals not previously vaccinated and who also have not previously been infected with SARS-CoV-2 (C. v(t) = 0.002).

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint Figure 5 : The total daily incidence (∆I+: purple) and strain-specific prevalence (I: red) individuals for a 16-strain epidemic burning through a population of size N = 100, 000 are plotted over a 3-year period for the cases where cross-immunity is absent (A. c = 0.0) or intermediate among nearest neighbors differing by one allele (B. c = 0.5). The remaining parameters are listed in Table 1 and for all strains include the same host transmission (β = 0.3) and virulence (pα = 0.02), waning half-life t half = 365, minor strain shedding ζ = 0.001 and within host mutation (µ = 0.001) parameter values. Note that each panel has its own vertical scale and that the horizontal scales are set by the J-RAMP and are 0-1100 days for the incidence curves and 0-1250 days for the prevalence curves, although all curves are drawn to 1100 days. Dominant or codominant strains are labeled in red (number and equivalane allele representation), while minor strains are labeled in light grey.

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint Figure 6 : Daily incidence (deep purple curves) and prevalence of individual strains (red curves; with strain number labels plotted in order of emergence over time) that exceed a maximum prevalence of 12 individuals at any time are plotted over a three year period (t ∈ [0, 1100] days), firstly, for three cases (A.-C. in a system with nearestneighbor only cross-immunity (labeled "simple", see Eq. 1) and, secondly, in two cases (D. and F.) in a system with cascading cross-immunity (Eq. A.7). In the first three cases we have: A., no vaccination (blue labels); an B. & C., a vaccination rate of v = 0.002 applied in years 2 and 3 (t ∈ [365, 110]) to individuals who have not previously been vaccinated, respectively for the cases of univalent vaccines (strain 0, green labels and strain 127, brown labels). In the fourth case, D., a quadravalent vaccine (strain 0, 7, 31 & 127; plum labels) is applied at a vaccination rate of v = 0.01 (i.e., 1%) at random, irrespective whether or not individuals selected have been previously vaccinated (or ill with COVID-19). In the fifth case, E ("robust immunity case"), a bivalent vaccine (strains 0 & 127; orange labels) is applied, as in cases B. & C., at a vaccination rate of v = 0.002 (i.e., 1%) at random to individuals selected that have been previously vaccinated. In this latter case, however immune waning has been reduced by a factor of 5 (t half = 1825 days) with cascading cross immunity increased from c = 0.5 to c = 0.9. The dashboards for all these cases, and miniature prevalence plots of all 128 strains for each of these four cases can be found in Figs. C.5-C.9 in the Appendix. The parameter values are the same as in Fig. 5B , apart from the immune parameters being greatly strengthen in the "robust immunity case" E., and in all cases strain specific transmission parameter values βj (j = 0, ..., 127) that were randomly assigned values between 0 and 0.3, except for the increasingly dominant strain sequence β0 = 0.3 (0 ≡ (0, 0, 0, 0, 0, 0, 0)), β1 = 0.325 (1 ≡ (0, 0, 0, 0, 0, 0, 1) ), β3 = 0.35 (≡ (0, 0, 0, 0, 0, 1, 1) ),..., β127 = 0.475 (127 ≡ (1, 1, 1, 1, 1, 1, 1) ) (particular values of βj are provided for the different strains depicted here; see Table in Fig. C.4 for all values of βj, j = 0, ..., 127). Note that the vertical axis are all on the same scale height for incidence, but are individually scaled for prevalence, while the time axis are all on the same scale for the prevalence curves which differs from the incidence curves, which among themselves are all on the same time scale.

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint Figure 7 : The daily incidences (vertical axis: number of individuals) are plotted over time (horizontal axis: number of days; plots have been squashed to ensure they are all on the same time scale) for two cases: A. adaptive vaccination for nearest-neighbor only cross immunity (Eq. 1) and B. adaptive vaccination for cascading cross immunity case (Eq. A.7), and should be compared with the non-vaccination case depicted in Fig. 6A and the non-adaptive vaccination case depicted in Fig. 6D . As in this latter non-adaptive case, vaccination is performed on non-infectious individuals, selected at random (i.e., re-vaccinations can occur) at a rate of 1% per day, but now the dominant strain is targeted, as well as a second strain if its prevalence exceeds 10 cases, with vaccine valency switches every 90 days. The realized adaptive vaccine valencies for these two runs are indicated by colored bars on the horizontal time axis (red, valency is (0,127); brown, valency is (127); orange, valency is (23); green, valency is (102)).

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. See Eq. A.6: note w(t) switches from 1 to 0 as immunity goes from complete to absent # This is "virulence" parameter of continuous-time SEIR models * * If α << 1 then pα j = 1 − e −α ≈ α * † See Eq. A.10. Mostly 0.002 is used. Setting p half I = 0 implements κ(t) = κ 0 , though κ(t) → κ 0 as p half I → ∞ * + Implies values of k and θ in Eq. A.13 are irrelevant * ‡ Strain independent-strain dependence requires more elaborate model * # Quantifies the pass-on rate of mutations rather than the mutational rate of a locus or gene * * * If α j << 1 then pα j = 1 − e −α j ≈ α j 27 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. probability is major strain when invades Eq. A.16 π ih,j (t) probability is major strain in A h Eq. A.17 when j is major strain in Ai 28 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The single infected individual at time zero will be designated Agent 1 (also 39 known as patient zero and denoted by A 1 ). Thus at time t it follows that

We note that individuals in set A(t) can be in a disease state E or I with 41 respect to pathogen j, but simultaneously can be in multiple immune states if 42 29 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint they have been infected with more than one pathogen strain in the past. We also 43 note that the distinction between symptomatic and asymptomatic individuals 44 in state I will not be considered here; and only need be incorporated if testing, 45 quarantining, and treatment processes are included in the model. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint 3. I j (t, τ ij ): An agent A i who is currently infectious with strain j, after being infected with strain j at time τ ij (this is the infectious stage that lasts for 84 σ I units of time). All these individuals belong to set I(t) ⊆ A(t) 85 4. V j (t, τ ij ): An agent A i who was infectious with strain j, having been 86 infected with strain j at time τ ij , but is now non-infectious with regard to 87 this strain-that is, recovered with some immunity to strain j, as well as 88 some cross immunity to strains closely related to j. All these individuals 89 belong to set A(t) Since an agent A i may be infected over time by more than one strain j, its • If A i (t) is infected, but not yet infectious, with pathogen strain j at time t but has not been infected with any other pathogen in its past history, then A i (t) = {∅, · · · , ∅, E j (t, τ j ), ∅, · · · , ∅}

• On the other hand if A i recovered from an infection with pathogen 0 at time τ 0 , and is now in infectious with pathogen j at time t, having become infected with this pathogen at time τ j then we write A i (t) = {V 0 (t, τ 0 ), ∅, · · · , ∅, I j (t, τ j ), ∅, · · · , ∅} As we shall see, an agent history may contain at most one instance of either E j 100 or I j , while possibly containing multiple instances of V j . (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint for newly infected susceptibles and indices removed for individuals that died, then by definition:

where the number of indices in the updated set I A (t + 1) is N A (t + 1) and the 109 updated number of dead is N D (t + 1) at time t + 1.

For mathematical convenience all susceptibles S will also be referred to as 111 A 0 :, i.e., there are N S (t) individuals referenced by A 0 at time t. It will be useful 112 to partition the set A(t) itself into three subsets at time t by identifying the sets 113 E(t) and I(t) which respectively contain all agents that are currently in a state 114 E j (t) or a state I j (t) at time t for some j = 0, ..., 2 J −1. We note the intersection 115 of these two sets is empty-i.e., E(t) ∩ I(t) = ∅-as will become apparent below 116 from the transmission process rules set up below. We will use the notation

to denote the set of agents in A(t) but not in E(t) or I(t). 118 We also identify the set of infectious agents with infectious strain j. If A j denotes an agent whose epidemiological state contains an entry I j (t, ), then

where the number of such agents is denoted by N Ij (t), and its index set by (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint consideration.

Waning immunity. Recall that we use A 0 to denote an anonymous (generic) 141 member of S and that A i for i > 0 refers to a specific individual with an 142 associated state list/vector. If some specific A i is in immune state V j having 143 been infected with this strain at time τ ij , we assume that the level of relative 144 susceptibility of agent A i to reinfection by strain j is given by (noting that the 145 existence of the value τ ij implies that infection of individual i by strain j at 146 time τ ij ensures that the V j is no longer "null")

We note the following: 1.) the first case implies that τ ij has yet to be defined; 2.) 148 the second case is equivalent to the statement that τ ij ≥ 0 now exists for strain 149 j, since this occurs at time t = τ ij (through the invocation of state E j (t, τ ij )); In addition to this, one may also designate certain changes in certain alleles as 163 "escape mutants" with respect to the progenitor strain. For these, the level of (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. In deriving a probability π inf ih,j of an agent A h being infected with strain 204 by and agent A i who is infectious with major strain j, we concatenate (i.e., 205 multiply together) several process, each of which involves nominal constants.

Thus, in all but one of these processes, the scaling of these constants can be 207 normalized and given a relative set of values across strains though one set of 208 34 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint constants though relative, will ultimately all be scaled by fitting the model to 209 real data. In our treatment below, constants associated with shedding and per-210 sistence will be scaled while those associated with within-host replication will be 211 kept unscaled to be ultimately fitted to data. In particular, the parametersβ j 212 associated with pathogen transmission (i.e., from contact to the start of within 213 host replication-see Fig. 1 ) will be scaled by fitting to epidemiological data, Pathogen shedding. We assume that shedding is affected by the immune state 219 of the infector A i and thus posit the shedding rates below for this individual 220 when its major infectious strain is j. In general, we have a matrix of shedding 221 ratesζ j before accounting for immunity and cross immunity that is specific to 222 agent A i . Immunity and cross-immunity act to reduce shedding rates through (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. 

If after receiving an initial infectious dose of pathogen, an individual is infected 258 primarily with strain , then we expect this strain to dominate unless intrinsic 259 mutational processes are high (which is not the case for COVID-19) or the 260 individual has some immunity to this dominant strain. In the latter case the 261 situation is ripe for an "escape mutation," that is one that evades the immune 262 system, to arise. Probability that infector A i with major strain j will result in infectee A h express 274 as its major strain is

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint

In the single strain case (J = 0), the waning immunity equation Eq. A.6 reduces 277 to (dropping the redundant index j = 0, and noting that the existence of a value 278 τ i implies A i has been infected at time τ i in the past)

if A i has never been infected

if t < τ i + σ I + σ E (A.18) and the modifying immunity functions φ ij (Eq. A.11) collapse to 1, which implies 280 that the pathogen shedding functionsζ i (Eq. A.12) collapse to 1. Without loss 281 of generality, we can also assume that single strain value η = 1 in Eq. A.13, 282 which implies that the probability of infection (Eq. A.15) reduces to

Further, since in the single strain case there are no mutations to consider, it 284 follows from Eq. A.16 that π inv h (t) = 1 for all h and we finally have that π ih (t) = (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. Of these, proportions

are expected to come from susceptibles in the sets S(t) and A S (t) 

We note that onlyN S i (t) andN AS i (t) are of interest because indi-333 vidual in states E and I cannot be reinfected. Also, we make the 334 assumption below that the first infection that an individual in set A 335 contracts in this contact loop, is the one that counts (i.e., there will 336 38 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. 

(c) Identify all infectious agents and their pathogens strains. Among all agents in the set A(t) (Eq. A.4), identify those that have an infectious strain I j for some j = 0, · · · , 2 J − 1. Thus, if the number of infectious agents with infectious strain j is N Ij (t) then consider the set

Initially, most of these sets will be empty, but will fill in over time. nomial drawing, we can now generate the number of newly exposed 357 individuals, N E+ 0 (t + 1) (the "+" is used to denote these are newly 358 added and the "0" that they are coming from the set S), with major 359 strain at time t + 1, have been infected by agent A i with major 360 pathogen strain j on the time interval [t, t + 1):

Multinomial N S * i (t); π i0,j0 (t), · · · , π i0,j 2 J −1 (t)

These individuals will be used to update list of currently infected 362 individuals in the sets A Ej , j = 0, ..., 2 J − 1 at time t+1, which is 363 computed in the outer loop computation, as presented below. We also 364 note that the probabilities in the above multinomial add to less than (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. A h ∈ A Ej for some j ∼ Multinomial 1; π ih,j0 (t), · · · , π ih,j 2 J −1 (t) (A.25)

We note here that since the agents A h , h ∈ I A\(E∪I) are drawn with (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint (a) Individuals in A S at time t. For the N S (t) individuals in A S at time t, we have N E+ 0j (t) enter set A Ej (t + 1) and we update

where Eq. A.23 ensures that N S (t + 1) ≥ 0 402 (b) Individuals in A S that are infected again over [t, t + 1). These in-403 dividuals can become reinfected as calculated in the contact loop.

Those that become reinfected with strain j, j = 0, ..., 2 J − 1 enter (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint Algorithm 1: Summary of simulation algorithm

Vaccination process Update N S and N A for agent in I do // 4) contact loop 4a) Numbers in various sets and associated index sets 4b) Infectious contacts with each group 4c) Identify all infectious agents and their pathogens strains 4d) Susceptible contacts 4e) Agent contacts Store updates for agent in E do // 5a) disease progression loop 5a) Individuals in A E at time t Store updates for agent in I do // 5b) disease progression loop 5b) Individuals in A I at time t Store updates Updates from loops 4), 5a) and 5b) 6a) Individuals in A S at time t // 6) updates in outer loop 6b) Individuals in A S that are infected again over [t, t + 1) 6c) Updating the immunity of individuals in A S 6d) Transfer from S to A 6e) Specify the valency of the vaccination

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Javascript or the R statistical platform. The API can also be accessed by an 443 onboard scripting interface that uses the Nashorn Javascript engine.

Additionally, using a novel design, elements of the internal algorithm are (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint Figure B .1: RAM frame shows the implementation of the nearest-neighbor crossimmunity formulation C of Eq. 1 as the default. In the red-bordered insets are the cascading cross-immunity with escape mutation formulations of Eqs. A.7 and A.9. Note the radio-buttons are set to 0 in the main figure and to 1 and 2 in the insets. Also note the "+" button which allows for an unlimited number of alternatives to be set up using consecutive integer numbers for the new radio-buttons that appear and pertain to the selection of each alternative. Note at the bottom the list of functions that can be altered at runtime. The "load default" button on the upper left hand side allows the user, when starting a new alternative, to insert the default code (which is immutable in radio-button 0) as a starting point. The frame also documents a list of terms in the upper panel that can be used to build any function. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint matrix function, as described by Eq. 1. Clicking the "+" button produced two new options, which appear in the insets. These option implement the alternative 484 cascading cross immunity schemes presented in Eqs. A.7 and A.9. Option 2 code is also more extensive documentation in a separate user guide (see Fig. B .2).

The platform duplicates a mini-development environment for building alter-497 native definitions. Once code has been entered the "Compile" button checks 498 the legality of the code and makes it available for use at runtime. Legally com-499 piled code will produce a "Compilation Successful" message. Errors will appear 500 with line numbers if they occur. Once the code is legal, the "Test" button can 501 be used with actual parameters entered into the small text fields to determine 502 correctness of the code. It is also possible to include print and println state- Instructions are comprised of opcodes (e.g., reset, step, get) followed by 0 524 or more arguments. Every BPL operation returns a result, even if empty, for 525 1 a bytecode is computer source code that is processed immediately by a program, usually referred to as an interpreter or virtual machine.

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Parameter assignment and retrieval. Every user-configurable element (in-530 cluding random number generator seeds) is addressed from BPL using a unique 531 three-letter "airport code" (see Table. B1). Additionally, pathogens are ad-532 dressed by their id number (0 to 2 J − 1) and agent states using identifiers S, E, 533 I, V, DI+ and DD (the latter two represent ∆I and ∆D, respectively). RAM 534 options are addressed in setOption and getOption using the name of the RAM 535 (e.g., "crossImmune"). Get and set operations can be used on each of these with 536 the exception of ENT (strain entropy), which is read-only. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Simulator operation. Simulation runs begin by executing the BPL reset instruction, followed by step, run for or run. The BPL interpreter operates 540 synchronously with the simulator by waiting to process subsequent commands 541 during a simulation run. Operational instructions can be interspersed with pa-542 rameter set/get or data retrieval to use in runtime decision-making. Note that 543 the reset operation restores the simulator to its state at the time of the last 544 reset, so that no parameter changes made during a run are persistent. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint Figure B .3: The list of Blackbox Programming Language (BPL) commands that can be used to write a simulation driver script, using the three-letter "airport codes" listed in Table A1 to access the parameters and variables in our coded algorithm. This list of commands can be accessed using the "Command Reference" button at the bottom of the Scripting Window.

As previously mentioned, the API supports remote control of the simulator 570 from independent platforms using the operating system's indigenous fifo chan-571 nels. Of particular interest is integration with the R statistical programming 572 48 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint Figure B .4: The JS scripting window accessed by selecting the "JS On" button in the main dashboard (see button second from left at bottom of Fig. 2A) . The script shown here was used to execute the adaptive vaccination strategy discussed in the main text.

environment. An R-package called "seiv" acts as a driver by synchronously is-573 suing BPL command strings and waiting for results. Consequently, a simulation 574 can be driven entirely from within the R platform, treating the simulator as a 575 "virtual package". (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint Figure B .5: Our SEIVD-IBM can be treated as an R-package called "seiv" and run as such in conjunction with other packages, such as ggplot2 and reshape to conduct multiple simulations and then carry out data and statistical analyses of the simulation results. The code shown here was used to produce the plots illustrated in Fig. C.2. 50 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint individuals) and dead (D(t fnl ) individuals) classes and the length of the epidemic (t fnl days) are provided here, together with summary statistics, for 10 simulations (runtime seed ranges from 0 to 9) in each case when the initial population sizes are N0 = 10, 000 and N0 = 30000 individuals respectively. The parameter values used are otherwise the same as those use to produce the simulation depicted in Fig. 2 in the main text.

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint Figure C.2: Plots of the mean (red) ± 1 standard deviation green and blue, the latter truncated at 0) of prevalence, incidence, and daily deaths for 100 runs (runtime seeds = 0,...,99) using the same parameter set used to produce the individual run (runtime seed = 0) depicted in Fig. 2 in the main text (also see Fig. C.1) . These figures where produced by running our SEIVD-IBM J-RAMP as an R-package as described in Appendix B, with the code documented in Fig. B .5.

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint (Table 1) , we set p half I = 0.005 then the total number of individuals that have been infected (disease class V ) at the end of one year is 10027/100, 000 ≈ 10%.

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint 54 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint Table 1 ). See Fig. 6 in main text for enlarged plots of predominant strains.

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint Table 1 ). See Fig. 6 in main text for enlarged plots of predominant strains.

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint Table 1 ). See Fig. 6 in main text for enlarged plots of predominant strains.

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint Table 1 ), except that the cross immunity matrix C is now defined by Eq. A.7. See Fig. 6 in main text for enlarged plots of predominant strains.

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint Table 1 ), except that the half-wining values has been increased fivefold t half = 1825 days the cross immunity matrix C is now defined by Eq. A.7 with c = 0.9. See Fig. 6 in main text for enlarged plots of the predominant strain 0. 59 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 11, 2021. ; https://doi.org/10.1101/2021.06.07.21258504 doi: medRxiv preprint

The impact of vaccination on the epidemiology of infectious