Vaccination compartmental epidemiological models for the delta and omicron SARS-CoV-2 variants

We explore the inclusion of vaccination in compartmental epidemiological models concerning the delta and omicron variants of the SARS-CoV-2 virus that caused the COVID-19 pandemic. We expand on our earlier compartmental-model work by incorporating vaccinated populations. We present two classes of models that differ depending on the immunological properties of the variant. The first one is for the delta variant, where we do not follow the dynamics of the vaccinated individuals since infections of vaccinated individuals were rare. The second one for the far more contagious omicron variant incorporates the evolution of the infections within the vaccinated cohort. We explore comparisons with available data involving two possible classes of counts, fatalities and hospitalizations. We present our results for two regions, Andalusia and Switzerland (including the Principality of Liechtenstein), where the necessary data are available. In the majority of the considered cases, the models are found to yield good agreement with the data and have a reasonable predictive capability beyond their training window, rendering them potentially useful tools for the interpretation of the COVID-19 and further pandemic waves, and for the design of intervention strategies during these waves.


I. INTRODUCTION
Over the last two and a half years, the COVID-19 pandemic has been deemed responsible, to date, for 760 million confirmed cases, and over 6.8 million deaths worldwide, as of this writing and according to the World Health Organization COVID-19 dashboard.As such, its emergence wreaked havoc in life as we knew it throughout the world and forced a dramatic modification of our social and economic activities during this interval.At the same time, it triggered a global mobilization of the scientific community to produce vaccines rapidly, especially through (thankfully, by that time, fairly mature) technology of mRNA-based methods.This effort led to the remarkable result of having a vaccine against SARS-CoV-2 within a year of its emergence.Nevertheless, this was far from the end of the story, as new variants of the SARS-CoV-2 virus kept emerging within 2020 and 2021.The so-called delta variant appeared in India in late 2020 and it had spread to 179 countries by November 2021.Subsequently, the delta variant was superseded by the so-called omicron variant that was reported in South Africa on November 2021, and subsequently became rapidly the predominant variant of SARS-CoV-2 thereafter.
The theoretical and mathematical modeling of infectious diseases such as COVID-19 has a long and time-honored history since the classic work of Kermack and McKendrick [1].Moreover, relevant efforts have been summarized in numerous venues in recent years, such as, e.g., [2][3][4], to mention only a few.The urgency and severity of the COVID-19 pandemic brought about an intense effort on the side of the mathematical and physical communities to develop analytical models and computational tools that could be used to examine the unprecedented volume of available data regarding the temporal (and spatial) evolution of the pandemic and to make predictions for the weeks (or in some cases month(s)) ahead.A notable example of comparison of such efforts can be seen in, e.g., websites such as [5].Relevant modeling efforts have now been summarized in a number of reviews such as [6,7], including ones of specialized modeling aspects such as the study of metapopulation network models [8], while other works summarized the challenges and difficulties of associated modeling [9,10].
Over the last year, a large portion of the focus of the modeling efforts has shifted towards the inclusion of vaccination in epidemniological models.While a lot of information is available regarding the effectiveness and efficacy of vaccines [11] (see also websites such as [12]) mathematical models can still be quite useful in a number of ways, including in guiding and informing distribution strategies thereof [13,14].It is in that light that numerous compartmental epidemiological models with vaccination strategies have arisen in the literature [15,16], including some specific to different geographical locations [17] and to different social infrastructures, such as nursing homes [18].While the relevant models feature different levels of complexity starting from SIRV (Susceptible-Infected-Recovered-Vaccinated) extensions of the classic SIR (Susceptible-Infected-Recovered) [15] model and progressively extending to multicomponent models such as, e.g., [19], our aim here is to build systematically on the earlier modeling attempt of [20] by considering SARS-CoV-2 variants that affect differently the vaccinated population.
More concretely, our aim is to present a model of the omicron variant (model A and its two implementations A1 and A2), in which its highly contagious nature allows for so-called breakthrough infections, whereby vaccinated individuals may still be infected.In that light, we account for the standard populations of our earlier work [20], including exposed, pre-symptomatics, and asymptomatics [21], as well as hospitalizations, recoveries and fatalities.In addition, we consider such populations in both the unvaccinated and vaccinated portions of the population and their interactions.The primary aim of the associated study is to explore the dynamical evolution of the omicron variant from the end of 2021 to early spring 2022.We also present a simpler model (model B) for the evolution of the earlier delta variant during the fall of 2021.In that case, vaccination was deemed to protect individuals from being infected, and the fraction of breakthrough infections was quite small, even in groups such as the potentially highly exposed group of healthcare workers [22].Accordingly, we assume that the vaccinated population may be effectively removed from the susceptible compartment.
We examine two versions of the proposed omicron model in Section II (models A1 and A2), their difference motivated by data availability.The population flows in Model A1 terminate at the fatalities compartment: as such, the model considers that fatalities in both the unvaccinated and vaccinated populations provide the most reliable data.Model A2 is motivated by the existence of systematic data for the total number of hospitalizations (conventional and ICU): here, population flows terminate at the hospitalizations compartment, i.e., they do not branch further to the fatalities compartment as in model A1.This is for a number of reasons: in primis, reporting of fatalities occasionally occurs retroactively (and less reliably).Our models are applied to two regions with similar populations (approximately 8 million inhabitants): Andalusia and Switzerland (including the Principality of Liechtenstein) The motivation for this choice arises, once again, from the availability of suitably stratified data, whereby both fatalities and hospitalizations are available for vaccinated and unvaccinated individuals.In Section III we propose, in addition to the more detailed model for the omicron variant presented in Section II, a model for the delta variant, model B. We use model B in the same spatial regions.We typically find that numerical results compare favorably to available data, both in terms of the comparison of the regression results and also in connection to testing beyond the end of the training period for the model parameters.Finally, in Section IV, we summarize our findings and present our conclusions.In the Appendix, we consider mathematically the question of structural identifiability of the models developed herein.

A. Model A1: Branches terminate at fatalities
The first model for the omicron variant, model A1, extends our compartmental epidemiological model used to examine the COVID-19 pandemic evolution in Mexico [20].Accordingly, the susceptible population S can turn to exposed (E) through interactions with either symptomatically infected (I), presymptomatic (P ), or asymptomatic (A) individuals.The exposed population E, in turn, can convert to either P , within a time scale 1/σ 1 , leading to different clinical stages of the disease or to A, a compartment that has been recognized to play a key role in the dynamical evolution of COVID-19 [21].The asymptomatic population A can only lead to undisclosed recoveries (denoted as U ), over a time scale 1/µ.On the other hand, the presymptomatic individuals P turn to infected with clinical symptoms I over a time scale 1/σ 2 .The addition of the latency period 1/σ 1 and the preclinical period 1/σ 2 constitute the incubation time scale of the disease, τ inc = 1/σ 1 + 1/σ 2 .Subsequently, the symptomatically infected can either turn to hospitalized H at a rate γ h , while the rest may recover (R) at a rate γ r .Finally, those in the hospitalized population of H can, again, branch into two populations: they either recover at a rate κ r , or they lead to fatalities (D) at a rate κ d .

Omicron Model A1
Population flows

Omicron Models A1 and A2
Susceptible interactions While these populations were also present in our earlier work [20], it is relevant to highlight the differences of the omicron-variant modeling.For the period under consideration (fall 2021 to spring 2022), vaccines had been deployed extensively in Andalusia and Switzerland.More importantly, breakthrough infections due to the omicron variant were substantial within the vaccinated population (contrary to the case of the delta variant considered in Section III).In light of that, we formulated two sets of populations: one representing unvaccinated individuals, denoted by (the subscript) u, and the other representing the substantial population of vaccinated individuals, denoted by (the subscript) v.Each population subset had its own set of parameters.This division of the whole population into two almost independent subgroups (they interact only through the vaccinated time series) is reflected in the presentation of the model schematic in Fig. 1.In fact, Fig. 1 summarizes population flows and interactions: the left panel illustrates the main population compartments and the corresponding flows, while the right panel shows the interactions of susceptibles with other compartments, both for vaccinated and unvaccinated populations.The model equations are: We made a number of simplifying assumptions to reduce the number of parameters and enhance the identifiability of the model (see also the relevant analysis in Appendix A).We consider a model with only four transmission rates β ij (i, j = u, v).We assumed that infectious contacts could only occur between four groups: between unvaccinated individuals (unvaccinatedunvaccinated contacts denoted by the superscript uu), between vaccinated and vaccinated individuals (denoted by the superscript vv) and across these two groups (denoted by uv -for vaccinated transmitting to unvaccinated and vu for the reverse path of infection).Notice that uv and vu are not a priori assumed to be equivalent.Within each subgroup of infection transmission (uu, vv, uv, and vu), the infections induced by the three infectious compartments P , A and I are assumed to occur at the same rate, i.e., the transmission rate is considered to be independent of whether the infectious individual exhibits symptoms (I) or not (A, P ).While we expect these transmission rates to differ (in fact, we know that even within a given population the S − I transmission rate differs from the S − A transmission rate, see for example, Ref. [20]), the identifiability analysis based on the available time series suggests that they would not be independently computable in a definitive way.In addition, we introduced a single constraint that requires that the incubation period of the disease τ inc = σ −1 1 + σ −1 2 be a value randomly sampled from a normal distribution with mean 3.42 and standard deviation 0.2755.This leverages information about the (shorter) incubation period associated with the omicron variant [23].
Vaccine efficiency is introduced via the parameter θ, whose variation bounds were set in the range 75%-95%.This factor multiplied by the time series of vaccinations V (t) effectively "transfers" individuals from the unvaccinated susceptible population to the vaccinated susceptible compartment.It is important to remark that H measures both conventional and critical hospitalizations together.While we recognize the relevance of the ongoing debate of distinguishing deaths "from COVID" vs. "with COVID" [24], unfortunately the data available herein do not allow for a definitive distinction between the two.
We obtained the best-fit parameters and initial conditions by minimizing an appropriately chosen norm.For both regions of interest, the time period used for the fits was from November 15, 2021 to March 1, 2022.The identification of the date a particular variant appeared in a geographical location is fraught with uncertainties.The choice of November 15, 2021 as the initial day of fittings stems from a number of indirect indications: Ref.
[25] reports a surge of cases in Germany at the beginning of November; Ref. [26] mentions that an omicron-variant case was reported on November 19; and the WHO site [27] mentions that in South Africa the first confirmed infection was reported on November 24, although arising in the sequencing of a sample collected on November 9. Additionally, inspection of the data shows a gradual increase starting at November 15, after a plateau.The effective parameter training period indicated above (till March 1, 2022), is followed by a prediction period (with the optimal parameters and initial conditions fixed, as determined in the training period).The predicted time series that terminates on March 29, 2022 is then compared to the reported data.Predictions do not go beyond that date because the measurement strategy in Andalusia changed, thereby rendering our fixed parameters of limited relevance to the new data.Moreover, around that time Spanish public policy also changed, and face masks were no longer required.In Switzerland some restrictions were removed in the middle of February.More details are presented in the appropriate results sections.
We perform two separate fits, i.e., we used two different norms to compare predictions to reported numbers depending on data availability.First, we fit the predicted total number of fatalities to the reported number by minimizing the norm N (i.e., the loss function) where the subscript "num" refers to predicted (calculated) numbers and "obs" to observation (reported numbers), and n refers to the number (days) of observations.This loss function effectively does not distinguish the compartmental origin of the fatalities, i.e., whether they arise from the u or v compartments: it only accounts for the cumulative number of fatalities.This will, inevitably, result in the determination of some parameters between the unvaccinated and the vaccinated populations that may not necessarily be epidemiologically meaningful (as we will see in the detailed comparisons of our predictions for Andalusia and Switzerland).Whenever we used norm (2), we also included the waning effect of the vaccines and booster vaccination effects, in addition to those fully vaccinated.It is well-documented that different vaccines have different waning immunities (see, for instance, the detailed analysis of Ref. [28]).However, to be able to account for these effects without adding a large number of additional coefficients, we assumed that vaccines are roughly effective for an interval of about 180 days.Consequently, we define V (t) as: where the subscripts f v refers to fully vaccinated, and b to booster.The above optimization via norm ( 2) is a point estimator, that is, a single set of parameters and initial conditions is obtained.To calculate their confidence intervals, and consequently the confidence interval of the predictions, we follow the bootstrapping method described in [29].The first step is to generate 250 random, synthetic, time series for the fatalities based on the reported data.To accomplish this, we first apply the optimization to find the best fit to the original data set: we refer to that optimization of the reported fatalities data as the "numerical truth".In this first optimization, we also included and H v (0) as fitting parameters: these parameters were fixed in the subsequent bootstrapping steps.Then, random noise of a prescribed level, empirically chosen to be 5%, was added to the "numerical truth" to generate 250 "polluted" (i.e., noisy) time series for the fatalities.Knowledge of the error in data collection may be helpful to select an appropriate noise level.In the second step, we apply the optimization procedure to find the best fit to each of the 250 synthetic fatalities time series to obtain 250 sets of parameters from which the confidence intervals for the parameters and predictions can be computed.The same bootstraping procedure was used for the hospitalization time series.
We also fitted separately, if the reported data allowed us, the vaccinated and unvaccinated fatalities time series using them as separate inputs to our minimization objective.In that case, the relevant norm is With this norm, we are genuinely treating the vaccinated compartment separately: we expect its fraction of fatalities (proportionally to the corresponding susceptible population) to be reflected in the obtained parameters.It should be added that in this case, given the way that the data are obtained, the vaccinated status corresponds to people who had received the full doses, independently of antibodies waning, boosting, or efficacy of vaccines.Consequently, θ has been fixed to 1 in every fit, and V (t) is defined as We mention here that an important consideration pertinent to the model concerns the identifiability of its coefficients (and initial conditions).This pertains to whether, based on the time series given, the unknown model parameters can be uniquely identified [30].In addition to the question whether all parameters can be uniquely identified (global identifiability) or some may have multiple possible values (local identifiability), there are also practical issues concerning whether different sets of parameters lead to similar (although not necessarily identical) observations; see, e.g., the discussion of [31].Here, following the approach presented in Appendix A (see also [32][33][34]), we find that all the parameters are globally identifiable, except For these eight parameters, the following combinations are globally identifiable: As concerns the initial conditions, H u (0) and H v (0) are not identifiable, in addition to the initial conditions for the terminal compartments FIG. 2. Omicron-variant model A1: Fit and prediction for the total number of fatalities in Andalusia (norm of Eq. ( 2)).The calculated curve is plotted in red along with its confidence/prediction intervals: red shade corresponds to the interquartile range, while the yellow shade corresponds to the 95% confidence interval comprised between the 2.5 and 97.5 percentiles.Reported data for the total number of fatalities are given by the black points.The vertical line, the beginning of the prediction interval, is March 1, 2022.
The identifiability analysis suggests that φ u , φ v and many other parameters are globally identifiable, i.e., they have a unique value given the functions D u (t) and D v (t).However, these theoretical-analysis results do not exactly transfer to numerical calculations for various reasons.A globally identifiable parameter may not necessarily have a sharp estimate due to the potential sloppiness [35] of the model.When the output functions are not sensitive to a particular parameter, a sharp estimate will not be expected, even though the parameter may be globally identifiable.Many other factors, e.g., the reliability and accuracy of the reported time series, may exacerbate the situation.The identifiability analysis assumes that both functions D u (t) and D v (t) are outputs, which includes much more information than the simple loss term Eq. ( 4) can provide.It is an open question how close these estimates are to the actual parameters.Interestingly, however, we find that these estimates still give fairly accurate predictions, even though certain parameter estimates may not be as sharp as desired.When the total death, D u (t) + D v (t), is the only output, we cannot obtain any identifiability results.The implications of there remarks on model identifiability are further elaborated in our comments of best-fit parameters.Hv(0) 370 -

Andalusia
As mentioned, the fitting time window we used to obtain the optimized parameters and initial conditions was from November 15, 2021 to March 1, 2022.The prediction interval ended on March 29, 2022.The time series for Andalusia is available from the Spanish Health Ministry, but we used the series compiled at [36].Note that in Andalusia the reported values of D(t), and the total number of hospitalizations J(t), the latter discussed in Section II B, correspond to the event day, whereas for Switzerland they correspond to the report day.
The vaccination data until April 29, 2021 were also taken from [36].After that date, they were extracted from the Regional Government of Andalusia (Junta de Andalucía, Ref. [37]).We ignored the vaccinations for kids under 12 years old, as there were many data anomalies, resulting in a time series that appears to be problematic.Irrespective of that, this population segment corresponds to only 4% of the total vaccinations.The fatalities time series we used did not report how many fatalities could be attributed to vaccinated or unvaccinated individuals.Therefore, for the region of Andalusia we used only norm (2), coupled to the modified vaccination time series as described in Eq. ( 3), to perform the optimizations.
Figure 2 shows the calculated fatalities time series (both fitting and prediction intervals) and the reported numbers.Table I (model A1 in column 3) presents the optimized parameters and initial conditions.The reported interquantile range arises from 250 fits in the bootstraping step, as discussed above.
We observe that the overall trend of the fatalities seems to be reasonably well captured by the model within its prediction intervals (and their associated uncertainty).We do note, however, a slight over-prediction towards the end of the time series, during March 2022.The transmission rates β ij (i, j = u, v) with at least one member of the unvaccinated population uu and uv are clearly higher than the vv rate between members of the vaccinated population.In fact, β uu is more that three times higher than β vv .The lowest transmission rate is predicted to be β vu , even lower that β vv .At this point it is important to recall our discussion about the use of the total number of deaths and the resulting inability to identify definitively the model parameters.Hence, the above numbers, even when they appear to be intuitively relevant, should be taken with a grain of salt.The latency period is approximately 3.5 days, as imposed by our constraint, and in agreement with [23].The role of asymptomatics, as reflected by the fraction φ i of exposed who become asymptomatics, is considerable, approximately 1/3 and independent of whether the population is vaccinated or not.The calculated fraction of asymptomatics is in reasonable agreement with Ref. [38] who reported a pooled fraction of asymptomatics for the omicron variant of 25.5% (95% confidence interval 17.0% -38.2%).The vaccinated and unvaccinated infectivity period for asymptomatic infections 1/µ i is approximately constant, at about three days, again independent of whether the u or v compartment is considered.
We also find that some parameters are more difficult to justify, Specifically, we find that the recovery rates of symptomatically infected individuals I i → R i and that of the hospitalized individuals H i → R i are almost independent of whether the population is vaccinated or not.The same holds for the transition rates I i → H i and the death rates H i → D i : all four of them are found to have weak variations.The independence of these rates on the administration of the vaccine might be related to the norm we used that does not distinguish between fatalities of vaccinated or unvaccinated individuals.We will return to this point in our analysis of the Switzerland data.

Switzerland
FIG. 3. Omicron-variant model A1: Fit and prediction for the total number of fatalities in Switzerland (norm Eq. ( 2)).The calculated curve is plotted in red along with its confidence/prediction intervals: red shade corresponds to the interquartile range, while the yellow shade corresponds to the 95% confidence interval comprised between the 2.5 and 97.5 percentiles.Reported data for the total number of fatalities are given by the black points.The vertical line, the beginning of the prediction interval, is March 1, 2022.
We chose to perform model calculations for a territory with a population similar in number to that of Andalusia, and for which adequate data are available.As such, we chose a region that contains Switzerland and the Principality of Liechtenstein (data are jointly reported) since the total population of this aggregate territory is 8.7M, (compared to 8.4M for Andalusia).Overall, we followed a procedure very similar to what we used for Andalusia, with a few minor changes.Identical fitting and prediction intervals are used as those for Andalusia.We do note, however, that starting February 17, 2022 most restrictions were lifted in Switzerland.We believe this is one of the reasons we observe a model under-prediction of the number of fatalities in Fig. 3 and in Figs. 4.
Case reporting was slightly different.The Swiss government through the Federal Office of Public Health provides daily the status (vaccinated, unvaccinated or unknown) of each hospitalized/deceased person at Ref. [39].In the absence of a concrete metric on how to partition unknown fatalities to vaccinated and unvaccinated individuals, we used the following procedure to convert these three time series into two, one associated with vaccinated D v (t i ) and the other to unvaccinated individuals D u (t i ).Let dv (t i ), du (t i ) and d (t i ) denote the number of daily reported fatalities with vaccinated, unvaccinated, and unknown state, respectively.We randomly sample an integer number δ i ∈ [0, d (t i )], following a uniform distribution, and then we define the daily number of vaccinated/unvaccinated deceased as ).Note that we followed the same procedure to generate the hospitalizations J v and J u used in model A2, in Section II B 2. As mentioned earlier, D(t) and J(t) for Switzerland correspond to the report day.
Given the reconstructed time series D i (t) we used norm (4), in addition to the norm (2) used in the case of Andalusia, to fit and predict the fatalities time series for the territory of Switzerland (and the Pricipality of Liechtenstein).We attempted to fit separately the vaccinated and unvaccinated deceased, using them as separate inputs to our minimization objective.As mentioned earlier, since we consider that vaccinated individuals have received the full dose (neglecting immunity waning, boosting, of vaccine efficiency) we take θ = 1 in every fit, and V (t) is defined as described in Eq. (5).Our results for the fit to the total number of fatalities are shown in Fig. 3, whereas those for the separate fits to vaccinated and unvaccinated deaths are presented in Fig. 4. Table II, columns two and four, presents the fitting parameters and initial conditions.FIG. 4. Omicron-variant model A1: Separate fits of vaccinated and unvaccinated fatalities in Switzerland (norm Eq. ( 4)).Left panel: Fatalities of vaccinated individuals.Right panel: Fatalities of unvaccinated individuals.Calculated curves are plotted in red along with its confidence/prediction intervals: red shade corresponds to the interquartile range, whereas the yellow shade presents the 95% confidence interval comprised between the 2.5 and 97.5 percentiles.Reported data for the total number of vaccinated and unvaccinated fatalities are given by the black points.The vertical line, the beginning of the prediction interval, is March 1, 2022.
We can see a clear model under-prediction of the fatalities (within the testing period), for both optimizations (norm (2) and ( 4)), despite an accurate following of the time-series trend throughout the period over which regression is performed.The under-prediction is more severe in the case of the total number of fatalities, Fig. 3, and in the vaccinated death time series of the right panel in Fig. 4. As mentioned earlier, we attribute the under-prediction to the fact that after the end of the fitting period, restrictions were considerably relaxed leading to more cases, and eventually more fatalities, a feature that was not explicitly factored in the model.
As regards the parameters of the model, we observe very similar trends to what we obtained for Andalusia.A notable exception is that in the total-deaths fit β uv is the highest transmission rate, retaining however β uu β vv in the case of norm (2).The latency period is well reproduced (as expected due to the constraint and Ref. [23]), and the fraction of asymptomatics is approximately 25% (again in agreement with [38]) irrespective of vaccination or not.The remaining parameters follow similar trends as reported in Table I for Andalusia.It is noteworthy that in both cases recovery, transmission, and death rates seem to depend relatively weakly on whether the vaccine had been administered or not.
A comparison of the parameters predicted by the two optimization is in order.When the two distinct populations are used in the regression, we observe, in the third column in Table II, that the transition rate β vv becomes the largest one with β vu being the smallest.While a calculated higher viral transmissivity of vaccinated individuals could, in principle, be attributed to taking fewer measures to limit pathogen transmission via behavioral changes, e.g., higher contact rates, negligence to use face masks, etc, it is not obvious that such an attribution is meaningful, rather than the potential outcome of the sloppiness of the model.Another surprising feature is that we do not find a significant dependence of the parameters on the norm used (apart from the noted difference in the transmission rates).The asymptomatic fraction is predicted to be slightly larger, approximately 30%, the H u → R u is slightly smaller, and the death rate is slightly larger.TABLE II.Optimal parameters and initial conditions for the omicron-variant models in Switzerland: Model fits to (a) the total number of fatalities (model A1, second column from the left), (b) the total number of hospitalizations (model A2, third column), (c) to the u, v total number of fatalities separately (model A1, fourth column), and (d) to the u, v total number of hospitalizations separately (model A2, fifth column).Population NCHL = 8.7M .Parameters descriptions are as defined in Table I.

Parameter
Median (interquartile range) Median (interquartile range) Median (interquartile range) Median (interquartile range) Fit: total number of deaths Fit: total hospitalizations Fit: u, v deaths separately Fit: u, v hospitalizations separately [Model A1, norm Eq. ( 2)] [Model A2, norm Eq. ( 7) ] [Model A1, norm Eq. ( 4)] [Model A2, norm Eq. ( 8 We now consider model A2, an omicron-variant model similar to A1, but where the population branches terminate at the total number of hospitalizations.The total hospitalization data appear, in our gauge, to be more reliable than fatalities, as the latter (at least in Andalusia) include deceased by any cause that may have recently generated a positive test.That is to say, we believe that numerous fatalities were attributed to COVID even though the primary reason for these events had not been COVID, but another occurrence, see, for example, [40].By considering the reported (total, conventional and critical) hospitalizations, this possible misattribution of fatalities to COVID-19 may be diminished.

Omicron Model A2
Population flows where J(t) in the total number of hospitalizations and γ ij is the rate symptomatically infected individuals u, v become hospitalized γ h,i or recover γ r,i .The optimal parameters (and initial conditions) are obtained by a procedure similar to what we used in model A1 with J(t) playing the role of D(t).Accordingly, the norms change: Eq. ( 2) becomes while Eq. ( 4) becomes 1. Andalusia The results of fitting the total hospitalizations (arising from both vaccinated plus unvaccinated populations using norm ( 7)) are shown in Figure 6.The optimal parameters and initial conditions are summarized in Table I, last column.
As in the case of the cumulative optimization of both the vaccinated and the unvaccinated fatalities, the results of Fig. 6 appear quite accurate, including the forward prediction for the month of March (despite a slight under-prediction of hospitalizations).Note that model A1 predictions were slightly above reported fatalities.However, a more careful inspection of the obtained parameters suggests that some of them may not be epidemiologically realistic.Inspection of the optimized transmission rates shows that the vv rate is larger than the uu rate (β vv > β uu and β vu > β uv and > β uu ).Whereas these inequalities may be related to changes in the behavior of vaccinated individuals, (for example, vaccinated individuals may take fewer precautions and socialize more) we believe instead that this aspect points to the non-identifiability of the model.It is also likely that the origin of these transmission rates stems from our regression's inability to expressly distinguish between the two u, v compartments,.More concretely, similarly to model A1, the appropriate identifiability analysis should consider that the only output is the cumulative number of hospitalizations J u (t)+J v (t).However, for such an output we could not obtain any identifiability results following the procedure described in Appendix A. The model is too complex and beyond the capability of the current identifiability-analysis packages.
FIG. 6. Omicron-variant model A2: Fit and prediction for the total number of hospitalizations in Andalusia (norm of Eq. ( 7)).The calculated curve is plotted in red along with its confidence/prediction intervals: red shade corresponds to the interquartile range, yellow shade to the 95% confidence interval comprised between the 2.5 and 97.5 percentiles.Reported data for the total number of hospitalizations are given by the black points.The vertical line, the beginning of the prediction interval, is March 1, 2022.

Switzerland
We followed the same procedure as that for Andalusia to generate the model A2 fits and predictions for Switzerland (including the Principality of Liechtenstein).As for model A1, we fitted both the total number of hospitalizations via norm (7) and the two vaccination-identified compartments via norm (8).We followed the same procedure as that used to determine D u (t) and D v (t) to obtain estimates for J u (t) and J v (t).
Figure 7 presents our results for the fit to total number of hospitalizations.Figure 8, instead, corresponds to the vaccinated and unvaccinated populations considered separately in the regression.Table II, third and fifth column, summarizes all the fitting parameters and initial conditions.FIG. 7. Omicron-variant model A2: Fit and prediction of the total number of hospitalizations in Switzerland (norm of Eq. ( 2)).The calculated curve is plotted in red along with its confidence/prediction intervals: red shade corresponds to the interquartile range, whereas the yellow shade presents the 95% confidence interval comprised between the 2.5 and 97.5 percentiles.Reported data for the total number of hospitalizations are given by the black points.The vertical line, the beginning of the prediction interval, is March 1, 2022.
We can see a clear model under-prediction of hospitalizations, despite an accurate following of the time-series throughout the period over which regression is performed.The under-prediction is more pronounced in the case of the fit to the total number of hospitalizations.In the case of the separate fittings, the hospitalizations of the vaccinated population are more under-predicted than the hospitalizations of unvaccinated individuals.We attribute this to the fact that, as also discussed in the context of fatalities, towards the end of the fitting period, restrictions were considerably relaxed leading to more cases, and eventually more fatalities.As regards the parameters of the model, we find the transmission rates to be more in line with what one might typically expect.In particular, β uu > β vv for both fittings (with the two different norms), but β uv > β vv for the fitting to the total hospitalization, whereas the reverse is true for the fitting to the two separate populations (although in both cases, the rates are fairly similar, even more so when considering the interquartile ranges).The fraction of asymptomatics varies from approximately 27% to 33%, again in reasonable agreement with [38].The remaining parameters do not seem to significantly depend on the norm chosen.
FIG. 8. Omicron-variant model A2: Separate fits of vaccinated and unvaccinated hospitalizations in Switzerland (norm Eq. ( 4)).Left panel: Hospitalizations of vaccinated individuals.Right panel: Hospitalizations of unvaccinated individuals.Calculated curves are plotted in red along with its confidence/prediction intervals: red shade corresponds to the interquartile range, whereas the yellow shade presents the 95% confidence interval comprised between the 2.5 and 97.5 percentiles.Reported data for the total number of vaccinated and unvaccinated hospitalizations are given by the black points.The vertical line, the beginning of the prediction interval, is March 1, 2022.
In our identifiability analysis of model A2 with the two vaccination compartments treated separately we considered, as in the case of model A1, that J u (t) and J v (t) are separately, and continuously, known.Moreover, we took θ = 1.In the case of model A2 and under the above assumptions, following the procedure described in Appendix A, we can show that all parameters and initial conditions are globally identifiable.This implies that, in principle, all initial conditions and parameters can be determined from the output J u (t) and J v (t).Practically speaking, however, only the discrete time series J u (t) and J v (t) are known: we do not have the full information of the continuous changes of J u (t) and J v (t) as the identifiability analysis supposes.Hence, the loss function used in the parameter estimation is based on a discrete time series reflecting a finite number of observations.Given the complexity of the model and the large number of parameters involved, the optimization package often fails to find a minimum.To alleviate the situation, we chose to fix certain initial conditions, even though such a choice is inconsistent with being globally identifiable.Combined with the possible sloppiness of the model, the result of such a choice may be that the estimates for some of the parameters may not be as sharp.We indicate the above to mitigate a potential impression (to the reader) that the mathematically obtained global identifiability of the model should be expected to translate into the most definitive model results.

III. DELTA VARIANT
Having explored the more elaborate model setting of the omicron variant we now turn to the simpler case of the delta variant.What simplifies the model considerably is that it is sufficient to consider a single susceptibles population since infection of vaccinated individuals was rare.Consequently, susceptibles who are vaccinated are added to a "withdrawn" population.An alternative option is to add them to the recovered population in the sense that this is a terminal compartment of the model.For the time period of 3-4 months for which the delta variant was dominant, the potential waning of immunity (either from recovery or from the vaccine) is not considered sufficient to allow these individuals to replenish the susceptibles compartment.For more persistent variants, replenishing the susceptibles population may be relevant.Furthermore, the SARS-CoV-2 vaccines were highly effective against the delta variant, leading to rather few breakthrough infections.This fact removes the need for a

𝜃𝑉(𝑡) 𝛽
A, P, I A schematic of the population flows (left panel) and the susceptible interactions is shown in Fig. 9.
The initial conditions are taken in a similar fashion as in the omicron variant, except for S(0), which is taken as a fitting parameter.We chose to render it a fitting parameter since susceptibles who became infected with previous variants are immune to the delta variant.However, their number is not definitively known.Moreover, as in our modeling of the omicron variant via models A1 and A2, we supposed that the transmission rate β is the same for all three infectious compartments: asymptomatics, presymptomatics, and for the symptomatically infected population.In addition, as in the case of the omicron-variant models, we imposed the single constraint on the incubation period τ inc = σ −1 1 + σ −1 2 to be equal to a value randomly sampled following a normal distribution whose mean now is 4.41 and standard deviation 0.3291, in line with what is reported in [23].
The norm associated with model B, and minimized during the optimization procedure is: Again, following the approach presented in Appendix A, we find that the three parameters σ 1 , φ, µ and the following combinations are globally identifiable.Thus σ 2 and the sum (γ h + γ r ) are locally identifiable.Initial conditions other than D(0), which is explicitly available through the fatality time series, are not identifiable.As discussed in the identifiability analyses of the two omicron models, this is just the theoretical result, based on the analysis of [32].As we have empirically observed, when we fix certain initial conditions or choose a bound for the parameters to be optimized, the identifiability properties of the model may change.In that light, the relevant parameter identifications should be considered with the associated practical "word of caution" indicated above.

A. Andalusia and Switzerland
In our simulations, the fitting window for the Andalusia calculation started on June 15, 2021, whereas it started on July 1, 2021 for the Switzerland simulations.The choice of the initial time was determined from the existence of a plateau in the associated time series.The fitting period ended on October 1, 2021 for both regions, and the prediction interval terminated on November 1, 2021.As mentioned in Section II the omicron variant appeared in November 2021.
Figure 10 shows the results of our model B simulations, both for the fitting and the prediction intervals.The left panel presents results for Andalusia, whereas the right panel for Switzerland.The best fitting parameters and initial conditions for both countries corresponding to the model B ODEs are presented in Table III.
In the case of Andalusia, we observe a high quality fit, not only for the regression interval but also for the prediction interval.Nevertheless, some parameters do not seem to be in agreement with current knowledge of the epidemiology of the delta variant of SARS-CoV-2.We believe, that the primary reason is the lack of identifiability of the model (both the local aspects thereof theoretically, as well as the practical aspect highlighted above in connection to data and initial condition choices).For example, the recovery time of asymptomatics, 1/µ is found to be ≈ 3.5 days, and an upper bound to the recovery time of symptomatically infected is 1/γ r ≈ 5.4 days.It may be expected that both time scales are likely to be longer than these predictions, although these numbers are in reasonable correspondence with findings, e.g., such as the ones of [41] for the delta variant.On the other hand, the fraction of asymptomatics, 8%, is in agreement with the review and analysis of [38] who found a considerably smaller fraction of asymptomatics associated with the delta than the omicron variants, again in agreement with our calculations.
Our model calculations for Switzerland in Fig. 10 provide a reasonable fit throughout the training interval.: calculations initially under-predict and later over-predict.Nevertheless, the predicted time series considerably under-predicts the number of fatalities over the testing period.Some optimized parameters for this territory differ significantly from those obtained for Andalusia.In particular, the transmission rate for the Switzerland data is higher than that for the Andalusia data, as are the H → R recovery rates.The Switzerland parameters, however, for the death rate and the vaccine efficiency are predicted to be lower than those for Andalusia.In this case, we do not have a definitive attribution of the relevant result (i.e., the underprediction of fatalities) in the case of Switzerland.The only change in policy that we could identify was that from September 13, 2021, access to most indoor public spaces like restaurants, bars, museums or fitness centres was permitted with a valid COVID certificate in Switzerland.No other restrictions were enforced on fully vaccinated and boosted people.

IV. CONCLUSIONS AND FUTURE CHALLENGES
In this work we presented a new class of compartmental epidemiological models that was motivated by the immunological properties of the delta and omicron variants of SARS-CoV-2.More generally, our aim was to present possibilities for settings where variants are highly transmissive (and hence relevant to consider vaccinated individuals and their epidemiological characteristics) as in the case of models A1-A2 for the omicron variant, as well as ones where breakthrough infections are more rare, and hence vaccination is tantamount to withdrawal from the susceptible population as in the case of model B for the delta variant.We, therefore, constructed model B, with the stipulation that vaccinated individuals were permanently withdrawn from the susceptible population based on the vaccination records and vaccine coverage rate.On the other hand, the epidemiology of the omicron variant suggests a substantial number of breakthrough infections, namely infections of vaccinated individuals.Accordingly, we developed models for both vaccinated and unvaccinated populations and analyzed their pairwise interaction and overall time evolution.Indeed, two classes of regression results were given.In the first (and more crude) regression, only the cumulative number of fatalities was accounted for in the optimization objective.This was done when the data did not allow the partitioning of fatalities (or hospitalizations) to vaccinated and unvaccinated components.In the second, more refined approach, fatalities (or hospitalizations) stemming from the two different (vaccinated or not) groups were separately considered.
We addressed the identifiability of the various models and considered mathematical issues (e.g., parameters globally and locally identifiable, given particular time series), we raised some practical considerations due to the finite nature of the available observations, and we considered the compatibility of the selection of some initial conditions.In the case where the time series associated with vaccinated and unvaccinated individuals are required, we identified the issue of how to handle the so-called "unknown" deaths if the individual vaccination status remains undeclared.We proposed a concrete approach to address such disparities, yet clearly these topics merit further investigation.
In our presentation, we focused on the region of Andalusia in Spain and the country of Switzerland (which included data from the Principality of Liechtenstein).These two territories have similar populations.In each territory, we presented studies of a regression effort involving the fatalities (model A1), as well as one terminating at the compartment of total (i.e., conventional plus critical) hospitalizations (model A2).Our models gave generally good agreement with the corresponding training sets, but also reasonable prediction intervals in comparison with the testing data for periods of about a month beyond the end of the training period (up to which the optimization is performed).In the cases where deviations from the predictions were more significant, plausible explanations were offered on the basis of, e.g., the relaxation of measures or other changes of policies.
Naturally, these models offer a starting point for further considerations and are intended as a stepping stone for further studies.On the one hand, it would be quite relevant to seek additional sources of data and other approaches to parameter estimation (than the regression and bootstrapping methodologies used here), to incorporate more accurately the measurement uncertainty and to improve the adequacy of the parameter estimation, in line with our expectations stemming from the analysis of the model identifiability.Another important direction is to add the spatial dimension to the proposed well-mixed ODE models, to incorporate the mobility of vaccinated individuals.This can be done either at the level of metapopulation models [8,42,43] or at that of PDE approaches [44][45][46][47].Finally, numerous additional dimensions of such modeling of vaccinations are relevant to consider such as, e.g., the age stratification of such effects [2,20,48].These directions are currently under consideration and will be reported in future publications.

V. DECLARATION OF COMPETING INTEREST
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
So far, the state variables for unvaccinated population H, I, P, E, S are expressed as functions of D and its derivatives.We have identical formulas for the vaccinated populations H v , I v , P v , E v , S v , except that for S v a negative sign should be added in front of θ.
The equation for E u , Eq. (1b), multiplied by (1 − φ u )σ 1 σ 2 γ h,u κ d,u becomes Note that the subscript u for the parameters k 3 , k 2 , k 1 , k 0 is still omitted.But for this equation they should be computed from the parameters associated with the unvaccinated population, whereas for the equation of D The next step is to scale variables as follows The above equation for D (5) u becomes The equation for A u , Eq. (1d), takes the form 4)  u + D u (κ s,u + γ s,u + σ 2 ) + D u (κ s,u γ s,u + κ s,u σ 2 + γ s,u σ 2 ) + D u (κ s,u γ s,u σ 2 ) − µ u A u .
Similarly, the equation of A v and the high order equation for D v are: Up to now, we rewrote the original system as a system of D (5) u , D v , A u , A v .It can be written as a first-order system (by using D, D , D , D (3) , D (4) ) so that the identifiability analysis package StructuralIdentifiability can be applied to find the identifiability property of the 18 parameters of this new system: β uu , β uv , β vu , β vv , σ 1 , σ 2 , φ u , φ v , γ s,u , γ s,v , µ u , µ v , κ s,u , κ s,v , θu , θv , α u , α v Then, the identifiability property of the 19 parameters of the original system β uu , β uv , β vu , β vv , σ 1 , σ 2 , θ, φ u , φ v , µ u , µ v , γ r,u , γ h,u , γ r,v , γ h,v , κ r,u , κ d,u , κ r,v , κ d,v can be derived.Furthermore, one may use the SIAN Webapp, with the globally identifiable parameters (from StructuralIdentifiability) as extra outputs, to find the identifiability property of the initial conditions.Identifiability results obtained from these calculations are reported in appropriate sections in the main text.

FIG. 1 .
FIG. 1. Schematic diagram of population flows according to model A1 (left panel) and susceptible interactions with other population compartments, for both models A1 and A2 (right panel).The symbol δij with i, j = u, v is the Kronecker delta.

FIG. 5 .
FIG. 5. Schematic diagram of the population flows according to model A2.The susceptible interactions are as in model A1, shown in the the right panel of Fig. 1.The symbol δij with i, j = u, v is the Kronecker delta.

FIG. 9 .σ 2 P
FIG. 9. Schematic diagram of population flows according to the delta-variant model B (left panel) and susceptible interactions with other population compartments (right panel).A single population is modeled, as we consider that the waning immunity time scale (either due to vaccine immunity or due to recovery) is much longer than the time scale of prevalence of the delta variant.Vaccinated individuals W are permanently removed from the susceptible compartment, their population becoming a terminal compartment of the model.

FIG. 10 .
FIG.10.Delta-variant model B: Fit and prediction of the total number of fatalities in Andalusia (left panel) and Switzerland (right panel).Norm(10) was used.The vertical line, the beginning of the prediction interval, is October 1, 2022.Note the significant difference in the number of fatalities.

TABLE I .
Optimal parameters and initial conditions for the omicron-variant models in Andalusia: Model fits to the total number of fatalities, discussed in Section II A 1 (model A1, third column) and to the total number of hospitalizations, discussed in Section II B 1 (model A2, right column).Population NAnd = 8.4M .

TABLE III .
(10)mal parameters and initial conditions for the delta-variant model, model B, in Andalusia (third column) and Switzerland (fourth column).Model fits to the total number of fatalities, norm(10).