A modeling-based approach to optimize COVID-19 vaccine dosing schedules for improved protection

While the development of different vaccines slowed the dissemination of SARS-CoV-2, the occurrence of breakthrough infections has continued to fuel the COVID-19 pandemic. To secure at least partial protection in the majority of the population through 1 dose of a COVID-19 vaccine, delayed administration of boosters has been implemented in many countries. However, waning immunity and emergence of new variants of SARS-CoV-2 suggest that such measures may induce breakthrough infections due to intermittent lapses in protection. Optimizing vaccine dosing schedules to ensure prolonged continuity in protection could thus help control the pandemic. We developed a mechanistic model of immune response to vaccines as an in silico tool for dosing schedule optimization. The model was calibrated with clinical data sets of acquired immunity to COVID-19 mRNA vaccines in healthy and immunocompromised participants and showed robust validation by accurately predicting neutralizing antibody kinetics in response to multiple doses of COVID-19 mRNA vaccines. Importantly, by estimating population vulnerability to breakthrough infections, we predicted tailored vaccination dosing schedules to minimize breakthrough infections, especially for immunocompromised individuals. We identified that the optimal vaccination schedules vary from CDC-recommended dosing, suggesting that the model is a valuable tool to optimize vaccine efficacy outcomes during future outbreaks.


S1. Mathematical model description
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 (1). Given that the timescale of drug delivery (intramuscular injection) and mRNA translation is much shorter (< 1 hour) (2) than that of the vaccine-induced immune response (days to weeks) (3), we assumed that the variable ሺ ሻ 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
where Dose indicates the dimensionless dose of the antigen administered via the vaccine. The concentration of the spike protein ሺ ሻ is described by the sum of Gaussians centered at , which represents the day on which a vaccine dose is injected out of the set of doses indicated by . NP is the characteristic time of clearance of the antigen-carrying nanoparticle (NP) from the body (4), estimated based on NP diameter of 100 nm for mRNA vaccines (5).
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 APC തതതതത 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 APC as they migrate towards the lymphoid tissue. The APC activation process is proportional to the antigen load (Agሺ ሻ), which can be derived either from the vaccine or natural infection and is either equal to ሺ ሻ or the viral load ሺ ሻ in the case of vaccination or infection, respectively, with 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 (
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. S2, and a death term determined by the death rate constant APC of activated APCs. Of note, in our model we have included a dimensionless coefficient ∈ [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 = 1 in healthy individuals, and < 1 in immunocompromised patients. Since our model is calibrated for patients who are immunocompromised due to anticancer therapy, immunosuppression in our model is characterized by T-and B-cell deficiency, which is one of the important immunological effects observed due to disruption of hematopoiesis leading to myelosuppression in patients undergoing anticancer therapy (7)(8)(9)(10)(11). 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 (12) 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 IL6 is the Michaelis constant for IL-6 effects.

Equation for the effector CD4+ T-cell density (
where is the death rate of effector T-cells.

Equation for the naïve B cell density ( ሺ ሻ)
where 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 BC is the differentiation rate of B cells into plasma cells.

Equation for the activated B cell density (
where BC is the differentiation rate of B cells into plasma cells.

Equation for the plasma cell density ( ሺ ሻ)
where 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 Ab , is proportional to the plasma cell density. The antibodies secreted into the plasma are then cleared at a rate Cl Ab , which is a lumped phenomenological parameter characterizing the various antibody clearance mechanisms.

Equation for the neutralizing antibody concentration ( ሺ ሻ)
Abሺ ሻ = Ab ⋅ ሺ ሻ ⏞ − Cl Ab ⋅ Abሺ ሻ ⏞ , Abሺ0ሻ = 0 (S11) 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 immunoregulatory elements (13,14). Each cytokine has a unique source and key role in the immune T-cells, effector CD8+ T-cells, and active APCs, tends to exhaust naïve CD4+ and CD8+ T-cell populations (12). The rate of change of cytokine concentration was modeled using production term and degradation terms, 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 (
where IFN1 is the production rate of type-I interferon.

Equation for the type-II interferon concentration (
where IFN2 is the production rate of type-II interferon.

Equation for the interleukin-6 concentration ( ሺ ሻ)
where, 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 natural infection. In the latter case, healthy susceptible cells are transformed 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 (15), as described by the equations below: Equation for the healthy respiratory epithelial cell density ( ሺ ሻ) where is the viral infectivity rate, ሺ ሻ is the viral load density in the respiratory tract, and 0 is the initial density of healthy cells.

Equation for the density of infected cells in the respiratory tract epithelium ( ሺ ሻ)
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 * ሺ ሻ is the density of effector CD8+ T-cells.

Equation for the viral load density in the respiratory tract ( ሺ ሻ)
where represents virion production rate, IFN1ሺ ሻ is the concentration of type-I interferons, IFN1 is the Michaelis constant of the virion production suppression factor, Ab is the antibody-mediated neutralization rate of viruses, Abሺ ሻ is the antibody concentration in the body, APCሺ ሻ is the density of naïve APCs in the respiratory tract, APC is the naïve APC-mediated clearance rate of viruses, and 0 is the initial viral load at the time of infection. Pearson correlation analysis (model calibration) Figure S3. Pearson correlation analysis between model fits and experimental data for immune response dynamics (comprising one or more of the following variables: viral load, naïve and effector CD4+ and CD8+ T-cells, type-I and type-II interferon, IL-6, and neutralizing antibody (IgG)) during A) moderate SARS-CoV-2 infection, vaccination in B) healthy individuals, C) cancer patients receiving chemotherapy, and D) cancer patients receiving immunotherapy. R value represents Pearson correlation coefficient.

Results
Pearson Correlation analysis (model validation) Figure S4. Pearson correlation analysis between model predictions and experimental data pooled from clinical studies on antibody (IgG) response dynamics following two doses of Moderna COVID-19 mRNA vaccine, and two and three doses of Pfizer-BioNTech COVID-19 mRNA vaccine. R value represents Pearson correlation coefficient. Figure S5. Local sensitivity analysis exhibiting correlation between parameter perturbation and sensitivity index. Parameters were perturbed linearly between ±50% around the baseline value, except that was perturbed between -50% to baseline.