(M1530-11-62) Mathematical Modeling to Optimize Dosing Schedules for COVID-19 Vaccines

Purpose: While the development of vaccines has slowed down the progression of the COVID-19 pandemic, the occurrence of breakthrough infections due to waning immunity and emergence of variants of concern (VOCs) continues to fuel the pandemic. The problem is further aggravated by logistic constraints preventing equitable global distribution of vaccines, especially to low-income countries, thereby jeopardizing the attainment of herd immunity. To impart at least partial protection with a single dose of vaccines to maximum possible fraction of the population, delayed administration of subsequent primary doses and booster shots has been suggested as a strategy by policy makers. To optimize vaccine dosing schedules, we have developed a mechanistic mathematical model of innate and adaptive immune response to vaccines.

Methods: We developed a mathematical model of innate and adaptive immune response to COVID-19 vaccines, as an adaptation of our previous model of immune response to viral infections [1]. The model incorporates key biological processes that are relevant to antigen presentation at the site of vaccination (i.e., muscle), the development of innate and adaptive immune responses in the lymphoid tissue, and the protection against infection in the respiratory tract, formulated as a system of ordinary differential equations. The model was calibrated and validated with clinical data for mRNA-based COVID-19 vaccines to study the effect of immune health status and dosing schedules on plasma antibody titers (a correlate protection against infection) leading to the estimation of population vulnerability to breakthrough infections. The model was then used to identify optimal schedules for the primary series and booster doses to minimize breakthrough infections in the population.

Results: We began by fitting the model to immune response kinetics of SARS-CoV-2 infection [2], which allowed us to estimate several unknown model parameters relevant to key immune response variables (Figure 1a). Subsequently, we calibrated the model with the clinical data for effector lymphocytes and neutralizing antibodies (Figure 1b), obtained from healthy individuals vaccinated with a double dose of mRNA vaccines (Pfizer) [3]. Finally, to accurately represent the vaccine-induced immune response in immunocompromised individuals, we also calibrated the immune health status parameter f with clinical data obtained from vaccinated cancer patients undergoing chemotherapy or immunotherapy (Figure 1c,d) [4]. Further, to evaluate the effects of dosing schedules and the immune health status, at the population level, on the variability in immune response to mRNA vaccines, we simulated the vaccination of a virtual population of 10,000 individuals with 3 doses and assessed the corresponding vulnerability to breakthrough infections. As shown in Figure 2a, the ‘average’ antibody kinetics across the 10,000 individuals stayed above the protective threshold for omicron (OM) and wildtype (WT) strains. However, for a significant fraction of the population, antibody levels remained below the OM threshold for a prolonged period (~5 months). Further, translating the antibody levels to vaccine efficacy, we observed that for a significant fraction of the 10,000 individuals, vaccine efficacy against OM fell below the 50% protection threshold (Figure 2b, 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, i.e., population at high risk of breakthrough infection (Figure 2c). Lastly, we intended to identify optimal vaccine dosing schedules to achieve continuous protection against the OM strain for prolonged periods. We generated three virtual cohorts of 10,000 individuals each to represent healthy, mildly immunocompromised, and highly immunocompromised individuals, and implemented several dosing schedules to identify optimal dosing times in each sub-population. As shown in Figure 3, the area under the curve (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), indicating optimal waiting periods between doses that minimize breakthrough infections. As shown in Figure 3a,d,g, as the immune status goes from healthy to highly immunocompromised, the position of the minima on the x-axis shows a right shift. In contrast, as shown in Figure 3b-i, the minima for the 1^st and 2^nd booster dose shows a left shift on the x-axis going from healthy to highly immunocompromised individuals. Through our analysis we highlighted critical waiting windows for immunocompromised individuals (26 days and 31 days after first dose for mildly and highly immunocompromised individuals, respectively) to ensure sufficient recovery of naïve B cells and minimize vulnerability to breakthrough infections in their sub-population.

Conclusion: The presented model is based on generalized innate and 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 optimize vaccine dosing schedules in case of future outbreak events.

References: 1. 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. 2020. doi: 10.1021/acsptsci.0c00183.
2. 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.
3. 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.
4. 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.

Acknowledgments: 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), National Science Foundation Grant DMS-1930583 (ZW, VC), the National Institutes of Health (NIH) Grants 1R01CA253865 (ZW, VC), 1R01CA226537 (ZW, VC), 1R01CA222007 (ZW, VC), and 1R01AI165372 (ZW). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Figure 1. Model calibration. Model calibration with literature-derived clinical data of immune system response kinetics during a) SARS-CoV-2 infection and b,c,d) vaccination in b) healthy individuals, c) cancer patients receiving chemotherapy, and d) cancer patients receiving immunotherapy. For consistency, all immunization data was based on double dose of Pfizer mRNA vaccine. Solid or dashed lines indicate model simulations; markers with errorbars represent mean±SD values of the data.



Figure 2. Effect of heterogeneity in vaccine dosing schedules and immune health status on breakthrough infections at the population scale. a) Average antibody levels in the plasma, b) corresponding vaccine or antibody efficacy, and c) population fraction vulnerable to breakthrough infections due to wildtype (WT, solid blue line) and omicron strain (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 one SD. Note that dose 1 was given on day 0 to each simulated individual, dose 2 was given between day 14 and day 56 after dose 1, and booster 1 was given between days 150 and 270 days after dose 2. Red and blue brackets on x-axis denote timing windows with respect to day 0 for second dose and first booster, respectively, used to design unique vaccine schedules in model simulations. Immune health status (f) of the simulated population varied between 0.5 to 1.



Figure 3. 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, and c,f,i) booster 2 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.