Gene and protein interaction network analysis in the epithelial-mesenchymal transition of Hertwig's Epithelial Root Sheath reveals periodontal regenerative drug targets - An in silico study.

Background and aim: Hertwig's Εpithelial Root Sheath (HΕRS) has a major function in the developing tooth roots. Earlier research revealed that it undergoes epithelial-mesenchymal transition, a vital process for the morphogenesis and complete development of the tooth and its surrounding periodontium. Few studies have demonstrated the role of HERS in cementogenesis through ΕMΤ. The background of this in-silico system biology approach is to find a hub protein and gene involved in the EMT of HERS that may uncover novel insights in periodontal regenerative drug targets. Materials and methods: The protein and gene list involved in epithelial-mesenchymal transition were obtained from literature sources. The protein interaction was constructed using STRING software and the protein interaction network was analyzed. Molecular docking simulation checks the binding energy and stability of protein-ligand complex.
Results: Results revealed the hub gene to be DYRK1A(Hepcidin), and the ligand was identified as isoetharine. SΤRIΝG results showed a confidence cutoff of 0.9 in sensitivity analysis with a condensed protein interaction network. Overall, 98 nodes from 163 nodes of expected edges were found with an average node degree of 11.9. Docking results show binding energy of -4.70, and simulation results show an RMSD value of 5.6 Å at 50 ns.
Conclusion: Isoetharine could be a potential drug for periodontal regeneration.

Introduction

HΕRS is a film or membrane with two epithelial layers that develops from the dental or enamel organ (an aggregation of cells lying above the dental papilla, observed in a developing tooth). It has a primary role in the formation of the roots of a tooth, and guides in determines shape, size and number of roots (Zeichner-David et al., 2003). It also participates in formation of cementum by epithelial–mesenchymal transition (EMT) and various other processes ranging from structural (dental ectomesenchymal tissues that are subdivided into dental follicle and dental papilla), regulating the root development timing, mesenchymal cell inducing differentiation into precursors of cementoblasts and odontoblasts (Chen et al., 2014). Malassez cell epithelial cells rests have been found to harbor clonogenic epithelial stem-cell populations that match functional and phenotypic characteristics with mesenchymal stromal/stem cells. HERS cells are believed to go through EMT to transform into functional cementoblasts through the enamel proteins synthesis in HERS cells, after which morphologic and phenotypic alterations, and transcripts for multiple cementum-associated proteins are expressed (Thomas, 1995). Most EMT interactions follow signalling pathways that transmit interactions between epithelial and mesenchymal cells.

EMT is an active process in which epithelial cells are differentiated into a mesenchymal phenotype. This phenomenon occurs during embryonic development and is integral to wound healing, tissue regeneration, and organ fibrosis. However, this process can also contribute to tumor progression and metastatic expansion. The shift in cell activity and differentiation are mediated by SΝAIL, zinc finger E-box-binding (ZΕB) transcription factor and the basic helix-loop-helix (HLH) transcriptional regulator, which are the key factors for transcription whose functions are regulated at various levels including the translational, post-translational, and transcriptional levels (Lamouille et al., 2014). Defects in these levels can arrest skeletal and tooth development and cause defects in multiple organs (Puthiyaveetil et al., 2016).

During isolation proteins do not function. The protein interactions mediate all biological processes such as signalling transduction and metabolic processes. Their central role in biological functions and control mechanisms can lead an organism to healthy and diseased conditions in case of derangement (Gonzalez and Kann, 2012). Analyzing proteins and gene interactions can aid in understanding the disease process, driving the discovery of novel biomarkers, drug targets, and regenerative materials. Proteomic knowledge of physiopathology basis of the disease, associated diagnostic biomarkers, potential drug targets, and hubs (=proteins with larger connectivity in-network) (Singh et al., n.d.) can be translated and it has an application on effective diagnosis, treatment, and immunotherapeutic strategies. (Ghamari et al., 2015, Keyvani et al., 2016, Atan et al., 2014, Rakshit et al., 2014, Rezaei-Tavirani et al., n.d..; Safari-Alighiarloo et al., 2014). Systems-based therapeutics provide accurate approaches that has a concern on the underlying biology than the conventional reductionist approaches. Here, hub targets, which indicate clinical state, are networks of interacting molecular entities that are computationally derived. They incorporate measurements from the expression of molecules obtained from clinically useful biological interactions (MacLellan et al., n.d.). Docking techniques are an effective tool in drug design, allowing virtual screening of libraries for drug targets. Structure-based drug discovery used molecular docking as it can accurately predict the conformation of ligands to the target binding sites (Meng et al., 2011). Molecular interactions networks can be utilized in gene expression signature modeling (Chuang et al., 2007, Martinez-Ledesma et al., 2015, Schramm et al., 2013, Winter et al., n.d..).

This paper presents an in-silico system biology approach to find the hub protein and gene involved in the EMT of HERS that may aid in establishing novel insights in periodontal regenerative drug targets.

Materials and methods

A protein list (Lamouille et al., 2014) (Table 1) and gene list (Åberg et al., 1997, Luan et al., n.d..; Sonoyama et al., n.d.) (Table 2) thought to be involved in the EMT of the HΕRS was obtained from various literature sources. The Search Tool for the Retrieval of Interacting Proteins/Genes popularly known as ‘STRING’ (version11.0) database helped estimate the probable Protein-Protein interaction network. This platform produces possible interactions between protein networks using predictions, literature pertaining to genome-related analysis and also information generated at a massive or exponentially fast rate, (Szklarczyk et al., 2011, STRING Protein-Protein Interactions Network, n.d.). All protein characteristics are given by the Universal Protein (“Universal Protein Resource UniProt. Available online: https://www.uniprot.org/ (accessed on 16 June 2020).,” n.d.). A score for node-node interaction is calculated in STRING and the output is analysed (“STRING Score Computation. Available online: https://version10.string-db.org/help/faq/ (accessed on 22 April 2021),” n.d.). The genes involved in the EMT of HERS were analyzed using GENEMANIA 3.6.0 for co-expression, identifying Hub genes, physical interactions, genetic interactions, functional annotation roles. The top 20 genes were selected.

Table 1

EMT Transcription factors, targets, and signaling pathways.

TranscriptionFactors	Downregulated Expression During EMT	Upregulated Expression During EMT	Regulatory Signaling Pathways	References
SNAIL1 & SNAIL2	ECAD, claudins, occludin, Crumbs 3, PALS 1, PATJ, Cytoceratines, desmoplakin, and plakophilin	Fibronectin, NCAD, collagen, MMP15, MMP2, MMP9, TWIST,I□1,I□2, ZEB1 & ZEB2	TGFβ- SMA□3, WNT –β- Catenin, Notch, P13K-AKT, NF-κ B,EGF & FGF	[26–33]
TWIST 1	ECAD, claudins, occludin, desmoplakin & plakoglobin	Fibronectin, NCAD, & α5 integrin	MAPK	[26][27][34–36]
ZEB1 & ZEB2	ECAD, ZO1, Crumbs 3 and plakophilin	NCAD, MMPs	TGFβ- SMA□3, WNT –β- Catenin, & RAS-MAPK	[26][27][37–42]
FOXT3	Unknown	unknown	β 1 integrin & laminin	[43]
FOXC2	ECAD,	Fibronectin, Vimentin, NCAD & α SMA	TGFβ- SMA□3	[44]
FOXF1	ECAD, claudin1, occludin, desmoglein 1β, desmoglein 2, □esmocollin 2, desmoplakin	Fibroectin, NCAD	unknown	[45]
FOXQ1	ECAD	Fibronectin, NCAD& Vimentin	unknown	[46]
FOXO3A	ECAD	SNAIL1	AKT	[47]
FOXA1	ECAD	Fibronectin, Vimentin, and SNAIL1	TGFβ, HGF and AKT	[48]
FOXA2	ECAD and ZO1	Fibronectin, Vimentin, NCAD, SNAIL1 & SNAIL2	TGFβ, HGF and AKT	[48]
Serpent, GATA4, AN□ GATA A6	ECAD, Xrumbs and claudins	NCAD & MMP1	unknown	[49]
HMGA2	ECAD	SNAIL1, SNAIL2 & TWIST	TGFβ, SΜA□3	[50]
SΟX9	Unknown	SNAIL2	BΜPs and PKA	[51]
KLF8	ECAD	ΜΜP9	unknown	[52]
CBFA-KAP1	Unknown	FSP1	unknown	[53]
ZNF703 (=Zeppo 1 in mice)	ECAD	Vimentin, NCAD, SNAIL1, and 27Cytokeratin	RHO-GTPase	[54]
PRΧ1	ECAD	Vimentin & Laminin	BΜP2 and TGHβ	[55]

Table 2

List of genes involved in EMT of HERS.

Gene	Reference
Gli1, Nfic, Fgf	[56] [57]
Tgfβ, Bmp, PTHrP	[58–60]
Shh, Msx2, Dlx2	[61][57][62]
Smad4	[63][57]
Sox2, Msx2, Tgfbr2	[21][57] [64–66]
Axin2, Wnt10a, Dkk1, and Sfrp1	[67][68]
Fgf10, Ptch2, Shh-Gli1	[69][70]

Docking

Selecting and preparing protein structure

DYRK1A (Hepcidin) was the target protein in this study. The central regulatory molecule responsible for systemic iron homeostasis is Hepcidin. Inhibition of Hepcidin is a promising strategy for the treatment of anaemia of chronic disease. Using the (Protein-Data-Bank), PDB ID for this protein's 3D structure 6A1G was extracted (https://www.rcsb.org/). Protein-ligand docking was done by fixing side chains, adding charges, atom bumps, energy minimization, and taking help of the virtual-screening-software (PyRx). This was followed by transforming the docking polypeptide structure into the required PDBQT-format, to be easily read by AutoDock Vina in the PyRx program (Dallakyan and Olson, 2015, Trott and Olson, 2010).

Selecting and preparing ligands

FINDSITEcomb 2.0 has better accuracy in virtual screening than the docking method under challenging conditions when no templates with greater than 30% sequence ID to the target are present in the ligand-binding databases. If sequence identity cutoff was fixed to 95%, FINDSITEcomb will be able to give the unprecedented mean AUC = 0.90 using modeled structures.

Using this tool, the ligand chosen for this study is isoetharine. PubChem database was used to identify the 3D structure of the ligand in the Spatial Data File (SDF) formatting style.

Molecular docking and selection of binding site

In-silico molecular docking, the software scores every interaction between the binding affinity of a macromolecule and a group of small molecules. In this investigation, docking was done using an open-source virtual screening software called AutoDock Vina in PyRx, which has advanced features that accurately render and predict the binding mode speed. The Vina Wizard control in the PyRx program selects the ligand molecules and protein that has to be docked. Blind docking changes the dimensions of the grid that shows on the protein. To begin the docking operation, select the “Run Vina” control. In the “Analyze Results” tab, the results are seen and can be exported to the working directory in CSV format.

Molecular dynamics simulation

Molecular dynamics were simulated for a period of 50 ns on Desmond v4.4 (Desmond, Schrödinger LLC, NY, USA). Docking experiments form the preliminary phase for molecular dynamics simulation, examining proteins and ligand complexes. In static settings, ligand-binding state can be predicted by Molecular Docking.

The physiological environment's ligand binding status was estimated using simulations. The MPPW (Maestro’s-Protein-Preparation-Wizard) with enhancement and minimization of the complex, was used to pre-process protein-ligand complex in addition to SBT (System-Builder-Tool) which was employed for building up the process. The solvent model chosen was a ΤIIP3P (Transferrable-Intermolecular-Interaction-Potential-3-Points) having a crystalline structure with 3 mutually perpendicular axes which were unequal in length. Force-field used was OPLS 2005 for simulation, with the aim of representing potential energy of this complex in terms of the bond, angle, and torsion. Neutrality was maintained by the addition of counter ions whenever necessary. Physiological circumstances were replicated through the addition of 0.15 M of salt (sodium chloride). The NPT ensemble was set at three hundred kelvin & one atm pressure for simulation. Before simulation the model was made relaxed. After every 50 ns, the analysis was done by using the saved trajectories, and simulation's stability was observed by recording the RMSD of the ligand and protein over a period of time.

Results

Analysis of protein-protein interaction

Following STRING results with the sensitivity of the analysis with a confidence cutoff of 0.9, a more condensed protein interaction network was obtained. Overall, 98 nodes from 163 nodes of expected edges were found with an average node degree of 11.9 (Fig. 1). We found a PPI enrichment value of e-16 (p <.01), which indicates the small size of the current set of proteins. The node degree average was 11.9, and the local clustering coefficient average was 0.586. The casted network shows the possible PPI among the proteins involved in the EMT of HΕRS Fig. 2.

Fig. 1

Analysis of protein interaction.

Fig. 2

Binding energy scores isoetharine with Hepcidin.

Genemania

DYRK1A was identified as the hub gene through this analysis.

Molecular docking

Both hydrogen and hydrophobic bonding contributed significantly to ligand binding. The best ligand binding affinity is seen in the 9th conformation of the protein with the binding energy score of −4.70 (Table 3). The ligand formed the hydrogen bonding with residues 239A, 241A, 291A (Supplemental Table 3). The ligand position is stable in the pocket which has contributed to hydrophobic residues found around the site along with 170A of PHE residue in a distance of 3.92 of ligand atom of 3401. Supplemental Table 3 lists the scores of docking, hydrogen bonding, and hydrophobic information.

Table 3

Binding energy scores isoetharine with Hepcidin/Docking Conformation score of isoetharine with Hepcidin.

Rank Pattern	Sub Rank	Run	Binding Energy	Cluster RMSD	Reference RMSD
1	1	9	−4.7	0	34.23
2	1	8	−4.64	0	53.24
3	1	5	−4.61	0	36.35
4	1	1	−4.52	0	32.92
5	1	10	−4.28	0	49.63
6	1	4	−4.25	0	38.38
7	1	7	−3.94	0	41.82
8	1	3	−3.67	0	56.29
9	1	2	−3.56	0	34.18
10	1	6	−3.42	0	35.25

Molecular-Dynamic-Simulation (MDS) results

Root-Mean-Square-Deviation (RMSD)

Analysis of the stimulated trajectories was done in Desmond (Desmond, Schrödinger, NY, USA). The ligand-protein bonds, Root mean square deviation (RMSD) and Root mean square fluctuation (RMSF) were determined by the analysis of the Molecular-Dynamics-trajectory.

The RMSD values of ‘C-α’ atoms of the ligand-bound polypeptide structure vary over time (Fig. 3). The stability of the protein is evident at 5 ns in the 6A1G-isoetharine RMSD plot. After that point, the RMSD values remain constant within 0.5 Å, which is within the tolerable limit. The RMSD values of the protein to ligand fit fluctuates within 1 Å. These features are indicative of the ligand’s stability when bound to the receptor within the observed simulation period.

Fig. 3

Root means square deviation (RMSD) of the C-alpha atoms of protein and the ligand complex (6a1g-isoetharine) with time. The left Y-axis depicts protein RMSD variation through time. The right Y-axis indicates ligand RMSD variation through time.

A specific particle when displaced, from its stable and constant location over a period of time, is indicated by the root means square fluctuation (RMSF) value. Fig. 4 depicts the residue-wise RMSF value of the protein bound to the ligand. The higher peaks of the residues visible in Fig. 4 correspond to the loop regions appreciable in Fig. 3 depicting the mean square deviation trajectories (=N-terminal and C-terminal zones). The lower RMSF values of the remnants of the binding-location are indicative the steadiness of a particular protein-binding ligand Fig. 5.

Fig. 4

Residue-wise Root Mean Square Fluctuation (RMSF) of ligand-bound protein 6a1g.

Fig. 5

Protein Secondary Structure element distribution by residue index throughout the protein structure (6a1g). Alpha helices are denoted by red columns, and beta-strands are denoted by blue columns.

The hydrophobic interactions and hydrogen bonding determine the most critical interactions of the ligand-proteins, as listed in Table 3. LYS_188, GLU_203, and ASP_307 are important in terms of Hydrogen-bonds (Fig. 6). Fig. 6 displays the stacked bar charts that have been normalized throughout the trajectory. For instance, in Fig. 6, a value of 1.0 signifies that the specific interaction was maintained for the complete simulation time. In Fig. 6, values above 1.0 are present because residual proteins may have multiple contacts of the same subtype with the ligand.

Fig. 6

Protein-ligand contact histogram.

The gyration radius as a measure of protein structural compactness is shown in Fig. 7. It depicts the stability of the secondary structures that are packed compactly into 3D protein structure. If folded correctly, protein will have a reasonably constant Rg value. The graph indicates that the protein is compactly packed with 3D space for most of the simulation time.

Fig. 7

The radius of Gyration.

Discussion

The most challenging part of periodontal regeneration is the formation of acellular extrinsic fiber cementum and alveolar bone. Investigators have long believed that the epithelial membrane cells have both cementoblasts as well as mesenchyme-like traits in vivo (Sonoyama et al., 2007, Zeichner-David et al., 2003) and in vitro (Bosshardt and Nanci, 2004, Huang et al., 2009, Lézot et al., 2000). Previous research has identified several molecules like enamel matrix proteins with the interaction with the role of Hertwig's epithelial root sheath (HERS) enhances periodontal acellular cementum structure formation. (Cancedda et al., 2017, Chen et al., 2010). Strategies adopted from studying the cells could help in creating drugs that activate or mimic physiological healing and repair. The spatial and temporal components and signals of root formation and development could be harnessed to help in the repair and regeneration of defects. Despite the importance of periodontal regeneration, extensive research on the resurgence of acellular cementum has so far proven inconclusive. In this study, we chose in silico approach for in-depth analysis of EMT of HERS genes.

Network biology is an emerging field of study that focuses on molecular and genetic interactions, disease network biomarkers, and therapeutic target discovery (Sonawane et al., 2019). High-throughput data collection allows for the interrogation of molecular interactions to discern the determinants of system-wide behavior (Barabási and Oltvai, 2004). Periodontal ligament (PDL) and cementum grow to parallel with the fragmentation of HERS and the differentiation of dental follicle mesenchymal cells during tooth root development. However, the significance of events after HERS fragmentation in dental epithelial cells of the PDL remains ambiguous.

Studies on whether functional cementoblasts are formed from the EMT of HERS have so far been inconclusive (Luan et al., n.d.; Thomas, 1995). Enamel proteins, E-cadherin, cytokeratin, cemento associated proteins, osteo associated proteins, vimentin expressed by HERS cells all have the belief that functional cementoblasts are formed via EMT of HERS cells (Kaneko et al., n.d., Obara et al., 1999, Thomas, 1995). Li et al. demonstrated that HERS has the potential to differentiate into cementoblasts and generate cementum-like tissue in vivo (Li et al., 2019).

Signalling mechanisms that respond to extracellular cues begin and control gene expression regulation during EMT and non-transcriptional alterations. The present study obtained the protein and gene list of HERS involved in EMT from various literature sources. Transforming growth factor-β (TGFβ) family signalling has an important function in the convergence of signalling pathways important in EMT. TGFβ1 induces EMT of Hertwig’s Epithelial Root Sheath (HERS) cells via PI3K/AKT-pathway stimulation (Sonoyama et al., n.d.). In-vitro research has shown that epithelial cells derived from HERS can form mesenchymal cells through EMT under TGF-β signalling (Sonoyama et al., n.d.; Akimoto et al., 2011, Lee et al., 2014).

Iron balance is highly regulated by Hepcidin. Excess Hepcidin can cause blockage of iron and in chronic conditions lead to anemia. Guo et al. reported that AD-associated bone loss can be prevented by correcting the potential hepcidin-FPN-iron axis. (Guo et al., 2021). Peng Zhang et al. found that hepcidin levels did not affect bone loss and were overexpressed (with ovariectomy) in murine models and suggested that Hepcidin may possibly be employed to prevent osteoporotic changes after menopause (Zhou et al., n.d.).

Using the system biology approach, we identified the hub gene and protein hepcidin. To identify a new drug for this target protein, we used the drug repurposing approach. The methodology we used to determine the ligand is based on FINDSITEcomb2.0 (Zhou et al., n.d.). The molecular docking technique helps in identifying accurate inhibitors for a particular protein based on the scores obtained by the binding affinity of different conformations of the docked poses. Tools for visualization like PLIP helps to locate the ligands in the binding pocket and bonds exhibited with the neighboring residues.

Based on an unprecedented mean (AUC = 0.90) obtained from modeled structures, isoetharine was the choice of ligand. Isoetharine binds to beta-2 adrenergic receptors present in the bronchi smooth muscle and activates intracellular adenyl cyclase, which in turn catalyzes ATP to cAMP. The CADD method used the interaction between protein and ligand. Structurally related ligands can evoke a range of receptor conformations, allowing multiple signaling states and providing an avenue to engineer and influence pharmacotherapy pharmacotherapy (Kenakin and Miller, 2010, Luttrell and Kenakin, 2011). The molecular dynamic simulation results show binding energy of −4.70 with an RMSD value of 34.23, implying that isoetharine could be a potential drug for regenerating acellular extrinsic fiber cementum and bone by targeting HERS hub protein hepcidin.

This research forms the basis for a new model of biological research on HERS related to root development and regeneration. These findings provide new insights into cellular events that could be the basis for future tissue engineering studies. Considering the infancy of this type of research, the main caveat of our study is the lack of in vitro data regarding the therapeutic potential of isoetharine. There is scant research in this area of study on drug activation of gene and tissue regeneration. Major investment and further studies on key signaling pathways and the molecules involved could elucidate the mechanisms and regulators of tissue regeneration. Future research involving more hub genes like DYRK1A can be identified in junctional epithelial cells in inflammatory conditions will give insights into periodontal soft tissue regeneration.

Conclusions

In conclusion, this research presents evidence, discovered using network biology, regarding previously unidentified drug targets involved in complex gene interactions in a developing tooth. Network analysis of genes involved in epithelial-mesenchymal transition of Hertwig’s Εpithelial Root Sheath (HΕRS) revealed DYRK1A as hub gene. DYRK1A is a potential periodontal drug target. Isoetharine may be a potential drug for periodontal regeneration.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.