Model

To address the challenging task of protein design, we performed thorough modelling to select the best candidates to be tested in the lab.

The first step in the process of engineering an intrinsic factor (IF) that can evade the autoantibodies is to identify the epitope on the IF and engineer an IF that can evade the clearance mediated by autoantibodies. Based on previous knowledge [1], the residues in the alpha helix region adjacent to the cobalamin (vitamin B12) binding pocket: 256Ser, 258Lys, 262Tyr and 265Val were to be mutated. Optimal mutations to substitute the residues with other or “similar” residues were decided and engineered in silico.

Then, we designed our mutants based on the rationale of preserving structure as much as possible. This principle was chosen as the antibody binding region is in close proximity to the B12 binding region. Thus, mutagenesis should prevent IF from binding to autoimmune antibodies while preserving ligand-receptor binding affinity. For this purpose, we used BLOSUM62 matrix [2] to detect closely related amino acids. These are the ones whose substitution is well tolerated in evolution.

Mutant creation

Previous Next

We performed protein structure predictions to verify that the designed mutants would preserve the properties of the wild type intrinsic factor. We utilised a well-known tool called AlphaFold2 [3] and used it as ColabFold [4]. Structures of all 21 mutants were obtained as PDB files, visualised and structurally aligned to the wild type IF crystallography structure (PDB ID: 2pmv) using PyMol [5]. This enabled us to refine the list of tested mutants excluding the improperly folded ones.

Docking and MD

Docking

The final step of the modelling consisted of screening for potent therapeutic candidates. We performed docking studies to ensure that the binding affinity of engineered protein is at least on par with the wild type IF. Docking studies were performed in 2 steps. First, we predicted the binding energy. AutoDock Vina [6] tool was selected for this task as it is a well-established tool for docking simulations. Energy refinement was performed using NNscore 2.0 neural network [7]. The NNScore 2.0 scoring function was used to estimate the binding energy and dissociation constant, as it is more accurate than the AutoDock Vina scoring function.

In addition, binding energies of the B12-IF complex with the cubilin receptor were assessed using MutaBind2 [8] since the strength of binding is a crucial factor for the B12-IF complex internalisation. However, we then found evidence that MutaBind2 has poor performance [9].

The most promising candidates were selected to be tested in wet lab experiments. We found that most proposed mutants have similar or better binding affinity to vitamin B12. Therefore, these results indicate that our mutant selection pipeline is adequate and applicable to the research question.

Molecular Dynamics

To have a closer look at the interactions at molecular level, attempts were made to model the IF-B12, protein-ligand complex for molecular dynamics study using GROMACS [10] software. Some of our initial tasks included pre-processing of the available crystallography structure data and topology generation. As the available structure is not fully resolved for the entire protein, a tool called Swiss-PDB viewer [11] was used to fill in the missing atoms (predict structure) and the structure was cleaned of glycosylation and other unwanted atoms were removed using Chimera [12]. For MD simulations, our force field of choice was CHARMM36 [13] which is well established for protein studies and widely accepted in the scientific community. However, a ligand such as B12 is not a recognized entity in the force field. Thus, it was important to process the topology files separately (for IF and B12) and it served as the bottleneck of our experiment.

One of the most challenging tasks in molecular simulation is the proper processing of ligands. Introducing new entities into the framework of the force field requires careful consideration and validation, often in the form of various quantum mechanical calculations [14]. This could be a very laboursome and tedious process. For preparing the ligand topology, an automated tool CGenFF [15] was used. The input file for the tool was generated using another software Avogadro [16]. An advantage of using CGenFF is that it returns penalty scores for each parameter and our penalty scores were larger than 50, which generally require manual reparametrization. The scores were exceptionally high around the corrin ring and metal ion - cobalt.

In the end, after trial and error, we were not able to implement successful parameterization which would accurately reflect a realistic model, due to lack of prior knowledge and time constraints.

Immunogenicity

Since we are modifying the protein, there is a high chance that artificial modifications could trigger an additional immunogenic response. Thus, we decided to screen our candidates for immunogenicity to ensure that the produced mutants do not trigger immunogenic responses themselves. Our focus lies primarily on possible responses from B lymphocytes [17].

We performed screens for continuous as well as discontinuous epitopes for possible antibodies of an immunogenetic response. The difference between them is that the discontinuous takes into consideration the folding while the linear only takes the ordered amino acids into account [18].

After performing an extensive literature review, we proceeded with the online tool IEDB (Immune Epitope Database and Analysis Resource) [19], more specifically the tools ElliPro for discontinuous and Antigen Sequence Properties for the linear epitopes. ElliPro employs a clustering algorithm combined with Thornton's method (focusing on protruding amino acid residues) for making predictions [20]. The tool Antigen Sequence Properties utilizes several methods that consider beta turns, surface accessibility probability, mobility of protein segments, physicochemical properties of residues and hydrophilicity. Furthermore, they apply a hidden Markov model, the propensity scale method, and a random forest algorithm.

We assumed that introduced mutations would lead to a detectable drop in immunogenicity. However, this hypothesis turned out to be invalid both for the linear epitope predictions and the discontinuous ones. We utilized the area under the per residue immunogenicity score curve (AUC) as a single-value metric to make immunogenicity prediction results comparable between mutants. For the discontinuous predictions, the wild type reached AUC of 202.1 and the mutant with the lowest immunogenicity (a lower score correlates to a lowered predicted chance of autoimmune response) at AUC 201.8. To conclude, the mutants were not predicted to have significantly altered immunogenicity. This could arise due to several reasons:

1. There is no significant difference in autoantibody binding meaning that the previously described epitope is not the actual one.
2. The previously described epitope is the actual one. However, our rationale chosen for the mutant design is not applicable to this case. Therefore, mutations with a stronger effect on protein structure should have been selected.
3. Single point mutations introduce very minor changes in structure that do not allow selected tools to detect any significant difference between them. However, it could still be enough to disrupt antibody binding in vitro.

References