Molecular dynamic simulation of Trastuzumab F(ab’)2 structure in corporation with HER2 as a theranostic agent of breast cancer

Trastuzumab as intact IgG are well researched for theranostic agent in HER2 overexpressed breast cancer. However, due to the relatively large of molecules it is slowly moved and weak penetration of the target cells. Fragmentation of trastzumab has been developed by pepsin cleavages to get the F(ab’)2 fragments. To observe the stability and accessibility of F(ab’)2 structure in corporation with HER2 (human epidermal growth factor receptor-2), the structure of antibody modeling had been developed with 1IGT as a template. Molecular dynamics (MD) of the F(ab’)2 structure simulation has been done in the aqueous phase with AMBER trajectories for 20 ns. Computational visualization by VMD (Visual Molecular Dynamics) were applied to identify binding site interaction details between trastuzumab F(ab’)2 and HER2 receptor. The results of MD simulations indicated that the fragmentation of trastuzumab F(ab’)2 did not change the structure and conformation of F(ab’)2 as a whole, especially in the CDR (Complementarity Determining Region) area. SASA (solvent accessibility surface area) analysis on lysine residues showed that formation of conjugate DOTA-F(ab’)2 predicted occur on outside of the CDR regions so its not interfered with binding affinity for the HER2 receptor. The molecular dynamic simulation of DOTA-F(ab’)2 with HER2 receptor in aqueous system generated ΔGbinding more highly (15.5066 kkal/mol) than positive control HER2-Fab (-45.1446 kkal/mol).


Introduction
The use of monoclonal antibody (mAb) in the diagnosis and therapy of cancer has been developed in two decades. Selectivity and specificity of antibodies which high enough against the target cells is one of reason why this techniques has been chosen [1][2] [3]. The mechanism of immunotherapy aimed directly at cancer cells, where there is an antigen or receptor site on cancer cells that are recognized by antibodies. The monoclonal antibodies capable of recognizing specific antigens on the surface of cancer cells. Antibody-antigen interaction typically occurs in CDR areas (Complementary Determining Region). One of the antibodies that has been commonly used for breast cancer therapy is Trastuzumab Mab (Herceptin) [4][5][6].
Trastuzumab, also known in commercial name Herceptin, currently used in a part of the immunotherapy of advanced breast cancers, expecially in extensive metastasis and in general, solid tumors over-expressing HER2 [7]. Different clinical trials confirm the efficiency of the antibody as an anticancer treatment. Trastuzumab binds to HER2-ECD domain. Interaction trastuzumab Fab to domain IV of the HER2-ECD categorizes this site as a possible target for anticancer therapies [5] [8].
Application of trastuzumab in breast cancer therapy still has many weaknesses, due to its size is relatively large (approximately 150 kD), thus resulting in the molecule is less effective in reaching target cells (low pharmacokinetics) and weak penetration of target cells (poor tumor penetration). This phenomena potentially generate an antibody response to the resistance of host cells [9]. To overcome these limitations, an approach that can be done is to reduce the size of the antibody to become smaller fragments by enzymatic hydrolysis or genetic engineering. Enzymatic fragmentation antibodies much easier to do because it does not require complicated steps and costs are relatively large and easy to do. Fragment bivalent F(ab') 2 can be produced from the breakdown of a disulfide bridge at the carbonyl side of the cysteine residues by pepsin hydrolisis which will produce fFc fragments and F(ab') 2 [10][11].
In this research, enzymatic fragmentation of trastuzumab has been conducted using pepsin to hydrolyze part of the antibody fragment that was not a CDR region (Complementary Determining regions). From the in silico results of these fragmentation obtained F(ab') 2 fragments with affinity expected not to change compared with the intact antibody. To study the mechanism of DOTA-F(ab')2-trastuzumab interaction with HER2 receptor conducted studies by in silico approach with molecular modeling and molecular dynamic (MD) of F(ab') 2 . The stability and accessibility of F(ab') 2 trastuzumab explored by molecular modeling of F(ab') 2 with 1IGT murine antibody as a template [12]. The structure modeling of F(ab') 2 results observed through molecular dynamics simulation in the aqueous phase with MD trajectories 20 ns. The Fab simulation used as a reference for the first phase conjugate formation of DOTA-F(ab') 2 which interacted to HER2 receptor and its affinity compared with native structure of complex-HER2 Fab (PDB file 1N8Z) [13]. The affinity interactions explored by MMGBSA between HER2 receptor and ligand F(ab') 2 respectively. In silico study of this immunoconjugates interaction expected to support the further development of a radiopharmaceutical over expressing HER2 positive breast cancer.

Structure modeling and CDR analysis the fragment F(ab') 2
Protein structure modeling and CDR analysis of F(ab') 2 fragment has been done by using antibody modeler in Biovia Discovery Studio 4.5 and Rosetta-based modeler antibody based on sequence data of trastuzumab IgG obtained from http://www.drugbank.ca/drugs / DB00072 [14]. Antibody modeling constructed with 1IGT as a template [15]. Amino acid sequences of F(ab') 2 heavy chain and light chain respectively combined and CDR analyzed through Rosetta software modeler antibody (antibody http://rosie.rosettacommons.org/) [16]. Type of intermolecular interaction visualized by Biovia Discovery Studio 4.5 and the value of RMSD from F(ab') 2 modeling determined by sequence alignment of the template structure. The validity of the F(ab') 2 structure were analyzed using Ramachandran plot Assessment (http: //www.RAMPAGE Ramachandran Plot Assessment.htm) and structure modeling results are stored in pdb format [12]. with several stages of energy minimization. Before the F(ab') 2 structure minimized prior modification of the structure by eliminating protons/H atom and change some amino acid residues of certain significant effect on the conformation of 3D protein and geometry optimization of molecules such as histidine (HID / HIE / HIP) and cysteine (CYS / CYX). Modified structured of F(ab')2 is saved in a new file name (fab2_ready.pdb). Furthermore pdb file run on Leap program (using AMBER 14) with force fields ff14SB mode with the grounds parameterization topology and coordinate files appropriate for protein structure.
Modified structure checked and neutralized by the addition of counter ions (10 Na + units) and solvated with water molecules using solvatebox fab TIP3PBOX10. F(ab') 2 solvation structure stored in a new file (fab2_ready_sol.pdb). Before the salvation parameters are stored, the parameters included in the molecular ion water by: loadamberparams frcmod.ionsjc_tip3p and parameterization file is stored in the form: saveamberparm c fab2_ready_sol.prmtop fab2_ready_sol.inpcrd [17].

Minimization.
Minimization in the SANDER module done by using two energy minimization functions, steepes descent and conjugate gradient. Stepees descent mode in the early stages minimization in which all atoms contained in the system with a constant force field restrain 100 kcal /mol except for water molecules and counter ions, carried 1,000 minimization stage. The second minimization is done by conjugate gradient mode where's all done restrain protein molecules except for a hydrogen atom with an energy restrain 5 kcal mol -1 as many as 2000 stage. Minimization of the third stage is done by providing harmonic restrain (protein backbone) with a 2000 phase minimization. The last stage minimization all the molecules is gradually unrestrain by 5000 stages minimization to minimize steric repulsion between atoms [18][19].

Heating.
The system heated gradually until it reaches physiological body temperature (310 o K) with phases (0-100K: 100-200K: 200-310K, respectively 20ps) in a state NVT using harmonic restraints on the main frame of the atoms (backbone).

Equilibrium.
After the heating process, the system reaches equilibrium at a temperature of physiological (310 o K) with weak restrain the protein, further equilibration in conditions of NTP by removing restrain gradually and move at a constant pressure with phase equilibration each for 100 ps, until the final stage where the molecule not restrained.

MD production.
Production trajectory run at NTP conditions. Temperature was controlled by Langevin thermostat with collision frequency 1 ps-1 and pressure controlled with Beredsen barostat with coupling constants 1 ps and a target pressure of 1 bar. Non bond cut-off value of 10 A is used and the electrostatic bond length is measured by Particle Mesh Ewald. Time step used during production is 2 fs and the time scale (running trajectory 20 ns), all hydrogen atoms are given constraint by using SHAKE algorithm. Short range repulsion-dispersion interactions were smoothly truncated at 10 Å. The particle mesh Ewald (PME) method was used to calculate long range electrostatic interactions, with a maximum grid spacing of 2.5 Å and using fourth-order (cubic) interpolation for the fast Fourier transforms [20].

Trajectory
Analysis. Trajectory analysis carried out by using ptraj module with the topology file simulation results (fab2ready_dry.prmtop). In order to be analyzed, all the coordinate files from MD simulation (in prod* .crd) combined with a script file (trajin/prod*.crd) into a trajout file (20ns_dry_fab2.binpos binpos). Furthermore, the combined results of the simulation file in the form file  2 . Structure accessibility of F(ab') 2 binding DOTA conjugates determined by calculating the SASA (solvent accessibility of surface area) parameters at all outside lysine residues in the F(ab')2 molecule. The degrees of freedom of lysine residues obtained were plotted based on the SASA parameters by using trajectory file from simulation. The conjugation accessibility visualized through the program VMD (Visual Molecular Dynamics) 1.9.1.

Molecular dynamic of HER2-Fab and HER2-DOTA-F(ab') 2 .
HER2-Fab complex structure (positive control) obtained from the crystal structure (1N8Z) by adding some missing residue using PRODAT databases based on homology to sequences obtained from NCBI.
Modeling and analysis of complex structures (HER2-Fab) conducted in PDBsum (http://www.chemogenomix.com/chemogenomix/ PDBsum.html). The structure modeling of HER2-Fab crosscheck by the stereochemistry mapping with Ramachandran plot to observe the favorable conformation of each amino acid residues in the region are permitted (allowed region). Molecular dynamics simulations of Fab with HER2 receptor is done through the program AMBER 14 in a state of solvated complex water at pH 7 and physiological temperature (37°C). The simulation results of HER2-Fab complex was observed through cpptraj program by plotting RMSD value over time and comparing the fluctualization residues on each structure (ligands and receptors). Inventory of hydrogen bonding performed on the CDR regions involved in the interaction with the receptor. Visualization distance and probability (occupancy) formation of hydrogen bonds between the donor and acceptor analyzed through cpptraj on each residue based on the simulation results MD trajectory file. The calculation of free energy value carried out by MMGBSA module where all frames conformation from MD simulation results (complex, ligand and receptor) were included in the calculation to obtain the free binding energy [21][22].

Structure Modeling of F(ab') 2 .
Before the molecular dynamic simulation , structure modeling of F(ab') 2 fragments carried out by using antibody modeller in Biovia Discovery Studio 4.5. Molecular dynamic (MD) simulation refers to the observation of conformational before and after conjugated with pSCN-DOTA molecules and predict accessibility for binding conjugate based on the SASA (solvent assesibility surface area) analysis. The results of structure modeling of F(ab') 2 as Figure 1.
Based on antibody modeling of F(ab') 2 structure with 1IGT as a template, it was appears that the F(ab') 2 fragment have a large similarity index sequences (78 408%) after superimpose with the 1IGT template by the RMSD value 1.5389 Å (Figure 2)., The overlay results indicate that the structure modeling of the F(ab') 2 results have a similar conformation to the template. To ensure the validity of the structure modeling, the sequence alignment carried out to the template 1IGT in SCR areas (structuraly conserve region) and SVR (structuraly variable region) through multiple sequence alignment (MSA) by the method of Psi-BLAST Biovia software contained in Discovery Studio 4.5. This aligment to ensure that the structure of the F (ab') 2 used in accordance with sequences that are desired. The MSA can also be known that the area is maintained (conserved region) of the amino acid sequence fragments including the CDR (Complementary Determining Region). The result of the sequence alignment F(ab') 2 tratuzumab are as follows: In the sequence alignment result (Figure 3) shows that the amino acid sequence of F(ab') 2 are retained (residues 1-220) naturaly conserve with relatively high homology especially in the CDR area (1-100). Although both include the type of F(ab') 2 trastuzumab and 1IGT (immunoglobulin) relative differ mainly in the area of non-CDR because of 1IGT is a chimeric monoclonal antibody types while F(ab') 2 is a type of humanized antibody. Chimeric monoclonal antibody is a combination of all parts of the Fv antibody mice with human antibody Fc (30% mouse antibody). Meanwhile, the humanized antibody (10% mouse antibody), did not take up the whole of a mouse antibody Fv part, but only part of the CDR only [23].
On the other hand, the mapping of the secondary structure of F(ab') 2 by using Ramachandran Plot Assessment (http: //www.RAMPAGE Ramachandran Plot Assessment.htm) (Figure 4) shows the stereochemical data structure was thermodynamically favorable where more 94% of the amino acid residues located in areas favorable and allowed regions and only 0.6% of residues located outside the region are allowed (disallowed region).  Determination of CDR regions is done by using Rosetta antibody modeler [16], where the analysis results for each CDR sequences of the light chain and heavy chain are as follows ( Figure 5).

Figure 5. CDR regions of Trastuzumab F(ab') 2 Light Chain and Heavy Chain
Based on analysis of CDR, obtained that CDR1 light chain regions of F(ab ')2 consist of 11 amino acid residues (residues 24-34), CDR2 region (residues 50-56) and CDR3 (residues 89-97) consists of 7 and 9 amino acid residues. In the heavy chain, CDR1 region consisting of 10 amino acid residues (residues 26-35) and CDR2 consists of 17 amino acid residues (residues 50-66) and CDR3 consists of 11 amino acid residues (residues 100-110). Each region is a part of the interface that is likely to interact strongly with domain area IV HER2 receptor. As predicted crystal structure of complex HER2-Fab, where residue contact amino acid Fab located at residue Trp99 that will interact with Phe573 through hydrophobic interactions, Asn30, Tyr92, Thr94, Arg50 and Gly103 to Gln602, Lys569, Asp560, Glu558, Asp560 and Lys593 through the interaction of hydrogen bonding and ionic interactions in the Asp560 residues Arg59 and Arg50 with Glu558 and Asp560. However, interaction need to be explored further through molecular dynamic simulation approach by mapping the interactions between the HER2 receptor structure that has been modeled Fab (PDB ID. 1N8Z) [13].
Accessibility and flexibility of the F(ab')2 structure were analyzed by observing changes in the value of RMSD and RMSF on the backbone structure of F(ab') 2 during the simulation period with 20 ns timescale at NTP conditions (temperature and constant pressure). The value of RMSD and RMSF is calculated based on the trajectory files that have been incorporated in the module cpptraj. RMSD (root mean sequare deviation) showed deviation in 3D structure during the simulation process F (ab')2 as a function of the coordinates of the starting position. On another hand RMSF (root mean square Fluctuation) is a measurement of the average mobility of atoms / residues.

Flexibility structure of F(ab') 2
Based on RMSD fluctualization (Figure 6), it appears that the conformation of the 3D structure of F(ab') 2 relatively stable even changes in the value of RMSD relatively large (5-20A). RMSD changes in F(ab') 2 in total residue due to the flexibility of the hinge region were quite free as a consequence of antibody fragment cuts off. Conversely, if compared with the value of RMSD in region 1 and region 2, which contains CDR regions, the structure and conformation of F (ab') 2 is relatively more stable (RMSD <2A). The stability of the F (ab') 2 structure is likely influenced by the formation of intermolecular disulfide bridges, hydrogen bonding and electrostatic interaction between aspartic acid residues, glutamic acid and lysine. Flexibility of F(ab') 2 structure in the hinge region became the reason why the structure intact antibodies are generally difficult to be crystallized [21] [24].
Flexibility and structure dynamics of F (ab') 2 at the residue levels can be explored more detail through the analysis of RMSF where mobility and flexibility of each residue is determined on the region 1 and region 2 area by comparing fluxtuality with the Fab structures (positive control). Fab structure derived from complex structure of HER2-Fab crystals (1N8Z) [13] [21]. Based on the results of RMSF seen that almost all amino acid residues F(ab') 2 has relatively small fluxtuality when compared with the structure of the Fab with RMSF values below 2 except for a few specific residues (Figure 7).

Accessibility of F(ab ')2 and the conjugate mechanism of formation of DOTA-F(ab') 2
Degrees of freedom of lysine residues and its accessibility to the pSCN-Bn-DOTA is determined by calculating the SASA on each lysine residue. From these calculations can predict where the all the most potential lysine residues exposed to the solvent so that thermodynamically allows for the conjugation reaction. SASA analysis of lysine residue results shown in Figure 8. Based on the SASA analysis, there are three lysine residues on each arm of the F(ab') 2 structure which has high exposure to solvents, ie. Lys169 (residual light chain / non-CDR), Lys190 (residual light chain / non-CDR ) and Lys208 (residue of heavy chain / non-CDR) (Figure 8). The structure of the F(ab') 2 consists of two arms (Fab' and Fab ''), so there are 6 lysine residues that have a high enough exposure and are outside the CDR regions [25].

HER2-Fab molecular dynamic simulations in the water phase
Molecular dynamic simulations of HER2-Fab conducted for 20 ns on NTP state (pressure and temperature constant). This simulation was performed with the same timescale as the simulation F(ab') 2 in a free state. The molecular dynamic simulation results show the RMSD value change as follows ( Figure 9):  203  304  405  506  607  708  809  910  1011  1112  1213  1314  1415  1516  1617  1718 1819 1920

RMSD of HER2-Fab dan Fab unbound
Reseptor Ligan Fab Fab unbound Based on the RMSD value (Figure 9), HER2-Fab complex structure is relatively more flexible with RMSD value larger than the Fab structures in a free state / not bound to the receptor. Flexibility HER2-Fab complex structure likely influenced by the adjustment phase during the formation of complex systems and their interaction loop in domain II of HER2 with CH region (heavy chain) of Fab molecules. This was confirmed by a shift in conformation observed during the simulation period, as shown in Figure 10.
Conformational change of HER2-Fab complex in the water phase (during 1-10 ns) is an adjustment of induced fit effect between the Fab molecules (ligands) with HER2 (receptors) where local loop (domain II) HER2 change the orientation of approaching the CH molecule ligand complex Fab resulting in more stable interactions. HER2-Fab complex formation at the interface region is mediated by the presence of non-covalent interactions between residue found on the contact area with several Fab CDR amino acid residues in the interface ( Figure 10). These interactions include intermolecular hydrogen bonding interactions, ionic interactions (salt bridges) and hydrophobic interactions (phi-phi aromatic). The third type of interaction is instrumental in the adjustment phase of the system so as to produce a stable compound [25][26]. Through the calculation of the value of binding energy by using MM-GBSA (Poisson-Boltzmann surface area) values obtained for G binding -45.1446 kcal / mol. If compared with the results of the wet experiment analysis of G binding (12.4-14.10 kcal / mol), the value of the free energy of binding ligands relatively larger than the experimental values, It is likely influenced by the value of the entropy of the system in which the in silico simulation calculations entropy values are not fully included in the calculation [13] [26].
Molecular dynamic simulation of the HER2 receptor with DOTA-F(ab') 2 prepared by the same procedures as the positive control (HER2-Fab). The binding energy result of complex interactions between HER2 and DOTA-F(ab') 2 calcultaed by MM-GBSA module generates the binding free energy (G binding ) amounted to 15.5066 kcal/mol. When compared with receptorligand binding affinity in the stimulation of HER2-Fab (positive control), the value of this G binding more positive, it means the relative binding affinity decreases. This is probably due to the presence of DOTA molecules in ligand structure F(ab') 2 which indirectly affect the stability of the complex conformational HER2-F(ab') 2 thus becomes a lower binding affinity.

Conclusions
The results of the molecular dynamics simulation F (ab ') indicates that fragmentation trastuzumab Fc molecules, the relative does not damage the structure and conformation of F(ab') 2 as a whole, especially in CDR areas that do not degrade binding affinity to HER2 receptors. The results of the SASA analysis on lysine residues of F(ab') 2 show that the reaction of formation of conjugate DOTA-F (ab') 2 is occured outside of the CDR regions so it does interfere with binding affinity for the receptor HER2. The results of the molecular dynamics simulation DOTAF(ab') 2 -HER2 in the aqueous phase to produce a relatively larger G binding (15.5066 kcal / mol) compared with HER2-Fab complex (-45.1446 kcal / mol) or in other words the binding affinity decreases.