2.1.1 Structural conservation fraction.

Using the structures for the trimeric spike of various coronaviruses (S1 Table), atoms were removed from each trimer so that only a monomer remained. Thus, the structural conservation fraction was calculated using a monomer rather than the full trimer, but the monomer retained the conformation that it had in the trimer. Using the monomer also removes redundant structures and is necessary to construct the structural sequence alignment because sequence alignments use single proteins. The STRAL-SV server [35] was used to align these monomeric spike protein structures. The STRAL-SV server addresses the fact that multi-domain proteins can adopt multiple conformations, such as the spike protein’s RBD being oriented down in one structure and up in another structure, which would have confounded a standard global alignment. In STRAL-SV, a sliding-window approach divided the monomer structures into overlapping fragments of 90 amino acids, and structural alignments were performed between the fragments. This process identified residues across different coronavirus spikes that structurally overlap in the alignment. This structural alignment process tolerates mutations that maintain backbone structure, as single amino acid substitutions often do, and it penalizes mutations that alter backbone structure. This data was used to construct a structure-based sequence alignment. A table is first constructed (see Fig 1) where the residues of the spike proteins are the rows, the column names correspond to coronaviruses (such as MHV, HKU1, and HKU2), and the entries of the table are the 20 amino acids or a placeholder entry of “.” The placeholder entry of “.” means that the coronavirus spike protein was not able to align at that residue.

The structural conservation fraction was calculated from the structure-based sequence alignment using a particular coronavirus as a reference. For each residue, the structural conservation fraction was the fraction of other coronaviruses that are able to structurally align at that residue. We illustrate this with an example from Fig 1‘s table, assuming for simplicity that the first row corresponds to residue 1. If we select residue 1 and MHV as a reference coronavirus, the other coronaviruses are HKU1 and HKU2. Only HKU1 is able to structurally align, so the structural conservation fraction is 0.5. It is also possible that the reference coronavirus did not align at the chosen residue, in which case the structural conservation fraction is 0. For example, consider residue 1 using HKU2 as a reference coronavirus. Since the entry is “.”, then the structural conservation fraction is 0. Note that Fig 1‘s table is truncated and is shown for illustration purposes only.

2.1.2 Biochemical conservation fraction from sequence analysis.

For each coronavirus, a set of several hundred spike sequences were collected from the NCBI Protein database and processed (see SI for details of collection and processing). All coronavirus alignments were then concatenated together into a single multiple sequence alignment as illustrated in S2 Fig.

The biochemical conservation fraction was calculated from the multiple sequence alignment using a particular coronavirus as a reference. For each residue, the biochemical conservation fraction was the fraction of sequences that have an amino acid in the same class (hydrophobic, polar, positive, or negative) as the reference coronavirus’s consensus (most common) amino acid. Thus, this calculation tolerates mutations if they are biochemically similar.

To illustrate this procedure using S2 Fig, assuming that the first row corresponds to residue 1 for simplicity, suppose we wish to calculate the biochemical conservation fraction of residue 1 using MHV as a reference coronavirus. The consensus amino acid from the MHV alignment is F, which is a hydrophobic amino acid. The biochemical conservation fraction is the fraction of amino acids in the first row of the concatenated alignment that are also hydrophobic.

2.1.3 Conservation fraction.

The conservation fraction was the average of the structural conservation fraction and the biochemical conservation fraction. For the design of antigens that may protect against SARS-CoV-2 variants, only severe acute respiratory syndrome coronavirus (SARS-CoV) and SARS-CoV-2 data were considered. Since SARS-CoV and SARS-CoV-2 are closely related, the threshold conservation fraction below which a residue was identified as variable was chosen to be 0.99.

For pan-coronavirus vaccine target identification, we considered 12 coronaviruses from different genera (S1 Table), which constituted an analysis of 12 spike protein structures and ~4000 spike protein sequences. We set a threshold conservation fraction of 0.8, which was chosen because it was the largest value at which we observed a conserved region comparable in size to an antibody footprint. The footprint size was estimated to be 7.8±0.4 nm² and was calculated as an average over ten RBD-antibody structures (see SI for details).

2.4.1 Binding free energy representation.

The BCR paratope and antigen epitope were each represented by a vector of residues. Both the paratope and epitope were 50 residues long, a length which was chosen based on the number of residues in the class 1 and class 2 epitopes [34,51–53]. For antibodies that did not have the epitope residues reported, the residues were found from the RBD-antibody structure using the PISA server [54–56]. This procedure allows identification of RBD residues at the RBD-antibody interface with buried surface areas greater than 0 Ų. In order to choose residues that were consistently part of the class 1 or class 2 epitope, only residues that were bound by 3 or more class 1 or class 2 antibodies were included. Using the conservation analysis described above, 9 of the 50 residues had conservation fractions above 0.99 and were assigned as conserved. The remaining 41 residues were assigned as variable.

For a seeding B cell, each residue, k, was associated with a number PAR(k), which was initially sampled from a uniform distribution between -0.18 and 0.90. After AM began, the value of PAR(k) changed due to mutations, and the bounds on these values were -1 and 1.5. The epitope residues were represented by a number EPI(k), with values of +1 for WT residues and negative values for mutated residues depending on how biochemically different the mutated residue is from the WT. For a mutation in the same class (such as hydrophobic to hydrophobic), the epitope residue was -1. For a mutation from hydrophobic to polar, polar to charged, or vice-versa, the epitope residue was -2. For a mutation from hydrophobic to charged or vice-versa, the epitope residue was -3. For a mutation from positively charged to negatively charged or vice-versa, the epitope residue was -4. To illustrate the construction of the epitope vector, we will use the Delta variant, which is defined using the mutations L452R and T478K, as an example. L452R is a hydrophobic to positive mutation (corresponding to a value of -3), and T478K is a polar to positive mutation (corresponding to a value of -2). Therefore, one residue of the epitope vector will have a value of -3, another will have a value of -2, and the remaining residues will have values of +1. The order of the residues does not matter in our simple model. Additionally, if the mutated residues are scaled by a factor of 0.5, meaning that the L452R and T478K mutations correspond to values of -1.5 and -1, then the results of the model do not change (S4 Fig).

The binding free energy was calculated according to (1) where PAR(k) was the paratope residue at position k and EPI(k) was the epitope residue at position k. Essentially, the binding free energy is a dot product between the paratope and epitope vectors. Larger values of E corresponded to stronger binding. This coarse-grained representation of the binding free energy was not based on structural considerations and had limitations: (1) the contribution of each residue to the binding free energy was weighted equally; (2) the order of the residues did not matter; (3) only RBD substitutions could be considered; (4) biochemically similar mutations were treated equivalently. The first limitation, that each residue is weighted equally, is a particularly significant assumption since it is known that certain residues more easily abrogate binding upon mutation than others [46]. The effect on our model is that our estimates of the potency of antibody responses to variants are likely to be conservative. Weighting all residues equally enables the virus to evade the response more easily because mutations at all sites, rather than a few, can abrogate binding. Additionally, we considered the possibility of reweighting residues according to their experimentally measured escape fractions. Such an approach is not appropriate because the escape fractions were measured using antibodies generated from WT infection, but our vaccine will generate different antibodies corresponding to different escape fractions. Reweighting the residues would then lead to binding free energies that are inaccurately biased toward WT escape residues. For this reason, we believe that our approach of equally weighting residues is more conservative.

That said, we emphasize that this limited representation of the binding free energy is not quantitative and is not meant to accurately reproduce experimental values. Rather, it is a qualitative tool meant to facilitate comparison between different antigens and obtain mechanistic insights into AM when used with the rest of the model. The utility of this approach has been demonstrated in previous work, which used a very similar model to provide insights into antibody evolution against HIV antigens that were then validated in animal models [22,39,57].

2.4.2 B cell expansion and mutation in the dark zone.

Upon initial immunization, 10 naive B cells seeded a new GC. Increasing the number of seeding cells does not change our qualitative results (S5 Fig). The B cells then expanded without mutation or selection until the population reached 5120 cells. In AM, activation-induced cytidine deaminase (AID) introduces somatic hypermutations (SHMs) into the BCR. This was modeled by mutating the paratope residues at a rate of 0.14 per sequence per division (B cells divided twice per GC cycle). SHMs had a 0.5 probability of being lethal, 0.3 probability of being silent, and 0.2 probability of modifying the binding free energy [58]. For mutations that modified the binding free energy, a random residue on the paratope was chosen, and the change in binding free energy for that residue was sampled from a shifted log-normal distribution: (2) where r was a standard normal random variable, ϵ was a shift parameter, μ was the mean of the log-normal distribution, σ was the standard deviation of the log-normal distribution, and δ limited the maximum ΔE. We used the same parameters as Sprenger et. al (ϵ = 3, μ = 1.9, σ = 0.5, δ = 1), which were set to match experimental distributions of changes in binding free energies between proteins due to single-residue mutations [59].

2.4.3 Selection in the light zone.

In the light zone of the GC, B cells underwent selection based on their binding affinity to antigen. First, each B cell had a probability of successfully internalizing antigen according to (3) where i was a sum over all antigens that the B cell encounters, c_i was the concentration of antigen i, E_i was the binding free energy to antigen i, E_act was a reference free energy, and e_scale was a scaling parameter. We used the same parameters as Sprenger et. al (e_scale = 0.08, E_act = 9k_bT). For a particular immunization scheme, the concentration that produces the highest antibody titers was chosen. The sum in Eq 3 did not apply for WT doses that only contain a single antigen. For cocktails, P_internalize depended on whether the B cell encounters one antigen at a time or all antigens on the FDC. If the B cell encountered one antigen at a time in each cycle, then a random antigen was chosen, and there was no summation. If the B cell encountered all antigens during each cycle, then the sum included all antigens i. In any given cycle, a B cell has a few chances to be positively selected. So, if the antigen concentration on FDCs is sufficiently high, then the sum over all antigens is likely to be more realistic.

B cells that did not internalize antigen underwent apoptosis and were removed from the population, while B cells that internalized antigen competed for T cell help. To model T cell selection if B cells encountered one antigen per cycle, the B cells were ranked according to their binding free energy to the encountered antigen in the last cycle, and the top F_cut,help fraction were positively selected. If B cells instead encountered all antigens per cycle, then the binding free energy ranking was based on all B cell–antigen pairs. The top F_cut,help fraction of pairs were selected, and each B cell survived with a probability equal to the frequency that it appeared in the top F_cut,help pairs. For example, suppose there are 6 antigens and thus 6 B cell–antigen pairs for a particular B cell. If there are 5000 B cells, then there are 30,000 possible B cell–antigen pairs in total. If 4 of the B cell–antigen pairs are in the top F_cut,help*30,000 pairs, then that B cell survives with a probability of . Following Sprenger et. al, we set F_cut,help = 0.7.

2.4.4 Recycling, exit, GC reaction termination, and seeding a new GC.

Following Sprenger et. al, positively selected B cells had a 0.7 probability of being recycled for subsequent cycles and a 0.3 probability of exiting the GC as either a plasma or memory cell. The GC reaction stopped if the B cell population exceeded the initial population of 5120 cells, if the number of cycles exceeded 250, or if all B cells died. The first condition is a proxy for the B cells internalizing all the antigen, the second is a proxy for antigen decay over time, and the third is GC extinction. Sprenger et al. set these parameters, along with the parameters from Section 2.4.3 and the bounds on paratope residues from Section 2.4.1, so that a simulation of single-antigen immunization matches experimentally measured GC dynamics of single-antigen immunization [60]. Upon a second immunization, a new GC was seeded using 10 B cells, which were a mixture of memory and naive cells. The fraction of memory cells was 1 (all memory cells), except in Section 3.2.4 where we studied the effect of varying this fraction.

2.4.5 Titer calculation.

The B cell population at the end of each GC reaction was analyzed for antibody titers. Identical B cells were grouped together as clones, and the binding free energy against an antigen sequence was calculated for each clone. The antibody titer was the total number of B cells that bind to the antigen with a free energy above E_th = 17 k_bT. Since the binding free energy representation is based on the class 1 and 2 epitopes, the titer calculated here theoretically corresponds to neutralizing antibody titer. While the titers quantitatively depended on E_th, the rank order of titers for different immunization schemes that we studied did not. So, the qualitative results did not depend on E_th (S6 Fig). The simulation was carried out N = 200 times, and the average antibody titer for antigen i was calculated as (4) where size(j) was the number of B cells belonging to clone j, and H represents the step function. Error bars were standard deviations across 10 sets of N simulations. The titers thus obtained may be considered to represent the titers generated on average in a typical vaccinated person. This model only considered titers upon completion of each GC reaction. In reality, the titers would be defined by the plasma cells exiting the GC secreting antibodies continuously during the GC reaction. We carried out simulations where plasma cells producing antibodies are continuously exported during the GC reaction, and the qualitative results remain the same (S7 Fig).

2.4.6 Mean panel titer calculation.

To analyze the coverage provided by antibodies against multiple different antigens, the binding free energy of each B cell clone with a test panel of 100 antigens was calculated. Increasing the panel size to 1000 antigens did not change the results (S8 Fig). The antibody titer against each panel antigen was calculated, and the mean panel titer was the mean titer across all 100 panel antigens.

Mean panel titers for several antigen panels were calculated. In the most general panel, residues took values of either +1 or -1, and the probability that the residue is +1 was the conservation fraction. This panel produced almost always resulted in mean panel titers of 0 because the panel antigens had too many mutations at the same time, which are unlikely to emerge in a variant antigen. So, we instead considered panels in which the antigens are restricted to N_mut mutations and varied N_mut. We considered cases in which the N_mut mutations occur in any variable residue and in which the mutations only occur in the 10 residues that are also mutated in the designed antigens. We also considered panels in which mutated residues take on different values (-1, -2, -3, and -4).

3.2.1 Reduced titers against SARS-CoV-2 variants following WT immunization.

We simulate AM following one WT immunization, which models infection or one vaccine dose, and two WT immunizations, which models two vaccine doses or infection followed by one vaccine dose. In Fig 4A, the anti-WT titers are high for one immunization (blue bars in various panels) and even higher for two immunizations (orange bars in various panels). In contrast, the titers against variants are reduced depending on the mutational distance between the variant and WT. For example, following two WT immunizations, anti-Beta/Gamma titers are reduced ~6-fold relative to anti-WT titers, anti-Delta titers are reduced ~2.5-fold, and anti-Alpha titers are nearly equivalent to anti-WT titers. This is due to the mutational distances of these variants as the Beta and Gamma variants have 3 RBD substitutions, the Delta variant has 2 RBD substitutions, and the Alpha variant has 1 RBD substitution. This also recapitulates experimental data demonstrating that the Alpha variant minimally reduces neutralizing titers while Beta, Delta, and Gamma significantly reduce titers [1,2,4–6,61–65]. These studies report differences between anti-Beta and anti-Gamma titers due to the importance of specific mutations or mutations outside of the RBD, but this coarse-grained model is not structurally detailed enough to recapitulate those differences.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Titers against certain variants are low following WT immunization but high following cocktail immunizations with designed antigens.

(a) Titers against WT and variants following different WT or cocktail immunization schemes. Sequences 2 and 3 are similar to sequence 1, so they are omitted for clarity. Sequences 5 and 6 are also omitted because they are similar to sequence 4. Beta and Gamma variants carry biochemically similar RBD mutations, so they are equivalent in this model. (b) Mean panel titers as a function of the number of mutations in panel antigens for different immunization schemes. Two types of panels are considered: one in which panel antigens can be mutated in any variable residue (All variable residues), and one in which panel antigen mutations can only occur in the 10 residues that are also mutated in sequences 1–6 (Seq1-6 mutated residues). It is possible for the panel antigens to have more mutations than in sequences 1–6, which have 5 mutations each. For example, for n = 8, every panel antigen has 8 mutations among the 10 residues that were mutated in the designed antigens.

https://doi.org/10.1371/journal.pcbi.1010563.g004

3.2.2 Titers are higher for immunization with a cocktail of designed antigens than WT and depend on the choice of antigens in the cocktail.

We then used our model to identify an optimal immunization scheme with our 6 designed antigens. We simulated a single immunization with all 6 antigens (Seq1-6), two immunizations with all 6 antigens (Seq1-6 | Seq1-6), a single immunization with sequences 1–3 (Seq1-3) and an immunization of sequences 1–3 followed by an immunization with sequences 4–6 (Seq1-3 | Seq4-6). We assume that the antigens are homogeneously distributed on the FDCs at high concentrations, so a B cell encounters all antigens in each cycle, but assuming that a B cell encounters one antigen in each cycle does not affect the results (see SI for details). The titers are reported in Fig 4A using different colored bars.

For anti-WT titers, WT immunizations provide greater protection than cocktail immunizations, which is expected because the cocktail antigens are mutated from the WT. However, the boosted cocktails (Seq1-6 | Seq1-6 and Seq1-3 | Seq4-6) still retain moderately high titers against WT, indicating that they may still be protective. The immunizations that start with Seq1-3 (Seq1-3 and Seq1-3 | Seq4-6) have slightly higher anti-WT titers than the analogous immunization starting with Seq1-6 (Seq1-6 and Seq1-6 | Seq1-6). This is because the antigens in Seq1-6 are mutated across 10 residues while the antigens in Seq1-3 are mutated across 5 residues, so the antibodies resulting from immunization with Seq1-3 will have more residues that resemble the WT sequence (more positive in the model). Since more positive residues increases the binding free energy with WT antigen, immunization with Seq1-3 results in higher anti-WT titers. The Alpha variant is similar to the WT, so the same reasoning applies to anti-Alpha titers.

Anti-Seq1 titers are low following WT immunization because sequence 1 contains numerous mutations. Seq1-6 immunization produces lower anti-Seq1 titers than Seq1-3 immunization because sequences 2 and 3 contain mutations in the same residues as sequence 1 while sequences 4–6 contain mutations in different residues from sequence 1. However, immunization with Seq1-6 | Seq1-6 versus Seq1-3 | Seq4-6 are nearly equivalent in anti-Seq1 titers. This is because the anti-Seq1 antibodies provided by the boost with Seq1-6 are much greater than the anti-Seq1 antibodies provided by the boost with Seq4-6, which compensates for the lower anti-Seq1 antibodies generated by priming with Seq 1–6 compared to Seq 1–3.

Anti-Seq2 and anti-Seq3 titers showed similar results to anti-Seq1 titers, so they are omitted for brevity. These titers are similar, even though sequence 2 and 3 have different mutations from sequence 1, because the mutations are in the same residues. Parenthetically, this phenomenon of titers mostly depending the site of mutation rather than the actual mutation is observed to some extent in experimental data, as escape residues often can abrogate antibody binding through numerous possible mutations [44], but there are some cases of biochemically similar mutations differentially affecting antibody binding. Our finding noted above illustrates the qualitative, but not quantitative, nature of our model.

For anti-Seq4 titers, immunization with Seq1-6 produces higher titers than immunization with Seq1-3 because sequences 1–3 are all mutated in different residues from sequence 4. Immunization with Seq1-6 | Seq1-6 produces higher anti-Seq4 titers than immunization with Seq1-3 | Seq4-6, which is somewhat unexpected because one might expect that a boost with Seq4-6 would produce higher anti-Seq4 titers than a boost with Seq1-6. The origin of this result is that, because Seq 4–6 are very different from Seq 1–3, the memory B cells produced by priming with Seq 1–3 are poorly adapted to Seq 4–6. Therefore, the boost with Seq4-6 results in many memory B cells in the secondary GC dying (S9 Fig), which reduces the anti-Seq4 titer from the boost. Similar to how anti-Seq2 and anti-Seq3 titers are similar to anti-Seq1 titers, anti-Seq5 and anti-Seq6 titers are similar to anti-Seq4 titers and are not shown for brevity.

For anti-Beta/Gamma and anti-Delta titers, cocktail immunization universally produces higher titers than WT immunization with the same number of immunizations (e.g. Seq1-6 is higher than WT and Seq1-6 | Seq1-6 is higher than WT | WT). The titers resulting from immunizations with different cocktails are roughly equivalent, such as the titers resulting from Seq1-6 | Seq1-6 versus Seq1-3 | Seq4-6.

Fig 4B shows the mean panel titer as a function of the number of mutations in the panel antigens from immunization with WT and cocktails. Mutated residues change from 1 to -4, which is chosen because the vaccine should ideally protect against emergent strains with biochemically dissimilar mutations. Using values such as -3, -2, or -1 instead of -4 does not change the qualitative results (S10 Fig). If the panel mutations occur in any variable residue (Fig 4B, left panel), then the mean panel titer decreases at a similar rate for the cocktails as for the WT immunization. This is because there are few conserved residues, so the contribution of conserved regions to the binding free energy is not sufficient to generate bnAbs that can protect against strains bearing many mutations at other residues. This is consistent with experiments showing that few antibodies targeting class 1 or class 2 epitopes of the RBD are broadly protective against all single RBD substitutions [46] let alone those bearing multiple such mutations.

We next studied whether immunization with our designed antigens can protect against panels of variants with mutations only in the residues that are also mutated in the designed antigens (Fig 4B, right panel). Note that the actual amino acids at the 10 residues can be different from those in our designed antigens. In this case, the mean panel titers do not significantly drop for immunization with Seq1-6 and Seq1-6 | Seq1-6. That is, vaccination using these schemes protects against variants with up to 10 mutations in particular variable residues. The mean panel titers for Seq1-3 | Seq4-6 do drop as the number of mutations increases, but less so than WT immunization. These results suggest that Seq1-6 | Seq1-6 may be an effective immunization scheme that could protect against current variants and some that may emerge in the future.

3.2.3 Titers against variants depend on previous exposure to WT and the number of mutations in the variant.

Many individuals have been previously exposed to WT SARS-CoV-2 either through vaccination or natural infection, which might affect the antibody response generated upon immunization with a cocktail of our designed antigens. This is because memory B cells generated during prior exposure may compete with naïve B cells that seed GCs. To study this process, we use our model to predict antibody response given previous WT immunizations. Fig 5 graphs antibody titer as a function of the number of previous immunizations with WT for different immunization schemes. Although naive cells are likely involved in seeding new GCs in vivo [66], here we assume that only memory cells seed new GCs to maximally illustrate the effects of memory. Section 3.2.4 studies the effect of including naive cells.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Titers against variants are affected by previous WT immunizations.

Titers against various antigens (indicated in panel title) as a function of the number of previous WT immunizations. Sequences 2 and 3 are similar to sequence 1, so they are omitted for clarity. Sequences 5 and 6 are also omitted because they are similar to sequence 4. Beta and Gamma variants carry biochemically similar RBD mutations, so they are equivalent in this model.

https://doi.org/10.1371/journal.pcbi.1010563.g005

First, we compare WT versus WT | WT and WT | WT | WT immunization as well as cocktail versus WT | cocktail and WT | WT | cocktail immunization. Intuitively, anti-WT and anti-Alpha titers steadily increase as the number of previous immunizations with WT increase. However, titers against variants with numerous mutations are not always higher for WT | cocktail and WT | WT | cocktail compared to just cocktail immunization. For example, the anti-Seq1 titers decrease from 0 to 1 previous WT immunization for Seq1-6 | Seq1-6 and Seq1-3 | Seq4-6. This is because the mutational distance between the previous WT immunizations and the cocktails induces GC extinction, thereby limiting antibody development. Nonetheless, the decrease is not large even though GCs are seeded entirely by memory cells.

Next, we compare WT followed by WT immunization versus WT followed by cocktail immunization (for example, WT | WT versus WT | cocktail). For variants with numerous mutations, immunization with the cocktails after the first exposure to WT show higher titers than boosting with WT. For example, anti-Seq1 and anti-Seq4 titers are significantly higher in the WT | Seq1-6 | Seq1-6 and WT | Seq1-3 | Seq4-6 immunizations than the WT | WT immunization. This is because Seq1 and Seq4 each have numerous mutations, so the WT immunizations are not protective against Seq1 and Seq4. Anti-Delta, anti-Alpha, and anti-WT titers are higher for WT | WT immunization than WT | cocktail immunization because of the low number of mutations in these variants. Since the Beta and Gamma variants have more mutations than the Delta and Alpha variants, the anti-Beta/Gamma titers are equivalent between WT | cocktail and WT | WT immunizations.

3.2.4 Sequential immunization of WT and cocktail produces significantly higher titers against current variants than sequential immunization of WT and WT when naive B cells also seed new GCs.

Anti-Beta/Gamma and anti-Delta titers are not higher for cocktail immunizations compared to WT immunizations under specific conditions (assuming 1–2 previous immunizations with WT have been given and new GCs are seeded entirely by memory cells). However, if naive cells also seed GCs, then the titers for cocktail immunization following immunizations with WT become significantly higher. Illustrating this, Fig 6 shows titers as a function of the memory cell fraction given 2 previous WT immunizations. At a memory cell fraction of 0.5, anti-Beta/Gamma titers are around 10 to 15-fold higher for cocktail immunizations than WT immunizations. This is a large difference compared to the memory cell fraction of 1, where the titers are nearly equivalent, suggesting that cocktail immunization is less sensitive to the memory cell fraction. Also, it is noteworthy that decreasing the memory cell fraction decreases the rank orders of the WT | WT | WT and WT | WT | WT | WT immunizations, while the rank orders of the WT | WT | cocktail immunizations are preserved. The reason for this phenomenon is as follows: in the WT | WT | cocktail immunization, the WT and cocktail antigens are sufficiently different so that naïve cells can better develop against Beta and Gamma variants than in the WT | WT | WT immunization. Thus, decreasing the memory cell fraction (increasing naïve cell fraction) increases the difference between WT | WT | cocktail immunization and WT | WT | WT immunization. A similar, but less pronounced, effect occurs for anti-Delta titers, as the titers resulting from WT | WT | cocktail immunization are around 1.5 to 2-fold higher than WT | WT | WT immunization. This decrease occurs but is less pronounced for anti-Delta titers because the Delta variant has 2 RBD mutations whereas the Beta and Gamma variants have 3 RBD mutations.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Titers against current variants are affected by the fraction of GC-seeding cells that are memory cells.

Titers against Beta/Gamma and Delta variants as a function of the fraction of GC-seeding cells that are memory cells. Beta and Gamma variants carry biochemically similar RBD mutations, so they are equivalent in this model.

https://doi.org/10.1371/journal.pcbi.1010563.g006

3.3 The S2’ cleavage site is structurally conserved and is a potential pan-coronavirus vaccine target.

Finally, we shift our focus from the designed RBD antigens and consider the problem of identifying a suitable target for a pan-coronavirus vaccine. To identify a pan-coronavirus vaccine target, we used the method outlined in Section 2.1 to calculate spike conservation across 12 coronaviruses from different genera. Fig 7A illustrates conservation using a threshold fraction of 0.8. None of the coronaviruses are conserved in the RBD with respect to all other coronaviruses, which is simply because different coronaviruses bind different receptors. For example, SARS-CoV-2 and SARS-CoV both bind to ACE2, leading them to have similar RBDs, while the MERS RBD is different because it binds to DPP4 (S1 Table).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. A conserved epitope in the S2 domain overlaps the S2’ cleavage site and fusion peptide.

(a) Spike protein structures colored by conservation fraction using each coronavirus as a reference. Green residues have conservation fractions above 0.8, blue residues have conservation fractions below 0.8 and are not in the RBD, and red residues have conservation fractions below 0.8 and are in the RBD. (b) Spike structure colored by conservation fraction using SARS-CoV-2 as a reference. The residues corresponding to the S2’ cleavage site and fusion peptide are outlined in yellow. This structure is shown in a separate panel in order to avoid covering up some residues with the outline.

https://doi.org/10.1371/journal.pcbi.1010563.g007

In contrast to this, the S2 domain of the spike possesses residues that are structurally conserved across many coronaviruses. In particular, 9 out of 12 of the coronaviruses share a conserved region of residues in the S2 subunit despite coming from diverse lineages (the residues in this region are provided in S3 Table). For example, SARS-CoV-2 and IBV both share this conserved region even though SARS-CoV-2 is a betacoronavirus that infects humans and IBV is a gammacoronavirus that infects birds.

We also considered different weights between the structural conservation fraction and the biochemical conservation fraction in S11 Fig. The effect of decreasing the structural conservation fraction weight makes some individual residues in the RBD become classified as conserved, which is because the variability across RBD structures is diminished. Nonetheless, the conserved patch in the S2’ site is always visible, even when the structural conservation fraction weight drops to 0.25.

These residues are conserved in many coronaviruses because they overlap with the S2’ cleavage site and fusion peptide (Fig 7B). Following S1/S2 cleavage, further cleavage at the S2’ site facilitates the release of the fusion peptide, which fuses the viral and host membranes so that the virus can infect the cell. Since the cleavage site and fusion peptide are functionally important, they must be conserved. It follows that if the virus is subject to evolutionary pressure by S2’-targeting antibodies, it is more difficult to develop escape mutations because potential mutations must also preserve S2’ site function.

Lastly, enzymes such as furin, TMPRSS2, and cathepsin [67–69] must be able to access cleavage sites, which means that the cleavage sites are exposed and accessible to proteins. Furthermore, this limits the virus’s ability to develop spike densities that sterically hinder antibodies from accessing the S2’ site. Electron tomography estimates that SARS-CoV-2 virions have 25–127 spikes and a diameter of 80 nm [70]. Even at the dense estimate of 127 spikes, this produces a density of 0.6 spikes per 100 nm², which is half as dense as influenza (1.2 spikes per 100 nm²) [25,71]. Thus, compared to the broadly neutralizing epitope on the stem of the influenza hemagglutinin spike, from a steric point of view, generating antibodies to the stem of the coronavirus spike may be more feasible. Taken together, these attributes suggest that targeting the S2’ site could serve as a potential pan-coronavirus vaccine strategy.