Multiple-Molecule Drug Design Based on Systems Biology Approaches and Deep Neural Network to Mitigate Human Skin Aging

Human skin aging is affected by various biological signaling pathways, microenvironment factors and epigenetic regulations. With the increasing demand for cosmetics and pharmaceuticals to prevent or reverse skin aging year by year, designing multiple-molecule drugs for mitigating skin aging is indispensable. In this study, we developed strategies for systems medicine design based on systems biology methods and deep neural networks. We constructed the candidate genomewide genetic and epigenetic network (GWGEN) via big database mining. After doing systems modeling and applying system identification, system order detection and principle network projection methods with real time-profile microarray data, we could obtain core signaling pathways and identify essential biomarkers based on the skin aging molecular progression mechanisms. Afterwards, we trained a deep neural network of drug–target interaction in advance and applied it to predict the potential candidate drugs based on our identified biomarkers. To narrow down the candidate drugs, we designed two filters considering drug regulation ability and drug sensitivity. With the proposed systems medicine design procedure, we not only shed the light on the skin aging molecular progression mechanisms but also suggested two multiple-molecule drugs for mitigating human skin aging from young adulthood to middle age and middle age to old age, respectively.


Introduction
Being the largest organ of the human body, the skin shows aging with biological age. Many people, especially female, like to spend money on cosmetics and pharmaceuticals regularly for preventing or reversing skin aging. Thus, changes in human skin caused by aging are important issues for both the pharmaceutical and cosmetic sectors worldwide [1]. Additionally, increasing life expectancy in developed countries reveals advancing age as the primary risk factor for numerous diseases [2]. Elder people tend to have dryness, itch, dyspigmentation, wrinkles, as well as benign and malignant tumors on skin. Under these worse conditions, they would feel sleep deprivation leading to having weakened immunity and getting infection. Hence, keeping our skin health promotes healthy aging [3]. Furthermore, identifying interventions, which are able to ameliorate skin aging progression, to delay, prevent or lessen age-related diseases is worth studying.
Human skin provides a primary protective barrier, routinely shielding us from allergens, microbes, and other environmental assaults, including solar ultraviolet (UV) irradiation, heat, infection, water loss, and injury. Skin aging is a complex process leading to the decrement of cutaneous functions and structures with time. Impaired epidermal barrier function, decline in resistance to infections and regenerative potential, and impairment of mechanical properties like loss of extensibility and elasticity are the essential biomarkers

Results
For the purpose of analyzing molecular progression mechanisms in human skin aging, extracting core signaling pathways from each core GWGEN becomes an essential issue. We defined three skin aging stages, including young-adult, middle-aged, and elderly human skin as shown in Figure 1. The research flowchart in Figure 2 shows how to construct the candidate GWGEN, real GWGENs, and core GWGENs so as to extract core singling pathways and investigate molecular progression mechanisms of human skin aging. By big database mining, the candidate GWGEN containing candidate PPIN and candidate GRN was constructed. With the help of the corresponding young-adult, middle-aged, and elderly skin microarray data, we applied system identification and system order detection methods to the candidate GWGEN for obtaining real GWGENs shown in Figures S1-S3 (Supplementary Materials), respectively. The Table 1 shows the number of nodes (e.g., proteins, TFs, miRNAs, and lncRNAs) as well as the edges standing for the interaction or regulation between two nodes for the candidate GWGEN and real GWGENs. According to Table 1, compared the nodes in candidate GWGEN to the nodes in real GWGENs, one could realize that the number of nodes diminish a lot in real GWGENs, reflecting that the false positives were removed successfully by the system order detection scheme. Since the real GWGENs were still too complicated to investigate molecular progression mechanisms of human skin aging, we applied principal network projection (PNP) method and selected the top-ranked 4000 nodes with significant projection values that could reflect 85% of the real GWGENs in three stages of skin aging to obtain core GWGENs (Figure 3a-c), respectively. In addition, for the genes in core GWGENs, we used the Database for Annotation, Visualization and Integrated Discovery (DAVID) Bioinformatics Resources version 6.8 to perform enrichment analyses for each stage of skin aging as shown in Tables S1-S3, respectively. Moreover, for investigating molecular progression mechanisms of skin aging conveniently, we denoted differential core signaling Molecules 2021, 26, 3178 4 of 28 pathways for young adult to middle-aged and middle-aged to elderly skin aging in respect of KEGG pathways, respectively. Based on skin aging molecular progression mechanisms, we identified essential biomarkers as drug targets for young-adult to middle-aged and middle-aged to elderly skin aging, respectively. Exploring candidate drugs toward our identified biomarkers, we trained a deep neural network of drug-target interaction in advance. We applied the trained model to predict the candidate drugs holding higher interaction probability with identified biomarkers. In order to narrow down the candidate drugs, we design two filters considering drug regulation ability and drug sensitivity. The more details will be discussed in the following sections. version 6.8 to perform enrichment analyses for each stage of skin aging as shown in Tables S1-S3, respectively. Moreover, for investigating molecular progression mechanisms of skin aging conveniently, we denoted differential core signaling pathways for young adult to middle-aged and middle-aged to elderly skin aging in respect of KEGG pathways, respectively. Based on skin aging molecular progression mechanisms, we identified essential biomarkers as drug targets for young-adult to middle-aged and middle-aged to elderly skin aging, respectively. Exploring candidate drugs toward our identified biomarkers, we trained a deep neural network of drug-target interaction in advance. We applied the trained model to predict the candidate drugs holding higher interaction probability with identified biomarkers. In order to narrow down the candidate drugs, we design two filters considering drug regulation ability and drug sensitivity. The more details will be discussed in the following sections.    Reseasrch flowchart of systems medicine design for human skin aging. Flowchart of using systems biology methods to construct the candidate GWGEN, real GWGENs, core GWGENs, and core signaling pathways to find skin aging progression mechanisms for identifying essential biomarkers. After obtaining the essential biomarkers, we applied trained a deep neural network of drug-target interactions to predict the potential candidate drugs holding higher probability. To narrow down the candidate drugs, we considered drug regulation ability by querying the CMap dataset and drug sensitivity by referring to the sensitivity dataset from DepMap portal. Consequently, we proposed two multiple-molecule drugs to mitigate the skin aging from young-adult to middle-aged and middle-aged to elderly-stage.

Figure 2.
Reseasrch flowchart of systems medicine design for human skin aging. Flowchart of using systems biology methods to construct the candidate GWGEN, real GWGENs, core GWGENs, and core signaling pathways to find skin aging progression mechanisms for identifying essential biomarkers. After obtaining the essential biomarkers, we applied trained a deep neural network of drug-target interactions to predict the potential candidate drugs holding higher probability. To narrow down the candidate drugs, we considered drug regulation ability by querying the CMap dataset and drug sensitivity by referring to the sensitivity dataset from DepMap portal. Consequently, we proposed two multiple-molecule drugs to mitigate the skin aging from young-adult to middle-aged and middle-aged to elderly-stage.

Differential Core Signaling Pathways from Young-Adult to Middle-Aged Skin Aging
The differential core signaling pathways from young-adult to middle-aged human skin were selected and analyzed as shown in Figure 4. According to our results, in core signaling pathways of young-adult skin aging only, receptor ESR1 receives microenvironment factor FASN to activate the TF SIRT6 through signaling transduction proteins PRR4 and LMNA. The TF SIRT6 could not only downregulate target gene RBBP8, which was modified by deacetylation, but also activate TF PARP1 to upregulate target gene XRCC1 to promote cell proliferation and DNA repair in young-adult skin aging only. The receptor ESR1 also regulates TF JUN through signaling transduction proteins GOT1 and CHEK2 to regulate target gene BRCA2 and KPNA2. TF JUN not only activates the target gene BRCA2, which was modified by phosphorylation, to promote DNA repair, but also downregulates the target gene KPNA2 to promote cell proliferation, DNA repair and cell cycle in young-adult skin aging only.
to promote cell proliferation and DNA repair in young-adult skin aging only. The receptor ESR1 also regulates TF JUN through signaling transduction proteins GOT1 and CHEK2 to regulate target gene BRCA2 and KPNA2. TF JUN not only activates the target gene BRCA2, which was modified by phosphorylation, to promote DNA repair, but also downregulates the target gene KPNA2 to promote cell proliferation, DNA repair and cell cycle in young-adult skin aging only. Figure 4. The core signaling pathways from young-adult to middle-aged skin aging. The green dot lines denote the signaling pathways in young-adult skin; the green dash lines represent the signaling pathways in middle-aged skin; the green solid lines indicate the signaling pathways in both stage; the green lines with arrow heads are upregulation (positive regulation); the green lines with circular heads are downregulation (negative regulation); the black solid lines with arrow heads mean activating cellular function; the black solid lines with circular heads mean inhibiting cellular function; the selected red target gene nodes indicate a higher gene expression in middle-aged skin . The core signaling pathways from young-adult to middle-aged skin aging. The green dot lines denote the signaling pathways in young-adult skin; the green dash lines represent the signaling pathways in middle-aged skin; the green solid lines indicate the signaling pathways in both stage; the green lines with arrow heads are upregulation (positive regulation); the green lines with circular heads are downregulation (negative regulation); the black solid lines with arrow heads mean activating cellular function; the black solid lines with circular heads mean inhibiting cellular function; the selected red target gene nodes indicate a higher gene expression in middle-aged skin compared with young-adult skin; the selected green target gene nodes indicate a lower gene expression in middle-aged skin compared with young-adult skin; the blue background shows young-adult skin; the brown background shows the overlap between young-adult and middle-aged skin; the skin color background shows middle-aged skin.
Next, in the core signaling pathways of both young-adult and middle-aged skin aging in Figure 4, the microenvironment factor FASLG was received by receptor FAS to activate TF E2F7 via signaling transduction proteins DAXX, MAP3K5 and AIFM1 to downregulate target gene APAF1 to promote apoptosis and cell-cycle arrest. In the next pathway, the receptor CXCR2 receive microenvironment factor CXCL1 to regulate TF E2F7, TF JUN, TF FOXM1 and miRNA MIR26B. First, the TF JUN was activated through signaling transduction proteins CENPJ and CDKN1A to activate TF FOXM1. As TF FOXM1, which was modified by phosphorylation, was activated, the target gene CAT was upregulated to inhibit ROS accumulation in both young-adult and middle-aged skin aging. The miRNA MIR26B was activated via signaling transduction proteins CENPJ, CDKN1A and AKT1 to inhibit target gene KPNA2 so as to promote cell proliferation. It is noted that, protein AKT1 was modified by phosphorylation. The TF E2F7 was also regulated by protein AIFM1 to downregulate target gene APAF1 in both young-adult and middle-aged skin aging.
In the core signaling pathways of middle-aged skin aging only, the TF FOXM1 was also activated through signaling transduction proteins CENPJ, CDKN1A, CDK4 and DCDC2 when receptor CXCR2 received microenvironment factor CXCL1 in middle-aged skin aging only as shown in Figure 4. With the activation of FOXM1, the target gene CCNB1 was upregulated to promote cell-cycle arrest and inhibit ROS accumulation. For the next pathway in middle-aged skin aging only, the microenvironment factor IGF1 was received by receptor IGF1R to regulate TF FOXO3 via signaling transduction proteins HMGCS2, ARRB1, PDK1 and AKT1. Additionally, the protein AKT1 and TF FOXO3 were modified by phosphorylation. The TF FOXO3 downregulates not only target gene SENS3 to promote DNA damage, apoptosis, and ROS accumulation, but also target gene GADD45A to promote DNA damage, apoptosis, and ROS accumulation and inhibit cell-cycle arrest in middle-aged skin aging only.
In summary, when the young-adult skin aging turned into middle-aged skin aging, DNA repair ability decreases and cell cycle starts to be arrested. Thereby, ROS accumulation increases and further promotes DNA damage and apoptosis in skin cells. Additionally, these molecular progression mechanisms from young-adult to middle-aged might potentially accelerate skin aging process in elderly skin aging. According to the core signaling analyses results and considering the overlap nodes between the GenAge and CMap datasets, we identified AIFM1, CAT, IGF1R, and LMNA as essential biomarkers for preventing skin aging from young adulthood to middle-age.

Differential Core Signaling Pathways from Middle-Aged to Elderly Skin Aging
The molecular progression mechanism based on differential core pathways from middle-aged skin to elderly skin aging is represented in Figure 5. In core pathways of middle-aged skin aging only, the TNF receptor superfamily member 1 alpha TNFRSF1A received microenvironment factor TNF to activate TF GATA2 through transduction proteins GABPA and STAT1 to upregulate target gene MMP9 so as to inhibit collagen stability and skin homeostasis in middle-aged skin aging only. Note that STAT1 was modified by phosphorylation. In the next pathway, the microenvironment factor NGF was accepted by neurotrophic receptor tyrosine kinase1 NTRK1 and then transmitted the signal through transduction proteins EME1, HSPB1, NEDD9 and CPNE2 to activate TF ETS1. TF ETS1 could downregulate target gene ERRFI1, which was modified by hypermethylation, through activating miRNA MIR573 to promote homeostasis in middle-aged skin aging only.
Next, we focus on the core pathways of both middle-aged and elderly skin aging. In the first pathway, the receptor NTRK1 receives the microenvironment factor NGF and then transmits the signal through transduction proteins KPNA2, KAT5, CST2 and HRAS to activate TF GATA2. The protein KAT5 was modified by phosphorylation. After GATA2 was activated, target gene BCL2 was upregulated to inhibit cell-cycle arrest and apoptosis in both middle-aged and elderly skin aging. For the second pathway, the receptor KIT could interact with the microenvironment factor KITLG to trigger TF AR through signaling transduction proteins FAM83H, HSPB1, PAX3 and H2AFB2. TF AR downregulated not only target gene TYR through triggering TF MITF to promote melanin synthesis, but also target gene CDH1, which was modified by phosphorylation, to promote cell-cycle, apoptosis and DNA damage in both middle-aged and elderly skin aging. In the third pathway, receptor LRP1 could receive the microenvironment factor CYR61 to trigger the TF ETS1 via signaling transduction proteins RAMP1, MCM2, GEMIN4 and NOP56. In both middle-aged and elderly skin aging, the TF ETS1 could negatively regulate target gene COL17A1 to promote melanin synthesis and inhibit collagen stability and skin homeostasis.
FOR PEER REVIEW 9 of 30 Figure 5. The core signaling pathways are obtained by projecting core GWGENs to KEGG pathways to investigate the aging progression mechanism from middle-aged to elderly skin aging. The green dotted lines denote the signaling pathways in middle-aged skin; the green dashed lines represent the signaling pathways in elderly skin; the green solid lines indicate the signaling pathways in both stages; the green lines with arrow heads are upregulation; the green lines with circle heads are downregulation; the black solid lines with arrow heads mean activating cellular function; the black solid lines with circle heads mean inhibiting cellular function; the selected red target gene nodes indicate a higher gene expression in elderly skin compared with middle-aged skin; the selected green target gene nodes indicate a lower gene expression in elderly skin compared with middle-aged skin; the blue background shows middle-aged skin; the brown background covers the overlap between middle-aged and elderly skin; the skin color background shows elderly skin.
Next, we focus on the core pathways of both middle-aged and elderly skin aging. In the first pathway, the receptor NTRK1 receives the microenvironment factor NGF and Figure 5. The core signaling pathways are obtained by projecting core GWGENs to KEGG pathways to investigate the aging progression mechanism from middle-aged to elderly skin aging. The green dotted lines denote the signaling pathways in middle-aged skin; the green dashed lines represent the signaling pathways in elderly skin; the green solid lines indicate the signaling pathways in both stages; the green lines with arrow heads are upregulation; the green lines with circle heads are downregulation; the black solid lines with arrow heads mean activating cellular function; the black solid lines with circle heads mean inhibiting cellular function; the selected red target gene nodes indicate a higher gene expression in elderly skin compared with middle-aged skin; the selected green target gene nodes indicate a lower gene expression in elderly skin compared with middle-aged skin; the blue background shows middle-aged skin; the brown background covers the overlap between middle-aged and elderly skin; the skin color background shows elderly skin.
In the core pathways of elderly skin aging only in Figure 5, CYR61/LRP1 could also trigger TF NDUFS4 through signaling transduction proteins CPNE2, MYH9 and ERCC6. The activated TF NDUFS4 might downregulate target gene CASP3 to inhibit apoptosis and promote DNA damage. For another pathway, the microenvironment factor IL6 was accepted by receptor IL6R to trigger TF YAP1 through signaling transduction proteins RHOB and CDK20. In the elderly skin aging only, TF YAP1 activated target gene CDC5L through inhibiting MIR126 to inhibit cell cycle and promote skin homeostasis and DNA damage.
In summary, for skin aging molecular progression mechanisms from middle age to old age, we found that the promotion of cell cycle process, the inhibition of apoptosis, and the damage of DNA arose in elderly skin. Furthermore, skin homeostasis and collagen stability were destroyed to cause lower immunity and epidermal thinning, that is, the increment of wrinkles. According to core signaling analyses and considering the overlap nodes between the GenAge and CMap datasets, we identified MMP9, IL6, BCL2, and CASP3 as essential biomarkers for preventing skin aging from middle age to old age. Moreover, by extracting differential core signaling pathways from young-adult to elderly skin aging, some cellular dysfunctions including proliferation, DNA repair and damage, cell-cycle arrest, apoptosis, ROS accumulation, collagen stability, skin homeostasis, and melanin synthesis are induced in the skin aging process shown in Figure 6.

The Application of Deep Neural Network of Drug-Target Interaction Prediction and the Design of Two Filters Considering Drug Regulation Ability and Drug Sensitivity
To explore the drug-target interaction toward our identified biomarkers, we trained a deep neural network for drug-target interaction prediction. The design framework is shown in Figure S4. The interaction dataset used for training are from BindingDB [35]. In total, there are 80,291 known drug-target interactions between 38,015 drugs and 7292 proteins. The number of unknown drug-target interactions is 19,966,109, which is greater than the known drug-target interactions. Considering the class imbalance problem, we randomly chose the number of unknown interactions and made them the same size as known interactions. We trained the model using 70% of data, including 10% of data as the validation set. The remaining 30% of data was used as the testing set. To the data preprocessing before training the model, we performed feature scaling by standardization. Assisted with principal component analyses (PCA) for dimensionality reduction, we had 1000 out of 1359 features. For the architecture of deep neural network of drug-target interaction, we used Adam as an optimizer (learning rate = 0.003) with binary cross-entropy loss. The input layer had 1000 neurons followed by 512, 256, 128, and 64 neurons of hidden layers, respectively. The output layer has one neuron. Except for using sigmoid function to the output layer, we set a nonlinear activation function ReLU for each hidden layer. Moreover, the dropout 0.5, 0.4, 0.3, and 0.1 was applied to each hidden layer, respectively. For the trained deep neural network of drug-target interaction prediction, the training accuracy, validation accuracy, and testing accuracy were 95.469%, 93.230%, and 93.077%, respectively. From the perspective of the deep neural network framework application, we used it to predict the potential candidate drugs for our identified biomarkers. When the score of candidate drug approaches one, it would be selected. In other words, the higher the score, the higher probability of interacting between the candidate drug and the identified biomarker. The upper horizontal part is the genetic and epigenetic progression mechanism from young-adult skin to middle-aged skin; the middle part indicates the genetic and epigenetic progression mechanisms from middle-aged skin to elderly skin; the red rectangle with orange background represents cellular functions; the yellow ellipse circles are microenvironment factors; the red dash lines surround the pathways and biomarkers that appear in two consecutive stages of skin; the black arrow lines represent the protein-protein interaction or transcriptional regulation; the black lines with circle head represent inhibit or downregulation; the red arrow lines represent the genes to induce cellular function; the red lines with circle head represent the genes to repress cellular function. The upper horizontal part is the genetic and epigenetic progression mechanism from young-adult skin to middle-aged skin; the middle part indicates the genetic and epigenetic progression mechanisms from middle-aged skin to elderly skin; the red rectangle with orange background represents cellular functions; the yellow ellipse circles are microenvironment factors; the red dash lines surround the pathways and biomarkers that appear in two consecutive stages of skin; the black arrow lines represent the protein-protein interaction or transcriptional regulation; the black lines with circle head represent inhibit or downregulation; the red arrow lines represent the genes to induce cellular function; the red lines with circle head represent the genes to repress cellular function.

The Application of Deep Neural Network of Drug-Target Interaction Prediction and the
In order to narrow down the candidate drugs predicted by the deep neural network framework based on the identified biomarkers, we designed two filters considering drug regulation ability and drug sensitivity. With the help of the CMap dataset, we could know whether a gene was upregulated or downregulated after treating the small molecule compound. The abnormal up or down gene expression could be found by comparing the gene expression of identified biomarkers to the later skin stage. The goal for the first filter is to select candidate drugs, which could reverse the abnormal gene expression. Afterwards, we used the drug sensitivity dataset (PRISM Repurposing Primary Screen) to consider drug sensitivity. The second filter aims to find the drugs with around zero values implying that they would not influence the cell line too much since we are not going to kill or proliferate cells toward the skin corresponding cell line. Consequently, we proposed niridazole, liothyroninr, decitabine, pinacidil, and allantoin as a multiple-molecule drug for mitigating skin aging from young adulthood to middle age; and allantoin, diclofenac, mepyramime, resveratrol, and azathioprine as multiple-molecule drug for mitigating skin aging from middle age to old age. The drug targets with their corresponding drugs are shown in Tables S4 and S5.

Investigating Skin Aging Molecular Progression Mechanisms by Differential Core Signaling Pathways from Young-Adult to Middle-Aged Human Skin Aging
In the first core pathway of young-adult skin aging only as shown in Figure 4, microenvironment factor fatty acid synthase FASN can promote cell proliferation, DNA repair, and cell-cycle arrest and interact with receptor ESR1 via crucial signaling transduction proteins PRR4, LMNA, GOT1, and CHEK2 to regulate TFs SIRT6, PARP1, and JUN. The signaling protein LMNA, which is an endogenous activator of TF SIRT6, could promote SIRT6-mediated downstream functions upon DNA damage. Moreover, protein LMNA could directly bind and activate TF SIRT6 toward histone deacetylation [36]. The TF SIRT6, which could control the longevity and regulation of DNA repair, could promote DNA repair and cell proliferation through the downregulation of the target gene RBBP8, which was mediated by deacetylation [37]. TF SIRT6 could also promote DNA repair under oxidative stress by activating TF PARP1 to upregulate target gene XRCC1 [38]. TF PARP1 serves as a genomic caretaker by participating in several molecular mechanisms such as DNA repair and cell-cycle regulation. Therefore, PARP1 was considered as a longevity assurance and aging-promoting factor [39]. The target gene XRCC1 upregulated by PARP1 was required for the viable and efficient repair for DNA single-strand breaks [40]. The TF JUN was also activated by the signaling transduction proteins GOT1 and CHEK2. The signaling transduction protein CHEK2 initiated by oxidative stress could regulate target gene BRCA2 and KPNA2 through interacting with TF JUN. In human cell, the serine kinase CHEK2 could induce the appropriate cellular response such as cell cycle checkpoint activation and DNA repair depending on the extent of damage, the cell type, and other factor. CHEK2 could participate in DNA repair by phosphorylating the target gene BRCA2 through TF JUN [41]. The karyopherin alpha2 KPNA2 expression had been reported to be induced in various proliferative skin disorders such as psoriasis and squamous cell carcinoma [42]. When the target gene KPNA2 was downregulated by TF JUN, cell proliferation, cell cycle and DNA repair induced by CHEK2 were indirectly promoted.
In the core pathways of both young-adult and middle-aged skin aging, the microenvironment factor FASLG binds receptor FAS to regulate TF E2F7 through signaling transduction proteins DAXX, MAP3K5 and AIFM1. Responding to ROS, the microenvironment factor FASLG was activated, then binding to death receptor FAS to promote apoptosis pathway [43]. Signaling transduction protein MAP3K5 known as apoptosis signal-regulating kinase 1 (ASK1) could respond to oxidative stress and be activated [44]. Target gene APAF1 is the core of the apoptosome, was activated by TF E2F7 to trigger the mitochondrial apoptotic pathway. Furthermore, target gene APAF1 was also involved in the maintenance of genomic stability by the cell-cycle arrest response elicited upon DNA damage and promoted apoptosis [45]. For other pathways in both young-adult and middle-aged skin aging, microenvironment factor CXCL1 bound the G-protein coupled receptor CXCR2 to activate signaling transduction protein CDKN1A through protein CENPJ. Protein CDKN1A known as CIP1, was a potent cyclin-dependent kinase inhibitor to regulate TFs E2F7, JUN, FOXM1, and miRNA MIR26B. First, protein CDKN1A transmits signal to AIFM1 so as to activate TF E2F7 to enhance apoptosis and cell-cycle arrest. TF JUN was also activated by protein CDNK1A to regulate FOXM1. Then TF FOXM1 upregulated target gene CAT, which was known as ROS detoxification enzyme and could defend the ROS accumulation [46]. After miRNA MIR26B inhibited the target gene KPNA2, the cell proliferation could be promoted [47].
In core pathways of middle-aged skin aging only, microenvironment factor CXCL1 also could regulate TF FOXM1 through signaling transduction proteins CENNPJ, CDKN1A, CDK4, and DCDC2 to trigger target gene CCNB1 as shown in Figure 4. Cyclin dependent kinase 4 CDK4, which was modified by phosphorylation, is a positive regulator of cell cycle entry and can stabilize and activate FOXM1, thereby promote cell cycle and suppress the levels of reactive oxygen species [48]. TF FOXM1 also had been reported to be essential for proper cell cycle progression via activating cell cycle gene CCNB1 for propelling specific cell cycle phase and inhibition ROS accumulation [49]. For another pathway, the microenvironment factor IGF1 was received by receptor IGF1R to activate FOXO3 via signaling transduction proteins HMGCS2, ARRB1, PDK1 and AKT1. The protein AKT1, which was modified by phosphorylation, could activate TF FOXO3. Protein phosphoinositide-dependent kinase PDK1 was one of the upstream kinases that activate AKT1. After AKT1, which is a key regulator of the PI3K/AKT1 signaling cascade controlling cell growth and survival, was activated and modified by phosphorylation, TF FOXO3 would be activated [50]. Moreover, it had been reported that the enhanced ROS production might further activate the signal of PI3K/AKT pathway, thus establishing a self-perpetuating cycle leading to further aging [51]. TF FOXO3, which was modified by phosphorylation, could downregulate target genes SESN3 and GADD45A [52]. TF FOXO3 could decline ROS rescue pathway through downregulating the peroxiredoxin gene SESN3, which is responsible for the biphasic ROS accumulation. Therefore, FOXO3-induced ROS was increased and then accelerated for apoptosis and DNA damage [53]. Furthermore, phosphorylated FOXO3 also inhibited proapoptotic activity such as cell-cycle arrest by downregulating GADD45A [54]. The cause of the pleiotropic action of GADD45 members, a decreased inducibility, might lead to far-reaching consequences such as DNA damage accumulation and disorder of cellular homeostasis and could eventually contribute to the aging process [55]. Therefore, we suggest that the downregulation of GADD45A also promotes ROS accumulation through cell-cycle arrest and the inhibition of proapoptotic activity.

Investigating Skin Aging Molecular Progression Mechanisms by Differential Core Signaling Pathways from Middle-Aged to Elderly Human Skin Aging
According to the core pathways of middle-aged skin aging only in Figure 5, the ligand TNF can inhibit collagen stability and skin homeostasis through receptor TNFRSF1A by transmitting the signal through significant signaling transduction proteins GABPA and STAT1 to TF GATA2. The proinflammatory cytokine tumor necrosis factor-alpha (TNF-A) inhibits collagen synthesis and enhances collagen degradation via increasing the production of target gene MMP9. It also increases the risk of cutaneous infections in the elderly by reducing skin immunity [56]. The activation of STAT1 is modified by phosphorylation. STAT1 has also been indicated as a potential target in the treatment of psoriasis, which is a chronic skin diseases [57]. TF GATA2 could upregulate target gene MMP9 to digest collagen type IV, which is an important component of the basement membrane in skin [58].
For the next pathway, ligand NGF can promote skin homeostasis through receptor NTRK1 to transmit the significant signal via signaling transduction proteins EME1, HSPB1, NEDD9, and CPNE2 to upregulate TF ETS1. In human skin, proliferating keratinocytes release NGF in an increasing amount. Receptor NTRK1, known as tyrosine kinase receptor (TrkA) is the high-affinity receptor for NGF. At the skin level, NTRK1 could mediate NGFinduced keratinocyte proliferation [59]. Note that protein NEDD9 could be modified by phosphorylation in human skin [60]. TF ETS1, which was regulated through the signaling pathway activated by ligand NGF, has been identified to be associate with skin aging [60]. The expression levels of MIR-573 were found to be lower in melanoma tissues and cell lines when compared to normal skin tissue. Moreover, MIR-573 reduction was demonstrated to be essential in melanoma initiation and progression [61]. Target gene ERRFI1, which was modified by hypermethylation, is required for proper epidermal homeostasis [62,63].
Focusing on core pathways in both middle-aged and elderly skin aging in Figure 5, the ligand NGF inhibits apoptosis and promotes cell-cycle arrest when received by receptor NTRK1 to activate TF GATA2 via signaling transduction proteins KPNA2, KAT5, CST2, and HRAS. Protein KAT5, which was modified by phosphorylation, has been presumed to serve as a potential biomarker for melanoma therapeutic target [64]. NGF can not only rescue human epidermal keratinocytes from spontaneous and UVB-induced apoptosis via NTRK1, but also protect keratinocytes from cell death via target gene BCL2 family of apoptosis inhibitors [59]. Antiapoptotic function of target gene BCL2 is regulated by phosphorylation. In addition, target gene BCL2 could not only regulate cell cycle progression, but also act as an antioxidant that may regulate intracellular ROS. Expression of target gene BCL2 has been observed to increase upon the induction of a senescence-like growth arrest or apoptosis by oxidative stress [65,66].
In the next core pathway of both middle-aged and elderly skin aging, the ligand KITLG promotes melanin synthesis, DNA damage, and inhibits cell-cycle arrest by modulating TFs AR and MITF via signaling transduction proteins FAM83H, HSPB1, PAX3, ATF5, and H2AFB2. The tyrosine kinase receptor KIT, its ligand KITLG, and TF MITF have been reported to play an important role of initiating and regulating signaling systems and transcription factors of melanin production. TF MITF also regulates melanocyte pigmentation by inducing target gene TYR [67]. Moreover, a previous study supposed that PAX3 and SOX10 could act together to induce the expression of MITF [68]. Target gene CDH1, which is downregulated by TF AR, has been reported to be regulated by phosphorylation [69]. It has been reported that cells lacking target gene CDH1 have a shortened G1 phase, accumulate DNA damage, and undergo apoptosis [70].
In the final core pathways of both middle-aged and elderly skin aging, the ligand CYR61 could modify skin homeostasis and melanin synthesis through receptor LRP1 to transmit signal by signaling transduction proteins RAMP1, MCM2, GEMIN4 and NOP56 for upregulating TF ETS1. Responding to oxidative stress, CYR61 was elevated in the dermis of chronologically aged human skin, promoting aberrant collagen homeostasis by downregulating collagen members, the major structural protein in skin, to promote collagen degradation [71,72]. The loss of target gene COL17A1 and MCM2 expression in advanced aged skin has been found to eventually cause epidermal thinning [73].
Focusing on the first core pathway of elderly skin aging only in Figure 5, through the signaling transduction starting from the ligand CYR61, TF NDUFS4 can promote apoptosis and DNA damage through signaling transduction proteins CPNE2, MYH9 and ERRC6. The ligand CYR61 interacting with receptor LRP1 has also been indicated to contribute to CCN1-induced ROS accumulation and CCN1/TNFA-induced apoptosis [74]. With the downregulating target gene CASP3 by TF NDUFS4, senescence fibroblast can resist apoptosis death [75].
In the final core pathway of elderly skin aging only, the ligand IL6 could be accepted by receptor IL6R and then the significant signal is transmitted through signaling transduction proteins RHOB and CDK20 to activate TF YAP1. Proinflammatory cytokine IL6 has been suggested to be a biomarker of health status in the elderly [76]. TF YAP1 has been identified to play a physiological role in skin homeostasis, which can promote cell proliferation in the basal layer [77]. Knockdown of target gene CDC5L induces mitotic arrest and DNA damage [78].

The Genetic and Epigenetic Molecular Progression Mechanisms from Young-Adult to Elderly Human Skin Aging
The overview of overall skin aging molecular progression mechanisms is shown in Figure 6. Microenvironments trigger corresponding ligand signals to initiate some cellular dysfunctions affecting skin aging progression. Thus, core signaling pathways with the genetic and epigenetic modifications play a significant role in cellular dysfunctions of signaling transductions for each stage of skin aging.
In Figure 6, the core pathways of young-adult skin aging only, ligand FASN (oxidative stress) binds to receptor ESR1 to mediate two pathways. Responding to DNA damage signal to cause of oxidative stress, LMNA directly binds and activates TF SIRT6 toward histone deacetylation [36]. Activated SIRT6 promotes DNA repair cell-cycle and proliferation through the downregulating gene RBBP8, which was modified by deacetylation [37]. In addition, TF SIRT6 also promotes DNA repair and cell cycle under oxidative stress by activating TF PARP1 to upregulate target gene XRCC1 [38]. Transduction protein CHEK2 also responds to oxidative stress from ligand FASN to activate TF JUN. TF JUN promotes DNA repair through phosphorylating target gene BRCA2 [41]. Moreover, TF JUN could downregulate target gene KPNA2 to promote cell proliferation, cell cycle and DNA repair [42].
In the core pathways of both young-adult and middle-aged skin aging, ligand CXCL1 (oxidative stress) binds to receptor CXCR2 to trigger protein CDKN1A. Activated protein CDKN1A not only upregulates target gene CAT to defend ROS accumulation through TF JUN and FOXM1, but also inhibits target gene KPNA2 to promote cell proliferation by activating MIR26B [46,47]. Next, responding to ROS induced from DNA damage, the ligand FASLG interacts with FAS to initiate apoptosis pathway. Signaling transduction protein MAP3K5, activated by oxidative stress, triggers TF E2F7 to downregulate target gene APAF1 and thereby involve in the maintenance of cell-cycle arrest upon DNA damage and promoting apoptosis [43][44][45]. Hence, in order to fight to the excessive accumulation of ROS upon decreasing the ability of DNA repair from young adulthood to middle-age, functions of apoptosis and cell-cycle arrest are raised.
In the core pathways of middle-aged skin aging only, the ligand CXCL1 (oxidative response) interacts with receptor CXCR2 and also activates signaling transduction protein CDK4. Phosphorylation of CDK4 positively regulates cell cycle entry and can stabilize and activate FOXM1 to upregulate target gene CCNB1 to promote cell cycle phase and suppress the level of ROS [48,49]. Moreover, the ligand IGF1 (oxidative response) is received by receptor IGF1R to activate PI3K/AKT signaling pathway. Transduction protein AKT1 can activate TF FOXO3 through the modification by phosphorylation. Furthermore, TF FOXO3, which is modified by phosphorylation, downregulates genes SESN3 and GADD45A. With the silence of target gene SESN3, ROS rescue pathway is declined, thus accelerating apoptosis with the increment of FOXO3-induced ROS. In addition, TF FOXO3 downregulates target gene GADD45A to promote ROS accumulation through cell-cycle arrest. The ligand TNF (proinflammatory cytokine) can inhibit collagen stability and skin homeostasis through activating TNFRSF1A and initiate the corresponding pathway. TF GATA2 was activated by phosphorylated transduction protein STAT1 to upregulate target gene MMP9. Increased gene MMP9 can inhibit collagen synthesis and enhance collagen degradation. The ligand NGF (proliferating keratinocytes) interacts with receptor NTRK1 to activate TF ETS1, which was identified to be associative with skin aging [59,60]. MIR573 is activated by TF ETS1 and then negatively regulates gene ERRFI1. Downregulated gene ERRFI1, which was modified by hypermethylation, can maintain proper epidermal homeostasis [62,63].
In the core pathways of both middle-aged and elderly skin aging, the ligand NGF/NTRK1 (signal of positive regulation of Ras signaling pathway) activates TF GATA2 so as to regulate gene BCL2 to protect keratinocytes from cell death [59]. Target gene BCL2 not only triggers antiapoptotic function through the modification by phosphorylation, but also regulates cell cycle progression to act as an antioxidant of intracellular ROS [65,66]. The ligand KITLG interacts with receptor KIT to promote melanin synthesis, DNA damage, and inhibits cell-cycle arrest by activating TF MITF through signaling transduction. TF MITF regulates melanocyte pigmentation by inducing gene TYR [67]. Due to the downregulation of gene CDH1, which is modified by phosphorylation, by TF AR, G1 phase is shortened to accumulate DNA damage and undergo apoptosis [69,70]. Next, the ligand CYR61 (oxidative stress) is received by receptor LRP1 to downregulate collagen members and promote collagen degradation [71,72]. After TF ETS1 is activated by signaling transduction, target gene COL17A1 can be downregulated to destroy the balance between collagen stability and skin homeostasis and eventually cause epidermal thinning [73].
In the core pathways of elderly skin aging only, the ligand CYR61 (oxidative stress) interacts with receptor LRP1 so as to contribute to the CCN1-induced ROS accumulation. When target gene CASP3 is downregulated by TF NDUFS4, senescence fibroblast could resist apoptosis death [74,75]. Moreover, since the ligand IL6 interacts with receptor IL6R to activate TF YAP1, which could promote cell proliferation in basal layer, it has been identified to play a physiological role in skin homeostasis [77]. Note that IL6 has been suggested as a biomarker of elderly health status [76]. After TF YAP1 downregulating target gene CDC5L through activating MIR126, functions of mitotic arrest and DNA damage were activated [78].

Two Multiple-Molecule Drugs Based on Identified Biomarkers to Mitigate Human Skin Aging
For mitigating the skin aging from young adulthood to middle age, we proposed a multiple-molecule drug including niridazole, liothyronine, decitabine, pinacidil, and allantoin. The drug targets were AIFM1, CAT, IGF1R, and LMNA as shown in Table 2. The black dot in Table 2 represents the proposed small molecules target to which identified biomarker (drug target). For instance, the niridazole has more potential to target to AIFM1 and CAT. Niridazole, an antiparasitic drug, could suppress delayed dermal hypersensitivity [79]. Studies have shown that combined therapy with liothyronine improved the treatment of hypothyroidism [80,81]. Decitabine, a DNA methyltransferase, induced changes in gene expression and cellular behavior associated with a regenerative response. Furthermore, wounds treated by decitabine were able to participate in regeneration [82]. Pinacidil is an effective antihypertensive drug for the treatment of mild to moderate essential hypertension [83]. In the meanwhile, according to the findings of one study, pinacidil may be utilized to prevent from UV-induced skin aging [84]. It is noted that allantoin, which is found in plants like chamomile, wheat sprouts, sugar beet, and comfrey, has been widely used in anti-aging serum [85,86]. Allantoin is also a well-known anti-irritating and hydrating agent as well as a peeling agent for skin [87,88]. For mitigating the skin aging from middle-aged to elderly, we proposed a multiplemolecule drug consisting of allantoin, diclofenac, mepyramine, resveratrol, and azathioprine. The drug targets were MMP9, IL6, BCL2, and CASP3 as shown in Table 3. In Table 3, the black dot shows the drug target to each specific small molecule. For example, the drug target for allantoin are MMP9 and IL6. Diclofenac is a nonsteroidal anti-inflammatory drug.
It has been used to treat actinic keratoses developing in fair-skinned individuals with a history of overexposure to ultraviolet light [89]. To mepyramine, it works by preventing the action of histamine, which is a compound produced by the body when getting venom from insect bites [90]. Moreover, one study mentioned the stimulation from histamine would upregulate matrix metalloproteinase 9 (MMP9), which is also our proposed drug target for mepyramine [91]. Resveratrol is abundant in grape skin and seeds [92]. Responding to infection, stress, injury, bacteria or fungal infections, and UV-irradiation, it a popular ingredient in skincare products [93]. In the field of dermatology, azathioprine is an effective immunosuppressant that is extremely valuable in treating pemphigoid, generalized eczematous disorders, and actinic dermatitis [94]. Taken together, most of the proposed small-molecule compounds are approved by the U.S. Food and Drug Administration (FDA). Drug repurposing for identifying new uses of old drugs with the proposed systems biology approaches might provide the alternative way to find the effective drugs for mitigating skin aging.

The Limitations and Advantages to the Proposed Systems Medicine Design Procedure for Human Skin Aging
Gene expression has been widely used to infer other molecular type measures, such as proteomics, copy number variation, and mutation. In this study, we used human skin microarray data processed with cubic spline interpolation to help us construct GWGENs by system identification method via solving constrained linear least-squares estimation problem. After that, we computed Akaike's information criterion (AIC) for each gene to prune false positives. Increasing samples through data interpolation and computing AIC for detecting real systems order, we conquered the overfitting issue. Even though we applied AIC and performed the data interpolation for increasing sample size in each skin aging stage, it is noted that the estimated real GWGENs are near-optimum solutions but not unique solutions. Furthermore, we include basal level in protein, gene, miRNA, and lncRNA systems modeling. These terms imply the unknown interaction or epigenetic modification, and mutation. If we found a basal level change, which was higher than a threshold, we inferred the corresponding node was influenced by epigenetic modification or mutation. These findings have to be verified by a literature survey. Based on the progression molecular mechanisms in each skin aging stage, we could identify essential biomarkers. For exploring the drug-target interaction to our identified biomarkers, we trained a deep neural network of drug-target interaction in advance. In the drug-target data which we used to train the prediction model, if pairs have not been mentioned as known interactions in the BindingDB, we would assign them in the group of negative samples, meaning no interaction. However, the negative samples in our study do not mean without interaction. They might just be lack of experimental evidence or record nowadays. Although the proposed system medicine design procedure exists the aforementioned limitations, it still provides another viewpoint to shed the light on the human skin aging progression based on system level. Moreover, drug repurposing strategy, giving new uses for old drugs, has been used in this study. Most of the suggested small molecules are approved by the FDA, which could shorten the time of clinical trials. Integrating systems biology approaches, deep learning framework and the design of two filters, we not only transferred biological knowledge into engineering interpretation but also applied them to drug discovery efficiently.

Overview of Systems Medicine Design Procedure of Human Skin Aging
In order to further understand skin aging molecular mechanisms from young adulthood to old age, we proposed a research flowchart as shown in Figure 2. At first, we collect several regulation and interaction databases including DIP [95], IntAct [96], BioGRID [97], BIND [98], MINT [99], HTRIdb [100], ITFP [101], Transfac [102], CircuitDB2 [103], and TargetScan [104] to construct the candidate GWGEN, which is composed of candidate protein-protein interaction network (PPIN) and candidate gene regulatory network (GRN). Moreover, the candidate GWGEN is a Boolean matrix. If two nodes have interaction, we would give one; if two nodes do not have interaction, we would give zero in it. With three-stage preprocessed microarray data, we then identify real GWGENs by system identification method and system order detection scheme. Since real GWGENs are still too complicated to investigate the skin aging progression mechanisms, we apply principal network projection (PNP) method to extract core GWGENs from real GWGENs based on the projection values. Subsequently, we denote the core signaling pathways in the style of KEGG pathways. According to the core signaling pathways, we investigate skin aging molecular mechanisms and identify essential biomarkers for young adulthood to middle age and middle age to old age, respectively. After that, we used the trained deep neural network of drug-target interaction to predict potential candidate drugs, which hold higher probability to have interactions with identified biomarkers. To narrow down the candidate drugs, we design two filters considering drug regulation ability and drug sensitivity by CMap [34] and PRISM Repurposing dataset [105]. Consequently, we propose two multiplemolecule drugs for slowing down human skin aging from young adulthood to middle age and from middle age to old age, respectively.

Data Preprocessing of Human Skin Microarray Data
We obtained human skin microarray data from GSE18876 containing the gene expression level of male skin. It included 50 ages in the range from 19 to 86 years old with 29,226 probes. One study has shown that OR52N2, SIRT6, CPT1B, TUBAL3, COL1A1 and MATN4 were significantly regulated with age. Furthermore, it also indicated that gene expressions of OR52N2, SIRT6 and CPT1B increased with age and gene expressions of TUBAL3, COL1A1 and MATN4 decreased with age [106]. Therefore, we sketched the changes of gene expression levels of these typical genes. Based on this line graph and gene expression trend in aforementioned study, we defined young-adult, middle-aged and elderly skin as 19 to 45 years old, 43 to 65 years old and 64 to 86 years old, respectively. That is, the averages of gene expressions of OR52N2, SIRT6 and CPT1B increased and the averages of gene expressions of TUBAL3, COL1A1 and MATN4 decreased from young adult stage to middle age, and then to old age in human male skin. In the estimation problem, one would easily face overfitting issue when the sample size is small and the feature size is big [107]. Hence, in this study, firstly, we increased the sample size to 500 for each skin aging stage by performing cubic spline data interpolation via splin, a MATLAB function [108][109][110]. Secondly, we utilized system order detection scheme by computing the AIC value to prune the false positives in the candidate GWGEN for finding the real GWGENs of the human skin aging systems. The more details would be discussed in the Section 4.5.

Dynamic Systems Modeling for the Candidate GWGEN
The candidate GWGWN consisting of PPIN and GRN. It is noted that GRN also includes miRNA regulation network and lncRNA regulation network. In the following contents, we would take PPIN and GRN as an example, and the rest of them could be found in Supplementary Materials. The PPIs of human-protein i in the candidate PPIN can be described as a dynamic equation shown as below: , for i = 1, . . . , I, −σ P i ≤ 0 and λ P i ≥ 0. (1) where p i (t), p j (t), and g i (t) indicate the expression levels of the ith protein, the jth protein, and the ith gene at time t, respectively; α ij denotes the interactive abilities between the ith protein with the jth protein in human skin cells; σ P i represents the degradation rate of the ith protein; λ P i indicates the translation effect from the corresponding mRNA to the ith protein; The basal level β P i signifies the regulations from other unknown regulators to the ith protein; I i denotes the number of human proteins interacting with the ith protein in the candidate GWGENs; P i (t) signifies the noise of the ith protein owing to model uncertainty or measurement noise at time t.
The k gene in the candidate GRN can be represented as a dynamic equation in the following: where g k (t), p i (t), m r (t), and o (t) indicate the expression level of the kth gene, the ith transcription factor(TF), the rth miRNA and the th lncRNA at time t, respectively; a G ki , −b G kr , and c G k represent the regulatory abilities of the ith TF, the repression ability of the rth miRNA, and the regulatory ability of the th lncRNA on the kth gene, respectively; −µ G k signifies the degradation rate of the gene expression of the kth gene; The basal level δ G k denotes the regulations from other unknown regulators to the kth gene such as phosphorylation; ω G k (t) signifies the noise of the kth gene owing to model uncertainty or measurement noise at time t; I k , R k , and L k mean the total number of TFs, miRNAs, and lncRNAs in the candidate GRN, respectively. Note that the biological regulatory mechanisms in skin cell in (2) involve TF transcription regulations by ∑ , the mRNA degradation by −µ G k g k (t), the basal level by δ G k , and the noise by ω G k (t). In this study, the effect of post-translational modification on skin aging is considered by the basal level term δ G k .

Systems Identification Approach in the Candidate GWGEN via Microarray Data
After systems modeling by Equations (1)-(4), we then perform the systems identification by solving the parameter estimation problems. The PPIN in Equation (1) can be rewritten in the following linear regression form: where ψ P i (t), represents the regression vector that can be obtained from the microarray data and θ P i denotes the unknown parameter vector to be estimated for the ith protein in PPIN. Furthermore, the Equation (3) of the ith protein can be augmented for Y i time points shown as below: . . .
which can also be simplified by where Therefore, the interaction parameters in the vector θ P i can be estimated by solving the following constrained least-squares estimation problem: where To estimate the interaction parameters in (1) by solving the parameter estimation problem in (6), we use an optimization toolbox function lsqlin in MATLAB. Simultaneously, we ensure the protein translation rate λ P i and the protein degradation rate −σ P i to always be non-negative and non-positive value, respectively; that is, λ P i ≥ 0 and −σ P i ≤ 0. Similarly, we rewrite the dynamic equation of GRN in the Equation (2) as the following linear regression form: where ψ G k (t), represents the regression vector that can be obtained from the microarray data and θ G k signifies the unknown parameter vector estimated for the kth gene in the GRN. Moreover, Equation (7) can be augmented for Y k time points in the following form: . . .

Pruning False Positives in Candidate GWGENs to Obtain Real GWGENs by System Order Detection Scheme
Due to the collected data, which we used for constructing the candidate GWGEN, come from different databases, the various experimental conditions and noises might result in getting many false-positive interactions and regulations after doing system identification. Thus, we have to apply system order detection scheme by computing AIC to detect the real system order of PPI model in Equation (1) and GRN model in Equation (2). According to Akaike's theory, the most accurate model has the smallest AIC value [111]. In other words, when the value of AIC achieves the minimum, the detected system order approaches to the real system order.
For PPI model in Equation (5), the AIC value of the ith protein can be defined in the following equation: whereθ P i denotes the estimated interactive parameters of the ith protein from the solutions of the parameter estimation problem in Equation (6), and the covariance of estimated In order to find out the real system order K * i of the ith protein in the PPI model so that AIC P i (K * i ), in Equation (11) can achieve the minimum value, we trade off the system order and the estimated residual error. By aforementioned system order detection method, PPIs with insignificant interaction abilities, which are out of K * i , could be regarded as false positives and be pruned away. For the GRN model in Equation (9), AIC value of the kth gene can be defined as the following equation: whereθ G k denotes the estimated regulatory parameters of the kth gene from the solutions of the parameter estimation problem in Equation (10), and the covariance of estimated In order to find out the real system order I * k , R * k , and O * k of the kth gene in GRN so that AIC G k (I * k , R * k , L * k ), in (12) can achieve the minimum value, we trade off the system order and to estimate residual error. In this way, to kth gene, the gene regulations with insignificant regulatory abilities, which are out of I * k , R * k , and O * k , can be treated as false-positives and be pruned away from the candidate GRN. It is noted that we apply the same system order detection scheme on the miRNA model and the lncRNA model, which could be found in the Section S1.3 of Supplementary Materials.
After performing system identification and system order detection scheme, which pruned away the insignificant interactions and regulations in the candidate GWGEN, we eventually obtained the real GWGENs for three stage of human skin aging. However, it is still quite difficult to investigate the progression mechanisms of skin aging from these real GWGENs due to their high complexity. Here, we introduce the principal network projection (PNP) method to extract the core networks from the real GWGENs as core GWGENs to solve this issue. The details are described in the following section.

Extracting Core Networks from Real GWGENs by the Principal Network Projection Method
The PNP method is a network structure projection approach based on the principal singular values so as to reduce network dimension via deleting insignificant structures. In order to use the PNP method to extract the core networks from the real GWGENs, we have to construct a network matrix H consisting all of the estimated interactions and regulations in the real GWGEN (with the ith row denoting the interactions or regulations on the ith node, i.e., protein, gene, miRNA or lncRNA of real GWGEN) in the following formation: whereα ij denotes the interactive abilities of the ith protein with the jth protein in the PPIN which could be obtained fromθ P i by solving parameter estimation problem in Equation (6) and pruning the false positives by AIC in Equation (11);â G ki ,b G kr , andĉ G kz represent transcriptional regulative abilities from the ith TFs, the rth miRNAs and the zth lncRNAs onto the kth protein-coding genes, respectively, which could be obtained fromθ G k by solving parameter estimation problem in Equation (10) and pruning the false positives by AIC in (12);â M ri ,b M rr , andĉ M rz indicate the transcriptional regulative abilities from the ith TFs, the rth miRNAs and the zth lncRNAs onto the rth miRNA's gene, respectively, which could be acquired fromθ M r by solving parameter estimation problem in Equation (S6) and pruning the false positives by AIC in Equation (S11);â L zi ,b L zr , andĉ L zz indicate the transcriptional regulative abilities from the ith TFs, the rth miRNAs and the zth lncRNAs onto the zth lncRNA's gene, respectively, which could be acquired fromθ L z by solving parameter estimation problem in Equation (S10) and pruning the false positives by AIC in Equation (S12). It is noted that if interactions or regulations do not exist in the candidate GWGEN via big data mining or already have been pruned by AIC, the corresponding components in matrix H are padded with zero.
As the H have been constructed, we thereby extract the core GWGEN from the real GWGEN by the PNP method shown as below. At first, the combined network matrix H can be a factorization of the following singular value decomposition (SVD) form as below: where U ∈ R (I+K+R+Z)×(I+R+Z) , V ∈ R (I+R+Z)×(I+R+Z) , and D = diag(d 1 , · · · , d I+R+Z ). D is composed of I + R + Z singular values of H and d 1 ≥ d 2 ≥ · · · ≥ d I+R+Z . The eigen expression fraction E h is defined in the following energy normalization: Then, we find out the minimum γ such that γ ∑ h=1 E h ≥ 0.85. That is, top γ singular vectors of matrix H contain 85% core network structure of the real GWGEN from the energy point of view. Additionally, we define the projections of H to the top γ singular vectors of V as N R (w, s) = h w,: × v T :,s for w = 1, 2, . . . , I * + R * + Z * and s = 1, 2, . . . , γ where h w:, and v T :,s denote the wth row of H and the sth column of V, respectively. Subsequently, for the top γ right-singular vectors, we define the 2-norm projection value of proteins, genes, lncRNAs, and miRNAs (i.e., the nodes) in the real GWGEN as below: [N R (w, s)] 2 1/2 for w = 1, 2, . . . , I * + R * + Z * and s = 1, 2, . . . , γ (17) If the projection value, D R (w), approaches to zero for the wth node, it means that the wth node is almost independent to the principal network structure. That is, the larger the projection value is, the greater the contribution of the corresponding node to the core network is. By doing so, we can extract the core GWGEN by collecting nodes with large projection values from the real GWGENs and denote them in the KEGG pathway style to investigate the progression mechanisms of human skin aging.

Data Preprocessing for Training Deep Neural Network of Drug-Target Interaction in Advance
The drug-target interaction dataset comes from BindingDB [35]. The descriptors of drugs and targets are transformed by PyBioMed [112]. We install this package and import PyMolecule module and PyProtein module to transform drugs and targets into their descriptors under python 2.7 environment. The PyMolecule module in PyBioMed is responsible to compute the commonly used structural and physicochemical descriptors to be drug features. The drug features include constitutional and geometrical descriptors. Furthermore, the PyProtein module in PyBioMed is responsible for calculating the widely used descriptors, including structural and physicochemical properties of proteins and peptides from amino acid sequences, to be target features. Subseqently, concatenating the drug descriptor and the target descriptor, we describe properties of a drug and its target by a feature vector shown in (18). Moreover, the total number of drug features and target features are 363 and 996, respectively.
where v drug−target indicates a feature vector of a drug-target pair; D denotes the feature vector of the drug; d i indicates the ith drug feature; T represents the feature vector of the target; t j is the jth target feature; I is the total number of drug features; J is the total number of target features. We conduct the same transformation for all the drug-target pairs to obtain their drug-target feature vectors.
Supplementary Materials: The following are available online, Figure S1: The real genome-wide genetic and epigenetic network (GWGEN) of young-stage skin, Figure S2: The real genome-wide genetic and epigenetic network (GWGEN) of middle-stage skin, Figure S3: The real genome-wide genetic and epigenetic network (GWGEN) of elder-stage skin, Figure S4: Deep neural network of drug-target interaction framework, Table S1: The pathway enrichment analysis of proteins through applying the DAVID in the core GWGEN of young-stage skin, Table S2: The pathway enrichment analysis of proteins through applying the DAVID in the core GWGEN of middle-stage skin, Table S3: The pathway enrichment analysis of proteins through applying the DAVID in the core GWGEN of elder-stage skin, Table S4: Drug targets with their corresponding small-molecule compounds, Table  S5: Drug targets with their corresponding small-molecule compounds.