Design of cross-reactive antigens with machine learning and high-throughput experimental evaluation

TL;DR

This paper describes using machine learning to design antigens that can protect against multiple variants of a pathogen, tested with fHbp from Neisseria meningitidis.

Contribution

A novel integrated strategy combining machine learning and experimental validation to design cross-reactive antigens.

Findings

Designed fHbp mutants that transfer conformational epitopes while maintaining cross-reactive epitope binding.

Confirmed structural integrity of top mutants using biophysical and x-ray crystallographic methods.

Developed an iterative antigen design platform to accelerate broadly protective vaccine development.

Abstract

Selecting an optimal antigen is a crucial step in vaccine development, significantly influencing both the vaccine’s effectiveness and the breadth of protection it provides. High antigen sequence variability, as seen in pathogens like rhinovirus, HIV, influenza virus, complicates the design of a single cross-protective antigen. Consequently, vaccination with a single antigen molecule often confers protection against only a single variant. In this study, machine learning methods were applied to the design of factor H binding protein (fHbp), an antigen from the bacterial pathogen Neisseria meningitidis. The vast number of potential antigen mutants presents a significant challenge for improving fHbp antigenicity. Moreover, limited data on antigen-antibody binding in public databases constrains the training of machine learning models. To address these challenges, we used computational models to predict fHbp properties and machine learning was applied to select both the most promising and informative mutants using a Gaussian process (GP) model. These mutants were experimentally evaluated to both confirm promising leads and refine the machine learning model for future iterations. In our current model, mutants were designed that enabled the transfer of fHbp v1.1 specific conformational epitopes onto fHbp v3.28, while maintaining binding to overlapping cross-reactive epitopes. The top mutant identified underwent biophysical and x-ray crystallographic characterization to confirm that the overall structure of fHbp was maintained throughout this epitope engineering experiment. The integrated strategy presented here could form the basis of a next-generation, iterative antigen design platform, potentially accelerating the development of new broadly protective vaccines.

Figures

Illustration of the mAb binding locations on the surface of fHbp. Aligned overlays of independent crystal structures, each containing a single Fab antibody fragment bound to fHbp v1.1 or fHbp v3.28; PDB: 2YPV, 5T5F, 5O14, 6H2Y, and 6XZW (Riahi et al., 2021; Desautels et al., 2024; Bianchi et al., 2019; Lopez-Sagaseta et al., 2018; Malito et al., 2013). (A) Bound Fab structures of the fHbp v1.1 specific mAbs are depicted in green and blue. (B) Bound Fab structures of the cross-reactive mAbs 4B3, 1A12, and 1E6 are depicted in orange, purple, and red respectively. Residues in the epitopes of fHbp v1.1-specific mAbs JAR5 (blue) and 12C1 (green) are indicated on the fHbp surface.

Round 1 summary of computational and experimental outcomes. (A) Interface ΔΔG calculated with FoldX for binding of fHbp v3.28-derived sequences to mAbs JAR5 and 12C1. ΔΔG values for binding to 1A12, 4B3, and 1E6 were also calculated (Supplementary Table S1). Of these sequences, and including v3.28, 41 were selected for production and are shown as red dots; sequences not selected for laboratory evaluation are shown as black dots. (B) For the selected sequences, experimentally validated binding is indicated with horizontal (pink, JAR5) or vertical (purple, 12C1) bars. Sequence m00006, which was selected as the starting sequence for Round 2, is at (−1.53, −2.11) in both panels.

Overview of the computational strategy for selection of diverse candidates with the GP model. All candidates within the defined search space are computationally modelled and multiple binding parameters are predicted using established biophysical tools, including STATIUM and FoldX. This data is combined in a machine learning model and the selection of candidates is driven by the training data provided. This allows the selection of a small set of candidates for experimental evaluation.

Sequences selected by the Gaussian process model for experimental evaluation. (A) Scatter plot of FoldX interface ΔΔG values (AnalyseComplex) for all 2,186 mutants evaluated computationally. The 108 mutants selected using the machine learning model and the decision rule are highlighted in red. Lines superimposed on the figure show effective selection thresholds in FoldX values corresponding to alternative sets of 108 mutants. If the 108 sequences with the best FoldX-predicted ΔΔG in binding 1A12 had been selected, points below the blue line fall in this set. Similarly, for if FoldX predicted binding of 12C1 is the only criterion, points left of the purple line are the resulting set, and for the equal sum of the two ΔΔGs for 1A12 and 12C1, points below and left of the orange line would constitute the resulting set. (B) The selected 108 points were tested in BLI. Among these, the best binders for 12C1 and 1A12 are marked, using cyan vertical and blue horizontal bars, respectively. A number of selected, strongly binding sequences do not lie in any of the sets constructed post hoc on the basis of FoldX binding predictions. (C) For each of group of sequences, sequence logos showing positions 221-223 and 249-252 are framed by the corresponding color. The sequence logo framed in black is all 2,187 mutants, including the parental sequence. (D) Pairwise distances among the 108 sequences in each selected or comparison set, using Blosum62. The black set is a size-matched, randomly selected set of 108 sequences from all 2,186 mutants considered, giving a fair comparison of intra-set sequence distances. Boxes show first quartile, median, and third quartile, while whiskers are 1.5 times the interquartile range or to the most extreme datum, whichever is narrower. Using the GP model, the decision rule traded off predicted binding performance for greater diversity in the selected set.

Alignment of the fHbp v3.28 mutant sequences with the highest affinity for mAb 12C1 and 1A12. The starting sequence selected for the machine learning design test, fHbp m000006, is shown in bold and deviations from this sequence are shown in green font. Residues identical in all sequences are in light grey. Locations allowed to mutate are highlighted.

Crystal structure of fHbp m002416. (A) Crystal structure of fHbp m002416 (blue) crystallized with JAR5 (black) PDB: 2YPV. (B) Crystal structure of fHbp m002416 bound to JAR5 (black) with fHbp m002416 colored by b-factors. (C) Crystal structure of fHbp v1.1 in complex with JAR5 (black) in space group C2221, fHbp v1.1 colored by b-factor. (D) Crystal structure of fHbp v1.1 in complex with JAR5 (black) in space group I4122, fHbp colored by b-factor. B-factor scale from high (red) to low (blue) and standardized across panels (B–D).

Tables

TABLE 1. Summary of measured affinities (Kd) for binding of five mAbs to select fHbp mutant sequences determined by BLI.

		mAb affinity
Antigen	Description	JAR5	12C1	1A12	1E6	4B3
fHbp v1.1	Wild type	200–700 pM	100–800 pM	Tight binding (no koff)	3–10 nM	1–5 nM
fHbp v3.28	Wild type	No binding	No binding	7–20 nM	2–4 nM	5–6 nM
m000019	Round 1 best mutant	254 pM	37 pM	Tight binding (no koff)	800 pM	2.4 nM
m000006	Round 2 starting sequence	Tight binding (no koff)	11.1 µM*	1.8 µM	870 pM	5.5 nM
m002416	Round 2 best mutant	631 pM	293 pM	544 pM	320 pM	3.8 nM

TABLE 2. Comparison of key sequences.

Residue Location	221	222	223	249	250	251	252
fHbp v1.1	Y	N	Q	T	V	N	G
fHbp v3.28	Y	G	S	I	G	E	K
m000006	Y	G	S	T	G	E	G
Allowed mutations based on StralSV	LTY	RGN	SGQ	STI	DVG	VNE	AGS
Consensus of top binders	Y	N	G/S	T/S	V	Any	G

TABLE 3. Thermostability of fHbp samples determined by differential scanning fluorimetry.

Sample	Melting temperatures
	Tm1	Tm2
fHbp v1.1	64.4	86.2
fHbp v3.28	59.0	84.1
fHbp m002416	47.2	82.1

Full text

Introduction

Advances in vaccine technologies, such as reverse vaccinology and structural vaccinology, have significantly enhanced our ability to develop vaccines for numerous diseases (Baden et al., 2021; Crank et al., 2019; Hsieh et al., 2020; Pizza et al., 2000; Polack et al., 2020; Rappuoli, 2001; Rappuoli et al., 2016). However, a significant remaining challenge in vaccinology is antigen sequence variation, where frequent mutations in infectious pathogens allow them to evade the immune system (McCulloch et al., 2017; Rappuoli, 2007). Consequently, developing broadly protective vaccines against highly mutating pathogens, such as influenza virus, HIV, rhinovirus, and meningococcus remains a difficult task (Rappuoli, 2007; Ng’uni et al., 2020).

Neisseria meningitidis factor H binding protein (fHbp) is a key component of vaccines that protect against meningitis and life-threatening sepsis caused by serogroup B meningococcus (Pace and Pollard, 2012; Rosenstein et al., 2001; Seib et al., 2015; Serruto et al., 2012; Zlotnick et al., 2015). Over 1,300 naturally occurring sequences of fHbp have been reported (Jolley and Maiden, 2010), classified into three variant groups (Bambini et al., 2009; Masignani et al., 2003; Murphy et al., 2009). Immunization with a single fHbp antigen typically induces an immune response specific to the administered variant group. Therefore, the licensed vaccines contain multiple antigens to ensure broad protection. This makes the fHbp antigen an ideal test case for the design of broadly protective antigens in variable pathogens. In addition, at least 20 crystal structures of fHbp variants and antibody-fHbp complexes are available in the protein data bank (Berman et al., 2016), providing a basis for rational design.

A significant challenge in rational antigen design is the vast search space. Each amino acid position has 20 different possibilities, resulting in an enormous combinatorial search space for a protein of 300 amino acids. This complexity makes it difficult for humans to navigate, so antigen design is typically driven by expert knowledge and focuses on the most promising protein domains. Advances in high-performance computing, simulation, and artificial intelligence (AI)/machine learning (ML) have paved the way for new approaches to rapidly design immunogens, including those not intuitive to human experts (Beede et al., 2020; Lalmuanawma et al., 2020; Ong et al., 2020). Unsupervised machine learning methods, including the use of protein sequence-based language models (Elnaggar et al., 2022; Prihoda et al., 2022; Verkuil et al., 2022), show substantial promise in developing proteins in silico that are fit, stable, or fold in vitro. These methods can also improve antibody binding behavior in a so-called “zero-shot” context (Hie et al., 2024). However, the need for large quantities of representative training data often makes these methods unsuitable for the scientific exploration of novel, idiosyncratic, or poorly characterized systems, such as vaccine antigen design (D'Amour et al., 2022). While substantial strides have been made, even the best ML structure prediction methods are uncertain or incorrect about many antibody-antigen co-structures, even when given high computational resources (Abramson et al., 2024), and structure-enabled ML antibody design remains challenging (Bennett et al., 2024). Conversely, there is substantial power and efficiency to be gained in iterative selection of designs using data obtained from measuring the quantity of interest in the system or protein of interest, including ML-assisted directed evolution approaches (Bachas et al., 2022; Shanehsazzadeh et al., 2024; Stanton et al., 2022; Yang et al., 2019).

Another strategy employed for large-scale in silico evaluation of protein-protein binding is physics-based calculations of binding energy. There are multiple simulation tools available to predict binding with differing levels of rigor. The fastest computational methods, such as FoldX (Schymkowitz et al., 2005) and STATIUM (DeBartolo et al., 2014) are the least rigorous. More elaborate methods, including Rosetta and molecular dynamics, are more rigorous and account for some flexibility in the protein backbone. While more rigorous approaches are expected to provide more accurate predictions, none of the available tools is completely accurate (Sirin et al., 2016).

In this work, we adopt a hybrid approach to manipulating antibody/antigen binding. This approach combines the advantages of physics-based simulations, machine learning, and experimentation. It has recently been shown that limited experimental data can be supplemented by molecular simulations to develop an ML tool for the prediction of protein stability (Jokinen et al., 2018) and that aggregation of a mixture of predictions can be a fruitful method for prioritizing mutations (Riahi et al., 2021). Another related application is redesigning prophylactic antibodies to counter viral escape (Desautels et al., 2024). Here, the motivating application is the creation of a vaccine antigen that would induce the production of broadly neutralizing antibodies in vivo, using in vitro binding to a panel of known, broadly neutralizing antibodies as a proxy for this desired shift in antigenicity. Collecting at large scale such binding data is still expensive and limited, as few data are collected in each system and each antigenic protein is distinct. Consequently, public data sets only contain thousands of relevant experimental observations (Sirin et al., 2016), while ML models typically employ millions or billions of datapoints. Simulation methods provide principled binding estimates, but effectively selecting antigen designs for subsequent in vitro testing requires a rigorous algorithmic framework for optimal performance.

Overall, this study establishes the capabilities required to implement machine learning as a primary component of an integrated strategy for vaccine antigen design. Specifically, this work attempts to expand the antigenicity of fHbp variant 3.28 (v3.28) by incorporating epitopes specific to fHbp variant 1.1 (v1.1), without compromising overlapping epitopes known to elicit cross-reactive antibodies. These two sequences have 60.2% identity and 75.7% similarity. Two phases of experimentation were conducted. In the set-up phase, the initial preparatory experiments focused on generating model-specific training data and confirmed that the target outcome was achievable. Next, the primary design campaign evaluated a Gaussian process (GP) model and Bayesian active learning method for candidate selection. The antigenicity of the candidate fHbp mutants was validated experimentally by measuring the affinity to a panel of antibodies (Bianchi et al., 2019; Lopez-Sagaseta et al., 2018; Malito et al., 2013; Malito et al., 2016; Veggi et al., 2020). Several fHbp v3.28 candidate antigens selected for experimental evaluation by the machine learning model gained binding to fHbp v1.1-specific antibodies, providing a proof of concept for further exploration of this approach.

Results

Preparatory experiments

The set-up phase focused on generating model-specific training data and confirming the feasibility of expanding the antigenicity of fHbp v3.28 to fHbp-specific antibodies. Three goals were considered in the design and selection of fHbp mutants in the set-up phase: (1) introduce substitutions into fHbp v3.28 or fHbp v1.1 to generate mutants with variable binding to two fHbp v1.1-specific monoclonal antibodies (mAbs) known as 12C1 and JAR5 (Figure 1), (2) select mutations in fHbp v3.28 that do not disrupt binding to cross-reactive mAbs 1A12, 4B3, and 1E6, and (3) determine the difficulty of introducing the 12C1 and JAR5 epitopes into fHbp v3.28 without disrupting binding to mAbs 1A12, 4B3, and 1E6. To accomplish these goals, residues within the 12C1 and JAR5 epitopes (as defined by their crystal structures with fHbp v1.1) (Malito et al., 2013; Malito et al., 2016) were identified and different combinations of these residues transferred from fHbp v1.1 into the fHbp v3.28 sequence (to test gain of binding) or vice-versa (to test loss of binding) (Figure 1). In locations where the amino acid residue is conserved between fHbp v1.1 and fHbp v3.28, a structurally compatible alternative residue for substitution at that location was selected by using StralSV, a computational tool that determines the probability of finding an amino acid given the local structure of the protein (Zemla et al., 2011).

Illustration of the mAb binding locations on the surface of fHbp. Aligned overlays of independent crystal structures, each containing a single Fab antibody fragment bound to fHbp v1.1 or fHbp v3.28; PDB: 2YPV, 5T5F, 5O14, 6H2Y, and 6XZW (Riahi et al., 2021; Desautels et al., 2024; Bianchi et al., 2019; Lopez-Sagaseta et al., 2018; Malito et al., 2013). (A) Bound Fab structures of the fHbp v1.1 specific mAbs are depicted in green and blue. (B) Bound Fab structures of the cross-reactive mAbs 4B3, 1A12, and 1E6 are depicted in orange, purple, and red respectively. Residues in the epitopes of fHbp v1.1-specific mAbs JAR5 (blue) and 12C1 (green) are indicated on the fHbp surface.

In total, 131 mutants of fHbp v1.1 and v3.28 were proposed. The binding affinity for mAbs 12C1, JAR5, 1A12, 4B3, and 1E6 to each mutant was estimated computationally. Specifically, 10 homology models were generated for each mutant and two approaches from the FoldX package, ‘BuildModel’ and ‘AnalyseComplex,’ were used to determine the change in binding free energy (ΔΔG) for the interaction of each mutant with the mAbs (Supplementary Table S1; see also Materials and Methods). Forty-eight mutants, eight fHbp v1.1 derivatives and forty fHbp v3.28 derivatives, were selected for experimental evaluation considering binding predictions for all five mAbs and the stated goals (Figure 2).

Round 1 summary of computational and experimental outcomes. (A) Interface ΔΔG calculated with FoldX for binding of fHbp v3.28-derived sequences to mAbs JAR5 and 12C1. ΔΔG values for binding to 1A12, 4B3, and 1E6 were also calculated (Supplementary Table S1). Of these sequences, and including v3.28, 41 were selected for production and are shown as red dots; sequences not selected for laboratory evaluation are shown as black dots. (B) For the selected sequences, experimentally validated binding is indicated with horizontal (pink, JAR5) or vertical (purple, 12C1) bars. Sequence m00006, which was selected as the starting sequence for Round 2, is at (−1.53, −2.11) in both panels.

In total, 50 fHbp proteins (48 mutants, plus fHbp v1.1 and fHbp v3.28) were produced in vitro on a small scale (4–12 mL) in Escherichia coli and evaluated for protein yield after one-step nickel affinity purification. From this group, 39 fHbp samples were successfully produced (final concentration >10 μg/mL), allowing the measurement of 195 antigen-antibody interactions (Supplementary Tables S2, S3) by biolayer interferometry (BLI). Affinity measurements (K_d_) fell between 750 pM and 5.5 nM for binding of all fHbp mutants to the cross-reactive mAb 4B3. This variation is expected given the accuracy of individual measurements by this technique. The consistent interaction with mAb 4B3, which recognizes a large discontinuous conformational epitope (Veggi et al., 2020), provided confidence that all fHbp mutants maintained the expected overall 3-dimensional fold. All mutants also bound the cross-reactive mAbs 1E6 and 1A12, though some mutants had an affinity reduced by up to 100-fold. Appreciable binding to mAb JAR5 was detected for 90% of the fHbp mutants produced, and 38% of the mutants also bound mAb12C1. Binding mutants were defined as having a measurable Octet K_d_ where the curve fit had an R ^2^ value of at least 0.8. The initial data for JAR5 and 12C1 binding to mutants includes 28 fHbp-mAb pairs with no measurable binding, four pairs with weak binding affinity (K_d_ > 50 nM), and 46 Ag-Ab pairs for which the binding of the fHbp mutant was similar to the Ag-Ab affinity observed for binding with fHbp v1.1. Mutants of fHbp v3.28 with the strongest binding to mAbs JAR5 and 12C1 (n = 2) had affinities equivalent to fHbp v1.1 (Table 1). Therefore, it was determined that the design problem selected for this project, gaining binding to mAbs JAR5 and 12C1 in fHbp v3.28 while maintaining binding to cross-reactive mAbs 1A12, 4B3, and 1E6, was achievable.

Finally, it was observed during these preliminary experiments that mouse 12C1 mAb performed poorly in the BLI affinity assay, with specifically low signal/noise ratio even when binding affinity was strong. Therefore, we humanized the mouse JAR5 and 12C1 mAbs and observed significant improvements in assay performance (Table 1). This was attributed to the substitution of anti-mouse mAb biosensors (AMC) for anti-human mAb biosensors (AHC). While data quality was improved, the Ag-Ab affinity measured was the same. The humanized JAR5 and 12C1 mAb reagents were used for all subsequent work.

Vaccine design pipeline integrating machine learning

The primary purpose of this work was to incorporate machine learning methods into our pipeline for vaccine design and selection. In the main experiment for this project, a machine learning method was implemented to select a diverse pool of fHbp mutant sequences after computational prediction of Ag-Ab binding. The design was focused on the 12C1 epitope (more difficult design problem than JAR5 epitope) and implemented a mutant design strategy that reduced the need for prior knowledge of fHbp v1.1. All mutants were made in fHbp m000006, a mutant of fHbp v3.28 identified in the preliminary work, that had good protein expression (88% yield when compared with wild-type fHbp v3.28). fHbp m000006 is primarily fHbp v3.28 with 51 residues mutated to their counterparts in fHbp v1.1. During the preliminary experiments, this mutant bound mAb JAR5 and did not bind mAb 12C1. The objective for the machine learning campaign was to computationally design and select mutants of fHbp m000006 that bound mAb 12C1 and improved binding to mAb 1A12, which shares an overlapping epitope with mAb 12C1.

Seven residue locations were selected for mutation, specifically in the region where the 12C1 and 1A12 epitopes overlap on the fHbp v3.28 surface, producing a design search space of 20^7^ (20 amino acids possible for each of seven locations). The search space was narrowed by selecting the top three structurally compatible amino acid substitutions using StralSV (Table 2) (Zemla et al., 2011). This resulted in a pool of 2,186 (3^7^-1) mutated sequences. This search space was probed with a hybrid computational approach as outlined conceptually in Figure 3. Calculations were performed using FoldX and STATIUM to calculate the ΔΔG for each mutant when binding all five mAbs (Supplementary Table S4). Machine learning was used to evaluate this data. Specifically, this work employed multi-objective Bayesian optimization and a GP model trained with published data from the AB-Bind database (Sirin et al., 2016) and project-specific data points collected in the preliminary experiments. Finally, a batch upper confidence bound decision rule (Desautels et al., 2014) was used to iteratively select 108 sequences, or 4.9% of the candidates for experimental evaluation (Figure 4). This decision rule ensured exploration (broad acquisition of information) with exploitation (selected mutants with the best predictors), while balancing the binding of fHbp mutants to mAbs 12C1 and 1A12 simultaneously. As constructed, the decision rule favored sequences with the best FoldX predictions, and the most important contribution of the GP model was to ensure diversity in the selected sequences for experimental evaluation. This stands in contrast to selection methods based on the output of individual computational predictors, such as FoldX alone, which might result in homogeneous, narrow groups of mutants (see Figure 5) that could fall victim to potential biases introduced by these tools.

Overview of the computational strategy for selection of diverse candidates with the GP model. All candidates within the defined search space are computationally modelled and multiple binding parameters are predicted using established biophysical tools, including STATIUM and FoldX. This data is combined in a machine learning model and the selection of candidates is driven by the training data provided. This allows the selection of a small set of candidates for experimental evaluation.

Sequences selected by the Gaussian process model for experimental evaluation. (A) Scatter plot of FoldX interface ΔΔG values (AnalyseComplex) for all 2,186 mutants evaluated computationally. The 108 mutants selected using the machine learning model and the decision rule are highlighted in red. Lines superimposed on the figure show effective selection thresholds in FoldX values corresponding to alternative sets of 108 mutants. If the 108 sequences with the best FoldX-predicted ΔΔG in binding 1A12 had been selected, points below the blue line fall in this set. Similarly, for if FoldX predicted binding of 12C1 is the only criterion, points left of the purple line are the resulting set, and for the equal sum of the two ΔΔGs for 1A12 and 12C1, points below and left of the orange line would constitute the resulting set. (B) The selected 108 points were tested in BLI. Among these, the best binders for 12C1 and 1A12 are marked, using cyan vertical and blue horizontal bars, respectively. A number of selected, strongly binding sequences do not lie in any of the sets constructed post hoc on the basis of FoldX binding predictions. (C) For each of group of sequences, sequence logos showing positions 221-223 and 249-252 are framed by the corresponding color. The sequence logo framed in black is all 2,187 mutants, including the parental sequence. (D) Pairwise distances among the 108 sequences in each selected or comparison set, using Blosum62. The black set is a size-matched, randomly selected set of 108 sequences from all 2,186 mutants considered, giving a fair comparison of intra-set sequence distances. Boxes show first quartile, median, and third quartile, while whiskers are 1.5 times the interquartile range or to the most extreme datum, whichever is narrower. Using the GP model, the decision rule traded off predicted binding performance for greater diversity in the selected set.

Alignment of the fHbp v3.28 mutant sequences with the highest affinity for mAb 12C1 and 1A12. The starting sequence selected for the machine learning design test, fHbp m000006, is shown in bold and deviations from this sequence are shown in green font. Residues identical in all sequences are in light grey. Locations allowed to mutate are highlighted.

Experimental evaluation of binding affinity was carried out for a selected batch of 102 candidate sequences, 83 of which could be produced in small scale E. coli cultures (8 mL) with sufficient yields (5%–50% yield compared with fHbp v3.28 wild type). Affinity measurements for the binding of 12C1, 1A12, and 1E6 were performed for all obtained mutants (Supplementary Table S5). Re-evaluation of the starting mutant, fHbp m000006, as a control in the main design campaign revealed unexpected weak binding of fHbp m000006 to mAb 12C1. This was attributed to improvements in the experimental assay. Therefore, final evaluation of the mutants focused on simultaneous improvement of binding to mAbs 12C1 and 1A12. All evaluated mutants bound mAb 1E6 with an affinity equivalent to the starting sequence (within 2-fold), or stronger, confirming that all fHbp mutants were folded. As expected, variation was observed among the binding of mutants to mAbs 12C1 and 1A12. In total, 40 mutants selected by machine learning had measurable affinity to mAb 12C1 and mAb 1A12. From this set, 17 fHbp mutants with the highest binding affinities were selected and it was confirmed that these samples also bound mAb 4B3 and mAb JAR5 with the expected affinity. Finally, six mutants were identified that bound all five antibodies with affinities similar to those obtained for fHbp v1.1 (affinity for mAbs 1A12 and 12C1 both stronger than 5 nM). When the respective sequences were aligned, they revealed a high degree of homology within the target mutation sites to fHbp v1.1 (Figure 5; Table 2) and a consensus sequence was derived (Table 2).

In total, 12 mutants, or 0.5% of the mutants in the decision set of 2,186 possibilities, conformed to this consensus sequence. The 108 mutants selected for experimental evaluation represented 4.9% of the total possibilities and included seven, or 58.3% of the potential fHbp candidates that matched the consensus sequence, a 12-fold enrichment. Only one experimentally evaluated mutant that conformed to this consensus showed weaker binding, namely, fHbp m002419, with an affinity to mAb 1A12 being 10 nM. Overall, the machine learning methods implemented here selected an enriched set of mutants for experimental testing, leading to multiple candidates capable of binding mAbs 12C1 and 1A12 with high affinity (stronger than 5 nM).

Biophysical characterization of a top mutant confirms protein integrity and antigenic design

The computational and experimental evaluation of fHbp mutants thus far focused solely on antigenicity (antibody binding) and did not confirm protein stability/integrity/folding beyond the inclusion of cross-reactive mAbs in the testing panel. The mutant from the ML design campaign with the highest measured affinities for 12C1 and 1A12, fHbp m002416, was produced in sufficient quantity to perform further biophysical characterization. Purified fHbp m002416 was monodisperse and monomeric as determined by size-exclusion chromatography (Supplementary Figure S1). The thermal stability of fHbp v1.1, fHbp v3.28, and fHbp m002416 was compared using differential scanning fluorimetry. As shown previously, fHbp v1.1 showed greater thermostability in the N-terminal domain (Tm1) relative to fHbp v3.28 (Malito et al., 2013; Johnson et al., 2012). By comparison, fHbp m002416 showed a reduction in thermostability, most significantly in the N-terminal domain, relative to both fHbp v1.1 and fHbp v3.28, suggesting a significant impact from the introduced mutations (Table 3).

To further explore the impact of mutations in fHbp m002416, a crystal structure was determined for this protein in complex with the Fab fragment of mAb JAR5 (Figure 6). JAR5 Fab was crucial for the crystallization of fHbp m002416, which did not crystallize in isolation. The JAR5 epitope/paratope interface adopts the same conformation as seen previously in the structure of JAR5 bound to fHbp v1.1 (Malito et al., 2016). The overall fold and domain structure reproduced the expected backbone for fHbp, including the 12C1 epitope where mutations were introduced. While >70% of the amino acid residues in fHbp m002416 were visible within the electron density, a section of the C-terminal domain could not be modelled due to flexibility and high local b-factors. When crystallized in space group C2, the solvent content of the crystal was 37%, the C-terminal domain of fHbp participated in few crystal contacts, and the disordered region was entirely exposed to the solvent channels. In contrast, crystals of fHbp v1.1 bound to JAR5 were previously grown in space groups C222_1_ (Malito et al., 2016) or I4_1_22, where the full C-terminal beta-barrel was modelled, despite above average b-factors relative to the N-terminal domain. Therefore, while fHbp m002416 did not maintain the thermostability of natural fHbp variants, it did maintain the overall fold of fHbp proteins while increasing the cross-reactivity to the tested panel of relevant antibodies. These results emphasized the importance of optimizing the computation-designed antigens for multiple parameters, including both antigenicity and stability, in future campaigns. Crystallographic statistics for new structures presented herein are summarized in Supplementary Table S6.

Crystal structure of fHbp m002416. (A) Crystal structure of fHbp m002416 (blue) crystallized with JAR5 (black) PDB: 2YPV. (B) Crystal structure of fHbp m002416 bound to JAR5 (black) with fHbp m002416 colored by b-factors. (C) Crystal structure of fHbp v1.1 in complex with JAR5 (black) in space group C2221, fHbp v1.1 colored by b-factor. (D) Crystal structure of fHbp v1.1 in complex with JAR5 (black) in space group I4122, fHbp colored by b-factor. B-factor scale from high (red) to low (blue) and standardized across panels (B–D).

Discussion

The type of interaction between antigen and antibody is distinct from general protein-protein interactions, where both binding partners have co-evolved over time (Akbar et al., 2021; Graves et al., 2020). Therefore, it is important to develop specific datasets and tools for predicting and optimizing antibody-antigen interactions. Protective epitopes can often be clustered close together (Neu et al., 2016) and few tools are available for optimizing the binding of antigens to multiple antibodies simultaneously. Machine learning methods employing GP-models have already been successfully applied to predict various protein properties, including thermostability (Jokinen et al., 2018; Pires et al., 2014; Romero et al., 2013), substrates for enzymatic reactions (Mellor et al., 2016), fluorescence (Saito et al., 2018), membrane localization, and peptide binding to MHC complexes (Ren et al., 2011). This work evaluated an integrated approach, incorporating both biophysical calculations of Ag-Ab affinity and machine learning strategies to optimize binding to multiple mAbs with overlapping binding sites.

The AB-bind database provided a limited number of datapoints (1,109) for generating a training dataset for immediate implementation of machine learning. Additionally, 45% of mutants in the AB-Bind database contain only alanine substitutions, and 37% represent a single point mutation to alanine (Sirin et al., 2016). Active ML strategies have shown small amounts of problem-specific data can significantly improve the performance of ML algorithms, especially when applied iteratively (Brochu et al., 2010; Liu, 2004; Shahriari et al., 2016). Therefore, preliminary experiments were designed to supplement the AB-bind data with fHbp-mAb datapoints, including a diversity of amino acid mutations and multi-residue mutations.

The primary design campaign focused on the 12C1 epitope and its overlap with the 1A12 epitope. While the preliminary experiments demonstrated that fHbp v3.28 can be mutated to bind fHbp v1.1-specific mAbs 12C1 and JAR5, the manual design strategy relied heavily on knowledge of the fHbp v1.1 sequence. By incorporating ML, the goal was to significantly reduce our reliance on information from fHbp v1.1. Additionally, the co-evolution of overlapping epitopes is a key point in vaccine design, making the interplay between 1A12 and 12C1 binding of specific interest.

All residue locations shared between the 12C1 and 1A12 epitopes were allowed to mutate, and the number of residues allowed in each position was down-selected. To provide an unbiased method for determining the search space, the algorithm StralSV (Zemla et al., 2011) was used to suggest the most compatible residues based on local structure. The top three residues for each location included known residues from fHbp v3.28, fHbp v1.1, and one to two additional residues. Thus, it was apparent that the known solution, all fHbp v1.1 residues in these locations, would be within the search space. After experimentally evaluating 108 computationally selected sequences, a consensus motif was determined based on the top six successful mutants. This sequence conflicted with the fHbp v1.1 sequence in only one position, residue 223, where a glycine or serine was preferred over glutamine. This substitution in the top candidates did not improve the binding of fHbp to either mAbs 12C1 or 1A12 over fHbp v1.1.

ML models provide a measure of the informational relationships and dependencies among a collection of many variant sequences. This allows the construction of batches of simultaneous experiments that strategically trade off exploitation of knowledge (prediction of strong performance) with the acquisition of valuable, non-redundant information on the overall design space. Balancing these objectives can result in practical gains in the quality of the selected sequences when taken as a group, including discovering effective designs not in the top tier of any ranking directly derived from computational tools (see Figure 5A), greater sequence diversity, and the balancing of competing design goals (Figures 5B,C). In this case, the search scope after down selection by StralSV was too narrow to allow the algorithm to propose a non-intuitive solution. Within the space allowed, the algorithm funneled toward the known fHbp v1.1 solution. To expand on this work in the future and find potentially novel solutions, the amount of down selection performed by StralSV will need to be reduced.

While antigenicity and binding affinity are key properties of antigens, developability and antigen stability are also crucial for developing an effective vaccine. Here, 18% of sequences selected by the ML algorithm were not successfully produced, reducing the efficiency of experimental evaluation. While the high-throughput production methods employed do not pin-point specific failures, the most likely issues include inadequate protein expression, protein aggregation, or reduced stability (Tokuriki and Tawfik, 2009). Given that the computational model only considers binding affinity, there was concern that stability and protein folding in successful designs had also been impacted. To further explore, the highest affinity binder, fHbp m002416, was selected for biophysical characterization. A crystal structure of this ML-derived mutant confirmed the expected 3D fold of the fHbp protein backbone. Interestingly, while this mutant folded correctly and had high affinity for both antibody targets, its thermostability was reduced. Therefore, an optimal ML strategy for vaccine antigen design must balance improved antigenicity with other protein properties, including protein expression and stability. This phenomenon has also been observed in the development of monoclonal antibody therapies and small molecule therapeutics using ML methods (Schneider et al., 2020).

Given the rapid advance of vaccine delivery platforms, including recent breakthroughs in using mRNA-encoded antigens, selecting and designing appropriate molecules is poised to become a rate-limiting step in vaccine development and future pandemic responses. The application of computational algorithms has massive potential to reduce the number of wet-lab experiments and the time required for this stage of development. Moving forward, an active ML approach that iterates between computational simulations and targeted wet-lab experimentation performed in a rapid high-throughput manner is expected to provide the most efficient tool for vaccine antigen development. Our work established the core capabilities required to develop this technology further and tackle more challenging design problems.

Materials and methods

Computational evaluation and selection of mutants for preparatory experiments

X-ray crystallographic structures of fHbp v1.1 bound to five specific antibodies, namely, 12C1, JAR5, 1A12, 1E6, 4B3, are available in the Protein Data Bank as 2YPV, 5T5F, 5O14, 6H2Y, and 6XZW (Bianchi et al., 2019; Lopez-Sagaseta et al., 2018; Malito et al., 2013; Malito et al., 2016; Veggi et al., 2020). All antibodies bind fHbp v1.1; antibodies 1A12, 1E6, and 4B3 also bind fHbp v3.28. Five homology models were created using the AS2TS system (Zemla et al., 2005) for fHbp v3.28 in complex with each antibody based on the available crystal structures. Models of fHbp v3.28 in complex with the antibody panel were used to define two putative epitopes each, including residues with at least one atom from fHbp within a distance ≤5 and/or ≤7 Å from the respective antibody (Supplementary Figure S2). Non-conserved amino acid locations within the proposed 12C1 and JAR5 epitopes were mutated to contain fHbp v1.1 residues. Overlapping positions were defined as those that are included in more than one epitope. Residues in overlapping regions were varied within the mutant panel to include residues present in fHbp v1.1, fHbp v3.28, or those suggested as conservative mutations by StralSV (Zemla et al., 2011). StralSV was used to select an alternative residue substitution at the positions where V1 and V3 had the same residue. The final mutants scored well according to ΔΔG assessment, and as anticipated showed binding to JAR5 experimentally confirmed. Sequences for all mutants evaluated are listed in Supplementary Document S1.

Five homology models were created using AS2TS (Zemla et al., 2005) for each sequence in the initial pool of candidates (5 × 131 = 655 total) to generate a pose in complex with each of the five antibodies. Each model was energy-minimized using molecular dynamics (MD) in GROMACS using the steepest descents method (Abraham et al., 2015). Both the energy-minimized, and original homology models were used in future calculations, for a total of 10 models per fHbp antigen sequence. Two approaches from the FoldX package were used to determine the change in binding free energy (ΔΔG); ‘BuildModel’ - “FOLDX total energies before and after introducing mutations” and ‘AnalyseComplex’ - “FOLDX energies in interfaces” (Schymkowitz et al., 2005). From this initial set of 131 mutants, 48 mutants were selected based on the best scores in all evaluated complexes and agreement in the two approaches for calculating ΔΔG values. All ΔΔG values calculated in this round are included in Supplementary Table S1.

Computational evaluation and selection of mutants using a GP model

From the training data, mutant m000006 was selected as the starting fHbp sequence. Seven residues from a region of the 12C1 epitope that overlaps with the 1A12 epitope were identified. StralSV was used to select the top three residue substitutions that would be compatible with these locations (Zemla et al., 2011). This yielded 2,186 (3^7, minus the unmutated m000006 sequence) sequences, and homology models were generated in AS2TS (Zemla et al., 2005) for each of them. FoldX (Schymkowitz et al., 2005) and STATIUM (DeBartolo et al., 2014) calculations were performed for each candidate to determine the change in binding free energy (ΔΔG) relatively to the starting sequence. This yielded 2,186 (3^7, minus the unmutated m000006 sequence) sequences (Supplementary Document S2) and homology models were generated in AS2TS (Zemla et al., 2005) for each of them. FoldX (Schymkowitz et al., 2005) and STATIUM (DeBartolo et al., 2014) calculations were performed for each candidate to determine the change in binding free energy (ΔΔG) relative to the starting sequence (Supplementary Table S4).

The relationships between likely experimental outcomes were modeled via a GP (Rasmussen et al., 2004). To represent mutant antigens to the GP model and make predictions about the effects of these mutations, a feature vector was constructed as follows: The contacts or potential contacts in these interfaces were identified by α-α distances of 10 Å or less. Each residue in such an interaction was assigned a chemical class and a size class (Supplementary Table S7); each identified interaction was then added to a tally of how many interactions exhibit particular chemical types (e.g., an aliphatic/aromatic interaction) and size types (e.g., a small/very large interaction). This set of all interaction features, for both the pre- and post-mutation versions of the antigen, was assembled into a feature vector. The outputs of the STATIUM and FoldX models were concatenated onto this feature vector and the whole feature vector was passed to the GP. The GP applies a normalizing feature transformation and Matérn kernel function to determine the similarity of any two such sequences’ mutational effects on the binding free energy. The model was implemented in Python in scikit-learn (Pedregosa et al., 2011). It was trained on the AB-Bind data set (Sirin et al., 2016) and antigen-specific results from round one experiments. For each interaction, a score is calculated, which is the model’s predicted mean, the model’s predicted variance, and a weighted combination of the FoldX and STATIUM scores. Because antigens must bind several target antibodies, separate predictions and separate scores were calculated for each of the target antibodies. The final score function was a sum over target antibodies:

$$[eqn]$$

Where $[eqn]$ was the antigen under consideration; $[eqn]$ and $[eqn]$ were the predicted (“posterior”) mean and standard deviation of the GP model for this antigen and the $[eqn]$ th target antibody, given selections 1 to t-1; $[eqn]$ was the set of simulation results (FoldX and STATIUM); and $[eqn]$ were weighting coefficients, respectively scalar-, scalar-, and vector-valued. We then selected the maximizer of $[eqn]$ , i.e.,

$$[eqn]$$

Where $[eqn]$ was the current decision set. Note that because desirable values for $[eqn]$ and $[eqn]$ were negative (i.e., mutation is desired to make binding free energy more negative), $[eqn]$ and $[eqn]$ were negative.

A batch upper confidence bound decision rule (Desautels et al., 2014) was adapted and applied to the posterior distribution to iteratively select the batch of 108 sequences from the decision set of 2,186. This decision rule balances exploration (broad acquisition of information) with exploitation (selecting sequences that are predicted to be promising in terms of free energy), while balancing the multiple target antibodies.

High-throughput cloning and sequencing

The fHbp mutant amino acid sequences were reverse-translated and codon-optimized for E. coli K12 expression using the Geneious R10 software (Version 10.1.3, Biomatters Ltd., San Diego, CA). Target sequences were synthesized as gBlocks (Integrated DNA technologies, IDT, Iowa, United States) and cloned into plasmid pBR322. Linearized plasmid was produced by PCR using Q5 Hot-start High-Fidelity DNA Polymerase (New England Biolabs) and site-specific primers (IDT), following the manufacturer’s protocol. Clones were assembled using NEBuilder HiFi DNA Assembly (E5520S, New England BioLabs) to perform Gibson Assembly in an automated 96-well plate format on a customized Tecan Fluent platform. This platform was used to resuspend dried DNA reagents, combine reaction components at the appropriate ratios, heat/cool the plate as required, transform competent bacteria, spread and incubate bacterial plates, and pick colonies for inoculation into LB media. New plasmid DNA constructs were isolated from bacterial culture with the MagJET Plasmid DNA Kit (ThermoFisher Scientific) automated in a 96-well format on a Kingfisher Flex (ThermoFisher Scientific). Sequence confirmation was performed by Sanger sequencing (Genewiz).

Production of antibody reagents

Plasmids for the expression of mAb 4B3, 1E6, 1A12, JAR5, and 12C1 have been described previously (Bianchi et al., 2019; Lopez-Sagaseta et al., 2018; Malito et al., 2013; Malito et al., 2016; Veggi et al., 2020). The human mAb 4B3, 1E6, and 1A12 and mouse antibodies JAR5 and 12C1 were prepared as previously reported (Veggi et al., 2020). In brief, plasmid DNA was prepared by maxi-prep and transiently transfected into Expi293 cell using an ExpiFectamine kit (ThermoFisher). Culture supernatant was harvest after a 4–5-day growth in a shake-flask culture at 37°C and 8% CO_2_. Antibodies were purified from the harvested supernatant by affinity chromatography (HiTrap MabSelect SuRe, Cytiva) followed by size-exclusion chromatography (Superdex 200 16/600, Cytiva). Antibodies were eluted from the affinity column with 0.1 M sodium citrate pH 3.0, which was immediately neutralized with 1 M Tris pH 9.0 (5:1 ratio). Buffer was exchanged into the storage buffer, 20 mM HEPES pH 7.0, 150 mM NaCl, during size-exclusion chromatography. The mouse antibodies had sub-optimal performance during initial binding studies (specifically low signal/noise ratio). Therefore, the variable domains from the mouse JAR5 and 12C1 mAb were subcloned into a human IgG framework containing the crystallizable domains (CH1, CH2, CH3, CL1). These chimeric antibodies, hJAR5 and h12C1, were produced using the same protocol and used for evaluating the mutants selected by the GP model.

Fab fragments for mouse JAR5 were sub-cloned from the original mouse vectors and produced in Expi293 cells using the same overall protocol as the mAb reagents. Culture supernatant was concentrated and diafiltered into Tris saline buffer (50 mM Tris pH 8.0, 150 mM NaCl) prior to purification. Fab fragments contained a Streptavidin tag fused to the C-terminus of the heavy chain which enabled purification by affinity chromatography (StrepTrap, Cytiva). Protein was eluted in 100 mM Tris pH 8.0, 150 mM NaCl, 1 mM EDTA, 2.5 mM desthiobiotin and further purified by size-exclusion chromatography (Superdex 200 Increase 10/300 GL, Cytiva) in 25 mM Tris pH 8.0, 150 mM NaCl. All antibody samples were concentrated using Millipore centrifugal filters. Antibody reagents were stored frozen at −80°C prior to characterization.

Production of fHbp reagents

Plasmids for the expression of wild-type fHbp variants 1.1 and 3.28 with a C-terminal 6-His fusion tag were reported previously (Masignani et al., 2003; Giuliani et al., 2005). Protein sample references and protein samples for X-ray crystallography were produced in low throughput on a large scale. Chemically competent E. coli BL21(DE3) star cells (ThermoFisher) were transformed with the respective fHbp plasmid and protein was expressed during growth in up to 1 L of Terrific Broth (TB) media containing 100 μg/mL carbenicillin. Cultures were incubated at 37°C for 16 h prior to harvest. Cells were lysed by sonication and insoluble debris was removed by centrifugation. Soluble protein was purified by nickel-affinity chromatography (HisTrap FF, Cytiva). Protein was eluted with a linear gradient of 20–250 mM imidazole in 50 mM Tris pH 8.0, 300 mM NaCl. Further purification was performed by size-exclusion chromatography (Superdex 200 Increase 10/300 GL, Cytiva) using 25 mM Tris-HCL pH 8.0 and 150 mM NaCl as the elution buffer. Protein samples were concentrated (Millipore-Sigma) and protein concentration was determined using optical density at 280 nm before storage at −80°C.

Mutant fHbp sequences were produced in parallel using a high-throughput approach. Protein was expressed using 4 mL cultures of TB media with 100 μg/mL carbenicillin in 24 deep-well plates shaken at 37°C. Harvested cells were lysed with BugBuster Master Mix reagent (Millipore Sigma). Mutant fHbp samples were purified using one-step nickel-affinity chromatography using the MagneHis protein purification system (Promega) and automated with a Kingfisher Flex (ThermoFisher). Proteins were eluted in 50 mM Tris pH 8.0, 300 mM NaCl, 500 mM imidazole and protein quality were assessed by SDS-PAGE.

Biolayer interferometry experiments (BLI)

The concentration of fHbp samples was determined using an Octet 384 Red Instrument (Sartorius FortéBio) operating at 30°C and anti-His biosensors. Highly pure fHbp v3.28 was used as the standard curve protein and all protein samples were diluted 1/10 in Octet buffer (1 x phosphate-buffered saline [PBS], 1% w/v bovine serum albumin). Binding kinetics were also measured using BLI at 30°C and anti-Human Fc biosensors. IgG antibodies at a concentration of 10 μg/mL were captured on the biosensor followed by binding of fHbp samples with concentrations ranging from approximately 150 nM–1 nM. For each fHbp-mAb pair, binding traces (average n = 4) were fit using the FortéBio Analysis 11.1 HT software package to determine K_D_, K_on_, and k_off_.

Differential scanning fluorimetry

Previous studies based on differential scanning calorimetry have demonstrated that fHbp proteins from variant groups 1, 2 and 3 all have two distinct melting points that correspond to the unfolding of the N- and C-terminal domains (Malito et al., 2013; Johnson et al., 2012). To measure both events, while using a minimal amount of protein, the stability of fHbp m002416 was determined using differential scanning fluorimetry. The first melting point (Tm1) was measured using the hydrophobic dye SYPRO orange and corresponds to the unfolding of the N-terminal domain. The second melting point (Tm2) was measured using intrinsic fluorescence and corresponds to the unfolding of the C-terminal domain.

fHbp samples were diluted in PBS to an approximate concentration of 0.5 mg/mL and a total volume of 20 μL. Samples were further diluted with 400 μL of buffer followed by concentration to 50 μL using 0.5 mL, 10 kDa cutoff devices from Millipore. This process was repeated twice to remove imidazole and other buffer components remaining from the purification. Final sample volume was approximately 80 μL. SYPRO Orange hydrophobic dye (ThermoFisher) was diluted from 1000X to 10X in PBS. Chemical DSF samples were prepared in duplicate by mixing 18 mL of protein with 2 mL of dye in a 96-well PCR plate and protecting from light. Samples were analyze using the ViiA 7 System (ThermoFisher) using a ramp speed of 0.02°C/s from 25°C to 99°C. Melting scans were analyzed using Protein Thermal Shift v1.4 and analyzed with default parameters. In addition, fHbp protein samples in PBS were loaded into high-sensitivity capillaries and analyzed using a Prometheus NT.48 system (nanoTemper). The melting scan was performed from 24°C to 110°C with a ramp of 1°C/min. Duplicate measurements were performed for each sample and melting temperature was determined using the manufacturer provided software.

X-ray crystallography

Purified fHbp 2,416 and JAR5 Fab were combined in a 1:1.1 ratio at room temperature in protein buffer (25 mM Tris pH 8.0, 150 mM NaCl) at a final protein concentration of 11.3 mg/mL. Sparse matrix screening was performed using a Mosquito automated drop setter (SPT Labtech) and drop volumes of 200 nL protein +200 nL precipitant. Crystals were identified in well B8 of the ProPlex screen (Molecular Dimensions); 0.1 M sodium citrate pH 5.0, 15% w/v PEG 4000, 0.1 M magnesium chloride and were harvested after an 8-day incubation at 20°C. 20% glycerol was added as cryoprotectant prior to plunge freezing in liquid nitrogen. Data were collected at the Advance Photon Source, beamline 22-ID equipped with an Eiger 16M detector. The X-ray wavelength was 1.0 Å and 800 images were collected with a frame width of 0.25° at a detector distance of 300 nM. Diffraction patterns were integrated and scaled using XDS (Kabsch, 2010). Data was truncated in XDS to remove frames from the end of the data collection that showed the most radiation damage. Molecular replacement was performed by placing individual models of fHbp v1.1 and JAR5 Fab extracted from PDB 5T5F using Phaser as implemented in Phenix (Adams et al., 2010; McC et al., 2007). The final model was built after iterative rounds of refinement in Phenix and manual corrections in Coot (Emsley et al., 2010). Visualization and analysis of the final model were performed with PyMOL 2.0 (Schrodinger LLC).