By applying a system identification method and system order detection scheme to two-sided microarray data (see materials and methods, supplementary methods and the flowchart in Figure 1), we identified the early-stage and late-stage GEINs of three biological replicates using the network visualizing software Cytoscape (9) (Figure 2). The numbers of identified nodes and edges are also shown in Tables 1 and 2, respectively. Among all three replicates, the node number of host transcription factors (TFs) at the early stage is higher than that at the late stage. In addition, the identified edges in Table 2 show significant differences in host-TFs to host-genes regulation between the two stages. These results reveal that the activities of host TFs are more abundant during the early stage.

To further characterize genes in Caco-2 cells according to their functional groups, we performed an enrichment analysis via the function annotation tool DAVID (10) on the conserved target genes among all three replicates based on the biological process categories of GO database and the protein information resources of the Swiss-Prot database (Table 3). The early stage of CDI was characterized by the disturbance of cell shape and epithelial cell barrier, as well as immune activation and metal binding, which plays an important role in the scramble for metallic nutrients between the host and pathogen. At the late stage, the results including inflammatory-related functions and molecule secretion/transport suggests that a strong inflammatory response is triggered to eliminate the pathogen.

Since GEINs are very complex, it is difficult to investigate the precise host–pathogen interaction process from these networks directly. We, therefore, performed the principal network projection method (PNP) method to extract the core nodes with high projection values, which compose the corresponding HPNs from early-stage and late-stage GEINs of Caco-2 cells during CDI.

Considering the crucial role of CD2787 (cwp84) in pathogenesis at the early stage of CDI, we recommend CD2787 as a potential drug target in the prevention of CDI. Moreover, other cell-wall proteins participating in the infectious process, such as CD0440, CD0237, CD0266, and CD1987 are also potential drug targets since the inhibition of these cell surface proteins may reduce the occurrence of cell adhesion. In Figures 3 and 4, four pathogen genes (CD0130, CD3256, CD1275, and CD2781) were previously predicted as essential survival/growth genes by applying flux balance analysis (FBA) and synthetic accessibility (SA) to a curated C. difficile metabolic network, as well as transposon-directed insertion site sequencing (TraDIS) to the C. difficile transposon mutant library (19, 20). However, CD0130 and CD3256 have important human-homologs, MAT2A and VARS, respectively (19). The inhibition of these proteins may result in unpredictable dysfunction of host cells. Therefore, CD1275 and CD2781 are more conservative choices for further drug design. Similarly, proteins involved in the defense mechanisms employed by C. difficile against reactive oxygen species (ROS), including CD0141 and CD2115 in the early stage of CDI, and proteins responsible for redox-state homeostasis (CD1690 and CD1631) in the late stage of CDI, have human-homologs (19). Therefore, the three remaining anti-ROS enzymes (CD2356, CD0171, and CD0179) are recommended as potential targets to repress the defense mechanisms of C. difficile. The heavy economic burden of CDI results from its high recurrence rate, suggesting that the inhibition of spore formation of C. difficile is a feasible therapeutic way to reduce the recurrence rate of CDI. In our results, the key sporulation pathway members included CD1214, CD2629, and CD2643. These proteins could become targets of further drug design.

In addition to the inhibition of pathogen mechanisms, another strategy of drug design includes the promotion of host defense mechanisms. It has been reported that low doses (0.2 ng/ml) of CD0660 is sufficient to induce significant cell rounding in vitro (21), suggesting that amplification of host defense appear to be required. In our results, the major pathogenic effects induced by C. difficile pathogenic factors include tight junction/cytoskeleton breakdown and the accumulation of endoplasmic reticulum (ER) stress, which result from the inactivation of Rho GTPases (RHOA, CDC42, and RAC1) and the dysfunction of chaperone proteins (HSP90B1, HSPA5, and HSP90B2P), respectively. Therefore, we aim to increase the expression of these genes to maintain the activities of these dysfunction components. The severe inflammation and apoptosis process, induced by NFKB1, REL, and IL-8, are responsible for cell death at the late stage of CDI. Though the activities of immune response and inflammation are important counter mechanisms against bacteria, the overexpression of these related genes cause tissue damage and apoptosis of host cells. Therefore, we suggest repressing the expression of NFKB1, REL, and IL-8 to relieve these activities.

Our results show that CD2356, CD0171, and CD0179 participate in the defense mechanisms of C. difficile against oxidative stress. The co-operation among these proteins provides a well-designed protection against human-produced ROS. The inhibition of these antioxidative proteins facilitates the host eliminating ability against pathogens, and the scattered ROS can induce the rapid necrosis of pathogen cells. Therefore, we recommend that CD2356, CD0171, CD1064, and CD0179 are potential drug targets for further drug design (see supplementary methods and Table S2 in Supplementary Material and Figure S5 in Supplementary Material).

During the pathogenesis of the bacterium, essential genes are required for the survival and growth of the pathogen. The absence of essential genes results in a decreased growth rate or death of the organism. The pathogen core networks extracted from GEINs display not only essential proteins required for the survival of the pathogen but also important enzymes that contribute to the pathogenesis and defense mechanism of the pathogen against various stress conditions. Here, we compare HPNs with the existing predicted essential proteins of C. difficile to show the difference and advantage of HPNs.

Recently, by applying FBA and SA to a curated C. difficile metabolic network, Larocque et al. predicted 76 essential C. difficile genes (19). Seven C. difficile proteins (CD2664, CD2335, CD3550, CD0198, CD1225, CD0130, and CD0123) in the HPN of the early stage of CDI (Figure S3 in Supplementary Material) and three C. difficile proteins (CD2588, CD1816, and CD0130) in the HPN of the late stage of CDI (Figure S4 in Supplementary Material), encoded by C. difficile genes appear to be required for the survival of the pathogen, have been identified in this study. Moreover, applying transposon-directed insertion site sequencing (TraDIS) to the C. difficile transposon mutant library led to the identification of 404 genes with no transposon insertion in the library, which can be considered as essential genes for C. difficile growth (20). Thirteen C. difficile proteins (CD2664, CD2335, CD0067, CD3550, CD3540, CD0198, CD1255, CD2714, CD3256, CD1316, CD0095, CD2739, and CD0123) in the HPN of the early stage of CDI (Figure S3 in Supplementary Material) and 15 C. difficile proteins (CD3170, CD2588, CD2744, CD2771, CD2462, CD3304, CD2781, CD3540, CD1145, CD2461, CD2793, CD0059, CD1275, CD0052, and CD1767) in the HPN of the late stage have been identified in this study as the products of likely essential genes for the growth of C. difficile. Overall, 15 pathogen genes in the HPN of the early stage of CDI, and 17 pathogen genes in the HPN of the late stage, were previously identified as likely essential genes of C. difficile (Table S1 in Supplementary Material).

Besides these previously identified essential proteins, HPNs also provide numerous crucial C. difficile enzymes that participate in the offensive and defensive mechanism utilized by the pathogen, such as well-known toxins (CD0660 and CD0663), toxin-regulators (CD0659, CD0661, and CD0664), cell-wall proteins (CD2787, CD1987, CD0237, and CD0440) that provide cell-adhesion abilities, defensive proteins (CD0141, CD0171, CD2115, CD1631, and CD1690) against ROS, and sporulation-related proteins (CD1214, CD2643, CD2629, and CD1511). The presence of these proteins in HPNs reveals that there are various cross-talk activities between the host and pathogen in the offensive and defensive mechanism.

The interactions between the small GTPase (CDC42) and downstream effectors (GRB2 and KBTBD7) suggest that the dynamic of cytoskeleton homeostasis can be affected by C. difficile toxins and host epigenetic activities to influence the change of cell shape displayed in Table 3. Similarly, at the late stage of CDI (Figure S4 in Supplementary Material), the subunits (NFKB1 and REL) of the NF-κB complex are regulated by acetyltransferase proteins (GCNT2 and B3GNT6) and the methyltransferase protein (PRDM14), resulting in the assembly of the NF-κB complex and the induction of apoptosis in host cells. The change in the methylation level of HSPA5 has been proposed in Hesson et al. (22) via NOME-Seq analysis of HSPA5 in intestinal disease, suggesting that altered nucleosome positioning induces differences in accessibility. In addition, in the late-stage HPN shown in Figure S4 in Supplementary Material, the significant difference in the basal level of the NF-κB subunit NFKB1 gene reflects the change of methylation level, which is consistent with a previous study showing that NFκB1 is hypomethylated in intestinal inflammation (23). These identified HPNs reveal the ability of acetylation and DNA methylation likely to alter cellular functions in order to adapt to bacterium infection.

Once the bacterium is attached to the surface of Caco-2 cells in the early stage of infection (Figure 3), the cell surface-associated cysteine protease CD2787 (cwp84) of bacterium interacts with FN1 on the host membrane. This negative interaction (CD2787, FN1) and the MIB2-induced ubiquitination of FN1 suggest that CD2787 could provide a degrading activity on FN1, which can be supported by the biological experiment (11). Since FN1 is involved in the maintenance of cell shape, it can bind cell surface and various compounds, including actin. We suggested that the degradation of FN1 results in the morphological change of Caco-2 cell or the disturbance of the actin cytoskeleton. The results also suggest that three pathogen proteins (CD0660, CD0663, and CD0478) secreted by C. difficile interact with four host receptors [FPR1, SCARA3, HSP90B1, and Arrestin Beta 2 (ARRB2)]. FPR1 has been reported as the receptor of CD0660 (TcdB) (24). Owing to its interaction with toxins (CD0660 and CD0663) and ubiquitination modified by ARIH2, we suggested that C. difficile toxins enter the host-cell cytoplasm via receptor-mediated endocytosis. For C. difficile toxins, another identified human cell surface receptor HSP90B1, which encodes a member of HSP 90 kDa family, plays a role in protein folding. It has been reported that the ER chaperones HSP90B1 and HSPA5 could be significantly downregulated by HDAC inhibition at the protein level (25). HSP90B1 has been identified as a cell surface receptor for C. difficile proteins, but HSP90B1 did not respond to CD0660 (13). The observation can be supported by the result in the GEIN that HDAC inhibition impairs the binding affinity of HSP90B1 for CD0660.

The interaction between CD0663 and CD46 (Figure 4) indicates that pathogen CD0663 could be the extracellular stimuli of CD46, and potentially responsible for CDI to recruit phagocytic cells such as neutrophils to clear toxins and whole microbes in the early stage. The result in Figure 4 suggests acetylation and interaction with toxins (CD0660 and CD0663) of HSP90B2P, which participate in impairing its chaperone-ability. Interestingly, we noticed that there are more host proteins influenced by CD0663 than CD0660 in the late stage of infection. By contrast, CD0660 binds more host proteins than CD0663 in the early stage (Figure 3). These results provide an explanation for the perennial argument about the cytotoxic responsibility of CD0660 (TcdB) and CD0663 (TcdA), suggesting that both toxins are likely essential for pathogenesis. The results suggested that this role change may result from CD2973-induced acetylation of CD0663 in the late stage (Figure 4). The rise of CD0663 activity could be the turning point for the progression of CDI. We observed that another host protein, FN1, could also change its role from the early stage to the late stage of CDI. FN1 is repressed by ubiquitination and CD2787-induced degradation in the early stage of CDI (Figure 3), but the CD46-triggered activation and NAT8L-induced acetylation of FN1 increased the corresponding expression levels (p-value < 8 × 10⁻³) despite the influence of CD2787, thus allowing FN1 to trigger downstream activities in the late stage of CDI in Figure 4. Therefore, we suggested that CD0660 functions prior to CD0663 and triggers rapid responses in the early stage of infection and CD0663 works actively in the late stage.

To further understand this, we separate the rearranged core pathways in Figures 3 and 4 based on different cellular responses and their corresponding pathways for further discussion. For the early stage of CDI, the genetic-and-epigenetic scheme of the core pathways (Figure 3) can be separated into three parts: the offensive mechanism of the pathogen and the corresponding pathogenesis of host cells (Figure 5A); the remedial actions employed by host cells in response to pathogen-induced injuries (Figure 5B); and the counter mechanisms of Caco-2 cells and the defense mechanisms of C. difficile (Figure 5C). During the late stage of CDI, the core pathways (shown in Figure 4) can be separated into two parts: the strong cellular responses triggered by epigenetic acetylation in host cells for depleting the tenacious pathogen at the late stage of CDI (Figure 6A) and the severe inflammation and apoptosis processes of Caco-2 cells as well as the endospore formation of C. difficile (Figure 6B).

Using a systems biology approach to construct genome-wide GEINs via big data mining and two-sided genome-wide data identification, we investigated the pathogenic model of CDI (Figure 7). The result suggests that the acetylation of HSP90B1 and the deacetylation of HSPA5, and the interaction of these HSPs with pathogen toxins (CD0660 and CD0663), impair the chaperone-activity of host cell, resulting in the formation and accumulation of misfolded proteins. In addition, C. difficile toxins (CD0660 and CD0663) display their cytotoxicity by inactivating Rho GTPases (RHOA and CDC42) via glucosylation, causing the disturbance of actin cytoskeleton homeostasis and tight junction breakdown (3). Another central epigenetic activity is ubiquitination. The ubiquitination-mediated endocytosis of FPR1 allows C. difficile toxins enter host cells and then CD0660 and CD0663 hijack a RAC1-mediated pathway to activate SYK-dependent ROS production (24, 28). ROS production is a primary response of the immune system, but the toxin-activated overproduction of ROS and free radicals could damage host lipids, proteins, and DNA. To repair the damaged DNA, SCARA3 and YY1 are activated to initiate the DNA damage response. The dysfunction of HSPs together with YY1 trigger the autophagy process, removing misfolded proteins, and ROS-induced damage (41). As a countermeasure, host miRNA can pass through the pathogen cell wall and then repress pathogen genes, including the antioxidative proteins CD3256 and CD0130 to promote the eliminating ability of host cells.

When the infection proceeds to the late stage (Figure 7B), several events together result in the turning point of CDI. The result suggests that the acetylation of CD0663 promotes the activity of this toxin and, thus, increases the probability of CD0663 to bind to CD46. CD46 can activate the complement system and enhance ROS production via NOX5. The complement system guides neutrophils, which can be recruited by interleukin-8 to remove toxins and pathogens (42). In addition, FN1 can be activated by CD46 and modified by acetylation, and these interactions could enhance the activity of FN1 and its downstream inflammatory responses. Furthermore, the newly produced CD1466 displays its immunogenicity by interacting with SNW1 and EGFR, thus inducing cytokine production and inflammation. Taken together, the presence of CD1466 and CD46, and the acetylation of CD0663 and FN1, result in enhanced oxidative stress and a severe inflammatory response. Meanwhile, the acetylation of HSPs (HSP90B1 and HSP90B2P) and their interactions with toxins (CD0663 and CD0660) impair the chaperone-activity, aggravating the accumulated ER stress (43). These activities in turn increase the cellular stress of the host cell. For example, neutrophil infiltration and the increased NADPH oxidase could cause tissue damage of the host cell, and the presence of C22orf28 reflects high ER stress (33). On the other hand, these processes are assuredly life threatening to C. difficile. To counteract these threatening stresses, C. difficile utilize numerous redox-related proteins, including superoxide dismutase CD1631, extracellular glutamate dehydrogenase CD0179, and thioredoxin CD1690 in the defense against oxidative stress. Finally, the accumulated oxidative and ER stress trigger the apoptosis process via the NF-κB complex. For C. difficile, the high levels of ROS and scattered neutrophil risk the survival of pathogens. Therefore, C. difficile transform to endospores by activating the sporulation pathway, and then lie dormant for the next infection.

As discussed above, CDI is characterized by cytoskeleton dysfunction, severe inflammation, and subsequent apoptosis. The actin cytoskeleton breakdown is mainly induced by CD0663 and CD0660, but the correlation between cytoskeleton dysfunction and apoptosis remains unclear. Here, we suggest that the main cause of apoptosis during CDI is the accumulation of cellular stress, including oxidative and ER stress. The emergence of oxidative and ER stress has been reported in studies using CD0660-infection models (21, 44) but generates few attention and further researches. We suggested that Caco-2 cells activate the DNA damage response and autophagy to counteract these stressors in the early stage of CDI. Unfortunately, the accumulated stress and tissue damage caused by severe inflammation induce host cells to undergo apoptosis in the late stage of infection.

There has been a long-lasted argument about whether CD0660 or CD0663 is responsible for the cytotoxicity in host cells. Some early studies reported that only CD0660 is essential for the virulence of C. difficile and a number of patients with CDI are caused by CD0663-negative and CD0660-positive strains (45). However, a later study using toxin knockdown technology on C. difficile 630 indicated that both CD0660 and CD0663 are responsible for disease (46). Our identified results also support this conclusion. During the early stage of CDI, CD0660 functions prior to CD0663 and triggers numerous host responses, such as ROS production and chaperone dysfunction. In the late stage, acetylation enhances the activity of CD0663, which takes the place of CD0660, inducing a severe inflammatory response and aggravating ER stress. The CD0663-activating ability of CD0660 identified in our toxin-regulation pathways is also consistent with this temporal relationship. Without CD0660, CD0663 may be not sufficient to initiate pathogenesis in host cells in the early stage of CDI and the low expression level of CD0660 cannot trigger the subsequent apoptosis in the absence of CD0663 during the late stage of infection.

By contrast, unlike pathogenic effects on the host, the intraspecies interactions and cellular mechanisms inside C. difficile are largely unknown. We found that C. difficile could generate redox-balancing proteins against oxidative stress and utilize toxin production and bacterial reproduction as offensive mechanisms at the early stage of CDI. During the late stage of infection, the result suggests that C. difficile utilized anti-ROS proteins, including CD0179, CD1631, and CD1690, as well as the DNA damage response to counteract the oxidative stress presented by host cells. The decreased activities of toxin production and bacterium-cell growth, and endospore formation also reveal that C. difficile actively transform to endospores to leave the infection site. Finally, the molecular mechanisms of progression from the early stage to the late stage of CDI are investigated, and the potential drug targets are also proposed for further drug design.

The crucial events of progression in CDI, such as the acetylation of CD0663 and FN1, the chaperone dysfunction of HSPs, and the pathogen silencing induced by host miRNA, are mainly caused by epigenetic regulations identified in our systems biology method. These results suggest that epigenetic regulation plays an important role in the progression of infection since these cellular activities could change the cellular functions of host cells in a more rapid and efficient manner than the adjustment of gene regulation.

This study has three main limitations. (1) The reported results are modeling outcomes of two-sided genome-wide expression data and are speculations and not proven. (2) Several new potential interactions, especially interspecies PPIs, based on sequence-homology between C. difficile and E. coli, and the interspecies PPIs between C. difficile and homo sapiens, need for further experimental evidence. (3) The findings relate to only on strain of C. difficile and only once specific cell line, Caco-2 cells. The true complexity of an entire gut significantly complicates interactions with other cells/bacteria/host. This is great and valuable work but certainly has limitations. To the best of our knowledge, there are few studies focusing on the epigenetic modulation on pathogenic and offensive mechanisms in the host cell infected by C. difficile. In addition, no existing whole-epigenomic data of the host cells can be considered as the basis of such studies. We show that epigenetic regulation will play a more important role in host/pathogen cross-talk mechanisms in CDI, which provides a novel direction for deeper studies on molecular mechanisms for drug target predictions and multiple drug discovery. With the identification of GEINs and HPNs, our knowledge of the bioinformatics of C. difficile and the core proteins is bound to increase, driving further experimental hypotheses and investigation directions for cross-talk mechanisms, including the protein–protein interaction network (PPIN) and gene regulation network (GRN) between host and C. difficile during infection. These further studies will complementarily complete our genetic-and-epigenetic network, providing a novel basis for developing the whole-genome cellular network of CDI.

A flowchart showing HPNs of the host and pathogen at the early and late stages of infection by big data mining, model construction, and network identification for investigating the cross-talk molecular mechanisms and inferring potential drug targets is shown in Figure 1. These processes can be divided into four steps: (1) big data mining and data preprocessing of host/pathogen gene/miRNA expression data (see supplementary methods); (2) construction of candidate GEIN, which consists of candidate host/pathogen intraspecies PPINs, candidate interspecies PPINs between host and pathogen, candidate host/pathogen gene/miRNA regulation networks, candidate miRNA regulation networks of host-miRNAs on host/pathogen-genes, and candidate lncRNAs regulation networks of host-lncRNAs on host-genes (see supplementary methods and Figure S6 in Supplementary Material); (3) identifying real GEINs of each stage from the candidate GEIN via system identification method and system order detection scheme, using the genome-wide microarray data of Caco-2 cells and C. difficile during infection; and (4) extraction of HPNs from the real GEINs using the PNP (see supplementary methods). We then investigated the crucial molecular mechanisms that contribute to the progression of CDI and inferred potential drug targets for further drug design.

The two-sided microarray raw data proposed by Janvilisri et al. (6) has two parts. The first one contains the mRNA/miRNA expression profiles of three biological replicates of the Caco-2 cell line at 0, 30, 60, and 120 min postinfection with C. difficile 630. Each biological replicate contains two technical replicates (GEO accession number GSE18407). The Caco-2 cell line was cultured in Dulbecco’s modified Eagle’s medium at 37°C prior to infection. The second part contains the mRNA expression profiles of three biological replicates of C. difficile 630 in Caco-2 cells at 0, 30, 60, and 120 min postinfection (GEO accession number GSE18407; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE18407). The platforms used in the host and pathogen were Phalanx Human OneArray and C. difficile 630/QCD32g58 array, respectively, which include 39,200 and 13,824 probes, respectively. The microarray data were validated using qRT-PCR. We took the average of the microarray data of the two technical replicates for further network identification. We applied one-way analysis of variance to calculate p-value of a gene using microarray data between early and late infection stages.

Since the candidate GEIN (see supplementary methods) was constructed by big data from numerous databases, experimental datasets, and literature, they contain some inevitable false-positive information. To address this, we built the dynamic model to characterize the molecular mechanisms of GEINs and to prune the false-positives in candidate GEIN, thus producing the real GEINs for Caco-2 cells and C. difficile during the infection process. For the PPIN of host proteins in candidate GEIN, the dynamic interaction model of the ith host protein can be described by the following dynamic equation: (1)piHt+1=piHt+∑f=1FiaifHpiHtpfHt+∑q=1QiciqHpiHtpqPt+αiHgiHt−γiHpiHt+κiH+ϖiHt, for i=1,2,…,I,αiH≥0 and−γiH≤0 where piH(t), pfH(t), giH(t) and pqP(t) represent the expression levels of the ith host protein, the f th host protein, the ith host gene and the qth pathogen protein at time t, respectively; aifH and ciqH denote the interactive ability between the ith and f th host protein and between the ith host protein and qth pathogen protein, respectively; F_i and Q_i signify the number of host proteins and pathogen proteins that interact with the ith host protein; αiH, −γiH, and κiH indicate the translation rate from the corresponding mRNA, the degradation rate, and the basal level of the ith host protein, respectively. In general, the basal level κiH in Eq. 1 represents an unknown activity affecting the expression of the ith host protein other than those mentioned above, such as epigenetic acetylation and ubiquitination. ϖiH(t) denotes the stochastic noise of the ith host protein at time t. The biological meaning of Eq. 1 is that the expression level of the ith host protein can be affected by various molecular mechanisms including the host intraspecies PPIs ∑f=1FiaifHpiH(t)pfH(t), interspecies PPIs ∑q=1QiciqHpiH(t)pqP(t), protein translation αiHgiH(t), protein degradation −γiHpiH(t), basal level κiH, and the corruption of stochastic noise ϖiH(t) In addition, the translation rate should be constrained to be non-negative and the protein degradation rate should be constrained to be non-positive in real PPIs.

For the GRN of host genes in the candidate GEIN, the dynamic model of the jth host gene can be described as follows: (2)gjHt+1=gjHt+∑i=1IjbjiHpiHt+∑n=1NjejnHlnHt+∑i′=1Ij′∑i″=1Ij″xjIj″i′−1+i″Hpi′Htpi″Ht−∑k=1KjdjkHgjHtmkHt−λjHgjHt+δjH+εjHt, for j=1,2,…,J,−djkH≤0 and −λjH≤0 where gjH(t), piH(t), mkH(t), and lnH(t) indicate the expression levels of the jth host gene, the ith host TF, the kth host miRNA, and the nth host lncRNA at time t, respectively; bjiH, −djkH and ejnH represent the regulation ability of the ith host TF, the kth host miRNA, and the nth host lncRNA on the jth host gene, respectively; I_j, K_j, and N_j denote the number of host TFs, host miRNAs and host lncRNAs, respectively, which regulate the expression level of the jth host gene; pi′H(t)pi″H(t) imply the ith host complex where pi′H(t) and pi″H(t) represent the subunit 1 and 2 of the ith host complex, respectively; Ij′ and Ij″ are the same, to denote the number of host complex, which regulate the jth host gene in the candidate GRN; xjIj″i′−1+i″H signifies the regulation ability of the ith host complex on the jth host gene; −λjH and δjH indicate the degradation rate and the basal level of the jth host gene, respectively. In general, the basal level δjH in Eq. 2 denotes an unknown regulation other than those mentioned above such as DNA methylation. εjH(t) represents the stochastic noise due to modeling residue at time t. Notably, in the case of the regulation ability xjIj″i′−1+i″H of the ith host complex on the jth host gene, the index Ij″i′−1+i″ guarantees the proper coordinate of the regulation ability xjIj″i′−1+i″H of the ith host complex pi′H(t)pi″H(t) in the host GRN system matrix of the jth host gene, i.e., the regulation abilities of host complexes on the jth host gene can be aligned to an one row matrix of xj1H,…,xjIj″H,xjIj″+1H,…,xj2Ij″H,xj2Ij″+1H,…,xjIj″i′−1+i″H,…,xjIj″Ij′H.

Therefore, the biological meaning of Eq. 2 is that the expression level of the jth host gene can be regulated by numerous molecular mechanisms, including the host TF regulations ∑i=1IjbjiHpiH(t) host lncRNA regulations ∑n=1NjejnHlnH(t), host complex regulations ∑i′=1Ij′∑i″=1Ij″xjIj″i′−1+i″Hpi′H(t)pi″H(t), host miRNA repressions −∑k=1KjdjkHgjH(t)mkH(t), mRNA degradation effect −λjHgjH(t), basal level δjH, and the corruption of stochastic noise εjH(t). Similar to protein model, the miRNA repression ability and gene degradation rate should be constrained to be non-positive. Since DNA methylation can directly influence the binding affinities of RNA polymerase to target genes (47), we assumed that the regulation by methyltransferase could cause the significant change of the basal level δjH and the change of δjH between the early stage and late stage of CDI in the dynamic model (2) implies the occurrence of methylation at the jth host gene in the infection process.

In Eq. 2, the expression of the kth host miRNA mkH(t) and the nth host lncRNA lnH(t) at time t can also be regulated by other regulators. Therefore, the dynamic regulatory equation of the kth host miRNA was modeled as follows: (3)mkHt+1=mkHt+∑i=1IkykiHpiHt−μkHmkHt+φkH+ςkHt,for k=1,2,…,K and −μkH≤0 where mkH(t) and piH(t) denote the expression levels of the kth host miRNA and the ith host TF at time t, respectively; ykiH represents the regulatory ability of the ith host TF on the kth host miRNA; I_k signifies the number of host TFs that regulate the expression level of the kth host miRNA; −μkH and φkH indicate the miRNA degradation rate and the basal level of the kth host miRNA, respectively; and ςkH(t) implies the stochastic noise at time t. The dynamic model of host miRNAs in Eq. 3 characterizes molecular regulatory mechanisms, including the transcription regulations ∑i=1IkykiHpiH(t), miRNA degradation effect −μkHmkH(t), basal level φkH, and the corruption of stochastic noise ςkH(t). In addition, the degradation rate should be constrained to be non-positive.

Similarly, the dynamic model of the nth host lncRNA in candidate GEIN can be described by the dynamic equation as follows: (4)lnHt+1=lnHt+∑i=1InzniHpiHt−χnHlnHt+ρnH+ϑnHt,for n=1,2,…,N and−χnH≤0 where lnH(t) and piH(t) represent the expression levels of the nth host lncRNA and the ith host TF at time t, respectively; zniH denotes the regulation ability of the ith host TF on the nth host lncRNA; I_n signifies the number of host TFs regulating the expression level of the nth host lncRNA; −χnH and ρnH indicate the degradation rate and the basal level of the nth host lncRNA, respectively; and ϑnH(t) is the stochastic noise due to the modeling residue. The dynamic model of host lncRNAs in Eq. 4 characterizes molecular regulatory mechanisms, including the transcription regulations ∑i=1InzniHpiH(t), degradation effect −χnHlnH(t), basal level ρnH, and the corruption of stochastic noise ϑnH(t). Furthermore, the constraint of this model is that the lncRNA degradation rate should be non-positive.

For the PPIN of pathogen proteins in candidate GEIN, the dynamic model of the qth pathogen protein can be described by the following equation: (5)pqPt+1=pqPt+∑o=1OqaqoPpqPtpoPt+∑i=1IqcqiPpqPtpiHt+αqPgqPt−γqPpqPt+κqP+ϖqPt,for q=1,2,…,Q,αqP≥0 and −γqP≤0 where pqP(t), poP(t), gqP(t), and piH(t) represent the expression level of the qth pathogen protein, the oth pathogen protein, the qth pathogen gene, and the ith host protein at time t, respectively; aqoP and cqiP denote the interactive ability between the qth pathogen protein and oth pathogen protein, and between the qth pathogen protein and ith host protein, respectively; O_q and I_q signify the number of pathogen proteins and host proteins that interact with the qth pathogen protein, respectively; αqP, −γqP and κqP indicate the translation rate, the degradation effect and the basal level of the qth pathogen protein, respectively; and ϖqP(t) denotes the stochastic noise of the qth pathogen protein at time t. The biological meaning of the Eq. 5 is that the expression level of the qth pathogen protein can be affected by various molecular interactive mechanisms, including the intraspecies PPIs ∑o=1OqaqoPpqP(t)poP(t), the interspecies PPIs ∑i=1IqcqiPpqP(t)piH(t), protein translation αqPgqP(t), protein degradation −γqPpqP(t), basal level κqP, and the corruption of stochastic noise ϖqP(t). Similar to the host protein dynamic model, the translation rate should be constrained to be non-negative and the protein degradation rate should be constrained to be non-positive.

For the GRN of pathogen genes in candidate GEIN, the dynamic model of the hth pathogen gene can be described as follows: (6)ghPt+1=ghPt+∑q=1QhbhqPpqPt−∑k=1KhdhkPghPtmkHt−λhPghPt+δhP+εhPt,for h=1,2,…,H, −dhkP≤0 and −λhP≤0 where ghP(t), pqP(t), and mkH(t) indicate the expression levels of the hth pathogen gene, the qth pathogen TF, and the kth host miRNA at time t, respectively; bhqP and −dhkP represent the regulatory ability of the qth pathogen TF and the kth host miRNA on the hth pathogen gene, respectively; Q_h and K_h denote the number of pathogen TFs and host miRNAs that regulate the qth pathogen gene; −λhP and δhP indicate the degradation rate and the basal level of the qth pathogen gene, respectively; and εhP(t) represent the stochastic noise due to the modeling residue at time t. The biological meaning of Eq. 6 is that the expression level of the hth pathogen gene can be regulated by various molecular mechanisms, including the pathogen TF regulations ∑q=1QhbhqPpqP(t), host miRNA repressions −∑k=1KhdhkPghP(t)mkH(t), mRNA degradation −λhPghP(t), basal level δhP, and the corruption of stochastic noise εhP(t). In addition, the host-miRNA repression rate and degradation rate should be constrained to be non-positive. By comparing to the well-proposed models in the cross-talk genome-wide GEINs (48), this study considered additional molecular mechanisms in the GEIN, including transcription regulations of host lncRNA and host protein complex in Eq. 2 and the dynamic models of host miRNA in Eq. 3 and host lncRNA in Eq. 4, to characterize molecular interaction during cell infection in more detail.

In order to identify the precise parameters, we applied a system identification method to the dynamic genetic-and-epigenetic Eqs 1–6 in the candidate GEIN. We rewrote the host PPIN dynamic Eq. 1 as the linear regression form below, (7)piHt+1=piHtp1Ht⋯piHtpFiHtpiHtp1Pt⋯piHtpQiPtgiHtpiHt1ai1H⋮aiFiHci1H⋮ciQiHαiH1−γiHκiH+ϖiHt≜ϕiHPtθiHP+ϖiHt,for i=1,2,…,I where ϕiHP(t) represents the regression vector that can be obtained from the microarray expression data and θiHP is the unknown parameter vector to be estimated for the ith host protein in host PPIN.

Equation 7 of the ith host protein can be augmented for T_i data points as the following form: (8)piH(t2)piH(t3)⋮ piHtTi+1=ϕiHPt1ϕiHPt2⋮ϕiHPtTiθiHP+ϖiHt1ϖiHt2⋮ϖiHtTi, for i=1,2,…,I which could be simply represented as follows: (9)PiH=ΦiHPθiHP+ΓiHP, for i=1,2,…I where PiH=piHt2piHt3⋮ piHtTi+1, ΦiHP=ϕiHPt1ϕiHPt2⋮ϕiHPtTi, ΓiHP=ϖiHt1ϖiHt2⋮ϖiHtTi.

Therefore, the parameters in the vector θiHP can be estimated by applying the following constrained least-squares estimation problem, (10)minθiHP ΦiHPθiHP−PiH22subject to 0⋯00⋯0−1000⋯00⋯0010 θiHP≤01

The parameters in the host PPIN dynamic Eq. 1 can be estimated by solving the constrained least-squares problem (10) via the help of lsqlin function in MATLAB optimization toolbox, and simultaneously the host protein translation rate αiH is guaranteed to be non-negative and the host protein degradation −γiH is guaranteed to be non-positive, i.e., αiH≥0 and −γiH≤0. Similarly, the constrained least-squares estimation problems of Eqs 2–6 are shown in supplementary methods.

As mentioned above, to avoid the overfitting problem in the parameter identification process, we have applied cubic spline to interpolate extra data points (five times number of the parameters in the parameter vector to be estimated, i.e., θIHP in host PPIN, θjHG in host GRN, θkHM in host miRNA dynamic model, θnHL in host lncRNA dynamic model, θqPP in pathogen PPIN, and θhPG in pathogen GRN). Therefore, with the microarray expression data, we could solve the constrained least-squares estimation problems in Eq. 10 and Eqs S3, S6, S9, S12, and S15 in Supplementary Material and identify the precise parameters in GEINs gene by gene (or protein by protein) via the lsqlin function of MATLAB optimization toolbox. Since the measurement technology of genome-wide protein expression of Caco-2 cells and C. difficile has not yet been realized, and about 73% variance of protein abundance can be explained by the corresponding mRNA abundance (49), the microarray data of gene expressions can replace protein expressions, providing sufficient information for solving above constrained least-squares parameter estimation problems.