Interactions between Inhibitors and 5-Lipoxygenase: Insights from Gaussian Accelerated Molecular Dynamics and Markov State Models

Inflammation is a protective stress response triggered by external stimuli, with 5-lipoxygenase (5LOX) playing a pivotal role as a potent mediator of the leukotriene (Lts) inflammatory pathway. Nordihydroguaiaretic acid (NDGA) functions as a natural orthosteric inhibitor of 5LOX, while 3-acetyl-11-keto-β-boswellic acid (AKBA) acts as a natural allosteric inhibitor targeting 5LOX. However, the precise mechanisms of inhibition have remained unclear. In this study, Gaussian accelerated molecular dynamics (GaMD) simulation was employed to elucidate the inhibitory mechanisms of NDGA and AKBA on 5LOX. It was found that the orthosteric inhibitor NDGA was tightly bound in the protein’s active pocket, occupying the active site and inhibiting the catalytic activity of the 5LOX enzyme through competitive inhibition. The binding of the allosteric inhibitor AKBA induced significant changes at the distal active site, leading to a conformational shift of residues 168–173 from a loop to an α-helix and significant negative correlated motions between residues 285–290 and 375–400, reducing the distance between these segments. In the simulation, the volume of the active cavity in the stable conformation of the protein was reduced, hindering the substrate’s entry into the active cavity and, thereby, inhibiting protein activity through allosteric effects. Ultimately, Markov state models (MSM) were used to identify and classify the metastable states of proteins, revealing the transition times between different conformational states. In summary, this study provides theoretical insights into the inhibition mechanisms of 5LOX by AKBA and NDGA, offering new perspectives for the development of novel inhibitors specifically targeting 5LOX, with potential implications for anti-inflammatory drug development.


Introduction
Inflammation is a complex physiological response of the body to tissue damage, infection, or other harmful stimuli [1].While the inflammatory response is crucial for protecting the body, its dysregulation can lead to chronic inflammatory diseases, such as rheumatoid arthritis [2], asthma, and inflammatory bowel disease [3].Arachidonic acid (AA) metabolism plays a key role in inflammation [4], with 5-lipoxygenase (5LOX) being the central enzyme in this pathway, catalyzing the conversion of arachidonic acid to leukotrienes [5].Leukotrienes (e.g., LTB4, LTC4, LTD4, and LTE4) are potent inflammatory mediators involved in leukocyte recruitment and activation, bronchoconstriction, and increased vascular permeability [6].The aberrant secretion of 5LOX and leukotrienes is closely associated with various inflammatory and allergic diseases, including asthma [7], rheumatoid arthritis [8], and psoriasis [9].Additionally, the products of the 5LOX pathway are linked to diseases such as cancer [10] and atherosclerosis [11].Therefore, 5LOX is considered a promising therapeutic target.
5LOX is composed of 673 amino acids (Figure 1A) and includes two structural domains.The N-terminal regulatory domain, consisting of residues 1-112, is mainly composed of beta sheets, while the C-terminal catalytic domain is significantly larger, spanning residues 125-673.The structural basis of 5LOX function involves non-heme iron in the catalytic site, which is essential for its enzymatic activity [12].The active site of the enzyme is deeply buried and must undergo conformational changes to access the substrate.Understanding the precise structure of 5LOX and its interactions with inhibitors is crucial for the rational design of effective therapeutic agents.
Although AKBA and NDGA exhibit significant therapeutic potential, the detailed mechanisms of their interactions with 5LOX remain unclear.Due to the presence of high-energy barriers, many biological events occur on the millisecond timescale [16].Conventional MD simulations often struggle to adequately sample the relevant conformational space of complex biomolecules [17].Gaussian accelerated molecular dynamics (GaMD) can be used to simulate long-timescale protein dynamics within relatively short computational times.Energy barriers are effectively overcome, significantly enhancing sampling efficiency, which allows for multiple protein conformational states to be explored in a relatively brief period [18,19].Numerous hundreds-of-nanosecond GaMD simulations have been effectively employed to capture significant conformational changes in various biomolecules.These applications include ligand binding [14,[20][21][22][23], protein kinases [24,25], fast-folding proteins [18,26], the CRISPR-Cas9 gene-editing system [27,28], and protein-protein interactions [29][30][31].The Markov state modeling (MSM) sampling method has emerged as an effective means of studying protein dynamics on millisecond timescales [32][33][34].The method uses a series of short molecular dynamics (MD) trajectories to sample the conformational space of a protein, identifying a discrete set of inter-transitioning conformational states.Each state is independent and unaffected by previous states.In this way, MSM can combine multiple independently generated MD trajectories to completely characterize the energy landscape of a protein [35,36].Accelerated molecular dynamics simulations combined with Markov state modeling (MSM) methods show significant potential in the fields of protein engineering and drug discovery [37].
In this study, GaMD was employed to investigate the conformational changes in 5LOX upon binding with the substrate arachidonic acid, the competitive inhibitor NDGA, and the allosteric inhibitor AKBA.Furthermore, the metastable states of proteins were classified using Markov models, revealing the transitions between different conformational states.The molecular mechanisms of 5LOX binding with different ligands were elucidated, providing new insights for the development of anti-inflammatory drugs.

Molecular Docking
The protein-ligand complex of 5LOX and its substrate, arachidonic acid, was subjected to molecular docking to simulate the most probable binding scenario.As Fshown in Figure 2, the substrate arachidonic acid binds tightly to the active pocket of the enzyme and has an alkyl interaction with the proteins LEU368, HIS372, LEU373, ALA410, and ILE415.The substrate binds to the active pocket and has van der Waals interactions with a large number of surrounding residues (Figure S2).The docking results were compared with the S663D Stable-5LOX in complex with arachidonic acid (PDB ID: 3V99), showing that arachidonic acid is docked in within the active cavity close to its position in the crystal structure (Figure S3).After evaluating the structural plausibility, we constructed the corresponding four systems for subsequent molecular dynamics simulations.

Structural Stability of the Four Systems
To evaluate the positional changes in molecular structure and ensure structural stability of each system, we calculated the root-mean-square deviation (RMSD) of the Cα atoms (Figure 3A,B).The average RMSD values for the Apo, 5LOX-AA, and 5LOX-NDGA systems are 1.56 Å, 1.50 Å, and 1.69 Å, respectively (Table 1), indicating minimal changes during the simulation.In contrast, the 5LOX-AKBA system exhibits larger RMSD fluctuations, with an average of 2.03 Å and a variance of 0.49 Å, suggesting significant conformational changes and adjustments in the protein upon binding with AKBA.Over time, the RMSD tends to stabilize, indicating that the protein gradually reaches a new stable state.The stability of ligand binding was assessed by calculating the RMSD values for the regions surrounding the binding sites of NDGA, AKBA, and AA (Figure S4).The RMSD values around AA consistently remained below 0.5 Å, indicating highly stable binding to 5LOX.In contrast, the regions around the binding sites of NDGA and AKBA exhibited greater fluctuations.Furthermore, to reveal changes in protein compactness during the simulation, we analyzed the radius of gyration (Rg).As shown in Figures 2C,D, the Rg value of the 5LOX-NDGA system slightly increases compared to the 5LOX-AA system, indicating that the binding of the competitive inhibitor NDGA induces greater changes in the 5LOX protein.Similarly, the Rg value for the 5LOX-AKBA system shows a noticeable increase compared to the Apo system, suggesting conformational changes upon binding with the allosteric inhibitor.In Figure 3A,C, 5LOX-AKBA showed a significant peak at approximately 280 ns.The protein exhibited a significant shift relative to its initial position, causing substantial fluctuations in the RMSD.After 300 ns, the protein returned to its original position.Further analysis revealed the formation of a new helix at residues 280-295 (Figure S5).To further understand the surface properties of the protein during the simulation, we analyzed the solvent-accessible surface area (SASA).As illustrated in Figure 3E,F, the SASA value of the free protein eventually stabilizes, with an average of 28,597.82Å 2 .The 5LOX-AKBA system shows a significant increase in the SASA compared to the Apo system, with an average value of 29,317.11Å 2 , indicating that the binding of the allosteric inhibitor enhances the protein's hydrophilicity.The SASA value for the 5LOX-NDGA system also significantly increases compared to the 5LOX-AA system, suggesting that the binding of the competitive inhibitor NDGA makes the protein conformation more hydrophilic.To ensure the reliability of our results and to prevent chance errors, we performed three independent simulations.We synthesized the results of these three simulations and found that none of them showed large fluctuations at 500 ns, indicating that the system reached equilibrium after this timepoint.We plotted the

Structural Stability of the Four Systems
To evaluate the positional changes in molecular structure and ensure structural stability of each system, we calculated the root-mean-square deviation (RMSD) of the Cα atoms (Figure 3A,B).The average RMSD values for the Apo, 5LOX-AA, and 5LOX-NDGA systems are 1.56 Å, 1.50 Å, and 1.69 Å, respectively (Table 1), indicating minimal changes during the simulation.In contrast, the 5LOX-AKBA system exhibits larger RMSD fluctuations, with an average of 2.03 Å and a variance of 0.49 Å, suggesting significant conformational changes and adjustments in the protein upon binding with AKBA.Over time, the RMSD tends to stabilize, indicating that the protein gradually reaches a new stable state.The stability of ligand binding was assessed by calculating the RMSD values for the regions surrounding the binding sites of NDGA, AKBA, and AA (Figure S4).The RMSD values around AA consistently remained below 0.5 Å, indicating highly stable binding to 5LOX.In contrast, the regions around the binding sites of NDGA and AKBA exhibited greater fluctuations.Furthermore, to reveal changes in protein compactness during the simulation, we analyzed the radius of gyration (R g ).As shown in Figure 3C,D, the R g value of the 5LOX-NDGA system slightly increases compared to the 5LOX-AA system, indicating that the binding of the competitive inhibitor NDGA induces greater changes in the 5LOX protein.
Similarly, the R g value for the 5LOX-AKBA system shows a noticeable increase compared to the Apo system, suggesting conformational changes upon binding with the allosteric inhibitor.In Figure 3A,C, 5LOX-AKBA showed a significant peak at approximately 280 ns.The protein exhibited a significant shift relative to its initial position, causing substantial fluctuations in the RMSD.After 300 ns, the protein returned to its original position.Further analysis revealed the formation of a new helix at residues 280-295 (Figure S5).To further understand the surface properties of the protein during the simulation, we analyzed the solvent-accessible surface area (SASA).As illustrated in Figure 3E,F, the SASA value of the free protein eventually stabilizes, with an average of 28,597.82Å 2 .The 5LOX-AKBA system shows a significant increase in the SASA compared to the Apo system, with an average value of 29,317.11Å 2 , indicating that the binding of the allosteric inhibitor enhances the protein's hydrophilicity.The SASA value for the 5LOX-NDGA system also significantly increases compared to the 5LOX-AA system, suggesting that the binding of the competitive inhibitor NDGA makes the protein conformation more hydrophilic.To ensure the reliability of our results and to prevent chance errors, we performed three independent simulations.We synthesized the results of these three simulations and found that none of them showed large fluctuations at 500 ns, indicating that the system reached equilibrium after this timepoint.We plotted the results of the three simulations as RMSD and Rg plots (Figure S6) and found no significant changes between the different simulation trajectories.Overall, after 500 ns of molecular dynamics simulation, all four systems exhibit good stability and are suitable for further study [38].
results of the three simulations as RMSD and Rg plots (Figure S6) and found no significant changes between the different simulation trajectories.Overall, after 500 ns of molecular dynamics simulation, all four systems exhibit good stability and are suitable for further study [38].In this analysis, we conducted mutations using FoldX 5.0 [39] on the residues surrounding the inhibitors NDGA and AKBA by substituting the active groups in their side chains with methyl groups, which have a lesser impact on the protein's structure.This method was employed to examine the influence of the active residues on the stability of the protein structure.The findings of this investigation are outlined in Tables 2 and 3.In Table 2, it is evident that the mutations of the residues around NDGA to alanine resulted in minimal changes in energy levels.This observation suggests a strong binding capacity of NDGA to the protein.Conversely, in Table 3, the mutations of the residues LEU66, ARG68, ARG101, HIS125, GLN129, LYS133, and ARG138 to alanine led to an increase in binding energies.This increase implies that these specific sites play a pivotal role in the binding of the inhibitors AKBA and 5LOX.In order to compare the binding ability of substrate AA and orthosteric inhibitor NDGA, the MM-PBSA were calculated as shown in Table 4.The binding energies of AA to 5LOX and NDGA to 5LOX were −28.50 ± 4.11 KJ/mol and −16.18 ± 2.04 KJ/mol, respectively.The substrates AA and 5LOX had smaller binding energies and were more tightly bound.

Flexibility Analysis of the 5LOX
To obtain a clearer understanding of protein volatility changes, RMSF (root-meansquare fluctuation) analysis was conducted for the four systems (Figure 4A).At residues 170-175, located near the active cavity of the protein (Figure 4B), the fluctuation of this segment was higher in the empty protein system but significantly reduced after the binding of the inhibitor AKBA.The increased structural rigidity of this protein segment may prevent substrate entry into the active site, thereby inhibiting protein activity.Conversely, upon binding of substrate AA, this amino acid residue exhibited less fluctuation and relative stability around the active site, ensuring efficient catalysis.At residues 176-195, which form the helix α2 in the protein (Figure 4C), the helix α2 of the 5LOX catalytic domain and the neighboring loop were elevated to create a flexible lid that regulates access to the active site.As shown in Figure 4A, the binding of both inhibitors AKBA and NDGA resulted in increased volatility of protein helix α2.This may directly affect substrate entry or product release, thereby impacting the enzyme's catalytic activity.At residues 280-300 (Figure 4D), where substrate arachidonic acid (AA) binds to the catalytic active site of 5LOX, a local conformational change was observed.The binding of the inhibitor NDGA did not significantly affect the distal structure.However, the fluctuation of the distal loop structure increased upon binding of the inhibitor AKBA, suggesting that AKBA may induce a conformational change in the protein through noncompetitive inhibition.At residues 406-425 (Figure 4E), volatility increased upon AKBA binding.This residue is also located in the helix α2 near the active site.These sites may be critical for reduced protein activity and warrant further investigation.
did not significantly affect the distal structure.However, the fluctuation of the distal loop structure increased upon binding of the inhibitor AKBA, suggesting that AKBA may induce a conformational change in the protein through noncompetitive inhibition.At residues 406-425 (Figure 4E), volatility increased upon AKBA binding.This residue is also located in the helix α2 near the active site.These sites may be critical for reduced protein activity and warrant further investigation.

Conformational Changes during MD Simulations
To compare the conformational changes among the four systems, we analyzed the Dictionary of Secondary Structure of Proteins (DSSP) for two specific regions: residues 150-200 and residues 280-300 of 5LOX.As shown in Figure 5, a significant conformational change occurs in the segment of residues 168-173 in the 5LOX protein.In the Apo system, this segment exists as a loop (Figure 5A), but in the AA-bound state, it gradually transforms into a stable α-helix (Figure 5B).Combined with the RMSF results, the volatility of the 168-173 segment was also significantly reduced (Figure 5D).The structure of the allosteric inhibitor AKBA upon binding is shown in red in Figure 5E.This indicates that, upon substrate binding, this sequence transforms into a more stable α-helix, stabilizing the substrate within the active cavity.The same conformational change was induced upon binding of the inhibitor AKBA.Compared to the null protein, this structure becomes rigid and inhibits the substrate from entering the active cavity, which may be an important reason for the inhibition.Residues 285-290 undergo a similar change (Figure S7).Instead, binding of the inhibitor NDGA occupies the active cavity and causes a change in the conformation of residues 190-195 from a loop to an α-helix (Figure 5C).Upon the binding of the AKBA

Conformational Changes during MD Simulations
To compare the conformational changes among the four systems, we analyzed the Dictionary of Secondary Structure of Proteins (DSSP) for two specific regions: residues 150-200 and residues 280-300 of 5LOX.As shown in Figure 5, a significant conformational change occurs in the segment of residues 168-173 in the 5LOX protein.In the Apo system, this segment exists as a loop (Figure 5A), but in the AA-bound state, it gradually transforms into a stable α-helix (Figure 5B).Combined with the RMSF results, the volatility of the 168-173 segment was also significantly reduced (Figure 5D).The structure of the allosteric inhibitor AKBA upon binding is shown in red in Figure 5E.This indicates that, upon substrate binding, this sequence transforms into a more stable α-helix, stabilizing the substrate within the active cavity.The same conformational change was induced upon binding of the inhibitor AKBA.Compared to the null protein, this structure becomes rigid and inhibits the substrate from entering the active cavity, which may be an important reason for the inhibition.Residues 285-290 undergo a similar change (Figure S7).Instead, binding of the inhibitor NDGA occupies the active cavity and causes a change in the conformation of residues 190-195 from a loop to an α-helix (Figure 5C).Upon the binding of the AKBA inhibitor, a more stable α-helix was observed in both residues, which may be a key factor in the inhibition of the protein's activity.Furthermore, in the AKBA system of replicate simulations, the same structural changes occurred at the same sites (Figure S8), further verifying the statistical significance of our results.
inhibitor, a more stable α-helix was observed in both residues, which may be a key factor in the inhibition of the protein's activity.Furthermore, in the AKBA system of replicate simulations, the same structural changes occurred at the same sites (Figure S8), further verifying the statistical significance of our results.

Dynamical Cross-Correlation Matrix Analysis
The dynamic cross-correlation matrix (DCCM) analysis of all Cα atoms is presented in Figure 6.The relative positions and structure of residues 168-173, 375-425, 375-400, and 285-290 in the 5LOX protein are shown in Figure S9.Significant changes in protein motility were caused by the binding of the allosteric inhibitor AKBA.Residues 168-173 formed a new α-helix upon AKBA binding, and the structure near this residue and the segment of residues 375-425 showed a significant positive correlation of motility (Figure 6D), suggesting that the inhibitor AKBA may indirectly affect 5LOX activity by influencing the conformational stability of the distal region.The covariance matrix of the segments

Dynamical Cross-Correlation Matrix Analysis
The dynamic cross-correlation matrix (DCCM) analysis of all Cα atoms is presented in Figure 6.The relative positions and structure of residues 168-173, 375-425, 375-400, and 285-290 in the 5LOX protein are shown in Figure S9.Significant changes in protein motility were caused by the binding of the allosteric inhibitor AKBA.Residues 168-173 formed a new α-helix upon AKBA binding, and the structure near this residue and the segment of residues 375-425 showed a significant positive correlation of motility (Figure 6D), suggesting that the inhibitor AKBA may indirectly affect 5LOX activity by influencing the conformational stability of the distal region.The covariance matrix of the segments of residues 285-290 and residues 375-400 became lighter in color (Figure 6D), and the negative correlation motions were significantly enhanced, which might have occurred due to relative motions between the structures of these two segments, which in turn affected the binding efficiency of the substrate molecules.This movement further validates the allosteric mechanism of AKBA.
of residues 285-290 and residues 375-400 became lighter in color (Figure 6D), and negative correlation motions were significantly enhanced, which might have occurred to relative motions between the structures of these two segments, which in turn affe the binding efficiency of the substrate molecules.This movement further validates th losteric mechanism of AKBA.

PCA Analysis
The principal component analysis (PCA) of carbon alpha (Cα) atoms was condu for the four systems to obtain relative free-energy maps (Figure 7).In Figure 7, the pro structure corresponds to the lowest energy state, and the active site cavity is highligh in different colors.The respective ratios of PC1 and PC2 are listed in Table 5.Table 6 the volume and depth of the active cavity at the time corresponding to the lowest ene state of the 5LOX protein during the simulation.The depth of the active cavity is measu by calculating the distance from the deepest point within the binding site to the bound of the binding site.

PCA Analysis
The principal component analysis (PCA) of carbon alpha (Cα) atoms was conducted for the four systems to obtain relative free-energy maps (Figure 7).In Figure 7, the protein structure corresponds to the lowest energy state, and the active site cavity is highlighted in different colors.The respective ratios of PC1 and PC2 are listed in Table 5.Table 6 lists the volume and depth of the active cavity at the time corresponding to the lowest energy state of the 5LOX protein during the simulation.The depth of the active cavity is measured by calculating the distance from the deepest point within the binding site to the boundary of the binding site.In the Apo system, the protein reaches its lowest energy states at 352.9 ns (PC1: −15.88,PC2: −44.05) and 398.9 ns (PC1: 41.57, PC2: −5.81), with corresponding active cavity volumes of 200.70 Å 3 and 215.04 Å 3 .These observations suggest that the protein cavity volume tends to stabilize over the course of the simulation.Conversely, in the 5LOX-AA system, the protein attains its lowest energy states at 116.6 ns (PC1: −23.60,PC2: −9.94) and 446.3 ns (PC1: 53.60, PC2: −18.47), with active cavity volumes of 559.10 Å 3 and 416.77Å 3 , respectively.At 446.3 ns, the protein cavity volume is reduced from the previous measurement, indicating that, over longer time scales, the protein may partially revert to a more compact structure while still maintaining an extended state to some extent.Similarly, in the 5LOX-NDGA system, the protein reaches its lowest energy states at 46.4 ns (PC1: −8.76, PC2: 36.05) and 179.1 ns (PC1: −44.44,PC2: −6.57), with active cavity volumes of 354.82 Å 3 and 601.60 Å 3 .These results indicate that NDGA binding increases the volume of the protein's active site cavity, occupying the active site and preventing effective substrate binding.Furthermore, in the 5LOX-AKBA system, the protein shows its lowest energy states at 63.1 ns (PC1: −87.54,PC2: −31.82) and 429.2 ns (PC1: 9.92, PC2: −10.41), with active cavity volumes of 272.90 Å 3 and 188.42 Å 3 , respectively.Upon AKBA binding, the active cavity volume decreases, suggesting that the protein may have closed certain structural domains, thereby preventing substrate access to the active site.This phenomenon aligns with the conformational change in residues 168-173 (from loop to α-helix) and the enhanced negative correlation motion of residues 285-290.These structural changes increase the rigidity of the active site, further reducing the volume of the active cavity and preventing substrate entry.In conclusion, the binding of AKBA resulted in a more stable and compact structure at the binding site and a significant reduction in the volume of the protein's active cavity, validating the mechanism by which AKBA inhibits 5LOX through allosteric inhibition.

Markov Model
To analyze the transitions between different conformational states, the entire Markov model calculation was performed using the Pyemma software package (version 2.5.7)[40].Trajectories of 5000 frames each, totaling 500 ns, were used for the analysis.The secondary structure information of residues 155-195 and 280-295 were selected for feature extraction, and the results were clustered into 10 classes using the K-Means algorithm.These features were then subjected to time-lagged independent component analysis (TICA).The TICA projection and secondary structure distribution with secondary structure coloring are shown in Figure 8A,B.
After binding with the orthosteric inhibitor NDGA and the allosteric inhibitor AKBA, the protein showed four major low-energy regions in TICA space, corresponding to different stable conformational states.To construct the Markov model, MSMs with different lag times were estimated, and the implied timescales (ITS) were calculated (Figure S10A,B).Lag times of 10 and 5 were chosen for MSM estimation.The validity of the MSM was checked by performing the Chapman-Kolmogorov test (Figure S10C,D).The predicted results matched well with the actual results, indicating that the Markov state model is valid.We identified four substates of the protein and calculated the mean first-passage time (MFPT) for transitions.PCCA (Perron-Cluster Cluster Analysis) [41] was used to color the TICA plots to show the distribution of each substate (Figure S11).The different feature spaces occupied by each state indicate good separation between the states.Finally, the constructed Markov state model was visualized.

Markov Model
To analyze the transitions between different conformational states, the entire Markov model calculation was performed using the Pyemma software package (version 2.5.7)[40].Trajectories of 5000 frames each, totaling 500 ns, were used for the analysis.The secondary structure information of residues 155-195 and 280-295 were selected for feature extraction, and the results were clustered into 10 classes using the K-Means algorithm.These features were then subjected to time-lagged independent component analysis (TICA).The TICA projection and secondary structure distribution with secondary structure coloring are shown in Figure 8A,B.After binding with the orthosteric inhibitor NDGA and the allosteric inhibitor AKBA, the protein showed four major low-energy regions in TICA space, corresponding to different stable conformational states.To construct the Markov model, MSMs with different lag times were estimated, and the implied timescales (ITS) were calculated (Figure S10A,B).Lag times of 10 and 5 were chosen for MSM estimation.The validity of the MSM was checked by performing the Chapman-Kolmogorov test (Figure S10C,D).The predicted results matched well with the actual results, indicating that the Markov state model is valid.We identified four substates of the protein and calculated the mean first-passage time (MFPT) for transitions.PCCA (Perron-Cluster Cluster Analysis) [41] was used to color the TICA plots to show the distribution of each substate (Figure S11).The different feature spaces occupied by each state indicate good separation between the states.Finally, the constructed Markov state model was visualized.In the NDGA-bound protein, the red segment of residues 275-295 remained stable without significant changes, while the blue α2 segment of residues 155-195 showed αhelix formation and disappearance.The transition from state S1 to states S2, S3, and S4 occurred more readily than the reverse, indicating that this segment more easily forms a stable α-helix (Figure 8C).Under the action of the allosteric inhibitor AKBA, the α-helix of residues 280-295 changed.The transition from states S1, S2, and S3 to state S4 took less time, indicating that this transition occurred more easily, and this segment readily formed an α-helix (Figure 8D).State S4 also tended to transition back to state S1, with the segment of residues 280-295 in state S1 forming a more stable α-helix structure.The Markov model analysis results were consistent with the previous PCA analysis results.In conclusion, AKBA binding made it easier for these two regions to transition to a stabilized α-helix, significantly affecting the dynamics and stability of 5LOX.Additional simulations were conducted to further investigate the stability of these substates under equilibrium conditions and their potential transitions.Conventional molecular dynamics simulations of the 5LOX-AKBA systems were performed for 300 ns for each substate structure separately.The results indicate that the RMSD values of each state show minimal fluctuations and stabilize (Figure S12), suggesting that each substate remains relatively stable under equilibrium conditions.Further analysis of the simulation trajectories revealed transitions between different structures within each trajectory.The durations of the different states were recorded during the simulations (Figure S13).These results suggest that transitions between substates are feasible under conventional simulation conditions and that these states can stably exist during the simulations.

Preparation of Four Systems
The 3D structures of 5LOX-AKBA (PDB ID: 6N2W) and 5LOX-NGDA (PDB ID: 6NCF) were downloaded from the RCSB protein Data Bank (www.rcsb.org,accessed on 1 March 2024) [42].Subsequently, the 5LOX protein was isolated from the protein-ligand complexes using PyMOL software (version number: PyMOL 2.6) [43].Since the 5LOX protein is a dimer, one of its chains was used in the simulations.Additionally, the structure of the protein residues surrounding NDGA was repaired using homology modeling through the Swiss-Model [44] online platform.Meanwhile, the substrate arachidonic acid (AA) was modeled using Discovery Studio Visualizer (version number: 21.1.0)[45].Furthermore, Gaussian 16 [46], a widely used software package for quantum chemistry calculations, was employed to optimize the molecular structure at the B3LYP/6-31G* level of theory.The optimized structure obtained from Gaussian 16 [46] then served as the starting point for molecular docking studies.Molecular docking [47], a computational method used to simulate interactions between protein molecules and small ligands, was subsequently performed using AutoDock Vina [48].In this study, the ligand AA was docked with 5LOX using AutoDock Vina.
In this study, four simulation systems were constructed, as shown in Table 7: Apo, 5LOX-AA, 5LOX-NGDA, and 5LOX-AKBA.Among the ligands, AA was the substrate arachidonic acid, NDGA was the orthosteric inhibitor NDGA, and AKBA was the allosteric inhibitor AKBA.The target protein in this study was 5LOX, which consists of 673 residues and 10,761 atoms.

Molecular Dynamics Simulations
The complexes of the four systems were performed by the PMEMD engine provided with AMBER22 [49].The forcefield of 5LOX was AMBER ff19SB forcefield [50], and the small molecule field was GAFF2 [51].The MCPB.py [52] tool was used to process Fe ions and construct the forcefield parameters for the metal ions The PME [53] method was used for treatment of long-range electrostatics.The distance between the solute surface and the box was set to 15 Å, and the box was filled with the OPC water model [54].Subsequently, sodium ions were randomly added to the simulation box to achieve neutralization of the system, and the steepest descent method was used to minimize the energy.After that, the 100 ps NVT and 100 ps NPT [55] were performed to maintain the system in a stable environment, and the temperature was kept at 310 K by the Langevin thermostat [56] and Langevin pressostat [57] methods.The simulations involved an initial short conventional molecular dynamics simulation of 50 ns to calculate the GaMD acceleration parameters and GaMD equilibration of the added boost potential for 50 ns.

Gaussian Accelerated Molecular Dynamics (GaMD) Simulation
Gaussian accelerated molecular dynamics (GaMD) simulation is a method independent of collective variables (CV), making it widely applicable and beneficial for studying complex biological systems [18].Its working principle involves adding a harmonic boost potential, which smooths the biomolecular potential energy surface and reduces energy barriers.GaMD significantly accelerates biomolecular simulations by several orders of magnitude [19].The GaMD boost potential follows a Gaussian distribution, allowing for energy reweighting through a cumulant expansion to the second order (i.e., "Gaussian approximation").This results in the accurate reconstruction of the biomolecular free-energy landscape.Finally, Gaussian accelerated molecular dynamics (GaMD) simulations with a timestep of 2 fs and a duration of 500 ns were conducted.The simulations were carried out with random initial atomic velocities for each system at the " dual-boost " level, where one boost potential was applied to the dihedral energetic term and the other to the total potential energy term.The reference energy was established at the minimum value E = V max , while the maximum standard deviation of the boost potential, σ 0 , was designated as 6.0 kcal/mol for both the dihedral and overall potential energy components.To calculate the true free energy, the PyReweighting toolkit [58] was utilized to reweight the trajectory data based on the GaMD simulations.A bin size of 1 Å and a cutoff of 500 frames in each bin were used for the calculation.

MM-PBSA
In this study, the molecular mechanics/Poisson-Boltzmann surface area (MM-PBSA) method was employed to evaluate the binding affinity between ligands and a specific protein.Initially, 5000 frames (500 ns) were selected from the molecular dynamics (MD) simulation, sampling every 10 frames, representing the conformations of both the proteinligand complex and its individual components for subsequent energy calculations.In the MM/PBSA analysis, the binding free energy (∆G bind ) is calculated using the following equations: In these equations, the change in gas-phase molecular mechanical energy is denoted as ∆E MM , while ∆G sol represents the variation in solvation free energy.Due to the minimal conformational alterations observed between the receptor and ligand before and after binding, the T∆S term is omitted.∆E MM comprises three components of molecular mechanics: covalent bonding energy changes (∆E int ), electrostatic energy (∆E ele ), and van der Waals interactions (∆E vdW ).Specifically, ∆E int accounts for changes in bond, angle, and torsion energies.In this study, only the ∆E ele and ∆E vdW components from Equation (3) were analyzed further.The ∆G sol term encompasses the combined effects of polar (∆G pol ) and non-polar (∆G nonpol ) solvation free energies.The polar solvation free energy (∆G pol ) was calculated by solving the linearized Poisson-Boltzmann equation using the PBSA tool within the AMBER 22 software.

Trajectory Analysis
Trajectory analyses, including RMSD, Rg, SASA, RMSF, and distance analyses, were performed using the cpptraj module of Amber22 [59].Cross-correlation matrix analysis can identify protein regions with significant conformational changes, aiding in the explanation of interactions between the inhibitor and the protein [60].During molecular dynamics simulations, the covariance matrix is constructed using the first two eigenvectors to describe the correlation of Cα atoms in the main chain.In this study, the dynamics cross-correlation matrix (DCCM) was calculated using R Studio (R version number: 4.3.3)[61].

PCA and Free-Energy Landscape
Principal component analysis (PCA) is an approach of dimensionality reduction, which can help us to identify those motion modes that have the most direct and profound effect on the overall dynamics of the protein [62,63].The PCA data obtained using the Amber cpptraj program were processed with the ddtpd (version number: ddtpd 1.3) [64] software to obtain the PC1 and PC2 coordinates corresponding to -LnP.These values were then used to plot a free-energy surface map in Origin, which represents the relative magnitude of the free energy for different conformations.

Markov Model
Markov state models (MSMs) play a crucial role in studying the conformational transitions and functional changes in allosteric proteins [65].They can identify and describe multiple conformational states, construct transition probability matrices to quantify the probabilities and times of transitions between these states, and build free-energy landscapes to illustrate the energy distributions of different conformations.This aids in understanding protein stability, transition mechanisms, and the kinetic pathways of allosteric proteins [66].

TICA Dimensionality Reduction Method
In this study, we used time-lagged independent component analysis (TICA) to reduce the dimensionality of high-dimensional datasets.TICA is a powerful technique that captures slow collective motions by maximizing the time-lagged covariance between data points, thus identifying the most relevant features that vary over time [66,67].Firstly, the MD trajectory data were processed to extract relevant features from the secondary structure and relative positions of residues 155-195 and 280-295.These features were then applied to the TICA algorithm.By choosing an appropriating lag time, we ensure that the TICA component captures important temporal correlations in the data.Subsequently, the TICA components were used to project the high-dimensional dataset into a low-dimensional space, effectively reducing the data complexity while preserving the most significant dynamic features.

K-Means Clustering Algorithm
The K-Means clustering algorithm is a widely used, unsupervised, machine learning method for partitioning a dataset into K distinct, non-overlapping clusters.K-Means operates by minimizing the variance within each cluster, ensuring that the datapoints in each cluster are as close as possible to the cluster centroid.In the context of Markov models, K-Means clustering is primarily used to cluster the conformational states in molecular simulation trajectories.Each cluster is defined as a microstate in the Markov model, resulting in a discrete trajectory of microstates.The Markov model can then compute the transition probabilities between the microstates, forming a transition probability matrix [68].The K-Means clustering algorithm was applied in our study to classify the protein conformations obtained from the molecular dynamics simulations.A total of 2-16 different numbers of clusters were subjected to cluster analysis, and the sum of square error (SSE) was recorded for each number of clusters (as shown in Figure S1).The SSE decreases as the number of clusters increases, but the rate of decrease slows down significantly at a cluster number of approximately 10.By analyzing the resulting 10 clusters, different conformational states were identified, and the dynamic behavior of the protein was described.

Determination of Lag Time
The lag time is crucial for constructing Markov state models (MSMs) from molecular dynamics (MD) simulation data.It defines the interval between frames used to calculate transition probabilities, ensuring the model captures the true dynamic behavior of the system.If the lag time is too short, there will not be enough time for transitions between states to decorrelate, causing the MSM to capture only rapid fluctuations rather than true long-term dynamics.Conversely, if the lag time is too long, important intermediate states and transitions may be missed.An appropriate lag time allows for the MSM to exhibit Markovian properties, accurately reproducing system dynamics [65].
To determine the optimal lag time, the Chapman-Kolmogorov (CK) test is used [40].This test evaluates the consistency of transition probabilities over different lag times.It compares the directly computed transition probabilities over a selected lag time with those obtained by multiplying transition probabilities over shorter times.If the MSM satisfies the CK equation, it indicates that the chosen lag time is appropriate, ensuring the model's reliability and accuracy.

Conclusions
In this study, molecular dynamics simulations were utilized to investigate four different systems, Apo, 5LOX-AA, 5LOX-NDGA, and 5LOX-AKBA, examining the protein changes induced by different ligands binding to 5LOX.The results are as follows.The substrate arachidonic acid was tightly bound to the enzyme's active pocket and interacted with specific residues (such as HIS600, PRO569, ALA603, and ALA424) through alkyl interactions.The orthosteric inhibitor NDGA was tightly bound in the protein's active pocket, occupying the active site through competitive inhibition, thereby preventing the effective binding of the substrate.The binding of NDGA increased the volume of the protein's active site cavity, causing structural changes near the active pocket.The binding of the allosteric inhibitor AKBA caused significant changes at the distal active site, resulting in the conformation of residues 168-173 shifting from a loop to an α-helix, and significant negative correlated motions between residues 285-290 and residues 375-400, reducing the distance between these segments.The stable conformation of the protein in the simulation showed a reduced active cavity volume, hindering the substrate's entry into the active cavity, thus inhibiting protein activity through allosteric effects.Finally, the change paths and transfer times were analyzed using Markov models to verify a more stable change structure.In conclusion, this study reveals the interaction mechanisms and secondary structure changes in 5LOX when bound to various ligands.These findings provide valuable insights for developing potential anti-inflammatory drugs.

Figure 1 .
Figure 1.(A) Sequence and the secondary structure diagram of 5LOX.(B) The active around the NDAG inhibitor binding to 5LOX; (C) the active residues around the AKBA i binding to the 5LOX-AKBA complex.

Figure 1 .
Figure 1.(A) Sequence and the secondary structure diagram of 5LOX.(B) The active residues around the NDAG inhibitor binding to 5LOX; (C) the active residues around the AKBA inhibitor binding to the 5LOX-AKBA complex.

Figure 2 .
Figure 2. The binding sites of the substrate arachidonic acid (AA) in 5LOX proteins.(The substrate AA is red, FE 2+ is orange and the active residues are yellow.)

Figure 2 .
Figure 2. The binding sites of the substrate arachidonic acid (AA) in 5LOX proteins.(The substrate AA is red, FE 2+ is orange and the active residues are yellow).

Figure 3 .
Figure 3. Analysis of structural stability.(A) The temporal evolution of the RMSD from their initial structure of the Apo, 5LOX-AA, 5LOX-NDGA, and 5LOX-AKBA systems.(B) Distribution of RMSD values in the four systems.(C) The temporal evolution of the Rg from their initial structure of the Apo, 5LOX-AA, 5LOX-NDGA, and 5LOX-AKBA systems.(D) Distribution of Rg values in the four systems.(E) The temporal evolution of the SASA from their initial structure of the Apo, 5LOX-AA, 5LOX-NDGA, and 5LOX-AKBA systems.(F) Distribution of SASA values in the four systems.The median (the horizontal line in the center), the mean (the black dot), and the interquartile range (the upper and lower edges of the box) are shown in the box plot for each set of data.

Figure 3 .
Figure 3. Analysis of structural stability.(A) The temporal evolution of the RMSD from their initial structure of the Apo, 5LOX-AA, 5LOX-NDGA, and 5LOX-AKBA systems.(B) Distribution of RMSD values in the four systems.(C) The temporal evolution of the R g from their initial structure of the Apo, 5LOX-AA, 5LOX-NDGA, and 5LOX-AKBA systems.(D) Distribution of R g values in the four systems.(E) The temporal evolution of the SASA from their initial structure of the Apo, 5LOX-AA, 5LOX-NDGA, and 5LOX-AKBA systems.(F) Distribution of SASA values in the four systems.The median (the horizontal line in the center), the mean (the black dot), and the interquartile range (the upper and lower edges of the box) are shown in the box plot for each set of data.

Figure 7 .
Figure 7.The free-energy landscape for the following four systems: (A) Apo, (B) 5LOX-AA, (C) 5LOX-NDGA, and (D) 5LOX-AKBA.The active cavity structures of representative conformations in the low-energy region are indicated by different colors.

Figure 7 .
Figure 7.The free-energy landscape for the following four systems: (A) Apo, (B) 5LOX-AA, (C) 5LOX-NDGA, and (D) 5LOX-AKBA.The active cavity structures of representative conformations in the low-energy region are indicated by different colors.

Figure 8 .
Figure 8. (A) Energy landscape colored by secondary structure percentage of 5LOX-NDGA system.(B) Energy landscape colored by secondary structure percentage of 5LOX-AKBA system.(C) Markov state models and mean first passage time of 5LOX-NDGA.(D) Markov state models and mean first passage time of 5LOX-AKBA.(Residues 155-195 are blue and residues 280-295 are red.)

Figure 8 .
Figure 8. (A) Energy landscape colored by secondary structure percentage of 5LOX-NDGA system.(B) Energy landscape colored by secondary structure percentage of 5LOX-AKBA system.(C) Markov state models and mean first passage time of 5LOX-NDGA.(D) Markov state models and mean first passage time of 5LOX-AKBA.(Residues 155-195 are blue and residues 280-295 are red).

Table 1 .
Mean and standard deviation of root-mean-square deviation (RMSD), radius of gyration (R g ), and solvent-accessible surface area (SASA) in the four systems.

Table 2 .
Alanine mutation of residues around NDGA binding.

Table 3 .
Alanine mutation of residues around AKBA binding.

Table 6 .
Active site information for the stabilization structure of the four systems.

Table 7 .
The four systems of molecular dynamics simulation.