This study investigated the potential targets of Dolastatin 16 based on network pharmacology approach. Several steps to obtain the most appropriate target of Dolastatin 16 included target identification using databases, PPI interaction, topological analysis using Cytoscape 3.9.1 software, and target accuracy. From 202 related targets, there were top 10 targets with higher degree resulted from topology analysis, then matched to the target prediction databases to determine the target with the highest accuracy. MMP9 was the target with the highest accuracy based on Super PRED and SWISS Target database, then it was decided as the target focus in the next analysis.
10 core targets were proceeded to disease analysis, gene ontology biological process, and KEGG pathway. Diabetic neuropathy (e.g., Diabetic Polyneuropathies, Diabetic Mononeuropathy, Asymmetric Diabetic Proximal Motor Neuropathy, Diabetic Asymmetric Polyneuropathy, Diabetic Autonomic Neuropathy, Diabetic Amyotrophy, Symmetric Diabetic Proximal Motor Neuropathy, and Diabetic Neuralgia) was the most possible disease based on the result in Fig. 4a. There are top 4 biological process involved MMP9 role including positive regulation of cell migration (GO:0030335), negative regulation of apoptotic process (GO:0043066), apoptotic process (GO:0006915), and positive regulation of protein phosphorylation (GO:0001934). These terms will associate with the process in the Lipid and atherosclerosis (hsa05417), Fluid shear stress and atherosclerosis (hsa05418), and TNF signaling pathway (hsa04668). Commonly, atherosclerosis and diabetes cases are relatively associated with chronic inflammation that involve TNF signaling pathway, but this study only focused on the diabetic condition.
Diabetic neuropathy is the most prevalent diabetes complications affecting sensory function. Over the years, around 50% of patients with diabetes develop diabetic neuropathy which is most affecting patients with type 2 diabetes mellitus (6,100/100,000 person per years) (Feldman et al. 2019). This chronic condition leads to issues with neuropathic pain, painless foot ulcers, and autonomic dysfunction (Bönhof et al. 2018). Among other types of diabetic neuropathies, diabetic polyneuropathy is most common condition that frequently leads to foot ulceration (De Gregorio et al. 2020). In the diabetic population, foot disorders contribute to around 85 percent of all amputation cases (Begum et al. 2022). In the next decades, the cases will exacerbate if no effective medications are applied.
TNF signaling pathway plays a vital role in the pathology of diabetic foot ulceration (DFU). Tumor necrosis factor alpha (TNF-α), a family of TNF, is a potent of proinflammatory cytokine produced by macrophages during acute inflammation. It is also evidenced in the patients with infected diabetic foot ulcer (Dhamodharan et al. 2019). In the signaling pathway, the inflammation path involves NF-κB signaling to promote MMP9 activation during diabetic ulceration. The five members of the NF-κB family include p65 (RelA), RelB, c-Rel, p105/p50, and p100/p52. When the NF-κB pathway is activated, the heterodimers RelA and p50 are the pivotal mediator for transcription of target genes. The illustration is shown in Fig. 9. The pathway mechanism was adopted from https://www.kegg.jp/pathway/hsa04668 and literature, then created in BioRender online tool (https://www.biorender.com).
In the inflammation cases, foot ulceration releases TNF-α extracellular signaling, then transferred to its receptor TNFR-1 (TNF receptor-1) containing DD (Death Domain). TNFR1 recruits TRADD (TNF receptor activated protein with death domain), TRAF2/5 (TNF receptor associated factor-2 and 5), and E3 ubiquitin ligase cIAP1/2 (cellular inhibitor of apoptosis 1 and 2). It also recruits RIP1 (Receptor interacting protein 1) (Hsu et al. 1995, 1996; Baker et al. 2011). Once the entire receptor complex is completed and activated, NF-κB pathway will be regulated.
Ubiquitination plays the important roles in the NF-κB activation. To begin with, the K63-linked ubiquitination is attached to RIP1 and recruited to NEMO to form TAK1-IKK complex. IKK comprises of two homologous catalytic subunits known as IKKα (IKK1), IKKβ (IKK2), and IKKγ (known as NF- κB essential modulator or NEMO). TAK1 phosphorylates and activates IKKβ. In the absence of external stimulus, NF-κB dimers are remain inactive due to sequestered in the cytoplasm by a family protein inhibitor IκBα which prevents the NF-κB dimers from entering the nucleus. To activate the NF-κB signaling, IKKβ mediates the phosphorylation of IκBα which is the crucial step in NF-κB activation. The phosphorylated IκBα associates with βTrCP protein to initiate the K48-linked ubiquitination of IκBα. The ubiquitinated IκBα is directed to the cytosolic proteosome and degraded, then NF-κB dimers (p50 and p65) translocate to the nucleus to bind to the DNA and drive gene expressions such as remodeling of extracellular matrix which involves MMP9, MMP3, and MMP14 (Hayden and Ghosh 2008; Israël 2010; Wertz and Dixit 2010).
In the inflammatory phase, the healing process in diabetic condition is different with the wound healing in normal individuals. In diabetic cases, the wounds are more vulnerable to bacterial infection due to impaired immune response, which results in the increased levels of pro-inflammatory cytokines, proteases, and reactive oxygen species (ROS) and the damage of structural components of extracellular matrix (ECM) (Cho et al. 2019;Demidova-Rice et al. 2012). ROS, combined with pro-inflammatory cytokines, can aggravate the wound by elevating the expression of matrix metalloproteinase (MMPs), leading to the destruction of the ECM and growth factors (Chang 2016).
Matrix metalloproteinases (MMPs), a family of zinc-dependent endopeptidases, are responsible for physiological and pathophysiological tissue remodeling. MMPs contribute to various wound healing phases by destroying almost all ECM protein components (Nagase et al. 2006; Yabluchanskiy et al. 2013). MMPs mainly classified as collagenases (such as MMP1, MMP8, MMP13), matrilysins (such as MMP7), stromelysins (such as MMP3, MMP10, and MMP11), gelatinases (such as MMP2 and MMP9), and membrane type metalloproteinases (such as MMP14) (Laronha and Caldeira 2020;Rohani and Parks 2015). In the diabetic foot ulcer condition, the higher level of MMPs containing protease activity causes protein degradation, ECM deterioration, granulation tissue remodeling, angiogenesis modulation, and growth factor activity regulation (Chang and Nguyen 2021;Park et al. 2014;Tardáguila-García et al. 2019).
MMP9 consists of the catalytic domain that contains two zinc ions, five calcium ions, and three repeats homologous to the type II module of fibronectin. MMP9 plays a role in the degradation of extracellular matrix (ECM) for both physiology and pathophysiology processes involving tissue remodeling (Yabluchanskiy et al. 2013). It is also involved in inflammatory process such as diabetes, arthritis, and cancer (Halade et al. 2013). In the diabetic foot condition, the higher level of MMP9 is the main cause of delaying wound to heal. Nguyen et.al reported that active MMP9 was detected with high levels in the severe infected wounds based on the affinity resin and proteomic analysis (Nguyen et al. 2018). This suggests that targeting MMP9 is a promising strategy for diabetic ulcer treatment. In this study, Dolastatin 16 may act as the inhibitor of MMP9 as the result from network pharmacology analysis. Figure 10 illustrates the role of Dolastatin 16 as the inhibitor of MMP9.
Figure 10 The role of Dolastatin 16 in diabetic wound healing process. The figure was created in Biorender.com (https://www.biorender.com) (Nguyen et al. 2018; Chen et al. 2023)
As shown in Fig. 10, when skin injury occurs, a mass of neutrophils will be drawn to the site of injury and release cytokines like MMP9 and ROS that kills pathogens/bacteria and regulates thrombus formation. The high levels of ROS trigger the NF-κB activation, which increase the production and overexpression of MMP9 leading to the chronic wound (Nguyen et al. 2018). Dolastatin 16 may act as the inhibitor of MMP9 which accelerates the wound healing process in diabetic foot patients.
To validate the MMP9 as the potential target of Dolastatin 16, this study also conducted molecular docking and molecular dynamics studies. In the molecular docking simulation, we used MMP9 in the catalytic (PDB ID: 1GKC) and hemopexin (PDB ID: 1ITV) domain. In the catalytic site, the inhibitor was not directed to bind to the catalytic Zn2+ ion as the metal site is the most conserved features across all MMPs. Targeting the catalytic site of non-zinc binding may offer a promising strategy to obtain more selective MMP9 inhibitor (Jacobsen et al. 2010). According to Fig. 7, the result showed that MMP9-catalytic site has low binding affinity value, at -9.7 kcal/mol. The binding affinity of Dolastatin 16-MMP9 also lower than a reverse hydroxamate which is the native ligand (-6.6 kcal/mol) and (R)-ND-336 (-8.9 kcal/mol) indicating that Dolastatin 16 could be more potent than the inhibitor in 1GKC and (R)-ND-336. (R)-ND-336 is the MMP9-catalytic site inhibitor designed by Nguyen et.al. (R)-ND-336 exhibited greater efficacy over Becaplermin which is FDA-approved drug for diabetic foot ulcer treatment (Nguyen et al. 2018). Therefore, it was used as another ligand control in this study.
Figure 11 Molecular docking result. (a) Dolastatin 16 bound in the binding site (yellow cavity) of MMP9 (PDB ID: 1GKC), (b) 3D view of Dolastatin 16 binding to the binding pocket, hydrogen bonds represented as dashed lines, and (c) 2D interaction of Dolastatin 16 with MMP9 (PDB ID: 1GKC)
Figure 11 demonstrates that Dolastatin 16 bind to the binding site of MMP9 which involves some key residues. Four hydrogen bonds were observed in the binding site of MMP9-catalytic site including Leu187, Leu188, Ala189, and Pro421. The hydrogen bond of Pro421 has the closest distance (1.799 Å) to the Dolastatin 16 structure which might enhance the binding to MMP9 to exhibit lower binding affinity. Other important residues in the binding pocket of MMP9-catalytic site are Val398, His401, His405, and Tyr423 involved in the hydrophobic contact. The results indicate that Dolastatin 16 could strongly bound to the catalytic domain of MMP9 via the hydrogen bonds and hydrophobic interaction. The detail interaction of Dolastatin 16 with the residues of MMP9-catalytic site and the other protein targets are attached in the SI Table 5.
The structural behavior, molecular flexibility, and stability of MMP9 docked with Dolastatin 16 compound were examined by 50ns of MD simulation employing YASARA software. This study analyzed RMSD, RMSF, Radius of Gyration, and MM-PBSA binding free energy. For clarity, only RMSD and binding energy are shown here in the main text, while the other figures can be found in the Supporting Information. The changes in Dolastatin 16 complex were evaluated and compared between the different systems. The conformational changes of protein-ligand complexes during the simulation can be analyzed by RMSD (Dong and Wang 2023). Figure 8a represents the RMSD-Cα plots for native protein MMP9 and protein-ligand complexes. As compared to the native protein, the complexes exhibit some fluctuations, indicating protein intervention resulted from its interaction with the ligands. However, the overall patterns are quite similar throughout the simulation. In the initial 10 ns, the complexes rose to approximately 3 Å, followed by a decrease to around 2 Å during the midterm 30 ns of simulation. Subsequently, the graph continues to increase until the end of simulation process. The major change of MMP9 structure due to the ligand binding may suggest that it is important for inhibition.
This study also examined the RMSD value of ligand movement to understand the movement of the ligand from the binding site during the simulation. A higher RMSD value indicates that the ligand has moved further away from the binding site. As illustrated in Fig. 8b, the plot of Dolastatin 16 complex remained stable at around 4 Å for the initial 35 ns of simulation, after which it increased to approximately 6 Å until the end of the simulation. This trend indicates that the ligand exhibits a significant movement from its initial binding pocket position.
To further evaluate the most favorable bindings of the complexes, the binding free energy was also assessed using MM-PBSA method. The results of MM-PBSA binding energy for two complexes are plotted over time as illustrated in Fig. 8c. A more positive value indicates a stronger binding of the complex (Mitra and Dash 2018). Based on the analysis, the Dolastatin 16-MMP9 complex exhibited a higher average positive binding energy of 154.098 kJ/mol, suggesting a strong binding affinity. On the other hand, (R)-ND-336-MMP9 demonstrated a lower binding energy of 74.467 kJ/mol, indicating a weaker binding. Even though Dolastatin 16 ligand experienced a significant movement from the binding site (referred to the RMSD ligand movement value), it could maintain the strong affinity with the binding pocket during simulation. The detail visualization of the MMP9 in complex with Dolastatin 16 in the period of simulation is provided in Fig. 12.
Figure 12a) Overall structure of native MMP9 in the catalytic domain before MD simulation, b) Molecular interaction between Dolastatin 16 at the MMP9 binding site: 10ns (top), 30 ns (mid), and 50 ns (bottom) during molecular dynamic simulation
To understand the intervention of MMP9 due to the ligand binding, Fig. 12 presents the comparison between structure of native MMP9 pre and post MD simulation. The figure reveals that there has been a conformational change in the C-terminus and N-terminus structure affected from the ligand binding to the S1’ pocket of MMP9. The result also aligned with the RMSD-Cα pattern that shows some fluctuations in the Dolastatin 16-MMP9 complex. Moreover, the figure demonstrates that the complex of Dolastatin 16-MMP9 remained to interact with the residue keys of MMP9 through hydrogen and hydrophobic interaction until the simulation finished. The continuous contribution of hydrogen bond interactions to the binding mode indicates that Dolastatin 16 have stable interactions with MMP9 protein. Hydrophobic interactions also play a crucial role in inhibition and help stabilize the ligand in the active site. The outcome of the MD simulation shows that Dolastatin 16 is well accommodated with the MMP9 pocket over the simulation period. Moreover, through MD simulation, we can conclude that Dolastatin 16-MMP9 is more stable than (R)-ND-336-MMP9 complex, indicating more potential as the MMP9 inhibitor in diabetic foot condition.