Introduction
Porcine Reproductive and Respiratory Syndrome (PRRS) is one of the most economically devastating infectious diseases currently facing the global swine industry, causing direct and indirect economic losses amounting to billions of dollars annually(Holtkamp et al., 2013). The main clinical features of the disease include: affected pigs exhibiting persistent high fever, progressive respiratory distress or even respiratory failure, significant deterioration in the semen quality of boars, sows experiencing abortions, and reproductive disorders such as stillbirths and mummified fetuses(Mil-Homens et al., 2025). Pathological studies have confirmed that the Porcine Reproductive and Respiratory Syndrome Virus (PRRSV) is the sole primary pathogen responsible for the disease(Zimmerman et al., 2019).
According to the newest classification scheme of the International Committee on Taxonomy of Viruses (ICTV), PRRSV can be divided into two genotypes, possessing great genetic diversity, the European type (PRRSV-1) and American type (PRRSV-2). Phylogenetic reconstruction of ORF5 gene sequences or entire-genome levels identified the prevalent PRRSV-1 isolates in China currently distributed in at least seven different phylogenetic subgroups with evident genetic backgrounds(Yuan et al., 2024), but PRRSV-2 is even more diverse in its genome, which in turn is phylogenetically separated into as many as nine lineages that constitute different evolutionally independent branches with numerous subclones possessing specific molecular features(Zimmerman et al., 2019).
Recently published epidemiological surveillance data indicate two serious facts, namely, 1) the pathogenicity of PRRSV-1 strains in Chinese pig herds keeps expanding continuously, including more serious clinical manifestations and prolonged post-infection pathogenesis; and 2) viral recombination (mainly as recombination of genomic fragments among different genotypes or lineages) have been repeatedly documented in laboratory reports(Pedro Mil-Homens et al., 2024). Whether these new molecular epidemiological features will lead to transformation in intensity and patterns of PRRS outbreaks especially, whether it will become critical risk factors triggering next large-scale PRRS epidemic is a major scientific problem of global interest in swine disease field(Merca et al., 2022).
The prevention and control measures for PRRS mainly involve the use of inactivated or attenuated vaccines. Literature reports indicated that there was only 30% - 50% amino acid homology among different PRRSV, resulting in weak cross-protection of vaccines(Zhang et al., 2022). The unique pathogenic mechanism of PRRS induces immunosuppression, making it difficult to effectively control the PRRSV and its variants with single treatment methods such as vaccines or antibiotics alone(Wang and Feng, 2024). Clinical practice showed that adding Traditional Chinese medicine (TCM) compound to feed could effectively prevent the outbreak of PRRS(Bello-Onaghise et al., 2020).
Guizhi Powder (GP), recorded in the Shang Han Lun, is composed of Cinnamon Twig, Peony Root, Licorice, Ginger, and Jujube . Modern pharmacological studies have shown that it has anti-inflammatory, analgesic, antiviral, antitussive, expectorant, and antiasthmatic effects; it also regulates the intestinal and cardiovascular systems; through mediating the action of endogenous pyrogens, it has a bidirectional regulatory effect on the body temperature center(Liu et al., 2024). GP could significantly inhibit lung damage caused by influenza virus pneumonia in mice and increase capillary permeability(Armour et al., 2021;Ling et al., 2023). In the compatibility system of GP, there are synergistic and balancing interactions among different medicinal ingredients(Wang et al., 2020).
Network pharmacology is a technology based on systems biology that constructs networks for biological systems, focusing on the synergistic effects of multi-component, multi-target, and multi-pathway signals through collaborative modulators of core targets and establishing connections between various components and targets(Noor et al., 2022;Zhao et al., 2023). Based on the identified component of GP, the application of original pseudo-time network pharmacology, molecular docking and molecular dynamics (MD) techniques, explores the main active components in GP that could effectively against PRRSV-induced lung injury and its potential mechanisms.
2. Materials and Methods
2.1 Materials
The medicinal materials, Cinnamon Twig, Peony Root, Licorice, Ginger, and Jujube were all purchased from Zhang Zhongjing Pharmacy(ZhengZhou). Chromatographic grade methanol, acetonitrile, ethanol and formic acid were all purchased from Merck, Germany. Purified water was prepared by Minipore (≥18.2Ω). Other reagents for the experiments were all of analytical purity and purchased from Sinopharm Group.
2.2 Methods
2.2.1 Extraction, identification, analysis of active components in GP and prediction of related targets
2.2.1.1 preparation of sample
According to the classic proportion of Cinnamon Twig, Peony Root, Licorice, Ginger, and Jujube in Shang Han Lun, the herbs were mixed and ground. 50.0g of the mixed sample was taken and distilled under rotation at 50℃ with 100.0 mL of ethanol-1% formic acid-water (55:45, v/v) for 2.0 h. Then, it was concentrated under reduced pressure until nearly dry. The sample was re-dissolved in 50ml of ethanol-1% formic acid-water (55:45, v/v), and ultrasonicated for 30 min to obtain the crude extracts. 5.0 mL of the crude extract was taken, centrifuged at 10000 rpm for 10 min, and the supernatant was collected. It was re-dissolved in 5.0 mL of ethanol-1% formic acid-water (55:45, v/v), ultrasonicated for 30 min, and centrifuged again at 10000 rpm for 10 min. The two supernatants were mixed and concentrated under reduced pressure. After drying, the refined extract was obtained. It was re-dissolved in 5.0 mL of methanol-0.1% formic acid-water (20:80, v/v), 1.0 mL of the re-dissolved solution was taken with a 1.0 mL syringe, and passed through a 0.22μm microporous membrane. The preparation of sample was ready for UPLC-Q/TOF-MS analysis. Meanwhile, the blank control sample was prepared.
2.2.1.2 UPLC-Q/TOF-MS analysis
UPLC conditions: Mobile phase A (aqueous phase) was 0.1% formic acid - water solution, and B (organic phase) was acetonitrile. The gradient elution program of the mobile phase: 0 - 5 min, 5% B; 5.1 - 60 min, 5% - 90% B; 60 - 68 min, 90% - 95% B; 68.1 - 70 min, 95% - 5% B. Flow rate was 0.3 mL/min, column temperature was 45℃, injection volume was 5.0 μL. The chromatographic column was Waters ACQUITY UPLC BEH C 18 (2.1×150 mm, 1.7 μm), and the pressure pump was Shimadzu LC-30AD system.
Mass spectrometry conditions: AB SCIEX X500R mass spectrometer was employed coupled with electrospray ionization(ESI) source. Both positive and negative ion mode scanning were designed. The essential parameters set as follows: The positive and negative ion voltages were 5500 V and -4500 V respectively, mass range m/z 100-1500, collision energy set at 45 eV, centroid data acquisition mode. Capillary voltage and cone voltage were set at 3.0 kV and 35 V, respectively. Ion source temperature and desolvation temperature were maintained in 550 °C and 350 °C, respectively. Nebulizer gas flow rate and cone gas flow rate were kept 3.0 L/min and 50 L/min, respectively.
2.2.1.3 Mass spectrometry data analysis
An efficient peak detection and identification strategy incorporating the NIST20 spectral libraries(Bittremieux et al., 2022;Kubo et al., 2023) was adopted. Firstly, the raw data obtained from UPLC-Q/TOF-MS were converted to mzML format using the Proteo Wizard MS Convert tool(Adusumilli and Mallick, 2017); subsequently, the blank sample data and GP’s sample data were imported into MZmine3 software for batch processing(Heuckeroth et al., 2024). In terms of parameter settings, the minimum peak height was set to 1000, the minimum number of scans required was 5, the MS 1 m/z tolerance was set to 10 ppm, the MS 2 m/z tolerance was also set to 10 ppm, the signal-to-noise ratio threshold was set to 10, and the peak width range was set between 0.1 and 0.3 minutes for deconvolution. Then, isotope grouping was performed, followed by alignment, gap filling, batch correction, and normalization using TIC method. After completing the above steps, feature peaks were extracted to form a feature list containing m/z (mass-to-charge ratio), retention time, and peak area. Third, the MS Search module in NIST20 spectral libraries was used for component matching, with a Score Threshold (score threshold) ≥700 as the screening criterion, ultimately obtaining the mass spectrometry information of the main active components in GP.
2.2.1.4 Prediction of targets of active components
The target prediction was performed using the PharmMapper(Wang et al., 2017) (https://www.lilab-ecust.cn/pharmmapper/), TCMSP(Ru et al., 2014) (https://www.tcmsp-e.com/), TCMID(Huang et al., 2018) (https://bidd.group/TCMID/) online databases. The SMILES structures of each active ingredient were used as input parameters, with the species set to Sus scrofa. Potential targets with a predicted score greater than or equal to 0.5 were screened. The prediction results were integrated, duplicate targets were removed, and an interaction network between the active ingredients of GP and potential target sites was constructed.
2.2.2 Mining of targets related to PRRS
Gene expression profile data of PRRSV infection of GEO databases (https://www.ncbi.nlm.nih.gov/geo/) with the precise search criterion as follow: search “PRRSV” as keyword and “Sus scrofa” as target species for screening. The GSE69515 data set was chosen finally that the temporal variation of gene expression in pigs after infection PRRSV were comprehensively recorded(Schroyen et al., 2015). Based on sample clinical data, 200 samples were classified into 4 experimental groups by 4 different phenotypic traits and further into 3 key time points (0 dpi, 4 dpi and 7 dpi) according to post infection time course. In the analysis, we applied an R language bioinformatic analysis pipeline which performs standardization on the raw gene expression data (preprocessing: Entrez ID convert, missing value imputation and normalizing) and statistical analysis with the DESeq2 package(Love et al., 2014). Using conservative screening criteria (|log2FoldChange|>1 and padj<0.05). Significantly differentially expressed genes (DEGs) in 4 phenotypic groups, respectively, were obtained in the 2 times’ comparison of 4 days vs 0 days and 7 days vs 0 days. Finally, by KEGG signaling pathway analysis, we systematically discovered the biological functional features of these DEGs in the viral infection process and the involved signal transduction pathways, thus comprehensively interpreting the molecular mechanism of PRRSV infection.
2.2.3 Recruitment of targets for lung injury of swine
Using ”lung injury” as the keyword, search and screen in the GeneCard database(Hauser et al., 2003) (https://www.genecards.org/), PharmGKB database(Thorn et al., 2013) (http://www.pharmgkb.org/), TTD database(Chen et al., 2002) (https://idrblab.net/ttd/), and DrugBank database (https://go.drugbank.com/) respectively. After summarizing the obtained targets and removing duplicates, verify the targets through the Uniprot database(Consortium, 2022) (https://www.uniprot.org/) and convert and standardize them with ”Sus scrofa”.
2.2.4 Pseudo-time cross-analysis of the active components targets, targets related to PRRS and lung injury
The active component targets, the DEGs of GSE69515, and the targets of porcine lung injury were subjected to time-series cross-analysis for the 0-4d and 0-7d time periods. The overlapping targets were presented using a Venn diagram(Gao et al., 2021).
2.2.5 Molecular docking
Retrieve and download the three-dimensional crystal structure of the target proteins and its corresponding amino acid sequence from the Uniprot. Homology modeling and preprocessing of the protein structure were carried out using Discovery Studio software (version 2019). This included supplementing missing amino acid residues, removing water molecules and ligand molecules from the structure, adding hydrogen to the protein, and adjusting the charge balance. Convert the processed protein structure into a PDBQT format file. Meanwhile, modify the structure of the small molecules and save it as a PDBQT format file. After setting the docking parameters, run the AutoDock Vina software(Trott and Olson, 2010) to perform a molecular docking simulation. Finally, screen out the active ingredient and target combinations with high binding affinity to the targets.
2.2.6 Molecular dynamics simulation and trajectories analysis
Using the GROMACS (v2020.06) software(Van Der Spoel et al., 2005), MD simulations were performed on the receptor-ligand complexes with better binding effects in molecular docking. First, the PyMOL software(Yuan et al., 2017) was used to process the ligand and receptor separately to obtain the optimal conformation files of the ligand and receptor. Then, the sobtop1.0 software (http://sobereva.com/soft/Sobtop) was used for processing to obtain the structural file and topology file of ligand, and subsequently, the relevant information of ligand was added to the protein topology file to construct the simulation system of the complex. Centered on this complex, a closed model with a 1.2 nm buffer distance was constructed. The Amber94 force field was used, the TIP3P water model was added, and Na⁺ and Cl⁻ were added to balance the charge, making the system electrically neutral.
We performed the constructed system with an energy minimization until the maximum force of the system converged. The detail parameters of the MD simulation were as follows : in any cases with the Hydrogen bonds, the constraints were performed by LINCS algorithm(Hess et al., 1998), the integration step size was 2 fs. And the electrostatics were computed by PME (Particle-mesh Ewald) method whose cutoff value was 1.2 nm.We defined the cutoff value for non-bonded interactions as 10 and updated every 10 steps. We used the V-rescale temperature coupling method to control the simulation temperature at 300 K and the Berendsen method to control the pressure at 1 bar. We performed 1st NVT (constant volume and temperature) and NPT (constant pressure and temperature) equilibrium simulations and then again 100 ns of simulation on the complex system.
2.2.7 Calculate the Gibbs free energy
Poisson-Boltzmann surface area based binding free energy based MD was applied to assess the amount of binding strength variation between conformational cluster, as a post processing of MD trajectories.Using the MM-PBSA method(Genheden and Ryde, 2015) to calculate the binding free energy of the ligand-receptor complex, it provides a detailed characterization of the contributions of each amino acid residue for the binding energy and investigates the dynamic mechanism of interaction between the ligand and receptor.
Results
3.1 Identification of Active Components in GP .
The mass spectrometry data analysis strategy we constructed has been successfully applied in qualitative identification of active components in GP. The detection results in positive and negative ion modes were integrated, and the detected spectral results were carefully annotated and merged through Origin softeware(Fig.1). It was indicated that many active components can be obtained in positive ion mode, and relatively less active components can be identified in negative ion mode.Combined with the compound-matching database in NIST, eventually more than 82 different chemical components were identified including flavonoids, triterpenoids, phenolic acids, coumarins and alkaloids were presented in Table 1.
A total of 82 reliable compounds were identified, with retention times(Rt) ranging from 3.21 - 66.45 min and measured m/z values differing from theoretical values by ≤ |8 ppm|, of which 90% had errors < 3 ppm, fully ensuring the accuracy of molecular formula assignment(Table 1). The polar segment (3-15 min) was dominated by flavonoid-O-glycosides and phenolic acids, with representative components being quercetin derivatives (No. 1, C 17 H 14 O 7, Rt 3.21 min, [M-H]⁻ 331.0838, error 7.76 ppm) and paeoniflorin (No. 13, C 23 H 28 O 11, Rt 13.72 min, [M+H] + 481.1704, error -0.08 ppm); the moderately polar segment (15-30 min) was concentrated with 31 isoflavones, flavan-3-ols, and coumarins, such as formononetin (No. 55, C 16 H 12 O 4, Rt 29.23 min, [M+H]⁺ 269.0808, error -0.15 ppm) and calycosin (No. 34, C 16 H 12 O 5, Rt 22.04 min, [M+H]⁺ 285.0757), whose MS² fragments m/z 269→253/237 and 285→270/255 showed typical pathways of methyl and methyl+water loss; the highly hydrophobic segment (>50 min) was occupied by triterpenoid saponin aglycones, sterols, and carotenoids, with 18α-hydroxyglycyrrhetic acid (No. 51, C 30 H 46 O 5, Rt 27.11 min, [M+H]⁺ 487.3418, error 0.1 ppm) and β-sitosterol (No. 82, C 29 H 50 O, Rt 66.45 min, [M+H]⁺ 415.3908, error 0.96 ppm) producing m/z 397.381, 273.2576, 213.23, 129.699 series fragments through consecutive dehydration (-18 Da) and side-chain cleavage (-C 2 H 4, -C 3 H 6 ); β-carotene (No. 80, C 40 H 56, Rt 60.12 min, [M+H]⁺ 537.4457, error 0.37 ppm) exhibited characteristic fragmentation of polyene chains (m/z 519.435, 429.328, 177.1633, 137.1327). The overall chromatographic-mass spectrometric characteristics revealed that GP contains a ”flavonoid-coumarin-triterpenoid-sterol” four-level polarity spectrum, combining high-polarity glycosides with low-polarity aglycones/sterols, providing precise chemical coordinates for subsequent research. Nevertheless, the isomeric selectivity here is low.
3.2 Application of A Novel Pseudo-time Network Pharmacology
3.2.1 GP Target Prediction, Lung Injury Target Prediction, and Mining of PRRSV-Related Targets
The targets information of the active components of GP was obtained from the the PharmMapper, TCMSP and TCMID databases, and 356 potential targets were screened using the Pharmmapper prediction tool (threshold ≥ 0.5). Further screening based on the criteria of OB ≥ 20% and DL ≥ 0.1, followed by standardization (species: Sus scrofa) and deduplication in the Uniprot database, yielded 307 action targets.
The databases GeneCard, PharmGKB, TTD, and DrugBank were searched to retrieve pig lung injury-related targets. After deduplication and standardization (species: Sus scrofa) in the Uniprot database, a total of 9334 targets were obtained, with GeneCard contributing the most and covering the other databases.
The GSE69515 datasets was restructured into four groups: HH0-4, LL0-4, HH0-7, and LL0-7. The whole blood gene expression after viral infection was analyzed (Fig.2 A1, B1, C1, D1). The DESeq2 package was used to screen for significantly DEGs (Fig.2 A2, B2, C2, D2). The results showed that at 4d vs 0d, the significantly different genes in the HH0-4 and LL0-4 groups were 70 and 120, respectively; at 7dpi, the significantly different genes in the HH0-7 group decreased to 17 (14 downregulated), and in the LL0-7 group to 31 (5 upregulated), both significantly reduced compared to 4dpi, possibly related to disease progression causing interference in the lungs and circulatory system (Fig.2 A3, B3, C3 and D3). The dynamic changes in these DEGs suggested their important roles in viral replication, immune response, and lung injury.
3.2.2 Analysis of Intersection Enrichment for Related Target Genes and KEGG Signaling Pathways
This study utilized the ”clusterProfiler v4.16.0” R package(Wu et al., 2021) and the ”org.Ss.eg.db” database (https://bioconductor.org) to perform KEGG analysis on 307 active component targets of GP, 9334 pig lung injury targets, and the DEGs in the HH 0-4, LL 0-4, HH 0-7, and LL 0-7 groups, identifying 157 common key pathways involving apoptosis, infection mechanisms, oxygen metabolism, and inflammatory responses, confirming the ”regulating gas( qi ) and activating blood” efficacy of GP and the important role of the lungs (Fig.3 B).
The number of intersecting common pathways between drugs and lung injury and the four groups of DEGs varied: 15 in the HH 0-4 group, 21 in the LL 0-4 group, none in the HH 0-7 group (possibly due to insufficient DEGs), and 21 in the LL 0-7 group(Fig.3 C 1-4). In the early stages of viral infection (HH 0-4 and LL 0-4), GP significantly affected key pathways, which weakened in the later stages (HH 0-7), suggesting its better early intervention effects; the consistently present 9 immune- and apoptosis-related pathways indicated its stable regulatory role(Fig.3 A).
Target gene intersection analysis showed that only the LL 0-4 group had 3 common target genes (JAK2, CASP3, and TLR2), the HH 0-4 group had 1 (CASP3), and the key gene shared by LL 0-4 and HH 0-4 was CASP3. GP had a stronger direct intervention effect on early-stage (LL 0-4) lung injury caused by PRRSV infection than in the later stages, with more significant effects on the low-virulence group (LL). JAK2, CASP3, and TLR2 were involved in cell proliferation, apoptosis, and pathogen recognition, respectively, indicating its multi-target regulation of immune and apoptotic processes, providing a theoretical basis for antiviral therapy.
3.3 Molecular docking and molecular dynamics simulation
3.3.1 Molecular docking results
Pseudo-time network pharmacology analysis showed that GP has potential intervention effects on early lung injury caused by PRRSV. A re-extraction and analysis of the GSE69515 found JAK2 was upregulated in the HH 0-4 group, CASP3 was highly expressed in the LL 0-4 group, and TLR2 was lowly expressed in both groups(Fig. 4). Based on the network, 7 active components related to these three targets were screened out: β-carotene, farnesol, DVF, kaempferol, naringenin, quercetin and β-sitosterol. Molecular docking showed that the binding energies of β-carotene with JAK2, CASP3, and TLR2 were -11.6, -11.6, and -15.5 kcal/mol (all below the -7 kcal/mol threshold), indicating strong affinity; the binding energy of β-sitosterol with JAK2 was -11.8 kcal/mol, and that of β-carotene with TLR2 was -9.8 kcal/mol, which may respectively inhibit viral proliferation though enhancing pathogen recognition via activating the expression level of TLR2.
We chosen 6 pairs of protein-ligand complex whose binding energies are smallest and mainly relative to CASP3, JAK2 and TLR2. All plots in Fig.5 indicated that the Van der Waals forces were the primary binding force for the complex stability and the two auxiliary binding forces were hydrogen bonds and hydrophobic interactions. Further, We carried on the MD simulations to validate the complex stability of CASP3-β-sitosterol/β-carotene -9.9/-10.2 kcal/mol, and −9.1/-15.7 kcal/mol for TLR2-farnesol / β-sitosterol complexes, respectively.
3.3.2 Molecular dynamics simulations results
Employing the GROMACS 2020.06 software to systematically investigate the dynamic behavior of 4 protein-ligand complex systems over a 100 ns timescale. The studied systems included: the CASP3-β-carotene complex (Complex I), the CASP3- β-sitosterol complex (Complex II), the TLR2-farnesol complex (Complex III), and the TLR2-β-sitosterol complex (Complex IV). During the simulations, MD trajectory conformations of each system were extracted at the 20 ns, 50 ns, 80 ns, and 100 ns time points(Fig.6). The overall analysis results showed that the ligand in all 4 complex systems could effectively bind within the corresponding protein binding pockets, which was entirely consistent with the trend observed in the preliminary molecular docking studies where β-carotene, β-sitosterol, and farnesol exhibited good binding energies with CASP3 and TLR2.
In order to precisely characterize the dynamic properties of each complex system, we calculated and compared some structural dynamics parameters, root-mean-squaredeviation (RMSD), root-mean-square fluctuation (RMSF), radius of gyration (Rg), solvent accessible surface area (SASA), principal component analysis (PCA), hydrogen-bond formation, of Complexes I to IV. Among them, the RMSD is the essential parameter characterizing system stability directly reflecting the intensity of molecules motion and showing that higher RMSD values are indicative of relatively harsh molecular motions and unfavorable system stabilities and that lower values of RMSD are indicative of relative system stabilities.In this work we utilized the initial protein structures of each system as baselines, and investigated the trends in RMSD dynamics of protein backbone atoms during the 100 ns simulation duration along with the protein conformation stabilities and ligand’s positional fluctuations with respect to the protein. Subplots A, B, C, and D of Fig.7 illustrated the RMSD dynamic evolution of Complexes I to IV for the 0-100ns duration.It can be seen from subplots A and B in Figure 7 that protein backbone of CASP3 has reached a steady state after the initial fluctuation and it was relatively stable, while Complexes I and II might have been bound unstably by the β-carotene and β-sitosterol, respectively, since their RMSD curves reflected typical transitional fluctuation state. Interestingly, throughout most of the simulation time, the protein RMSD of Complexes I and II were approximately overlapped with the overall system RMSD, indicating a good stability of the overall system. It might be due to the ligands binding to surface pocket of CASP3 and they had a lower stability, which was also speculated from the trajectory conformation plots A and B. In short, they kept stable in their binding states. As can be seen from plots C and D in Fig.7, in the subsequent relaxation, the RMSD fluctuation magnitudes of the farnesol / β-sitosterol-TLR2 complexes substantially reduced (resulting in more stable structure binding).
Subplots E, F, G and H of Figure 7 showed the results of RMSF analysis for Complexes I, II, III and IV, respectively. The RMSF distributions can reveal the sensitively the local vibration features of protein atoms near the ligand-binding sites. In our analysis, it indicated that the β-carotene caused greater RMSF fluctuations in the regions of 1-40 and 160-200 amino acids of the CASP3 protein (Fig.7 E), and β-sitosterol caused visible fluctuations peaks in the region of 1-40 of the CASP3 protein (Fig.7 F).These results complement the stability feature of RMSD calculated Complex I and Complex II, since the effects of ligands perturbation were identified to drive the dynamical pattern of the corresponding part of their RMSF. However, the RMSF fluctuation peaks of farnesol and β-sitosterol in the TLR2 protein were relatively small but the perturbation frequency of its each individual amino acid residues was more pronounced(Fig.7 G and H).
In the effects of CASP3-TLR2 secondary structure interactions, proteins’ overall spatial tightness has very obviously changed. To determine the structural variation quantitatively, we used the Rg parameters that are useful physics quantities reflecting the protein’s structural tightness of the Complexes I to IV. Fig7 subplots I, J, K and L indicated that CASP3 and TLR2 protein in Complexes I to IV during a 100 ns simulation maintained relatively unchanged folded conformations although there were some radial movements along the X, Y and Z spatial axes, the Rg values of all complexes stayed generally unchanged.
The SASA is an important parameter for research into the conformational dynamics of protein in solvent environment. By SASA analysis, we systematically analyzed the variations of the protein solvent-exposed surface area. The analysis results indicated that the total SASA of the Complex I presented a convergent trend of variation between 130-180 nm; the varying pattern of the Complex II was nearly the same as Complex I. Contrarily, the overall SASA of Complexes III and IV did fluctuate between the limits of 245-265 nm and did not change throughout the simulation(Fig.8 A, B, C, and D).
In order to further understand the internal dynamics features of Complexes I to IV, we build up a PCA cross-correlation matrix composed by the fluctuations of Cα backbone atoms of CASP3 and TLR2. The 2D free energy landscape(FEL) topography was constructed. Here, red areas reflect a relatively high-energy and weak stability state of the system, and blue areas reflect a relatively low-energy and high stability conformations. Comparable with results of RMSD analysis of Complexes I and II, the FEL plots of Complexes I and II exhibited more red regions while the PCA plots of Complexes III and IV showed larger and more gathered blue regions, implying that overall two complex systems were more structurally stable(Fig.8 E, F, G and H).
The secondary interaction – hydrogen bonds play an important role in protein-ligand complex system for keeping protein-ligand binding interface stable, and the dynamic change characteristics of hydrogen bond directly indicate the stability and specificity of protein-ligand binding interface. Hydrogen bond number formed in simulation was calculated. It was found that there was no obvious hydrogen bond formation in complex I in the process of hydrogen bond analysis(Fig.8 I).The hydrogen bond number distribution for the Complexes II, III, and IV was displayed in Fig8 subplots J, K and L, and the results of analysis showed that Complex II built a more extensive hydrogen-bonding network, which suggests that the stability of the complex largely originates from the hydrogen-bond interaction, while possibly there was another sort of non-covalent interaction between the Complexes III and IV, and these synergic effects are to sustain the well-binding stability of the systems.
3.3.3 Gibbs free energy analysis
In order to explore the variation of Gibbs binding free energy for the complex II-IV formation processes, the MM-PBSA approach was applied here to evaluate the binding free energy (Table 3) and contributions of different energy terms throughout a 100 ns MD simulations.
The calculation results showed that the total binding free energy (ΔGtotal) of Complex I is -34.31±2.53 kcal/mol, that of Complex II is -39.46±0.64 kcal/mol, while those of Complex III and IV are -28.34±1.49 kcal/mol and -48.60±0.57 kcal/mol, respectively. These values are consistent with the trends obtained from molecular docking studies, thereby confirming the rationality of the MD simulations performed for the ComplexⅠ II, III, and IV systems.
References
Adusumilli, R., and Mallick, P. (2017). Data Conversion with ProteoWizard msConvert. Methods Mol Biol 1550 , 339–368. DOI: 10.1007/978-1-4939-6747-6_23.Armour, M., Al-Dabbas, M.A., Ee, C., Smith, C.A., Ussher, J., Arentz, S., Lawson, K., and Abbott, J. (2021). The effectiveness of a modified Gui Zhi Fu Ling Wan formulation (Gynoclear™) for the treatment of endometriosis: a study protocol for a placebo-controlled, double-blind, randomised controlled trial. Trials 22 , 299. DOI: 10.1186/s13063-021-05265-x.Bello-Onaghise, G., Wang, G., Han, X., Nsabimana, E., Cui, W., Yu, F., Zhang, Y., Wang, L., Li, Z., Cai, X., and Li, Y. (2020). Antiviral Strategies of Chinese Herbal Medicine Against PRRSV Infection. Front Microbiol 11 , 1756. DOI: 10.3389/fmicb.2020.01756.Bittremieux, W., Wang, M., and Dorrestein, P. (2022). The critical role that spectral libraries play in capturing the metabolomics community knowledge. DOI: 10.26434/chemrxiv-2022-hjtt2.Chen, X., Ji, Z.L., and Chen, Y.Z. (2002). TTD: Therapeutic Target Database. Nucleic Acids Res 30 , 412–415. DOI: 10.1093/nar/30.1.412.Consortium, T.U. (2022). UniProt: the Universal Protein Knowledgebase in 2023. Nucleic Acids Research 51 , D523–D531. DOI: 10.1093/nar/gkac1052.Gao, C.-H., Yu, G., and Cai, P. (2021). ggVennDiagram: An Intuitive, Easy-to-Use, and Highly Customizable R Package to Generate Venn Diagram. Frontiers in Genetics Volume 12 - 2021. DOI: 10.3389/fgene.2021.706907.Genheden, S., and Ryde, U. (2015). The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities. Expert Opin Drug Discov 10 , 449–461. DOI: 10.1517/17460441.2015.1032936.Hauser, E.R., Mooser, V., Crossman, D.C., Haines, J.L., Jones, C.H., Winkelmann, B.R., Schmidt, S., Scott, W.K., Roses, A.D., Pericak-Vance, M.A., Granger, C.B., and Kraus, W.E. (2003). Design of the Genetics of Early Onset Cardiovascular Disease (GENECARD) study. Am Heart J 145 , 602–613. DOI: 10.1067/mhj.2003.13.Hess, B., Bekker, H., Berendsen, H., and Fraaije, J. (1998). LINCS: A Linear Constraint Solver for molecular simulations. Journal of Computational Chemistry 18. DOI: 10.1002/(SICI)1096-987X(199709)18:123.0.CO;2-H.Heuckeroth, S., Damiani, T., Smirnov, A., Mokshyna, O., Brungs, C., Korf, A., Smith, J.D., Stincone, P., Dreolin, N., Nothias, L.-F., Hyötyläinen, T., Orešič, M., Karst, U., Dorrestein, P.C., Petras, D., Du, X., Van Der Hooft, J.J.J., Schmid, R., and Pluskal, T. (2024). Reproducible mass spectrometry data processing and compound annotation in MZmine 3. Nature Protocols 19 , 2597–2641. DOI: 10.1038/s41596-024-00996-y.Holtkamp, D.J., Kliebenstein, J.B., and Neumann, E.J. (2013). Assessment of the economic impact of porcine reproductive and respiratory syndrome virus on United States pork producers. Journal of Swine Health and Production 21 , 72–84. https://www.aasv.org/jshap/jshap-abstract/?v21n2p72Huang, L., Xie, D., Yu, Y., Liu, H., Shi, Y., Shi, T., and Wen, C. (2018). TCMID 2.0: a comprehensive resource for TCM. Nucleic Acids Res 46 , D1117–d1120. DOI: 10.1093/nar/gkx1028.Jain, M., Muthukumaran, J., and Singh, A.K. (2021). Structural and functional characterization of chitin binding lectin from Datura stramonium: insights from phylogenetic analysis, protein structure prediction, molecular docking and molecular dynamics simulation. Journal of Biomolecular Structure and Dynamics 39 , 1698–1716. DOI: 10.1080/07391102.2020.1737234.Kubo, A., Kubota, A., Ishioka, H., Hizume, T., Ubukata, M., Nagatomo, K., Satoh, T., Yoshida, M., and Uematsu, F. (2023). Construction of a Mass Spectrum Library Containing Predicted Electron Ionization Mass Spectra Prepared Using a Machine Learning Model and the Development of an Efficient Search Method. Mass Spectrometry 12 , A0120–A0120. DOI: 10.5702/massspectrometry.A0120.Kumari, G., Nigam, V.K., and Pandey, D.M. (2023). The molecular docking and molecular dynamics study of flavonol synthase and flavonoid 3’-monooxygenase enzymes involved for the enrichment of kaempferol. Journal of Biomolecular Structure and Dynamics 41 , 2478–2491. DOI: 10.1080/07391102.2022.2033324. Li, X., Chen, S., Zhang, L., Niu, G., Zhang, X., Yang, L., Ji, W., and Ren, L. (2022). Coinfection of Porcine Circovirus 2 and Pseudorabies Virus Enhances Immunosuppression and Inflammation through NF-κB, JAK/STAT, MAPK, and NLRP3 Pathways. International Journal of Molecular Sciences 23 , 4469. https://www.mdpi.com/1422-0067/23/8/4469.Liang, Q., Jiang, Y., Liang, L., and Liang, H. (2025). Factors Influencing the Efficacy of Chinese Medicinals. Chinese medicine and natural products 05 , e73–e85. DOI: 10.1055/s-0045-1809610.Ling, W., Li, X., Liu, X., Wu, Q., Wang, W., Meng, J., Sun, B., and Lv, B. (2023). Mechanism of the Anti-Influenza Functions of Guizhi Granules Based on Network Pharmacology, Molecular Docking, and in Vitro Experiments. Chemistry & Biodiversity 20 , e202201228. DOI: 10.1002/cbdv.202201228. Liu, X., Chen, L., Sun, P., Zhan, Z., and Wang, J. (2024). Guizhi Fuling Formulation: A review on chemical constituents, quality control, pharmacokinetic studies, pharmacological properties, adverse reactions and clinical applications. Journal of Ethnopharmacology 319 , 117277. DOI: 10.1016/j.jep.2023.117277.Love, M.I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15 , 550. DOI: 10.1186/s13059-014-0550-8. Merca, C., Lindell, I.C., Ernholm, L., Selling, L.E., Nunes, T.P., Sjölund, M., and Dórea, F.C. (2022). Veterinary syndromic surveillance using swine production data for farm health management and early disease detection. Preventive Veterinary Medicine 205 , 105659. DOI: 10.1016/j.prevetmed.2022.105659. Mil-Homens, M., Silva, G., Holtkamp, D., Linhares, D., Osemeke, O., Dion, K., Baker, K., Robbins, R., Sparks, J., Jensen, R., Arruda, A., Corzo, C., Vanderwaal, K., Kikuti, M., Yeske, P., Glowzenski, L., Gillespie, T., Petznick, T., and Lopez, W. (2025). Proposing a clinical case definition for porcine reproductive and respiratory syndrome virus outbreaks in sow herds based on key productivity indicators. Journal of Swine Health and Production 33 , 153–160. https://www.aasv.org/shap/issues/v33n4/v33n4p153.html.Noor, F., Asif, M., Ashfaq, U.A., Qasim, M., and Tahir Ul Qamar, M. (2023). Machine learning for synergistic network pharmacology: a comprehensive overview. Briefings in Bioinformatics 24. DOI: 10.1093/bib/bbad120.Noor, F., Tahir Ul Qamar, M., Ashfaq, U.A., Albutti, A., Alwashmi, A.S.S., and Aljasir, M.A. (2022). Network Pharmacology Approach for Medicinal Plants: Review and Assessment. Pharmaceuticals 15 , 572. https://www.mdpi.com/1424-8247/15/5/572. Pedro Mil-Homens, M., Jayaraman, S., Rupasinghe, K., Wang, C., Trevisan, G., Dórea, F., C. L. Linhares, D., Holtkamp, D., and S. Silva, G. (2024). Early Detection of PRRSV Outbreaks in Breeding Herds by Monitoring Productivity and Electronic Sow Feed Data Using Univariate and Multivariate Statistical Process Control Methods. Transboundary and Emerging Diseases 2024 , 9984148. DOI: 10.1155/2024/9984148.Ru, J., Li, P., Wang, J., Zhou, W., Li, B., Huang, C., Li, P., Guo, Z., Tao, W., Yang, Y., Xu, X., Li, Y., Wang, Y., and Yang, L. (2014). TCMSP: a database of systems pharmacology for drug discovery from herbal medicines. Journal of Cheminformatics 6 , 13. DOI: 10.1186/1758-2946-6-13. Sánchez-Carvajal, J.M., Ruedas-Torres, I., Carrasco, L., Pallarés, F.J., Mateu, E., Rodríguez-Gómez, I.M., and Gómez-Laguna, J. (2021). Activation of regulated cell death in the lung of piglets infected with virulent PRRSV-1 Lena strain occurs earlier and mediated by cleaved Caspase-8. Veterinary Research 52 , 12. DOI: 10.1186/s13567-020-00882-x.Schroyen, M., Steibel, J.P., Koltes, J.E., Choi, I., Raney, N.E., Eisley, C., Fritz-Waters, E., Reecy, J.M., Dekkers, J.C., Rowland, R.R., Lunney, J.K., Ernst, C.W., and Tuggle, C.K. (2015). Whole blood microarray analysis of pigs showing extreme phenotypes after a porcine reproductive and respiratory syndrome virus infection. BMC Genomics 16 , 516. DOI: 10.1186/s12864-015-1741-8.Shao, C., Yu, Z., Luo, T., Zhou, B., Song, Q., Li, Z., Yu, X., Jiang, S., Zhou, Y., Dong, W., Zhou, X., Wang, X., and Song, H. (2022). Chitosan-Coated Selenium Nanoparticles Attenuate PRRSV Replication and ROS/JNK-Mediated Apoptosis in vitro. International Journal of Nanomedicine 17 , 3043–3054. DOI: 10.2147/IJN.S370585.Singh, S., Kumar, R., Payra, S., and Singh, S.K. (2023). Artificial Intelligence and Machine Learning in Pharmacological Research: Bridging the Gap Between Data and Drug Discovery. Cureus 15 , e44359. DOI: 10.7759/cureus.44359.Thorn, C.F., Klein, T.E., and Altman, R.B. (2013). PharmGKB: the Pharmacogenomics Knowledge Base. Methods Mol Biol 1015 , 311–320. DOI: 10.1007/978-1-62703-435-7_20.Trott, O., and Olson, A.J. (2010). AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J Comput Chem 31 , 455–461. DOI: 10.1002/jcc.21334.Van Der Spoel, D., Lindahl, E., Hess, B., Groenhof, G., Mark, A.E., and Berendsen, H.J. (2005). GROMACS: fast, flexible, and free. J Comput Chem 26 , 1701–1718. DOI: 10.1002/jcc.20291.Wang, H., and Feng, W. (2024). Current Status of Porcine Reproductive and Respiratory Syndrome Vaccines. Vaccines (Basel) 12. DOI: 10.3390/vaccines12121387. Wang, J., Li, X., Bello, B.K., Yu, G., Yang, Q., Yang, H., Zhang, W., Wang, L., Dong, J., Liu, G., and Zhao, P. (2022). Activation of TLR2 heterodimers-mediated NF-κB, MAPK, AKT signaling pathways is responsible for Vibrio alginolyticus triggered inflammatory response in vitro. Microbial Pathogenesis 162 , 105219. DOI: 10.1016/j.micpath.2021.105219.Wang, Q., Xu, B., Wang, F., Xu, F.F., Zhang, X., Zhang, Y.C., Du, H., Xia, C.Y., Bao, L.W., Wang, Z.Z., Qiao, Y.J., and Xiao, W. (2020). [Predictive model for hygroscopicity of contents in Guizhi Fuling Capsules]. Zhongguo Zhong Yao Za Zhi 45 , 242–249. DOI: 10.19540/j.cnki.cjcmm.20191219.302.Wang, X., Shen, Y., Wang, S., Li, S., Zhang, W., Liu, X., Lai, L., Pei, J., and Li, H. (2017). PharmMapper 2017 update: a web server for potential drug target identification with a comprehensive target pharmacophore database. Nucleic Acids Res 45 , W356–w360. DOI: 10.1093/nar/gkx374. Wang, X., Wang, Z.-Y., Zheng, J.-H., and Li, S. (2021). TCM network pharmacology: A new trend towards combining computational, experimental and clinical approaches. Chinese Journal of Natural Medicines 19 , 1–11. DOI: 10.1016/S1875-5364(21)60001-8. Wei, P., Huang, S., Yang, J., Zhao, M., Chen, Q., Deng, X., Chen, J., and Li, Y. (2024). Identification and characterization of chemical constituents in Mahuang Guizhi Decoction and their metabolites in rat plasma and brain by UPLC-Q-TOF/MS. Chinese Herbal Medicines 16 , 466–480. DOI: 10.1016/j.chmed.2024.01.006.Wu, T., Hu, E., Xu, S., Chen, M., Guo, P., Dai, Z., Feng, T., Zhou, L., Tang, W., Zhan, L., Fu, X., Liu, S., Bo, X., and Yu, G. (2021). clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation 2 , 100141. DOI: 10.1016/j.xinn.2021.100141.Yuan, S., Chan, H.C.S., and Hu, Z. (2017). Using PyMOL as a platform for computational drug design. WIREs Computational Molecular Science 7 , e1298. DOI: 10.1002/wcms.1298.Yuan, Z., Sun, Y., Niu, X., Yan, Q., Zeng, W., Du, P., Xie, K., Fang, Y., Wang, L., Ding, H., Yi, L., Zhao, M., Fan, S., Zhao, D., and Chen, J. (2024). Epidemiologic Investigation and Genetic Variation Analysis of PRRSV, PCV2, and PCV3 in Guangdong Province, China from 2020 to 2022. Viruses 16 , 1687. https://www.mdpi.com/1999-4915/16/11/1687.Zamzami, M.A. (2023). Molecular docking, molecular dynamics simulation and MM-GBSA studies of the activity of glycyrrhizin relevant substructures on SARS-CoV-2 RNA-dependent-RNA polymerase. Journal of Biomolecular Structure and Dynamics 41 , 1846–1858. DOI: 10.1080/07391102.2021.2025147. Zhang, H., Luo, J., Wan, Q., Wang, X., Wu, Z., Yang, M., and Wang, Y. (2025). Effect of freeze-pressure regulated extraction technology on the physicochemical properties and pharmacological activities of guizhi extract. Frontiers in Chemistry Volume 13 - 2025. DOI: 10.3389/fchem.2025.1581429.Zhang, H., Xiang, L., Xu, H., Li, C., Tang, Y.-D., Gong, B., Zhang, W., Zhao, J., Song, S., Peng, J., Wang, Q., An, T., Cai, X., and Tian, Z.-J. (2022). Lineage 1 Porcine Reproductive and Respiratory Syndrome Virus Attenuated Live Vaccine Provides Broad Cross-Protection against Homologous and Heterologous NADC30-Like Virus Challenge in Piglets. Vaccines 10 , 752. DOI: 10.3390/vaccines10050752.Zhang, P., Zhang, D., Zhou, W., Wang, L., Wang, B., Zhang, T., and Li, S. (2024). Network pharmacology: towards the artificial intelligence-based precision traditional Chinese medicine. Briefings in Bioinformatics 25. DOI: 10.1093/bib/bbad518. Zhao, L., Zhang, H., Li, N., Chen, J., Xu, H., Wang, Y., and Liang, Q. (2023). Network pharmacology, a promising approach to reveal the pharmacology mechanism of Chinese medicine formula. Journal of Ethnopharmacology 309 , 116306. DOI: 10.1016/j.jep.2023.116306.Zimmerman, J.J., Dee, S.A., Holtkamp, D.J., Murtaugh, M.P., Stadejek, T., Stevenson, G.W., Torremorell, M., Yang, H., and Zhang, J. (2019). ”Porcine Reproductive and Respiratory Syndrome Viruses (Porcine Arteriviruses),” in Diseases of Swine .), 685–708. https://onlinelibrary.wiley.com/doi/abs/10.1002/9781119350927.ch41
Information & Authors
Information
Version history
Copyright
This work is licensed under a Non Exclusive No Reuse License.