Model-based estimation of transmissibility and reinfection of SARS-CoV-2 P.1 variant.

Background: The SARS-CoV-2 variant of concern (VOC) P.1 (Gamma variant) emerged in the Amazonas State, Brazil, in November 2020. The epidemiological consequences of its mutations have not been widely studied, despite detection of P.1 in 36 countries, with local transmission in at least 5 countries. A range of mutations are seen in P.1, ten of them in the spike protein. It shares mutations with VOCs previously detected in the United Kingdom (B.1.1.7, Alpha variant) and South Africa (B.1.351, Beta variant).
Methods: We estimated the transmissibility and reinfection of P.1 using a model-based approach, fitting data from the national health surveillance of hospitalized individuals and frequency of the P.1 variant in Manaus from December-2020 to February-2021.
Results: Here we estimate that the new variant is about 2.6 times more transmissible (95% Confidence Interval: 2.4-2.8) than previous circulating variant(s). Manaus already had a high prevalence of individuals previously affected by the SARS-CoV-2 virus and our fitted model attributed 28% of Manaus cases in the period to reinfections by P.1, confirming the importance of reinfection by this variant. This value is in line with estimates from blood donors samples in Manaus city. Conclusions: Our estimates rank P.1 as one of the most transmissible among the SARS-CoV-2 VOCs currently identified, and potentially as transmissible as the posteriorly detected VOC B.1.617.2 (Delta variant), posing a serious threat and requiring measures to control its global spread.

Introduction

The Japanese National Institute of Infectious Diseases identified the new P.1 SARS-CoV-2 variant from travelers returning from Amazonas State, Brazil, on 6-January-2021[1]. P.1 was eventually reported in Manaus city (Amazonas state capital), on 11 January-2021[2]. Later, it was identified in samples collected since 6-Dec-2020 from Manaus[3]. According to phylogenetic studies, P.1 likely emerged in the Amazonas state in early[3] or late[4] November 2020. This variant shares mutations with other variants of concern (VOCs) previously detected in the United Kingdom and South Africa (B.1.1.7 and B.1.351, respectively)[2]. Mutations of these two other variants are associated with greater transmissibility and immune evasion[5,6], which confer them the status of variant of concern. However, information, data, and analyses on the epidemiology of P.1 are still incipient.

The Coronavirus disease 2019 (COVID-19) outbreak in Manaus (April–May 2020) was followed by a period of high but stable incidence, after which the proportion of individuals who were infected by the SARS-CoV2 virus may have reached 42%[7] to 76%[8] by November 2020. From December 2020 to February 2021 the city was devastated by a new outbreak that caused a collapse in the already fragile health system[9], with shortages of oxygen supply[10], while the frequency of P.1 increased sharply from 0% in November 2020 to 73% in January 2021[4]. The pathogenicity of P.1 variant is still unknown, although recent studies point to increased viral load in individuals infected with the new variant[4], suggesting it could be higher than the one from previous circulating strain. We analyzed Brazilian national health surveillance data on COVID-19 hospitalizations and the frequency of P.1 among sequences from residents of Manaus city using a model-based approach (an extended SEIR compartmental model—see Fig. 1) to estimate the relative transmissibility in comparison to the previous local variant(s), and relative force of reinfection of the P.1 variant, i.e. the ratio between the force of infection by P.1 on previously infected individuals (reinfections) and the force of infection by P.1 on susceptible ones (new infections). We estimate that the P.1 variant is about 2.6 times more transmissible than previous circulating variant(s), and 28% of Manaus cases in the period were due to reinfections.

Fig. 1

Diagram of the extended deterministic compartmental model (SEAIHRD).

The model compartments and the respective connections between them are summarized in this diagram, and they are named as S: Susceptible, E: Exposed (pre-symptomatic), H: Hospitalized (severe infected individuals), I: Infected (symptomatic individuals, not hospitalized), A: Asymptomatic. D: Deceased, R: Recovered. Compartments are subdivided into three age categories, not represented here for simplicity. Compartments with subindex 1 represent the wild-type variant, subindex 2 refers to the VOC P.1. Continuous lines represent flux between each compartment; dashed lines, infection pathways. Small arrows indicate force of reinfection and transmissibility. λ = force of infection. β = relative transmission rate. p = relative force of reinfection. γ = average time between being infectious and presenting symptoms. σ = proportion of severe cases that require hospitalization. α = proportion of asymptomatic cases. ν = average time between being infectious and recovering for severe cases. ν = average time between being infectious and recovering for mild/asymptomatic cases. μ = in-hospital mortality ratio.

Methods

In order to estimate key parameters of the variant of concern (VOC) P.1, we developed a model and fitted it to time-series data of number of hospitalized cases and frequency of the P.1 variation. The fitting approach used here can be applied to other regions where as soon as data on number of cases and frequency of a new variant are available. It primarily requires surveillance data to determine proper model initial conditions. In Brazil, these are the hospitalized cases data. Stratification by age allows the model to also consider the different death rates, asymptomatic and hospitalized proportions of each age class, important features for SARS-CoV-2. Contact levels between different age classes, which may vary from one place to another, can also be considered. For special cases in which information such as contact between ages classes and age distribution are not available (or even unnecessary for some other disease), the model can be easily simplified. In this sense, the method proposed here demands low-detailed data and relies on the structure of a simple compartmental model to measure quantities of interest, such as transmissibility and relative force of reinfection. More information about the methodology here applied is available in the Supplementary Information (SI) in sections 1–3.

Model

A deterministic compartmental model (Fig. 1) was developed to model the infection of Susceptible individuals moving to the Exposed (pre-symptomatic) compartment, which can progress to three alternative compartments: Hospitalized (severely ill), Infected (symptomatic but non-hospitalized), and Asymptomatic. Eventually, individuals move to Recovered or Deceased. Two variants are considered: 1-non-P.1 (“wild-type”) and 2-new/P.1. The latter is assumed to infect Recovered individuals previously infected by the wild-type, and no reinfections of wild-type due to waning immunity occur. Compartments were stratified into three age categories: young (<20 years old), adult (≥20 and <60 years old) and elderly (≥60 years old), with different rates for outcomes. The key parameters of relative transmissibility and relative force of reinfection were estimated by a maximum likelihood fitting to the weekly number of new hospitalizations and to genomic surveillance data. Three additional model parameters with unknown values were also estimated. The remaining parameters (24 out of 29) were fixed, using current values from the literature. Sensitivity to different pathogenicity of the P.1 variant was explored by repeating the fit assuming IHR as a free parameter (SA1). The sensitivity to the period analysed was also explored by another fit excluding the health system collapse period (SA2) (see below. Further model details about the model are available in the Supplementary Methods (Section 1 of the SI).

Dataset

We used the Brazilian epidemiological syndrome surveillance system for influenza, SIVEP-Gripe (https://opendatasus.saude.gov.br), to track COVID-19 hospitalized cases. All hospitalized patients with Severe Acute Respiratory Illness (SARI) are reported to SIVEP-Gripe with symptom onset date and SARS-CoV-2 test results. SIVEP-Gripe, due to its universal coverage and mandatory notification of SARI cases, has an homogeneous testing effort to diagnose SARS-CoV-2 infections, and is currently the best source for Brazilian data at the national level. As the data are publicly available by the Brazilian National Health System, no ethical approval was needed to perform the analysis, according the the National Ethical Commission (CONEP) of the National Health Council, Resolution Number 510 of April 7, 2016. Hospitalization data provides the most accurate basis to infer incidence in Manaus, because mild cases are vastly under-reported and testing capacity fluctuates, while mortality data are harder to relate to total number of cases, since the city’s health system endured a prolonged stress even before the collapse, with large variation in the in-hospital fatality rate over time[11]. Time-series of frequency of sequenced genomes identified as P.1 in Manaus were extracted from published datasets[3,12]. Because of the uncertainty of the mortality rate, hospitalization data are more reliable to monitor the Brazilian situation. However, hospitalization suffers an important dampening when the health system collapses, presenting a false spreading control precisely because of the overload on the health services, which cannot admit more patients than their capacity. We take this into account making a sensitivity analysis (SA2—see below). Likewise, when the health system is overwhelmed, the IHR can be higher (as it can be only a trait of the P.1 variant), leading to greater mortality rates than the parameters used in this work. For this reason, we also taken the differing IHR into account for P.1 using a sensitivity analysis (SA1—see below). Data and code are available in the following link: https://zenodo.org/record/5594600[13], and in the Github repository https://github.com/covid19br/model-P1-variant.

Nowcasting

Data for hospitalized COVID-19 cases among residents in Manaus from 01 November 2020 to 31 January 2021 was obtained from SIVEP-Gripe database of 15 February 2021. The hospitalized cases of the last 10 weeks in the time-series were nowcasted[14] to correct for notification delay Data used in parameters estimation were collected from the SIVEP-Gripe In this system, reporting of cases can be delayed for several reasons, including the notification system itself and confirmation of RT-PCR test results. The nowcasting procedure estimates, based on the past delay distribution, the number of cases that already occurred but were not yet reported. A window of 10 weeks is the acting window on the series, since delays greater than this are rare.

Nowcasting requires a pair of dates: (i) onset date of the event and (ii) report date of the event. The delay distribution is modeled as being best described as a Poisson distribution for days since the onset date to the report date. We considered the first symptoms date as the onset date. For the report date, we used the latest between the test result date and the clinical classification date. The nowcasting algorithm were developed by ref. [14], and implemented in the NobBS (Nowcasting by Bayesian Smoothing) package in [15].

Model parametrization and initial condition estimation

The model requires appropriate mid-epidemic initial conditions in order to give relevant results. In the model, the number of new hospitalizations at a given time—h, is directly proportional to the number of exposed individuals at that time, therefore data was used to get an approximation of the number of exposed people. Also, to quantify the number of people belonging to the recovered class, prevalence was used. In Table 1 we present the parameters considered for the wild-variant are described below. The parameters for the P.1 variant are the same except for those considered in the model fitting and in the Sensitivity Analysis 1 (SA1). See more detailed for the analysis of the initial conditions in the Supplementary Methods (Section 3 of the SI).

Table 1

Epidemiological parameters.

Parameter	Description	Value	Source
γ	Average time in days between being infected and developing symptoms	5.8	[25]
νi	Average time in days between being infectious and recovering for asymptomatic and mild cases	9.0	[26]
νs	Average time between being infectious and recovering/dying for severe cases	8.4	SIVEP-Gripe for São Paulo State
ξ	Reduction on the exposure of symptomatic cases (due to symptoms/quarantining)	0.1	Assumed
ξsev	Reduction on the exposure of severe cases (due to hospitalization)	0.9	Assumed
ω	Relative infectiousness of pre-symptomatic individuals	1.0	Assumed
α	Proportion of asymptomatic cases	[0.67, 0.44, 0.31]	Juvenile[27] Adult and Elderly[28]
σ	Proportion of symptomatic cases that require hospitalization	[0.001, 0.012, 0.089]a	[29]
μ	In-hospital mortality ratio	[0.417, 0.188, 0.754]	[11]
χ	Case report probability	1.0	Assumed

aThe proportion is weighted by the age distribution of the population with each age category.

Maximum likelihood estimation

Given the cumulative daily curves of hospitalization for wild-type variant (C1), and P.1 variant, (C2) we can obtain the daily variation of each curve (namely and ). Those curves are summed up to give the total number of weekly new cases:where τ is a discrete index given in weeks.

To calculate the frequency of P.1 in a given time period T, we use the proportion of new cases in this period from the wild-type and P.1 variant as follows:where is a discrete index given in T periods. The time period T depends on the dataset of genome sequences (weekly[3] and monthly[12]).

Using the maximum likelihood method, we fitted the model by estimating five parameters, namely, the relative transmissibility (), the relative force of reinfection of P.1 (p), initial total prevalence (ρ0 = [R/N]), initial fraction of cases that were caused by the new variant (P0), and intrinsic growth rate of the wild-type variant (r). The initial fraction of P.1 cases (P0) accounts for the uncertainty in the time of emergence of the new variant: the simulation starts at beginning of November, so this initial value is below 1 individual, and only reaches this threshold by mid to late November, depending on the value of P0. The parameter r incorporates effects related to contact rates for the wild-type variant, such as non-pharmacological interventions relaxation, elections, and others; it affects the transmissibilities of both variants (β1 and β2) in the same way, and so is independent of Δβ.

Number of hospitalization cases were assumed to follow a Poisson distribution, with expected value given by Eq. (1). The recorded number of P.1 in genome samples was assumed to follow a binomial distribution with an expected value equal to the product of the total number of genome sequences sampled in each date and the proportion of P.1 cases (Eq. (2)). The log-likelihood function for the model fitting was then:where Pois is a Poisson distribution with parameter λ, x is the number of recorded hospitalizations in week i, Bin is a Binomial distribution with parameters N (total number of trials) and π (probability of success at each trial), n is total number of sequences in clinical samples in week or day j, y is the number of P.1 sequences in each of these samples, and θ(. ) is the logit function.

The model was then fitted by finding the values of the five above mentioned parameters that minimize the negative of the log-likelihood function (Eq. (3)), using the function mle2, from the R package bbmle[16]. To find starting values for the optimization performed by mle2 we calculated the log-likelihood function for one million combinations of parameters values in a regular reticulate within reasonable ranges. The 100 sets of parameters that were local minima (that is, with highest log-likelihood values) were used as starting values for the computational minimization.

The confidence intervals for the expected number of cases and frequency were estimated from 20,000 parametric bootstrap samples assuming that the estimated parameters follow a multivariate normal distribution. The parameters of these multivariate distributions were the estimated values and estimated variance-covariance matrix of the parameters. We then calculate the 2.5% and 97.5% quantiles of each parameter to obtain the confidence intervals of our estimates. We check the reliability of these estimates by verifying that the log-likelihood profiles satisfy required conditions for identifiability, as detailed in the Supplementary Methods (Section 4 of the SI and Fig. S1).

Sensitivity analysis

The model fitting assumed a constant infection hospitalization rate (IHR, parameter σ) for each age group over time for both variants. Recent evidence suggests that prior SARS-CoV-2 infection protects most individuals against reinfection[17], so reinfections might have lower IHR. Because the pathogenicity of the P.1 variant is unknown, the model fitting was repeated assuming that the odds ratio of the IHR in each age class for P.1 infections compared to wild-type variant infections (SA1) is a free parameter. Moreover, as the collapse of Manaus health system hindered hospitalizations of new severe cases and may have affected case recording in surveillance databases, the model fitting was repeated considering only the period prior to the collapse (10 January 2021) (SA2). Sensitivity analysis, latin hypercube explorations and likelihood profiles characterization are important methods which can be applied in this kind of model, where fitting depends on a set parameters, which were obtained from the literature and one has no further information about confidence intervals.

Table 2

Summary of the fitted parameters and respective confidence intervals considering the entire period, November 1, 2020 to January 31, 2021 maintaining the same pathogenicity of the previous variant.

Parameter	Main fitting			SA1			SA2
	Estimate	2.5%	97.5%	estimate	2.5%	97.5%	Estimate	2.5%	97.5%
Relative transmission rate for the new variant	2.61	2.45	2.76	2.52	2.28	2.76	2.95	2.70	3.20
Relative force of reinfection of P.1	0.032	0.026	0.040	0.053	0.044	0.065	0.000	0.000	0.000
Prevalence of previous infection (2020-11-01) (%)	78	73	83	73	67	78	71	69	74
Initial fraction of the new variant (2020-11-01) (×10−5)	30.4	8.2	112.9	8.5	1.4	50.8	17.6	5.0	62.4
Intrinsic growth rate (days−1)	0.029	0.024	0.035	0.045	0.037	0.052	0.030	0.026	0.034
Relative IHR odds ratio	1a	–	–	0.74	0.63	0.85	1a	–	–

Sensitivity analyses were performed considering different pathogenicity of the P.1 variant (SA1) and data censuring after the collapse of the healthcare system (SA2) in Manaus, Brazil, on January 10, 2021.

aParameter was fixed, not estimated, in this analysis.