Quentin Clairon

Mélanie Prague

Delphine Planas

Timothée Bruel

Laurent Hocqueloux

Thierry Prazuck

Olivier Schwartz

Rodolphe Thiébaut
email: rodolphe.thiebaut@u-bordeaux.fr

Jérémie Guedj

Modeling the evolution of the neutralizing antibody response against SARS-CoV-2 variants after several administrations of Bnt162b2

Because SARS-CoV-2 constantly mutates to escape from the immune response, there is a reduction of neutralizing capacity of antibodies initially targeting the historical strain against emerging Variants of Concerns (VoC)s. That is why the measure of the protection conferred by vaccination cannot solely rely on the antibody levels, but also requires to measure their neutralization capacity. Here we used a mathematical model to follow the humoral response in 26 individuals that received up to three vaccination doses of Bnt162b2 vaccine, and for whom both anti-S IgG and neutralisation capacity was measured longitudinally against all main VoCs. Our model could identify two independent mechanisms that led to a marked increase in humoral response over the successive vaccination doses. In addition to the already known increase in IgG levels after each dose, we identified that the neutralization capacity was significantly increased after the third vaccine administration against all VoCs, despite large inter-individual variability. Consequently, the model projects that the mean duration of detectable neutralizing capacity against non-Omicron VoC is between 366 days (Beta variant, 95% Prediction Intervals PI [323; 366]) and 606 days (Alpha variant, 95% PI [555; 638]). Despite a very low protection after three doses, the mean duration of detectable neutralizing capacity against Omicron variants varies between 184 days (BA.5 variant, 95% PI [155; 215]) and 268 days (BA.1 variant, 95% PI [238; 302]). Our model shows the benefit of incorporating the neutralization capacity in the follow-up of patients to better inform on their level of protection against the different SARS-CoV-2 variants as well as their optimal timing of vaccine administration. December 14, 2022 1/19

Developed vaccines against SARS-CoV-2 have been a turning point against the ongoing Covid-19 pandemic. When the Wuhan virus was dominant, they help to dramatically reduce the number of severe cases as well as infection and transmission rates. For mRNA vaccines, it was in great part explained by the high level of induced antibodies a few weeks/months after injection and linked to high neutralizing capacity, the ability to prevent viruses to enter and infect target cells. However, decreasing antibody concentration over time and apparition of variants escaping their neutralizing action dramatically reduced the initial vaccine efficacy. As a countermeasure, additional injections were used to re-establish significant antibody population and ensure a long-term neutralizing activity against emerging variants. To infer if this multi-dose strategy fulfils such task, we construct a model of the evolution of the induced antibodies and their neutralizing capacity against different variants. This model helps us to quantify the gain brough by each new injection on both antibody population and their neutralizing ability against all tested variants as well as the dramatic differences between them. We also predict the long-term evolution of neutralizing activity, years after last injection, and thus discuss the longevity of the induced protection by vaccine.

Introduction

The discovery and the rapid availability of several vaccines against SARS-CoV-2 virus has been a turning point in the combat against Covid-19 [START_REF] Watson | Global impact of the first year of COVID-19 vaccination: a mathematical modelling study[END_REF]. Although their efficacy may vary to some extent, it is undisputable that large scale vaccination campaigns have dramatically reduced both the risk of severe diseases [START_REF] Thomas | Safety and efficacy of the BNT162b2 mRNA Covid-19 vaccine through 6 months[END_REF][3][START_REF] Baden | Efficacy and safety of the mRNA-1273 SARS-CoV-2 vaccine[END_REF] and, to a lesser extent, the rates of transmission and disease acquisition [START_REF] Harris | Effect of vaccination on household transmission of SARS-CoV-2 in England[END_REF][START_REF] Pritchard | Impact of vaccination on new SARS-CoV-2 infections in the United Kingdom[END_REF][START_REF] Eyre | Effect of Covid-19 vaccination on transmission of alpha and delta variants[END_REF], resulting in millions of saved lives [START_REF] Watson | Global impact of the first year of COVID-19 vaccination: a mathematical modelling study[END_REF][START_REF] Gupta | Vaccinations Against COVID-19 May Have Averted Up To 140,000 Deaths In The United States: Study examines role of COVID-19 vaccines and deaths averted in the United States[END_REF][START_REF] Meslé | Estimated number of deaths directly averted in people 60 years and older as a result of COVID-19 vaccination in the WHO European Region, December 2020 to[END_REF].

However vaccine efficacy has been jeopardized by the apparition of various Variants of Concern (VoCs) that partially escape immune protection. A clear decrease in the neutralization capacity has been observed [START_REF] Cele | Omicron extensively but incompletely escapes Pfizer BNT162b2 neutralization[END_REF][START_REF] Planas | Reduced sensitivity of SARS-CoV-2 variant Delta to antibody neutralization[END_REF] which has translated to a substantial reduction of efficacy against transmission and disease acquisition with Delta and Omicron variants, and, to a lesser extent, to a decrease of efficacy against severe COVID19 disease [START_REF] Bernal | Effectiveness of Covid-19 vaccines against the B. 1.617. 2 (Delta) variant[END_REF][START_REF] Andrews | Covid-19 vaccine effectiveness against the Omicron (B. 1.1. 529) variant[END_REF]. The concern caused by a potential loss of protection against VoC has been further enhanced by the natural waning immunity and the progressive reduction in antibody levels over time [START_REF] Mizrahi | Correlation of SARS-CoV-2-breakthrough infections to time-from-vaccine[END_REF][START_REF] Goldberg | Waning immunity after the BNT162b2 vaccine in Israel[END_REF][START_REF] Levin | Waning immune humoral response to BNT162b2 Covid-19 vaccine over 6 months[END_REF]. This has supported boosting strategies with one or two additional vaccine doses to maintain a high level of protection. However the optimal time to administer boosters, and how these times may vary for different VoC, remains unclear.

To characterize in detail the duration of protection against SARS-CoV-2, it is therefore essential to measure not only total anti-S IgG antibodies over time, as typically done in large observational studies, but also how this translates in terms of neutralization capacity. The latter requires intensive in vitro measurements, but it provides a much more accurate description of the level of protection present in the sera of COVID-19 vaccine recipients [START_REF] Cromer | Prospects for durable immune control of SARS-CoV-2 and prevention of reinfection[END_REF][START_REF] Zhu | WHO international standard for SARS-CoV-2 antibodies to determine markers of protection[END_REF]. Then, a detailed characterization of the immunological or virological factors modulating the duration of protection can be obtained by using mathematical models of immune marker dynamics [START_REF] Korosec | Long-term predictions of humoral immunity after two doses of BNT162b2 and mRNA-1273 vaccines based on dosage, age and sex[END_REF].

Here we propose to use for the first time a mathematical model to analyze the joint kinetics of anti-S IgG antibodies and neutralization capacity after repeated vaccine injections against the main VoCs. For that purpose we relied on data from a cohort of Bnt162b2 vaccine recipients, in which both antibody kinetics and neutralizing activity were measured longitudinally [START_REF] Planas | Reduced sensitivity of SARS-CoV-2 variant Delta to antibody neutralization[END_REF][START_REF] Planas | Sensitivity of infectious SARS-CoV-2 B. 1.1. 7 and B. 1.351 variants to neutralizing antibodies[END_REF][START_REF] Planas | Duration of BA. 5 neutralization in sera and nasal swabs from SARS-CoV-2 vaccinated individuals, with or without Omicron breakthrough infection[END_REF]. We built on previous models of antibody kinetics [START_REF] Pasin | Dynamic of the humoral immune response to a prime-boost Ebola vaccine: quantification and source of variation[END_REF][START_REF] Balelli | A model for establishment, maintenance and reactivation of the immune response after two-dose vaccination regimens against Ebola virus[END_REF] to develop a novel approach to quantify the kinetics of neutralization December 14, 2022 2/19 capacity, and we use this model to characterize the effects of repeated vaccine administrations on the neutralizing activity. We finally use the model to discuss the duration of protection conferred by Bnt162b2 against VoCs.

Materials and methods

Data

Population Study

Data originate from a cohort of N=29 subjects who received up to three injections of Bnt162b2 (ClinicalTrials.gov:NCT04750720 and ClinicalTrials.gov:NCT05315583). In brief SARS-CoV-2 naive patients were recruited in Orléans, France between August 27, 2020 and May 24, 2022. Individuals were followed for up to 483 days after their first vaccine injections (see more details on the data in [START_REF] Planas | Reduced sensitivity of SARS-CoV-2 variant Delta to antibody neutralization[END_REF][START_REF] Planas | Sensitivity of infectious SARS-CoV-2 B. 1.1. 7 and B. 1.351 variants to neutralizing antibodies[END_REF][START_REF] Planas | Duration of BA. 5 neutralization in sera and nasal swabs from SARS-CoV-2 vaccinated individuals, with or without Omicron breakthrough infection[END_REF]). Two patients without longitudinal follow-up and 1 immunocompromised individual were not included in our analysis. In total, N = 26 individuals were analyzed (see

Longitudinal markers of immune response

Two types of measurements were available at each visit: 1) anti-spike binding IgGs, measured in BAU/mL) neutralization titers of sera provided in ED 50 , which is the effective dilution required to neutralize 50% of an arbitrary viral load of reference (eg, the higher the ED 50 the larger the protection level). Neutralization capacity was assessed against historical strain (D614G), Alpha, Beta, Delta, and Omicron variants (strains BA.1, BA.2, and BA.5).

In brief IgGs markedly increased after each dose, but rapidly declined over time, with a rate that did not substantially differ after the second or the third dose (Fig In contrast the kinetics of neutralizing activity was much more heterogeneous, and was characterized by large differences against the different VoCs. Further the neutralizing activity was markedly increased after the third dose against all Omicron variants, albeit remaining at much lower levels than against the other VoC (Fig 1 -bottom).

Mechanistic model of B cells response

Mechanistic model for antibody kinetics

We rely on a simplified and rescaled version of a previously published model in the context of vaccine against Ebola infection [START_REF] Balelli | A model for establishment, maintenance and reactivation of the immune response after two-dose vaccination regimens against Ebola virus[END_REF]. In brief, after each dose, cells transfected with Bnt162b2 generate antigen, noted V , which triggers the constitution of a memory December 14, 2022 4/19 compartment, noted M , at a rate ρ. This memory compartment is a general one accounting for all cell populations able to differentiate into secreting cells upon antigen presence. These can be activated or memory B-cells either circulating or present in germinal centers. So, M can differentiate into secreting plasma cells, noted S, at a rate µV . These cells then produce antibodies, noted Ab, at a rate θ. V , S, and Ab are degraded at rates δ V , δ S , and δ Ab , respectively, leading to the following ODE system:

V = -δ V V Ṁ = ρV -µV M ˙ S = µV M -δ S S Ȧb = θ S -δ Ab Ab. (1)
Assuming that individuals are naive of infection, and noting t 1 the time of first injection, the initial conditions are given by: M (t 1 ) = S(t 1 ) = Ab(t 1 ) = 0.

To model the effect of repeated doses, we consider that V is a function presenting discontinuity at time of first, second and third injection (t

Ṡ = f M k e -δ V (t-t k ) -δ S S Ȧb = ϑS -δ Ab Ab (2)
where f M k = M k M1 is the fold-change for steady-state memory compartment after k th injection compared to the first one (by definition f M1 = 1). Of note, [START_REF] Balelli | A model for establishment, maintenance and reactivation of the immune response after two-dose vaccination regimens against Ebola virus[END_REF] initial model contains two populations of secreting cells S and L, differing by their life expectancy. In our case, preliminary statistical analysis conclude that there was no statistical differences between model adjustements when accounting for S and L or S only (results no shown). This allows us to reduce the number of unknown parameters. This is crucial for parameters related to cell kinetics known to be very different for newly developed mRNA vaccines comparing to viral vector ones and for which no values have been previously inferred. Thus, the retained model ( 2) is complex enough to account for the effect of multiple injections on antibody concentration evolution while avoiding identifiability issues. We define η ODE = (f M2 , f M3 , δ V , δ S , ϑ, δ Ab ) the vector of model parameters defining the dynamics of the system.

Functional model for neutralizing activity

After modeling antibody concentration evolution in the previous section, we aim to model their neutralizing activity. This means in our case proposing a model describing the evolution of ED ν 50 with respect to Ab. We consider the following linear model:

ED ν 50 (t) = F (ν, t)Ab(t).
The function F (ν, t) represents the relationship between the concentration of binding antibodies in BAU/mL and its neutralization capacity against the VoC ν. It is variant-specific and time-varying, let us first derive its expression for the strain D614G before moving to any arbitrary VoCs. After t 1 , we assume a proportional relationship December 14, 2022 5/19

between Ab and neutralizing activity against D614G i.e F (D614G, t) = γ (equivalently ED D614G 50 (t) = γAb(t)). After additional injections, we assume there is a neutralization gain quantified by the fold-change f 2 after t 2 and f 3 after t 3 i.e. F (D614G, t) = γf 2 when t ∈ [t 2 ; t 3 ] and F (D614G, t) = γf 3 for t ≥ t 3 . Now, we account for VoCs specific neutralizing activity by modifying baseline value γ by the fold-changes f ν such that

F (ν, t) = F (D614G, t)f ν = γf ν when t ∈ [t 1 ; t 2 ] and F (ν, t) = F (D614G, t)f ν = γf ν f 2 for t ∈ [t 2 ; t 3 ].
We assume that the relative gain brought by third injection can be also VoC-specific. That is why we introduce the fold-changes g ν to quantify this gain i.e. F (ν, t) = F (D614G, t)f ν g ν = γf ν f 3 g ν for t ≥ t 3 . This piece-wise constant function can be then expressed in a general form:

F (ν, t) = γf ν (1 t<t2 + f 2 1 t∈[t2;t3] + f 3 g ν 1 t≥t3 ).
The choice of this model is the result of exploration based on the minimization of an adjustment criteria. In particular, the current model only quantifies the effect of the repetition of injections on affinity enhancement. Other factors can play a role as the elapsed time since antigen presentation, for example to account for the progressive Memory B-cells repertoire expansion [START_REF] Muecksch | Increased potency and breadth of SARS-CoV-2 neutralizing antibodies after a third mRNA vaccine dose[END_REF][START_REF] Wang | Memory B cell repertoire from triple vaccinees against diverse SARS-CoV-2 variants[END_REF]. An alternative neutralization model only considering the time factor has been developed. This supplemental analysis is described in Appendix B in S1 Appendix but lead to a less accurate model (in terms of AIC). A general model accounting for both factors, the number of injections and the elapsed time, has been also tested leading to non-significant improvements over the retained model and at the expense of identifiability problems (results not shown). We also investigate the possibility of a variant-specific fold-change after second injection. This was discarded due to practical identifiability issues.

We define η = (η ODE , γ, f ν , f 2 , f 3 , g ν ) the vector of model parameters that have to be estimated from the observed data. Description of the model parameters can be found in Table 2 Observation model

The structural model used to describe the log-transformed concentration of binding antibodies in BAU/mL for the i th individual (i = 1, . . . , N ) at the j th time point (j = 1, . . . , n i ) is:

Y BAU ij = log 10 (Ab(η i , t ij )) + e BAU ij ,
where e BAU ij is the residual additive error which follow a normal distribution of mean zero and constant standard deviation σ BAU . The vector η i is the specific value for individual i of vector η.

We also consider a log-transformation of ED ν 50 raw measurements for the variant in the list {D614G, Alpha, Beta, Delta, BA.1, BA.2, BA.5}. For the i th individual at the j th time point, we have:

Y ν ij = log 10 (ED ν 50 (η i , t ij )) + e ν ij ,
where e ν ij is the residual additive error for variant ν which follow a normal distribution of mean zero and constant standard deviation σ ν .

Statistical model for parameters over time and injections

Fixed parameters. Because not all parameters can be estimated when only concentration of binding antibodies and antibody neutralizing activity are measured, some parameters had to be fixed based on the current literature. Here, we just fix one parameter, the clearance of antibodies found to be stable across studies and infections [START_REF] Fatehi | Comparing antiviral strategies against COVID-19 via multiscale within-host modelling[END_REF] with a half life of 21 days. Further, the model predictions were found largely insensitive to the choice of the degradation rate of V and S. Using a profiled likelihood approach [START_REF] Murphy | On profile likelihood[END_REF], we fixed their half-life to 0.25 and 51 days, respectively.

Inter-individual variability. In the vector η, some parameters have to be individual-specific to account for inter-individual variability. It is the case for ψ i = (ϑ, f M2 , γ). We suppose it follows a log-normal distribution such that:

ψ i = ψ 0 exp(u i ),
where ψ 0 is the fixed effect and average mean value in the population. The vector u i is individual random effects, which follow a normal distribution of mean zero and standard deviation Ω, and account for heterogeneity across individual. We assume that other parameters in vector η except error measurements are also estimated in log-transformation and are common to all individuals in the population. Altogether, the vector of parameters to estimate is given by θ = (η, Ω, σ BAU , σ ν ).

Estimation procedure Parameters were estimated (and named θ in the following) with the SAEM algorithm implemented in MONOLIX software version 2022R1 [START_REF] Lixoft | Monolix version 2020R1[END_REF] allowing to handle left censored data [START_REF] Samson | Extension of the SAEM algorithm to left-censored data in nonlinear mixed-effects model: Application to HIV dynamics model[END_REF]. Likelihood was estimated using the importance sampling method and standard error were obtained by asymptotic approximation and inversion of the Fisher Information Matrix. Graphical and statistical analyses were performed using R version 3.4.3.

Simulation of long-term humoral response

Next, we used the model to predict the long-term evolution of Ab and ED ν 50 over time. To account for uncertainty in our predictions, we used a Monte-Carlo sampling method, where K = 1000 replicates of parameters values θ (k) were sampled in the posterior distribution of the parameter estimates to derive 95% prediction intervals (PI) of the predicted trajectories. December 14, 2022 8/19

Finally we used these predictions to calculate the time to reach a given threshold value. To take into account between-subjects variability, we added a second layer to our Monte-Carlo sampling method and we sampled N = 100 replicates in the population parameter distribution. We used these predictions to derive the probability of having a concentration of binding antibodies greater than given thresholds, in particular higher than 264 BAU/mL, which corresponds to the standard threshold of protection defined by Feng et al. [START_REF] Feng | Correlates of protection against symptomatic and asymptomatic SARS-CoV-2 infection[END_REF] and adopted by WHO. The level of neutralizing activity has been identified as a correlate of protection for vaccine efficacy against the historical strain [START_REF] Gilbert | Immune correlates analysis of the mRNA-1273 COVID-19 vaccine efficacy clinical trial[END_REF][START_REF] Alexandre | Modelling the response to vaccine in non-human primates to define SARS-CoV-2 mechanistic correlates of protection[END_REF]. However, to date, no threshold for ED ν 50 value ensuring protection has been isolated for D614G, let alone for the new VoCs. So, for a range of threshold values, we calculated the probability that the neutralizing activity against each VoC would be higher than these values over time, especially if this activity was still detectable at a given time. In this way, we can compare the longevity of neutralizing activity between VoCs even in the absence of a clear threshold of protection for each of them.

Results

Mechanistic model for humoral response

We first aimed to investigate whether there is a proportional relationship between the evolution of concentration of binding antibodies and its neutralization capacity. for each VoC after each injection. First, we notice that these ratios are different for the variants. Then, using a Student's t test, we compared the evolution of these ratios with respect to the previous vaccination. In most cases, the ratios improved significantly, indicating an intrinsic gain in neutralisation that cannot be explained by the variation in antibody concentration alone, justifying the need to quantify this phenomenon precisely. Estimation of model parameters can be found in Table 2. This estimation indicates that multiple injections both increase antibody concentration and intrinsic affinity per constant antibody unit. Regarding antibody concentration, estimation of mechanistic parameters indicates a significant increase in the size of the memory compartment. It increased by f M2 = 4.5 (95% Confidence interval CI [2.9; 7.0]) after the second injection and by f M3 = 12.4 (95% CI [11.0; 14.0]) after the third injection compared to the first one.

Regarding neutralization per constant antibody concentration unit, we found that there are two main influencing factors: the repetition of the injections and the VoC.

Regarding repeated injections effect for the original strain D614G, the second dose increases neutralization by a factor f 2 = 8.3 (95% CI [5.6; 12.3]) and the third one by f 3 = 17 (95% CI [9.7; 28.4]) compared to the first injection. Now regarding the neutralization capacities for emerging VoCs, they are significantly decreased compared to the original strain, with the exception of Alpha, where there is no significant change in neutralization compared to D614G. It ranges from a reduction of 70% (95% CI [60%; 80%]) for Delta to a dramatic reduction of 98.4% (95% CI [97.7%; 99%]) for BA.5. Still, we find that the sequential injection strategy confers a gain in long-term neutralizing capacities for all VoCs. The second injection increases neutralization against all VoCs by the same factor f 2 = 8.3 (same as D614G). The third injection increases neutralization in a VoC-specific manner, given by f 3 g ν . It ranges from an increase in fold change of 21 (95% CI [4; 94]) for Delta to 230 (95% CI [40; 1200]) for BA.1 times higher for the third injection than for the first injection. For comparison with D614G, the neutralization is for an average individual (i.e random effects u i set to 0). As expected, the level of the response is higher after a repeated number of injections for both binding antibody concentration and neutralization for all variants. Interestingly, the neutralization curves for BA.1, BA.2, and BA.5 are significantly lower than for the other variants, with no overlap in prediction intervals. The first and second doses elicit a neutralization response for Omicron (BA.1, BA.2, and BA.5) that remains below the detection limit in most individuals (which is consistent with the observed data) but is dramatically enhanced by the third injection. Regarding the ED ν 50 /BAU ratio (See S1 Fig), we find that for the same concentration of binding antibodies, neutralization is significantly increased after each new injections for all variants and is significantly different for Alpha, D614G, {Delta, Beta} and {BA.1, BA.2, BA.5} variants. [START_REF] Feng | Correlates of protection against symptomatic and asymptomatic SARS-CoV-2 infection[END_REF] of 264 BAU/mL each days after the last injection in the counterfactual scenario where subjects only received one, two or three doses. The same is done for neutralizing activity again the VoCs (bottom, left: one, middle: two or right: three). It is possible to see the drastic effect of repeated injections on the levels reached by both binding antibodies concentration and neutralization for all variants. Strikingly, the full response duration is similar in length for the binding antibodies concentration after two or three injections. However, whereas response higher than 264 BAU/mL is reached in 100% (95% PI [99%; 100%]) of the population after three injections, it is only reached in 75% (95% PI [60%; 90%]) of the population after two doses, and not even reached in 2% (95% PI [0%; 3%]) of the population after the first injection. Table 3 provides the time needed for a proportion of a vaccinated population to return under a certain threshold. It explores multiple thresholds (100 BAU/mL, 264 BAU/mL, and 1000 BAU/mL for antibodies concentrations; and undetectability, 100 and 1000 for neutralization) that could be investigated when and if a clear level of correlate of protection is found. For all markers, there is a systematically and significantly higher duration of humoral activity after three compared to two injections. After three injections, duration of neutralization against Omicron variants (BA.1,BA.2 and BA.5) is significantly lower than for other variants for all thresholds. December 14, 2022 12/19

Long-term predictions

As already shown in

Discussion

We proposed here a modeling framework to characterize the kinetics of antibodies to successive doses of Bnt162b2 vaccine. The originality of our approach is that we relied on both the kinetics of anti-S IgG binding antibodies and their neutralization against the major VoCs that have emerged since 2021. Our model quantifies the benefit of successive injections and can be used to predict the optimal duration against each VoC.

The model shows a large and significant benefit of the third dose against all VoCs in terms of neutralizing activity. However, both the maximum level achieved and the rate of decline could vary greatly between VoCs. Accordingly, the mean duration of detectable neutralizing activity after the third dose of vaccine was 20, 12.5, 9, 8.5 and 6 months for Alpha, Delta, and Omicron BA.1, BA.2 and BA.5 respectively. Our results also highlight the wide variability in patient response, with at least 5% of patients with undetectable neutralizing activity against Omicron BA.5 only 3 months after the third December 14, 2022 13/19 injection.

These results were obtained based on a number of hypotheses, which we summarize below. First, the model of antibody concentration dynamics remains simplified, with the memory compartment simply represented by a piecewise constant function over successive doses. In addition, the model assumes only one type of secreting cell population and thus overlooks the complexity of the B-cell response mechanism. A piecewise constant function was also used to model the effects of successive doses on neutralization activity and did not account for other potential phenomena, such as time since last injection. Finally, the model assumed that the second dose would result in a similar change in protection for all variants. In the future, application of such approaches to larger populations of individuals may allow some of these hypotheses to be relaxed without compromising the identifiability of the parameters.

We restrict our analysis to vaccination induced humoral response and discarded patient data after breakthrough infection. To this date, 11 out of 26 followed subjects were infected with Omicron and this proportion is likely to increase as it is the case in the global vaccinated population. Thus, it is of great interest to model the hybrid protection induced by vaccination followed by natural infection. Still, it requires to deeply modify our model to integrate two different antibody populations, one coming December 14, 2022 14/19 from vaccination and targeting the historical strain and the other one targeting the Omicron variant. That is why this analysis is left to future works.

One of the main advantages of the model is its flexibility to easily incorporate information on new VoCs and to use the strength of information obtained on other viral variants to update the model as data become available. In fact, we continuously updated the model to include successive Omicron variants. Interestingly, despite the small number of samples available, a high degree of precision was achieved for all variants. For example, although patients had on average only two data points with detectable neutralization against the BA.2 variant, this was sufficient to achieve a good precision of the model parameters (Table 2). Moreover, by using all available data (eg, by analyzing anti-S IgG and neutralizing activity of all patients simultaneously), the model reveals some signals of kinetics that were not visible when analyzing the individual markers separately. For example, we identified different slopes of antibody decline that directly affect the prediction of protection duration. Using the same data set and a simpler single-slope model, the time to undetectable neutralizing levels after the third dose of vaccine for D416G was estimated to be 11.5 months [START_REF] Planas | Duration of BA. 5 neutralization in sera and nasal swabs from SARS-CoV-2 vaccinated individuals, with or without Omicron breakthrough infection[END_REF], which is much shorter than our estimate of 15 months (derived from Fig 4). In their approach, Planas et al. [START_REF] Planas | Duration of BA. 5 neutralization in sera and nasal swabs from SARS-CoV-2 vaccinated individuals, with or without Omicron breakthrough infection[END_REF] chose to adjust the anti-S IgG and ED ν 50 decline separately for the different VoCs without considering causal relationships between them. On the contrary, our model assumes an influence of antibody concentration on the development of ED ν 50 . Similar results are shown for the Delta variant (11.5 vs. 12.5 months, respectively), demonstrating the importance of a model-based approach to predict neutralizing activity in the long term. Predictions for the Omicron strains were similar for BA.1 and BA.2, respectively 9 and 8.5 months, and show a reduction for BA.5 strains with 6 months. Interestingly, our results also suggest a longer duration of detectable neutralizing activity than what has been directly extrapolated from other observational studies [START_REF] Brockman | Reduced magnitude and durability of humoral immune responses to COVID-19 mRNA vaccines among older adults[END_REF][START_REF] Mwimanzi | Older adults mount less durable humoral responses to two doses of COVID-19 mRNA vaccine, but strong initial responses to a third dose[END_REF][START_REF] Mwimanzi | Brief Report: Impact of age and SARS-CoV-2 breakthrough infection on humoral immune responses after three doses of COVID-19 mRNA vaccine[END_REF], although this difference may also be due in part to the different experimental procedure used to measure neutralization.

The large and VoC-dependent variability in neutralization values argues for the use of individualized approaches to identify patients most at risk. Although it should be acknowledged that such an approach is hampered by the lack of an established neutralizing activity threshold as correlates of protection, its level was found to be associated with the risk of breakthrough infection. In a cohort of elderly nursing home residents, none of those the individuals with ED50 above 2136 had Omicron BA.1 breakthrough infection [START_REF] Bruel | Neutralising antibody responses to SARS-CoV-2 omicron among elderly nursing home residents following a booster dose of BNT162b2 vaccine: A community-based, prospective, longitudinal cohort study[END_REF]. A model-based study found that a threshold of 1000 dramatically reduced peak viral load, suggesting that such a threshold may be a good indicator of protection against infection (unpublished results). Interestingly, our results show that neutralizing levels for all Omicron variants remain largely below this value (Fig 5 ), consistent with the current understanding that BNT162b2 is poorly effective against disease acquisition in the Omicron era [START_REF] Bruel | Neutralising antibody responses to SARS-CoV-2 omicron among elderly nursing home residents following a booster dose of BNT162b2 vaccine: A community-based, prospective, longitudinal cohort study[END_REF][START_REF] Lapointe | Serial infection with SARS-CoV-2 Omicron BA. 1 and BA. 2 following three-dose COVID-19 vaccination[END_REF]. Fortunately, the vaccine has shown high efficacy against severe disease to date [START_REF] Hacisuleyman | Vaccine breakthrough infections with SARS-CoV-2 variants[END_REF][START_REF] Kuhlmann | Breakthrough infections with SARS-CoV-2 omicron despite mRNA vaccine booster dose[END_REF].

To date, the use of a fourth dose of vaccine to increase efficacy in France has been limited to high-risk patients who were not represented in this cohort. Nevertheless, we used the model to predict the neutralization levels that could be achieved after a fourth vaccine dose. Under the conservative assumption, yet consistent with available observational study [START_REF] Regev-Yochay | Efficacy of a fourth dose of COVID-19 mRNA vaccine against omicron[END_REF], that this injection does not increase affinity or maturation parameters, our model predicts a similar duration of detectable neutralization as after the third dose, ranging from 183 to 268 days for the Omicron variants. Assuming that the fourth dose allows a similar increase in maturation and affinity as after the third dose, the model predicts that the duration of detectable neutralization could be much longer, ranging from 596 to 680 days for Omicron variants (see this supplementary analysis in Appendix C in S1 Appendix). December 14, 2022 15/19

Fig 1 .

 1 Fig 1. Top: longitudinal evolution of the binding antibody concentration of anti-S IgG. Bottom: longitudinal evolution of the neutralizing activity against VoCs after the first (left), second (middle) and third (right) vaccination dose. Squares represent median values, and plain horizontal lines represent the minimal and maximal encountered values among subjects. The lower limit of detection (LOD) is equal to 6 BAU/mL for IgG and 30 for ED 50 . Given the limited number of samples available, data were grouped, using a one week sliding window after the first dose, 20 days in the first 100 days following the second or third infection, and 50 days for the other data points.

Fig 2 displays

 2 the observed relationship from data between antibody concentration and ED ν 50

Fig 2 . 19 Fig 3 .

 2193 Fig 2. Evolution of the predicted ratio ED ν 50 /BAU for each VoC after successive vaccine doses. Each circle represents a ratio ED ν 50 /BAU computed when both measurements for ED ν 50 and BAU where available for a given patient at a given observation time. Most of patients contribute several times due to the repeated measurements made over time after each dose. All predictions below the limit of detection for ED ν 50 were removed to avoid overoptimistic ED ν 50 /BAU ratio when replacing ED ν 50 values by detection threshold. This explains why very few values are available for Beta and Delta and none for Omicron strains for one dose case. Comparison between vaccine dose were done using a Student test, "**" indicate p-values lower than 0.01, "*" lower than 0.1 and ">0.1" higher than 0.1.

  Fig 4, we can use the estimated models to predict the long-term trajectories of markers corresponding to the mean parameter values as well as 95% December 14, 2022 11/19

Fig 4 .

 4 Fig 4. Left: Predicted evolution of binding antibody concentration. The horizontal line corresponds to the value of 264 BAU/ml considered as a threshold against symptomatic infection. Right: Predicted kinetics of ED ν 50). The horizontal line corresponds to the LOD. In all panels, the shaded area is the 95% prediction interval.

Fig 5 (

 5 Fig 5 (top) displays the probability of having antibody concentration higher than the protection threshold established by[START_REF] Feng | Correlates of protection against symptomatic and asymptomatic SARS-CoV-2 infection[END_REF] of 264 BAU/mL each days after the last

Fig 5 .

 5 Fig 5. Top: Predicted probability of having predicted antibody concentration (anti-S IgG) greater than 264 BAU/mL. Bottom: Probability of having detectable neutralizing activity against VoCs after after the first (left), second (middle) or third (right) vaccination dose. Simulations were performed assuming that the second and third vaccination doses occurred at day 27 and 269, respectively.

Table

Table 1 .

 1 ). Briefly, all subjects received at least 2 doses, administered on average 27 days after the first injection. Characteristics of the analysed population

	Characteristics	Median	Median Time of vaccination
		[Min; Max] or n (%)	[Min; Max]
			since first dose since second dose
	Men	14 (54%)
	Age	59 [33; 95]
	Follow-up duration after first-dose (days)	368 [168; 483]
	Number of follow-up visits	14 [2; 18]
	Number of vaccination doses
	1st	26 (100%)	-	-
	2nd	26 (100%)	22 [17; 60]	-
	3nd	23 (88%)	243 [175; 385]	221 [154; 361]

N = 22 subjects received a third injection, administered on average 269 days after the first injection. During the follow-up N = 12 had a positive PCR, and only data prior to infection were analyzed, leaving an average follow-up of 11 visits and a median follow-up time of 362 days.

Ethics statement

This study was approved by the Ethics Committee ILE DE FRANCE IV. The cohort was approved by the national external committee (CPP Ile-de-France-IV IRB No. 00003835). Study participants did not receive any compensation.

  1 -top).

	December 14, 2022	3/19

  1 , t 2 , t 3 ). By denoting k = 1, . . . , 3 the dose number, on each interval [t k , t k+1 ], solving previous ODE for V gives us V k (t) = V 0 e -δ V (t-t k ) where V 0 is the initial antigen concentration, assumed equal from all doses.

	Because this model is not identifiable when only Ab are measured, we derived a
	structurally identifiable approximated model described in Equation 2; see Appendix A
	in S1 Appendix for a description of this simplification. Briefly, it consists of rescaling
	the model for S = µV 0 M 1	-1 S and assuming that M can be replaced by its
	steady-state value M k if equilibrium is reached quickly after each injection:

Table 2 .

 2 . Model parameters and estimation. Fcn:fold-change in neutralisation

	December 14, 2022

Acknowledgments

We thank Isabelle Staropoli, Florence Guivel-Benhassine, Françoise Porrot and all the members of the Virus and Immunity Unit for discussion and help, as well as Fabienne Peira, Vanessa Legros, Barbara De Dieuleveult, Aurelie Theillay, Sandra Pallay and Daniela Pires Roteia (CHR Orléans) for their help with the cohorts. Part of the experiments presented in this paper were carried out using the PlaFRIM experimental testbed, supported by Inria, CNRS (LABRI and IMB), Université de Bordeaux, Bordeaux INP and Conseil Régional d'Aquitaine (see https://www.plafrim.fr). We thank Simulations Plus, Lixoft division for the free academic use of the MonolixSuite.


Prof Roberto Lorusso
email: roberto.lorusso@mumc.nl

Maria Elena De Piero

MD Silvia Mariani

Michele Di Mauro

Thierry Folliguet

Fabio Silvio Taccone

Luigi Camporota

Justyna Swol

Dominik Wiedemann

Mirko Belliato

PhD Maria Elena De Piero

Lars Mikael Broman

Alain Vuylsteke

Yigal Kassif

Anna Mara Scandroglio

Vito Fanelli

Philippe Gaudard

Stephane Ledot

Julian Barker

Udo Boeken

Sven Maier

Alexander Kersten

Bart Meyns

Matteo Pozzi

M Finn

Peter Pedersen

Kaan Schellongowski

Nicholas Kirali

Barrett

Jordi Riera

Thomas Mueller

Jan Belohlavek

MD M Di

Mauro Phd

a multicentre, prospective observational study

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Since the end of 2019, the SARS-CoV-2 pandemic has affected more than 600 million people worldwide, with more than 6•5 million deaths (as of [START_REF] Who | COVID-19 Dashboard[END_REF].1 COVID-19 represents an unprecedented cause of hospitalisations, requiring intensive care unit (ICU) admission in 5-32% of admitted patients.1-3 Almost 90% of patients admitted to an ICU in the first wave of the pandemic required mechanical ventilation,2,3 with the most severe cases of acute respiratory distress syndrome (ARDS) requiring extracorporeal membrane oxygenation (ECMO).2-6 Despite the lessons learned from the 2009 H1N1 influenza pandemic,7 scepticism and uncertainty regarding appropriate selection of patients for ECMO were common during the initial phase of the first wave of COVID-19 from February to March, 2020. These doubts were exacerbated by scarcity of resources, differences in health-care systems, variability in institutional protocols, and cultural factors, such as heterogeneity concerning population characteristics (race and ethnicity, genetic background, and socioeconomic status).8 Once again, the patient selection process and ECMO management optimisation became crucial. Nevertheless, more than 2 years after the beginning of the COVID-19 pandemic, much remains unknown about the determinants of mortality, factors that could serve as a guide to prevent the futile initiation of extracorporeal support, or modifiable factors that could be targeted to improve outcomes. Even less evidence is available regarding outcomes after discharge;9-11 however, this information will be essential to understand the disease burden in patients with COVID-19 who survived ECMO. Furthermore, the post-COVID-19 condition, also known as long COVID, has been shown to have consequences for post-discharge quality of life and even delayed mortality after acute disease, which could play a part in the health trajectory after ECMO.

The EuroECMO-COVID study was developed by the European Extracorporeal Life Support Organization (ELSO) to provide weekly reports on patients undergoing ECMO for refractory COVID-19 in Europe and adjacent countries, including information on demographics, support type, and health status (deceased, alive, or still on ECMO).12, [START_REF] Lorusso | ECMO for COVID-19 patients in Europe and Israel[END_REF] The aim of this analysis of data from the EuroECMO-COVID study was to investigate inhospital outcomes in addition to survival and health status at 6 months after ECMO initiation-including respiratory, cardiac, and neurocognitive symptoms, and employment status-of patients with COVID-19 who were supported by ECMO during the first wave of the COVID-19 pandemic.

Methods

Study design and participants

The EuroECMO-COVID study is a prospective, international, multicentre observational study that was started in March, 2020, to follow the status of patients undergoing ECMO for refractory COVID-19 in Europe and adjacent countries, including Belarus, Israel, Norway, Russia, and Switzerland.12,13 Data collection is still ongoing. Currently, the study includes 204 participating centres and releases near real-time information on ECMO use in patients with COVID-19 every week throughout the European chapter of the ELSO (Euro-ELSO). Patients aged 16 years or older who received ECMO support for refractory COVID-19 from March 1, 2020, to Sept 13, 2020the first wave-were included in this analysis. Diagnosis of COVID-19 was confirmed at each centre through a PCR test on a nasopharyngeal specimen. All patients receiving ECMO for other concomitant causes were excluded.

This study complied with the Declaration of Helsinki, and the use of data was approved by the Ethics Committee of the Maastricht University Medical Centre, Maastricht, Netherlands (METC 2020-1343), which issued a consent exemption for the current analysis of anonymised data due to the increased benchmark risks during the pandemic. All participating centres had the study approved by the local ethics committee.

Procedures

The EuroECMO-COVID study collects data prospectively on patients undergoing ECMO for refractory COVID-19. Participating hospitals regularly deliver a predefined dataset to the coordinating centre and data collection is done on the basis of a predefined protocol.

Variables collected as part of the study dataset reported here addressed pre-ECMO profile, ECMO management, and outcomes until 6 months after initiation or death (appendix pp [START_REF] Melhuish | Outcomes of extracorporeal membrane oxygenation support for patients with COVID-19: a pooled analysis of 331 cases[END_REF][START_REF] Barbaro | Extracorporeal membrane oxygenation support in COVID-19: an international cohort study of the Extracorporeal Life Support Organization registry[END_REF][START_REF] Barbaro | Extracorporeal membrane oxygenation for COVID-19: evolving outcomes from the international Extracorporeal Life Support Organization Registry[END_REF][START_REF] Lebreton | Extracorporeal membrane oxygenation network organisation and clinical outcomes during the COVID-19 pandemic in Greater Paris, France: a multicentre cohort study[END_REF]. Data on in-hospital complications and causes of death were collected, including end-of-life decisions by the attending health-care personnel.

Follow-up continued until March 14, 2021, through direct patient contact or institutional charts review. Data on 6-month status included hospitalisation status, respiratory symptoms, neurocognitive complications (physical disability, muscle pain, attention disorders, and headache), and cardiovascular symptoms (chest pain, peripheral oedema, fatigue, and tachycardia). Finally, 6-month survivors were asked to report about their return to work in part-time or full-time employment.

Outcomes

In-hospital mortality and mortality 6 months after ECMO initiation were the primary outcomes. Secondary inhospital outcomes were the median duration of ECMO support, median duration of ICU and hospital stay, inhospital complications (eg, renal failure, major bleeding, neurological complications, gastrointestinal complications, and respiratory complications) from ECMO initiation to hospital discharge or death, cause of death, in-hospital time of death (on ECMO or after ECMO weaning), and heart or lung transplantation. Secondary 6-month outcomes were any kind of hospitalisation, work status, and respiratory, cardiac, and neurocognitive symptoms (appendix pp [START_REF] Melhuish | Outcomes of extracorporeal membrane oxygenation support for patients with COVID-19: a pooled analysis of 331 cases[END_REF][START_REF] Barbaro | Extracorporeal membrane oxygenation support in COVID-19: an international cohort study of the Extracorporeal Life Support Organization registry[END_REF][START_REF] Barbaro | Extracorporeal membrane oxygenation for COVID-19: evolving outcomes from the international Extracorporeal Life Support Organization Registry[END_REF][START_REF] Lebreton | Extracorporeal membrane oxygenation network organisation and clinical outcomes during the COVID-19 pandemic in Greater Paris, France: a multicentre cohort study[END_REF].

Statistical analysis

Demographic and clinical variables were expressed as n (%) for categorical variables and median (IQR) for continuous variables after evaluation for normality. Categorical data were compared with a χ² test. Continuous variables were analysed with a t test or a Mann-Whitney U test, as appropriate.

First, we described the population characteristics, ECMO details, and in-hospital outcomes for the whole cohort and separately for patients who died in hospital (in-hospital non-survivors) and those who did not die in hospital (in-hospital survivors). In-hospital survival was investigated with the Kaplan-Meier method and log-rank test for comparison. A patient's discharge or transfer to another facility was treated as a censoring event. To estimate the association between variables and mortality, we fitted a frailty mixed-Cox proportional hazards model that included both fixed and random effects. The random effect was used to consider differences between centres. [START_REF] Balan | A tutorial on frailty models[END_REF] We considered different sets of variables that we deemed to be important in terms of an association with inhospital mortality. Each variable set was developed on the basis of clinical practice and the authors' search of the literature. The mixed-Cox model estimated the hazard ratio (HR) for death in hospital as a function of the best discriminatory set of variables. Specifically, we included two models: one accounting only for pre-ECMO initiation variables along with initial ECMO mode, and another model with pre-ECMO and post-ECMO initiation variables.

The variables considered were as follows: age categories (<60 years, 60-69 years, or ≥70 years); sex; body-mass index (BMI) categories (≤29•9 kg/m² or obesity class 1 [30•0-34•9 kg/m²], class 2 [35•0-39•9 kg/m²], or class 3 [≥40•0 kg/m²] according to the WHO classification); race; history of chronic obstructive pulmonary disease (COPD), chronic renal failure, arterial hypertension, diabetes, or pre-existing cardiovascular disease; use of antiviral drugs, vasopressors, or inotropes; ECMO indication; pre-ECMO prone positioning; duration of pre-ECMO endotracheal intubation (<4 days or ≥4 days); initial ECMO mode (venovenous vs venoarterial); ECMO duration; ECMO configuration change; renal failure with and without renal replacement therapy; lung complications (pneumothorax, embolism, and superinfection); ischaemic stroke; major bleeding (intracranial, gastrointestinal, and lung bleeding); and bowel ischaemia. The cutoff of 4 days for pre-ECMO endotracheal intubation was found using time receiver operating characteristic package to identify the best predicting cutoff at different timepoints (0, 30, 60, 90, 120, 150, and 180 days). The 4-day cutoff showed the best combination of predictive parameters: sensitivity 50•7, specificity 54•1, positive predictive value of 56•2, negative predictive value of 48•7 at 180 months (see appendix p 26). The final models were then assessed for goodness of fit using Akaike's information criterion (AIC) and Bayesian information criterion (BIC) as reported. Discrimination power was assessed with Harrell's C index. A calibration plot was drawn using bootstrap resampling in 500 new samples and mean absolute prediction error was reported (appendix pp [START_REF] Urner | Venovenous extracorporeal membrane oxygenation in patients with acute covid-19 associated respiratory failure: comparative effectiveness study[END_REF][START_REF] Hall | Multi-institutional analysis of 505 patients with coronavirus disease-2019 supported with extracorporeal membrane oxygenation: predictors of survival[END_REF]. We chose only variables with missing data of 20% or less to be included in the model; for selected variables that still had missing data, we applied a multiple imputation process. Briefly, we used fully specified chained equations in the R package (version 1.1.463).15 Mechanisms underlying missing data were investigated using the Little test. Because missing data were not missing completely at random (MCAR) and monotonic pattern was not found, chained equations or the Markov chain Monte Carlo (MCMC) method was used, checking for healthy convergence. Ten imputed datasets were created and pooled (appendix pp 20-21). [START_REF] Van Buuren | mice: multivariate imputation by chained equations in R[END_REF] We plotted marginal distribution of observed data and ten imputed datasets per variable to verify whether the imputations were reasonable; in all cases, we achieved a reasonable imputation.

The parameters of interest from a mixed-Cox model was the hazard ratio (HR; 95% CI; p value) for death in hospital. Survival curves were obtained with Kaplan-Meier and log-rank tests for comparison. A p value of less than 0•0001 was considered statistically significant. A subgroup analysis focused on hospital survivors to describe their status 6 months after ECMO initiation. Each patient's follow-up period was defined as the time from the start of ECMO to death, 6-month follow-up, or last clinical contact. The 6-month outcomes are reported as counts and percentages of available data. Data completeness regarding the 6month follow-up is reported in the appendix (p 22).

All data were merged from separate deidentified files into SPSS (version 26.0), Prism (version 8.0), and R studio, for data management and statistical analysis. The study was registered with ClinicalTrials.gov, NCT04366921.

Role of the funding source

There was no funding source for this study.

Results

During the first COVID-19 pandemic wave, 133 centres from 21 countries (appendix pp [START_REF] Emanuel | Fair allocation of scarce medical resources in the time of Covid-19[END_REF][START_REF] Supady | Respiratory indications for ECMO: focus on COVID-19[END_REF] provided in-hospital and 6-month follow-up information for patients with COVID-19 who received ECMO. An additional 44 centres participated in the study but did not implant any ECMO in patients with COVID-19 during the first wave (figure 1 During the index hospitalisation, 852 (74%) of 1149 patients had major complications, with a different distribution between the 613 (50%) in-hospital survivors and the 602 (50%) inhospital non-survivors (table 3). The most frequent cause of death was multiorgan failure, followed by irreversible respiratory compromise, sepsis, neurological complication, acute pulmonary embolism, and cardiac arrest (figure 2). Patients who died had a median ICU stay of 24 days (IQR 14-39). Patients who did not die in hospital had a median ICU stay of 37 days (23-58) and an overall hospital stay length of 52 days (34-81; table 3).

Two multivariable models (one based on pre-ECMO parameters and the other on pre-ECMO and post-ECMO initiation parameters) were created with the best possible performance. The first model reported the association between pre-ECMO initiation parameters and in-hospital mortality: age of 60 years or older, chronic renal failure, use of inotropes or vasopressors before ECMO initiation, time from intubation to ECMO initiation of 4 days or more were found to be independent predictors (figure 3A, appendix p 28). The model showed good performance (AIC 149; BIC 31; Harrell's C 0•67; and mean absolute error 0•12). In the model that used pre-ECMO and post-ECMO initiation parameters, in addition to the associations shown in the first model, the need for ECMO configuration change, renal failure with or without renal replacement therapy, major bleeding, ischaemic stroke, and bowel ischaemia were also associated with higher in-hospital mortality (figure 3B, appendix p 29). The second model also showed good performance (AIC 201; BIC 41; Harrell's C 0•69; and mean absolute error 0•020).

Discussion

We report patient characteristics, in-hospital outcomes, and the associations between variables, such as patient demographics, comorbid status, pre-ECMO settings, ECMO characteristics and complications, and in-hospital mortality in 1215 patients with COVID-19 who were undergoing ECMO support, and we describe the largest 6-month follow-up during the first wave of the pandemic in this patient setting. Our study had four main findings. First, in patients who required ECMO for SARS-CoV-2 infection during the first pandemic wave, in-hospital mortality was high (50%), with death occurring after a median ICU stay of 24 days. Second, age of 60 years or more, use of inotropes and vasopressors before ECMO initiation, chronic renal failure, and a time from intubation to ECMO initiation of 4 days or more were associated with higher risk of death in hospital. Third, of the patients who were discharged, 5% died within 6 months of ECMO initiation. Finally, a substantial number of patients still had respiratory, cardiac, or neurocognitive symptoms and required hospitalisation or rehabilitation at 6-month follow-up. Only 24% of in-hospital survivors returned to full-time employment at 6 months.

Scepticism initially characterised the use of ECMO for patients with COVID-19 with respiratory failure. However, it soon became clear that selected patients might benefit from extracorporeal support.12,13,16-21 A first analysis of the ELSO Registry, reporting on 1035 patients with COVID-19 undergoing ECMO during the first pandemic wave, showed a 90-day mortality of 38%, with 68 patients still in hospital at the time of data analysis.17 Despite these early promising results, the first pandemic wave was characterised by uncertainty regarding appropriate COVID-19 therapies and optimal ECMO selection criteria. Moreover, the ethical circumstances surrounding the allocation of scarce resources during the first COVID-19 wave tended to favour a utilitarian approach, minimising the use of ECMO.22,23 In addition, significant differences were noticed in local and regional mortality rates, duration of ECMO support, and resource constraints.8,18 Our study highlights the substantial amount of resources required to manage a patient receiving ECMO during the first wave of the pandemic and the complex in-hospital outcomes of these patients. Almost half of the patients in this study died after a median of 14 days on ECMO support, an ICU stay of 24 days, and an overall hospital stay of 30 days. Patients who did not die were supported on ECMO for a median of 16 days; they remained in the ICU for 37 days and in hospital for 52 days.

In times of restricted resources and challenges in the ethical allocation of treatments, the early recognition of determinants of mortality and the identification of safe and effective treatments became fundamental to patient care. Initially, ELSO provided well defined indications and patient selection criteria for the use of ECMO, based on previous experience.4 However, COVID-19 is a new disease, and although several lessons from the past can be adapted to its management, many other aspects of care need special consideration.24 Time from intubation to ECMO start has previously been shown to be a determinant of patient outcome, leading to the ELSO guideline of 7 days as the upper time limit for ECMO consideration.7 Although a single-centre experience25 and a multi-institutional analysis26 questioned such a time threshold, our findings confirm that an even shorter time from intubation to ECMO initiation (<4 days) might be a determinant of death in patients with COVID-19 receiving ECMO. Similar results have been described by Hall and colleagues,27 who estimated that an ECMO patient intubated on day 14 after the diagnosis of COVID-19 versus day 4 had a relative hazard of survival of 0•65 (95% CI 0•44-0•96). In our analysis, only 54% of patients had a tracheostomy. More patients who did not die in hospital underwent a tracheostomy, suggesting that it was a practice reserved for patients with signs of positive outcomes.28 According to the recommendations by several major societies, tracheostomy was done 2-3 weeks after intubation only in selected patients because of the initial alarming high death rate in ICUs and the procedural risks for the attending staff.29 Finally, our findings confirmed that renal failure represents risk factors for in-hospital mortality;17 bowel ischaemia, neurological complications, and major bleeding also characterised a less favourable prognosis for these patients. By contrast, our data do not provide confirmation for chronic respiratory disease and cardiac arrest before ECMO as a risk factor for in-hospital mortality, as shown by data in the ELSO Registry.17

The post-discharge outcomes of patients treated for COVID- [START_REF] Lebreton | Extracorporeal membrane oxygenation network organisation and clinical outcomes during the COVID-19 pandemic in Greater Paris, France: a multicentre cohort study[END_REF] are not yet recorded in sufficient detail, particularly in patients who have received ECMO.31,32 Overall, survivors of ECMO are reported to have adverse changes in quality of life, with 12% of patients losing their jobs and almost 30% having significant disabilities affecting their survival.33,34 In patients with COVID-19, the effects on quality of life might be worsened by long COVID,31,35,36 which is characterised by severe fatigue, muscle pain, breathlessness, heart palpitations, and problems with attention, memory, and cognition.23-28 Our study confirms the hypothesis of long-lasting symptoms in patients who receive ECMO, despite a very low post-discharge mortality rate at 6 months for hospital survivors. In line with the results of previous smaller reports,9 only 24% of patients in our study returned to full-time work and 13% of patients returned to part-time work at 6 months. Furthermore, a substantial number of patients still required oxygen therapy or had shortness of breath and cardiac or neurocognitive symptoms. Whether these symptoms can be attributed to ECMO, COVID-19 itself, or a combination of both still needs to be clarified with larger specific studies. Nevertheless, our findings highlight the importance of a multidisciplinary post-ECMO follow-up programme tailored to address physical and psychological rehabilitation in patients who have received ECMO, especially in the COVID-19 era.9,34,37

The structured and systematic prospective data collection, regular weekly reports, and participation of large-volume centres in the EuroECMO-COVID study guarantee data robustness and granularity. Nevertheless, our study is observational by nature, limiting causal inferences. Moreover, under-reporting of ECMO implants or complications is plausible. Furthermore, specific data on ECMO details, such as selection criteria or local ECMO protocols, are not captured by the database and could not be included in this analysis. Indeed, each country might have followed different guidelines according to the local health-care system, resources, and ethical guidance.8 Perception of futility or the legal possibility of withdrawing support might have varied between centres. Nevertheless, this issue was addressed in our mixed-Cox analysis by using random effect to consider the differences between centres. Patients who did not die were either discharged home or transferred to another hospital, but this specification is missing, and this information might have been additionally informative. Furthermore, data on baseline functional status, previous physical limitations or working status, and an in-depth analysis of the patient's functional status at 6 months-including the 6 min walk test or spirometry, and an extended evaluation of psychological status, cognitive impairment, and quality of life through the use of specific tools-were either not available or not possible. Due to the particularity of this cohort in pandemic conditions, findings could differ in analyses of subsequent waves of COVID-19, although the EuroELSO study has reported a provisional comparison of in-hospital survival between the first and second waves (lower survival rate during the second wave).12 A partial overlapping of the cohort reported here with previously reported series cannot be excluded.17,18 In particular, we estimate a 15-20% overlap between this study and the ELSO Registry analysis of the first COVID-19 wave.17 Finally, a comparison between patients with COVID-19 and those with other causes of hypoxaemic respiratory failure-with potential differences in intrinsic pathophysiological characteristics, clinical presentation, and management-is not provided here due to the well-known differences that distinguish them.

In our study, in-hospital mortality after initiation of ECMO for COVID-19 during the first wave of the pandemic reached almost 50%. On the basis of their association with in-hospital mortality, patient's age, timing of cannulation (<4 days or ≥4 days from intubation), and use of inotropes or vasopressors represent essential factors to be considered in the care of patients, and the role of these factors needs to be confirmed in analyses of subsequent waves. Furthermore, identifying evolving cardiorespiratory compromise, guided by appropriate monitoring, represents an essential step in targeting a patient's physiological needs and a consequently adopting a more efficient ECMO configuration that could improve in-hospital results. The post-discharge prognosis in our study was favourable in terms of survival, but dedicated post-ECMO follow-up programmes with an emphasis on the psychological care of patients are warranted to support patients' recovery in view of the expected effects of persistent, disabling symptoms on health, employment status, and quality of life.

Contributors

Data sharing

The data dictionary and EuroELSO policies are available online. All submitting institutions have unrestricted access to their own data and routine reports generated using these data. Reports are produced by EuroELSO, which include standard outcome reports, data analyses as a result of specific requests to EuroELSO and scientific publications using submitted data. Any scientific study, abstract, presentation, or manuscript resulting from use of data requested from EuroELSO must be presented to the EuroELSO Scientific Committee or its delegate for review before presentation or submission for publication and must comply with the EuroELSO Publication Policy and Data Access. Data and reports will be provided through a secure transfer mechanism and comply with the GDPR EU regulation.

Figure 1 :

 1 Figure 1: Study profile

Figure 2 :

 2 Figure 2: Reasons for in-hospital deathAt 6-month follow-up, for the 613 in-hospital survivors, data collection completeness ranged from 70% (n=428) for data on returning to full-time work, to 94% (n=577) for mortality (appendix p 22). 30 (5%) of the 577 patients with data on 6-month survival died after hospital discharge. Lung transplantation was required for six (2%) of 309 hospital survivors and heart transplantation was required for three (1%) of 305 hospital survivors with reported data. One patient died soon after lung transplantation. 26 (5%) of 520 patients were still hospitalised or admitted to a nursing home facility at 6-month follow-up. Post-discharge respiratory rehabilitation was still required in 88 (17%) of 522 patients. Two patients (<1%) were still on mechanical ventilation, and 31 (6%) of 521 patients needed oxygen therapy. The most common residual symptoms at follow-up were dyspnoea (185 [35%] of 523 patients), cardiac symptoms (52 [10%] of 514), and neurocognitive symptoms (66 [13%] of 512) symptoms. 102 (24%) of 431 discharged patients had returned to full-time work; 57 (13%) of 428 were employed in part-time work (figure4).

Figure 3 :

 3 Figure 3: Mixed-Cox regression analysis of determinants of in-hospital mortality

Figure 4 :

 4 Figure 4: Status (A) and functional symptoms (B) of survivors 6 months after extracorporeal membrane oxygenation initiation

  ). Data were available for 1215 patients who were supported with ECMO (942 [78%] men and 267 [22%] women of 1209 patients for whom sex data were available; median age 53 years [IQR 46-60]) and diagnosed with SARS-CoV-2 infection. Patients' geographical distribution is presented in the appendix (p 23). Baseline patient characteristics and comparisons between in-hospital survivors and non-survivors are reported in table 1. Pre-ECMO use of inotropes was required in 361 (32%) of 1119 patients and vasopressors in 934 (79%) of 1178 patients (table 1). Time from hospital presentation to ICU admission was 7 days (IRQ 4-10), and the median time from intubation to ECMO initiation was 4 days (2-8). 1105 (91%) of 1215 patients received isolated venovenous ECMO; the remaining patients had combined cardiorespiratory or isolated cardiac support (table 2). The median ECMO support duration was 15 days (IRQ 8-27). Overall, 109 (11%) of 1039 patients required an ECMO reconfiguration: 23 (26%) of 87 due to refractory hypoxaemia, five (6%) of 87 due to left ventricular failure, seven (8%) of 87 due to right ventricular failure, and 19 (22%) of 87 due to biventricular failure. 590 (54%) of 1097 patients had tracheostomy; the median time from intubation to tracheostomy was 16 days (IQR 8-25; table 2).

  RL, MEDP, JBe, FST, LMB, TM, JS, JR, and NB conceived the study. SM and MDM conceived the analysis plan. RL, MEDP, SM, MDM, VLC, LMB, TM, and JBe designed the study. RL, MEDP, SM, and MDM drafted the initial manuscript. RL, MEDP, SM, MDM, TF, ST, LC, JS, DM, MB, LMB, AV, YK, AMS, VF, PG, SL, JBa, UB, SM, AK, BM, MP, FMP, PS, KK, NB, JR, TM, JBe, and all members of the EuroECMO-COVID Study Group analysed and interpreted the results, critically edited the manuscript, approved the final work, and agreed to be accountable for the accuracy and integrity of the work. RL, MEDP, SM, and MDM accessed and verified the data. The Chief Investigator and the Coordinating Scientific Committee had complete access to data. Other study members had access to the data upon approval of the Coordinating Scientific Committee and all participating centres, in compliance with the EuroELSO Publication Policy and Data Access and verified the data, RL, NB, TM, and JBe were responsible for the decision to submit the manuscript. The corresponding author confirms that all authors and contributors have seen and approved the final work.Julian Barker, Martin O Schmiady, Udo Boeken, David Schibilsky, Sven Maier, Alexander Kersten, Bart Meyns, Matteo Pozzi, Mariusz Kowalewski, Finn M Pedersen, Peter Schellongowski, Luis F Pinto, Pedro E Silva, Igor Kornilov, Kaan Kirali, Aaron Blandino Ortiz, Leen Vercaemst, Simon Finney, Peter Paul Roeleveld, Matteo Di Nardo, Felix Hennig, Marta Velia Antonini, Mark Davidson, Tim J Jones, Nicholas Barrett, Jordi Riera, Thomas Mueller, Jan Belohlavek. Affiliations for group members are provided in the appendix (pp 4-5).RL and JBe are past chairmen of the European Extracorporeal Life Support Organization (ELSO) Steering Committee. NB is the current chairman of the European ELSO Steering Committee. MVA, NB, JBe, MB, MDM, TM, JR, and JS are European ELSO Steering Committee members. RL reports personal fees from Abiomed and Xenios and honoraria payed institution from Medtronic, LivaNova, Getinge and Eurosets, unrelated to the submitted works. FST is advisory member for Eurosets unrelated to the submitted work. DW is proctor for Abbott and scientific advisor for Fresenius. NB is advisory board for ALung. JR reports lecture and advisory fees from Medtronic and Werfen outside the submitted work. PS reports advisory fees from Getinge but not related to the submitted work. AV received consultant fees from Inspira Technologies. MB is an advisory member for Eurosets, congress speaker for Hamilton Medical, and course speaker for Estor, all unrelated to the submitted work. LMB is advisory board for Eurosets and Xenios and is board member of EuroELSO Enterprise, unrelated to the present work. JBe reports consultancy and lecture fees from Abiomed and Getinge and is advisory board member of AstraZeneca and Boehringer Ingelheim. SL is advisory member of Inspira Technologies, unrelated to the present work. All other authors declare no competing interests.

	Declaration of interests
	EuroECMO-COVID Study Group
	Members of the EuroECMO-COVID Study Group (authors and additional contributors to this

study) were Roberto Lorusso, Maria Elena De Piero, Silvia Mariani, Michele Di Mauro, Valeria Lo Coco, Thierry Folliguet, Fabio Silvio Taccone, Luigi Camporota, Iwan C C van der Horst, Bas C T van Bussel, Ronny M Schnabel, Thijs Delnoij, Justyna Swol, Dominik Wiedemann, Mirko Belliato, Lars Mikael Broman, Alain Vuylsteke, Gil Bolotin, Yigal Kassif, Anna Mara Scandroglio, Vito Fanelli, Philippe Gaudard, Luca Lorini, Stephane Ledot,

Acknowledgments

We thank all contributors for their continuous commitment to sharing information and supporting each other during all phases of the COVID-19 pandemic. Additional contributors and collaborators are listed in the appendix (pp 4-12). There were no funding sources for this study.