Diabetes mellitus caused by mutations in human insulin: analysis of impaired receptor binding of insulins Wakayama, Los Angeles and Chicago using pharmacoinformatics

Several naturally occuring mutations in the human insulin gene are associated with diabetes mellitus. The three known mutant molecules, Wakayama, Los Angeles and Chicago were evaluated using molecular docking and molecular dynamics (MD) to analyse mechanisms of deprived binding affinity for insulin receptor (IR). Insulin Wakayama, is a variant in which valine at position A3 is substituted by leucine, while in insulin Los Angeles and Chicago, phenylalanine at positions B24 and B25 is replaced by serine and leucine, respectively. These mutations show radical changes in binding affinity for IR. The ZDOCK server was used for molecular docking, while AMBER 14 was used for the MD study. The published crystal structure of IR bound to natural insulin was also used for MD. The binding interactions and MD trajectories clearly explained the critical factors for deprived binding to the IR. The surface area around position A3 was increased when valine was substituted by leucine, while at positions B24 and B25 aromatic amino acid phenylalanine replaced by non-aromatic serine and leucine might be responsible for fewer binding interactions at the binding site of IR that leads to instability of the complex. In the MD simulation, the normal mode analysis, rmsd trajectories and prediction of fluctuation indicated instability of complexes with mutant insulin in order of insulin native insulin < insulin Chicago < insulin Los Angeles < insulin Wakayama molecules which corresponds to the biological evidence of the differing affinities of the mutant insulins for the IR.


Introduction
Diabetes mellitus (DM) is a heterogeneous group of metabolic diseases characterised by chronic hyperglycaemia resulting from decreased insulin secretion and/or increased insulin resistance. Insulin is a globular protein molecule critical for hormonal regulation of mammalian metabolism and plays a key role in the regulation of blood glucose levels and fat storage. Insulin consists of two chains, A and B with 21 and 30 amino residues, respectively, constrained by one intra-and two inter-chain disulphide bonds (Adams et al., 1969). The A chain (InsA) contains two α-helices ranging from amino acid 1 to 8, and 12 to 18, while in case of B chain (InsB) the helix range is amino acids 9 to 19. The cellular activities of insulin are initiated by insulin binding to the cell surface insulin receptor (IR), which is a disulphide-linked (αβ) 2 heterodimer. The extracellular portion of each αβ protomer contains six domains (L1, CR, L2, FnIII-1, FnIII-2 and FnIII-3) and an insert domain (ID) within FnIII-2. The α-chain component of the ID is terminated by a segment termed αCT, spanning residues 704-719. Defective binding of structurally abnormal insulin or its precursor, pro-insulin can result in diabetes mellitus. Mutations in naturally occurring insulin molecules affect insulin activity through alterations in binding to the IR. There are several naturally occurring mutations in the insulin molecule which have been described in humans and are associated with diabetes mellitus. These include insulin Wakayama (Val 3A → Leu 3A ) where valine is mutated to leucine at A3, insulin Los Angeles (Phe 24B → Ser 24B ) and insulin Chicago (Phe 25B → Leu 25B ) in which phenylalanine is substituted by serine and leucine at B24 and B25, respectively. These substitutions occur within the crevice adjoining the IR binding surface and impair receptor binding by many fold in affected individuals (Wan et al., 2005;Wollmer et al., 1981). The orientations of amino acids Val 3A , Leu 3A , Phe 24B , Ser 24B , Phe 25B and Leu 25B in natural and mutant insulin molecules are depicted in Figure 1.
These mutations on the natural insulin molecule especially at A3, B24 and B25 reduce the binding affinity for IR manyfold (Hayashi et al., 2011). Moreover, there are substantial differences in the biological activity between individual mutant molecules. There has been no published pharmacoinformatics analysis of this reduction in binding affinity and the aim of this study was to correlate the biological activity with the structural analysis.
It has been reported that Val A3 , present within the N-terminal α-helix of A-chain plays an important role in the insulin fold and is recognised as an essential element (Hua et al., 1996;Nakagawa & Tager, 1992;Weiss et al., 2000). The side chains of Val A3 pack within the non-polar crevice between InsA and InsB which adjoin the classical receptor binding surface (Baker et al., 1988;Blundell et al., 1971;Nakagawa & Tager, 1992;Pullen et al., 1976). The crevice is surrounded by sequence-continuous side chains of Ile A2 and Leu B11 , the internal disulphide bridge between Cys A6 and Cys A11 , Tyr A19 and Leu B11 in the hydrophobic core and two amino residues of InsB such as Tyr B26 and Pro B28 which project from an overlying B-chain. It was expected that substitution of Val A3 by Leu A3 would not change conservation and stability because leucine demonstrates a greater intrinsic α-helical propensity than valine (Lyu, Liff, Marky, & Kallenbach, 1990;O'Neil & DeGrado, 1990;Padmanabhan, Marqusee, Ridgeway, Laue, & Baldwin, 1990). However, it was observed that such substitutions reduced activity by 500-fold relative to that of normal insulin molecule (Nakagawa & Tager, 1992). Sequences of natural insulin are remarkable for an invariant phenylalanine present at positions B24 and B25 and for the conservation of phenylalanine and tyrosine at positions B25 and B26, respectively (Conlon, 2001). Phe B24 , Phe B25 and Tyr B26 amino residues play crucial role for the stability, self-assembly and activity of the insulin molecule (Mirmira, Nakagawa, & Tager, 1991;Mirmira & Tager, 1989;Zakova et al., 2013). In 1969, it was shown that (Adams et al., 1969;Baker et al., 1988) on the basis of the crystal structure of the zinc insulin hexamer that residues at the positions B24 to B26 participate in an aromatic-rich dimer interface and this interface occurs three times in the hexamers. Owing to substitution of aromatic phenylalanine at positions B24 and B25 by non-aromatic serine and leucine respectively in insulin Los Angeles and Chicago, binding affinity towards IR is reduced significantly (Weitzel, Bauer, & Eisele, 1978;Wollmer et al., 1981).
Apart from crystallography, molecular docking approaches which take two or more structures as input and predict the structure of their complex are useful for analysing protein-protein interactions. The 3D complex of protein-protein interaction can be used to estimate the effect of mutations (Benedix, Becker, de Groot, Caflisch, & Bockmann, 2009;Kamisetty, Ramanathan, Bailey-Kellogg, & Langmead, 2011;Moretti et al., 2013), and thus for protein design (Ahmad et al., 2012;Baker, 2006;Chen & Keating, 2012;Kortemme & Baker, 2004) and determining the functional consequences of mutations associated with disease (Nisius, Sha, & Gohlke, 2012). Aspects of the conformational variability of a protein molecule can be modelled using the molecular dynamics (MD) simulation method. In MD simulation, the detailed time evolution of a system of particles may be evaluated by integrating equation of motion for the particles on their potential energy surface. The MD method may provide effective modelling of protein-protein interactions at a reasonable computing cost. In order to understand the mechanism of deprived binding affinity of the insulin molecules for the IR, the present study evaluated proteinprotein molecular docking by considering mutant insulin molecules as ligands and recently crystallised IR as receptor. Furthermore, MD simulations were performed to analyse the stability of complexes and exploration to provide insight into the binding interactions between mutant insulins and binding residues of IR.

Materials and methods
In order to understand impaired binding interactions of mutant insulin molecules to the IR, protein-protein docking and MD study were performed using the online ZDOCK server (Pierce et al., 2014) and Amber 14 (http://ambermd.org/#Amber14), respectively. The crystal complex of the IR with insulin molecule (Protein Data Bank [PDB] ID: 3W14  and natural insulin (PDB: 3IR0) were obtained from RCSB PDB (Parasuraman, 2012). The DeepView (Guex & Peitsch, 1997) program was used to mutate the residues on natural insulin and the spatial orientation was maintained. The online server contact map analysis (Sobolev et al., 2005) was used to analyse the binding interactions between insulin and IR.

Protein-protein docking
In the current research, ZDOCK server (Pierce et al., 2014) was used for the protein-protein docking study. The mutant insulin molecules and IR were considered as ligand and receptor molecule, respectively. The ZDOCK is the rigid-body protein-protein docking program which uses the Fast Fourier Transform algorithm (Chu & Alan, 2000) that enables an efficient global docking search on a 3D grid. This program ranks top hundred probable predictions by a combination of shape complementarily, electrostatic and statistical potential terms for scoring. It is reported that ZDOCK achieves high predictive accuracy with more than 70% success in the top thousand predictions (Pierce, Hourai, & Weng, 2011). Jobs of the server are run on a dedicated computer server at the University of Massachusetts Medical School (Pierce et al., 2014) which provides interfaces of input two protein structures and email.

Molecular dynamics
In order to check stability and binding orientation MD simulations were performed for crystal complex of insulin bound to the IR, and best docked complexes between mutant insulins and IR. The MD was carried out using GPU-based PMEMD program and the protein structures were parameterised and represented with amber force field (leaprc.ff99SB) as implemented in Amber 14 (Gotz et al., 2012). The generalised Born (GB) solvation model was used for the implicit solvent with the igb set to 5 which has been shown to agree better with the Poisson-Boltzmann treatment in calculating the electrostatic part of the solvation free energy (Feig et al., 2004). In the production MD of the GB simulations, the Langevin thermostat was used with a collision frequency of 2.0 ps-1 and used an improved GB-Neck of the GB-OBC model that correspond to mbondi2 radii (Onufriev, Bashford, & Case, 2004) as our starting radii set. It has been shown that GB-OBC worked quite well with mbondi2 radii (Nguyen, Roe, & Simmerling, 2013). It is also reported that slight increases in structural stability of the complexes can be achieved in the performance of GB model (HCT or OBC GB models) by setting the atomic intrinsic radii used to define the dielectric boundary of hydrogens bound to nitrogen, H(N) from their Bondi radii of 1.2-1.3 Å (forming the mbondi2 radius set) (Mongan, Simmerling, McCammon, Case, & Onufriev, 2007). Initial minimisation of all the systems was carried out for 1000 steps of steepest descent followed by 1000 steps of heating up of the system to 300.00 K and equilibration. The MD production was carried out for at least 10 ns for each of the system. The frames were recorded at every 200 steps of simulation. The analysis of the conformational changes in the trajectory of the protein was analysed using bio3d (Grant, Rodrigues, ElSawy, McCammon, & Caves, 2006) and the secondary structures were predicted using DSSP (Kabsch & Sander, 1983). The visualisation of the results was done using VMD (Humphrey, Dalke, & Schulten, 1996), chimera (Pettersen et al., 2004) and Pymol Molecular Graphics System, Version 1.7.2 Schrödinger, LLC.

Crystal complex of IR and native insulin
The crystal structure (PDB ID: 3W14) of IR bound with native insulin is depicted in Figure 2(a), which explains that helix of InsB (B7-B19) engages with the C-terminal end of L1 domain but InsA has no interaction with domain L1. The domain αCT interacts with both InsA and InsB of the insulin molecule. It was found that two important amino residues, His710 and Phe714 present in the αCT domain were found to be critical for interaction with both chains of the insulin molecule. His710 inserts into a pocket surrounded by insulin residues Val A3 , Gly B8 , Ser B9 and Val B12 , while Phe714 occupies a hydrophobic crevice formed by amino residues Gly A1 , Ile A2 , Tyr A19 , Leu B11 , Val B12 and Leu B15 (Ward, Menting, & Lawrence, 2013). Another amino acid, Asn711 was directed towards Gly A1 , Val A3 and Glu A4 . Two amino residues of InsB such as Val B12 and Tyr B16 were found to be critical to interact with L1 domain. Val B12 was connected to Phe39, Phe64 and Arg65, while Tyr B16 clashed with Phe39 (Ward et al., 2013). It was observed that αCT domain is an important domain for interactions with both chains of native insulin molecule.

Molecular docking
The molecular docking study was performed with the mutant insulin molecules, Wakayama, Los Angeles and Chicago in the receptor cavity of IR. The native insulin molecule bound in the IR in crystal structure (PDB ID: 3W14) was deleted to create the space for the variant insulin molecules. In order to remove any structural strain or conformational changes the mutated insulin molecules were considered for 1-ns MD simulation analysis. The rmsd vs time along with number of frame vs. total potential and kinetic energies are given in Supplementary Figure S1. The simulated conformers were aligned with respective insulin molecule to calculate rmsd values. The rmsd values for insulin Wakayama, insulin Los Angeles and insulin Chicago were found to be 1.001, 1.043 and 1.005 Å, respectively. The aligned structures are portrayed in Supplementary Figure S2. The low rmsd values clearly explained that docking using simulated and non-simulated mutant insulin molecules at the IR will give comparable binding interactions. Further both simulated and non-simulated mutant insulin molecules were considered for molecular docking study and were found to have similar binding interactions. The best docked pose of each docking study is given in Figure 2(b)-(d).

Insulin Wakayama
The best docked pose insulin Wakayama molecule (Figure 2(b)) explains that Tyr A19 and Asn A21 are able to create connections with L1 domain of IR. Tyr A19 is connected to Arg65, while Asn A21 is surrounded by Leu37, Phe39, Phe64, Phe96 and Glu97. Both InsA and InsB of the insulin Wakayama interact with αCT domain of IR. Asn A21 , Glu B21 , Arg B22 and Phe B25 were found to be important for interactions with Glu706, Val713, Phe714 and Val715. InsB of the variant insulin molecule also bound to L1 domain through Phe39, Lys40, Arg65 and Glu97. Interestingly, the Leu A3 in the mutant insulin molecule was not found to create any connection with IR. With the exception of the length of the main chain of valine and leucine the structures of both amino acids are almost similar. However, the orientation and position of the methyl group present in the both amino acids are different and this might be a possible explanation for no interaction with IR via Leu A3 in the variant insulin molecule.

Insulin Los Angeles
The binding interactions between insulin Los Angeles and IR are depicted in the Figure 2(c). The complex suggests that both InsA and InsB of variant insulin Los Angeles molecule interact with domains L1 and αCT. The best docked pose explains that Tyr A19 , Cys A20 , Asn A21 , Ser B9 , Val B12 , Leu B15 , Glu B21 , Ser B24 and Phe B25 are crucial amino residues for connecting to domain L1 of IR through Phe39, Lys40, Phe64, Arg65 and Phe96. Two important amino residues Phe714 and Val715 belonging to αCT domain successfully interact with Val A3 and Glu A4 , respectively. Glu706, Asp707, His710 and Val713 along with Phe714 and Val715 of the domain αCT were found to be critical to establish binding interactions with InsB of the insulin Los Angeles through Cys B7 , Ser B9 , His B10 , Glu B13 , Tyr B26 and Thr B30 . It is important to note that Ser B24 in the variant insulin Los Angeles only connects with L1 domain but not with αCT domain. Owing to reduction of surface volume at positions 24 of InsB (Serine is smaller in size compared to phenylalanine) as well as lack of aromatic ring in Ser B24 , there was a lack of interaction with IR αCT domain.
Insulin Chicago Insulin Chicago was considered as ligand and IR as receptor for molecular docking study. The best docked complex is shown in the Figure 2(d). The complex shows that Asn16 and Lys40 are important amino residues of L1 domain of IR for binding interactions with Tyr A19 , Cys A20 and Asn A21 , while Val A3 , Glu A4 , Cys A7 and Thr A8 are connected to αCT through His710, Asn711, Val713, Phe714 and Val715. InsB of the variant insulin Chicago was unable to establish any binding interaction with L1 domain. However, two amino acids belonging to InsB, Glu B21 and Arg B22 clashed with αCT domain via Val713, Phe714 and Val715. It has been reported that Phe B24 , Phe B25 and Tyr B26 are important amino residues for interaction with IR (Mirmira et al., 1991). However, the docking results of insulin Chicago shows that none of the three amino acids are involved in binding interactions. The possible reason may be lack of aromatic ring in Leu B25 and variation on orientation of the functional groups in the insulin Chicago compared to native insulin.

Molecular dynamics
MD simulations are likely to play an increasingly important role in the development of novel pharmacological therapeutics. Thus, it can be expected that MD simulation would provide more reliable structural information on the mutations in insulin. In order to compare the dynamic stability of the best docked poses of insulin molecules inside the IR, the MD simulations were carried out based on protocols described in the materials and methods section to predict the impact of mutation on insulin molecules. MD gives information on the stability of binding interactions between the ligand and receptor cavity. The published crystal structure of IR, bound with native insulin (PDB ID: 3W14) and best docked complexes of IR with mutant insulin molecules were taken into consideration for MD study.
The production of the MD of the complexes were at least for 10 ns for each of the complexes as shown in the energy plots of Figure 3. The energy plot is shown from 7000 ps to the last ps in order to show clearly the fluctuation in the energy at the points of energy equilibrium. The geometry stability was obtained from 12,000 steps of the simulation as shown by the rmsd plot of Figure 3. The rmsd of the complexes from their respective initial geometries during the trajectory is in the order of complexes of insulin Chicago > insulin Wakayama > insulin Los Angeles > native insulin.
Further trajectory analysis of the MD was done using the frames from 10,000 steps (4 ns) to the last step of each of the complexes. The sample of the geometries of the insulin complexes of native insulin, insulin Chicago, insulin Los Angeles and insulin Wakayama with the IR at the first step of MD and at the 10 ns of the MD as shown in Figure 4. The correlation of the common IR part of the complexes with that of the native insulin before the MD was observed to be .098, .070 and .062 Å while at the 10 ns of the MD the correlation changes to 25.241, 21.964 and 20.087 Å, respectively. Even though the change in the rmsd of the IR of the insulin Wakayama is lower than that of insulin Chicago and insulin Los Angeles, its insulin complex was found to be unstable during the MD simulation. At 10 ns of the simulation, the insulin chains of Wakayama have completely dissociated from the IR (measured to a distance of 95.036 Å from closest atom of other aligned insulin complexes as shown in Figure 4) but those of native insulin, insulin Chicago and insulin Los Angeles still maintain a stable interaction with the IR.
The geometrical analysis of the complexes before the trajectory also shows some important distinguishing features in the complexes. The normal mode analysis (NMA) of the complexes shows that insulin Los Angeles has similar directions of vector fields with that of complex of native insulin while insulin Chicago and insulin Wakayama are in opposite directions. Using the NMA simulation techniques to probe large-scale motions in complexes, the highest predicted functional motions in complexes are shown to take place in the mutated region of insulin Wakayama ( Figure 5). Generally in all the complexes, the larger fluctuations were predicted to occur at the loop regions in the fluctuation plot of Figure 5). The prediction of fluctuation shows a gradual increase in the fluctuation of the insulin part of the complexes in the order of native insulin < insulin Chicago < insulin Los Angeles < insulin Wakayama. Based on the NMA prediction for insulin Wakayama, it would expected in the course of dynamics to have high fluctuation in the insulin region which we eventually observed to be true.

Conformational changes in the complexes
The rmsd of the sampled conformation during the MD of complexes with native insulin, insulin Chicago, insulin Los Angeles and insulin Wakayama was found to be predominantly around 6, 7, 13 and 9, respectively, from the starting geometries as shown in the histogram of distribution plots in Figure 6. In order to examine the relationship between different conformations sampled during the trajectory, the principal component analysis (PCA) of the trajectory was carried out. The conformational differences are described by the axes of maximal variance of the distribution of structures using principal components (orthogonal eigenvectors). The three dimensions (i.e. first three PC) sufficiently capture 73.1, 78.2, 75.8 and 77.0% of the total variance in MD trajectory of the complexes of native insulin, insulin Chicago, insulin Los Angeles and insulin Wakayama, respectively, as shown in Figure 7. There are periodic jumps between the different Figure 4. The trajectory of the complexes at step one before the md and at time step of 10 ns during the md. Note: The complexes as distinguished with colour: NAT (blue), CHI (red), LOS (green) and WAK (yellow).TBSD conformers of each complex throughout the trajectory has clearly shown by the continuous colour scale of the PC from blue to white to red as shown in Figure 7. In order to focus on the relationships between individual structures in terms of their major structural displacements, these distinct conformers were further highlighted by clustering structures in PC space as shown in Figure 8.
From the results of the clustering, there are two major groups of conformation for each of the complexes. Figure 8 shows the distinct groups of conformations along the PC1 plane for complexes with native insulin (the larger grouping centred around −50 and a second grouping at +150), insulin Chicago (the larger grouping centred around −50 and a second grouping at +230), insulin Los Angeles (a very small grouping at −330 and a very large grouping centred around +100) and insulin Wakayama (one grouping centred around −200 and a longer grouping at +40). The contribution of the residual fluctuation to the conformational changes along the PC1 and PC2 is shown for each complex in Figure 9. The highest contribution of the insulin residues to the observed conformational variances observed along both PC1 and PC2 is found in the complex of insulin Wakayama and the next significant contribution is observed in the complex of native insulin. The figure also shows the groups of the IR residues that significantly contribute to the conformational changes in the complexes of native insulin , insulin Chicago (500-600, 1100-1200), LOS (200-300, 1100LOS (200-300, -1200 and insulin Wakayama (300-400, 750-850).

Discussion
Naturally occurring mutations of the human insulin molecule are associated with diabetes mellitus and this has given insight into structure-function relationships. The present research was focused on three important mutant insulin molecules, Wakayama, Los Angeles and Chicago which have very low binding affinity for IR compared to the native insulin molecule. In the case of insulin Wakayama, the molecular docking results suggested that the amino residue at position A3 was unable to connect L1, while valine at the same position in native insulin molecule successfully formed interactions with αCT domain of IR (Figure 2). MD study of insulin Wakayama and docked complex between insulin Wakayama and IR indicated instability of the molecule. RMSD and RMSF of both insulin Wakayama and docked complex explained that due to mutation at A3 the stability was disturbed. The branched-chain aliphatic residues changed within the non-polar crevice as valine was replaced by leucine. Such alteration is easy to tolerate in the protein molecule and mutant insulin molecule is almost similar in structure to the native insulin molecule. Observation of several groups (Derewenda et al., 1991;Hua, Shoelson, Kochoyan, & Weiss, 1991;Ludvigsen, Olsen, & Kaarsholm, 1998;Wan, Xu, Chu, Katsoyannis, & Weiss, 2003;Xu et al., 2002) and computational output of the current work explained that crevice related to A3 opens due to binding with IR which enables docking of Val A3 within non-polar pocket of receptor. Due to steric hindrance between Leu A3 and IR, the stability of the complex is reduced and this leads to impaired binding affinity of insulin Wakayama towards IR. In the MD study of insulin Los Angeles and Chicago, it was observed that the average RMSD and RMSF of docked complexes mutant insulin molecules are higher than crystal complex of native insulin. Moreover, the molecular docking results suggested that InsB of insulin Chicago failed to connect with L1 domain of the IR. Amino residue at B24 in insulin Los Angeles was unable to connect to the domain αCT. It is reported that the aromatic side chain is favourable for interactions with IR. In both mutant insulin Los Angeles and Chicago, aromatic amino residue phenylalanine substituted by non-aromatic serine and leucine, respectively, leads to changes in binding interaction with IR. The charged or polar side chains reduce the hormone binding (Pandyarajan et al., 2014) and in insulin Los Angeles, the number of polar side chains is increased and this may reduce the binding affinity towards IR.

Conclusion
In this study, we examined the binding of mutant insulins Wakayama, Los Angeles and Chicago using pharmacoinformatics analysis of the reduced binding affinity to the IR. The molecular docking study of all three mutant insulin was carried out to analyse the binding interaction with the binding amino residues of the IR. The study explained that due to small size and less surface area of the replaced amino acids, the binding interactions were reduced and this altered the binding affinity towards IR. In order to confirm the possible reason for this, MD study was also performed. The rmsd of the protein backbone of the complexes during the trajectory were in the order of insulin Chicago > insulin Wakayama > insulin Los Angeles > native insulin. Aligning the trajectory of complexes with insulin Chicago, insulin Los Angeles and insulin Wakayama with the native insulin at 10 ns of simulation, insulin Chicago has the highest rmsd followed by insulin Los Angeles while insulin Wakayama is least. It was also observed that, the insulin part of native insulin, Chicago and Los Angeles still maintain the stable interaction with IR but a complete dissociation of the insulin from IR was observed for insulin Wakayama during the MD simulation. The PC analysis also shows that the highest contribution of the insulin residues to the two distinct conformational variances observed along both PC1 and PC2 is found in the complex of insulin Wakayama and the next significant contribution is observed in the complex with native insulin. Therefore, the trajectories of the MD successfully explained that possible reasons may be the smaller size and variation in the size chains of the amino acids which affects the binding interactions as well as the stability of the complex. The outcome of the current research can assist in identifying new possible insulin analogues for the effective treatment of diabetes mellitus.

Supplementary material
The supplementary material for this paper is available online at http://dx.doi.org/10.1080/07391102.2016.1160258.

Disclosure statement
No potential conflict of interest was reported by the authors.