Distinct binding mode of BAFF antagonist antibodies belimumab and tabalumab, analyzed by computer simulation

B cell-activating factor (BAFF) can bind with specific receptors to activate signalling pathways associated with the B cell activation. Belimumab and tabalumab are anti-BAFF (B cell depleting) monoclonal antibodies, with therapeutic efficacy demonstrated for the treatment of autoimmune disorders, while belimumab was approved by FDA in 2011 as a targeted therapy for systemic lupus erythematosus (SLE) and exhibited better clinical outcome than tabalumab. In this investigation, the combination modes of BAFF-belimumab and BAFF-tabalumab complexes were simulated in silico to better understand the reason for the comparative inhibitory difference between belimumab and tabalumab. The structures of belimumab and tabalumab were constructed through homology modelling. The combination mode of BAFF-belimumab complex was analyzed by molecular dynamics simulation, while that of BAFF-tabalumab complex was analyzed by protein-protein docking following the molecular dynamics simulation. Both belimumab and tabalumab were bound with BAFF at the same hydrophobic center to which the natural receptors of BAFF bind as well. Belimumab heavy chain components I51, F54, K58, D100, D101, L102, L103, and P105 and R27, Y30, K49, and S65 of belimumab light chain contribute to the BAFF-belimumab interaction mainly via hydrogen bonds, salt bridges, and hydrophobic interactions. More importantly, belimumab could bind to L83 of BAFF and produce steric hindrance with the adjacent BAFF trimers, while tabalumab could not. Therefore, our results indicated that belimumab has a better clinical outcome compared with tabalumab mainly because belimumab could bind to L83 of BAFF and interfere the formation of a BAFF 60-mer, besides mediating inhibition of the interaction of BAFF with its receptors.

BAFF is a type of membrane-bound protein, but released in a soluble form by the action of a metalloprotease [27][28][29]. The soluble BAFF (sBAFF) is formed from the hydrolysis of membrane-type BAFF and can improve the activity of B cells, CD4+ T cells, and NK cells to ramp up the immune responses [30]. Three kinds of BAFF receptor have been reported; they are B cell maturation antigen (BCMA), BAFF receptor 3 (BR3) along with the trans-membrane activator and calcium modulator, and cyclophilin ligand interactor (TACI) [2,[31][32][33][34][35][36][37][38][39][40][41]. Although BAFF often functions in the homo-trimeric form, a more active BAFF 60-mer has also been reported [42]. Twenty soluble BAFF trimers oligomerize into a 60-mer through the trimer-trimer interaction via a long DE loop (KVHVFGDELS), and this 60-mer is considered to be more active than the trimer due to the clustered accumulation of these receptors on the B cell surface [43,44]. Since BAFF can bind to the B lymphocytes, it can subsequently induce their proliferation, differentiation, and secretion of immunoglobulins [45,46]. Some researchers have demonstrated that the elevated levels of BAFF are closely related to SLE pathogenicity [26]. Therefore, BAFF and its downstream signal transduction factors have become therapeutic targets for the treatment of many autoimmune disorders.
Antibodies targeting BAFF can inhibit the proliferation and activation of B cells by blocking the binding of BAFF and its receptors [47]. Anti-BAFF monoclonal antibodies can reduce the elevated levels of immunoglobulin and lupus-related damage factors in the peripheral blood of SLE patients [48]. Belimumab (tradename Benlysta) is a fully human monoclonal IgG1λ antibody, which can neutralize the soluble BAFF [12,[49][50][51][52][53][54][55]. It was approved by FDA in 2011 as a targeted therapy for SLE [56][57][58]. Clinical trials of SLE patients with belimumab exhibit reduced number of circulating naive B cells, activated B cells, and enhanced clinical efficacy [48,59,60]. Tabalumab (formerly LY2127399) is a fully human IgG4 monoclonal antibody. Similar to the natural BAFF receptors, tabalumab binds both membranebound and soluble BAFF, while belimumab binds only the soluble fraction of BAFF [61,62]. Tabalumab improved the signs and symptoms in patients of rheumatoid arthritis (RA) who had previously become naive to biological diseasemodifying antirheumatic drugs (bDMARDs) such as TNF inhibitors in an initial dose range study [63]. Telitacicept is a fusion protein constructed from the extracellular part of TACI and Fc part of human IgG1. TACI has a strong affinity for BAFF and a proliferation-inducing ligand (APRIL), which is a type II membrane protein that can stimulate B cell proliferation and tumor cell growth [64,65]. Therefore, telitacicept can block the interaction of BAFF and APRIL with their respective receptors to block their biological activity. This bi-specific blocking effect may be more effective than blocking BAFF alone, consequently inhibiting the proliferation of B lymphocytes and treat autoimmune diseases more effectively. On March 12, 2021, telitacicept was approved by China's National Medical Products Administration (NMPA).
In this study, we employed homology modelling, molecular docking, and molecular dynamics (MD) simulation to develop a three-dimensional (3D) structural model of BAFF-belimumab complex, to understand the mechanism of interaction of belimumab in atomic detail. After that, we calculated the binding free energy and energy of decomposition of each residue to figure out the binding mechanism and alanine mutational analyses to validate our findings. Then, we built a 3D structure model of tabalumab and performed MD simulation for BAFF-tabalumab complex, to find out a possible reason why tabalumab failed in the clinical trials [66]. These results provide an insight into the mechanism of interaction between BAFF and its antibody antagonists and further provide a foundation for designing novel BAFFtargeted inhibitors in future.

Construction of molecular systems for BAFF and belimumab
The 3D crystal structure of BAFF in complex with belimumab was downloaded from the website of Protein Data Bank (http:// www. rcsb. org/) with an access code of 5Y9J. Since the continuous peptide structure was needed in the subsequent MD simulation, and there was a breakpoint in the downloaded 3D structure (S136 to G142 in heavy chain), the complete and continuous structure of belimumab heavy chain was obtained through homologous construction by importing its primary sequence into the SWISS-MODEL (https:// swiss model. expasy. org/) and using the downloaded structure as a template. To build BAFF-belimumab complex, the constructed heavy chain of belimumab was matched to the downloaded complex by MatchMaker of Chimera Ver1.14.

Energy minimization and MD simulation of BAFF-belimumab complex
During the molecular dynamics simulation and after building the starting structure, a simulation was carried out in parallel to obtain an equilibrated system. The SANDER module of the AMBER 16 package along with the AMBER parm99SB force field was used [67]. A chloride ion was added to neutralize the charge for the system. The system was soaked in the TIP3P water box with the truncated octahedron periodic box having a size of 1 nm [68]. Generally, the Particle Mesh Ewald (PME) method was used to deal with the long-range electrostatic interactions [69,70]. The SHAKE algorithm which supports atomic vibration in calculation simulation system was applied to constrain hydrogens. In the simulation, 0.8 nm was set as the cut-off distance, if/ when nonbonded interactions were encountered [71].
The solvated system which performed energy minimization was conducted with the help of MD simulation. The water-solvent molecules and ions were optimized first, while the atoms of the residues of the complex were restrained.
Then, for all residues, the main chain residues were kept frozen while relaxing the side chains. Finally, all atoms were allowed to move freely while 2 fs was set up for each integration step and the ensemble type was set to NPT (a kind of isothermal isobaric system). The steepest descent method was executed to go on the energy minimization for the first 5000 steps. Then, the conjugated gradient method was used for the next and subsequent 2500 steps. The temperature range was set from 0 to 300 K within the period of 100 ps after the structure was minimized to continue with the MD simulation. The systems were equilibrated in the NPT ensemble (1 atm and 300 K) for 100 ps when the proteins were kept frozen. Next, the complexes were equilibrated for 500 ps without restraint. At last, the whole system was subjected to 200 ns of molecular dynamics under the abovementioned NPT ensemble. To analyze the structure in detail, the PTRAJ module was used to collect the atomic coordinates of the system.

Binding free energy calculation on BAFF-belimumab system, per-residue free energy decomposition analysis and hydrogen bond analysis
After MD simulation, the interactions of the protein complex system were investigated. They were validated by calculating the binding free energy and decomposing the binding free energy into per-residue contributions. The molecular mechanics MM-GBSA algorithms in AMBER 16 were implemented to calculate the binding free energies between BAFF and belimumab. The total binding free energies in condensed phase can be calculated by the following formula: where ∆G bind is the binding free energy in solution state, ∆E MM is the molecular mechanics energy which is comprises an electrostatic interaction (∆E ele ), a van der Waals interaction (∆E vdw ) and an internal energy (∆E int ) for a ligand or a protein, whereas ∆G sol is the solvation free energy, comprising a polar component (∆G polar ) and a nonpolar one (∆G NP ). The calculation for ∆G polar can be done through GB method [72]. The nonpolar component ∆G NP is calculated by the solvent-accessibility surface area (SASA) equation ∆G NP = γSASA + β, where SASA can be obtained from the MSMS program, and the values of γ and β were set as 0.72 kcal/(mol*nm 2 ) and 0.00 kcal/mol, respectively [73]. Previous literature reports that entropy does not contribute much to the binding affinities of similar ligands [74].
The free binding energy decomposition of residue inhibitor interaction for each of BAFF and belimumab residues can also be calculated by alanine mutation. The free energy changes (∆∆G bind ) were estimated by comparing the delta free energy of the wild type (∆G WT ) with the delta free energy of the mutant into alanine (∆G mut ), ∆∆G bind = ∆G WT − ∆G mut , while the negative ∆∆G bind indicates that the wild-type protein is more favorable as the ∆G WT is more negative.
The hydrogen bond analysis was carried out by "hbond," which is a command of cpPTRAJ and allows hydrogen bond calculation using geometric criteria. Automatic determination of hydrogen bond donors/acceptors used the simplistic criterion that "hydrogen bonds are FON," which is, hydrogens bonded to F, O, and H atoms are considered donors, while F, O, and N atoms are considered acceptors.

Construction of BAFF-tabalumab complex
Due to the absence of a 3D structure of tabalumab, we constructed the 3D structure of tabalumab through SWISS-MODEL using 4NPY for heavy chain (86.04% homology) and 5WT9 for the light chain (97.98% homology) as a template which was recommended by the search template results. The constructed 3D structure of tabalumab was used for running the MD simulation using the same parameters as BAFF-belimumab complex to attain stable structure with low energy. Finally, the interaction between BAFF and tabalumab was analyzed and information used for making schematic diagrams via Chimera (Ver 1.14).

Homology construction of BAFF-belimumab system
Homologous construction was done by using the SWISS-MODEL to smoothen the remedy of the break point. To rank the constructed structure, Ramachandran plots were calculated to assess the geometric properties of the backbone conformations. This step indicated that all (100%) of the residues fell within the most favored region for both the heavy and light chains of belimumab. The constructed system was further refined by performing an optimized geometric calculation of the mechanics using the SANDER module.

MD simulation and free binding energy calculation
To ensure the proper conformational space sampling, we performed a 200 ns MD simulation for the BAFF-belimumab complex. The atomic root mean square deviation (RMSD) fluctuations of the backbone atoms (CA, C, and N atoms) of BAFFs and belimumabs were separately calculated and outlined into a function of time (Fig. 1).
It can be concluded that the BAFF and two chains of belimumab have relatively smooth curves in the last 20 ns period with the average RMSD values of 0.200 nm for heavy chain, 0.127 nm for light chain of belimumab, and 0.135 nm for BAFF. The system becomes equilibrated as judged the by RMSDs.
The individual energy terms of the binding free energies of the BAFF-belimumab complex were calculated by the MM-GBSA method for the last 20 ns data which was relatively stable according to the RMSD value. Different terms of energy values are listed in Table 1. The enthalpy of the BAFF-belimumab complex was estimated at − 48.68 kcal/ mol by the MM-GBSA approach.
In order to identify which of the interactions have a significant impact on the system, the binding free energy was divided into polar and nonpolar interactions. It was clear that the nonpolar interaction function was the dominant force, and the polar energy provided an unfavorable contribution to the formation of BAFF-belimumab complex by comparing values of the two types of energy. The absolute value of the nonpolar energy is about twice as much as the polar one. In addition, the polar interaction in the gas phase provided a favorable energy component, but it is mainly unfavorable in the solution phase for belimumab-BAFF complex.
In order to confirm the binding ability of belimumab and BAFF, the ∆G bind of belimumab-BAFF complex and the complex of BAFF and a natural receptor B cell maturation antigen (BCMA) were compared ( Table 2). The MD simulation of BCMA-belimumab complex was also carried out by the same method, and the individual energy terms of the binding free energies of BCMA-belimumab complex were obtained by the last 10-ns simulation structures, which were relatively stable according to the RMSD value to compare with belimumab-BAFF complex. It showed that nonpolar forces contribute much more to the complex binding for both the systems; this is consistent with the interaction between BAFF and its other natural receptors we reported earlier [71, Fig. 1 Time dependence of RMSD values along the MD trajectories for the complex of BAFF and belimumab. The heavy chain of belimumab is shown in red, while the light chain of belimumab is shown in blue, and BAFF is shown in black Table 1 The binding free energies of the BAFFbelimumab complex (kcal/mol) *Abbreviations of the components in expanded form: E vdw the van der Waals contribution from molecular mechanics; E ele the electrostatic energy as calculated by the molecular mechanics force field; G polar the electrostatic energy to the solvation free energy; G NP the energy of nonpolar components contributing to the solvation free energy; E MM the free energy of the complex in vacuum which equals to E vdw , E ele , and E int ; G sol the free energy of the complex in solvation states which is equal to G polar and G NP ; total, the whole binding free energy including E MM and G sol ; delta, the sum of BAFF and belimumab minus complex  75]. Compared with the natural receptor BCMA, the binding energy of belimumab with BAFF is higher in absolute value, which means that belimumab-BAFF complex interacts more strongly and binds more closely with each other than the BCMA-BAFF complex. This is the reason why belimumab can inhibit the binding of BAFF and their natural receptors. Thus, the function of the belimumab binding with BAFF to inhibit the function of BAFF can be realized.

The binding mode between BAFF and belimumab
To locate the binding site and investigate the interaction mechanism, we analyzed the binding mode between BAFF and belimumab (Fig. 2). BAFF can interact with both the heavy and light chains of belimumab primarily through hydrophobic interactions, hydrogen bonds, and salt bridges. The major binding site was found to consist of eight residues (I51, F54, K58, D100, L101, L102, L103, and P105) in the heavy chain of belimumab. The hydrophobic part of belimumab heavy chain was in close contact with Y65, L70, I92, and P123 of BAFF while heavy K58/ BAFFE 125, and heavy D100/ BAFF R124 formed the salt bridge (Fig. 2a). BAFF formed three hydrogen bonds with the light chain of belimumab; they were formed by the backbone of light R27 with backbone of BAFF S84, hydroxyl of light S65 with backbone of BAFF G80, and backbone of light G93 with backbone of BAFF S21, whereas K49 of belimumab light chain and D81 of BAFF formed the salt bridge, and the hydrophobic effect of the light chain and BAFF is mainly due to the Y30 of light chain and L83 and L85 of BAFF (Fig. 2b).

The binding free energy decomposition of per residue and computational alanine scanning
To further analyze the binding mechanism, binding free energy decomposition was performed on a per-residue contribution spectrum. The criterion of |∆G bind | >1 kcal/ mol was employed to identify the important residues that contribute to the free energy. According to this standard, seven residues in heavy chain of belimumab (I51, F54, D100, D101, L102, L103, and P105), four residues in light chain of belimumab (R27, Y30, K49, and S65), and eight residues in BAFF (Y65, L70, G80, L83, L85, I92, P123, and R124) were selected. Total free energy of each key residue was resolved and shown in Fig. 3. For belimumab, D100 in heavy chain and K49 in light chain mainly provided electrostatic interaction, while the side chains of L101, L102, and L103 in heavy chain formed a hydrophobic site and mainly contributed for the van der Waals interactions. These results of these interactions are consistent with our previous docking model; however, K58 of belimumab heavy chain and E125 was found on the surface of the protein. As for the desolvation step; it led to a decrease in the absolute value of binding energy although they (K58 and E125) formed a salt bridge. The ∆G bind value of heavy K58 and BAFF E125 was − 0.90 kcal/mol and 1.59 kcal/mol, respectively. And the same situation existed in D81 of BAFF with the ∆G bind value  of 0.19 kcal/mol. The hydrogen bond interaction between G94 of belimumab light chain and G20 of BAFF was weak and unstable in MD simulation, which can be concluded by their ∆G bind value of − 0.80 kcal/mol and − 0.82 kcal/mol, respectively. There was a hydrogen bond between R27 of belimumab light chain and S84 of BAFF, but the binding energy of BAFF S84 was only − 0.16 kcal/mol, which is far from meeting the set criterion of |∆G bind | >1 kcal/mol. According to the results of energy decomposition, the skeleton energy decomposition was − 1.07 kcal/mol, but the side chain energy decomposition was 0.91 kcal/mol. It indicated that the adverse effects caused by the side chain led to the failure of BAFF S84 to meet the criterion and cannot be considered as a key residue for belimumab-BAFF complex, although there was the existence of a hydrogen bond.
According to the results of hydrogen bond analysis (Table 3), the salt bridge formed by D100 of belimumab heavy chain and R124 of BAFF was very stable, while the salt bridges formed by belimumab heavy chain K58-BAFF E125 and belimumab light chain K49-BAFF D81 were not stable, with the occupancy only at 12.39% and 3.97%, respectively. The hydrogen bonds that were formed separately by R27 and S65 of belimumab light chain with S84 and G80 of BAFF were relatively stable, with high hydrogen bond occupancy, but it was not stable for the hydrogen bond between G94 of belimumab light chain and G20 of BAFF. These were consistent with the results of energy decomposition.
In order to confirm whether the side chain of specific residues play an important role in the binding of BAFF-belimumab complex or not, seventeen residues were selected to perform the computational alanine scanning ( Table 4). P105 of belimumab heavy chain, G80 and P123 of BAFF were not mutated into alanine because the skeletal structure changed significantly when these side chains were replaced. In addition to R27 of belimumab light chain and L85 of BAFF, the binding energies of other residues became less negative after computational alanine scanning. It indicated that if any one of these amino acids mutated to alanine, they would weaken the binding affinity of BAFF-belimumab complex. Among

Comparison of binding modes between BAFF-tabalumab, BAFF-belimumab, and BAFF-BR3 complexes
Tabalumab could interact with BAFF mainly by heavy chain through hydrophobic interactions, hydrogen bonds, and salt bridges. BAFF formed two hydrogen bonds with the heavy chain of tabalumab; they were formed by heavy N58 with BAFF E125 and heavy Y100 with backbone of BAFF R124. In addition, DXL motif (D101, I102, L103) of tabalumab heavy chain is involved in binding with BAFF, heavy D101 and BAFF R90 formed salt bridge. The hydrophobic effect of tabalumab and BAFF was primarily due to I102, L103 of tabalumab heavy chain and L70, V86 of BAFF (Fig. 4).
As shown in Fig. 5, the hydrophobic center of BAFF was the active site which could bind with its natural receptors such as BCMA and BR3 at the DxL motif [76], and a structure similar to the DxL motif could also be found in CDR3 of heavy chains in both belimumab and tabalumab (Fig. 6). As a result, the mechanisms of inhibition of BAFF by belimumab and tabalumab were considered to be "competitive" with the natural receptors.  The comparison of protein-protein docking schematic between BAFF-tabalumab complex and BAFF-belimumab complex for a the section included in the red square was the steric hindrance between the light chain of belimumab and the adjacent BAFF trimer after the combination of belimumab and BAFF, b the section included in the red square was the DxL motif where belimumab, tabalumab, and BR3 bound and acted with BAFF, c the interaction between Y30 of beli-mumab light chain and L83 of BAFF; tan, local structure of BAFF60mer; sky blue, heavy chain of belimumab; plum, light chain of belimumab; light green, tabalumab; red, the nature receptor of BAFF (BR3); the section included in the red square was the steric hindrance between the light chain of belimumab and the adjacent BAFF trimer after the combination of belimumab and BAFF It was reported that the dissolved BAFF trimer could aggregate to form a 60-mer, and this 60-mer was thought to be more active than the BAFF trimer. In addition, L83 of BAFF specifically contributed to the formation of the BAFF 60-mer [77]. As a result, whether the antibody could interact with L83 of BAFF became the key point of consideration for the inhibition of BAFF. For the BAFF-belimumab complex, the heavy chain of belimumab could bind with the L83 of BAFF via hydrophobic interactions. This hydrophobic effect brought the light chain close to BAFF and caused hindrance that prevented the DE loop of BAFF to interact with the adjacent BAFF homo-trimer. This hindrance consequently blocked the formation of BAFF 60-mer, and inhibited soluble BAFF more specifically and effectively, compared with its natural receptors. But for the BAFF-tabalumab complex, the light chain of tabalumab could not combine with L83 of BAFF, so that no steric hindrance existed between the light chain of tabalumab and the adjacent BAFF trimer. As a result, tabalumab could not prevent the formation of BAFF 60-mer, and this seemed to be the very reason why tabalumab failed in clinical trials.
By alignment of the sequences of belimumab and tabalumab (Fig. 6), it was noted that R30 was in light CDR1 of tabalumab, instead of Y30 in belimumab. It was interesting to note that the difference in single amino acid of CDR1 caused the two antibodies to function differently. The amino acids corresponding to R27 in light CDR1 and G94 in light CDR3 could not be found in tabalumab, which made the light chain of tabalumab unable to form hydrogen bond with BAFF. Both the heavy CDR3 of belimumab and tabalumab had DxL motif. They could bind to the active site of BAFF and competitively inhibit the effect of BAFF, which showed that the heavy chain of belimumab and tabalumab shared a similar mechanism in the process of inhibiting BAFF activity, except prevention of BAFF 60-mer formation.

Conclusions
This investigation successfully constructed, in silico, the theoretical complexes of BAFF with belimumab and tabalumab respectively, through homologous construction. The binding model between BAFF and belimumab was analyzed through MD simulation. The amino acids that played a key role in the interaction were determined, and the corresponding binding energy values were calculated through the energy decomposition. Contribution of I51, F54, K58 D100, D101, L102, L103, and P105 of belimumab heavy chain and R27, Y30, K49, S65 of belimumab light chain, to the BAFF-belimumab interaction, was primarily by hydrogen bonds, salt bridges, and hydrophobic effects. Both belimumab and tabalumab had a DxL motif, which was the main binding site of the natural receptors with BAFF, could cause competitive inhibition with the natural receptors. The light chain of belimumab could contact with L83 of BAFF which was the key residue for the formation of BAFF 60-mer, and then the steric hindrance between belimumab light chain and adjacent BAFF trimers potentiated its inability to form BAFF 60-mer. Tabalumab contact with BAFF at the similar site as belimumab and the natural receptors while its light Fig. 6 The alignment result for both light chain and heavy chain of belimumab and tabalumab for a the alignment result of light chain, b the alignment result of heavy chain; the section included in the yellow square was the CDR1/2/3 of belimumab and tabalumab according to the label, the section included in the yellow square was the DxL motif chain could not contact with L83 of BAFF could not affect the formation of BAFF 60-mer and had a less inhibitory effect. Therefore, it is reasonable that belimumab showed better clinical therapeutic efficacies than tabalumab and was approved by FDA for SLE in 2011. In summary, our computation studies may provide the foundation for designing novel therapeutic antibodies in future.
Author contribution Sun J and Wei J conceived and designed the experiments. Jiang Y performed the experiment. Wei J and Jiang Y analyzed the binding modes. Jiang Y and Sun J prepared the manuscript. All authors have read and approved the final manuscript to be submitted for publication.
Funding This study was supported by the National Natural Science Foundation of China (Grant No.81273308).

Data availability Not applicable.
Code availability Not applicable.

Conflict of interest
The authors declare no competing interests.