A modeling-based approach to optimize COVID-19 vaccine dosing schedules for improved protection =============================================================================================== * Prashant Dogra * Carmine Schiavone * Zhihui Wang * Javier Ruiz-Ramírez * Sergio Caserta * Daniela I. Staquicini * Christopher Markosian * Jin Wang * H. Dirk Sostman * Renata Pasqualini * Wadih Arap * Vittorio Cristini ## Abstract While the development of different vaccines has slowed the dissemination of SARS-CoV-2, the occurrence of breakthrough infections continues to fuel the pandemic. As a strategy to secure at least partial protection, with a single dose of a given COVID-19 vaccine to maximum possible fraction of the population, *delayed* administration of subsequent doses (or boosters) has been implemented in many countries. However, waning immunity and emergence of new variants of SARS-CoV-2 suggest that such measures may jeopardize the attainment of herd immunity due to intermittent lapses in protection. Optimizing vaccine dosing schedules could thus make the difference between periodic occurrence of breakthrough infections or effective control of the pandemic. To this end, we have developed a mechanistic mathematical model of adaptive immune response to vaccines and demonstrated its applicability to COVID-19 mRNA vaccines as a proof-of-concept for future outbreaks. The model was thoroughly calibrated against multiple clinical datasets involving immune response to SARS-CoV-2 infection and mRNA vaccines in healthy and immunocompromised subjects (cancer patients undergoing therapy); the model showed robust clinical validation by accurately predicting neutralizing antibody kinetics, a correlate of vaccine-induced protection, in response to multiple doses of mRNA vaccines. Importantly, we estimated population vulnerability to breakthrough infections and predicted tailored vaccination dosing schedules to maximize protection and thus minimize breakthrough infections, based on the immune status of a sub-population. We have identified a critical waiting window for cancer patients (or, immunocompromised subjects) to allow recovery of the immune system (particularly CD4+ T-cells) for effective differentiation of B-cells to produce neutralizing antibodies and thus achieve optimal vaccine efficacy against variants of concern, especially between the first and second doses. Also, we have obtained optimized dosing schedules for subsequent doses in healthy and immunocompromised subjects, which vary from the CDC-recommended schedules, to minimize breakthrough infections. The developed modeling tool is based on generalized adaptive immune response to antigens and can thus be leveraged to guide vaccine dosing schedules during future outbreaks. Keywords * booster * COVID-19 * SARS-CoV-2 * cancer * immunocompromised * mathematical modeling * breakthrough infection * omicron * vaccines * variants of concern ## 1. Introduction Since December 2019, the COVID-19 pandemic caused by SARS-CoV-2 has afflicted more than 655 million individuals and caused more than 6.67 million deaths worldwide [1]. Global vaccination programs along with public health measures such as social distancing and masking are anticipated to be the most effective approaches to attain herd immunity and curb the pandemic [2, 3]. Herd immunity represents a scenario where a virus cannot spread due to a dearth of susceptible hosts and can be achieved through natural infection and/or vaccination of the population. In December 2020, the first COVID-19 vaccine obtained Emergency Use Authorization from the United States Food and Drug Administration (FDA), and as of December 2022, 50 vaccines have obtained regulatory approval in at least one country [4]. As a result, over 63.4% of the world population is fully vaccinated and ∼69% of the population has received at least a single dose of a COVID-19 vaccine. However, due to the inequitable allocation of vaccines, only ∼26% of the people in low-income countries have received at least a single dose [5, 6], which cans facilitate the emergence of new variants of SARS-CoV-2 and thus resurgence of the pandemic. According to a meta-analysis, seroconversion rates related to the development of neutralizing antibodies in the sera of individuals doubly vaccinated with COVID-19 vaccines have been found to be dependent on patient immunological health status; seroconversion positivity in immunocompetent individuals can be up to 99%, while in immunosuppressed patients the efficacy of vaccination varies for different diseases (e.g., solid tumors ∼92%, immune-mediated inflammatory diseases ∼78%, hematological cancers ∼64%, and organ transplant recipients ∼27%) [7, 8]. Due to limited protection, immunocompromised individuals are more vulnerable to infection and are at a higher risk of developing severe or lethal COVID-19. Thus, immunizing the majority of the population is a means to additionally protect individuals who are susceptible or unable to receive a vaccine. However, the emergence of breakthrough infections in previously infected or vaccinated individuals is a major challenge. The key biological reasons for breakthrough infections are attributed to: (i) waning immunity over time, and (ii) emergence of mutant variants of SARS-CoV-2, referred to as variants of concern (VOCs) [9, 10]. Depending on demographics and the type of vaccine administered, the humoral response (i.e., neutralizing antibodies) against SARS-CoV-2 has been found to be substantially reduced within about six months after two-dose vaccination [11-13]. Thus, vaccines with an initial effectiveness of 90% are only ∼30-70% effective after six months [14-16]. Further, coronaviruses tend to have high genetic diversity due to their large genome size (26.4 – 31.7 kb), high mutation rate caused by a low-fidelity viral polymerase (∼10−4 substitutions per site per year), and high recombination frequency (up to 25% for the entire genome *in vivo*) [17]. As a result of selection pressure imposed by neutralizing antibodies on viral surface proteins, particularly the receptor binding domain (RBD) and the N-terminal domain (NTD) of the spike protein, which are the targets of most of the COVID-19 vaccine-induced neutralizing antibodies, SARS-CoV-2 show clusters of mutations as documented in the genomes of VOCs [18]. Mutations that confer greater fitness such as increased transmission rates and improved antibody escape are positively selected, leading to antigenic drift that makes the vaccination-induced neutralizing antibodies partially ineffective against the mutant strains [17]. This predisposes the vaccinated or previously infected individuals to breakthrough infections [19] (though the severity of symptoms tends to be milder) [20]. Currently, additional (booster) doses of COVID-19 vaccines are being used to reinforce protection and minimize breakthrough infections [21-24]. Boosters are being administered to fully vaccinated individuals since ∼June 2021, except in low-income countries [25], and prioritized for high-risk populations such as the elderly and immunocompromised patients [26]. According to the Centers for Disease Control and Prevention (CDC), a two-dose schedule (3- to 8-week gap) followed by a third dose (5-month gap) of mRNA vaccine (Pfizer-BioNTech or Moderna) is recommended for immunocompetent adults, while a three-dose schedule (3- to 4-week gap between doses 1, 2, and 3) followed by a fourth dose (12-week gap) is recommended for immunocompromised adults [27]. These scheduling recommendations are based on clinical trials, which are generally limited to healthy volunteers, thereby may require optimization, especially for special populations, to achieve better protection at the population scale. A mathematical modeling approach, which is data-driven and based on first principles of physiology, immunology, and biophysics can be a valuable tool to simulate population-scale heterogeneity in immune health status and immune response to vaccines, thereby supporting rational design of dosing schedules. In addition, given the disparities in global vaccine allocation, optimization of dosing schedule to extend the gaps between doses with no major effect on efficacy could allow for improved distribution of vaccines to countries without the capacity to provide for themselves, reduce costs, and promote vaccine compliance, thereby benefiting the overall population, but especially patients in critical care. Using a mathematical modeling approach, we identified optimal vaccine dosing schedules of mRNA COVID-19 vaccines for immunocompetent and immunocompromised individuals to minimize breakthrough infections at the population scale. Clinical evidence that demonstrates acceptable vaccine effectiveness, despite delayed follow-up doses, sets the premise for our theoretical investigation [28-30]. Previous mathematical models that have been developed to identify optimal vaccine allocation and dosing schedules to minimize hospitalizations and deaths due to COVID-19 are primarily age-structured compartmental models, based on epidemiological principles (e.g., susceptible, exposed, infectious, and removed (SEIR) models), which focus on the transmission of the virus under different vaccination scenarios and the analysis of strategies to reduce the rate of infection [31-37]. These models, however, lack mechanistic details relevant to virus-host interaction, the immune response to vaccines, and the time-dependent variation in vaccine efficacy due to inter-individual variability, vaccine efficacy against VOCs, and other biological/physiological factors. To this end, as an adaptation of our previous mechanistic models of complex biological systems [38-43], we have developed a mathematical model that accurately simulates the adaptive immune response to COVID-19 vaccines at the individual scale. The model was calibrated and validated with clinical data for mRNA-based COVID-19 vaccines to conduct analysis in virtual cohorts comprising immunocompetent and immunocompromised digital twins (virtual cancer patients undergoing chemotherapy and/or immunotherapy). The model identified optimal schedules for vaccination doses that minimize vulnerability to breakthrough infections, especially against VOCs (specifically Omicron), while retaining vaccine efficacy above the protection threshold in populations with different health statuses. ## 2. Methods ### 2.1 Model development Based on our previous mathematical modeling of the immune response to SARS-CoV-2 infection [44], we developed a model of the adaptive immune response to COVID-19 vaccines. As shown in **Figure 1**, the model incorporates key biological processes that are relevant to antigen presentation at the site of vaccination (i.e., muscle), the development of adaptive immune responses in the lymphoid tissue, and protection against infection in the respiratory tract. The model was formulated as a system of ordinary differential equations (ODEs, **Equations 1-17**), which describe the kinetics of key immune response variables following vaccination or infection. The equations were solved numerically as an initial value problem in MATLAB R2018a. While some of the model parameters were known a priori (**Table 1**), the remainder were estimated by non-linear least squares fitting of the model to multiple clinical datasets obtained from the literature [45-47]. The model was then used to simulate the immune response to mRNA-based COVID-19 vaccines in healthy and immunocompromised populations and was implemented to identify optimal vaccine dosing schedules to minimize breakthrough infections. The model equations are described in detail below. View this table: [Table 1.](http://medrxiv.org/content/early/2023/01/31/2022.09.14.22279959/T1) Table 1. List of model parameters. ![Figure 1.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/01/31/2022.09.14.22279959/F1.medium.gif) [Figure 1.](http://medrxiv.org/content/early/2023/01/31/2022.09.14.22279959/F1) Figure 1. Model schematic. Diagram shows key variables and system interactions incorporated into the mathematical model. Upon respiratory tract infection by SARS-CoV-2 or intramuscular administration of mRNA vaccines, antigen presenting cells (e.g., macrophages) engage the adaptive immune system to produce antibodies and activate T-lymphocytes to build immunity against infection. Cytokines secreted by infected cells (e.g., IFN-I) and immune cells (e.g., IFN-II, IL-6) in the process have modulatory effects on the immune system. Abbreviations: IFN-I, type-I interferon; IFN-II, type-II interferon; IL-6, interleukin 6. At the site of vaccination, nanoparticles carrying the mRNA of SARS-CoV-2 spike protein are endocytosed into myocytes, leading to the translation and expression of spike protein on myocytes [48]. Given that the timescale of drug delivery (intramuscular injection) and mRNA translation is much shorter (< 1 hour) [49] than that of the vaccine-induced immune response (days to weeks) [50], we assumed that the variable *C**a*(*t*) represents the concentration of vaccine-induced spike protein in the muscle cells that can trigger the immune response via antigen-presenting cells (APCs). #### Concentration kinetics of the exogenously administered antigen (via vaccine) in muscle cells (*C**a* (*t*)) ![Formula][1] where Dose indicates the dimensionless dose of the antigen administered via the vaccine. The concentration of the spike protein *C**a*(*t*) is described by the sum of Gaussians centered at τ*i*, which represents the day on which a vaccine dose is injected out of the set of doses indicated by *S**T*. *T*NP is the characteristic time of clearance of the antigen-carrying nanoparticle (NP) from the body [39], estimated based on NP diameter of 100 nm for mRNA vaccines [51]. The population of naïve (or immature) APCs is maintained through continuous regeneration and presumably maintained at a steady state. Thus, we used a logistic growth term to include this contribution, where *γ*APC is the exponential growth rate, and ![Graphic][2] is the carrying capacity of the APC population. Naïve APCs at the site of expression of spike proteins recognize, process, and present the antigen via major histocompatibility complex (MHC) during differentiation into activated APC (APC*) at a rate *T*APC as they migrate towards the lymphoid tissue. The APC activation process is proportional to the antigen load (Ag(*t*)), which can be derived either from the vaccine or natural infection and is either equal to *C**a*(*t*) or the viral load *V*(*t*) in the case of vaccination or infection, respectively, with *K**ν* being the Michaelis constant for antigen-induced activation of naïve APCs. #### Equation for the naïve APC density at the site of vaccination or natural infection (APC(*t*)) ![Formula][3] Activated APCs are primarily responsible for the induction of the adaptive immune response, and their population is determined by the activation of naïve APCs, which we discussed in Eq. 2, and a death term determined by the death rate constant, *δ*APC of activated APCs. #### Equation for the activated APC density (APC*(*t*)) ![Formula][4] Activated APCs migrate from the site of vaccination or natural infection to the lymphoid tissue to interact with naïve T-cells (CD8+ or CD4+) and transform them into their active or effector forms. Alternatively, naïve B-cells are activated by the binding of soluble antigens, which however in the current model is replaced by binding to active APCs, given that the density of active APCs is dependent on antigen load in the body. For the naïve cells, population density is determined by cell regeneration and cell activation, where we used a logistic growth term with *γ*CD4, *γ*CD8, and *γ*B as the growth rates of naïve forms of CD4+ T-cells, CD8+ T-cells, and B cells, respectively; ![Graphic][5], ![Graphic][6], and ![Graphic][7] are the carrying capacities of the corresponding cell the product of active APC density and the corresponding naïve cell density; *T*CD4,*T*CD8 and *T*B populations, respectively. The activation term has second-order kinetics and is proportional to are the activation rates of lymphocytes indicated by the subscript. Activation of T-cells is amplified by the presence of type-II interferons (IFN2(*t*)) secreted by activated T-cells [52], with possible saturation effects. Thus, we used a Michaelis-Menten term to model this process in which *K*IFN2 is the Michaelis constant of type-II interferon effects. Of note, in our model we have included a dimensionless coefficient *f* ∈ [0, 1] that represents an immunosuppression factor to modulate the carrying capacity (i.e., homeostasis levels) of the naïve immune cell population to model immunocompromised subjects, such that *f* = 1 in healthy individuals, and *f* < 1 in immunocompromised patients. Note that since our model is calibrated for patients who are immunocompromised due to anticancer therapy, therefore immunosuppression in our model is characterized by T- and B-cell immunodeficiency, which is one of the important immunological effects observed due to disruption of hematopoiesis leading to myelosuppression in patients undergoing anticancer therapy [53-57]. Also, in the case of naïve CD4+ and CD8+ T-cells we have included the ability of interleukin-6 (IL-6) to cause T-cell exhaustion [58] by including an additional term that limits the carrying capacity of these cells. This term uses the concentration of IL-6 in a Michaelis-Menten function, where, *K*IL6 is the Michaelis constant for IL-6 effects. #### Equation for the naïve CD4+ T-cell density (CD4(*t*)) ![Formula][8] #### Equation for the effector CD4+ T-cell density (CD4*(*t*)) ![Formula][9] where *δ**T* is the death rate of effector T-cells. #### Equation for the naïve CD8+ T-cell density (CD8(*t*)) ![Formula][10] #### Equation for the effector CD8+ T-cell density (CD8*(*t*)) ![Formula][11] #### Equation for the naïve B cell density (*B*(*t*)) ![Formula][12] where *T*B is the transition rate of naïve B cells into their activated form. Of note, the activated B cells differentiate into antibody-secreting plasma cells upon interaction with effector CD4+ T-cells. We modeled this interaction using second-order kinetics, where is *T*BC the differentiation rate of B cells into plasma cells. #### Equation for the activated B cell density (*B**(*t*)) ![Formula][13] where *T*BC is the differentiation rate of B cells into plasma cells. #### Equation for the plasma cell density (*P*(*t*)) ![Formula][14] where *δ**P* is the death rate of plasma cells. Virus-neutralizing antibodies are secreted by plasma cells, such that their rate of production, characterized by the first-order rate constant *P*Ab, is proportional to the plasma cell density. The antibodies secreted into the plasma are then cleared at a rate ClAb, which is a lumped phenomenological parameter characterizing the various antibody clearance mechanisms. #### Equation for the neutralizing antibody concentration (Ab(*t*)) ![Formula][15] Following vaccination or natural infection, the immune system produces different cytokines to regulate cellular activation and differentiation, as discussed above. In the specific case of SARS-CoV-2, it has been shown that type-I and type-II interferons, and IL-6 are the relevant immune response. For instance, type-I interferons (IFN1(*t*)), secreted by virus-infected cells or immunoregulatory elements [59, 60]. Each cytokine has a unique source and key role in the vaccine-affected cells, lowers the production of new virions by infected cells [52]; type-II interferon (IFN2(*t*)), produced by effector CD4+ and effector CD8+ T-cells, accelerates the differentiation of naive T-cells into their effector form in a positive feedback loop fashion [52]; and IL-6, secreted by effector CD4+ T-cells, effector CD8+ T-cells, and active APCs, tends to exhaust naïve CD4+ and CD8+ T-cell population [58]. The rate of change of cytokine concentration was modeled using a production term and a degradation term, where production and degradation are modeled as first-order processes, with degradation characterized by a common degradation rate constant *δ*cyt. #### Equation for the type-I interferon concentration (IFN1(*t*)) ![Formula][16] where *P*INF1 is the production rate of type-I interferons. #### Equation for the type-II interferon concentration (IFN2(*t*)) ![Formula][17] where *P*INF2 is the production rate of type-II interferons. #### Equation for the interleukin-6 concentration (IL6(*t*)) ![Formula][18] where, *P*IL6 is the production rate of IL-6. The entire immune cascade can be triggered either by a vaccine (as we have already elaborated), or through an infection caused by the SARS-CoV-2 virus. In the latter case, the infection is characterized by the transformation of healthy susceptible cells into infected cells by the virus, followed by production of new viral particles by the infected cells. With the intent to develop a generalized mathematical model capable of simulating immune response to vaccines as well as infections, we incorporate the infection process into our model, with the respiratory tract as a representative site. For this, we used a target cell limited model of acute viral infection [61], as described by the equations below: #### Equation for the healthy respiratory epithelial cell density (*H*(*t*)) ![Formula][19] where *β* is the viral infectivity rate, *V*(*t*) is the viral load density in the respiratory tract, and *H* is the initial density of healthy cells. #### Equation for the density of infected cells in the respiratory tract epithelium (*I*(*t*)) ![Formula][20] where *δ* represents the cytopathic death rate of infected cells, *δ*C is the death rate of infected cells mediated by effector CD8+ T-cells, and CD8*(*t*) is the density of effector CD8+ T-cells. #### Equation for the viral load density in the respiratory tract (*V*(*t*)) ![Formula][21] where *P**ν* represents virion production rate, IFN1(*t*) is the concentration of type-I interferons, *K*IFN1 is the Michaelis constant of the virion production suppression factor, *k*Ab is the antibody-mediated neutralization rate of viruses, Ab(*t*) is the antibody concentration in the body, APC(*t*) is the density of naïve APCs in the respiratory tract, *k*APC is the naïve APC-mediated clearance rate of viruses, and *V* is the initial viral load at the time of infection. ### 2.2 Model calibration and validation Using the built-in MATLAB function *lsqcurvefit*, non-linear least squares regression was performed to fit the model to literature-derived clinical data to estimate the unknown model parameters (**Table 1**). The datasets used for model calibration included: (i) viral load and immune response kinetics following a SARS-CoV-2 infection [45], (ii) immune response kinetics following vaccination with mRNA vaccines in healthy individuals [46], and (iii) cancer patients undergoing chemotherapy or immunotherapy [47]. Further, to test the predictive ability of our model to accurately reproduce the immune response to mRNA vaccines, we simulated two and three doses of the Pfizer-BioNTech and Moderna vaccines in healthy individuals using the parameters obtained from model calibration for healthy population (**Table 1**), and compared it to published clinical data [46, 62, 63]. Specifically, for calibration of immune response to infection, we digitized average, longitudinal viral load and immune variable data (naïve CD4+ T-cells, naïve CD8+ T-cells, effector CD4+ T-cells, effector CD8+ T-cells, type-I interferons, type-II interferons, IL-6, and neutralizing antibody titer) for moderately infected COVID-19 patients (N = 80 patients) from Lucas et al. [45]. Moderately infected patients were chosen over severely infected ones for the lack of pharmacological intervention in the former, thereby allowing the calibration of purely immune response effects on containing infection. Further, for the calibration of immune response to mRNA-based COVID-19 vaccines in healthy subjects, we extracted average, longitudinal immune response data (neutralizing antibody titer, effector CD4+ T-cells, effector CD8+ T-cells) to two doses of Pfizer-BioNTech COVID-19 mRNA vaccine (N = 31 subjects) from Collier et al. [46]. The two doses were given 21 days apart, with immune response measured 2 to 4 weeks (3 weeks average) after second dose, 6 months after first dose, and 8 months after first dose. From the same study, average immune response data (neutralizing antibody titer) for two doses of Moderna COVID-19 mRNA vaccine (21 days apart) was also extracted for model validation (N = 22 subjects, see below). For calibration of vaccine immune response in cancer patients undergoing antineoplastic treatment, longitudinal antibody titer data (N = 63 patients receiving chemotherapy, N = 16 patients receiving immunotherapy) following two doses of Pfizer-BioNTech COVID-19 mRNA vaccine was extracted from Peeters et al. [47]. The two doses were given 21 days apart and immune response was measured at the time of second dose and 7 and 28 days after second dose. Additional data for model validation (antibody titer kinetics) following two and three doses of Pfizer-BioNTech COVID-19 mRNA vaccine in healthy subjects was obtained from Bayart et al. (N = 158 subjects) [62] and Papazisis et al. (N = 110 subjects) [64], respectively. In Bayart et al., the two doses were given 21 days apart and the immune response was measured 14, 28, 42, 56, 90 and 180 days following the first dose. In Papazisis et al., two doses were given 21 days apart, following which a third dose was given 9 months after the second dose, and immune response was measured 2 weeks after the first dose, 2 weeks after the second dose, 3 months, 6 months, 9 months, and 12 months after the second dose. To account for the uncertainty in parameter estimation during model calibration, model predictions were accompanied by 90% prediction interval obtained through 10,000 simulation runs of the model, where each simulation was obtained for a unique set of parameter values drawn through Latin hypercube sampling (LHS). Note that for LHS, all parameter values were chosen from within a ±10% range of the baseline values estimated during model calibration with healthy subjects. ### 2.3 Vaccine efficacy estimation In accordance with the literature [65, 66], we used the plasma levels of neutralizing IgG anti-spike antibodies against SARS-CoV-2 as predictors of vaccine efficacy (i.e., correlate of protection against SARS-CoV-2). For this, we characterized an empirical correlation between neutralizing antibody titer (Ab(*t*)) and vaccine efficacy (*V*eff(*t*)) based on clinical data from the literature [67]. The following Michaelis-Menten function was thus used: ![Formula][22] where ![Graphic][23] is the maximum possible efficacy of antibody and *K*eff is the Michaelis constant for vaccine efficacy. As shown in **Figure S1** (solid blue curve), the above function is in excellent agreement with the clinical data, giving an estimate for *K*eff = 18.95 U/mL (i.e., 50% efficacy) and ![Graphic][24]. While only 50% efficacy is necessary to obtain approval for clinical use of vaccines [68], our analysis is based on a more stringent threshold to ensure protection in majority of recipients. According to Goldblatt et al. [67], the average plasma antibody titer for various COVID-19 vaccines to be protective against wildtype strain (WT) of SARS-CoV-2 is 154 U/mL, which corresponds to a vaccine efficacy of ∼82.3% on the Michaelis-Menten curve. Therefore, 154 U/mL was used as a threshold to differentiate protected versus non-protected individuals in our analyses. *Note that in this work we assume that individuals with antibody titer above the protection threshold are fully protected, while the ones below are fully at risk of infection* [69]. Of note, for the VOCs, the protective threshold was corrected for by using the binding score (Abescape) of the antibodies obtained from the literature [70]. For this, the previous function was modified to: ![Formula][25] where Abescape is a dimensionless binding score (Abescape ∈ [0, 1]) obtained from Greaney et al. that for the VOC [70] that quantifies antibody escape, i.e. the inability of antibodies to neutralize the virus due to inefficient binding. As per Greaney et al., the value of Abescape for WT is 1 and that for the VOC studied here (i.e., omicron (OM)) is 0.2, indicating that the strength of binding of antibodies to OM is five times lesser than that for WT. As a result, 770 U/mL (i.e., 5 times of 154 U/mL) was the estimated threshold of protection against OM, corresponding to 82.3% vaccine efficacy (based on dotted orange curve in **Figure S1**). Of note, the above calculations assume that the mutations in the RBD- or NTD-domain of SARS-CoV-2 spike protein negatively affect the binding affinity of antibodies [71], which implies that to obtain a similar protection against OM, or other VOCs, a higher antibody titer is necessary. ### 2.4 Sensitivity analysis To evaluate the relative effect of model parameters on antibody titer following an injection of mRNA-based COVID-19 vaccine, we performed global sensitivity analysis (GSA) and local sensitivity analysis (LSA) with parameters of interest. For this, model parameters were perturbed from their baseline values and the effect of parameter perturbation on model output of interest quantified (i.e., antibody titer). Firstly, to rank order the parameters for their relative importance in determining antibody titer following vaccination, GSA was performed where all model parameters were *simultaneously* perturbed over a uniformly distributed range of ±50% around the baseline parameter values (obtained from model calibration), except parameter *f* that was perturbed between the values of 0.1–1 (left half-Gaussian distribution); area under the antibody concentration kinetics curve (AUC) from zero to 15 days post injection was calculated for each simulation (i.e., for a given combination of parameter values). Note, to comprehensively investigate the vast multiparameter space (22 parameters), yet to minimize the number of simulations, Latin Hypercube Sampling (LHS) was used to obtain 10,000 combinations of parameter values, and 10 such replicates were obtained, based on our previously developed workflows. [72, 73] Multivariate linear regression analysis was then performed on every replicate, and regression coefficients were determined as a quantitative measure of parameter sensitivity index (SI). A distribution of regression coefficients (or SI) was obtained for each parameter, and one-way ANOVA with Tukey’s test was used to rank order the parameters in terms of their sensitivity, such that a higher SI represents a greater influence on model output (i.e., AUC of antibody titer). Next, to evaluate the correlation between parameter perturbations and change in antibody titer, LSA was performed, where parameters were perturbed individually at 100 linearly spaced levels over a uniformly distributed range of ± 50% around the baseline parameter values (obtained from model calibration). The corresponding change in AUC of antibody titer (from zero to 15 days post injection) with respect to parameter perturbation was calculated with the following formula for SI: ![Formula][26] where, AUC0−15*d* is the AUC of antibody titer under baseline conditions from 0-15 days following injection, ![Graphic][27] is the AUC of antibody titer under parameter perturbation condition, *P**i* is the baseline value of parameter *i*, and ![Graphic][28] is the perturbed value of parameter *i*. ### 2.5 Virtual patient cohort design To perform population-scale numerical experiments, two types of patient populations were generated, namely Cohort A and Cohort B, as described below. #### Cohort A A virtual cohort of 10,000 individuals was generated using Latin hypercube sampling (LHS) [74-76] from twelve parameter distributions (**Figure S2**), such that each individual of the cohort varied in terms of their immune health status defined by *f*, underlying biology (characterized by the higher-ranking parameters of GSA, i.e., top 10 parameters), and vaccination schedule. The chosen range for the parameter values was such that the *f* parameter varied between 0.5 and 1 (left half-Gaussian distribution with mode equal to 1 and standard deviation equal to 15% of the mode, **Figure S2c**), while the other biological parameters were normally distributed with a mean equal to the baseline parameter value and one standard deviation equal to 5% of the mean value (**Figure S2d-l**); the dosing schedules varied between two weeks and eight weeks for the second dose (continuous uniform distribution, **Figure S2a**), and between five months and nine months for the first booster dose (i.e., third dose; continuous uniform distribution, **Figure S2b**). Note that to reflect the distribution of immune health status and biological variability in a population realistically, i.e., *majority* of the population comprising healthy individuals, Gaussian distributions with limited variance were chosen over unform distributions for LHS. #### Cohort B Three virtual sub-cohorts of 10,000 individuals each to represent healthy, mildly immunocompromised, and highly immunocompromised individuals were generated through LHS. The range of *f* values used to represent immune health status was *f* = 0.9 to 1 for healthy (continuous uniform distribution), *f* = 0.7 to 0.9 for mildly immunocompromised (continuous uniform distribution), and *f* = 0.5 to 0.7 for highly immunocompromised individuals (continuous uniform distribution). Also, biological variability was included in each sub-cohort by LHS of the relevant biological parameters (identified through GSA, same as used for cohort A), assumed to be normally distributed with a mean equal to the baseline parameter value and one standard deviation equal to 5% of the mean value. For each sub-cohort, we tested 100 dosing schedules ranging from two to eight weeks (after the first dose) for the second dose (continuous uniform distribution), 0.5 to nine months (after the second dose) for the first booster (i.e., third dose) (continuous uniform distribution), and one to nine months (after the first booster) for the second booster (i.e., fourth dose; continuous uniform distribution) for their effect on continuity of protection against OM (i.e., vulnerability to breakthrough infection due to mutants). ### 2.6 Vulnerability kinetics and vaccine dosing schedule optimization To study the temporal evolution and quantify the vulnerability to breakthrough infections at the population scale, we calculated a vulnerability kinetics curve in our numerical experiments (as shown in **Figure 6c**). From the vaccine efficacy calculation (based on **Figure S1**), on a given day, the fraction of simulated individuals below the protection threshold for OM or WT (i.e., <82.3% efficacy) was calculated to obtain the population fraction that is at a high risk of infection. Performing this calculation daily for the entire simulation period gave us the curve shown in **Figure 6c**, referred to as the vulnerability kinetics curve. Subsequently, we calculated the area under the curve (AUC) as a measure of total vulnerability to breakthrough infections, which was then used as a metric for optimizing dosing schedules to impart prolonged protection against OM, as discussed below. Note that we do not use epidemiological principles to model the spread of infection or risk of exposure amongst the individuals of virtual cohorts, given that the goal of our work is to only evaluate the effect of dosing schedule optimization in maintaining protection in already exposed individuals. To optimize the timing of the second dose, immune response kinetics for each virtual individual (Cohort B) was simulated for up to 150 days after the first dose (given on day 0). From the corresponding antibody concentration kinetics, the vaccine efficacy kinetics for OM were computed using **Eq. 19**. Subsequently, we estimated the vulnerability to breakthrough infections over time. From the vulnerability kinetics plot, the area under the curve (AUC0-150d) was calculated using the trapezoidal method. After calculating the AUC0-150d for 100 dosing schedules (ranging from two to eight weeks for the three cohorts), we identified the schedules that led to a minimum in the three sub-populations, which translates to a minimized vulnerability to breakthrough infections resulting from OM. Next, using the optima found in the previous step, we repeated the process to identify the optimal timing for the first booster (third dose) in the three cohorts. In this case, the total simulated time was 500 days. Thus, the AUC0-500 d was calculated from the breakthrough infection vulnerability kinetics plots to identify the minima. Finally, using the optimal dosing schedules for the second dose and first booster (third dose), we estimated the optimal timing for the second booster (fourth dose) using the same process as described before. In this case, the total simulation time was 700 days. ## 3. Results and Discussion ### 3.1 Model calibration The focus of this work was to mechanistically model the individual-scale immune response to COVID-19 vaccines and apply it to optimize vaccine dosing schedules to maximize protection against SARS-CoV-2 and thus minimize breakthrough infections in the population. For this purpose, we began by fitting the model to immune response kinetics of SARS-CoV-2 infection [45], which allowed us to estimate several unknown model parameters relevant to key immune response variables that were otherwise difficult to compute from vaccination data alone (**Table 1**). This enabled the reliable simulation of immune response kinetics following infection. As shown in **Figure 2a**, the numerical solutions of the model are in agreement with the clinical data for viral load and immune response kinetics following SARS-CoV-2 infection [45]. This is also indicated by the strong Pearson correlation between the observations and the model fits (**Figure S3a**; *R* > 0.99). The computed kinetics of viral load in the respiratory tract predict an incubation period of eight to nine days, which is in accordance with values established in the literature [77]. Moreover, the simulations closely approximated the kinetics for eight additional cellular and molecular immune response variables, including naïve and effector lymphocytes, antibodies, interferons, and interleukins. This suggests that the model predictions are within physiological limits and thus the estimated parameter values are reliable. The results also showed that the viral load peaks around day 10, reaching a level of ∼107GE·mL−1 while adaptive immunity variables (lymphocytes, neutralizing antibodies) peaked at around day 15, which led to clearing of the infection within five weeks without any pharmacological intervention. ![Figure 2.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/01/31/2022.09.14.22279959/F2.medium.gif) [Figure 2.](http://medrxiv.org/content/early/2023/01/31/2022.09.14.22279959/F2) Figure 2. Model calibration. Model calibration with literature-derived clinical data of immune system response kinetics during **a**) SARS-CoV-2 infection in moderately infected subjects, vaccination in **b**) healthy individuals, **c**) cancer patients receiving chemotherapy, and **d**) cancer patients receiving immunotherapy. For consistency, all immunization data was based on two simulations; markers with error bars represent mean ± standard deviation values. doses of the Pfizer-BioNTech COVID-19 mRNA vaccine. Solid or dashed lines indicate model simulations; markers with error bars represent mean ± standard deviation values. Subsequently, as shown in **Figure 2b**, we calibrated the model with the clinical data obtained from healthy individuals vaccinated with mRNA vaccines (specifically Pfizer-BioNTech) [46]. For this purpose, a double dose of the vaccine was simulated in accordance with the schedule used for the individuals in the study [46]. A Gaussian function described the kinetics of antigen load following injections on days 0 and 28 (**Eq. 1**). The solutions for the various immune response variables were computed over a period of eight months and fitted to the available clinical data for effector T-cells (CD4+ and CD8+) and neutralizing antibodies. Our results showed a high degree of correlation between the model fits and clinical measurements (**Figure S3b**; *R* > 0.96). To ensure that the model can reproduce immune responses elicited by the vaccines over long time periods, some of the parameters were refitted (**Table 1**). Since during the previous calibration, the characteristic time of simulation is a few weeks unlike the current simulation where the simulated time is a few months, we recalibrated some parameters to ensure long term accuracy of the simulation. Also, to adjust the model for it to be able to capture any fundamental differences between response to infection and vaccines, we performed the recalibration. In addition, some parameters required recalibration because of the variation in units of measurement between experiments. An important observation is the gradually waning levels of neutralizing antibodies and effector lymphocytes, which suggests that protection conferred by mRNA vaccines is temporal, warranting the use of boosters. To accurately represent the vaccine-induced immune response in immunocompromised individuals, we also calibrated the model with clinical data obtained from vaccinated cancer patients undergoing chemotherapy or immunotherapy (**Figure 2c,d**) [47]. In both cases, we assumed that due to the underlying pathophysiology and associated treatment, the levels of some immune system parameters were only a fraction (0 < *f* < 1) of their values in healthy individuals (*f* = 1). Therefore, keeping all other model parameters from the previous two fits as constants, we fitted the model to two datasets [47] to estimate the parameter *f*, which resulted in a value of *f* = 0.517for chemotherapy-treated patients and *f* = 0.588 for immunotherapy-treated cancer patients. The model fits were also in good agreement with clinical data (**Figure S3c,d**; *R* > 0.96). ### 3.2 Model validation To test the ability of our model to accurately reproduce the immune response to mRNA vaccines, we simulated two and three doses of the Pfizer-BioNTech and Moderna COVID-19 vaccines in healthy individuals (data not used for calibration). As shown in **Figures 3** and **S4** (*R* = 0.92), the computed neutralizing antibody (IgG) kinetics closely resembles the literature-derived clinical data following two doses of the Pfizer-BioNTech COVID-19 vaccine [62], two doses of the Moderna COVID-19 vaccine [46], and three doses of the Pfizer-BioNTech COVID-19 vaccine [78]. The dosing schedules were obtained from the respective clinical studies, and the parameter values were based on the values calibrated for healthy individuals in the previous section (**Table 1**). The ability of the model to accurately predict the response to the third dose, despite not using the third dose data during model calibration, highlights the biological and physiological robustness of our mechanistic model. Having established the validity of our model to reliably reproduce neutralizing antibody kinetics with various mRNA vaccines and dosing schedules, we proceeded to perform numerical experiments to explore the heterogeneity in immune responses and optimize dosing schedules to minimize breakthrough infections. ![Figure 3.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/01/31/2022.09.14.22279959/F3.medium.gif) [Figure 3.](http://medrxiv.org/content/early/2023/01/31/2022.09.14.22279959/F3) Figure 3. Model validation. Validation of the mathematical model with antibody kinetics data derived from the literature for healthy individuals vaccinated with two doses of Pfizer-BioNTech COVID-19 mRNA vaccine (red squares), two doses of Moderna COVID-19 mRNA vaccine (black circles), and three doses of Pfizer-BioNTech COVID-19 mRNA vaccine (blue triangles). Solid line indicates model predictions, grey bands represent 90% prediction intervals, and markers with error bars represent mean ± standard deviation values of clinical data. Yellow diamonds on the x-axis denote timing of injection (i.e., first dose given on day 0, second dose given on day 21 after first dose, and third dose given 9 months after second dose). ### 3.3 Sensitivity analysis To study the relative significance of model parameters in governing humoral response to mRNA-based COVID-19 vaccines, quantified as area under the curve of neutralizing antibody titer kinetics following a unit dose of the vaccine, global sensitivity analysis (GSA) and local sensitivity analysis (LSA) were performed with 22 model parameters that characterize the key immunological interactions, processes, and immune states considered to be important for vaccine-induced immune response generation. As shown in **Figure 4**, the immunosuppression factor *f* stands out as the most relevant parameter in determining antibody response to vaccines, which highlights the importance of immune health status (i.e., immune cell sufficiency) in governing vaccine-induced protection. Further, as shown in the inset of **Figure 4**, LSA reveals a positive monotonic correlation between change in *f* and sensitivity index (SI) for antibody titer, which signifies better antibody response in immunologically competent subjects, and thus warrants the need for additional doses or optimized dosing frequency in immunocompromised patients [8]. Following this, biological parameters that characterize antibody production *P*Ab, death of antibody secreting plasma cells *δ**P*, and antigen-induced activation of naïve APCs, *T*APC were observed to also influence antibody response strongly (**Figure 4**), thereby indicating the relevance of antigen presentation, plasma cell population, and antibody secretion from plasma cells in humoral immunity development. While the relationship between change in antibody titer and perturbations of parameters *P*Ab or *T*APC is monotonically increasing within the studied parameter range, that between antibody titer and *δ**P* is monotonically decreasing (**Figure 4, inset**). Further, additional parameters belonging to CD4+ T-cell activation *T*CD4, growth of B cells *γ**B* death of activated APCs *δ*APC, differentiation of B cells into antibody secreting plasma cells *T*BC, potency of the injected antigen to activate naïve APCs, v, and activation of naïve B cells, are also observed to have a moderately significant effect on antibody titer (**Figure 4**). This further confirms the relevance of previously identified processes, in addition to the levels of activated CD4+ T-cells and B cells in governing antibody response. While the parameters identified in this analysis do not contribute directly to immune health status in our model, except *f*, their significance in antibody response warrants their inclusion in virtual patient cohort generation to study the effect of the underlying population-scale biological variability on vaccine-induced protection. ![Figure 4.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/01/31/2022.09.14.22279959/F4.medium.gif) [Figure 4.](http://medrxiv.org/content/early/2023/01/31/2022.09.14.22279959/F4) Figure 4. Parameter sensitivity analyses. Violin plot showing parameters ranked in descending order (from left to right) for their sensitivity, as indicated by the sensitivity index (SI) obtained through global sensitivity analysis (GSA). Correlation of parameter perturbation and its effect on antibody titer (quantified as SI) obtained through local sensitivity analysis (LSA) for the top ten ranking parameters of GSA is shown in the inset. Note that parameters were perturbed in the range of ± 50% around their baseline values for both GSA and LSA. ### 3.4 Heterogeneity in immune response to vaccines at the individual and population scale To study the influence of (i) vaccine dosing schedules and (ii) the immune status of an individual on neutralizing antibody levels and vaccine efficacy, we simulated immune responses under different dosing schedules in representative healthy and immunocompromised subjects. Based on the dosing schedules used across various countries, we considered three vaccination regimens: *rapid, intermediate*, and *delayed*. In all cases the first dose was given on day 0; (i) Rapid: second dose is given two weeks after the first dose, and the first booster (third dose) is given five months after the second dose; (ii) Intermediate: second dose is given four weeks after the first dose, and the first booster is given seven months after the second dose; (iii) Delayed: second dose is given eight weeks after the first dose, and the first booster (third dose) is given nine months after the second dose. Here, the immune health status was defined by the non-dimensional, empirical parameter *f*, such that healthy individuals have *f* = 1, mildly immunocompromised subjects have *f* = 0.75, and highly immunocompromised individuals have *f* = 0.55. As previously discussed, *f* = 0.55 corresponds to cancer patients undergoing chemotherapy or immunotherapy, whereas *f* = 0.75 simulates individuals with underlying conditions that may also affect the immune system, but usually to a lesser degree, (e.g., autoimmune diseases). As per evidence in the literature, plasma antibody titer is a correlate of protection against infection [65, 66]. Therefore, we used the computed neutralizing antibody levels as predictors of vaccine efficacy (i.e., protection against SARS-CoV-2; see section 2.3 in **Methods**). We used the model to predict the humoral response to mRNA vaccines following three dosing schedules in representative healthy or immunocompromised individuals for a 600-day period. As shown in **Figure 5a,d,g** (upper subplot in each panel), irrespective of the dosing schedule, the antibody levels remain above the protection threshold for both OM (770 U/mL) and WT strain (154 U/mL) for a much longer duration in healthy individuals (as indicated by the shaded grey area and quantified as the *T*safe value = 383–443 days) than in mildly immunocompromised subjects (**Figure 5b,e,h**; *T*safe = 162–228 days). In contrast, in cancer patients undergoing antineoplastic treatment (i.e., highly immunocompromised subjects), *T*safe was zero days across all dosing schedules (**Figure 5c,f,i**). This suggests that highly immunocompromised subjects are vulnerable to infection with OM throughout the 600-day simulation period, however protection against WT is intermittently present depending upon the dosing schedule. Of note, within both the healthy and mildly immunocompromised subjects, the intermediate dosing schedule leads to higher *T*safe values (443 days if healthy, 228 days if mildly immunocompromised; **Figure 5d,e**) than the rapid dosing schedule (383 days if healthy, 162 days if mildly immunocompromised; **Figure 5a,b**) or the delayed dosing schedule (396 days if healthy, 182 days if mildly immunocompromised; **Figure 5g,h**). Nonetheless, the protection window in these cases is not continuous for the chosen dosing schedules, and an intermediate ‘gap’ is observed between the second dose and first booster (third dose) that highlights the period when antibody levels temporarily fall below the protective threshold for OM and/or WT. The duration of this gap varies according to immune health status and dosing schedule. ![Figure 5.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/01/31/2022.09.14.22279959/F5.medium.gif) [Figure 5.](http://medrxiv.org/content/early/2023/01/31/2022.09.14.22279959/F5) Figure 5. Effect of vaccine dosing schedule and immune health status on antibody levels and vaccine efficacy. Simulations in representative (**a,d,g**) healthy and (**b,c,e,f,h,i**) immunocompromised individuals show antibody levels and vaccine efficacy against wildtype strain (WT) and Omicron variant (OM) of SARS-CoV-2 following (**a,b,c**) rapid, (**d,e,f**) intermediate, and (**g,h,i**) delayed vaccine dosing schedules. Yellow diamonds on the x-axes, in the upper sub-panel, indicate injection timepoints. In each upper sub-panel, the black solid line represents antibody levels, with the dashed blue and red lines indicating protective threshold against WT and OM, respectively. The lower sub-panel shows vaccine efficacy (colored solid line for WT and shaded area for OM), with the dashed black line indicating the 82.3% threshold of protection. Note: The value *T*safe indicated in every upper sub-panel represents the number of days when antibody levels are above the protective threshold for both WT and OM. The corresponding vaccine efficacy kinetics are shown in the lower subplots in **Figure 5**. The shaded area represents the vaccine efficacy against OM, and the solid-colored line indicates the vaccine efficacy against WT. The continuous color mapping assigns blue to efficacies above the protection threshold (>82.3%) and red to efficacies equal to or below the protection threshold (≤ 82.3%). As visible from the bluish region of the shaded area, for any given dosing schedule, protection threshold healthy individuals have greater vaccine efficacy against OM than immunocompromised individuals. In highly immunocompromised subjects the shaded area always remained below the protective threshold (82.3%), indicating a high risk of becoming infected with OM (**Figure 5c,f,i**). As expected, due to limited *antibody escape* [79], the vaccine efficacy against WT is greater than that against OM in all individuals under all dosing schedules (as indicated by the colored solid line). Further, in healthy individuals, the three dosing schedules produced antibody titers above the WT protection threshold for the majority of the simulation period (**Figure 5a,d,g**). In mildly immunocompromised individuals, protection against WT does not persist continuously (**Figure 5b,e,h**). For example, in the delayed dosing schedule shown in **Figure 5h**, the period between day 236 and day 330 (∼three months) indicates a vaccine efficacy of less than 82.3%. In highly immunocompromised cases, the three dosing schedules provide limited protection against WT, with prolonged periods of lapse in immunity. Though we only considered representative individuals, these observations collectively highlight the importance of optimizing the dosing schedule based on the immune health status of a sub-population to achieve continuous, long-term protection against both WT and other VOCs (e.g., OM). To evaluate the effects of dosing schedules and immune health status on the variability in immune response to mRNA vaccines at the population level, we simulated the vaccination of a virtual population of 10,000 individuals with three doses (**Cohort A**; see section 2.5 in **Methods** for details of dosing schedule) and assessed the corresponding vulnerability to breakthrough infections. Note that the dosing schedule for each simulated individual was obtained randomly from a continuous time interval (red and blue brackets on x-axis of **Figure 6**) to replicate the real-world heterogeneity in dosing time intervals. As shown in **Figure 6a**, the average antibody kinetics across the 10,000 individuals remained above the protective threshold for OM and WT. However, for a significant fraction of the population, antibody levels remained below the OM threshold for a prolonged period (∼five months). This is evident from the shaded area representing 90% prediction interval. Further, translating the antibody levels to vaccine efficacy using **Eq. 19**, we observed that for a significant fraction of the 10,000 individuals, vaccine efficacy against OM fell below the 82.3% protection threshold (see **Figure 6b**, orange shaded area). Subsequently, we quantified the fraction of the virtual population that presented a vaccine efficacy below the protective threshold for OM and WT (**Figure 6c**, see section 2.6 in **Methods)**. This population fraction can alternatively be interpreted as the fraction of vaccinated individuals in a population that is vulnerable to breakthrough infections, i.e., becoming infected despite being vaccinated. As observed in **Figure 6c**, this fraction increases to about 0.5 (or ∼50% of the population) for OM in vaccinated individuals (two doses), and then declines rapidly following administration of the first booster (third dose). However, due to waning antibody levels, which translate into declining efficacy, the vulnerable fraction begins to increase again and becomes 1 (i.e., 100% of population) in about six months after the booster window. In contrast, for WT, the vulnerable fraction of the population peaks at about 0.1 (or ∼10% of the population) in vaccinated individuals (two doses) and then decreases again after administration of the first booster (third dose), suggesting effective protection against WT in vaccinated individuals for up to ∼1.5 years, irrespective of the dosing schedule or immune health status. Of note, in the population-scale simulation, immune health status is non-uniformly distributed across the population, as defined by the left half-Gaussian distribution (**Figure S2c**); this indicates that a major proportion of the population is healthy. It is worth mentioning that the sharp rise in **Figure 6c** of the population fraction several months after the first booster warrants the administration of a second booster to curb the vulnerability to VOCs and WT. Given that the proposed dosing schedules do not warrant continuous protection against VOCs and/or WT, it is imperative to optimize the schedules to achieve long-term protection in the population without lapses. ![Figure 6.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/01/31/2022.09.14.22279959/F6.medium.gif) [Figure 6.](http://medrxiv.org/content/early/2023/01/31/2022.09.14.22279959/F6) Figure 6. Effect of heterogeneity in vaccine dosing schedules and immune health status on breakthrough infections at the population scale. **a**) Average antibody levels in plasma, **b**) corresponding vaccine or antibody efficacy, and **c**) population fraction vulnerable to breakthrough infections due to wildtype strain (WT, solid blue line) and Omicron variant (OM, dotted orange line) of SARS-CoV-2 over time. Solid and dotted lines in **a,b**) represent average behavior of 10,000 simulated individuals and shaded bands indicate 90% prediction interval. Note that the first dose was administered on day 0 to each simulated individual, the second dose was administered between day 14 and day 56, and the third dose (i.e., first booster) was administered between day 150 and day 270. Red and blue brackets on x-axis denote timing windows with respect to day 0 for second dose and third dose, respectively, used to design unique vaccine schedules in model simulations. Immune health status (*f*) of the simulated population varied between 0.5 to 1. ### 3.5 Vaccine dosing schedule optimization Following the previous numerical experiments, we intended to identify optimal vaccine dosing schedules to achieve continuous protection against OM (as a representative example) for prolonged periods. We generated three virtual cohorts of 10,000 individuals (**Cohort B**) each to represent healthy, mildly immunocompromised, and highly immunocompromised individuals, and implemented several dosing schedules to identify optimal times for the second dose, the third dose (first booster), and the fourth dose (second booster) in each sub-cohort (see section 2.6 in **Methods**). As shown in **Figure 7**, the AUC of vulnerability kinetics curves follows a non-linear relationship with respect to dosing schedules, and a minimum is visible for each dose and population sub-type (highlighted by a red circle). As shown in **Figure 7a,d,g**, as the immune status changes from healthy to highly immunocompromised, the position of the minima on the x-axis shows a right shift, such that the optimal time for the second dose in healthy, mildly immunocompromised, and highly immunocompromised individuals is 18 days, 25 days, and 30 days after the first dose, respectively. In contrast, as shown in **Figure 7b,e,h**, the minima for the first booster shows a left shift on the x-axis from healthy to highly immunocompromised individuals, such that the optimal time for first booster is 164 days (∼5.5 months), 115 days (∼4 months), and 36 days (1.2 months) after the second dose for healthy, mildly immunocompromised, and highly immunocompromised individuals, respectively. Similarly, as shown in **Figure 7c,f,i**, the minima for the 2nd booster shows a left shift from healthy to highly immunocompromised individuals, such that the optimal schedule for the second booster is 223 days (∼7.5 months), 195 days (∼6.5 months), and 126 days (∼4 months) after the first booster for healthy, mildly immunocompromised, and highly immunocompromised individuals, respectively. ![Figure 7.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/01/31/2022.09.14.22279959/F7.medium.gif) [Figure 7.](http://medrxiv.org/content/early/2023/01/31/2022.09.14.22279959/F7) Figure 7. COVID-19 vaccine dosing schedule optimization. Area under the curve (AUC) of breakthrough infection vulnerability kinetics curve obtained from simulation of 10,000 individuals from different population sub-types under unique dosing schedules and immune health status. Estimated AUC versus dosing schedules for **a,d,g**) dose 2, **b,e,h**) booster 1 (i.e., dose 3), and **c,f,i**) booster 2 (i.e., dose 4) for **a,b,c**) healthy, **d,e,f**) mildly immunocompromised, and **g,h,i**) highly immunocompromised individuals, obtained through model simulations. Each black dot represents one AUC value. Red dot in each plot represents the corresponding minima for each dose and population sub-type. It is intuitive to expect inter-vaccination periods to be longer for healthy individuals than for immunocompromised patients. This is evidenced by data presented in **Figure 5**, where the antibody titer stays above the OM protection threshold for a longer period in healthy individuals than in their immunocompromised counterparts, thereby allowing the possibility to delay subsequent doses. Although this is true for the first and the second boosters (**Figure 7b** vs. **Figures 7e,h**; **Figure 7c** vs. **Figures 7f,i**), the trend is reversed for the second dose (**Figure 7a** vs. **Figures 7d,g**), where healthy individuals seem to require the second dose sooner than immunocompromised individuals to ensure continuity of protection against OM. This observation can be explained in light of a key mechanistic assumption of our model. Recall that the immune health status parameter *f* scales the homeostasis level of naïve immune cells (![Graphic][29]). In immunocompromised individuals, *f* ranges from 0.5 to 0.9; therefore, the homeostasis level of naïve immune cells is less than that in healthy individuals (**Figures 2b,c,d**). As a result, when the second dose is given too soon after the first dose in immunocompromised individuals, due to reduced levels of CD4+ T-cells and therefore slower activation of B cells, the production of neutralizing antibodies from plasma cells may be thwarted, thereby rendering an individual vulnerable to infection. Therefore, permitting the CD4+ T-cell and B cell population to regenerate after the first dose will allow antibody titers to rise to levels associated with adequate protection. Of note, since healthy individuals produce or activate immune cells more quickly (given *f* = 1), they are ready to receive a second dose sooner than immunocompromised individuals. However, in the case of healthy individuals, as shown in **Figure 7a**, the AUC0-150 d values are smaller than those of immunocompromised patients (**Figure 7d,g**) for up to ∼six weeks of delay after the first dose. A six-week delay after the first dose predisposes ∼30% (obtained from the ratio of AUC0-150 d value at six weeks (i.e., ∼45) to maximum possible value of AUC0-150 d (i.e., 150)) of the healthy population to a breakthrough infection over 150 days under no public health restrictions. This indicates that although an optimal waiting period for healthy individuals is two and a half weeks after the first dose (which predisposes only ∼15% healthy population over 150 days), if required due to logistic constraints, waiting longer (up to six weeks) will still allow the healthy individuals to be more protected than immunocompromised individuals. ### 3.6 Testing model-predicted optimal dosing schedules Finally, to demonstrate the impact of the previously identified optimal dosing schedules (for the second dose and the two boosters) in reducing vulnerability to breakthrough infections, we simulated a vaccination regimen with four doses in 10,000 virtual individuals per group, belonging to the three cohorts of interest (healthy, mildly immunocompromised, and highly immunocompromised; see **Cohort B** in section 2.5 in **Methods**), and measured the vaccine efficacy and corresponding level of vulnerability to infection over a period of two years. Note that to simulate a more realistic test scenario and add variability to the optimal dosing schedules identified previously (**Figure 7**), we sampled the dosing schedules from within a ± 10 % uniform distribution around the optimal values. As shown in **Figures 8a,c,e**, the average vaccine efficacy for WT and OM is above the protection threshold in all sub-populations for an extended period of time, the duration of which is dependent on the viral strain and population sub-type. Therefore, the corresponding vulnerability to breakthrough infections for OM and WT remains at almost zero for most of the two-year period in healthy individuals and shows only two intermittent windows of ∼two months each where the vulnerability is as high as ∼0.065 (**Figure 8b**). In mildly immunocompromised individuals, the optimized protocol exhibits similar results, although the vulnerability to infection after fourth dose begins to rise sooner in comparison to the healthy population (**Figure 8d**). Furthermore, as shown in **Figure 8f**, in the highly immunocompromised cohort, the same trend continues; although complete protection against OM and WT is observed for a shorter duration, the results are nonetheless notably more promising compared to the observed findings in **Figures 5c,f,i**, where vaccine efficacy remained below the OM protection threshold throughout the 600-day window under conventional dosing schedules. ![Figure 8.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/01/31/2022.09.14.22279959/F8.medium.gif) [Figure 8.](http://medrxiv.org/content/early/2023/01/31/2022.09.14.22279959/F8) Figure 8. Testing model-predicted optimal dosing schedules. **a,c,e)** Vaccine efficacy and **b,d,f**) vulnerability to breakthrough infections due to wildtype strain (WT, solid blue line) and Omicron variant (OM, dotted orange line) of SARS-CoV-2 in **a,b**) healthy, **c,d**) mildly immunocompromised, and **e,f**) highly immunocompromised individuals. For each population sub-type, testing was done on 10,000 simulated individuals with unique *f* and dosing schedule values. Colored bands represent 90% prediction intervals. Finally, the optimal dosing schedules identified above are summarized in **Figure 9** (green bands), with a comparison made to the CDC-recommended dosing schedules being currently implemented for the Pfizer-BioNTech vaccine (blue bands). The ongoing CDC guidelines for COVID-19 vaccination for healthy individuals (not moderately or severely immunocompromised and <50 years of age) include 3 doses with intervals of 3-8 weeks between the first and second dose (represented as 21 days) and 5 months between the second and third dose (represented as 140 days). The model-predicted schedule closely recapitulates the CDC guidelines with the inclusion of a fourth dose to prolong immunity for 385 days (>1 year). Although the model distinguishes between two immunocompromised cancer populations (mildly and highly), the CDC guidelines suggest a schedule of 4 doses for patients who are moderately or severely immunocompromised (with intervals of 21, 21, and 84 days, respectively). According to the model-predicted optimal dosing schedule, longer gaps between doses (or boosters) would not compromise the immunity of healthy and immunocompromised patients that could represent a solution to logistic constraints. ![Figure 9.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/01/31/2022.09.14.22279959/F9.medium.gif) [Figure 9.](http://medrxiv.org/content/early/2023/01/31/2022.09.14.22279959/F9) Figure 9. Model-predicted optimal dosing and CDC-recommended dosing schedules for the Pfizer-BioNTech vaccine in healthy and immunocompromised sub-populations. The ongoing CDC guidelines for dosing schedules are represented by the blue bands and those predicted by the model are shown in green. ## 4. Conclusions In summary, we developed a mechanistic mathematical model of adaptive immune response to COVID-19 vaccines and viral infection in healthy and immunocompromised individuals. The model was formulated as a system of ODEs to account for key biological interactions leading to the development of antigen-induced humoral and cellular immunity. Following the calibration and validation of the model with published clinical data, numerical experiments were performed to study the effects of immune health status and vaccine dosing schedules on plasma antibody titers (a correlate of protection against infection), leading to the estimation of population vulnerability to breakthrough infections. Through simulations of virtual individuals, the model was then applied to identify optimal dosing schedules of the vaccines to minimize breakthrough infections in the population. Through our analysis, we highlighted critical waiting windows for immunocompromised individuals (25 days and 30 days after first dose for mildly and highly immunocompromised individuals, respectively) to ensure sufficient time for the development of immune recall responses and minimize vulnerability to breakthrough infections in their sub-populations. In the case of healthy individuals, while the optimal waiting period after first dose was found to be two and a half weeks, we proposed that it can be extended (without much compromise to protection) up to six weeks. Thereby, we make the case for longer waiting period between doses without compromising the immunity at the population scale. The presented model is based on generalized adaptive immune response to antigens and can thus be adapted to investigate different infections or vaccines, given appropriate data for model calibration. Through our proof-of-concept study, we have thus presented a novel approach to optimizing vaccine dosing schedules in case of future outbreaks. While the current model is thoroughly calibrated and validated, the mechanistic underpinnings of immunosuppression and innate immune response need to be considered in greater detail in future studies. Also, model adaptations relevant to the other types of COVID-19 vaccines may need to be considered as well. Lastly, our results also suggest the need for follow-up boosters (more frequently for immunocompromised subjects due to rapidly waning immunity) to ensure continued immunity against breakthrough infections and reinfections, especially given the emergence of novel VOCs. ## Supporting information Supplementary Information [[supplements/279959_file03.docx]](pending:yes) ## Data Availability All data produced in the present study are available upon reasonable request to the authors. [https://www.nature.com/articles/s41586-020-2588-y](https://www.nature.com/articles/s41586-020-2588-y) [https://www.nejm.org/doi/10.1056/NEJMc2115596](https://www.nejm.org/doi/10.1056/NEJMc2115596) [https://www.esmoopen.com/article/S2059-7029(21)00236-2/fulltext](https://www.esmoopen.com/article/S2059-7029(21)00236-2/fulltext) [https://www.mdpi.com/2076-393X/9/10/1092](https://www.mdpi.com/2076-393X/9/10/1092) [https://www.mdpi.com/2076-393X/10/6/876](https://www.mdpi.com/2076-393X/10/6/876) ## Author contributions PD conceived and supervised the study. PD and CS curated the data, developed the mathematical model, and performed modeling analysis and simulations. PD, CS, ZW, JRR, SC, and VC designed numerical experiments. PD, CS, DIS, CM, JW, HDS, RP, and WA helped with interpretation of modeling results. All authors contributed to manuscript writing and editing. ## Competing Interests DIS, RP, and WA are listed as inventors on a patent application related to immunization strategies (International Patent Application PCT/ US2020/053758, entitled Targeted Pulmonary Delivery Compositions and Methods Using Same). DIS, CM, RP, and WA are inventors on International Patent Application PCT/US2021/040392, entitled Enhancing Immune Responses Through Targeted Antigen Expression, which describes immunization technology adapted for COVID-19. PhageNova Bio has licensed these intellectual properties and DIS, CM, RP, and WA may be entitled to standard royalties. RP and WA are founders and equity stockholders of PhageNova Bio. RP is Chief Scientific Officer and a paid consultant of PhageNova Bio. RP and WA are founders and equity shareholders of MBrace Therapeutics; RP is a Board Member and paid consultant, and WA is a Scientific Advisor at MBrace Therapeutics. These arrangements are managed in accordance with the established institutional conflict-of-interest policies of Rutgers, The State University of New Jersey. The remaining authors declare no competing interests. ## Acknowledgements This study was conducted under the umbrella of the International Academic Affiliation Agreement between the Houston Methodist Academic Institute (Houston, TX, USA) and the University of Naples Federico II (Naples, Italy). The work was supported in part by the Cockrell Foundation (PD, VC), the National Institutes of Health (NIH) Grants 1R01CA253865 (ZW, VC), 1R01CA222007 (ZW, VC), and 1R01CA226537 (ZW, VC, RP and WA), Rutgers Cancer Institute of New Jersey (NCI Cancer Center Support Grant number P30CA072720), and by awards from the Levy-Longenbaugh Donor-Advised Fund (to RP and WA). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. We thank Prof. Alan S. Perelson for his valuable criticism of our work. PD and CS also acknowledge Luca Messina and Maria J. Peláez for helpful scientific discussions, Rachael E. Whitehead for her contributions to the model schematic, and Life Science Editors for editing support. ## Footnotes * This version of the manuscript has been revised to update the following: include additional analysis (Figure 4), updated the previous figure 3, and improved the description of the model and results. * Received September 14, 2022. * Revision received January 31, 2023. * Accepted January 31, 2023. * © 2023, Posted by Cold Spring Harbor Laboratory This pre-print is available under a Creative Commons License (Attribution-NonCommercial 4.0 International), CC BY-NC 4.0, as described at [http://creativecommons.org/licenses/by-nc/4.0/](http://creativecommons.org/licenses/by-nc/4.0/) ## References 1. 1.WHO. WHO COVID-19 Dashboard Geneva: World Health Organization; 2020 [updated 202004/12/2022]. Available from: [https://covid19.who.int/](https://covid19.who.int/) 2. 2.Mukaigawara M, Hassan I, Fernandes G, King L, Patel J, Sridhar D. An equitable roadmap for ending the COVID-19 pandemic. Nature Medicine. 2022. doi: 10.1038/s41591-022-01787-2. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41591-022-01787-2&link_type=DOI) 3. 3.Dogra P, Koay EJ, Wang Z, Vahidy FS, Ferrari M, Pasqualini R, et al. Is the worst of the COVID-19 global pandemic yet to come? Application of financial mathematics as candidate predictive tools. Translational Psychiatry. 2021;11(1):299. doi: 10.1038/s41398-021-01429-0. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41398-021-01429-0&link_type=DOI) 4. 4.Team VGCVT. COVID19 Vaccine Tracker VIPER Group COVID19 Vaccine Tracker Team; 2022 [12/2/2022]. Available from: [https://covid19.trackvaccines.org/](https://covid19.trackvaccines.org/). 5. 5.Our World in Data: Coronavirus (COVID-19) Vaccinations: Our World in Data; [4/13/2022]. Available from: [https://ourworldindata.org/covid-vaccinations?country=OWID_WRL](https://ourworldindata.org/covid-vaccinations?country=OWID_WRL). 6. 6.Pilkington V, Keestra SM, Hill A. Global COVID-19 Vaccine Inequity: Failures in the First Year of Distribution and Potential Solutions for the Future. Frontiers in Public Health. 2022;10. doi: 10.3389/fpubh.2022.821117. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.3389/fpubh.2022.821117&link_type=DOI) 7. 7. Lee Aryb, Wong SY, Chai LYA, Lee SC, Lee MX, Muthiah MD, et al. Efficacy of covid-19 vaccines in immunocompromised patients: systematic review and meta-analysis. BMJ. 2022;376:e068632. doi: 10.1136/bmj-2021-068632. [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6MzoiYm1qIjtzOjU6InJlc2lkIjtzOjE5OiIzNzYvbWFyMDJfNC9lMDY4NjMyIjtzOjQ6ImF0b20iO3M6NTA6Ii9tZWRyeGl2L2Vhcmx5LzIwMjMvMDEvMzEvMjAyMi4wOS4xNC4yMjI3OTk1OS5hdG9tIjt9czo4OiJmcmFnbWVudCI7czowOiIiO30=) 8. 8.Parker EPK, Desai S, Marti M, Nohynek H, Kaslow DC, Kochhar S, et al. Response to additional COVID-19 vaccine doses in people who are immunocompromised: a rapid review. The Lancet Global Health. 2022;10(3):e326–e8. doi: 10.1016/S2214-109X(21)00593-3. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S2214-109X(21)00593-3&link_type=DOI) 9. 9.Aschwanden C. Five reasons why COVID herd immunity is probably impossible. Nature. 2021:520–2. 10. 10.Lin D-Y, Gu Y, Wheeler B, Young H, Holloway S, Sunny S-K, et al. Effectiveness of Covid-19 Vaccines over a 9-Month Period in North Carolina. New England Journal of Medicine. 2022;386(10):933–41. doi: 10.1056/NEJMoa2117128. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1056/NEJMoa2117128&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=35020982&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 11. 11.Levin EG, Lustig Y, Cohen C, Fluss R, Indenbaum V, Amit S, et al. Waning Immune Humoral Response to BNT162b2 Covid-19 Vaccine over 6 Months. New England Journal of Medicine. 2021;385(24):e84. doi: 10.1056/NEJMoa2114583. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1056/NEJMoa2114583&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=34614326&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 12. 12.Peng Q, Zhou R, Wang Y, Zhao M, Liu N, Li S, et al. Waning immune responses against SARS-CoV-2 variants of concern among vaccinees in Hong Kong. eBioMedicine. 2022;77. doi: 10.1016/j.ebiom.2022.103904. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.ebiom.2022.103904&link_type=DOI) 13. 13.Pérez-Alós L, Armenteros JJA, Madsen JR, Hansen CB, Jarlhelt I, Hamm SR, et al. Modeling of waning immunity after SARS-CoV-2 vaccination and influencing factors. Nature Communications. 2022;13(1):1614. doi: 10.1038/s41467-022-29225-4. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41467-022-29225-4&link_type=DOI) 14. 14.Khoury DS, Cromer D, Reynaldi A, Schlub TE, Wheatley AK, Juno JA, et al. Neutralizing antibody levels are highly predictive of immune protection from symptomatic SARS-CoV-2 infection. Nat Med. 2021;27(7):1205-11. Epub 2021/05/19. doi: 10.1038/s41591-021-01377-8. PubMed PMID: 34002089. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41591-021-01377-8&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=34002089&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 15. 15.Tartof SY, Slezak JM, Fischer H, Hong V, Ackerson BK, Ranasinghe ON, et al. Effectiveness of mRNA BNT162b2 COVID-19 vaccine up to 6 months in a large integrated health system in the USA: a retrospective cohort study. The Lancet. 2021;398(10309):1407–16. doi: [https://doi.org/10.1016/S0140-6736(21)02183-8](https://doi.org/10.1016/S0140-6736(21)02183-8). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S0140-3946736(21)02183-8&link_type=DOI) 16. 16.Irizarry RA, Robles-Fontán MM, Nieves EG, Cardona-Gerena I. Time-Varying Effectiveness of Three COVID-19 Vaccines in Puerto Rico. Available at SSRN 3957118. 2021. 17. 17.Singh J, Pandit P, McArthur AG, Banerjee A, Mossman K. Evolutionary trajectory of SARS-CoV-2 and emerging variants. Virology Journal. 2021;18(1):166. doi: 10.1186/s12985-021-01633-w. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1186/s12985-021-01633-w&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=34389034&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 18. 18.Markosian C, Staquicini DI, Dogra P, Dodero-Rojas E, Lubin JH, Tang FHF, et al. Genetic and Structural Analysis of SARS-CoV-2 Spike Protein for Universal Epitope Selection. Molecular Biology and Evolution. 2022;39(5). doi: 10.1093/molbev/msac091. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/molbev/msac091&link_type=DOI) 19. 19.Hacisuleyman E, Hale C, Saito Y, Blachere NE, Bergh M, Conlon EG, et al. Vaccine Breakthrough Infections with SARS-CoV-2 Variants. New England Journal of Medicine. 2021;384(23):2212–8. doi: 10.1056/NEJMoa2105000. PubMed PMID: 33882219. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1056/NEJMoa2105000&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=33882219&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 20. 20.Bergwerk M, Gonen T, Lustig Y, Amit S, Lipsitch M, Cohen C, et al. Covid-19 Breakthrough Infections in Vaccinated Health Care Workers. New England Journal of Medicine. 2021;385(16):1474–84. doi: 10.1056/NEJMoa2109072. PubMed PMID: 34320281. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1056/NEJMoa2109072&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=34320281&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 21. 21.Bar-On YM, Goldberg Y, Mandel M, Bodenheimer O, Freedman L, Kalkstein N, et al. Protection of BNT162b2 Vaccine Booster against Covid-19 in Israel. New England Journal of Medicine. 2021;385(15):1393–400. doi: 10.1056/NEJMoa2114255. PubMed PMID: 34525275. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1056/NEJMoa2114255&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=34525275&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 22. 22.Bar-On YM, Goldberg Y, Mandel M, Bodenheimer O, Amir O, Freedman L, et al. Protection by a Fourth Dose of BNT162b2 against Omicron in Israel. New England Journal of Medicine. 2022. doi: 10.1056/NEJMoa2201570. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1056/NEJMoa2201570&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=35381126&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 23. 23.Regev-Yochay G, Gonen T, Gilboa M, Mandelboim M, Indenbaum V, Amit S, et al. Efficacy of a Fourth Dose of Covid-19 mRNA Vaccine against Omicron. New England Journal of Medicine. 2022;386(14):1377–80. doi: 10.1056/NEJMc2202542. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1056/NEJMc2202542&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=35297591&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 24. 24.Lusvarghi S, Pollett SD, Neerukonda SN, Wang W, Wang R, Vassell R, et al. SARS-CoV-2 BA.1 variant is neutralized by vaccine booster-elicited serum, but evades most convalescent serum and therapeutic antibodies. Science Translational Medicine. ():eabn8543. doi: doi:10.1126/scitranslmed.abn8543. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1126/scitranslmed.abn8543&link_type=DOI) 25. 25.Ritchie H, Mathieu E, Rodés-Guirao L, Appel C, Giattino C, Ortiz-Ospina E, et al. Coronavirus Pandemic (COVID-19): Our World in Data; 2020 [04/15/2022]. Available from: [https://ourworldindata.org/coronavirus](https://ourworldindata.org/coronavirus). 26. 26.Williamson EJ, Walker AJ, Bhaskaran K, Bacon S, Bates C, Morton CE, et al. Factors associated with COVID-19-related death using OpenSAFELY. Nature. 2020;584(7821):430–6. doi: 10.1038/s41586-020-2521-4. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41586-020-2521-4&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=32640463&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 27. 27.Control CfD. Interim COVID-19 Immunization Schedule for Ages 5 Years and Older: Center for Disease Control; 2022 [4/15/2022]. Available from: [https://www.cdc.gov/vaccines/covid-19/downloads/COVID-19-immunization-schedule-ages-5yrs-older.pdf](https://www.cdc.gov/vaccines/covid-19/downloads/COVID-19-immunization-schedule-ages-5yrs-older.pdf). 28. 28.Amirthalingam G, Bernal JL, Andrews NJ, Whitaker H, Gower C, Stowe J, et al. Serological responses and vaccine effectiveness for extended COVID-19 vaccine schedules in England. Nature Communications. 2021;12(1):7217. doi: 10.1038/s41467-021-27410-5. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41467-021-27410-5&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=34893611&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 29. 29.Wong CKH, Xiong X, Lau KTK, Chui CSL, Lai FTT, Li X, et al. Impact of a delayed second dose of mRNA vaccine (BNT162b2) and inactivated SARS-CoV-2 vaccine (CoronaVac) on risks of all-cause mortality, emergency department visit, and unscheduled hospitalization. BMC Medicine. 2022;20(1):119. doi: 10.1186/s12916-022-02321-4. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1186/s12916-022-02321-4&link_type=DOI) 30. 30.Hall VG, Ferreira VH, Wood H, Ierullo M, Majchrzak-Kita B, Manguiat K, et al. Delayed-interval BNT162b2 mRNA COVID-19 vaccination enhances humoral immunity and induces robust T cell responses. Nature Immunology. 2022;23(3):380–5. doi: 10.1038/s41590-021-01126-6. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41590-021-01126-6&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=35115679&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 31. 31.Silva PJS, Sagastizábal C, Nonato LG, Struchiner CJ, Pereira T. Optimized delay of the second COVID-19 vaccine dose reduces ICU admissions. 2021;118(35):e2104640118. doi: doi:10.1073/pnas.2104640118. [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NDoicG5hcyI7czo1OiJyZXNpZCI7czoxODoiMTE4LzM1L2UyMTA0NjQwMTE4IjtzOjQ6ImF0b20iO3M6NTA6Ii9tZWRyeGl2L2Vhcmx5LzIwMjMvMDEvMzEvMjAyMi4wOS4xNC4yMjI3OTk1OS5hdG9tIjt9czo4OiJmcmFnbWVudCI7czowOiIiO30=) 32. 32.Han S, Cai J, Yang J, Zhang J, Wu Q, Zheng W, et al. Time-varying optimization of COVID-19 vaccine prioritization in the context of limited vaccination capacity. Nature Communications. 2021;12(1):4673. doi: 10.1038/s41467-021-24872-5. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41467-021-24872-5&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=34344871&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 33. 33.Matrajt L, Eaton J, Leung T, Brown ER. Vaccine optimization for COVID-19: Who to vaccinate first? 2021;7(6):eabf1374. doi: doi:10.1126/sciadv.abf1374. [FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6MzoiUERGIjtzOjExOiJqb3VybmFsQ29kZSI7czo4OiJhZHZhbmNlcyI7czo1OiJyZXNpZCI7czoxMjoiNy82L2VhYmYxMzc0IjtzOjQ6ImF0b20iO3M6NTA6Ii9tZWRyeGl2L2Vhcmx5LzIwMjMvMDEvMzEvMjAyMi4wOS4xNC4yMjI3OTk1OS5hdG9tIjt9czo4OiJmcmFnbWVudCI7czowOiIiO30=) 34. 34.Moore S, Hill EM, Dyson L, Tildesley MJ, Keeling MJ. Modelling optimal vaccination strategy for SARS-CoV-2 in the UK. PLOS Computational Biology. 2021;17(5):e1008849. doi: 10.1371/journal.pcbi.1008849. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1371/journal.pcbi.1008849&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 35. 35.Gonzalez-Parra G. Analysis of Delayed Vaccination Regimens: A Mathematical Modeling Approach. 2021;2(3):271–93. PubMed PMID: doi:10.3390/epidemiologia2030021. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.3390/epidemiologia2030021&link_type=DOI) 36. 36.Qi L, Ryan W, Sean TM, Mateo D, Emmanuel T, William B, et al. Optimization of vaccination for COVID-19 in the midst of a pandemic. Networks and Heterogeneous Media. 2022;17(3):443–66. doi: 10.3934/nhm.2022016. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.3934/nhm.2022016&link_type=DOI) 37. 37.González-Parra G, Cogollo MR, Arenas AJ. Mathematical Modeling to Study Optimal Allocation of Vaccines against COVID-19 Using an Age-Structured Population. 2022;11(3):109. PubMed PMID: doi:10.3390/axioms11030109. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.3390/axioms11030109&link_type=DOI) 38. 38.Dogra P, Butner JD, Ruiz Ramírez J, Chuang Y-l, Noureddine A, Jeffrey Brinker C, et al. A mathematical model to predict nanomedicine pharmacokinetics and tumor delivery. Computational and Structural Biotechnology Journal. 2020;18:518–31. doi: [https://doi.org/10.1016/j.csbj.2020.02.014](https://doi.org/10.1016/j.csbj.2020.02.014). 39. 39.Dogra P, Ramírez JR, Butner JD, Peláez MJ, Chung C, Hooda-Nehra A, et al. Translational Modeling Identifies Synergy between Nanoparticle-Delivered miRNA-22 and Standard-of-Care Drugs in Triple-Negative Breast Cancer. Pharmaceutical Research. 2022;39(3):511–28. doi: 10.1007/s11095-022-03176-3. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1007/s11095-022-03176-3&link_type=DOI) 40. 40.Anaya D, Dogra P, Wang Z, Haider M, Ehab J, Jeong D, et al. A Mathematical Model to Estimate Chemotherapy Concentration at the Tumor-Site and Predict Therapy Response in Colorectal Cancer Patients with Liver Metastases. Cancers 2021, 13, 444. s Note: MDPI stays neutral with regard to jurisdictional claims in published …; 2021. 41. 41.Goel S, Zhang G, Dogra P, Nizzero S, Cristini V, Wang Z, et al. Sequential deconstruction of composite drug transport in metastatic breast cancer. 2020;6(26):eaba4498. doi: doi:10.1126/sciadv.aba4498. [FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6MzoiUERGIjtzOjExOiJqb3VybmFsQ29kZSI7czo4OiJhZHZhbmNlcyI7czo1OiJyZXNpZCI7czoxMzoiNi8yNi9lYWJhNDQ5OCI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIzLzAxLzMxLzIwMjIuMDkuMTQuMjIyNzk5NTkuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 42. 42.Goel S, Ferreira CA, Dogra P, Yu B, Kutyreff CJ, Siamof CM, et al. Size-Optimized Ultrasmall Porous Silica Nanoparticles Depict Vasculature-Based Differential Targeting in Triple Negative Breast Cancer. 2019;15(46):1903747. doi: [https://doi.org/10.1002/smll.201903747](https://doi.org/10.1002/smll.201903747). 43. 43.Staquicini DI, Barbu EM, Zemans RL, Dray BK, Staquicini FI, Dogra P, et al. Targeted phage display-based pulmonary vaccination in mice and non-human primates. Med. 2021;2(3):321-42.e8. doi: [https://doi.org/10.1016/j.medj.2020.10.005](https://doi.org/10.1016/j.medj.2020.10.005). 44. 44.Dogra P, Ruiz-Ramírez J, Sinha K, Butner JD, Peláez MJ, Rawat M, et al. Innate Immunity Plays a Key Role in Controlling Viral Load in COVID-19: Mechanistic Insights from a Whole-Body Infection Dynamics Model. ACS Pharmacology & Translational Science. 2021;4(1):248–65. doi: 10.1021/acsptsci.0c00183. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1021/acsptsci.0c00183&link_type=DOI) 45. 45.Lucas C, Wong P, Klein J, Castro TBR, Silva J, Sundaram M, et al. Longitudinal analyses reveal immunological misfiring in severe COVID-19. Nature. 2020;584(7821):463–9. doi: 10.1038/s41586-020-2588-y. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41586-020-2588-y&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=32717743&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 46. 46. Collier A-rY, Yu J, McMahan K, Liu J, Chandrashekar A, Maron JS, et al. Differential Kinetics of Immune Responses Elicited by Covid-19 Vaccines. New England Journal of Medicine. 2021;385(21):2010–2. doi: 10.1056/NEJMc2115596. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1056/NEJMc2115596&link_type=DOI) 47. 47.Peeters M, Verbruggen L, Teuwen L, Vanhoutte G, Vande Kerckhove S, Peeters B, et al. Reduced humoral immune response after BNT162b2 coronavirus disease 2019 messenger RNA vaccination in cancer patients under antineoplastic treatment. ESMO Open. 2021;6(5):100274. Epub 2021/10/02. doi: 10.1016/j.esmoop.2021.100274. PubMed PMID: 34597941; PubMed Central PMCID: PMCPMC8423808. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.esmoop.2021.100274&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=34597941&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 48. 48.Liang F, Lindgren G, Lin A, Thompson EA, Ols S, Röhss J, et al. Efficient Targeting and Activation of Antigen-Presenting Cells In Vivo after Modified mRNA Vaccine Administration in Rhesus Macaques. Molecular therapy : the journal of the American Society of Gene Therapy. 2017;25(12):2635–47. Epub 2017/09/30. doi: 10.1016/j.ymthe.2017.08.006. PubMed PMID: 28958578; PubMed Central PMCID: PMCPMC5768558. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.ymthe.2017.08.006&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=28958578&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 49. 49.Milo R, Phillips R. Cell biology by the numbers: Garland Science; 2015. 50. 50.Cox RJ, Brokstad KA. Not just antibodies: B cells and T cells mediate immunity to COVID-19. Nature Reviews Immunology. 2020;20(10):581–2. doi: 10.1038/s41577-020-00436-4. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41577-020-00436-4&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=32839569&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 51. 51.Schoenmaker L, Witzigmann D, Kulkarni JA, Verbeke R, Kersten G, Jiskoot W, et al. mRNA-lipid nanoparticle COVID-19 vaccines: Structure and stability. International Journal of Pharmaceutics. 2021;601:120586. doi: [https://doi.org/10.1016/j.ijpharm.2021.120586](https://doi.org/10.1016/j.ijpharm.2021.120586). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.ijpharm.2021.120586&link_type=DOI) 52. 52.Abbas AK, Lichtman AH, Pillai S. Basic immunology e-book: functions and disorders of the immune system: Elsevier Health Sciences; 2019. 53. 53.Ahlmann M, Hempel G. The effect of cyclophosphamide on the immune system: implications for clinical cancer therapy. Cancer Chemother Pharmacol. 2016;78(4):661–71. doi: 10.1007/s00280-016-3152-1. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1007/s00280-016-3152-1&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 54. 54.Epstein RS, Aapro MS, Basu Roy UK, Salimi T, Krenitsky J, Leone-Perkins ML, et al. Patient Burden and Real-World Management of Chemotherapy-Induced Myelosuppression: Results from an Online Survey of Patients with Solid Tumors. Adv Ther. 2020;37(8):3606–18. Epub 2020/07/10. doi: 10.1007/s12325-020-01419-6. PubMed PMID: 32642965; PubMed Central PMCID: PMCPMC7340862. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1007/s12325-020-01419-6&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=32642965&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 55. 55.Epstein RS, Nelms J, Moran D, Girman C, Huang H, Chioda M. Treatment patterns and burden of myelosuppression for patients with small cell lung cancer: A SEER-medicare study. Cancer Treatment and Research Communications. 2022;31:100555. doi: [https://doi.org/10.1016/j.ctarc.2022.100555](https://doi.org/10.1016/j.ctarc.2022.100555). 56. 56.Atwal D, Joshi KP, Ravilla R, Mahmoud F. Pembrolizumab-Induced Pancytopenia: A Case Report. Perm J. 2017;21:17–004. Epub 2017/07/27. doi: 10.7812/tpp/17-004. PubMed PMID: 28746020; PubMed Central PMCID: PMCPMC5528855. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.7812/TPP/17-024&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=28746020&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 57. 57.Ueki Y, Suzuki M, Horikawa Y, Watanabe H, Yamaguchi Y, Morita C, et al. Pembrolizumab-induced pancytopenia in a patient with squamous cell lung cancer. Thoracic Cancer. 2020;11(9):2731–5. doi: [https://doi.org/10.1111/1759-7714.13582](https://doi.org/10.1111/1759-7714.13582). 58. 58.Diao B, Wang C, Tan Y, Chen X, Liu Y, Ning L, et al. Reduction and Functional Exhaustion of T Cells in Patients With Coronavirus Disease 2019 (COVID-19). 2020;11. doi: 10.3389/fimmu.2020.00827. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.3389/fimmu.2020.00827&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 59. 59.Galbraith MD, Kinning KT, Sullivan KD, Araya P, Smith KP, Granrath RE, et al. Specialized interferon action in COVID-19. 2022;119(11):e2116730119. doi: doi:10.1073/pnas.2116730119. [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NDoicG5hcyI7czo1OiJyZXNpZCI7czoxODoiMTE5LzExL2UyMTE2NzMwMTE5IjtzOjQ6ImF0b20iO3M6NTA6Ii9tZWRyeGl2L2Vhcmx5LzIwMjMvMDEvMzEvMjAyMi4wOS4xNC4yMjI3OTk1OS5hdG9tIjt9czo4OiJmcmFnbWVudCI7czowOiIiO30=) 60. 60.Vatansever HS, Becer E. Relationship between IL-6 and COVID-19: to be considered during treatment. 2020;15(12):817–22. doi: 10.2217/fvl-2020-0168. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.2217/fvl-2020-0168&link_type=DOI) 61. 61.Perelson AS, Ke R. Mechanistic Modeling of SARS-CoV-2 and Other Infectious Diseases and the Effects of Therapeutics. Clin Pharmacol Ther. 2021;109(4):829–40. Epub 2021/01/08. doi: 10.1002/cpt.2160. PubMed PMID: 33410134; PubMed Central PMCID: PMCPMC8142935. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/cpt.2160&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=33410134&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 62. 62.Bayart J-L, Douxfils J, Gillot C, David C, Mullier F, Elsen M, et al. Waning of IgG, Total and Neutralizing Antibodies 6 Months Post-Vaccination with BNT162b2 in Healthcare Workers. Vaccines. 2021;9(10):1092. PubMed PMID: doi:10.3390/vaccines9101092. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.3390/vaccines9101092&link_type=DOI) 63. 63.Kontopoulou K, Nakas CT, Papazisis G. Significant Increase in Antibody Titers after the 3rd Booster Dose of the Pfizer–BioNTech mRNA COVID-19 Vaccine in Healthcare Workers in Greece. 2022;10(6):876. PubMed PMID: doi:10.3390/vaccines10060876. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.3390/vaccines10060876&link_type=DOI) 64. 64.Kontopoulou K, Nakas CT, Papazisis G. Significant Increase in Antibody Titers after the 3rd Booster Dose of the Pfizer–BioNTech mRNA COVID-19 Vaccine in Healthcare Workers in Greece. Vaccines. 2022;10(6):876. PubMed PMID: doi:10.3390/vaccines10060876. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.3390/vaccines10060876&link_type=DOI) 65. 65.Cromer D, Steain M, Reynaldi A, Schlub TE, Wheatley AK, Juno JA, et al. Neutralising antibody titres as predictors of protection against SARS-CoV-2 variants and the impact of boosting: a meta-analysis. The Lancet Microbe. 2022;3(1):e52–e61. doi: 10.1016/S2666-5247(21)00267-6. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S2666-5247(21)00267-6&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=34806056&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 66. 66.Khoury DS, Cromer D, Reynaldi A, Schlub TE, Wheatley AK, Juno JA, et al. Neutralizing antibody levels are highly predictive of immune protection from symptomatic SARS-CoV-2 infection. Nature Medicine. 2021;27(7):1205–11. doi: 10.1038/s41591-021-01377-8. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41591-021-01377-8&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=34002089&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 67. 67.Goldblatt D, Fiore-Gartland A, Johnson M, Hunt A, Bengt C, Zavadska D, et al. Towards a population-based threshold of protection for COVID-19 vaccines. Vaccine. 2022;40(2):306–15. doi: [https://doi.org/10.1016/j.vaccine.2021.12.006](https://doi.org/10.1016/j.vaccine.2021.12.006). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.vaccine.2021.12.006&link_type=DOI) 68. 68.WHO. Vaccine efficacy, effectiveness and protection: WHO; 2021 [1/6/2023]. Available from: [https://www.who.int/news-room/feature-stories/detail/vaccine-efficacy-effectiveness-and-protection](https://www.who.int/news-room/feature-stories/detail/vaccine-efficacy-effectiveness-and-protection). 69. 69.Siber GR. Methods for estimating serological correlates of protection. Developments in biological standardization. 1997;89:283–96. Epub 1997/01/01. PubMed PMID: 9272362. [PubMed](http://medrxiv.org/lookup/external-ref?access_num=9272362&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=A1997BJ43T00033&link_type=ISI) 70. 70.Greaney AJ, Starr TN, Bloom JD. An antibody-escape estimator for mutations to the SARS-CoV-2 receptor-binding domain. Virus Evolution. 2022;8(1). doi: 10.1093/ve/veac021. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/ve/veac021&link_type=DOI) 71. 71.Harvey WT, Carabelli AM, Jackson B, Gupta RK, Thomson EC, Harrison EM, et al. SARS-CoV-2 variants, spike mutations and immune escape. Nature Reviews Microbiology. 2021;19(7):409–24. doi: 10.1038/s41579-021-00573-0. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41579-021-00573-0&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=34075212&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 72. 72.Dogra P, Ramírez JR, Butner JD, Peláez MJ, Chung C, Hooda-Nehra A, et al. Translational Modeling Identifies Synergy between Nanoparticle-Delivered miRNA-22 and Standard-of-Care Drugs in Triple-Negative Breast Cancer. Pharmaceutical Research. 2022. doi: 10.1007/s11095-022-03176-3. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1007/s11095-022-03176-3&link_type=DOI) 73. 73.Dogra P, Butner JD, Ramírez JR, Chuang Y-l, Noureddine A, Brinker CJ, et al. A mathematical model to predict nanomedicine pharmacokinetics and tumor delivery. Computational and Structural Biotechnology Journal. 2020;18:518–31. 74. 74.Wang Z, Bordas V, Deisboeck T. Identification of critical molecular components in a multiscale cancer model based on the integration of Monte Carlo, resampling, and ANOVA. Frontiers in physiology. 2011;2:35. 75. 75.Wang Z, Deisboeck TS, Cristini V. Development of a sampling-based global sensitivity analysis workflow for multiscale computational cancer models. IET systems biology. 2014;8(5):191–7. 76. 76.Butner JD, Dogra P, Chung C, Ruiz-Ramírez J, Nizzero S, Plodinec M, et al. Dedifferentiation-mediated stem cell niche maintenance in early-stage ductal carcinoma in situ progression: insights from a multiscale modeling study. Cell death & disease. 2022;13(5):485. Epub 2022/05/22. doi: 10.1038/s41419-022-04939-x. PubMed PMID: 35597788; PubMed Central PMCID: PMCPMC9124196 declare no competing interests. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41419-022-04939-x&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=35597788&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 77. 77.Lauer SA, Grantz KH, Bi Q, Jones FK, Zheng Q, Meredith HR, et al. The Incubation Period of Coronavirus Disease 2019 (COVID-19) From Publicly Reported Confirmed Cases: Estimation and Application. Annals of internal medicine. 2020;172(9):577–82. Epub 2020/03/10. doi: 10.7326/m20-0504. PubMed PMID: 32150748; PubMed Central PMCID: PMCPMC7081172. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.7326/M20-0504&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=32150748&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 78. 78.Kontopoulou K, Nakas C, Papazisis G. Significant increase of antibody titers after the 3rd booster dose of the Pfizer–BioNTech mRNA COVID-19 vaccine in healthcare workers in Greece. Available at SSRN. 2021. 79. 79.Planas D, Saunders N, Maes P, Guivel-Benhassine F, Planchais C, Buchrieser J, et al. Considerable escape of SARS-CoV-2 Omicron to antibody neutralization. Nature. 2022;602(7898):671–5. doi: 10.1038/s41586-021-04389-z. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41586-021-04389-z&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=35016199&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F01%2F31%2F2022.09.14.22279959.atom) 80. 80.Chitu V, Yeung YG, Yu W, Nandi S, Stanley ERJCPiI. Measurement of Macrophage Growth and Differentiation. 2011;92. 81. 81.De Boer RJ, Homann D, Perelson AS. Different Dynamics of CD4+ and CD8+ T Cell Responses During and After Acute Lymphocytic Choriomeningitis Virus Infection. 2003;171(8):3928–35. doi: 10.4049/jimmunol.171.8.3928 %J The Journal of Immunology. [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6ODoiamltbXVub2wiO3M6NToicmVzaWQiO3M6MTA6IjE3MS84LzM5MjgiO3M6NDoiYXRvbSI7czo1MDoiL21lZHJ4aXYvZWFybHkvMjAyMy8wMS8zMS8yMDIyLjA5LjE0LjIyMjc5OTU5LmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 82. 82.Maizel AL, Morgan JW, Mehta SR, Kouttab NM, Bator JM, Sahasrabuddhe CG. Long-term growth of human B cells and their use in a microassay for B-cell growth factor. 1983;80(16):5047–51. doi: 10.1073/pnas.80.16.5047 %J Proceedings of the National Academy of Sciences. [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NDoicG5hcyI7czo1OiJyZXNpZCI7czoxMDoiODAvMTYvNTA0NyI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIzLzAxLzMxLzIwMjIuMDkuMTQuMjIyNzk5NTkuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) [1]: /embed/graphic-3.gif [2]: /embed/inline-graphic-1.gif [3]: /embed/graphic-4.gif [4]: /embed/graphic-5.gif [5]: /embed/inline-graphic-2.gif [6]: /embed/inline-graphic-3.gif [7]: /embed/inline-graphic-4.gif [8]: /embed/graphic-6.gif [9]: /embed/graphic-7.gif [10]: /embed/graphic-8.gif [11]: /embed/graphic-9.gif [12]: /embed/graphic-10.gif [13]: /embed/graphic-11.gif [14]: /embed/graphic-12.gif [15]: /embed/graphic-13.gif [16]: /embed/graphic-14.gif [17]: /embed/graphic-15.gif [18]: /embed/graphic-16.gif [19]: /embed/graphic-17.gif [20]: /embed/graphic-18.gif [21]: /embed/graphic-19.gif [22]: /embed/graphic-20.gif [23]: /embed/inline-graphic-5.gif [24]: /embed/inline-graphic-6.gif [25]: /embed/graphic-21.gif [26]: /embed/graphic-22.gif [27]: /embed/inline-graphic-7.gif [28]: /embed/inline-graphic-8.gif [29]: /embed/inline-graphic-9.gif