Drug information processing method, system and equipment based on HPV L2 transmembrane region hydrophobic pocket and medium

Abstract

The invention provides a drug information processing method, system and equipment based on an HPV L2 transmembrane region hydrophobic pocket and a medium, and belongs to the technical field of drug information processing. Rebuilding a physiological environment in the lipid bilayer; optimizing the structure through molecular docking and dynamic simulation; setting and operating molecular dynamics simulation; analyzing trajectory data and extracting dynamic features of a hydrophobic pocket; identifying a key residue; and calculating combined free energy and combining experimental verification. According to the method, the characteristics of the hydrophobic pocket are accurately represented through dynamic simulation, calculation and experimental verification are combined, the dynamic characteristics of the hydrophobic pocket are accurately revealed through molecular dynamics simulation, and MM-PBSA calculation and SPR experimental verification are combined, so that small molecule drugs targeting an HPV L2 transmembrane region can be more accurately screened and optimized, and the drug research and development efficiency is improved.

Description

Technical Field

[0001] This invention belongs to the field of drug information processing technology, specifically relating to a drug information processing method, system, device and medium based on the hydrophobic pocket of the transmembrane region of HPV L2. Background Technology

[0002] Human papillomavirus (HPV) infection is a major cause of various diseases, including cervical cancer and genital warts. Currently available vaccines primarily target the L1 protein of HPV, but their effectiveness against existing infections is limited. The L2 protein, another important structural protein of HPV, plays a crucial role in the viral invasion process.

[0003] However, in related technologies, when X-ray crystallography is used to resolve the transmembrane region structure of HPV L2, a buffer solution is added, which does not simulate the lipid bilayer environment of the cell membrane. The conformation of transmembrane proteins is highly dependent on the hydrophobic interactions of lipid molecules. For example, the orientation of the transmembrane helix is ​​determined by the hydrophobic compression of the lipid tail. In the absence of lipids, proteins are prone to problems such as deviations in the tilt angle of the transmembrane helix (>10°), closure of hydrophobic pockets due to lack of lipid support, and abnormal orientation of polar residue side chains.

[0004] Inappropriate molecular dynamics simulation parameters and incomplete capture of dynamic features are prevalent issues in related technologies. Specifically, in traditional molecular dynamics (MD) simulations of the HPV L2 transmembrane region, improper system parameter selection and excessively short simulation durations (<50 ns) fail to cover slow dynamic processes such as pocket opening and closing and lipid-protein interaction adjustments. Such parameter deficiencies cause the simulation system to deviate from its physiological state, resulting in distorted dynamic characteristics. Furthermore, short-duration simulations cannot capture critical dynamic processes, leading to subsequent analyses based on incomplete dynamic data and an inability to determine the optimal conformational window for drug binding.

[0005] The identification of key residues in related technologies is based on sequence conservation, making the amino acid sequence conserved. In the transmembrane region of HPV L2, some non-conserved residues, although having different sequences, all form hydrophobic interactions with drugs in dynamics; some conserved residues have side chains facing outwards from the pocket and do not actually participate in drug binding. Summary of the Invention

[0006] This invention provides a drug information processing method based on the hydrophobic pocket of the transmembrane region of HPV L2. The method analyzes the structure and dynamic characteristics of the HPV L2 hydrophobic pocket, providing reliable target structure, dynamic information, and energy basis for targeted drug screening, thereby improving screening efficiency and hit rate and reducing invalid experiments. The methods include: Step S101: Obtain the high-resolution structure of the HPV L2 transmembrane region by X-ray crystallography; Step S102: Use a structural reconstruction tool to reconstruct the physiological environment of the HPV L2 transmembrane protein corresponding to the high-resolution structure described in step S101 in the lipid bilayer; Step S103: Optimize the transmembrane protein structure of HPVL2 after reconstructing the physiological environment in step S102 through molecular docking and molecular dynamics simulation; Step S104: Set the molecular dynamics simulation parameters, including using specific molecular dynamics software and controlling the temperature, pressure and simulation step size to preset parameters under canonical isobaric and isothermal ensemble conditions; Step S105: Based on the preset parameters set in step S104, run the molecular dynamics simulation and record the trajectory data during the simulation process; Step S106: Analyze the trajectory data recorded in step S105, extract the dynamic change characteristics of the hydrophobic pocket in the HPV L2 transmembrane region, and use the pocket volume calculation tool to calculate the volume change information of the hydrophobic pocket. Step S107: Identify key residues in the hydrophobic pocket of step S106 by hydrogen bond analysis and hydrophobicity analysis; Step S108: The binding free energy of the small molecule to the HPV L2 transmembrane protein was calculated using the molecular mechanics-Poisson-Boltzmann surface area method. The binding affinity between the two was determined by surface plasmon resonance experiments, and the good agreement between the calculated results and the experimental results was compared.

[0007] It should be further explained that step S102 specifically includes the following methods: Based on the amino acid sequence characteristics of HPV L2 transmembrane proteins, matching lipid bilayer components were selected; The membrane construction module of the structure reconstruction tool generates the initial topology of the lipid bilayer; The high-resolution structure of the HPV L2 transmembrane protein obtained in step S101 is translated to the central region of the lipid bilayer. The protein is rotated so that the axis of the transmembrane helix is ​​perpendicular to the lipid bilayer plane. The initial spatial position of the protein and lipid is optimized by minimizing energy using the steepest descent method. Adding ions with opposite charges to the lipid bilayer system and calculating the system charge density according to the Poisson-Boltzmann equation neutralizes the overall net charge of the protein and lipids to achieve electroneutrality.

[0008] It should be further explained that step S103 specifically includes the following methods: Step S1031: Assign force field parameters to candidate small molecules and calculate partial atomic charges using a graph theory topology generation algorithm to form a small molecule ligand library compatible with the force field of HPV L2 transmembrane region proteins; Step S1032: Using the ligand library obtained in S1031 as the object, the Fast Fourier Transform correlation algorithm is used to search for shape complementarity and chemical scores between each small molecule and HPV L2 transmembrane protein, generating and sorting multiple binding conformations for stability testing. Step S1033: Perform restricted energy minimization and short-time equilibrium on the highest-scoring binding conformation in S1032 in an explicit solvent environment based on potential energy surface gradient descent; Step S1034: Using the trajectory of the stabilized complex after S1033, calculate the root mean square deviation and root mean square fluctuation, and select the binding conformation with the lowest structural deviation as the optimal representative.

[0009] It should be further explained that step S105 specifically includes the following methods: Step S1051: Verify the gro, top, and mdp files using GROMACS, expand the cubic periodic cell, assign Maxwell-Boltzmann initial velocities, and generate the initial simulation input file; Step S1052: The structure obtained in step S1051 is subjected to the steepest descent and conjugate gradient second-order energy minimization in sequence to reduce the maximum atomic force to below 2kJ / (mol·nm). The system stability is confirmed by the 100ps short-range simulated main chain RMSD≤0.2nm. Otherwise, it is back to S1051 for re-verification. Step S1053: Based on the stable structure in step S1052, first perform 5ns NVT equilibration and then 10ns NPT equilibration, with temperature control at 310±2K, pressure control at 1±0.2atm and main chain position restriction retained, and complete the equilibration after the density reaches about 0.997g / cm³. Step S1054: Remove all position restrictions from step S1053, and use the leap-frog integrator, Verlet neighbor list, 1.0nm intercept and PME long-range processing to write the trajectory every 2ps, the energy every 10ps, and the checkpoint every 100ps, thereby continuously acquiring the production trajectory for 200ns with a step size of 2fs. Step S1055: Verify the frame integrity of the trajectory obtained in step S1054 using gmxcheck. Automatically complete if there are ≤5 missing frames, and rerun the corresponding time period if there are >5 missing frames. Ensure that the main chain RMSD fluctuation of the entire trajectory is within 0.1–0.3nm. Finally, name and archive the trajectory, energy, and checkpoint files and generate a simulation report to complete data storage and quality confirmation.

[0010] It should be further explained that step S106 specifically includes the following methods: Step S1061: The trajectory from step S105 was dehydrated, deionized, and fitted with rmsd using the GROMACStrjconv tool, retaining only the protein atom coordinates of the HPV L2 transmembrane region to obtain the corrected pure protein trajectory. Step S1062: Based on the purified protein trajectory in step S1061, MDAnalysis is used to extract the side chain heavy atoms of leucine, isoleucine, valine, phenylalanine, tryptophan, and methionine as pocket marker atoms to form a set of hydrophobic pocket boundary atoms. Step S1063: Calculate the root mean square fluctuation of the pocket-marked atom set in step S1062 and perform principal component analysis to obtain the pocket flexibility quantization value and main motion modes; Step S1064: Based on the pocket-marked atom set in step S1062, the volume of the hydrophobic pocket in each frame is calculated using the POVME3.0 convex hull algorithm, and the volume time series is output to realize continuous monitoring of the pocket opening and closing degree; Step S1065: Using the volume time series from step S1064 as the object, perform independent copy correlation test and key residue mutation resimulation to verify the reproducibility and causality of the volume calculation and ensure the reliability of the dynamic volume results of the hydrophobic pocket.

[0011] It should be further explained that step S107 specifically includes the following methods: Step S1071: Using PyMOL and GROMACS together, export the numbers of the hydrophobic pocket boundary residues and generate a dedicated trajectory and index file. After verifying the coverage rate by set operation, the pocket residue set required for subsequent analysis is locked. Step S1072: Using the pocket residue set obtained in S1071 as the object, hydrogen bonds are detected frame by frame using the VMD plugin according to the standard that the donor-acceptor distance is less than 3.5 Å and the included angle is greater than 120°. The hydrogen bond occupancy rate is calculated and a network diagram is drawn. Residues with an occupancy rate greater than 20% are screened to form potential hydrogen bond key residues. Step S1073: Based on the pocket residue set in S1071, Hopp-Woods scaling is used for assignment and solvent-accessible surface area calculation is combined with probe radius of 1.4 Å. Residues with high hydrophobicity values ​​and low surface area are included in the hydrophobic core candidate set, which together with hydrogen bond residues constitute a dual-attribute candidate pool. Step S1074: Correlation analysis is performed on the hydrogen bond occupancy rate of S1072 and the hydrophobicity value of S1073 with the RMSF and volume change of S1064. Intersecting residues with hydrogen bond occupancy rate greater than 20%, RMSF less than 0.2nm, and mutual information greater than 0.4 are retained to form a preliminary screening list of key residues. Step S1075: For the residues in the initial screening list of S1074, the key residue list is finally determined and classified and output by PyMOL orientation check, GROMACS minimum distance less than 4Å screening and homologous subtype conservation comparison, with hydrophobicity value and hydrogen bond occupancy rate quantification parameters.

[0012] It should be further explained that step S108 specifically includes the following methods: S1081: The binding free energy was calculated using the gmxmmpbsa tool of GROMACS. The dielectric constant ε=80 was set to simulate the polarity of the water environment. The accessible surface area of ​​the nonpolar solvent was approximated by the LCPO model. The total energy difference ΔGca between the small molecule and the protein complex was calculated. S1082: SPR experiments were performed using a Biacore 8K instrument. HPV L2 protein was immobilized on the surface of a CM5 sensor chip, and small molecule solutions of varying concentrations were injected. The real-time binding response signal was monitored, and the dissociation constant Kde was obtained by fitting a 1:1 Langmuir binding model. S1083: Construct the computation-experimental data correlation matrix, perform thermodynamic conversion between ΔGca output from step S1081 and Kde output from step S1082, calculate the absolute error Δ=|ΔGca-ΔGe| and the relative error δ=(Δ / ΔGe)×100%, and generate a consistency evaluation report. S1084: Based on the error analysis results of step S1083, establish a nonlinear error compensation model ΔGco=ΔGca+α×(ΔGe-ΔGca)^2+β×ln(Kde); Where α and β are compensation coefficients optimized through machine learning. The model is trained using historical data from steps S1081-S1083 to generate the compensated binding free energy ΔGco. S1085: Perform cross-validation loop, compare the ΔGco generated in step S1084 with the Kde measured in step S1082 through thermodynamic conversion again, and calculate the compensating consistency index.

[0013] This application also provides a drug information processing system based on the hydrophobic pocket of the HPV L2 transmembrane region, the system comprising: The crystal structure acquisition module is used to acquire the high-resolution structure of the HPV L2 transmembrane region by X-ray crystallography. The membrane environment reconstruction module is used to reconstruct the physiological environment of HPV L2 transmembrane proteins corresponding to high-resolution structures in a lipid bilayer using structural reconstruction tools. The structure optimization module is used to optimize the structure of HPV L2 transmembrane protein after reconstructing the physiological environment through molecular docking and molecular dynamics simulation; The simulation parameter setting module is used to set molecular dynamics simulation parameters, including using specific molecular dynamics software and controlling the temperature, pressure and simulation step size to preset parameters under canonical isobaric and isothermal ensemble conditions. The molecular dynamics sampling module runs molecular dynamics simulations based on preset parameters and records trajectory data during the simulation process. The pocket dynamic analysis module is used to analyze the recorded trajectory data, extract the dynamic change characteristics of the hydrophobic pocket in the HPV L2 transmembrane region, and use the pocket volume calculation tool to calculate the volume change information of the hydrophobic pocket. The critical residue identification module is used to identify critical residues in hydrophobic pockets through hydrogen bond analysis and hydrophobicity analysis. The free energy calculation and verification module is used to calculate the binding free energy of small molecules to the HPV L2 transmembrane protein using the molecular mechanics-Poisson-Boltzmann surface area method, determine the binding affinity between the two through surface plasmon resonance experiments, and compare the goodness of agreement between the calculated results and the experimental results.

[0014] According to another embodiment of this application, an electronic device is provided, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the steps of the drug information processing method based on the hydrophobic pocket of the HPV L2 transmembrane region.

[0015] According to another embodiment of this application, a storage medium is also provided, on which a computer program is stored, which, when executed by a processor, implements the steps of the drug information processing method based on the hydrophobic pocket of the HPV L2 transmembrane region.

[0016] As can be seen from the above technical solutions, the present invention has the following advantages: The drug information processing method based on the hydrophobic pocket of the HPV L2 transmembrane region provided by this invention obtains high-resolution structures of 1.8-2.5 Å using X-ray crystallography, clearly revealing the atomic-level details of the HPV L2 transmembrane region. It reconstructs the physiological environment within a lipid bilayer, simulating the lipid composition and ion gradient of the cervical epithelial cell membrane. This avoids problems such as transmembrane helical orientation deviation and hydrophobic pocket closure caused by the lack of lipid hydrophobic interactions in in vitro single protein structures, making the protein conformation closer to the actual in vivo state.

[0017] This invention screens high-affinity binding modes and then eliminates spatial conflicts in static docking through 50ns MD simulation, verifying structural stability. The NPT ensemble and 2fs step size are set to ensure the MD simulation conforms to physiological thermodynamic conditions. A 200ns long-duration simulation covers key dynamic processes of the protein, avoiding the limitation of short-duration simulations in capturing pocket dynamic changes. Volume fluctuations, conformational clusters, and opening widths are extracted to quantify the frequency and amplitude of pocket opening-closing transitions, providing dynamic evidence for drug embedding conditions. Combining hydrogen bond analysis and hydrophobicity analysis, and through triple verification using side chain orientation, ligand distance, and trans-subtype conservation, it is ensured that the identified residues are the direct sites of drug binding. The MM-PBSA method decomposes the binding free energy, quantifying the energy nature of drug-protein binding; the SPR experiment calculates the dissociation constant using a 1:1 Langmuir model, providing experimental affinity data; and the good agreement analysis verifies the reliability of the calculation method. Attached Figure Description

[0018] To more clearly illustrate the technical solution of the present invention, the accompanying drawings used in the description will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0019] Figure 1 The flowchart shows a drug information processing method based on the hydrophobic pocket of the HPV L2 transmembrane region. Figure 2 A schematic diagram of a drug information processing system based on the hydrophobic pocket of the L2 transmembrane region of HPV; Figure 3 This is a schematic diagram of an electronic device. Detailed Implementation

[0020] The following describes in detail the drug information processing method based on the hydrophobic pocket of the HPV L2 transmembrane region according to this application. Specific details, such as particular system structures and techniques, are presented for illustrative purposes and not for limitation, in order to provide a thorough understanding of the embodiments of this application. However, those skilled in the art will understand that this application can also be implemented in other embodiments without these specific details.

[0021] It should be understood that, when used in this specification, terms include indicating the presence of a described feature, integral, step, operation, element, and / or component, but do not exclude the presence or addition of one or more other features, integrals, steps, operations, elements, components, and / or collections thereof. The terms include, encompass, have, and variations thereof mean including but not limited to, unless otherwise specifically emphasized.

[0022] The statements such as "one embodiment" or "some embodiments" described in this application mean that one or more embodiments of this application include the specific features, structures, or characteristics described in that embodiment. Therefore, the statements such as "in one embodiment," "in some embodiments," "in other embodiments," and "in still other embodiments" in this application do not necessarily refer to the same embodiment, but rather mean one or more, but not all, embodiments, unless otherwise specifically emphasized.

[0023] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0024] Please see Figure 1 The diagram shows a flowchart of a drug information processing method based on the hydrophobic pocket of the HPV L2 transmembrane region in a specific embodiment. The method includes: S101: Obtain the high-resolution structure of the HPV L2 transmembrane region by X-ray crystallography.

[0025] In some embodiments, during the structural analysis by X-ray crystallography, recombinant expression and purification of HPV L2 transmembrane protein are performed: Escherichia coli BL21(DE3) strain is selected as the expression host, the gene encoding the HPV L2 transmembrane region (amino acid residues 45-175) is cloned into the pET-28a vector, the concentration of the inducing agent IPTG is set to 0.5 mM, and the bacterial cells are collected after induction at 37°C for 4 hours; the target protein is purified by Ni-NTA affinity chromatography, and then further removed by gel filtration chromatography to ensure that the protein purity is >95%.

[0026] This embodiment also included crystal culture: using the seated drop gas-phase diffusion method, the protein concentration was adjusted to 10 mg / mL, mixed with crystallization buffer at a 1:1 volume ratio, and incubated at 20°C for 14 days to form crystals; diffraction data were collected using a synchrotron radiation source, with the diffraction resolution set to 1.8-2.5 Å; data processing was performed using HKL-2000 software, and a complete diffraction dataset was obtained through integration, scaling, and merging; the molecular substitution method was used for structural analysis, and the final structure had an R working factor < 0.25 and an R free factor < 0.30. The dihedral angles of the protein backbone were verified using PROCHECK software to ensure structural reliability, and the structure was then saved.

[0027] It should be noted that protein crystals are composed of a large number of regularly arranged protein molecules. X-rays interact with the electrons of atoms in the crystal to produce diffraction. By recording the intensity and position of the diffraction points, the electron density map is reconstructed using Fourier transform. Then, combined with the amino acid sequence, an atomic model is constructed, and finally a high-resolution three-dimensional structure is obtained.

[0028] Step S102: Use a structural reconstruction tool to reconstruct the physiological environment of the HPV L2 transmembrane protein corresponding to the high-resolution structure described in step S101 in the lipid bilayer.

[0029] In some embodiments, the HPV L2 transmembrane protein is naturally present in the lipid bilayer of the host cell membrane in vivo. Its conformation, the openness of its hydrophobic pockets, and the exposure of key residues are highly dependent on the lipid environment and solvent conditions. The structural reconstruction tool recreates the physiological microenvironment of the protein by simulating the lipid composition, solvent properties, and ion gradient of the cell membrane. The force field parameter (CHARMM36m) describes the interactions (van der Waals, electrostatics, hydrogen bonding) between protein, lipid, and solvent molecules, ensuring that the molecular motion in the reconstruction system conforms to physical laws. Energy minimization eliminates unreasonable repulsion between atoms, bringing the system to a low-energy stable state.

[0030] In this embodiment, the system with minimized energy has no atomic conflicts, ensuring the stability of subsequent molecular dynamics simulations; the generated gro and top files can be directly used for GROMACS simulations, reducing parameter adjustment work in subsequent steps.

[0031] Step S103: Optimize the transmembrane protein structure of HPVL2 after reconstructing the physiological environment in step S102 through molecular docking and molecular dynamics simulation to improve its reliability in drug screening.

[0032] In some embodiments, the molecular docking stage can use the AutoDockVina tool to preprocess the protein structure of S102, remove solvents and ions, retain lipid molecules that form stable hydrophobic interactions with the transmembrane region of the protein, add Gasteiger charges using AutoDockTools, screen 200 HPV inhibitor analogs as ligands, add hydrogen using OpenBabel and calculate AM1-BCC charges to generate three lowest energy conformations.

[0033] In this embodiment, the hydrophobic pocket center is located using PyMOL, a docking grid is set, and after docking is performed, the complexes are sorted by binding energy. The top 10% of complexes are selected, and complexes in which the ligands are not embedded in the pocket or do not interact with the key residues are removed using PyMOL, leaving 30 optimal complexes.

[0034] In the molecular dynamics optimization stage of this embodiment, GROMACS can be used to simulate 30 complexes. The gro / top file of S102 is loaded, the protein is replaced with the docking complex structure, and the force field is CHARMM36m. First, energy minimization is performed, then 10ns NVT equilibration and 10ns NPT equilibration are performed, and finally 50ns production simulation is run.

[0035] In this embodiment, the RMSD of the complex was monitored in the simulation, and complexes with fluctuations >0.3 nm were removed. For stable complex analysis, GROMACS was used to calculate the RMSF of residues, POVME3.0 was used to calculate the pocket volume, and VMD was used to analyze the ligand-protein hydrogen bond occupancy rate. Finally, the protein structure was extracted and superimposed on the original S102 structure. It was confirmed that the positional deviation of key residues was <0.2 nm, and the optimized structure was saved. This structure can ensure conformational stability and intact pocket features, which meets the requirements of structural reliability for drug screening.

[0036] Step S104: Set the molecular dynamics simulation parameters, including using specific molecular dynamics software and controlling the temperature, pressure and simulation step size to preset parameters under canonical isobaric and isothermal ensemble conditions.

[0037] In some embodiments, GROMACS2023.2 is selected as the molecular dynamics software, and simulation parameters are set based on the protein structure optimized by S103.

[0038] A canonical isobaric and isothermal ensemble is used to simultaneously control temperature and pressure. The temperature is maintained at 310 K by a velocity-rescale temperature controller with a temperature control time constant of 0.1 ps. The atomic velocities are adjusted at each step to match the system's kinetic energy with the target temperature. The pressure is maintained at 1 standard atmosphere by a Parrinello-Rahman pressure controller with a pressure control time constant of 2 ps. Isotropic pressure coupling is used, and pressure stabilization is achieved by adjusting the volume of the simulation chamber.

[0039] In this embodiment, the integration step size is set to 2 femtoseconds. Based on the vibration period of protein atoms, the step size is selected as 1 / 5 of the vibration period to ensure the stability of numerical integration and avoid energy non-conservation due to excessively large step size.

[0040] The output is set to save one frame of trajectory data every 2 ps, an energy file every 10 ps, ​​and a checkpoint file every 100 ps, ​​to facilitate restarting after a simulation interruption.

[0041] The reasonable step size and interaction parameters in this embodiment ensure the stability and accuracy of the simulation. The output settings not only fully record dynamic data but also control the amount of data, taking into account the needs of subsequent analysis.

[0042] Step S105: Based on the preset parameters set in step S104, run the molecular dynamics simulation and record the trajectory data during the simulation process.

[0043] In some embodiments, before the simulation starts, the parameter file of S104, the optimized structure file of S103, and the topology file are compiled into an executable tpr file using the gmx grompp tool of GROMACS, while verifying the integrity of the files.

[0044] The preset parameters can be a control temperature of 310 Kelvin, a pressure of 1 standard atmosphere, and a simulation step size of 2 femtoseconds.

[0045] During the simulation, energy data is extracted every 10 ns using the gmxenergy tool to check whether the energy is stable, whether the temperature is maintained at 310±2K, and whether the pressure is maintained at 1±0.2 atm. If any abnormality is found, the simulation is paused, parameters are adjusted, and then restarted. After the simulation ends, the integrity of the trajectory file is verified using gmxcheck to ensure that the total number of frames is 100,000. If the number of missing frames is ≤5, linear interpolation using gmxtrjcat is used to complete the simulation. If the number of missing frames is >5, the simulation for the corresponding time period is rerun.

[0046] As can be seen, the 200 nanosecond simulation duration in this embodiment is sufficient to cover the key dynamic processes of HPV L2 transmembrane proteins, avoiding incomplete capture of dynamic features due to insufficient duration. Real-time monitoring of energy, temperature, and pressure ensures that the simulation process is in thermodynamic equilibrium, preventing the system from deviating from physiological conditions.

[0047] Step S106: Analyze the trajectory data recorded in step S105, extract the dynamic change characteristics of the hydrophobic pocket in the HPV L2 transmembrane region, and use the pocket volume calculation tool to calculate the volume change information of the hydrophobic pocket.

[0048] In some embodiments, trajectory preprocessing can call the gmxtrjconv tool of GROMACS to remove solvents and ions from alignedtrajectory.xtc, retaining only the protein structure; set the -skip10 parameter to extract data every 10 frames, reducing the amount of data while retaining dynamic features; and create a pocket residue index file using gmxmakendx for subsequent targeted analysis.

[0049] The dynamic identification and volume calculation of hydrophobic pockets can be performed using the POVME 3.0 tool. Import the preprocessed trajectory and pocket residue index, and set the parameters: probe radius 1.4 Å, grid resolution 0.2 Å, and search range extending 10 Å outward from the pocket center. Enable dynamic boundary mode, identify pocket regions frame by frame, and determine the internal space of the pocket based on the non-collision principle of the probe spheres. Calculate the volume using a grid-filling method, counting the number of 0.2 ų grids completely inside the pocket. Set a 5-frame sliding window to smooth the volume data, eliminate instantaneous noise, and output a dataset showing the volume change over time, including each frame's time point and corresponding volume value.

[0050] The drug-binding capacity of the HPV L2 transmembrane hydrophobic pocket in this embodiment depends not only on its static structure but also on its volume fluctuations, shape changes, and opening degree. These characteristics are driven by protein conformational motion and reflect the functional state of the pocket under physiological conditions. Trajectory analysis focuses on the target region through preprocessing, POVME3.0 quantifies the pocket volume based on the probe ball principle, and VMD and GROMACS tools extract multi-dimensional features such as conformation, opening, and flexibility to jointly construct a dynamic behavioral map of the pocket.

[0051] Step S107: Identify key residues in the hydrophobic pocket of step S106 by hydrogen bond analysis and hydrophobicity analysis.

[0052] In some embodiments, hydrogen bond analysis uses the Hydrogen Bonds plugin of VMD to analyze the trajectory of pocket residues after S106 preprocessing. The hydrogen bond determination criteria are set as follows: donor-acceptor distance < 3.5 Å and donor-hydrogen-acceptor angle > 120°. The analysis objects include hydrogen bonds inside the pocket residues and hydrogen bonds between the residues and surrounding water / lipid molecules. The hydrogen bond formation is detected frame by frame, and the hydrogen bond donor, acceptor and duration are recorded. The hydrogen bond occupancy rate of each residue is calculated by a custom script, and residues with an occupancy rate > 20% are screened as potential key hydrogen bond residues. A hydrogen bond network snapshot is generated, high-frequency hydrogen bond connections are marked, and transient non-specific hydrogen bonds are excluded.

[0053] Hydrophobicity analysis can be performed using the Hopp-Woods hydrophobicity scale to calculate the hydrophobicity values ​​of pocket residues. For each residue, a value is assigned based on the type of side-chain amino acid, and the average hydrophobicity value of the residue is calculated by weighting by atomic mass. The solvent accessible surface area (SASA) of each pocket residue in GROMACS is calculated using the gmxsasa tool with a probe radius of 1.4 Å. The mean SASA of the residues is statistically analyzed (the smaller the SASA, the stronger the hydrophobicity). A scatter plot of hydrophobicity values ​​versus SASA is plotted, and quadrants are defined: residues with high hydrophobicity values ​​> 2.0 and low SASA are included in the hydrophobic core candidate set, while surface hydrophilic residues with high SASA are excluded.

[0054] The key residue screening and verification method involves performing intersection analysis between hydrogen bond candidate residues and the hydrophobic core candidate set, retaining residues that simultaneously satisfy hydrogen bond occupancy >20% and high hydrophobicity and low SASA; the side chain orientation of residues is verified using PyMOL. The minimum distance between the residue and the optimal ligand for docking with S103 is calculated using GROMACS' gmxmindist, and residues with a distance <4 Å are retained.

[0055] In this embodiment, the binding of the small drug molecule to the hydrophobic pocket of the HPV L2 transmembrane region is based on two key interactions: hydrogen bonding and hydrophobic interaction. Hydrogen bonding analysis screens for residues that stably participate in hydrogen bonding by occupancy rate, excluding transient interactions; hydrophobicity analysis combines hydrophobicity values ​​with solvent exposure to identify residues that constitute the hydrophobic core of the pocket.

[0056] Step S108: The binding free energy of the small molecule to the HPV L2 transmembrane protein was calculated using the molecular mechanics-Poisson-Boltzmann surface area method. The binding affinity between the two was determined by surface plasmon resonance experiments, and the good agreement between the calculated results and the experimental results was compared.

[0057] In some embodiments, binding free energy is an energy index measuring the strength of intermolecular binding. The MM-PBSA method quantifies the binding energy of protein-small molecules from simulated trajectories by decomposing gas-phase interactions, solvation effects, and entropy contributions. Surface plasmon resonance (SPR) is based on the physical phenomenon that changes in surface refractive index cause resonance angle shifts when small molecules bind to proteins on the chip surface, and calculates experimental affinity through kinetic parameters. Both methods describe binding strength from two dimensions: calculated energy and experimental interaction. The consistency analysis verifies the reliability of the calculation methods. High consistency indicates that the calculation results can replace some experiments, reducing drug screening costs; moderate or no consistency suggests the need to optimize calculation parameters to improve prediction accuracy.

[0058] In some specific embodiments, step S108 specifically includes the following methods: Step S1081: Use the gmxmmpbsa tool of GROMACS to perform the binding free energy calculation, set the dielectric constant ε=80 to simulate the polarity of the water environment, and use the LCPO model to approximate the accessible surface area of ​​the nonpolar solvent. Calculate the total energy difference ΔGca between the small molecule and the protein complex, and output the calculation results as a standardized data file for subsequent steps.

[0059] The ΔGca value generated in the steps of this embodiment will be used as the core calculated value for the calculation of experimental consistency evaluation in step S1083.

[0060] Step S1082: SPR experiment was performed using a Biacore 8K instrument. HPV L2 protein was immobilized on the surface of a CM5 sensor chip, and small molecule solutions of varying concentrations were injected. The real-time binding response signal was monitored, and the dissociation constant Kde value was obtained by fitting a 1:1 Langmuir binding model.

[0061] The Kde value generated in this embodiment will serve as the core experimental value for the calculation-experimental consistency evaluation in step S1083, and will also serve as the input parameter for the error compensation model in step S1084.

[0062] Step S1083: Construct the computation-experimental data correlation matrix, perform thermodynamic conversion between ΔGca output in step S1081 and Kde output in step S1082, calculate the absolute error Δ=|ΔGca-ΔGe| and the relative error δ=(Δ / ΔGe)×100%, and generate a consistency evaluation report.

[0063] Step S1084: Based on the error analysis results of step S1083, establish a nonlinear error compensation model ΔGco=ΔGca+α×(ΔGe-ΔGca)^2+β×ln(Kde); Where α and β are compensation coefficients optimized through machine learning. The model is trained using historical data from steps S1081-S1083 to generate the compensated binding free energy ΔGco.

[0064] Step S1085: Execute the cross-validation loop, compare the ΔGco generated in step S1084 with the Kde measured in step S1082 through thermodynamic conversion again, and calculate the compensating consistency index (R²=0.95 or above is considered qualified); if the index is not met, the force field parameter adjustment in step S1081 is triggered, and finally the experimentally verified binding free energy dataset is output.

[0065] In some embodiments, a dynamic correlation between calculated values ​​and experimental values ​​is established: initial calculated values ​​are obtained through the MM-PBSA method, and experimental baseline values ​​are obtained through SPR experiments; then, error analysis is used to locate the source of calculation deviation, and a compensation model containing Kde terms is constructed; finally, cross-validation is used to ensure that the compensated calculated values ​​and experimental values ​​reach a preset degree of agreement.

[0066] As can be seen, the method in this embodiment is suitable for the initial screening of high-throughput drug screening; the SPR experiment directly measures the binding affinity between molecules, and the results have high reliability and can be used as a verification standard for the calculation method; the goodness-of-match analysis clarifies the reliability boundary of the calculation method and avoids drug screening errors caused by calculation deviations; when the goodness-of-match is high, candidate small molecules can be screened first by MM-PBSA and then verified by SPR, which greatly reduces the experimental workload and improves the accuracy of the overall drug screening system.

[0067] In one embodiment of the present invention, based on step S102, the following will provide a possible embodiment and describe its specific implementation in a non-limiting manner. Step S102 specifically includes the following methods: S1021: Based on the amino acid sequence characteristics of HPV L2 transmembrane proteins, matching lipid bilayer components were selected. The lipid components included phosphatidylcholine, sphingomyelin, and cholesterol, and their proportions were determined with reference to the lipidomics data of natural epithelial cell membranes. S1022: The membrane construction module of the structure reconstruction tool generates the initial topology of the lipid bilayer, which includes the atomic coordinates, bonding relationships and stress field parameters of the lipid molecules. S1023: The high-resolution structure of the HPV L2 transmembrane protein obtained in step S101 is translated to the central region of the lipid bilayer. The protein is rotated so that the axis of the transmembrane helix is ​​perpendicular to the lipid bilayer plane. The initial spatial position of the protein and lipid is optimized by minimizing energy through the steepest descent method to reduce the repulsive effect between them. S1024: Add ions with opposite charges to the lipid bilayer system, calculate the system charge density according to the Poisson-Boltzmann equation, and neutralize the overall net charge of the protein and lipid to achieve electroneutrality.

[0068] In step S102 of this embodiment, the hydrophobicity of the transmembrane helix of the HPV L2 transmembrane region protein is matched with the lipid composition of the natural cell membrane to ensure that the lipid environment is consistent with the native state of the protein. Next, the topological structure of the lipid bilayer is generated using a membrane construction tool to determine the atomic type and interaction parameters of each lipid molecule. Then, the protein is extracted from the crystal structure and positioned at the center of the lipid bilayer. The protein orientation is adjusted so that the transmembrane helix is ​​perpendicular to the membrane plane. The relative positions of the protein and lipid are optimized by minimizing energy to avoid structural overlap or excessive tension. The net charge of the system is calculated and ions are added for neutralization to ensure the electroneutrality of the simulated system and avoid abnormal electrostatic interactions.

[0069] In one embodiment of the present invention, based on step S103, the following will provide a possible embodiment and its specific implementation will be described in a non-limiting manner. Step S103 specifically includes the following methods: Step S1031: Prepare a small molecule ligand library, convert the three-dimensional structure of the candidate drug molecule into a simulation-compatible format, and use a force field assignment tool to assign appropriate atom types and charge parameters to it.

[0070] This embodiment employs a graph-based molecular topology generation algorithm, which allocates force field parameters based on atom type and bond connectivity, and calculates partial atomic charges using chemical sensing rules and charge balance formulas, ensuring the accuracy of the physicochemical properties of small molecules.

[0071] Step S1032: Using molecular docking software, the small molecule ligand is docked with the reconstructed HPV L2 transmembrane protein structure, and a docking method based on a fast search algorithm is used to generate multiple possible binding conformations.

[0072] It should be noted that this embodiment uses a fast Fourier transform correlation algorithm for spatial matching. The mathematical model is based on shape complementarity and chemical scoring function. It maximizes the complementarity between ligands and pockets through rotation and translation search, and uses empirical methods combined with free energy formulas for conformation ranking.

[0073] Step S1033: Perform molecular dynamics simulation preprocessing on the docking conformation with the highest score, including energy minimization and short-time equilibrium simulation in an explicit solvent environment to eliminate local atomic conflicts and stabilize the complex structure.

[0074] It should be noted that numerical integration is performed using Newton's equations of motion, a restricted energy minimization algorithm is employed, and the mathematical model is based on the principle of potential energy surface gradient descent. The total energy of the system is minimized by iteratively adjusting the atomic positions, and harmonic potential constraints are used to maintain the protein's main structure.

[0075] Step S1034: By analyzing the simulated trajectory, the stability of the combined conformation is evaluated. The root mean square deviation and root mean square fluctuation are used as discrimination indicators to screen out the optimized conformation with the lowest structural deviation for subsequent drug screening.

[0076] It should be noted that structural bias is calculated based on time series analysis algorithms, and the root mean square bias formula is used to quantify structural changes. Conformation stability is assessed by comparing changes in atomic coordinates over time, and a clustering algorithm is used to identify representative structures. The most stable binding conformations are then selected to form an optimized structure set.

[0077] This embodiment optimizes the protein-ligand complex structure through computational simulation to more closely approximate the actual biological binding state. By combining the spatial search capabilities of molecular docking with the relaxation effect of molecular dynamics simulation, it can be seen that by rapidly sampling and identifying potential binding sites, and then using physical force fields to refine the structure, the accuracy of structure optimization and the reliability of drug screening can be improved.

[0078] In one embodiment of the present invention, based on step S105, the following will provide a possible embodiment and describe its specific implementation in a non-limiting manner. Step S105 specifically includes the following: S1051: Call the GROMACS verification step S104 associated with the gro structure file, topology file, and mdp parameter file to check the atom type matching and the rationality of bond lengths and bond angles; expand the system box by 0.5nm using the gmx editconf tool and set the periodic boundary conditions to cubic boxes; initialize the atom velocities using the Maxwell-Boltzmann distribution and generate the initial simulation input file using gmxgrompp.

[0079] It should be noted that the type of each atom is matched with the force field parameter library using a hash table. If the atom type does not match in the library, it is marked as an error and the initialization is terminated to avoid calculation errors caused by missing force field parameters.

[0080] In this embodiment, velocity initialization ensures that the initial energy of the system conforms to thermodynamic conditions, periodic boundary conditions eliminate the boundary effects of the finite simulation box, and file verification ensures the integrity and consistency of the input data, providing a reliable initial system for subsequent simulations.

[0081] S1052: Perform secondary energy minimization: First, use the steepest descent method to iterate for 2000 steps to reduce the maximum atomic force of the system to below 5 kJ / (mol·nm); then switch to the conjugate gradient method to iterate for 1000 steps to reduce the maximum force to below 2 kJ / (mol·nm); Use the gmxenergy tool to extract the potential energy curve after energy minimization. If the curve does not converge, increase the number of iterations to 5000 and recalculate. After completion, run a 100ps short-range simulation and calculate the protein backbone RMSD using gmxrms. If the RMSD increases by more than 0.2nm within the initial 50ps, the system is considered unstable and returns to S1051 to re-verify the file.

[0082] It should be noted that minimizing the secondary energy gradually reduces the system's potential energy and eliminates unreasonable spatial conflicts between atoms; short-range simulation and RMSD monitoring verify whether the system can remain stable under dynamic conditions, avoiding the waste of resources caused by directly entering long simulations.

[0083] S1053: Staged processing, in which the first stage: 5ns NVT equilibration is performed, and a velocity-rescale temperature controller is used to maintain 310K. Positional constraints are imposed on the protein backbone atoms, allowing lipids and solvents to move freely. One frame of trajectory is output every 1ps, and temperature fluctuations are monitored by gmxtemp to ensure that the fluctuation range is within 310±2K.

[0084] Phase 2: Perform 10ns NPT balancing, retain main chain constraints, enable Parrinello-Rahman pressure controller to maintain 1 atm, output one frame trajectory every 1ps, monitor pressure fluctuations using gmxpressure to ensure fluctuation range is within 1±0.2 atm; after balancing, use gmxdensity to calculate system density and verify whether it is stable around 0.997 g / cm³.

[0085] It should be noted that the phased equilibration first fixes the protein backbone, allowing lipids and solvents to adapt to the environment (NVT phase), and then adjusts the box volume to physiological density (NPT phase). Temperature and pressure control work together to maintain physiological conditions, while positional constraints prevent drastic conformational changes in the protein during the initial equilibration phase.

[0086] S1054: Remove all atomic position constraints and use the leap-frog integrator to solve Newton's equations of motion; set the neighbor search algorithm to Verlet list, set the cutoff distance for short-range van der Waals and Coulomb interactions to 1.0 nm, and use the ParticleMeshEwald (PME) method for long-range Coulomb interactions; save a trajectory frame every 2 ps, an energy file every 10 ps, ​​and a checkpoint file every 100 ps.

[0087] It should be noted that the leap-frog integrator efficiently solves atomic motion trajectories, the PME method accurately calculates long-range electrostatic interactions, and the Verlet list reduces redundant calculations, ensuring that trajectory data covers important dynamic changes in protein conformation.

[0088] S1055: After the simulation, use gmxcheck to verify the integrity of the trajectory file, check if the total number of frames is 100,000, and use gmxtrjcat to complete the missing frames if the number of missing frames is ≤5. If the number of missing frames is >5, rerun the simulation for the corresponding time period. Calculate the RMSD of the entire trajectory protein backbone using gmxrms to ensure that the overall fluctuation range is within 0.1-0.3nm. Name the trajectory file, energy file, and checkpoint file according to the simulation date-protein ID-200ns and generate a simulation report.

[0089] It should be noted that trajectory completion ensures data continuity, checksum guarantees data integrity, RMSD trend analysis verifies the reliability of simulation results, and the final stored high-quality trajectory data provides complete and accurate raw material for subsequent dynamic feature extraction of hydrophobic pockets.

[0090] As can be seen, step S105 calculates the position, velocity, and force changes of each atom in continuous time by solving Newton's equation of motion for each atom, thereby simulating the dynamic movement of HPV L2 transmembrane protein under physiological conditions. The duration of 200 ns is sufficient to cover the opening and closing of the hydrophobic pocket of the protein conformation and the side chain swinging of key residues. The trajectory data records these dynamic characteristics, providing information in the time dimension for pocket volume changes and residue interactions. The phased balancing and strict parameter control ensure that the simulation system gradually transitions from the initial state to a stable physiological state.

[0091] In one embodiment of the present invention, based on step S106, the following will provide a possible embodiment and describe its specific implementation in a non-limiting manner. Step S106 specifically includes the following methods: S1061: The trajectory data recorded in step S105 is preprocessed. The trjconv tool of GROMACS is used to remove water molecules, ions and irrelevant ligand molecules in the system, and only the atomic coordinates of the HPV L2 transmembrane protein are retained. At the same time, the trajectory is fitted with rmsd with the protein structure optimized in step S103 as a reference to correct the coordinate offset caused by the overall translation and rotation of the protein.

[0092] S1062: Hydrophobic pocket boundaries were defined based on the side chain atoms of hydrophobic residues in the HPV L2 transmembrane region protein. The side chain heavy atoms of leucine, isoleucine, valine, phenylalanine, tryptophan, and methionine were selected as pocket marker atoms, and the spatial coordinates of these marker atoms in each frame trajectory were extracted using the MDAnalysis tool.

[0093] S1063: Perform dynamic feature analysis on the trajectory of the pocket-marked atoms, calculate the root mean square fluctuation of the pocket in each frame to quantify the flexibility, and use the formula RMSFi=√[(1 / T)*Σ(t=1 to T)(xi(t)-xiavg)²+(yi(t)-yiavg)²+(zi(t)-ziavg)²], where T is the number of simulated frames and xi(t) is the coordinate of the i-th marked atom at time t; use principal component analysis to extract the principal components of the pocket-marked atom coordinates and identify the main motion patterns of the pocket.

[0094] S1064: The volume of the hydrophobic pocket in each frame trajectory is calculated using the POVME3.0 tool. Based on the convex hull algorithm of pocket marker atoms, the marker atoms are regarded as vertices of the convex hull. The volume of the space enclosed by the convex hull is calculated by the formula V=(1 / 6)*Σ(|(ai×aj)·ak|), where ai, aj, and ak are the three vertex vectors of the convex hull tetrahedron.

[0095] S1065: Verify the volume calculation results by selecting 2-3 independent 200-nanosecond MD simulation copies and comparing the time-series correlation of pocket volumes in each copy; mutate key residues in the pocket, rerun the 50-nanosecond MD simulation, and observe whether the pocket volume change meets expectations to ensure the reliability of the results.

[0096] It can be seen that the dynamic characteristics and volume changes of the hydrophobic pocket in step S106 are key to drug binding; the drug needs to adapt to different pocket conformations to bind effectively. Trajectory preprocessing removes irrelevant molecules, focusing on the structural changes of the protein itself; the pocket boundary is defined based on hydrophobic residues because the formation of the hydrophobic pocket is based on the spatial arrangement of these residues; RMSF and PCA quantify the pocket's flexibility and main motion modes, giving the pocket size at different time points; result validation ensures that these characteristics are not accidental results of simulation errors.

[0097] In conjunction with step S106, irrelevant molecules in the trajectory can be cleaned up and the coordinates corrected; pockets can be labeled with hydrophobic residues; dynamic features can be extracted using statistical methods; the volume can be calculated using geometric algorithms; and the results can be verified through multi-copy comparison and mutant simulation.

[0098] In one embodiment of the present invention, based on step S107, the following will provide a possible embodiment and its specific implementation will be described in a non-limiting manner. Step S107 specifically includes the following methods: S1071: Load the list of hydrophobic pocket boundary residues determined in step S1062 using PyMOL and export the residue numbers; use the gmxtrjconv tool of GROMACS to extract the atomic coordinates of the residues from the preprocessed trajectory to generate pocket residue-specific trajectories; also extract the topological information of the corresponding residues and save it as an index file; verify the integrity of the residue set using gmxmakendx to ensure that no non-pocket residues are omitted or incorrectly included.

[0099] This embodiment can extract the coordinate data of the corresponding residues from the entire trajectory based on the residue index file, maintain the integrity of the time series through coordinate mapping, and remove irrelevant atoms to reduce the amount of data.

[0100] In this embodiment, the integrity verification uses set operations to compare manually labeled residues with extracted residues. The coverage rate is calculated as (number of extracted residues / number of labeled residues) × 100%, which must be ≥95% to ensure that no critical residues are missed.

[0101] Thus, accurately defining the pocket residue set is the foundation for subsequent analysis, trajectory extraction focuses on the target region and reduces interference from irrelevant data, while the verification step ensures the integrity and accuracy of the analyzed object.

[0102] S1072: The Hydrogen Bonds plugin of VMD was used to analyze the trajectory of pocket residues: the hydrogen bond determination criteria were set as donor-acceptor distance < 3.5 Å and donor-hydrogen-acceptor angle > 120°. The analysis objects included hydrogen bonds inside the pocket residues and hydrogen bonds between the residues and surrounding water molecules and lipid molecules; the hydrogen bond formation was detected frame by frame, and the hydrogen bond donor residues, acceptor residues and duration were recorded.

[0103] The hydrogen bond occupancy rate of each residue is calculated using a custom script, and residues with an occupancy rate >20% are selected as potential key hydrogen bond residues; a dynamic hydrogen bond network graph is generated, and frequently occurring hydrogen bond connections are marked.

[0104] The hydrogen bond determination criteria in this embodiment ensure the authenticity of the interaction, the occupancy rate quantifies the stability of the residues participating in the hydrogen bond, and the network analysis reveals the overall position of the residues in the hydrogen bond interaction, providing a basis for identifying specifically binding residues.

[0105] S1073: The average hydrophobicity value of pocket residues was calculated using the Hopp-Woods hydrophobicity scale: For each residue, a value was assigned according to the type of its side chain amino acid, and the weighted average of all atoms was taken as the residue hydrophobicity value; the solvent-accessible surface area of ​​each pocket residue was calculated using the gmx sasa tool of GROMACS, with the probe radius set to 1.4 Å, and the average SASA value of the residues was calculated; a scatter plot of hydrophobicity value and SASA was plotted, and quadrants were divided: residues with high hydrophobicity value and low SASA were included in the hydrophobic core candidate set.

[0106] In this embodiment, the hydrophobicity value quantifies the chemical properties of the residues, and the SASA reflects their spatial exposure state. The combination of the two can accurately identify the residues that constitute the core of the hydrophobic pocket, which are the main binding sites of the drug's hydrophobic groups.

[0107] S1074: Integrate the hydrogen bond occupancy rate of S1072, the hydrophobicity data of S1073, and the dynamic characteristics of S1064 to calculate the hydrogen bond occupancy rate and the Pearson coefficient of the residue RMSF, and screen hydrogen bond residues with an occupancy rate > 20% and RMSF < 0.2nm. The mutual information (MI) between hydrophobicity value and volume change is calculated, and hydrophobic core residues with MI > 0.4 are screened. Intersection analysis is performed on the two types of candidate residues, and residues that simultaneously satisfy hydrogen bonding and hydrophobicity conditions are retained to form a preliminary list of key residues.

[0108] This embodiment combines static interaction characteristics with dynamic stability in its association analysis, excluding residues that only participate in the action transiently or are dynamically unstable. The initial screening set can better reflect the sustained contribution to drug binding.

[0109] S1075: Structural verification of the initial screening residues using PyMOL: checking whether the side chains of the residues face the inside of the pocket; The minimum distance between the residue and the optimal ligand for docking in step S103 was calculated using gmxmindist from GROMACS, and residues with a distance < 4 Å were retained. Non-conservative residues were removed by referring to the conserved residue database of the homologous HPV subtype L2 protein. The list of key residues was classified according to their function type and included quantitative parameters such as residue number, average hydrophobicity value, and hydrogen bond occupancy rate.

[0110] It should be noted that structural validation ensures that the residues can interact with the drug in space, ligand distance screening focuses on directly acting residues, conservation analysis ensures the broad applicability of the residues, and multi-dimensional validation ensures that the key residues finally identified have functional necessity and reliability.

[0111] In step S107 of this embodiment, the binding of the hydrophobic pocket of the HPV L2 transmembrane region to the drug molecule is based on two types of key interactions: hydrogen bonding and hydrophobic interaction. Key residues are residues that play a dominant role in these interactions and remain stable during dynamic changes. By analyzing the formation mode of hydrogen bonds and the hydrophobic characteristics of residues in the trajectory, combined with the dynamic stability of residues, residues that are crucial to drug binding can be identified. The identification of these residues provides target sites for subsequent drug design, ensuring that the designed small molecule can accurately match the interaction characteristics of the pocket.

[0112] It should be noted that the Cytoscape network analysis tool treats each frame of hydrogen bond network as a residue-node hydrogen bond edge structure, calculates the degree and betweenness centrality of nodes, and filters residues with a degree > 5 and betweenness centrality > 0.3 to identify residues that play a pivotal role in the hydrogen bond network, rather than relying solely on static occupancy.

[0113] Temperature factor-corrected hydrophobicity analysis combines the B factor of X-ray crystallography to correct the calculated hydrophobicity value and reduce the weight of highly flexible residues.

[0114] In this embodiment, virtual mutations were performed on the initial screening residues, and a 10 ns molecular dynamics simulation was conducted on the mutants and the optimal ligands of step S103 to calculate the binding energy change (ΔΔG) before and after the mutation. Residues with ΔΔG > 2 kJ / mol were retained, and the reliability of identification was enhanced through functional verification.

[0115] It should be understood that the sequence number of each step in the above embodiments does not imply the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of the present invention.

[0116] The following are embodiments of the drug information processing system based on the hydrophobic pocket of the HPV L2 transmembrane region provided in this disclosure. This system and the drug information processing method based on the hydrophobic pocket of the HPV L2 transmembrane region in the above embodiments belong to the same inventive concept. For details not described in detail in the embodiments of the drug information processing system based on the hydrophobic pocket of the HPV L2 transmembrane region, please refer to the embodiments of the drug information processing method based on the HPV L2 transmembrane region.

[0117] like Figure 2 As shown, the system includes: Crystal structure acquisition module 201 is used to acquire the high-resolution structure of the HPV L2 transmembrane region by X-ray crystallography. Membrane environment reconstruction module 202 is used to reconstruct the physiological environment of HPV L2 transmembrane proteins corresponding to high-resolution structures in a lipid bilayer using structural reconstruction tools. The structure optimization module 203 is used to optimize the structure of the HPVL2 transmembrane protein after reconstructing the physiological environment through molecular docking and molecular dynamics simulation. The simulation parameter setting module 204 is used to set molecular dynamics simulation parameters, including using specific molecular dynamics software and controlling the temperature, pressure and simulation step size to preset parameters under canonical isobaric and isothermal ensemble conditions. The molecular dynamics sampling module 205 runs a molecular dynamics simulation based on preset parameters and records trajectory data during the simulation process. The pocket dynamic analysis module 206 is used to analyze the recorded trajectory data, extract the dynamic change characteristics of the hydrophobic pocket in the HPV L2 transmembrane region, and use the pocket volume calculation tool to calculate the volume change information of the hydrophobic pocket. The critical residue identification module 207 is used to identify critical residues in the hydrophobic pocket through hydrogen bond analysis and hydrophobicity analysis. The free energy calculation and verification module 208 is used to calculate the binding free energy of the small molecule to the HPV L2 transmembrane protein using the molecular mechanics-Poisson-Boltzmann surface area method, determine the binding affinity between the two through surface plasmon resonance experiments, and compare the goodness of agreement between the calculated results and the experimental results.

[0118] like Figure 3 As shown, this application also provides an electronic device, including a display module 103, a memory 102, a processor 101, and a computer program stored in the memory and executable on the processor 101. When the processor 101 executes the program, it implements the steps of a drug information processing method based on hydrophobic pockets in the transmembrane region of HPV L2.

[0119] In embodiments of the present invention, electronic devices include, but are not limited to, laptop computers, desktop computers, workstations, personal digital assistants, servers, blade servers, mainframe computers, and other suitable computers. Electronic devices may also represent various forms of mobile devices, such as personal digital processors, cellular phones, smartphones, wearable devices, and other similar computing devices. The components shown herein, their connections and relationships, and their functions are merely examples and are not intended to limit the implementation of the embodiments described and / or claimed herein.

[0120] In this embodiment, processor 101 may be implemented using at least one of an application-specific integrated circuit, a programmable logic device, a field-programmable gate array, a processor, a controller, a microcontroller, a microprocessor, or an electronic unit designed to perform the functions described herein. In some cases, such an implementation may be implemented within a controller. For software implementation, implementations such as processes or functions may be implemented with separate software modules that allow the performance of at least one function or operation. Software code may be implemented by a software application (or program) written in any suitable programming language, and the software code may be stored in memory and executed by the controller.

[0121] The display module 103 is used to display information input by the user or information provided to the user. The display module 103 may include a display panel, which may be configured in the form of a liquid crystal display, an organic light-emitting diode, or the like.

[0122] The memory 102 can be used to store software programs and various data. The memory 102 may include high-speed random access memory, and may also include non-volatile memory, such as at least one disk storage device, flash memory device, or other volatile solid-state storage device.

[0123] This application also provides a storage medium storing a computer program thereon, which, when executed by a processor, implements the steps of a drug information processing method based on hydrophobic pockets in the HPV L2 transmembrane region.

[0124] The storage medium may be any combination of one or more readable media. A readable medium may be a readable signal medium or a readable storage medium. A readable storage medium may be, for example,, but not limited to, an electrical, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any combination thereof. More specific examples (a non-exhaustive list) of readable storage media include: electrical connections having one or more wires, portable disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), optical fiber, portable compact disk read-only memory (CD-ROM), optical storage devices, magnetic storage devices, or any suitable combination thereof.

[0125] In a storage medium, a readable signal medium may include a data signal propagated in baseband or as part of a carrier wave, carrying readable program code. This propagated data signal may take various forms, including but not limited to electromagnetic signals, optical signals, or any suitable combination thereof. A readable signal medium may also be any readable medium other than a readable storage medium, capable of sending, propagating, or transmitting a program for use by or in conjunction with an instruction execution system, apparatus, or device.

[0126] The above description of the disclosed embodiments enables those skilled in the art to make or use the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A drug information processing method based on the hydrophobic pocket of the HPV L2 transmembrane region, characterized in that the method... include: S101: Obtain the high-resolution structure of the HPV L2 transmembrane region by X-ray crystallography; S102: Using a structural reconstruction tool, reconstruct the physiological environment of the HPV L2 transmembrane protein corresponding to the high-resolution structure described in step S101 in the lipid bilayer; S103: Optimize the structure of HPV L2 transmembrane protein after reconstructing the physiological environment in step S102 through molecular docking and molecular dynamics simulation; S104: Set the molecular dynamics simulation parameters, including using specific molecular dynamics software and controlling the temperature, pressure and simulation step size to preset parameters under canonical isobaric and isothermal ensemble conditions; S105: Based on the preset parameters set in step S104, run the molecular dynamics simulation and record the trajectory data during the simulation process; S106: Analyze the trajectory data recorded in step S105, extract the dynamic change characteristics of the hydrophobic pocket in the HPV L2 transmembrane region, and use the pocket volume calculation tool to calculate the volume change information of the hydrophobic pocket. S107: Key residues in the hydrophobic pocket of step S106 were identified by hydrogen bond analysis and hydrophobicity analysis. S108: The binding free energy of the small molecule to the HPV L2 transmembrane protein was calculated using the molecular mechanics-Poisson-Boltzmann surface area method. The binding affinity between the two was determined by surface plasmon resonance experiments, and the good agreement between the calculated and experimental results was compared.

2. The drug information processing method based on the hydrophobic pocket of the HPV L2 transmembrane region according to claim 1, characterized in that, Step S102 specifically includes the following methods: Based on the amino acid sequence characteristics of HPV L2 transmembrane proteins, matching lipid bilayer components were selected; The membrane construction module of the structure reconstruction tool generates the initial topology of the lipid bilayer; The high-resolution structure of the HPV L2 transmembrane protein obtained in step S101 is translated to the central region of the lipid bilayer. The protein is rotated so that the axis of the transmembrane helix is ​​perpendicular to the lipid bilayer plane. The initial spatial position of the protein and lipid is optimized by minimizing energy using the steepest descent method. Adding ions with opposite charges to the lipid bilayer system and calculating the system charge density according to the Poisson-Boltzmann equation neutralizes the overall net charge of the protein and lipids to achieve electroneutrality.

3. The drug information processing method based on the hydrophobic pocket of the HPV L2 transmembrane region according to claim 1, characterized in that, Step S103 specifically includes the following methods: Step S1031: Assign force field parameters to candidate small molecules and calculate partial atomic charges using a graph theory topology generation algorithm to form a small molecule ligand library compatible with the force field of HPV L2 transmembrane region proteins; Step S1032: Using the ligand library obtained in S1031 as the object, the Fast Fourier Transform correlation algorithm is used to search for shape complementarity and chemical scores between each small molecule and HPV L2 transmembrane protein, generating and sorting multiple binding conformations for stability testing. Step S1033: Perform restricted energy minimization and short-time equilibrium on the highest-scoring binding conformation in S1032 in an explicit solvent environment based on potential energy surface gradient descent; Step S1034: Using the trajectory of the stabilized complex after S1033, calculate the root mean square deviation and root mean square fluctuation, and select the binding conformation with the lowest structural deviation as the optimal representative.

4. The drug information processing method based on the hydrophobic pocket of the HPV L2 transmembrane region according to claim 1, characterized in that, Step S105 specifically includes the following methods: Step S1051: Verify the gro, top, and mdp files using GROMACS, expand the cubic periodic cell, assign Maxwell-Boltzmann initial velocities, and generate the initial simulation input file; Step S1052: The structure obtained in step S1051 is subjected to the steepest descent and conjugate gradient second-order energy minimization in sequence to reduce the maximum atomic force to below 2kJ / (mol·nm). The system stability is confirmed by the 100ps short-range simulated main chain RMSD≤0.2nm. Otherwise, it is back to S1051 for re-verification. Step S1053: Based on the stable structure in step S1052, first perform 5ns NVT equilibration and then 10ns NPT equilibration, with temperature control at 310±2K, pressure control at 1±0.2atm and main chain position restriction retained, and complete the equilibration after the density reaches about 0.997g / cm³. Step S1054: Remove all position restrictions from step S1053, and use the leap-frog integrator, Verlet neighbor list, 1.0nm intercept and PME long-range processing to write the trajectory every 2ps, the energy every 10ps, and the checkpoint every 100ps, thereby continuously acquiring the production trajectory for 200ns with a step size of 2fs. Step S1055: Verify the frame integrity of the trajectory obtained in step S1054 using gmxcheck. Automatically complete if there are ≤5 missing frames, and rerun the corresponding time period if there are >5 missing frames. Ensure that the main chain RMSD fluctuation of the entire trajectory is within 0.1–0.3nm. Finally, name and archive the trajectory, energy, and checkpoint files and generate a simulation report to complete data storage and quality confirmation.

5. The drug information processing method based on the hydrophobic pocket of the HPV L2 transmembrane region according to claim 1, characterized in that, Step S106 specifically includes the following methods: Step S1061: The trajectory from step S105 was dehydrated, deionized, and fitted with rmsd using the GROMACStrjconv tool, retaining only the protein atom coordinates of the HPV L2 transmembrane region to obtain the corrected pure protein trajectory. Step S1062: Based on the purified protein trajectory in step S1061, MDAnalysis is used to extract the side chain heavy atoms of leucine, isoleucine, valine, phenylalanine, tryptophan, and methionine as pocket marker atoms to form a set of hydrophobic pocket boundary atoms. Step S1063: Calculate the root mean square fluctuation of the pocket-marked atom set in step S1062 and perform principal component analysis to obtain the pocket flexibility quantization value and main motion modes; Step S1064: Based on the pocket-marked atom set in step S1062, the volume of the hydrophobic pocket in each frame is calculated using the POVME3.0 convex hull algorithm, and the volume time series is output to realize continuous monitoring of the pocket opening and closing degree; Step S1065: Using the volume time series from step S1064 as the object, perform independent copy correlation test and key residue mutation resimulation to verify the reproducibility and causality of the volume calculation and ensure the reliability of the dynamic volume results of the hydrophobic pocket.

6. The drug information processing method based on the hydrophobic pocket of the HPV L2 transmembrane region according to claim 1, characterized in that, Step S107 specifically includes the following methods: Step S1071: Using PyMOL and GROMACS together, export the numbers of the hydrophobic pocket boundary residues and generate a dedicated trajectory and index file. After verifying the coverage rate by set operation, the pocket residue set required for subsequent analysis is locked. Step S1072: Using the pocket residue set obtained in S1071 as the object, hydrogen bonds are detected frame by frame using the VMD plugin according to the standard that the donor-acceptor distance is less than 3.5 Å and the included angle is greater than 120°. The hydrogen bond occupancy rate is calculated and a network diagram is drawn. Residues with an occupancy rate greater than 20% are screened to form potential hydrogen bond key residues. Step S1073: Based on the pocket residue set in S1071, Hopp-Woods scaling is used for assignment and solvent-accessible surface area calculation is combined with probe radius of 1.4 Å. Residues with high hydrophobicity values ​​and low surface area are included in the hydrophobic core candidate set, which together with hydrogen bond residues constitute a dual-attribute candidate pool. Step S1074: Correlation analysis is performed on the hydrogen bond occupancy rate of S1072 and the hydrophobicity value of S1073 with the RMSF and volume change of S1064. Intersecting residues with hydrogen bond occupancy rate greater than 20%, RMSF less than 0.2nm, and mutual information greater than 0.4 are retained to form a preliminary screening list of key residues. Step S1075: For the residues in the initial screening list of S1074, the key residue list is finally determined and classified and output by PyMOL orientation check, GROMACS minimum distance less than 4Å screening and homologous subtype conservation comparison, with hydrophobicity value and hydrogen bond occupancy rate quantification parameters.

7. The drug information processing method based on the hydrophobic pocket of the HPV L2 transmembrane region according to claim 1, characterized in that, Step S108 specifically includes the following methods: S1081: The binding free energy was calculated using the gmxmmpbsa tool of GROMACS. The dielectric constant ε=80 was set to simulate the polarity of the water environment. The accessible surface area of ​​the nonpolar solvent was approximated by the LCPO model. The total energy difference ΔGca between the small molecule and the protein complex was calculated. S1082: SPR experiments were performed using a Biacore 8K instrument. HPV L2 protein was immobilized on the surface of a CM5 sensor chip, and small molecule solutions of varying concentrations were injected. The real-time binding response signal was monitored, and the dissociation constant Kde was obtained by fitting a 1:1 Langmuir binding model. S1083: Construct the computation-experimental data correlation matrix, perform thermodynamic conversion between ΔGca output from step S1081 and Kde output from step S1082, calculate the absolute error Δ=|ΔGca-ΔGe| and the relative error δ=(Δ / ΔGe)×100%, and generate a consistency evaluation report. S1084: Based on the error analysis results of step S1083, establish a nonlinear error compensation model ΔGco=ΔGca+α×(ΔGe-ΔGca)^2+β×ln(Kde); Where α and β are compensation coefficients optimized through machine learning. The model is trained using historical data from steps S1081-S1083 to generate the compensated binding free energy ΔGco. S1085: Perform cross-validation loop, compare the ΔGco generated in step S1084 with the Kde measured in step S1082 through thermodynamic conversion again, and calculate the compensating consistency index.

8. A drug information processing system based on the hydrophobic pocket of the HPV L2 transmembrane region, characterized in that, The system is used to implement the drug information processing method based on the hydrophobic pocket of the HPV L2 transmembrane region as described in any one of claims 1 to 7; the system includes: The crystal structure acquisition module is used to acquire the high-resolution structure of the HPV L2 transmembrane region by X-ray crystallography. The membrane environment reconstruction module is used to reconstruct the physiological environment of HPVL2 transmembrane proteins corresponding to high-resolution structures in lipid bilayers using structural reconstruction tools. The structure optimization module is used to optimize the structure of HPV L2 transmembrane protein after reconstructing the physiological environment through molecular docking and molecular dynamics simulation; The simulation parameter setting module is used to set molecular dynamics simulation parameters, including using specific molecular dynamics software and controlling the temperature, pressure and simulation step size to preset parameters under canonical isobaric and isothermal ensemble conditions. The molecular dynamics sampling module runs molecular dynamics simulations based on preset parameters and records trajectory data during the simulation process. The pocket dynamic analysis module is used to analyze the recorded trajectory data, extract the dynamic change characteristics of the hydrophobic pocket in the HPV L2 transmembrane region, and use the pocket volume calculation tool to calculate the volume change information of the hydrophobic pocket. The critical residue identification module is used to identify critical residues in hydrophobic pockets through hydrogen bond analysis and hydrophobicity analysis. The free energy calculation and verification module is used to calculate the binding free energy of small molecules to the HPV L2 transmembrane protein using the molecular mechanics-Poisson-Boltzmann surface area method, determine the binding affinity between the two through surface plasmon resonance experiments, and compare the goodness of agreement between the calculated results and the experimental results.

9. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the steps of the drug information processing method based on the hydrophobic pocket of the HPV L2 transmembrane region as described in any one of claims 1 to 7.

10. A storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the steps of the drug information processing method based on the hydrophobic pocket of the HPV L2 transmembrane region as described in any one of claims 1 to 7.

Images