Modelling the response to vaccine in nonhuman primates to define SARSCoV2 mechanistic correlates of protection
Abstract
The definition of correlates of protection is critical for the development of nextgeneration SARSCoV2 vaccine platforms. Here, we propose a modelbased approach for identifying mechanistic correlates of protection based on mathematical modelling of viral dynamics and data mining of immunological markers. The application to three different studies in nonhuman primates evaluating SARSCoV2 vaccines based on CD40targeting, twocomponent spike nanoparticle and mRNA 1273 identifies and quantifies two main mechanisms that are a decrease of rate of cell infection and an increase in clearance of infected cells. Inhibition of RBD binding to ACE2 appears to be a robust mechanistic correlate of protection across the three vaccine platforms although not capturing the whole biological vaccine effect. The model shows that RBD/ACE2 binding inhibition represents a strong mechanism of protection which required significant reduction in blocking potency to effectively compromise the control of viral replication.
Editor's evaluation
This work should be of interest to a broad readership in infectious diseases, especially those people interested in modeling of infections. It combines statistical and mechanistic modeling to find assayable correlates of immunity for vaccines. This method could be relevant to many diseases or vaccines, although the particular markers identified here likely will be limited in their generalizability.
https://doi.org/10.7554/eLife.75427.sa0Introduction
There is an unprecedented effort for SARSCoV2 vaccine development with 294 candidates currently evaluated (World Health Organization, 2021). However, variants of concern have emerged before the vaccine coverage was large enough to control the pandemics (Cobey et al., 2021). Despite a high rate of vaccine protection, these variants might compromise the efficacy of current vaccines (Kuzmina et al., 2021; Planas et al., 2021; Lustig et al., 2021; Zhou et al., 2021). Control of the epidemic by mass vaccination may also be compromised by unknown factors such as longterm protection and the need of booster injections in fragile, immunocompromised, elderly populations, or even for any individual if protective antibody levels wane. Furthermore, the repeated use of some of the currently approved vaccine could be compromised by potential adverse events or by immunity against vaccine viral vectors (Greinacher et al., 2021). Finally, the necessity to produce the billions of doses required to vaccinate the world’s population also explains the need to develop additional vaccine candidates.
The identification of correlates of protection (CoPs) is essential to accelerate the development of new vaccines and vaccination strategies (Koch et al., 2021; Jin et al., 2021). Binding antibodies to SARSCoV2 and in vitro neutralization of virus infection are clearly associated with protection (Khoury et al., 2021; Yu et al., 2020; Earle et al., 2021; Feng et al., 2021). However, the respective contribution to virus control in vivo remains unclear (Zost et al., 2020), and many other immunological mechanisms may also be involved, including other antibodymediated functions (antibodydependent cellular cytotoxicity [ADCC], antibodydependent complement deposition [ADCD], antibodydependent cellular phagocytosis [ADCP]; Yu et al., 2020; Mercado et al., 2020; Tauzin et al., 2021), as well as T cell immunity (McMahan et al., 2021). Furthermore, CoP may vary between the vaccine platforms (Plotkin, 2013; Plotkin, 2020; Bradfute and Bavari, 2011; Dagotto et al., 2020).
Nonhuman primate (NHP) studies offer a unique opportunity to evaluate early markers of protective response (MuñozFontela, 2020; Eyal and Lipsitch, 2021). Challenge studies in NHP allow the evaluation of vaccine impact on the viral dynamics in different tissue compartments (upper and lower respiratory tract) from day 1 of virus exposure (Yu et al., 2020; Mercado et al., 2020; Corbett et al., 2020). Such approaches in animal models may thus help to infer, for example, the relation between early viral events and disease or the capacity to control secondary transmissions.
Here, we propose to apply a modelbased approach on NHP studies to evaluate (i) the immune mechanisms involved in the vaccine response and (ii) the markers capturing this/these effect(s) leading to identification of mechanisms of protection and definition of mechanistic CoP (Plotkin and Gilbert, 2012). First, we present a mechanistic approach based on ordinary differential equation (ODE) models reflecting the virushost interaction inspired from models proposed for SARSCoV2 infection (Gonçalves et al., 2020; Kim et al., 2021; Gonçalves et al., 2021; Wang et al., 2020; Marc et al., 2021; Ke et al., 2021) and other viruses (Myers et al., 2021; Baccam et al., 2006; Goyal et al., 2019; Goyal et al., 2017). The proposed model includes several new aspects refining the modelling of viral dynamics in vivo, in addition to the integration of vaccine effect. A specific inoculum compartment allows distinguishing the virus coming from the challenge inoculum and the virus produced de novo, which is a key point in the context of efficacy provided by antigenspecific preexisting immune effectors induced by the vaccine. Then, an original data mining approach is implemented to identify the immunological biomarkers associated with specific mechanisms of vaccineinduced protection.
We apply our approach to a recently published study (Marlin et al., 2021) testing a proteinbased vaccine targeting the receptorbinding domain (RBD) of the SARSCoV2 spike protein to CD40 (αCD40.RBD vaccine). Targeting vaccine antigens to dendritic cells via the surface receptor CD40 represents an appealing strategy to improve subunitvaccine efficacy (Flamar et al., 2012; Zurawski et al., 2017; Cheng et al., 2018; Godot et al., 2020) and for boosting natural immunity in SARSCoV2 convalescent NHP.
We show that immunity induced by natural SARSCoV2 infection, as well as vaccineelicited immune responses contribute to viral load control by (i) blocking new infection of target cells and (ii) by increasing the loss of infected cells. The modelling showed that antibodies inhibiting binding of RBD to ACE2 correlated with blockade of new infections and RBDbinding antibodies correlate with the loss of infected cells, reflecting importance of additional antibody functionalities. The role of RBD/ACE2binding inhibition has been confirmed in two other vaccine platforms.
Results
A new mechanistic model fits the in vivo viral load dynamics in nasopharyngeal and tracheal compartments
The mechanistic model aims at capturing the viral dynamics following challenge with SARSCoV2 virus in NHP. For that purpose, we used data obtained from 18 cynomolgus macaques involved in the vaccine study reported by Marlin et al., 2021, and exposed to a high dose (1 × 10^{6} pfu) of SARSCoV2 administered via the combined intranasal and intratracheal route. The viral dynamics during the primary infections were characterized by a peak of genomic RNA (gRNA) production 3 days postinfection in both tracheal and nasopharyngeal compartments, followed by a decrease toward undetectable levels beyond day 15 (Figure 1—figure supplement 1). At the convalescent phase (median 24 weeks after the primary infection), 12 macaques were challenged with SARSCoV2 a second time, 4 weeks after being randomly selected to receive either a placebo (n=6) or a single injection of the αCD40.RBD vaccine (n=6) (Figure 1A). A third group of six naïve animals were infected at the same time. Compared to this naïve group, viral dynamics were blunted following the second challenge of convalescent animals with the lowest viral load observed in vaccinated animals (Figure 1B, Figure 1—figure supplement 2).

Figure 1—source data 1
 https://cdn.elifesciences.org/articles/75427/elife75427fig1data1v2.xlsx

Figure 1—source data 2
 https://cdn.elifesciences.org/articles/75427/elife75427fig1data2v2.xlsx

Figure 1—source data 3
 https://cdn.elifesciences.org/articles/75427/elife75427fig1data3v2.xlsx

Figure 1—source data 4
 https://cdn.elifesciences.org/articles/75427/elife75427fig1data4v2.xlsx

Figure 1—source data 5
 https://cdn.elifesciences.org/articles/75427/elife75427fig1data5v2.xlsx

Figure 1—source data 6
 https://cdn.elifesciences.org/articles/75427/elife75427fig1data6v2.xlsx

Figure 1—source data 7
 https://cdn.elifesciences.org/articles/75427/elife75427fig1data7v2.xlsx

Figure 1—source data 8
 https://cdn.elifesciences.org/articles/75427/elife75427fig1data8v2.xlsx

Figure 1—source data 9
 https://cdn.elifesciences.org/articles/75427/elife75427fig1data9v2.xlsx

Figure 1—source data 10
 https://cdn.elifesciences.org/articles/75427/elife75427fig1data10v2.xlsx

Figure 1—source data 11
 https://cdn.elifesciences.org/articles/75427/elife75427fig1data11v2.xlsx
We developed a mathematical model to better characterize the impact of the immune response on the viral gRNA and subgenomic RNA (sgRNA) dynamics, adapted from previously published work (Gonçalves et al., 2020; Kim et al., 2021; Baccam et al., 2006), which includes uninfected target cells (T) that can be infected (I_{1}) and produce virus after an eclipse phase (I_{2}). The virus generated can be infectious (V_{i}) or noninfectious (V_{ni}). Although a single compartment for de novo produced viruses (V) could be mathematically considered, two distinct ODE compartments were assumed for a better understanding of the model. We completed the model by a compartment for the inoculum to distinguish between the injected virus (V_{s}) and the virus produced de novo by the host (V_{i} and V_{ni}). In both compartments of the upper respiratory tract (URT), the trachea and nasopharynx, viral dynamics were distinctively described by this model (Figure 2A). Viral exchange between the two compartments was tested (either from the nasopharynx to the trachea or vice versa). However, as described in the literature (Gonçalves et al., 2021; Ke et al., 2020; Pinky et al., 2021) and demonstrated by the additional modelling work in Appendix 1 ‘Model building’, viral transport within the respiratory tract plays a negligible role in viral kinetics compared with viral clearance. Consequently, no exchange was considered in the model. Using the gRNA and sgRNA viral loads, we jointly estimated (i.e., shared random effects and covariates) the viral infectivity (β), the viral production rate (P), and the loss rate of infected cells (δ) in the two compartments. We assumed that gRNA and sgRNA were proportional to the free virus and the infected cells, respectively. This modelling choice relied on both biological and mathematical reasons (see section Materials and methods for more details). Due to identifiability issues, the duration of the eclipse phase (1/k), the clearance of free viruses from the inoculum (c_{i}) and produced de novo (c) were estimated separately by profile likelihood and assumed to be identical in the two compartments of the URT. In addition, infectious and noninfectious viruses were assumed to be cleared at the same rate. We estimated the viral infectivity at 0.95 × 10^{–6} (CI_{95%} [0.18 × 10^{–6}; 4.94 × 10^{–6}]) (copies/mL)^{–1} day^{–1} in naïve animals, which is in the range of previously reported modelling results whether in the case of SARSCoV2 virus (Kim et al., 2021; Wang et al., 2020) or influenza (Myers et al., 2021; Baccam et al., 2006). We found estimates of the loss rates of infected cells of 1.04 (CI_{95%} [0.79; 1.37]) day^{–1}, corresponding to a mean halflife of 0.67 day. This estimation was consistent with previously published results obtained on SARSCoV2 virus showing the mean value of this parameter ranging from 0.60 to 2 day^{–1} (i.e., halflife between 0.35 and 1.16 days) (Gonçalves et al., 2020; Kim et al., 2021; Gonçalves et al., 2021; Wang et al., 2020; Marc et al., 2021). The eclipse phase (3 day^{–1}) was found similar to the values commonly used in the literature (Gonçalves et al., 2020; Marc et al., 2021; Myers et al., 2021; Baccam et al., 2006). Here, we distinguished the clearance of the inoculum which was much higher (20 virions day^{–1}) as compared to the clearance of the virus produced de novo (3 virions day^{–1}). While the halflife of the virus de novo produced usually approximates 1.7 hr (i.e., c=10 day^{–1}) (Gonçalves et al., 2020; Gonçalves et al., 2021; Marc et al., 2021; Myers et al., 2021), because of this distinction, our model provided a higher estimation of 5.5 hr which remained in accordance with the estimations obtained by Baccam et al., 2006, on influenza A. Furthermore, the viral production by each infected cells was estimated to be higher in the nasopharyngeal compartment (12.1 × 10^{3} virions cell^{–1} day^{–1}, CI_{95%} [3.15 × 10^{3}; 46.5 × 10^{3}]) as compared to the tracheal compartment (0.92 × 10^{3} virions cell^{–1} day^{–1}, CI_{95%} [0.39 × 10^{3}; 2.13 × 10^{3}]). These estimations are in agreement with the observation of the intense production of viral particles by primary human bronchial epithelial cells in culture (Robinot et al., 2021). In particular, they are in the range of estimates obtained within the URT, either in NHP (Gonçalves et al., 2021) or in humans (Wang et al., 2020), with the product p × T_{0} equals to 15.1 × 10^{8} (CI_{95%} [3.98 × 10^{8}; 58.1 × 10^{8}]) and 0.21 × 10^{8} (CI_{95%} [0.088 × 10^{8}; 0.48 × 10^{8}]) virions mL^{–1} day^{–1} in the nasopharynx and the trachea, respectively. By allowing parameters to differ between animals (through random effects), the variation of cell infectivity and of the loss rate of infected cells captured the observed variation of the dynamics of viral load. The variation of those parameters could be partly explained by the group to which the animals belong reducing the unexplained variability of the cell infectivity by 66% and of the loss rate of infected cells by 54% (Supplementary file 1). The model fitted well the observed dynamics of gRNA and sgRNA (Figure 2B).

Figure 2—source data 1
 https://cdn.elifesciences.org/articles/75427/elife75427fig2data1v2.xlsx

Figure 2—source data 2
 https://cdn.elifesciences.org/articles/75427/elife75427fig2data2v2.xlsx

Figure 2—source data 3
 https://cdn.elifesciences.org/articles/75427/elife75427fig2data3v2.xlsx
Modelling of the dynamics of viral replication argues for the capacity of αCD40.RBD vaccine to block virus entry into host cells and to promote the destruction of infected cells
We distinguish the respective contribution of the vaccine effect and postinfection immunity on the reduction of the cell infection rate and the increase of the clearance of infected cells. Because blocking de novo infection and promoting the destruction of infected cells would lead to different viral dynamics profile (Figure 2—figure supplement 1), we were able to identify the contribution of each mechanism by estimating the influence of the vaccine compared to placebo or naïve animals on each model parameter. The αCD40.RBD vaccine reduced by 99.6% the infection of target cells in the trachea compared to the naïve group. The estimated clearance of infected cells was 1.04 day^{–1} (95% CI 0.75; 1.45) in naïve macaques. It was increased by 80% (1.86 day^{–1}) in the convalescent macaques vaccinated by αCD40.RBD or not.
The mechanistic model allows predicting the dynamics of unobserved compartments. Hence, a very early decrease of the target cells (all cells expressing ACE2) as well as of the viral inoculum which fully disappeared from day 2 onward were predicted (Figure 2C). In the three groups, the number of infected cells as well as infectious viral particles increased up to day 2 and then decreased. We show that this viral dynamic was blunted in the vaccinated animals leading to a predicted maximum number of infectious viral particles in the nasopharynx and the trachea below the detection threshold (Figure 2C). The number of target cells would be decreased by the infection in the naïve and the convalescent groups, whereas it would be preserved in vaccinated animals.
The RBDACE2binding inhibition is the main mechanistic CoP explaining the effect of the αCD40.RBD vaccine on new cell infection
In our study (Marlin et al., 2021), an extensive evaluation of the immunological response has been performed with quantification of spikebinding antibodies, antibodies inhibiting the attachment of RBD to ACE2, antibodies neutralizing infection, SARSCoV2specific CD4^{+} and CD8^{+} T cells producing cytokines and serum cytokine levels (Figure 3, Figure 1—figure supplements 3–5). Therefore, based on our mechanistic model, we investigated if any of these markers could serve as a mechanistic CoP. Such a CoP should be able to capture the effect of the natural immunity following infection, associated or not to the vaccine (group effect) estimated on both the rate of cell infection and the rate of the loss of infected cells. To this aim, we performed a systematic screening by adjusting the model for each marker and we compared these new models with the model without covariates and with the model adjusted for the groups. In particular, our approach allowed us to benefit from all the information provided by the overall dynamics of the immunological markers after the exposure by integrating them as timevarying covariates (see the Materials and methods section for a detailed description of the algorithm). We demonstrate that the RBDACE2binding inhibition measure is sufficient to capture most of the effect of the groups on the infection of target cells (Figure 4A and B). The integration of this marker in the model explains the variability of the cell infection rate with greater certainty than the group of intervention, reducing the unexplained variability by 87% compared to 66% (Supplementary file 1). The marker actually takes into account the variation between animals within the same group. Hence, it suggests that the levels of antiRBD antibodies induced by the vaccine that block attachment to ACE2 are highly efficient at reflecting the neutralization of new infections in vivo. Furthermore, when taking into account the information provided by the RBDACE2binding inhibition assay, the effect of the group of intervention was no longer significant (Supplementary file 1). Finally, we looked at the estimated viral infectivity according to the binding inhibition assay in each animal. A positive dependence was found between the viral infectivity and the RBDACE2binding inhibition measure, linking an increase of 10^{3} AU of the marker, whether over time or between animals, with an increase of 1.8% (95CI% [1.2%; 2.3%]) of the viral infectivity (see Supplementary file 4). Accordingly, the values at the time of exposure were not overlapping at all, distinguishing clearly the vaccinated and unvaccinated animals (see Figure 4C).
In the next step, several markers (IgGbinding antiRBD antibodies, CD8^{+} T cells producing IFNγ) appeared to be associated to the rate of loss of infected cells (Figure 4—figure supplement 1A). Both specific antibodies and specific CD8^{+} T cells are mechanisms commonly considered important for killing infected cells. We retained the antiRBD binding IgG Ab that were positively associated to the increase of the loss of infected cells. For unknown reason the IFNγ response was high in unstimulated conditions in the naïve group. Thus, although this marker was associated with a decrease of the loss rate of infected cells, it appears essentially here as an indicator of the animal group. Further studies would be needed to fully confirm the place of IFNγ response as a mechanistic marker.
A large part of the variation of the infection rate (71%) and loss rate of infected cells (60%) were captured by the two markers of CoP: the RBDACE2binding inhibition and the antiRBDbinding Ab concentration. Using the estimated parameters, the effective reproduction number could be calculated (R) which is representing the number of cells secondarily infected by virus from one infected cell (Figure 4D). When looking at this effective reproduction number according to the groups, the vaccinated animal presented from the first day of challenge an effective R below 1 meaning that no propagation of the infection started within the host. These results were consistent when taking the value of RBDACE2binding inhibition at the time of the challenge without considering the evolution of the inhibition capacity over time (Figure 4—figure supplement 1B). This means that the dynamics of the viral replication is impacted very early during the infection process in immunized (i.e., both convalescent and vaccinated) animals and that vaccinated animals were protected from the beginning by the humoral response. Then, we looked at the threshold of the markers of interest leading to the control of the withinhost infection (as defined by R<1) which was around 30,000 AU for the RBDACE2binding inhibition assay. For the animals in the naïve and the convalescent groups, the observed values of binding inhibition measured by ECL RBD (the lower the better) and of IgG antiRBDbinding antibodies (the higher the better) led to R>1, whereas in vaccinated animals, the value of ECL RBD led to R<1. Therefore, our modelling study shows that the inhibition of binding of RBD to ACE2 by antibodies is sufficient to control initial infection of the host (Figure 4E). According to the observed value of ECL RBD in vaccinated animals (e.g., 66 AU in Figure 4E), a decrease of more than 2 log_{10} of the inhibition capacity (to reach 81,000 AU), due to variant of concern (VoC) or waning of immunity, would have been necessary to impair the control of the withinhost infection. Moreover, a decrease of the neutralizing activity (i.e., increased ECL) could be compensated by an increase of cell death as measured by an increase of binding IgG antiRBD as a surrogate. As an example, increasing IgG antiRBD from 2.5 to 10 in the animal MF7 of the convalescent group would lead to a control of the infection.
In conclusion, the αCD40.RBD vaccineelicited humoral response leads to the blockade of new cell infection that is well captured by measure of the inhibition of attachment of the virus to ACE2 through the RBD of the spike protein. Hence, the inhibition of binding of RBD to ACE2 is a promising mechanistic CoP. Indeed, this CoP fulfills the three criteria of leading to the best fit (lower BIC), the best explanation of interindividual variability, and fully captured the effect of the group of intervention.
The model revealed the same CoP related to another proteinbased vaccine but not with mRNA1273 vaccine
We took the opportunity of another study testing a twocomponent spike nanoparticle proteinbased vaccine performed in the same laboratory and using the same immune and virological assays (Brouwer et al., 2021), measured only at the time of exposure, for applying the proposed model and methodology. In this study, six animals were vaccinated and compared to four naïve animals (Figure 5A and B). The good fit of the data (Figure 5C and D) allows for estimating the effect of the vaccine that appeared here also to decrease the infectivity rate (by 99%) and increase the clearance of the infected cells by 79%. Looking at the best mechanistic CoP following the previously described strategy, we ended here again with the inhibition of RBD binding to ACE2 as measured by ECL RBD. In fact, this marker measured at baseline before challenge fulfilled the three criteria: (i) it led to the best model in front of a model adjusted for group effect, (ii) it rendered the group effect nonsignificant, and (iii) it explained around 71% of the infectivity rate variability, compared to 65% of variability explained by the groups. Interestingly, here again, the inhibition assay led to a clear separation of the estimated rate of infectivity between vaccinees and the placebo group (Figure 5E).
Finally, we applied our approach to a published NHP study performed to evaluate several doses of mRNA1273 vaccine (Corbett et al., 2020). Using available data, we compared the viral dynamics in the 100 µg, 10 µg, and placebo groups, enrolling a total of 12 rhesus macaques in a 1:1:1 ratio. Similar to the previous study, only immune markers measured at the time of exposure were available in this study, in addition to viral dynamics. We started from the same model as defined previously. We estimated a reduction of the infection rate by 97% but we did not find any additional effect. Looking at potential mechanistic CoP, we retained neutralization as measured on live cells with Luciferase marker. Although this marker led to the best fit and replaced the group effect (which was nonsignificant after adjustment for the marker), it explained only 15% of the variability of estimated viral infectivity, while 19% were explained by the groups.
In conclusion, we demonstrated, based upon challenge studies in NHP vaccinated with two different proteinbased vaccine platforms, that both vaccines lead to the blockade of new cell infection. Neutralizing antibodies likely represent a consistent mechanistic correlate of protection (mCoP). This could change across vaccine platforms especially because mechanisms of action are different.
Discussion
We explored the mechanistic effects of three SARSCoV2 vaccines and assessed the quality of markers as mCoP. This model showed that neutralizing and binding antibodies elicited by a nonadjuvanted proteinbased vaccine targeting the RBD of spike to the CD40 receptor of antigen presenting cells are reliable mCoP. Interestingly, we found the simpler and easier to standardize and implement binding inhibition assay may be more relevant to use as a CoP than cellculture neutralization assays. This result has been replicated in another study testing a nanoparticle spike vaccine. The model was able to capture the effect of the vaccines on the reduction of the rate of infection of target cells and identified additional effects of vaccines beyond neutralizing antibodies. This latter consisted of increasing the loss rate of infected cells which was better reflected by the IgGbinding antibodies and CD8^{+} Tcell responses in the case of the CD40targeting vaccine. One limitation of our study is that the prediction potential of our model relies on the range of the immune markers measured. However, our approach would allow a full exploitation of the data generated as in systems serology where nonneutralizing Ab functions, such as ADCC, ADCP, ADCD, and Abdependent respiratory burst (ADRB) are explored (Chung et al., 2015). The role of ADCC in natural infection has been previously shown (Dufloo et al., 2021), ADCD in DNA vaccine recipients (Yu et al., 2020) and with Ad26 vaccine (Alter et al., 2021). Here, we extended significantly these data by modelling the viral dynamic, showing that two other proteinbased vaccines exert an additional effect on infected cell death which relied on the level of IgG antiRBDbinding antibodies especially for the CD40.RBDtargeting vaccine. Measurements of other nonneutralizing Ab functions would probably also capture this additional effect.
The next question after determining which marker is a valid mCoP is to define the concentration that leads to protection, looking for a threshold effect that will help to define an objective (Khoury et al., 2021; Jin et al., 2021). In the context of SARSCoV2 virus, several emerged variants are leading to a significant reduction of viral neutralization as measured by various approaches. However, a 20fold reduction of viral neutralization might not translate in 20fold reduction of vaccine efficacy (Emary et al., 2021). First, there are many steps between viral neutralization and the reduction of viral infectivity or the improvement of clinical symptoms. Second, the consequences of a reduction of viral neutralization could be alleviated by other immunological mechanisms not compromised by the variant. In the context of natural immunity, when the level of neutralizing antibodies was below a protective threshold, the cellular immune response appeared to be critical (McMahan et al., 2021; Chandrashekar et al., 2020). We showed with our model that an improvement of infected cell destruction could help to control the withinhost infection and is quantitatively feasible.
The control of viral replication is the key for reducing infectivity (Leung et al., 2020; Marks et al., 2021) as well as disease severity (Néant, 2021; Gutmann et al., 2021). According to our nonlinear model linking the neutralization to the viral replication, a decrease of 4 to 20fold in neutralization as described for the variants of concern (Planas et al., 2021; Zhou et al., 2021) is not enough, especially in the context of the response to CD40.RBDtargeting vaccine, to compromise the control of viral replication. The results showing a conserved effectiveness of mRNA vaccines in humans infected by the alpha or beta variants (Charmet et al., 2021), although a decrease of neutralization has been reported (Planas et al., 2021), are consistent with this hypothesis. However, this is highly dependent upon the mode of action of currently used vaccines and upon the VoC that may much more compromise the neutralization but being also intrinsically less pathogenic such as Omicron (Nyberg et al., 2022).
The analysis performed extended significantly the observation of associations between markers as previously reported for SARSCoV2 vaccine (Yu et al., 2020) and other vaccines (Kester et al., 2009) because it allows a more causal interpretation of the effect of immune markers. However, our modelling approach requires the in vivo identification of the biological parameters under specific experimentations. On the other hand, the estimation of parameters included in our model also provided information on some aspect of the virus pathophysiology. Notably, we found an increased capacity of virion production in nasopharynx compared to the trachea which could be explained by the difference in target cells according to the compartment (Travaglini et al., 2020). This result needs to be confirmed as it may also be the consequence of a different local immune response (Pizzorno et al., 2020). The choice of the structural model defining the hostpathogen interaction is a fundamental step in the presented approach. Here, it was well guided by the biological knowledge, the existing models for viral dynamics (Goyal et al., 2019; Gonçalves et al., 2021; Smith et al., 2018), and the statistical inference allowing the selection of the model that best fit the data. As the number of observations was relatively small in regard to the number of model parameters, we investigated overfitting issues. This was done using a bootstrap approach to evaluate the stability of confidence intervals of the estimated parameters. Results are provided in Appendix 2 ‘BICc as selection criteria and multiple testing adjustment’. Many modelling choices for the statistical model were made in this approach and more theoretical work evaluating the robustness of the results in their regards may be relevant for future works. In particular, we could relax the constraint of linear interpolation of marker dynamics by using simple regression models, allowing in the same time the integration of error model to account for measurement error for timevarying covariates (Dafni and Tsiatis, 1998; Carroll et al., 2006; Wu, 2009). Moreover, by construction, we assumed similar interindividual variability and effects of covariates within the two URT compartments as well as similar values for the viral infectivity and the loss rate of infected cells. Viral load dynamics measured in lungs being different from those in the URT (Lui et al., 2020; Goyal et al., 2020), the relaxation of this hypothesis of homogeneous physiological behavior in the URT may be pertinent to extend the model to the LRT. Finally, it should be underlined that the dynamics of the immune response has not been modelled as suggested for instance for Bcell response (Balelli et al., 2020). This clearly constitutes the next step after the selection of the markers of interest as done in the present work.
In conclusion, the modelling of the response to two new promising SARSCoV2 vaccines in NHP revealed a combination of effects with a blockade of new cell infections and the destruction of infected cells. For these two vaccines, the antibody inhibiting the attachment of RBD to ACE2 appeared to be a very good surrogate of the vaccine effect on the rate of infection of new cells and therefore could be used as a mechanistic CoP. This modelling framework contributes to the improvement of the understanding of the immunological concepts by adding a quantitative evaluation of the contributions of different mechanisms of control of viral infection. In terms of acceleration of vaccine development, our results may help to develop vaccines for ‘hardtotarget pathogens’, or to predict their efficacy in aging and particular populations (Pollard and Bijker, 2021). It should also help in choosing vaccine dose, for instance at early development (Rhodes et al., 2018) as well as deciding if and when boosting vaccination is needed in the face of waning protective antibody levels (Gaebler et al., 2021; Vanshylla et al., 2021), at least in NHP studies although the framework could be extended to human studies using mixed approaches of within and between hosts modelling (Goyal et al., 2022) providing that enough information is collected.
Materials and methods
Experimental model and subjects details
Request a detailed protocolCynomolgus macaques (Macaca fascicularis), aged 37–66 months (18 females and 13 males) and originating from Mauritian AAALAC certified breeding centers were used in this study. All animals were housed in IDMIT facilities (CEA, Fontenayauxroses), under BSL2 and BSL3 containment when necessary (Animal facility authorization #D9203202, Préfecture des Hauts de Seine, France) and in compliance with European Directive 2010/63/EU, the French regulations and the Standards for Human Care and Use of Laboratory Animals, of the Office for Laboratory Animal Welfare (OLAW, assurance number #A582601, US). The protocols were approved by the institutional ethical committee ‘Comité d’Ethique en Expérimentation Animale du Commissariat à l’Energie Atomique et aux Energies Alternatives’ (CEtEA #44) under statement number A20011. The study was authorized by the ‘Research, Innovation and Education Ministry’ under registration number APAFIS#244342020030216532863v1.
Evaluation of antispike, antiRBD, and neutralizing IgG antibodies
Request a detailed protocolAntispike IgG were titrated by multiplex bead assay. Briefly, Luminex beads were coupled to the spike protein as previously described (Fenwick et al., 2021) and added to a BioPlex plate (BioRad). Beads were washed with PBS 0.05% tween using a magnetic plate washer (MAG2x program) and incubated for 1 hr with serial diluted individual serum. Beads were then washed and antiNHP IgGPE secondary antibody (Southern Biotech, clone SB108a) was added at a 1:500 dilution for 45 min at room temperature (RT). After washing, beads were resuspended in a reading buffer 5 min under agitation (800 rpm) on the plate shaker then read directly on a Luminex Bioplex 200 plate reader (BioRad). Average MFI from the baseline samples were used as reference value for the negative control. Amount of antispike IgG was reported as the MFI signal divided by the mean signal for the negative controls.
AntiRBD and antinucleocapside (N) IgG were titrated using a commercially available multiplexed immunoassay developed by Mesoscale Discovery (MSD, Rockville, MD) as previously described (Johnson et al., 2020). Briefly, antigens were spotted at 200–400 μg mL^{–1} in a proprietary buffer, washed, dried, and packaged for further use (MSD Coronavirus Plate 2). Then, plates were blocked with MSD Blocker A following which reference standard, controls, and samples diluted 1:500 and 1:5000 in diluent buffer were added. After incubation, detection antibody was added (MSD SULFOTAGTM AntiHuman IgG Antibody) and then MSD GOLDTM Read Buffer B was added and plates read using a MESO QuickPlex SQ 120 MM Reader. Results were expressed as arbitrary unit (AU) mL^{–1}.
AntiRBD and antiN IgG were titrated by ELISA. The nucleocapsid and the spike RBD (Genbank # NC_045512.2) were cloned and produced in Escherichia coli and CHO cells, respectively, as previously described (Flamar et al., 2012). Antigens were purified on Ctag column (Thermo Fisher) and qualitycontrolled by SDSPAGE and for their level of endotoxin. Antigens were coated in a 96well plates Nuncimmuno Maxisorp (Thermo Fisher) at 1 μg mL^{–1} in carbonate buffer at 4°C overnight. Plates were washed in TBS Tween 0.05% (Thermo Fisher) and blocked with PBS 3% BSA for 2 hr at RT. Samples were then added, in duplicate, in serial dilution for 1 hr at RT. Noninfected NHP sera were used as negative controls. After washing, antiNHP IgG coupled with HRP (Thermo Fisher) was added at 1:20,000 for 45 min at RT. After washing, TMB substrate (Thermo Fisher) was added for 15 min at RT and the reaction was stopped with 1 M sulfuric acid. Absorbance of each well was measured at 450 nm (reference 570 nm) using a Tristar2 reader (Berthold Technologies). The EC50 value of each sample was determined using GraphPad Prism 8 and antibody titer was calculated as log (1/EC50).
The MSD pseudoneutralization assay was used to measure antibodies neutralizing the binding of the spike protein to the ACE2 receptor. Plates were blocked and washed as above, assay calibrator (COVID 19 neutralizing antibody; monoclonal antibody against S protein; 200 μg mL^{–1}), control sera, and test sera samples diluted 1:10 and 1:100 in assay diluent were added to the plates. Following incubation of the plates, an 0.25 μg mL^{–1} solution of MSD SULFOTAGTMconjugated ACE2 was added after which plates were read as above. Electrochemiluminescence (ECL) signal was recorded.
Viral dynamics modelling
Request a detailed protocolThe mechanistic approach we developed to characterize the impact of the immune response on the viral gRNA and sgRNA dynamics relies on a mechanistic model divided in three layers: first, we used a mathematical model based on ODEs to describe the dynamics in the two compartments, the nasopharynx and the trachea. Then, we used a statistical model to take into account both the interindividual variability and the effects of covariates on parameters. Finally, we considered an observation model to describe the observed log_{10} viral loads in the two compartments.
For the mathematical model, we started from previously published models (Gonçalves et al., 2020; Kim et al., 2021; Baccam et al., 2006) where the nasopharynx and trachea were respectively described by a target cell limited model, with an eclipse phase, as model of acute viral infection assuming targetcell limitation (Baccam et al., 2006). We completed the model by adding a compartment for the inoculum that distinguishes the injected virus (V_{s}) from the virus produced de novo (V_{i} and V_{ni}). To our knowledge, this distinction has not been proposed in any previous work. Two main reasons led us to make this choice. First, it allowed us to study the dynamics of the inoculum, in particular during the early phase of viral RNA load dynamics. Second, as described in more detail below, it gave us the opportunity to use all the information provided by the preclinical studies, such as the known number of inoculated virions, to define the initial conditions of the ODE model rather than estimating or randomly fixing them for V_{i} and V_{ni}, as is usually done. Consequently, for each of the two compartments, the model included uninfected target cells (T) that can be infected (I_{1}) either by infectious viruses (V_{i}) or inoculum (V_{s}) at an infection rate β. After an eclipse phase, infected cells become productively infected cells (I_{2}) and can produce virions at rate P and be lost at a per capita rate δ. The virions generated can be infectious (V_{i}) with proportion µ while the (1−µ) remaining proportion of virions is noninfectious (V_{ni}). Mathematically, a single compartment (V) for de novo produced virions could be considered in the model, with µV and (1−µ)V representing the respective contributions of infectious and noninfectious viruses to the biological mechanisms. However, to have a better visual understanding of the distinction between the two types of viruses, we wrote the model with distinct compartments, V_{i} and V_{ni}.
Finally, virions produced de novo and those from the inoculum are cleared at a rate c and c_{i}, respectively. Distinct clearances were considered to account for the effects of experimental conditions on viral dynamics. In particular, it is hypothesized that, animals being locally infected with large numbers of virions, a large proportion of it is assumed to be rapidly eliminated by swallowing and natural downstream influx, in contrast to the de novoproduced virions. However, it is important to keep in mind that this distinction was possible because of the controlled experimental conditions performed in animals, (i.e., exact timing and amount of inoculated virus known, and frequent monitoring during the early phase of the viral dynamics). Because of identifiability issues, similar clearances for infectious and noninfectious viruses were used. Accordingly, the model can be written as the following set of differential equations, where the superscript X denotes the compartment of interest (N, nasopharynx or T, trachea):
where ${T}^{X}\left(t=0\right)$, ${I}_{1}^{X}\left(t=0\right)$, ${I}_{2}^{X}\left(t=0\right)$, ${V}_{i}^{X}\left(t=0\right)$, ${V}_{ni}^{X}\left(t=0\right)$, and ${V}_{s}^{X}\left(t=0\right)$ are the initial conditions at the time of exposure. The initial concentration of target cells, that are the epithelial cells expressing the ACE2 receptor, is expressed as ${T}_{0}^{X}=\frac{{T}_{0}^{X,nbc}}{{W}^{X}}$, where ${T}_{0}^{X,nbc}$ is the initial number of cells and W^{X} is the volume of distribution of the compartment of interest (see the subsection ‘Consideration of the volume of distribution’). Each animal was exposed to 1 × 10^{6} pfu of SARSCoV2 representing a total of 2.19 × 10^{10} virions. Over the total inoculum injected (5 mL), 10% (0.5 mL) and 90% (4.5 mL) of virions were respectively injected by the intranasal route and the intratracheal route leading to the following initial concentrations of the inoculum within each compartment: $V}_{S,0}^{N}=\frac{0.10\times Ino{c}_{0}}{{W}^{N}$ and ${V}_{S,0}^{T}=\frac{0.90\times {\mathrm{I}\mathrm{n}\mathrm{o}\mathrm{c}}_{0}}{{W}^{\mathrm{T}}}$, with Inoc_{0} the number of virions injected via the inoculum.
Using the gRNA and sgRNA viral loads, we estimated the viral infectivity, the viral production rate, and the loss rate of infected cells within each of the two compartments of the URT (Supplementary file 2). To account for interindividual variability and covariates, each of those three parameters was described by a mixedeffect model and jointly estimated between the two compartments as follows:
where ${\beta}_{0},\mathrm{log}\left({\delta}_{0}\right)$, and $\mathrm{log}\left({P}_{0}\right)$ are the fixed effects, $\left\{{\varphi}_{conv}^{\theta}\mid \theta \in \left\{\beta ,\text{}\delta ,P\right\}\right\}$ and $\left\{{\varphi}_{CD40}^{\theta}\mid \theta \in \left\{\beta ,\delta ,P\right\}\right\}$ are respectively the regression coefficients related to the effects of the group of convalescent and αCD40.RBDvaccinated animals for the parameters β, δ, and P, and ${u}_{i}^{\theta}$ is the individual random effect for the parameter θ, which is assumed to be normally distributed with variance ${\omega}_{\theta}^{2}$ . A logtransformation was adopted for the parameters δ and P to ensure their positivity while a log_{10}transformation was chosen for viral infectivity to also improve the convergence of the estimation. Because of the scale difference between the parameter β and the other parameters (see Supplementary file 2), the mere use of the logtransformation for this parameter led to convergence issues. The use of a log_{10}transformation allowed to overcome this problem. Moreover, as shown in Equation 2, a joint estimation of the parameters β, δ, and P between the two compartments of the URT was considered. In this regard, a homogeneous interindividual variability within the URT was assumed as well as a similar contribution of the covariates to the value of the parameters. Parameters in the trachea were then either equal or proportional to those in the nasopharynx. This modelling choice, resulting in a smaller number of parameters to be estimated, was made mainly to address identifiability issues and to increase the power of the estimation. All other parameters included in the targetcell limited models were assumed to be fixed (see the subsection ‘Parameter estimation’ for more details).
In practice, after the selection of the optimal statistical model (see Appendix 1 ‘Model building’), random effects were added only to the parameters β and δ (i.e., ω_{β} ≠0, ω_{δ} ≠0, and ω_{P}=0), and the estimation of multiple models identified the viral production rate P as the only parameter taking different values between the trachea and nasopharynx. (i.e., β^{N}=β^{T} with f_{β}^{T}=0, δ^{N}= δ^{T} with f_{δ}^{T}=0, while P^{N}≠P^{T}). Finally, the adjustment of the model for the categorical covariates of groups of treatment, natural infection, and/or vaccination identified β and δ as the parameters with a statistically significant effect of these covariates (i.e., ${\varphi}_{conv}^{P}=0$ and ${\varphi}_{CD40}^{P}=0$).
For the observation model, we jointly described genomic and subgenomic viral loads in the two compartments of the URT. We defined genomic viral load, which characterizes the total viral load observed in a compartment (nasopharynx or trachea), as the sum of inoculated virions (V_{s}), infectious (V_{i}), and noninfectious virions (V_{ni}). The sgRNA was described as proportional to the infected cells (I_{1} + I_{2}). This choice was driven by two main reasons. First, sgRNA is only transcribed in infected cells (Sawicki et al., 2007). Second, as described by Miao et al., 2011, to overcome identifiability issues between the parameters β and P typically observed in targetcell limited models. The comparison of the two observation models describing sgRNA as either proportional to virions produced de novo (V_{i} +V_{ni}) or proportional to infected cells (I_{1} + I_{2}) confirmed this conclusion. In addition to a better BICc value (–25 points) compared with the first model, the second one allowed the estimation of both β and P by counteracting identifiability problems faced with the first model (results not shown). Accordingly, the log_{10}transformed gRNA and sgRNA of the ith animal at the jth time point in compartment X (nasopharynx or trachea), denoted ${gRNA}_{ij}^{X}$ and ${sgRNA}_{ij}^{X}$, respectively, were described by the following equations:
where $\mathrm{\Theta}}_{i}^{X$ is the set of parameters of the subject i for the compartment X and ε are the additive normally distributed measurement errors.
Consideration of the volume of distribution
Request a detailed protocolTo define the concentration of inoculum within each compartment after injection, nasopharyngeal and tracheal volumes of distribution, labelled W^{N} and W^{T}, respectively, were needed. Given the estimated volumes of the trachea and the nasal cavities in four monkeys similar to our 18 macaques (Figure 2—figure supplement 2A–C) and the welldocumented relationship between the volume of respiratory tract and animal weights (Asgharian et al., 2012), the volume of distribution of each compartment was defined as a step function of NHP weights:
where weight_{i} is the weight of the monkey i in kg. Using Equation 4 and weights of our 18 NHPs (mean = 4.08; [Q1; Q3] = [3.26; 4.77]), we estimated WT = 2 and WN = 4 mL for a third of them (n=12) (Figure 2—figure supplement 2D), leading to the initial concentration of target cells ${T}_{0}^{X}$ (see ‘Viral dynamics modelling’ for equation) fixed at 3.13 × 10^{4} cells mL^{–1} and 1.13 × 10^{4} cells mL^{–1} in nasopharynx and trachea, respectively. Similarly, their initial concentrations of challenge inoculum ${V}_{S,0}^{X}$ were fixed at 5.48 × 10^{8} copies,mL^{–1} and 9.86 × 10^{9} copies,mL^{–1} in nasopharynx and trachea respectively. For the last third of NHPs (n=6), WT = 3 and WN = 5.5 mL leading to ${T}_{0}^{X}$ fixed at 2.27 × 10^{4} cells mL^{–1} in nasopharynx and 7.50 × 10^{3} cells mL^{–1} in trachea while ${V}_{S,0}^{X}$ was fixed at 3.98 × 10^{8} copies mL^{–1} in nasopharynx and 6.57 × 10^{9} copies mL^{–1} in trachea. Through this modelling, we assumed a homogenous distribution of injected virions and target cells within nasopharyngeal and tracheal compartments. In addition, the natural downward flow of inoculum toward lungs, at the moment of injection, was indirectly taken into account by the parameter of inoculum clearance, c_{i}.
Parameter estimation
Request a detailed protocolAmong all parameters involved in the three layers of the mechanistic model, some of them have been fixed based on experimental settings and/or literature. That is the case of the proportion of infectious virus (µ) that has been fixed at 1/1000 according to previous work (Gonçalves et al., 2021) and additional work (results not shown) evaluating the stability of the model estimation according to the value of this parameter. The initial number of target cells, that are the epithelial cells expressing the ACE2 receptor, ${T}_{0}^{X,nbc}$ was fixed at 1.25 × 10^{5} cells in the nasopharynx and 2.25 × 10^{4} cells in trachea (Gonçalves et al., 2021; Supplementary file 2). The duration of the eclipse phase (1/k), the clearance of the inoculum (${c}_{i}$) and the clearance of the virus produced de novo (c) were estimated by profile likelihood. The profile likelihood consists in defining a grid of values for the parameters to be evaluated and sequentially fixing these parameters to one of these combinations of values. The model and all the parameters that are not fixed are then estimated by maximizing the loglikelihood. In this process, all parameters that are assumed to be fixed in the model (i.e., μ and the initial conditions) are held fixed. Finally, the optimal set of parameters is chosen as the one optimizing the loglikelihood. Although the available data did not allow the direct estimation of these three parameters, the use profile likelihood enabled the exploration of various potential values for k, c, and ${c}_{i}$ . In a first step, we explored the 18 models resulting from the combination of three values of k∈{1, 3, 6} day^{–1} and six values for c∈{1, 5, 10, 15, 20, 30} day^{–1}, assuming that the two parameters of virus clearance were equal, as first approximation. As shown in Supplementary file 3a, an eclipse phase of 8 hr (k=3) and virus clearance higher than 15 virions per day led to lowest values of –2loglikelihood (–2LL, the lower the better). In a second step, we fixed the parameter k at 3 day^{–1} and estimated the 70 models resulting from the combination of 10 values for c∈{1, 2, 3, 4, 5, 10, 15, 20, 25, 30} day^{–1} and 7 values for∈{1, 5, 10, 15, 20, 25, 30} day^{–1} (Supplementary file 3b). The distinction of the two parameters of free virus clearance enabled to find much lower halflife of inoculum (~50 min) than halflife of virus produced de novo (~5.55 hr), with c=3 day^{–1} compared to c_{i} = 20 day^{–1}.
Once all these parameters have been fixed, the estimation problem was restricted to the determination of the viral infectivity β, the viral production rate P, the loss rate of infected cells δ for each compartment, the parameter ${\alpha}_{vlsg}$ in the observation model, regression coefficients for groups of intervention $\left({\varphi}_{conv},{\varphi}_{CD40}\right)$, and standard deviations for both random effects ($\omega $) and error model (σ). The estimation was performed by maximum likelihood estimation using a stochastic approximation EM algorithm implemented in the software Monolix (http://www.lixoft.com). The Fisher information matrix was calculated by stochastic approximation, providing for each estimated parameter its variance, from which we were able to derive its 95% confidence interval. Selection of the compartment effect on parameters (β, δ, P) as well as random effects and covariates on the statistical model (Equation 2) was performed by the estimation of several models that were successively compared according to the corrected Bayesian information criterion (BICc) (to be minimized). After the removal of random effect on the viral production (${\omega}_{P}=0$) allowing the reduction of the variance on the two other random effects, all combinations of compartment effects were evaluated, leading to the final selection of a single effect on P $\left({f}_{\beta}^{T}={f}_{\delta}^{T}=0\right)$. Then, the effect of group intervention was independently added on model parameters among β, δ, P, and c. Once the group effect on the viral infectivity identified as the best one, the addition of a second effect on the remaining parameters was tested, resulting in the selection of the loss rate of infected cells. Finally, the irrelevance of the addition of a third effect was verified.
The possibility of migration of free plasma virus between the nasopharynx and the trachea was tested. However, as widely described in the literature, the transport of viral particles within the respiratory tract is negligible in the viral dynamics and is difficult to estimate. The reader can refer to Appendix 1 ‘Model building’ for an additional modelling work conducted to estimate this exchange and provided the same conclusion. Accordingly, the two compartments of the URT were assumed are distinct in our model.
Algorithm for automatic selection of biomarkers as CoP
Request a detailed protocolAfter identifying the effect of the group of intervention on both the viral infectivity (β) and the loss rate of infected cells (δ), we aimed at determining whether some immunological markers quantified in the study could capture this effect. Nowadays, many methods for selecting constant covariates already exist (Chowdhury and Turin, 2020) and are implemented in software like Monolix. However, these latter do not allow timevarying covariates. In this section, we present the algorithm we implemented to select timevarying covariates. We proposed a classical stepwise datadriven automatic covariate modelling method (Figure 4—figure supplement 2). However, initially implemented to select covariates from more than 50 biomarkers, computational time restricted us to consider only a forward selection procedure. Nevertheless, the method can be easily extended to classical stepwise selection in which both forward selection and backward elimination are performed sequentially. Although the method was developed for timevarying covariates, it can also be applied to constant covariates.
At the initialization step (k=0) (see Figure 4—figure supplement 2), the algorithm requests three inputs: (World Health Organization, 2021) a set of potential $M$ covariates, labelled Marker m for $m\in \left\{1,\cdots ,M\right\}$ (e.g., immunological markers); (Cobey et al., 2021) a set of P parameters on which covariates could be added, labelled θ_{p} for $p\in \left\{1,\cdots ,P\right\}$(e.g., β and δ); and (Kuzmina et al., 2021) an initial model (e.g., the model without covariates), labelled M^{0}, with $\theta}_{p}^{0$ being the definition of the parameter θ_{p}. At each step k>0, we note M^{k}^{−1} the current model resulting in the model built in the step k−1. Then, each combination of markers and parameters that have not already been added in M^{k}^{−1}, labelled $r(r\in \{\mathrm{M}\mathrm{a}\mathrm{r}\mathrm{k}\mathrm{e}\mathrm{r}\phantom{\rule{thinmathspace}{0ex}}m\phantom{\rule{thinmathspace}{0ex}}\u2a02$ $\theta}_{p}\notin {M}^{k1}\mid \phantom{\rule{thinmathspace}{0ex}$ $m\in \left\{1,\text{}\dots M\right\},\text{}p\in \left\{1,\dots P\right\}\})$, are considered and tested in an univariate manner (each relation r is independently added in M^{k}^{−1} and ran). To this end, the parameter θ_{p} involved in this relationship r is modified as ${\theta}_{p}^{k}\left(t\right)={\theta}_{p}^{k1}\left(t\right)\times exp({\varphi}_{m}^{p}\times Marke{r}_{m}\left(t\right))$, where ${\varphi}_{m}^{p}$ is the regression coefficient related to the marker and $Marke{r}_{m}\left(t\right)$ being the trajectory of the marker over time, while other parameters remain unchanged $\left(\forall {\theta}_{q}\notin r,{\theta}_{q}^{k}\left(t\right)={\theta}_{q}^{k1}\left(t\right)\right)$. Once all these models evaluated, the one with the optimal value of a given selection criterion defining the quality of the fits (e.g., the lowest BICc value) is selected and compared to the model M^{k}^{−1}. If the value of the criterion is better than the one found for M^{k}^{−1}, then this model is defined as the new current model, M^{k}, and the algorithm moves to the step k+1. Otherwise, the algorithm stops. The algorithm can also be stopped at the end of a fixed number of step k.
The objective of this algorithm being to identify mechanistic CoP, at each step, the selected model should respect, in addition to the best fits criterion, the two other criteria defining mCoP meaning the ability to capture the effect of the group of intervention and the ability to better explain the variability on individual parameters than the model adjusted for the group effect. To this end, we verify that in the selected model additionally adjusted for the group of intervention, the group effect appears as nonsignificantly different from 0 using a Wald test. Then, we check that the variances of random effects in the selected model are lower or equal to the ones obtained in the model adjusted only for the group effect.
Modelling hypothesis for timedependent covariates in our application
Request a detailed protocolUsing a populationbased approach to estimate our mechanistic model and similar to the adjustment of the model for constant covariates (e.g., groups of intervention), timevarying covariates are incorporated into the statistical model as individualspecific explanatory variables in the mixedeffects models. To implement the algorithm for selecting the timevarying covariates, many modelling choices were made. First, targeting covariates able to fully replace the group of intervention, we kept a similar mathematical relationship between parameters and immune markers than the one used with the constant covariate (see Equation 2). Accordingly, we adjusted the model parameters additively in logarithmic scale. In this regard, at each step k (k>0), the parameter θ_{p} was defined as $\mathrm{log}\left({\theta}_{p}^{k}\left(t\right)\right)$ $=\mathrm{log}\left({\theta}_{p}^{k1}\left(t\right)\right)$ $+{\varphi}_{m}^{p}\times Marke{r}_{m}\left(t\right)$. However, this choice may affect the results and other choices may be more relevant under different conditions. Second, because immune markers are observed only at discrete time points, whereas the estimation of the model is performed in a continuous way, we introduced immune markers as timevarying covariates using linear interpolation. Lets denote Marker_{i,j} the value of the marker observed for the ith animal at the jth time point, with $i\in \{1,\dots ,n\}$ and $j\in \{1,\dots ,J\}$. By linear interpolation, the timecontinuous marker was defined as, $\mathrm{\forall}t>0$,
As previously described in the Results section, three different studies were considered in this work: a main study reported by Marlin et al., 2021, testing the αCD40.RBD vaccine, and two additional studies (Corbett et al., 2020; Brouwer et al., 2021) evaluating a twocomponent spike nanoparticle vaccine and the mRN1273 vaccine, respectively. In the main study, the method was applied with both timevarying covariates and constant covariates for which only baseline value was considered, such that Marker_{i}(t)=Marker_{i}(t=0) (see Supplementary file 1). For the other two studies, only the baseline values were considered as covariates, the dynamics being not available. To assess the robustness of the results, several selection criteria were tested: AIC, BIC, loglikelihood, the percentage of explained interindividual variability, and similar results were obtained for all (results not shown). Moreover, as presented in Appendix 2 ‘BICc as selection criteria and multiple testing adjustment’, we verified the robustness of the use of BIC as selection criteria despite the multiplicity of the tests. The identification of antibodies inhibiting the attachment of the RBD to the ACE2 receptor (ECLRBD) as the first timevarying CoP led to the definition of the timevarying viral infectivity for the ith animal as described in Equation 5, while the selection antiRBD IgGbinding antibodies led to the elimination rate of infected cells given in Equation 6.
Quantification and statistical analysis
Request a detailed protocolIn each of the three studies used in this work, no statistical tests were performed on the raw data (i.e, observations), whether for viral load or for immune marker measurements, to identify statistical differences between treatment groups, as the statistical analyses were already been performed in the respective papers. Statistical significance of the effect of groups in model estimation is indicated in the tables by stars: *, p<0.05; **, p<0.01; ***, p<0.001 and were estimated by Wald tests (Monolix software version 2019R1).
Model parameters were estimated with the SAEM algorithm (Monolix software version 2019R1). Graphics were generated using R version 3.6.1 and Excel 2016 and details on the statistical analysis for the experiments can be found in the accompanying figure legends. Horizontal red dashed lines on graphs indicate assay limit of detection.
Appendix 1
Model building
In the model presented in the manuscript, we considered the two compartments of the URT, trachea, and nasopharynx, as two distinct compartments (i.e., without transfer of virus between them), as described by Equation AE1. In each of them, the viral dynamics are described by a targetcell limited model augmented with a compartment describing the dynamics of the inoculated virus (V_{s}). Moreover, in the statistical model describing the model parameters, the three parameters β, δ, and P were assumed as jointly estimated between the two compartments, with shared random effects and covariates and considering that parameters β and δ are equal in both trachea and nasopharynx (β^{T} = β^{N}, δ^{T} = δ^{N}).
Initially, random effects were added on the three parameters. However, taken into consideration identifiability issues that are usually encountered between the viral infectivity (β) and the viral production (P), we decided to remove the possibility of interindividual variability on the parameter P. This choice was also driven by multiple model estimations showing less robust estimations when variability was allowed in both parameters β and P. In particular, the estimate of the viral production was impacted by a ratio between the parameter and its standard error (RSE) higher than 100%.
Comparison of the parameters between the tracheal and the nasopharyngeal compartments
To decide which of these three parameters were assumed to be equal between the two compartments, all possibilities were tested and compared, using the BICc as selection criteria. As shown in Appendix 1—table 1, we started with the model in which all parameters were equal between the two compartments and we progressively relaxed this hypothesis. During this step, no exchange of virions between the two compartments of the URT was possible (g=0). Once all models estimated, we kept the one with the lowest value of BICc, meaning with the highest negative difference of BICc compared to the initial model. We identified the model with only the viral production varying between the two compartments as the best one to fit the data.
Identification of group effects
Once the structure of the statistical model defined, we tried to identify on which parameters an effect of the group of treatment could be identified and by extension on which biological mechanisms. In this step, we were interested in four parameters: β, δ, P, and c, the latter being the clearance of de novo produced virions. In the study, three groups of treatments were considered as constant categorical covariates: naïve, convalescent, and convalescent vaccinated. We performed a forward selection approach using the BICc as selection criteria to find the best model, using the model without covariate as initial model. At each step the model decreasing the most the value of the BICc is selected and the procedure stops once the BICc does not decrease anymore. At each step of the procedure, the statistical significance of covariate added into the model was verified via a Wald test. As shown in Appendix 1—table 2, the selected model identified a group effect on the viral infectivity and the loss rate of infected cells.
Based on all these results, the optimal statistical model with adjustment for groups of treatment was defined as follows:
Exchange of viruses between the nasopharyngeal and tracheal compartments
Afterward, we tested the possibility of an exchange of free plasma virus from between the two compartments of the URT. We made the hypothesis of a constant firstorder exchange and we tested the addition a transfer of virions from nasopharyngeal to tracheal compartments and vice versa, with a migration rate g_{NT} and g_{TN}, respectively. To this end, equations of infectious (V_{i}) and noninfectious (V_{ni}) viruses in Equation AE1 between the two compartments were linked as follows:
with the arrow symbolizing the modification of the equations defined in Equation AE1 and g_{NT} and g_{TN} being two positive rates. As a first step, we tried to estimate either bidirectional or one of the two unidirectional transfers using the data from the 18 NHPs of the first study described in the main paper. However, data were too spare to bring enough information to get estimations. Consequently, as a second step, additional data were used: two naïve macaques were exposed to the same dose (1 × 10^{6} pfu) of SARSCoV2 than the 18 NHPs of the main study. However, instead of being inoculated via intratracheal (4.5 mL) and intranasal (0.5 mL) routes, these latter received inoculum via intragastric (4.5 mL) and intranasal (0.5 mL) routes. Similar to the main study, the viral gRNA dynamics in both tracheal and nasopharyngeal compartments were repeatedly measured during the 20 days following the challenge (Figure 2—figure supplement 2E).
These two additional macaques having not received intratracheal inoculum, viral dynamics measured in this same compartment was expected to come from (at least partially) an exchange with the nasopharynx and thus bring information about it. However, having only two macaques without virions inoculated via intratracheal route, no enough information were available to totally estimate the model with exchanges. Consequently, these two additional NHPs having similar characteristics than the 18 NHPs involved in the main study, we made the assumption that the viral dynamics in nasopharynx after inoculation and the viral dynamics in the trachea, once the transfer initiated, should be described by the same model (without inoculum in trachea) and those by the same parameters. We expected that the difference of dynamics in trachea between these two set of macaques could allow an estimation of the parameters g_{TN} and/or g_{NT}. For that reason, we estimated the model in Equation AE1 using data from the 18 NHPs of the main study. Then using the data from the two additional NHPs, and assuming all parameters of the model resulting from Equation AE2 as fixed (see Appendix 1—table 2), except g_{TN} and g_{NT}, we tried to quantify the transfers of virions.
The estimation of multiple models on those two animals tended to conclude that only a unidirectional transfer of viruses from the nasopharyngeal to the tracheal compartment should be explored, with an estimation of g_{NT} ranging from 0.9 to 2.5 day^{–1}. Once these values quantified, we tried to update/reestimate the model, initially estimated on the 18 NHPs, using only a unidirectional transfer from nasopharynx to trachea and fixing the value of the migration rate at the different values aforementioned. However, all tested values of g_{NT} led irremediably to a degradation of the model with an increase of at least two points of BICc.
An estimation of the parameter g_{NT} by profile likelihood (results not shown) led to a strictly increasing profile of the likelihood (the lower the better) and was thus no more conclusive. Consequently, no exchange of virions were assumed in the final model and the parameters g_{NT} and g_{TN} were fixed at 0 day^{–1}.
Appendix 2
BICc as selection criteria and multiple testing adjustment
In the case of classic covariate selection approaches using pvalues as selection criteria, particular attention must be paid to take into account the dependence of the results on the number tests performed.
Over the years, multiple corrections have been proposed to adjust results for test multiplicity (e.g., Bonferroni correction, Benjamini and Hochberg correction among others).
Although we verified the significance of the covariate selected in our model, our covariate selection approach relies on the BICc. To ensure the robustness of the BICc as selection criterion despite the multiplicity of the tests, we performed an additional simulation work.
We simulated M=25 longitudinal variables for 18 individuals and with similar time points than those found on our data, meaning at days 0, 4, 9, and 20 postinfection. Variables were simulated as whitenoise random variables such that for the ith subject at the jth time point, the mth variable was defined as ${X}_{ij}^{m}\sim \mathcal{N}\left(0,{\sigma}^{2}\right)$ , with m=1, …, M. In our simulations, we tested five values for the variance σ^{2} ranging from 1% to 10% (five variables simulated for each value of σ).
Assuming these variables as our timevarying covariates, we applied the forward selection approach used in our method by testing each of them in a univariate manner of both β and δ.
As shown in Appendix 2—figure 1, the 50 models built to evaluate the adjustment of either β or δ for the simulated variables provide similar results in terms of BICc, and thus whatever the value of the standard deviation σ used. Consequently, these results appear as quite robust to the multiplicity of the test. Moreover, as expected, adjustments for whitenoise random variables depict the degradation of the model in comparison to the model without covariates.
Evaluation of the robustness of the estimation
To evaluate the robustness of the parameter estimates obtained on our models, despite the small number of independent observations, we performed a bootstrap procedure with replacement (Thai et al., 2014), for B=50 iterations. The bootstrap parameter estimate was calculated as the median of the parameter estimates from the B bootstrap samples while the standard error of each parameter was calculated according to the definition of Thai et al., 2014, which means with the SE of the lth component of the vector of parameters given by:
with ${\widehat{\theta}}_{b}^{*\left(l\right)}$ being its estimate obtained at the bth iteration of the bootstrap and ${\widehat{\theta}}_{B}^{\left(l\right)}$ the bootstrap parameter estimate. For each bootstrap sample, we paid attention to keep the 1:1:1 ratio between the three groups of treatment, with six animals selected within each group. Results are reported in Appendix 2—table 1 and Appendix 2—table 2.
Data availability
No unique reagents were generated for this study. Data that support the findings of this study are provided in the source data files of this paper and gather data from (1) the study [Marlin, Nature Com 2021] used in this analysis, which are also directly available online in the section Source data of this related paper (https://www.nature.com/articles/s41467021253820#Sec17); (2) the study [Brouwer, Cell 2021] used in this analysis, which are also available from the corresponding authors of the related paper and (3) the study [Corbett, NEJM 2020] used in this analysis, which are also available online in the section Supplementary Material of the related paper, excel file labelled ("Supplementary Appendix 2"). Data from the main study [Marlin, Nature Com 2021] can also be found in the openaccess repository Dryad using the following DOI: https://doi.org/10.5061/dryad.1zcrjdfv7. The original code (mlxtran models and R) as well as model definition files including the full list of parameters used are available and freeofcost on github (Inria SISTM Team) at the following link: https://github.com/sistm/SARSCoV2modelingNHP, (copy archived at swh:1:rev:a704c80daebc949434694d3f4441e48293c461cc).

Dryad Digital RepositoryViral loads and antibody, cytokine and Tcell responses in NHPs following vaccination targeting SARSCoV2 RBD domain to cells expressing CD40.https://doi.org/10.5061/dryad.1zcrjdfv7
References

Kinetics of influenza A virus infection in humansJournal of Virology 80:7590–7599.https://doi.org/10.1128/JVI.0162305

A model for establishment, maintenance and reactivation of the immune response after vaccination against Ebola virusJournal of Theoretical Biology 495:110254.https://doi.org/10.1016/j.jtbi.2020.110254

BookMeasurement Error in Nonlinear ModelsChapman and Hall/CRC.https://doi.org/10.1201/9781420010138

TLR3 agonist and CD40targeting vaccination induces immune responses and reduces HIV1 reservoirsThe Journal of Clinical Investigation 128:4387–4396.https://doi.org/10.1172/JCI99005

Variable selection strategies and its importance in clinical prediction modellingFamily Medicine and Community Health 8:e000262.https://doi.org/10.1136/fmch2019000262

Concerns about SARSCoV2 evolution should not hold back efforts to expand vaccinationNature Reviews. Immunology 21:330–335.https://doi.org/10.1038/s41577021005449

Evaluation of the mRNA1273 Vaccine against SARSCoV2 in Nonhuman PrimatesThe New England Journal of Medicine 383:1544–1555.https://doi.org/10.1056/NEJMoa2024671

Approaches and Challenges in SARSCoV2 Vaccine DevelopmentCell Host & Microbe 28:364–370.https://doi.org/10.1016/j.chom.2020.08.002

How to test SARSCoV2 vaccines ethically even after one is availableClinical Infectious Diseases: An Official Publication of the Infectious Diseases Society of America 73:ciab182.https://doi.org/10.1093/cid/ciab182

SARSCoV2 viral dynamics in nonhuman primatesPLOS Computational Biology 17:e1008785.https://doi.org/10.1371/journal.pcbi.1008785

Withinhost mathematical models of hepatitis B virus infection: Past, present, and futureCurrent Opinion in Systems Biology 18:27–35.https://doi.org/10.1016/j.coisb.2019.10.003

Multiscale modelling reveals that early superspreader events are a likely contributor to novel variant predominanceJournal of the Royal Society, Interface 19:20210811.https://doi.org/10.1098/rsif.2021.0811

Thrombotic Thrombocytopenia after ChAdOx1 nCov19 VaccinationThe New England Journal of Medicine 384:2092–2101.https://doi.org/10.1056/NEJMoa2104840

Immunological surrogate endpoints of COVID2019 vaccines: the evidence we have versus the evidence we needSignal Transduction and Targeted Therapy 6:48.https://doi.org/10.1038/s4139202100481y

Evaluation of a novel multiplexed assay for determining IgG levels and functional activity to SARSCoV2Journal of Clinical Virology 130:104572.https://doi.org/10.1016/j.jcv.2020.104572

Viral dynamics of SARSCoV2 across a spectrum of disease severity in COVID19The Journal of Infection 81:318–356.https://doi.org/10.1016/j.jinf.2020.04.014

Neutralizing Response against Variants after SARSCoV2 Infection and One Dose of BNT162b2The New England Journal of Medicine 384:2453–2454.https://doi.org/10.1056/NEJMc2104036

Transmission of COVID19 in 282 clusters in Catalonia, Spain: a cohort studyThe Lancet. Infectious Diseases 21:629–636.https://doi.org/10.1016/S14733099(20)309853

ON IDENTIFIABILITY OF NONLINEAR ODE MODELS AND APPLICATIONS IN VIRAL DYNAMICSSIAM Review. Society for Industrial and Applied Mathematics 53:3–39.https://doi.org/10.1137/090757009

Quantifying dose, strain, and tissuespecific kinetics of parainfluenza virus infectionPLOS Computational Biology 17:e1009299.https://doi.org/10.1371/journal.pcbi.1009299

Nomenclature for immune correlates of protection after vaccinationClinical Infectious Diseases 54:1615–1617.https://doi.org/10.1093/cid/cis238

Complex correlates of protection after vaccinationClinical Infectious Diseases 56:1458–1465.https://doi.org/10.1093/cid/cit048

A guide to vaccinology: from basic principles to new developmentsNature Reviews. Immunology 21:83–100.https://doi.org/10.1038/s41577020004797

A contemporary view of coronavirus transcriptionJournal of Virology 81:20–29.https://doi.org/10.1128/JVI.0135806

Influenza Virus Infection Model With Density Dependence Supports Biphasic Viral DecayFrontiers in Microbiology 9:1554.https://doi.org/10.3389/fmicb.2018.01554

Evaluation of bootstrap methods for estimating uncertainty of parameters in nonlinear mixedeffects models: a simulation study in population pharmacokineticsJournal of Pharmacokinetics and Pharmacodynamics 41:15–33.https://doi.org/10.1007/s109280139343z

Modeling the viral dynamics of SARSCoV2 infectionMathematical Biosciences 328:108438.https://doi.org/10.1016/j.mbs.2020.108438

BookMixed Effects Models for Complex DataChapman and Hall/CRC.https://doi.org/10.1201/9781420074086
Decision letter

Frederik GrawReviewing Editor; Heidelberg University, Germany

Miles P DavenportSenior Editor; University of New South Wales, Australia
Our editorial process produces two outputs: i) public reviews designed to be posted alongside the preprint for the benefit of readers; ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.
Decision letter after peer review:
Thank you for submitting your article "SARSCoV2 mechanistic correlates of protection: insight from modelling response to vaccines" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Miles Davenport as the Senior Editor. The reviewers have opted to remain anonymous.
The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission. The reviewers appreciated the developed approach and analyses, but had several concerns regarding the robustness of the methods and estimates, as well as the conclusions drawn from it. In particular, they identified the following main points that would need to be addressed in order to advance the study:
Essential revisions:
(1) A more detailed explanation of the used methods. This especially applies to the selection of biomarkers, as well as addressing the potential risk of overfitting in case of the mechanistic model.
(2) A revision of the title to specify that the analysis is limited to nonhuman primates in the light of the specific comments by the reviewers below.
(3) A more detailed discussion of the limitations of the study and reconsideration of the conclusions with regard to the interpretability of the nonhuman data for the human situation and possible caveats in the analysis. An inclusion of human data sets as suggested by reviewer 2 is not necessarily warranted. But a more detailed/extended discussion of the literature on correlates of protection should be included. Please also have a look at the detailed suggestions and links provided by reviewer 2.
We would welcome a substantially revised version of the manuscript addressing all of the issues mentioned by the reviewers. Please also see the detailed suggestions that are made by the reviewers below, and we hope you will find them helpful in this regard.
Reviewer #1 (Recommendations for the authors):
The model used seems adequate, but there are several modeling options that could be better justified. This would help the reader understand the modeling approach better. In some cases, this is simply a case of explaining earlier why you are making the choice that you make. For example, why do you model infection and noninfectious virus (Line 119)? I think this is to fit both virus and infected cells, but this is not clear when you present this choice, which then makes it a little unclear in line 124 when you say that sgRNA is proportional to infected cells. Again, why this choice? Wouldn't it also be reasonable to consider gRNA plus sgRNA to be proportional to infected cells? Another modeling choice that could be better explained is why you consider the inoculum virus separately (its own compartment and its own clearance rate)? In Line 401 you say, "to be able to distinguish", but why is that needed?
When you describe the model in the Results section, it would also be important to say something about the coupling of URT and trachea – this is only at the end of the model description in the Methods, but since you mention the two compartments in the results, the reader (who may not read the Methods first) may be confused.
Some more questions about the model. Is the clearance rate the same for infectious and noninfections virus? Line 470, you use profile likelihood for these parameters, what about the values of the other parameters when you are calculating these profiles? Line 511, this seems a little bit circular: you estimate all parameters with gs zero, then fix all parameters and try to estimate the two gs. But perhaps the other parameters would be different with g not zero?
The process for automatic selection of biomarkers also needs a little more clarification. It seems to be a generalization of stepwise analyses, but the issue of the combination of two lower ranked markers potentially being better than a higher ranked one is not mentioned or discussed. This is an issue because each parameter is tried individually. Also, you say (line 526) that it is possible to add timevarying covariates in this methodology, but you don't say anything about how that is done (it is not trivial), and in fact you almost don't say anything about results with timevarying covariates. Some of these aspects may also be relevant for discussion or limitation of the approach.
Line 238, even though you say that the mechanistic model was better, that is not what is presented in Table S1, where the group effect model is better. Please clarify. Also, in line 246 you state "no additional effect", but Table S1 presents results for δ. Also, Table S1 is not mentioned in this section of the manuscript.
I may have missed it, but it seems that in the "Results" nothing is said about time varying markers…
Something seems strange in some graphs in figure 2b and S8. In some cases, the mean value appears to be below every individual value; or sometimes above every individual value (for example, in the naïve). Also, are the thin lines individual predictions as stated, why do these lines have kinks instead of being smooth and do they really go through every data point in every animal (this is clear for example in the naïve, where the thin lines are clearly visible, are not smooth and go through every data point)?
Reviewer #2 (Recommendations for the authors):
– What the authors state the main contribution as (new framework of modeling) and what I found to be the main conclusion (vaccineinduced antibody binding indicates protectiveness) are different. Because the model and analytical approach is not particularly novel, I would suggest redesigning the manuscript to highlight the data and biological conclusions.
– Getting into the business of claiming there is some magical "correlate" is very tricky. I don't feel that this is what we need to understand the infection or response to vaccination, and the work in this manuscript does not assess sufficient data to make broad claims about any clinicallymeaningful correlate. With as much data as there is in the literature, a validation using human data is warranted.
– In the introduction, the authors reference studies that suggest "binding antibodies to SARSCoV2 and in vitro neutralization of virus infection are clearly associated with protection". This seems to be the same conclusion that the authors came to and makes me question: what is new here?
– A more robust introduction that highlights the current conclusions and limitations of models in the literature is needed, particularly because the model presented is only a minor modification of those models. In addition, the model lacks immune dynamics, which would make its usefulness limited.
– In addition to the point above, while the model is a standard, previously published viral dynamics model with minor modifications, it was not adequately described for the broad readership of eLife. Perhaps I missed it, but I couldn't find where the terms were statistically justified?
– Several models have shown that viral load data, even when defined between compartments, is insufficient to distinguish transport between the nasal passages and trachea (e.g., Khan et al. Viruses, Ke et al. MedRxiv, Pinky et al. PCB, among others). In light of these studies, this part of the model seems unnecessary and one unsupported by the data. In addition, these other studies have also shown that viral and immune parameters are distinct in the nasopharynx and trachea. In the authors' work, the data were joined and no effective comparison was made, which may cloud the parameter estimates and conclusions. Minimally, a discussion and reference of these published works is needed.
– Another limitation of the model and the data is that viral loads in the nasopharynx can be similar when disease is not. Can the approach here be used to assess vaccine efficacy (in terms of reducing disease)?
– Curiously, CD8 T cells were measured in the serum, but this type of measurement tends to not reflect the events in the tissue where resident T cells are thought to be important. There were some comments made about the infected cell clearance rate, but it did not seem that the data were used to evaluate the parameter estimates or model conclusions.
– Macaques, unfortunately, are not a great model to assess the individual heterogeneity observed in humans. Dissecting the heterogeneity from different sources would be greatly beneficial, but it seems naïve to assume that this can be done using animal models. I suggest using a human dataset, where there are many to choose from in the literature, or artificially creating one from human data and see if the model can still perform.
Possible data sets and information:
(i) CDC list of studies of correlates of protection, reinfection studies, etc.:
https://www.cdc.gov/coronavirus/2019ncov/science/sciencebriefs/vaccineinducedimmunity.html
(ii) These quantify antibody binding + other immunology: https://www.science.org/doi/full/10.1126/science.abm3425
https://www.nature.com/articles/s41586021037382?r=artikellink
https://www.sciencedirect.com/science/article/pii/S0092867421007066
(iii) One that has nice temporal data of lots of different antibodies and cells: https://www.science.org/doi/10.1126/science.abf4063
(iv) One that shows antibody levels + whether the macaque was protected: https://www.nature.com/articles/s41586020030416
– In reference to Table S1, it is stated that vaccinated and unvaccinated animals could be distinguished. Why is this the goal? We would know who is/isn't vaccinated – a better question is who among the vaccinated would not be protected from infection and/or disease.
– In Line 189191, the authors state "both specific antibodies and specific CD8^{+} T cells are mechanisms commonly considered important for killing infected cells. We retained the antiRBD binding IgG Ab that were positively associated to the increase of the loss of infected cells.". Unfortunately, this statement is incorrect as antibodies cannot kill infected cells. Antibodies neutralize virus so that cells are not infected, but they cannot kill cells.
– Prior studies have found that antibodies have a limited role during infection (e.g., see Goyal et al. Viruses). How does this fit into the results?
– Figure 2 is difficult to assess what is data and what is the model, and the model fit is difficult to see. It seems a though the 95% confidence intervals (gray shading) are also not indicative the model. Why are there multiple peaks and does the model capture that?
– The discussion could be improved. It lacked discussion of how the work fits into the literature, particularly other published SARSCoV2 models, how the correlate might be used in the clinic, etc. In addition, a discussion about the differences in dynamics within the nasopharynx and trachea in vaccinated individuals would have been interesting.
– It's unclear to me why the reproduction number R was under 1 for the alphaCD40.RBD group when there was replication for several days in this group.
Reviewer #3 (Recommendations for the authors):
First, a suggestion regarding my overfitting concern:
• The number of nonhuman primates in each arm of these studies is understandably limited, but crossvalidation could be used to reduce the extent of overfitting even with limited replicates.
Next, two suggestions to enhance clarity:
• The rationale for the use of profile likelihood should be made explicit. Why do a few parameters using this approach and then later do rest with full maximum likelihood? Is it a limitation of the optimization software?
• In the model, why separately track virions from the inoculum versus virions generated in host? If it is just to allow for different clearance rates, why do we expect these to be different?
https://doi.org/10.7554/eLife.75427.sa1Author response
Essential revisions:
1.) A more detailed explanation of the used methods. This especially applies to the selection of biomarkers, as well as addressing the potential risk of overfitting in case of the mechanistic model.
As pointed out by the reviewers, the first version of the article was somewhat poor in explanations about the methods used. Accordingly, we have described the modeling by adding some more precision. In particular, we have added more explanations about:
1. The choice of the mechanistic model used and more specifically about
a. The distinction of infectious and noninfectious viruses (see Reviewer 1, remarks 1a, 3a)
b. The distinction of the inoculated and the produced de novo virions (see Reviewer 1, remark 1b)
c. About our modeling choice to integrate subgenomic data in the observation model (see Reviewer 1, remark 1a)
d. About the coupling of URT compartments (the trachea and nasopharynx) and the exchange of virions between these two compartments (see Reviewer 1, remarks 2 and 3c and reviewer 2, remark 6).
2. The algorithm of covariate selection that we implemented, and more specifically about
a. The type of algorithm used and its robustness (see Reviewer 1 remark 4a, Reviewer 3 remark 1)
b. The difference of the algorithm compared to classical algorithms to incorporate timevarying variables (see Reviewer 1, remark 4b)
c. The limitation of the proposed algorithm (see Reviewer 1, remark 4b)
2.) A revision of the title to specify that the analysis is limited to nonhuman primates in the light of the specific comments by the reviewers below.
As requested, we modified the title of the article to specify its current limitation to nonhuman primates. The proposed new title is “Modelling the response to vaccine in NonHuman Primates to define SARSCoV2 mechanistic correlates of protection”.
3.) A more detailed discussion of the limitations of the study and reconsideration of the conclusions with regard to the interpretability of the nonhuman data for the human situation and possible caveats in the analysis. An inclusion of human data sets as suggested by reviewer 2 is not necessarily warranted. But a more detailed/extended discussion of the literature on correlates of protection should be included. Please also have a look at the detailed suggestions and links provided by reviewer 2.
No application on human data has been included in this revised version of the article. However, a particular attention has been paid to information provided by reviewer 2 to enrich the discussion about the proposed approach, whether about its limitations or its potential extensions.
Reviewer #1 (Recommendations for the authors):
The model used seems adequate, but there are several modeling options that could be better justified. This would help the reader understand the modeling approach better. In some cases, this is simply a case of explaining earlier why you are making the choice that you make. For example, why do you model infection and noninfectious virus (Line 119)? I think this is to fit both virus and infected cells, but this is not clear when you present this choice, which then makes it a little unclear in line 124 when you say that sgRNA is proportional to infected cells. Again, why this choice? Wouldn't it also be reasonable to consider gRNA plus sgRNA to be proportional to infected cells?
About infectious and noninfectious virus:
The SARSCoV2 virus exists in both infectious and noninfectious form (1,2) with a proportion μ of virus able to infect new target cells relatively low compared to noninfectious one. According to (3), this ratio can be quantified with a lower bound ranging from 10^{5} and 10^{4}. Expert opinion as well as the evaluation of this same ratio on similar NHP than those used on our main study (unpublished results) led us to the value of 10^{3}. The final choice of the value of the parameter μ was guided by the comparison of the results obtained for μ=10^{3} and 10^{4}. Indeed, estimations appeared as more stable with the highest value.
Nevertheless, we agree with the fact that outside the ODE system, compartments of virions are always used as the sum of V_{i} and V_{ni}. Consequently, the model could be rewritten assuming a single compartment for plasma free viruses (V), paying attention to replace the terms βV_{i}T by βμVT to keep the ratio between infectious and noninfectious viral particles. We decided to write the model with distinct compartments between V_{i} and V_{ni} for a better visual understanding of the model and of the parameter μ.
We add the justification about the modelling of both infectious and noninfectious viruses in the manuscript (see Results page 7, Lines 121122 and Methods, page 21, Lines 453457).
[Addition applied in section Results, page 7, Lines 121122]
“Although a single compartment for de novo produced viruses (V) could be mathematically considered, two distinct ODE compartments were assumed for a better understanding of the model.”
[Addition applied in section Methods, page 21, Lines 453457]
“… remaining proportion of virions is noninfectious (V_{ni}). Mathematically, a single compartment (V) for de novo produced virions could be considered in the model, with μV and (1 μ)V representing the respective contributions of infectious and noninfectious viruses to the biological mechanisms. However, to have a better visual understanding of the distinction between the two types of viruses, we wrote the model with distinct compartments, V_{i} and V_{ni}.”
About subgenomic RNA:
In our model, we assume the subgenomic RNA (sgRNA) as proportional to infected cells. First, sgRNA is only transcribed in infected cells (4). Second, the use of this additional marker allowed us to improve the identifiability of model parameters. According to the paper of Miao et al. (2011) (5), the use of observations for both virus and infected cell compartments allows to counteract structural identifiability issue in this type of targetcell limited models, in particular between the parameters of infectivity (β) and viral production (P). In accordance with this paper, we tested different parametrizations to jointly describe gRNA and sgRNA in the observation model. In particular, we compared the use of the two following observation models:
(a) with gRNA defined by (V_{i} + V_{ni} + V_{s}) as total viral load and sgRNA, resulting from viral replication, as proportional to (V_{i} + V_{ni}),
(b) with gRNA defined by (V_{i} + V_{ni} + V_{s}) and sgRNA as proportional to (I_{1} + I_{2}).
As expected, the observation model (a) was impacted by identifiability issues compelling us to fix either β or P, while the estimation of these two parameters was possible with the model (b). Moreover, this latter led to better values of BIC (25 points of BICc).
As proposed by the reviewer, another possibility could have been to define the sum of gRNA and sgRNA as proportional to infected cells assuming gRNA characterized by (V_{i} + V_{ni} + V_{s}) and gRNA + sgRNA by (I_{1} + I_{2}). But we did not try this solution.
To better explain this choice of modelling, we modified the manuscript as follows (section Results, page 7, Lines 133135 and section Methods, page 2425, Lines 508521).
[Addition applied in section Results, page 7, Lines 133135]
“We assumed that gRNA and sgRNA were proportional to the free virus and the infected cells, respectively. This modeling choice relied on both biological and mathematical reasons (see the section methods for more details).”
[Addition applied in section Methods, page 2425, Lines 508521]
“For the observation model, we jointly described genomic and subgenomic viral loads in the two compartments of the URT. We defined genomic viral load, which characterizes the total viral load observed in a compartment (nasopharynx or trachea), as the sum of inoculated virions (V_{s}), infectious (V_{i}) and noninfectious virions (V_{ni}). The sgRNA was described as proportional to the infected cells (I_{1} + I_{2}). This choice was driven by two main reasons. First, sgRNA is only transcribed in infected cells (76). Second, as described by Miao et al. (2011) (77), to overcome identifiability issues between the parameters β and P typically observed in targetcell limited models. The comparison of the two observation models describing sgRNA as either proportional to virions produced de novo (V_{i} + V_{ni}) or proportional to infected cells (I_{1} + I_{2}) confirmed this conclusion. In addition to a better BICc value (25 points) compared with the first model, the second one allowed the estimation of both β and P by counteracting the identifiability problems faced with the first model (results not shown). Accordingly, the log_{10}transformed gRNA and sgRNA of the ith animal at the jth time point in compartment X (nasopharynx or trachea), denoted gRNA_{ij}^{X} and sgRNA_{ij}^{X} respectively, were jointly described by the following equations …”
Another modeling choice that could be better explained is why you consider the inoculum virus separately (its own compartment and its own clearance rate)? In Line 401 you say, "to be able to distinguish", but why is that needed?
In our model, we distinguished between the inoculated virus and the virus produced de novo by the infected cells by including a compartment V_{s}. This choice was guided by requests of immunologists for a better understanding of the dynamics of the inoculum. Moreover, this allowed us to use all the information provided by preclinical studies and avoided to estimate or randomly fixe the initial conditions for the V_{i} and V_{ni} compartments, as usually done. Finally, this modeling choice allowed to differentiate the elimination rate of the two virus types. In particular, the inoculum corresponds to a high number of locally inoculated virions, a large fraction of which can be eliminated faster than de novo produced virions due to the experimental conditions (swallowing and natural downstream influx). As described in the paper, using profile likelihood, we showed an improvement in the estimation of model (gain of 5 points of likelihood) when distinct values for the clearance of inoculated and de novo produced virions were considered, leading to a faster elimination of the inoculum, as expected.
To better explain the introduction of the inoculum as ODE compartment, the manuscript was modified as follows (section Methods, page 20, Lines 440448 and pages 21, Lines 458465).
[Addition applied in section Methods, page 20, Line 440448]
“We completed the model by adding a compartment for the inoculum that distinguishes the injected virus (V_{s}) from the virus produced de novo (V_{i} and V_{ni}). To our knowledge, this distinction has not been proposed in any previous work. Two main reasons led us to make this choice. First, it allowed us to study the dynamics of the inoculum, in particular during the early phase of viral RNA load dynamics. Second, as described in more detail below, it gave us the opportunity to use all the information provided by the preclinical studies, such as the known number of inoculated virions, to define the initial conditions of the ODE model rather than estimating or randomly fixing them for V_{i} and V_{ni}, as is usually done.”
[Addition applied in section Methods, pages 21, Line 458465]
“Finally, virions produced de novo and those from the inoculum are cleared at rates c and c_{i}, respectively. Distinct clearances were considered to account for the effects of experimental conditions on viral dynamics. In particular, it is hypothesized that, animals being locally infected with large numbers of virions, a large proportion of it is assumed to be rapidly eliminated by swallowing and natural downstream influx, in contrast to the de novo produced virions. However, it is important to keep in mind that this distinction was possible because of the controlled experimental conditions performed in animals, (i.e., exact timing and amount of inoculated virus known, and frequent monitoring during the early phase of the viral dynamics).”
When you describe the model in the Results section, it would also be important to say something about the coupling of URT and trachea – this is only at the end of the model description in the Methods, but since you mention the two compartments in the results, the reader (who may not read the Methods first) may be confused.
We have reformulated the presentation of the two compartments (trachea and nasopharynx) in the Results section to better insist on the coupled modelling earlier in the paper. Moreover, we have considered the comment of Reviewer 2 (see comment 6 below) to justify the absence of viral exchange between the two compartments. The manuscript was modified accordingly (section Results, page 6, Lines 107109 and page 7, Lines 124133, and section Methods page 28, Lines 591596).
[Addition/modification applied in section Results, page 6, Lines 107109]
“The viral dynamics during primary infections were characterized by a peak in genomic RNA (gRNA) production three days postinfection in both tracheal and nasopharyngeal compartments, followed by …”
[Addition/modification applied in section Results, page 7, Lines 124133]
“In both compartments of the upper respiratory tract (URT), the trachea and nasopharynx, viral dynamics were distinctively described by this model (Figure 2A). Viral exchange between the two compartments was tested (either from the nasopharynx to the trachea or vice versa). However, as described in the literature (28, 42, 43) and demonstrated by the additional modeling work in Appendix 1 “Model building”, viral transport within the respiratory tract plays a negligible role in viral kinetics compared with viral clearance. Consequently, no exchange was considered in the model. Using the gRNA and sgRNA viral loads, we jointly estimated (i.e., shared random effects and covariates) the viral infectivity (β), the viral production (P), and the loss rate of infected cells (δ) in the two compartments.”
[Addition applied in section Methods, page 28, Lines 591596]
“The possibility of migration of free plasma virus between the nasopharynx and the trachea was tested. However, as widely described in the literature, the transport of viral particles within the respiratory tract is negligible in the viral dynamics and is difficult to estimate. The reader can refer to the Appendix 1 “Model building” for an additional modelling work conducted to estimate this exchange and provided the same conclusion. Accordingly, the two compartments of the URT were assumed as distinct in our model.”
Some more questions about the model. Is the clearance rate the same for infectious and noninfections virus?
The clearance rate for infectious and noninfectious viruses is the same, in both URT compartments (c=20 day^{1}). Due to identifiability issues, we were not able to identify distinct clearances, although this could be biologically relevant under certain conditions. Because this point was unclear to the reviewer, we modified the manuscript to clearly mention it in the paper (section Results, page 7, Lines 135139 and section Methods page 22, Line 468).
[Addition applied in section Results, page 7, Lines 135139]
“Due to identifiability issues, the duration of the eclipse phase (1/k), the clearance of free viruses from the inoculum (c_{i}) and produced de novo (c) were estimated separately by profile likelihood and assumed to be identical in the two compartments of the URT. In addition, infectious and noninfectious viruses were assumed to be cleared at the same rate.”
[Addition/modification applied in section Methods, page 21, Line 466]
“Because of identifiability issues, similar clearances for infectious and noninfectious viruses were used. Accordingly, …”
Line 470, you use profile likelihood for these parameters, what about the values of the other parameters when you are calculating these profiles?
Because of identifiability issues, we performed profile likelihood to determine appropriate values for the three parameters k, c, and c_{i}. Profile likelihood consists in sequentially testing different combinations of values of the parameters of interest by fixing them in the model and optimizing their values by choosing the combination that maximizes the likelihood. Accordingly, for each tested combination, all unfixed parameters were estimated here by the SAEM algorithm. Since we are only interested in the maximum likelihood estimates, we did not look at the value of the estimated parameters for each combination, but only for the one that yields the highest loglikelihood (see section Methods, page 26 lines 554559).
[Addition applied in section Methods, page 26, Lines 554559]
“The profile likelihood consists in defining a grid of values for the parameters to be evaluated and sequentially fixing these parameters to one of these combinations of values. The model and all the parameters that are not fixed are then estimated by maximizing the loglikelihood. In this process, all parameters that are assumed to be fixed in the model (i.e. μ and the initial conditions) are held fixed. Finally, the optimal set of parameters is chosen as the one optimizing the loglikelihood. Although the available data ….”
Line 511, this seems a little bit circular: you estimate all parameters with gs zero, then fix all parameters and try to estimate the two gs. But perhaps the other parameters would be different with g not zero?
We understand the reviewer’s concerns on this point. Nevertheless, the estimation was not completely circular as different data were used to estimate the parameters of the model in the absence of virus exchange and to estimate g. While the 18 NHPs from the main study were used in the first case, similar data obtained for two additional naive NHPs that had not received intratracheal inoculum were used to quantify exchanges. Accordingly, we assumed that their viral loads observed in the trachea should come from an exchange with the nasopharynx and thus should provide information to estimate g. However, because of the small number of new NHPs, we did not have enough observations to fully reestimate the model with exchange. Consequently, hypothesizing that the two naïve NHPs had similar physiological characteristics to the 18 others, we fixed the model parameters to the value found without exchange and we estimated g only. In a second time, we tried to estimate the global model directly with the 20 NHPs but results were not more conclusive, with g tending to zero.
In response to this comment, to the previous comment 2 and to the comment 6 of the Reviewer 2, we transferred the subsection “Methods/Exchange of viruses between nasopharynx and trachea compartments” in Appendix 1 “Model building” and we mentioned this file in section Methods (page 28, Lines 594).
The process for automatic selection of biomarkers also needs a little more clarification. It seems to be a generalization of stepwise analyses, but the issue of the combination of two lower ranked markers potentially being better than a higher ranked one is not mentioned or discussed. This is an issue because each parameter is tried individually.
The algorithm implemented for biomarker selection is not exactly a stepwise algorithm as it involves only a forward selection procedure. As mentioned by the reviewer, the case in which two lower ranked markers are better than one higher ranked marker was not mentioned in the manuscript. The use of stepwise selection (i.e. a combination of forward selection and backward elimination) could reduce the number of these cases. Although we could theoretically have implemented stepwise selection, we did not integrate backward elimination in the present work for two practical reasons. The first one is the computational time required to run the algorithm, since a large number of timevarying markers (more than 50) were tested univariately on all parameters of interest. Second, the proposed method for selecting covariates relies not only on a single selection criterion (here, the best fit criterion), as is usually the case, but also on the other two other criteria that define a mCoP (better capture of interindividual variability and replace the group of intervention). In our specific applications, using BIC as selection criterion (rather than pvalues) and the selection stopping after two steps, backward elimination would have been useless. Nevertheless, we manually checked the statistical significance of all added covariates at each step. Finally, it is important to note that we are not interested in finding “the” CoP (it is not a biological paper), but we want to find a robust CoP with the good statistical properties.
For sake of simplicity, we merged the modifications related to this comment with those related to the next comment. Both comments relate to the algorithm for selecting covariates.
Also, you say (line 526) that it is possible to add timevarying covariates in this methodology, but you don't say anything about how that is done (it is not trivial), and in fact you almost don't say anything about results with timevarying covariates. Some of these aspects may also be relevant for discussion or limitation of the approach.
We agree that the explanations about timevarying variables and the application of the covariate selection method to them were not entirely clear and sufficiently emphasized in the manuscript. We modified the paper to take these points more into account. In particular, we mentioned:
a) the lack of algorithms for selecting timevarying covariates already implemented in software such as Monolix (see section Methods page 29, lines 602611).
b) the use of linear interpolation to incorporate the dynamics of the markers as a timevarying covariate in the statistical model (see section Methods, pages 3031, lines 638672).
We also added a table in the supplementary files with the parameter estimates of the model in which viral infectivity was adjusted for the timevarying dynamics of antibodies inhibiting the attachment of ACE2 and RBD (see Table 2 of this document, referred as Supplementary file 4 in the manuscript). Finally, we extended the discussion to mention the limitations of the proposed algorithm (see section Discussion, page 16, lines 344349).
[Addition/modification applied in section Methods, page 29, Lines 602611]
“… capture this effect. Nowadays, many methods for selecting constant covariates already exist (80) and are implemented in software such as Monolix. However, these latter do not allow timevarying covariates. In this section, we present the algorithm we implemented to select timevarying covariates. We proposed a classical stepwise datadriven automatic covariate modeling method (Figure S10). However, initially implemented to select covariates from more than 50 biomarkers, computational time restricted us to consider only a forward selection procedure. Nevertheless, the method can be easily extended to classical stepwise selection in which both forward selection and backward elimination are performed sequentially. In particular, this extension could contribute to reduce the occurrence of the case in which the combination of two lower ranked markers might be better than a single higher ranked marker. Although the method was developed for timevarying covariates, it can also be applied to constant covariates.”
[Addition applied in section Methods, pages 3031, Lines 638672]
“Modelling hypothesis for timedependent covariates in our application
Using a populationbased approach to estimate our mechanistic model and similar to the adjustment of the model for constant covariates (e.g., groups of intervention), timevarying covariates are incorporated into the statistical model as individualspecific explanatory variables in the mixedeffects models. The identification of antibodies inhibiting the attachment of the RBD domain to the ACE2 receptor (ECLRBD) as the first timevarying CoP led to the definition of the timevarying viral infectivity for the ith animal as described in Equation (5), while the selection antiRBD IgGbinding antibodies (IggRBD) led to the elimination rate of infected cells given in Equation (6).
[Addition/modification applied in section Discussion, pages 16, Lines 344349]
…that best fit the data. Nevertheless, many modeling choices for the statistical model were made in this approach and more theoretical work evaluating the robustness of the results in their regards may be relevant for future works. In particular, we could relax the constraint of linear interpolation of marker dynamics by using simple regression models, allowing in the same time the integration of error model to account for measurement error for timevarying covariates (6365). Moreover, by construction, ….”
Line 238, even though you say that the mechanistic model was better, that is not what is presented in Table S1, where the group effect model is better. Please clarify.
We thank the reviewer for this relevant comment. Indeed, some information is missing in the table to fully understand the results. For the first criterion, we have reported the BICc values for the model without covariates (i.e., the initial model in the algorithm), for the model adjusted for group effects on β and δ, except for the third study, and for the model adjusted only for the marker, i.e., without group effect on δ. We added the BICc value of the model with β adjusted for the marker and δ adjusted for the groups of treatment in Table S1 (see below; labelled now Supplementary file 1) to clarify the results. In addition, as described in the algorithm, and in response to comment (4a), the algorithm is based on three criteria, and criterion 1 (best fit) aims at selecting the marker with the optimal value compared with the model without covariates or adjusted for immune markers. The main purpose of the comparison with the model adjusted for groups of treatment is to quantify the gap between the two models. As this is not the case here, the selection of two markers might be necessary to perform better than group effects.
The groups of treatment essentially intervene for the other two criteria: replacement of the effect and explanation of interindividual variability.
Also, in line 246 you state "no additional effect", but Table S1 presents results for δ. Also, Table S1 is not mentioned in this section of the manuscript.
An effect of the group of treatment was well identified only for the parameter β for the third dataset. The information on the parameter δ provided in Table S1 (criterion 3) refers to the interindividual variability of the parameter, random effects being considered on both β and δ. We included this information in the Table (see below) to verify that adjusting the model for β does not severely degrade the estimation of the variance of random effects on δ.
I may have missed it, but it seems that in the "Results" nothing is said about time varying markers…
In the Results section, timevarying markers were indeed not mentioned and we thank the reviewer for this relevant comment. We have corrected the manuscript accordingly (see section Results page 10, lines 195200, and page 1011, lines 209214).
[Addition/modification applied in section Results, pages 10, Lines 195200]
“To this aim, we performed a systematic screening by adjusting the model for each marker, and we compared these new models with the model without covariates and with the model adjusted for the groups. In particular, our approach allowed us to benefit from all the information provided by the overall dynamics of the immunological markers after the exposure by integrating them as timevarying covariates (see supplemental information for a detailed description of the algorithm) …”
[Addition applied in section Results, pages 1011, Lines 209214]
“… Finally, we looked at the estimated viral infectivity according to the binding inhibition assay in each animal. A positive dependence was found between the viral infectivity and the RBDACE2 binding inhibition measure, linking an increase of 10^{3} AU of the marker, whether over time or between animals, with an increase of 1.8% (95CI% [1.2%; 2.3%]) of the viral infectivity (see Supplementary file 4). Accordingly, the values at the time of exposure were not overlapping at all, distinguishing clearly the vaccinated and unvaccinated animals (see Figure 4C).”
Something seems strange in some graphs in figure 2b and S8. In some cases, the mean value appears to be below every individual value; or sometimes above every individual value (for example, in the naïve). Also, are the thin lines individual predictions as stated, why do these lines have kinks instead of being smooth and do they really go through every data point in every animal (this is clear for example in the naïve, where the thin lines are clearly visible, are not smooth and go through every data point)?
We thank the reviewer for these helpful comments. Indeed, these figures were very confusing.
In these figures, it is important to note that observations were represented by thick dashed lines (mean value) and shaded areas (interindividual variability), while the predicted data were represented by thin solid lines. Individual predictions were given by thin solid lines. These latter were not smooth because only predictions at time of observations were reported. We have corrected the figures accordingly (see Figure 2b and Figure S8 in Supplementary Materials).
Reviewer #2 (Recommendations for the authors):
– What the authors state the main contribution as (new framework of modeling) and what I found to be the main conclusion (vaccineinduced antibody binding indicates protectiveness) are different. Because the model and analytical approach is not particularly novel, I would suggest redesigning the manuscript to highlight the data and biological conclusions.
We think that the use of such mechanistic modelling for defining a mechanistic correlate of protection constitutes a novelty because, to our knowledge, mechanistic CoP has been defined with counterfactual approaches only (6).
However, we do acknowledge that the important message here is the result in the context of COVID19, therefore we have better highlighted the biological conclusion and dampen the methodological aspect (see section Introduction page 5, lines 77, 8082 and section Discussion page 14, lines 289290).
[Modification applied in section Introduction, page 5, line 77 and lines 8082]
“Here, we propose to apply a modelbased approach on NHP studies to evaluate i) […]. First, we present a mechanistic approach based on ordinary differential equation (ODE) models reflecting the virushost interaction inspired from models proposed for SARSCoV2 infection (26–31) and other viruses (3235).”
[Modification applied in section Discussion, page 14, lines 289290]
“We explored the mechanistic effects of three SARSCoV2 vaccines and assessed the quality of markers as mechanistic CoP (mCoP).”
– Getting into the business of claiming there is some magical "correlate" is very tricky. I don't feel that this is what we need to understand the infection or response to vaccination, and the work in this manuscript does not assess sufficient data to make broad claims about any clinicallymeaningful correlate. With as much data as there is in the literature, a validation using human data is warranted.
We certainly started this work with the hypothesis that neutralization is a likely candidate. However, in this work we try to respond to the following unknown questions:
– How much the neutralization could explain the dynamics of viral load?
– Which neutralization assay best explain the dynamics of the viral load?
– Is there any additional marker that participates to the control of the viral load and how it works?
The validation with human data is not straightforward as the quantification of the effects obtained in this study is helped by the experimental conditions of the studies in NHP: the date of infection is known and the viral load is precisely tracked.
As recommended, we have modified the manuscript to focus more on the results in these NHP experiments than the framework for defining a mCoP (see section Conclusion, page 17, lines 358360)
[Addition applied in section Conclusion, page 17, lines 358360]
“In conclusion, the modelling of the response to two new promising SARSCoV2 vaccines in NHP revealed a combination of effects with a blockade of new cell infections and the destruction of infected cells.”
– In the introduction, the authors reference studies that suggest "binding antibodies to SARSCoV2 and in vitro neutralization of virus infection are clearly associated with protection". This seems to be the same conclusion that the authors came to and makes me question: what is new here?
Two original aspects should be seen here:
1. We provide an indication on the biological effect of markers suspected as correlates of protection
2. We provide a quantification of the effect of the vaccine and more precisely of the markers that play the role of surrogate markers
– A more robust introduction that highlights the current conclusions and limitations of models in the literature is needed, particularly because the model presented is only a minor modification of those models. In addition, the model lacks immune dynamics, which would make its usefulness limited.
We have revised the Discussion section to better underline the existing literature on such dynamical models. We have also included a point on the limitation due to absence of modelling of the dynamics of the immune response (see Introduction, page 5, lines 8082 and section Discussion, pages 1617, lines 342351 and 356358).
[Modification applied in section Introduction, page 5, lines 8082]
“First, we present a mechanistic approach based on ordinary differential equation (ODE) models reflecting the virushost interaction inspired from models proposed for SARSCov2 infection (26–31) and other viruses (3235).”
[Modification applied in section Discussion, pages 1617, lines 340349 and 354356]
“This result needs to be confirmed as it may also be the consequence of a different local immune response (60). The choice of the structural model defining the hostpathogen interaction is a fundamental step in the presented approach. Here, it was well guided by the biological knowledge, the existing models for viral dynamics (34, 61, 62) and the statistical inference allowing the selection of the model that best fit the data. Nevertheless, many modeling choices for the statistical model were made in this approach and more theoretical work evaluating the robustness of the results in their regards may be relevant for future works. In particular, we could relax the constraint of linear interpolation of marker dynamics by using simple regression models, allowing in the same time the integration of error model to account for measurement error for timevarying covariates (63–65). […]. Finally, it should be underlined that the dynamics of the immune response has not been modelled as suggested for instance for B cell response (68). This clearly constitutes the next step after the selection of the markers of interest as done in the present work.”
– In addition to the point above, while the model is a standard, previously published viral dynamics model with minor modifications, it was not adequately described for the broad readership of eLife. Perhaps I missed it, but I couldn't find where the terms were statistically justified?
Many corrections were made in the manuscript to better justify our modelling choices. The reviewer can refer to the comments (1a), (1b) and (3a) of the Reviewer 1.
We also added justifications about our modeling choices for the statistical model describing the three parameters β, δ and P:
(a) A logtransformation was considered for δ and P to ensure their positivity while a log10transformation was used for β to also tackle some convergence issues and enable a good estimation of the fisher matrix.
(b) The two URT compartments were jointly considered in the statistical model for the three parameters. Although these latter can take different values between the compartments, this modeling choice allowed for the sharing of random effects and covariates. This choice was driven by the hypothesis of homogenous interindividual variability within the URT. Moreover, this modeling of the parameters appeared us as more mechanistically correct rather than considering totally distinct parameters, and provided more power for their estimation.
To account for these different points, the manuscript was modified as follows (see section Methods page 2324, Lines 487506). Furthermore, to better explain our final choice for the statistical model we added an additional file, entitled “Model Building” describing the procedure applied.
Finally, we added a discussion in the manuscript about how these modeling choices could influence the selection of covariates (see section Discussion, page 16, Lines 349353).
[Addition applied in section Methods, page 2324, Lines 487506]
“A logtransformation was adopted for the parameters δ and P to ensure their positivity while a log_{10}transformation was chosen for viral infectivity to also improve the convergence of the estimation. […] Finally, the adjustment of the model for the categorical covariates of groups of treatment, natural infection and/or vaccination, identified β and δ as the parameters with a statistically significant effect of these covariates (i.e., φ^{P}_{conv} = 0 and φ^{P}_{CD40}=0).”
[Addition applied in section Discussion, page 16, Lines 349353]
“…measurement error for timevarying covariates. Moreover, by construction, we assumed similar interindividual variability and effects of covariates within the two URT compartments as well as similar values for the viral infectivity and the loss rate of infected cells. Viral load dynamics measured in lungs being different from those in the URT (66, 67), the relaxation of this hypothesis of homogeneous physiological behavior in the URT may be pertinent to extend the model to the LRT. Finally, it should be underlined that the dynamics of the immune response ….”
– Several models have shown that viral load data, even when defined between compartments, is insufficient to distinguish transport between the nasal passages and trachea (e.g., Khan et al. Viruses, Ke et al. MedRxiv, Pinky et al. PCB, among others). In light of these studies, this part of the model seems unnecessary and one unsupported by the data. In addition, these other studies have also shown that viral and immune parameters are distinct in the nasopharynx and trachea.
We thank the reviewer for these relevant comments and the papers. Considering this remark and the remark 2 of the Reviewer 1, we modified the Results section to integrate papers mentioned in this comment and we transferred the section in which we tried to quantify virus exchanges in Appendix 1 “Model building”. Moreover, we added in this file the procedure applied on our model to build our statistical model and identify similar viral infectivity and loss rate of infected between the trachea and the nasopharynx. Nevertheless, although parameters were jointly estimated between the two compartments, distinct values of viral production were found.
In the authors' work, the data were joined and no effective comparison was made, which may cloud the parameter estimates and conclusions. Minimally, a discussion and reference of these published works is needed.
The parameter estimates obtained in the first study with the model adjusted for the groups of intervention were compared with previously reported modelling results. This comparison was included in the Results section (see pages 78, Lines 139158).
[Modification / addition applied in section Results, pages 78, Lines 139158]
“We estimated the viral infectivity at 0.95x10^{6} (CI_{95%} [0.18x10^{6}; 4.94x10^{6}]) (copies/ml)^{1} day^{1} in naïve animals, which is in the range of previously reported modelling results whether in the case of SARSCoV2 virus (27, 29) or influenza (32, 33). […] In particular, they are in the range of estimates obtained within the URT, either in NHP (28) or in humans (29), with the product pxT_{0} equals to 15.1x10^{8} (CI_{95%} [3.98x10^{8}; 58.1x10^{8}]) and 0.21x10^{8} (CI_{95%} [0.088x10^{8}; 0.48x10^{8}]) virions/ml/day in the nasopharynx and the trachea, respectively.…”
– Another limitation of the model and the data is that viral loads in the nasopharynx can be similar when disease is not. Can the approach here be used to assess vaccine efficacy (in terms of reducing disease)?
The mathematical model can be extended to incorporate clinical progression as it has been proposed in other contexts (8). However, it was not possible to fit such model here because of the reduce sample size of the first trial and clinical data were not available in the two others.
– Curiously, CD8 T cells were measured in the serum, but this type of measurement tends to not reflect the events in the tissue where resident T cells are thought to be important. There were some comments made about the infected cell clearance rate, but it did not seem that the data were used to evaluate the parameter estimates or model conclusions.
In agreement to the reviewer comment, a great CD8 T cell response to the αCD40.RBD vaccine has been demonstrated in mice splenocytes but was not found in macaques’ blood (9). We systematically analyzed the available markers but we were not surprised that we could not estimate an association between the CD8 T cell response in the blood and the clearance of infected cells. This is what we have quickly discussed in the current paper.
– Macaques, unfortunately, are not a great model to assess the individual heterogeneity observed in humans. Dissecting the heterogeneity from different sources would be greatly beneficial, but it seems naïve to assume that this can be done using animal models. I suggest using a human dataset, where there are many to choose from in the literature, or artificially creating one from human data and see if the model can still perform.
It is true that results may differ from macaques and humans although the importance of neutralizing antibodies has been demonstrated in the two species. One advantage of the experimentations in macaques was the availability of repeated measures of viral load and the control of the date of infection which were very helpful for the estimation of model parameters. The translation of such model in humans is certainly an interesting perspective that we are now mentioning in the Discussion section (see page 17, lines 368370).
[Addition applied in section Discussion, page 17, Lines 368370]
“… waning protective antibody levels (71, 72), at least in NHP studies although the framework could be extended to human studies using mixed approaches of within and between hosts modelling (73) providing that enough information is collected.”
– In reference to Table S1, it is stated that vaccinated and unvaccinated animals could be distinguished. Why is this the goal? We would know who is/isn't vaccinated – a better question is who among the vaccinated would not be protected from infection and/or disease.
The objective was not to distinguish which animal was vaccinated or not but to use this information to identify which part of the model would be impacted and thus on which biological mechanisms vaccination and/or natural infection seem to play a role. In the Table S1, information about groups is used for the criterion 2 to verify that the immune marker selected by the selection method is able to capture information provided by the group: once adjusted for the maker, the covariate of group is no more significant.
– In Line 189191, the authors state "both specific antibodies and specific CD8^{+} T cells are mechanisms commonly considered important for killing infected cells. We retained the antiRBD binding IgG Ab that were positively associated to the increase of the loss of infected cells.". Unfortunately, this statement is incorrect as antibodies cannot kill infected cells. Antibodies neutralize virus so that cells are not infected, but they cannot kill cells.
The antibodies can prevent target cell infection and concur to the opsonization of viral particles but there are also playing an important role for the destruction of infected cells by several (Fcmediated) pathways including antibodydependent cellular phagocytosis, antibodydependent cellular cytotoxicity or complementdependent cytotoxicity (see for instance Figure 2 in (11)).
– Prior studies have found that antibodies have a limited role during infection (e.g., see Goyal et al. Viruses). How does this fit into the results?
It is true that cellular response is expected to play a greater role during infection. However, antibodies are still playing an important role to control ongoing infection as demonstrated by neutralizing monoclonal antibodies used for treatment of COVID19 (see (11)).
– Figure 2 is difficult to assess what is data and what is the model, and the model fit is difficult to see. It seems a though the 95% confidence intervals (gray shading) are also not indicative the model. Why are there multiple peaks and does the model capture that?
We thank the reviewer for this relevant comment. The reviewer can refer to the comment 7 of the reviewer 1 for our answer about this figure. A new version of Figure 2 is now proposed.
– The discussion could be improved. It lacked discussion of how the work fits into the literature, particularly other published SARSCoV2 models, how the correlate might be used in the clinic, etc. In addition, a discussion about the differences in dynamics within the nasopharynx and trachea in vaccinated individuals would have been interesting.
The Discussion section has been revised by adding:
– the perspective on extension of the model inspired by other models proposed in the literature
– one more sentence on the differences in dynamics within the nasopharynx and trachea
– a perspective on clinical implications
[Addition/modification applied in section Discussion]
The reviewer can refer to the comments 4, 5 and 10 of the reviewer 2 for the related modifications of the section Discussion.
– It's unclear to me why the reproduction number R was under 1 for the alphaCD40.RBD group when there was replication for several days in this group.
The reproductive ratio is reflecting the trend of the replication of the virus which is tending to extinct over time when R<1. However, at a given time, R could be < 1 although there are still infected cell and replicating virus. R is estimated to be <1 in the nasopharynx from the beginning and drop below 1 in the trachea from day 2. Hence, the dynamics of the observed gRNA is a contribution of the virus inoculated in the first days and some replication (below the detection threshold in the nasopharynx and with one observed blip in the trachea).
Reviewer #3 (Recommendations for the authors):
First, a suggestion regarding my overfitting concern:
• The number of nonhuman primates in each arm of these studies is understandably limited, but crossvalidation could be used to reduce the extent of overfitting even with limited replicates.
We thank the reviewer for this relevant comment. Indeed, the number of NHPs is quite limited. As recommended, we performed a bootstrap procedure with replacement (B=50 times) to obtain more robust estimates. Due to time constraints, the procedure was performed only for the model adjusted for group of treatment and the one with β adjusted the selected marker. The bootstrap parameter estimate was calculated as the median of the parameter estimates from the B bootstrap samples while the standard error of each parameter was calculated according to the definition of Thai et al. (2013) (12), which means with the SE of the lth component of the vector of parameters given by:$\hat{\text{SE}}}_{B}^{\left(l\right)}=\sqrt{\frac{1}{B1}\sum _{b=1}^{B}{\left({\hat{\theta}}_{b}^{\ast \left(l\right)}{\hat{\theta}}_{B}^{\left(l\right)}\right)}^{2}$ with $\hat{\theta}}_{b}^{\ast (l)$ being its estimate obtained at the bth iteration of the bootstrap and $\hat{\theta}}_{B}^{(l)$ the bootstrap parameter estimate. For each bootstrap sample, we paid attention to keep the 1:1:1 ratio between the 3 groups of treatment, with 6 animals selected within each group. Results were then reported in Table 1 and Table 2 below, respectively for the two models aforementioned.
In addition, we performed a convergence assessment in Monolix for these two models to evaluate the robustness of the convergence, with 15 runs. The results are shown in Author response image 1 and Author response image 2.
[Addition applied in section Methods page 27, Lines 578580]
“Furthermore, we performed a bootstrap procedure with replacement (79) (50 samples) on the optimal model to obtain and verify the robustness of the estimates (see Supplementary file 2 for the results), and to reduce potential overfitting that could result from the small number of NHPs.”
[Addition/modification applied in supplementary information file – Table S2 (now defined as Supplementary file 2)]
[Addition in the supplementary files of the Table S5 (Supplementary file 4) presenting parameter estimates when viral infectivity is adjusted for timevarying pseudoneutralization (ECLRDB) and loss rate of infected cells adjusted for group effects]
Next, two suggestions to enhance clarity:
• The rationale for the use of profile likelihood should be made explicit. Why do a few parameters using this approach and then later do rest with full maximum likelihood? Is it a limitation of the optimization software?
As mentioned by the Reviewer 1 in the comment (3b), we performed profile likelihood to determine the value of the three fixed parameters c, c_{i} and k. We used profile likelihood for these parameters because of identifiability issues. Similar to Gonçalves et al. (2021) (3), data did not provide enough information to estimate the duration of the eclipse phase and virus clearances with full maximum likelihood (e.g., not enough time points or data not enough specific to the mechanisms of interest). Profile likelihood (also referred as sensitive analyses in some papers) was then performed to roughly identify values that could be biologically relevant and lead to better estimations of the model. Other parameters are then estimated considering c, c_{i} and k fixed at the values found by profile likelihood. The reviewer can refer to comments (3a) and (3b) to see the modifications we applied to the manuscript to better explain this point.
• In the model, why separately track virions from the inoculum versus virions generated in host? If it is just to allow for different clearance rates, why do we expect these to be different?
We distinguished virions that have been inoculated and those produced by infected cells. The reviewer can refer to the comment 1. (b) of reviewer 1 for our feedback about this comment.
References
1. Czuppon P, Débarre F, Gonçalves A, Tenaillon O, Perelson AS, Guedj J, et al. Success of prophylactic antiviral therapy for SARSCoV2: Predicted critical efficacies and impact of different drugspecific mechanisms of action. PLOS Computational Biology. 1 mars 2021;17(3):e1008752.
2. Krambrich J, Akaberi D, Ling J, Hoffman T, Svensson L, Hagbom M, et al. SARSCoV2 in hospital indoor environments is predominantly noninfectious. Virology Journal. 2 juin 2021;18(1):109.
3. Gonçalves A, Maisonnasse P, Donati F, Albert M, Behillil S, Contreras V, et al. SARSCoV2 viral dynamics in nonhuman primates. PLOS Computational Biology. 17 mars 2021;17(3):e1008785.
4. Sawicki SG, Sawicki DL, Siddell SG. A Contemporary View of Coronavirus Transcription. Journal of Virology. janv 2007;81(1):20‑9.
5. Miao H, Xia X, Perelson AS, Wu H. On Identifiability of Nonlinear ODE Models and Applications in Viral Dynamics. SIAM Review. janv 2011;53(1):3‑39.
6. Gilbert PB, Fong Y, Carone M. Assessment of Immune Correlates of Protection via Controlled Vaccine Efficacy and Controlled Risk. arXiv:210705734 [stat] [Internet]. 12 juill 2021 [cité 8 févr 2022]; Disponible sur: http://arxiv.org/abs/2107.05734
7. Robinot R, Hubert M, de Melo GD, Lazarini F, Bruel T, Smith N, et al. SARSCoV2 infection induces the dedifferentiation of multiciliated cells and impairs mucociliary clearance. Nature communications. 2021;12(1):4354.
8. Guedj J, Thiébaut R, Commenges D. Joint Modeling of the Clinical Progression and of the Biomarkers’ Dynamics Using a Mechanistic Model. Biometrics. 2011;67(1):59‑66.
9. Marlin R, Godot V, Cardinaud S, Galhaut M, Coleon S, Zurawski S, et al. Targeting SARSCoV2 receptorbinding domain to cells expressing CD40 improves protection to infection in convalescent macaques. Nat Commun. 1 sept 2021;12(1):5215.
10. Goyal A, Reeves DB, Schiffer JT. Multiscale modelling reveals that early superspreader events are a likely contributor to novel variant predominance. Journal of The Royal Society Interface. 19(189):20210811.
11. Taylor PC, Adams AC, Hufford MM, de la Torre I, Winthrop K, Gottlieb RL. Neutralizing monoclonal antibodies for treatment of COVID19. Nat Rev Immunol. juin 2021;21(6):382‑93.
12. Thai HT, Mentré F, Holford NHG, VeyratFollet C, Comets E. Evaluation of bootstrap methods for estimating uncertainty of parameters in nonlinear mixedeffects models: a simulation study in population pharmacokinetics. Journal of Pharmacokinetics and Pharmacodynamics. févr 2014;41(1):15‑33.
https://doi.org/10.7554/eLife.75427.sa2Article and author information
Author details
Funding
Agence Nationale de la Recherche (ANR10LABX7701)
 Yves Levy
 Rodolphe Thiébaut
Agence Nationale de la Recherche (ANR11 1018 INBS0008)
 Roger Le Grand
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We would like to thank J Guedj and O Terrier for fruitful discussions on the model definition. We thank S Langlois, J Demilly, N Dhooge, P Le Calvez, M Potier, JM Robert, T Prot, and C Dodan for the NHP experiments; L Bossevot, M Leonec, L MoenneLoccoz, M CalpinLebreau, and J Morin for the RTqPCR, ELISpot and Luminex assays, and for the preparation of reagents; AS Gallouët, M GomezPacheco, and W Gros for NHP Tcell assays and flow cytometry; B Fert for her help with the CT scans; M Barendji, J Dinh, and E Guyon for the NHP sample processing; S Keyser for the transports organization; F Ducancel and Y Gorin for their help with the logistics and safety management; I Mangeot for her help with resources management and B Targat contributed to data management. The monkey and syringe pictures in Figure 1 was created with BioRender.com. This work was supported by INSERM and the Investissements d’Avenir program, Vaccine Research Institute (VRI), managed by the ANR under reference ANR10LABX7701. MA has been funded by INRIA PhD grant. The Infectious Disease Models and Innovative Therapies (IDMIT) research infrastructure is supported by the 'Programme Investissements d’Avenir', managed by the ANR under reference ANR11INBS0008. The Fondation Bettencourt Schueller and the Region IledeFrance contributed to the implementation of IDMIT’s facilities and imaging technologies used to define volume of respiratory tract. The NHP study received financial support from REACTing, the Fondation pour la Recherche Medicale (FRM; AMCoVPath). We thank Lixoft SAS for their support. Numerical computations were in part 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 Miles Davenport and Frederik Graw as Senior Editor and Reviewing Editor of our paper, respectively, and the three anonymous reviewers for their time and their constructive comments.
Ethics
Cynomolgus macaques (Macaca fascicularis), aged 37–66 months (18 females and 13 males) and originating from Mauritian AAALAC certified breeding centers were used in this study. All animals were housed in IDMIT facilities (CEA, Fontenayauxroses), under BSL2 and BSL3 containment when necessary (Animal facility authorization #D9203202, Préfecture des Hauts de Seine, France) and in compliance with European Directive 2010/63/EU, the French regulations and the Standards for Human Care and Use of Laboratory Animals, of the Office for Laboratory Animal Welfare (OLAW, assurance number #A582601, US). The protocols were approved by the institutional ethical committee “Comité d’Ethique en Expérimentation Animale du Commissariat à l’Energie Atomique et aux Energies Alternatives” (CEtEA #44) under statement number A20011. The study was authorized by the “Research, Innovation and Education Ministry” under registration number APAFIS#24434–2020030216532863 v1.
Senior Editor
 Miles P Davenport, University of New South Wales, Australia
Reviewing Editor
 Frederik Graw, Heidelberg University, Germany
Version history
 Preprint posted: November 1, 2021 (view preprint)
 Received: November 9, 2021
 Accepted: June 22, 2022
 Accepted Manuscript published: July 8, 2022 (version 1)
 Version of Record published: July 14, 2022 (version 2)
Copyright
© 2022, Alexandre et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 830
 Page views

 242
 Downloads

 3
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Immunology and Inflammation
The coexpression of inhibitory receptors (IRs) is a hallmark of CD8+ Tcell exhaustion (Tex) in people living with HIV1 (PLWH). Understanding alterations of IRs expression in PLWH on longterm antiretroviral treatment (ART) remains elusive but is critical to overcoming CD8+ Tex and designing novel HIV1 cure immunotherapies. To address this, we combine highdimensional supervised and unsupervised analysis of IRs concomitant with functional markers across the CD8+ Tcell landscape on 24 PLWH over a decade on ART. We define irreversible alterations of IRs coexpression patterns in CD8+ T cells not mitigated by ART and identify negative associations between the frequency of TIGIT+ and TIGIT+ TIM3+ and CD4+ Tcell levels. Moreover, changes in total, SEBactivated, and HIV1specific CD8+ T cells delineate a complex reshaping of memory and effectorlike cellular clusters on ART. Indeed, we identify a selective reduction of HIV1 specificCD8+ Tcell memorylike clusters sharing TIGIT expression and low CD107a that can be recovered by mAb TIGIT blockade independently of IFNγ and IL2. Collectively, these data characterize with unprecedented detail the patterns of IRs expression and functions across the CD8+ Tcell landscape and indicate the potential of TIGIT as a target for Tex precision immunotherapies in PLWH at all ART stages.

 Genetics and Genomics
 Immunology and Inflammation
Ageassociated DNA methylation in blood cells convey information on health status. However, the mechanisms that drive these changes in circulating cells and their relationships to gene regulation are unknown. We identified ageassociated DNA methylation sites in six purified bloodborne immune cell types (naive B, naive CD4^{+} and CD8^{+} T cells, granulocytes, monocytes, and NK cells) collected from healthy individuals interspersed over a wide age range. Of the thousands of ageassociated sites, only 350 sites were differentially methylated in the same direction in all cell types and validated in an independent longitudinal cohort. Genes close to ageassociated hypomethylated sites were enriched for collagen biosynthesis and complement cascade pathways, while genes close to hypermethylated sites mapped to neuronal pathways. In silico analyses showed that in most cell types, the ageassociated hypo and hypermethylated sites were enriched for ARNT (HIF1β) and REST transcription factor (TF) motifs, respectively, which are both master regulators of hypoxia response. To conclude, despite spatial heterogeneity, there is a commonality in the putative regulatory role with respect to TF motifs and histone modifications at and around these sites. These features suggest that DNA methylation changes in healthy aging may be adaptive responses to fluctuations of oxygen availability.