Molecular Dynamics Simulation

Background

Smart responsive hydrogels, as a class of "intelligent" materials that undergo physicochemical property changes in response to external stimuli (such as temperature, pH, and light), have shown immense potential in cutting-edge biomedical fields like drug controlled-release, tissue engineering, and biosensing. Among them, thermosensitive hydrogels have garnered significant attention for wound dressing applications. However, most current thermosensitive polymers are synthetic materials, and their biocompatibility and biodegradability remain key challenges for clinical application.

Chitosan, as a natural polymer with a wide range of sources, excellent biocompatibility and antimicrobial peptide function, are an ideal cornerstone for constructing novel biomedical materials. This project innovatively selects hydroxybutyl chitosan (L-HBC) as the core matrix, aiming to utilize its unique low-temperature liquid to body-temperature gel phase transition property to construct an injectable engineered bacteria encapsulation system that can perfectly conform to irregular wounds. Although experiments have preliminarily confirmed the macroscopic thermosensitive properties of L-HBC, the underlying molecular-level driving mechanism—that is, how temperature changes trigger conformational changes in polymer chains, and how the interaction networks between polymer and solvent, and between polymer and polymer, reorganize—remains unclear.

To investigate the microscopic mechanism of the thermosensitive properties of our L-HBC hydrogel material, we plan to perform molecular dynamics simulations on the entire gel system to explore its microscopic processes at different temperature fields. The thermosensitive property of L-HBC can be summarized as being liquid at 4°C and below, and solid at room temperature. Therefore, we have chosen 4°C (277.15K) and 37°C (310.15K) as two different simulation temperatures. We hope to verify that the thermosensitivity of L-HBC originates from an entropy-driven dehydration process, and its macroscopic phase transition is co-driven by the disruption of the polymer-water hydrogen bond network and the enhancement of inter-polymer chain interactions upon heating.

In addition, this study aims to elucidate the feasibility of L-HBC hydrogel as an encapsulation carrier for living cells (such as engineered yeast) from a molecular level. An ideal cell encapsulation matrix must achieve "gentle immobilization" of the cells—that is, provide sufficient physical anchoring to prevent cell leakage. We speculate that the molecular basis for this affinity primarily stems from the dynamic hydrogen bond network formed between the gel matrix and the cell surface (rich in proteins). To this end, we selected the antimicrobial peptide Pexiganan as a key model for the cell surface. Through molecular dynamics simulations, we will systematically investigate the formation, stability, dynamic characteristics, and thermodynamic driving forces of the hydrogen bond network between L-HBC and the peptide chain. This part of the research aims to reveal the intrinsic mechanism of "gentle immobilization" of bioactive molecules by the hydrogel at the atomic scale, providing a solid theoretical foundation for its application as a high-performance cell encapsulation platform.

Model Assumptions

Before conducting the aforementioned molecular dynamics simulations, we must establish a series of clear assumptions and approximations. This helps us understand the scope and limitations of the results.

1. Core Scientific Hypothesis

We hypothesize that the macroscopic properties of the L-HBC hydrogel observed in wet experiments (thermosensitive gelation, encapsulation ability) are determined by the dynamic competition of non-covalent interaction networks at the atomic scale. The purpose of the simulation is to verify:

1. (1) Thermosensitivity originates from a thermodynamic switch between an enthalpy-driven process (polymer-water hydrogen bonds stabilizing the system at low temperature) and an entropy-driven process (dehydration leading to polymer chain aggregation at high temperature).
2. (2) Bioaffinity arises from the formation of an energetically and entropically favorable, stable hydrogen bond network between L-HBC and the Pexiganan peptide chain through an interfacial desolvation process.

2. Approximations in Force Field and Physical Model

We employ the GAFF (General Amber Force Field) to describe the entire L-HBC and Pexiganan peptide chain composite system. The core assumption behind this choice is that for questions like "why they bind" and "how they bind," accurately describing the interaction modes (e.g., hydrogen bonds) and relative affinity trends between molecules is more crucial than calculating the absolute binding free energy. Although GAFF is not as optimized for specific protein secondary structures as dedicated protein force fields, it ensures consistency and compatibility in describing the forces between the two different chemical entities—the polymer and the peptide chain. This is vital for studying the interactions at their interface, effectively avoiding computational artifacts that might arise from mixing different force fields, thereby providing a reliable microscopic mechanism for the efficient loading observed in wet experiments.

3. Assumptions on System Scale and Representativeness

Our simulation system consists of 6 L-HBC chains with a degree of polymerization of 100 and 10 Pexiganan peptide chains, which is not a true macroscopic hydrogel. We assume this oligomeric system is a "minimal functional unit." That is, the underlying physicochemical principles that drive the collapse or expansion of these 6 chains at different temperatures are the same as those that drive the formation or dissolution of the macroscopic gel network. Similarly, the binding mode of the 10 peptide chains with this "microgel" core can represent the anchoring mechanism of biomolecules within the pores of a real gel. This assumption allows us to capture the core mechanistic insights necessary to explain the phenomena observed in wet experiments at a feasible computational cost, rather than pursuing a quantitative reproduction of the material's macroscopic properties.

4. Assumption on Timescale

Our simulations are conducted on the nanosecond (ns) timescale, which is far shorter than the actual time for macroscopic gelation. We assume that the key molecular events upon which our scientific hypothesis relies—such as the rapid reorganization of hydrogen bond networks, the initial collapse trend of chain conformations, and the rapid adsorption of peptide chains—will all manifest clearly within the nanosecond scale. The purpose of the simulation is not to reproduce the complete equilibrium process but to capture the initial kinetic trends of the system in response to external stimuli (temperature, addition of peptide chains). For example, observing a continuously non-converging RMSD at 37°C is, in itself, strong support for the "entropy-driven continuous collapse" hypothesis, perfectly explaining the macroscopic sol-gel transition process. Similarly, the formation of a stable hydrogen bond network between the peptide and sugar chains within 4 ns is sufficient to demonstrate their high affinity, providing direct theoretical evidence for the "gentle immobilization" observed in wet experiments.

Molecular Dynamics Simulation to Investigate the Thermosensitive Properties of L-HBC

To investigate the thermosensitive mechanism of the L-HBC material at the atomic level, we employed all-atom molecular dynamics (MD) simulations. This computational approach allows us to track the dynamic trajectories of every atom in the system, enabling the direct observation of conformational changes in the L-HBC polymer chains, intermolecular interactions, and solvation behavior at different temperatures. The entire simulation workflow follows a rigorous, multi-stage protocol designed to ensure that the final data used for analysis accurately represents the system's properties at thermodynamic equilibrium.

Before commencing the "production simulation" for data collection, the initially constructed molecular system must undergo a series of meticulous "preprocessing and equilibration" steps. This critical preparatory phase includes energy minimization to resolve unfavorable atomic contacts, followed by a stepwise relaxation under different ensembles (such as NVT and NPT) to allow the system's temperature, pressure, and density to reach and stabilize at their target values. Only after the system is fully equilibrated, removing any artifacts from the initial configuration, do we proceed with the extended production simulation to ensure the collected molecular trajectories are accurate for subsequent property analysis.

Preprocessing

Initial Structure Construction

L-HBC is obtained by randomly modifying HBC monomers (degree of polymerization 100) with L-DOPA at a 30% ratio. We first used the chemical drawing software Chemdraw to draw a polysaccharide with a degree of polymerization of 10, randomly modified the amino sites with L-DOPA, and generated a 3D structure, as shown below:

Considering the molecular weight limitation of Chemdraw, we proceeded to use Materials Studio to generate L-HBC with a degree of polymerization of 100. First, we converted the ".mol" file generated by Chemdraw into a ".mol2" file using Openbabel^[1]. Then, we used Materials Studio to build the 100-degree-of-polymerization L-HBC, added hydrogens, and optimized the 3D structure, as shown below:

Topology File Generation and Single-Chain Pre-equilibration

To generate the topology file for a single L-HBC chain, we first used L-HBC (in ".mol2" format) as input for Sobtop^[2]. The program was set to automatically recognize atom types for the general GAFF force field^[3]. For atoms not covered by the GAFF force field, their types were automatically supplemented using the Universal Force Field (UFF). For the calculation of atomic charges, considering the extremely high cost of high-precision quantum chemistry calculations for the entire long chain, we adopted an efficient alternative: using the built-in interface of Sobtop, we employed the Openbabel program to calculate MMFF94 atomic charges for all atoms on the L-HBC chain. According to the Sobtop author's evaluation, although MMFF94 charges are inferior to RESP charges in reproducing electrostatic potential, their quality is acceptable for simulations of organic systems and they have the advantage of not being dependent on a specific conformation.

Furthermore, all bond, angle, dihedral, and improper dihedral parameters were directly matched and assigned from the parameter library that accompanies the GAFF force field, built into Sobtop. For the few parameters missing from the library, Sobtop estimated them based on standard chemical rules.

Through the above process, we generated a complete and self-consistent set of GROMACS topology files and corresponding coordinate files for a single L-HBC chain. To obtain a lower-energy, more relaxed single-chain structure to serve as a "building block" for the subsequent multi-chain system, we performed a thorough pre-equilibration of a single L-HBC chain in a vacuum.

This process was carried out in a huge cubic box with a side length of 200 nm to prevent the chain from interacting with its periodic images. The pre-equilibration procedure included:

• 1. Energy Minimization: 50,000 steps of energy minimization using the steepest descent method to eliminate severe spatial clashes in the initial model.
• 2. Simulated Annealing: Subsequently, a 200 ps simulated annealing dynamics run was performed. During this process, the system temperature was linearly increased from 300 K to 600 K and then linearly cooled back down to 300 K. This "high-temperature quenching" process effectively helps the polymer chain overcome local energy barriers and explore a wider conformational space.
• 3. Constant Temperature Equilibration: Finally, the dynamics simulation was continued for another 100 ps at a constant temperature of 300 K to ensure the stability of the chain conformation.

All simulations were performed using the GROMACS 2023.5 software package^[4]. We also sincerely thank Teacher Ze Liu for providing the computing resources: an Ubuntu-18.04 server equipped with two NVIDIA RTX 3090 GPUs. The final pre-equilibrated single-chain structure is shown in PyMOL^[5] as follows:

Constructing the Six-Chain Polysaccharide Aggregate Molecular System

To simulate the aggregation behavior of L-HBC in solution, we first needed to construct an initial conformation containing multiple polymer chains. This step was accomplished using the Packmol^[6] software package. Packmol is capable of randomly packing pre-defined molecules into a simulation box according to specified spatial constraints.

Our specific construction process was as follows: First, we used 6 copies of the pre-equilibrated single-chain L-HBC structure as input. Then, we randomly placed these 6 chains with arbitrary orientations inside a cubic box with a side length of 260 nm. We set a minimum tolerance distance of 4.5 Å between molecules to avoid excessively severe atomic overlaps in the initial structure. This parameter was determined after multiple tests to be optimal; a container volume that is too small or a minimum tolerance that is too low can easily cause the system to crash during simulation. Through this process, we obtained a file containing 6 L-HBC chains with a relatively uniform initial distribution, which served as the starting point for all subsequent GROMACS simulations. The structure of the 6 L-HBC chains is shown below:

Production Simulation

To safely relax the relatively crowded initial system constructed by Packmol, we employed a procedure that included energy minimization and a two-stage equilibration with position restraints:

• 1. Energy Minimization: 50,000 steps of energy minimization using the steepest descent method until the maximum force on any atom was less than 100.0 kJ mol⁻¹ nm⁻¹, to eliminate severe atomic overlaps.
• 2. NVT Equilibration: In the constant particle number, volume, and temperature (NVT) ensemble, the system temperature was stabilized at the target temperature using the V-rescale thermostat. During this stage, a position restraint of 1000 kJ mol⁻¹ nm⁻² was applied to all heavy atoms of the L-HBC chains for 100 ps. This step aimed to allow water and ion molecules to find their equilibrium positions around the fixed polymer backbones.
• 3. NPT Equilibration: In the constant particle number, pressure, and temperature (NPT) ensemble, the system pressure was stabilized at 1 bar using the more gentle C-rescale barostat for 100 ps. This stage continued to maintain position restraints on the polymer chains, with the goal of allowing the system's density to reach a reasonable equilibrium value while avoiding disruption of the system structure due to drastic changes in box volume.

After the system was fully equilibrated, all position restraints were removed. We then conducted two independent sets of production simulations, each for 1 ns and 10 ns, at two different temperatures: 4°C (277.15 K) and 37°C (310.15 K). During the production phase, the Parrinello-Rahman barostat was used to obtain the correct NPT ensemble. A time step of 2 fs was used, and all bonds involving hydrogen were constrained using the LINCS algorithm. Long-range electrostatic interactions were handled by the PME method. However, due to computational and storage limitations, the 10 ns simulation was actually run for 8 ns. The following are demonstration videos of the 1 ns and 8 ns simulations at 4°C and 37°C, respectively:

Results and Analysis

To quantify the molecular dynamics simulation process of the interaction of the 6 L-HBC chains, we plotted the changes in the system's RMSD and radius of gyration (Rg) during the 1 ns and 8 ns simulations at 4°C and 37°C. The RMSD change is shown in the schematic diagram below:

From the visualization process and image analysis, it is evident that a 1 ns simulation time is insufficient for the system to converge at either temperature. Further analysis of the 8 ns RMSD plots shows that the simulation at 4°C tends to plateau after 5 ns, but at 37°C, it has not yet plateaued at 8 ns (and may still be slowly increasing).

This indicates that L-HBC has the following thermosensitive properties:

• 1. Structural Stability at Low Temperature: At 4°C, the RMSD quickly reaches a plateau, indicating that the six-chain polysaccharide aggregate rapidly finds one or several low-energy, relatively stable conformational states and primarily undergoes small-scale thermal motion around this stable state for the remainder of the simulation. This suggests that at low temperatures, intermolecular interactions (especially hydrogen bonds formed with water molecules) are strong, restricting large-scale movement of the chains.
• 2. Conformational Exploration at High Temperature: At 37°C, the RMSD continuously fails to plateau, indicating that the system has not yet reached its final equilibrium conformation. The molecular chains possess higher kinetic energy, sufficient to cross energy barriers, and are actively undergoing conformational sampling. They are constantly folding, rearranging, and searching for a new, more stable state at the higher temperature. This itself is a strong signal of thermosensitive behavior, as it shows that the increase in temperature significantly alters the system's dynamics and energy landscape.

Since the system clearly did not converge in the 1 ns simulations, we only plotted the change in the radius of gyration for the 8 ns simulations:

Observing the plots, it is clear that for the 4°C simulation, the change in Rg is gentle, whereas at 37°C, Rg shows a downward trend. This indicates that at 4°C, the overall size of the aggregate did not undergo major changes, maintaining a relatively loose and extended state. This is likely due to the stable hydrogen bond network formed between the polysaccharide chains and a large number of water molecules. At 37°C, however, the continuous decrease in Rg indicates that the six-chain polysaccharide aggregate is undergoing collapse—the molecular chains are moving closer to each other, and the entire aggregate is becoming increasingly compact.

The simulation visualizations and the analysis of RMSD and Rg plots effectively explain the macroscopically observed thermosensitive properties of the hydrogel:

• 1. 4°C: Enthalpy-driven, polymer-water interactions dominate. The chains are fully solvated by water molecules, remaining separate from each other, and the system is in a liquid state.
• 2. 37°C: Entropy-driven, polymer-polymer interactions dominate. To release ordered water molecules and increase the total entropy of the system, the polymer chains undergo dehydration and collapse, cross-linking with each other to form a network, and the system transitions to a gel state.

Summary

The dynamics simulation of L-HBC at different temperature fields has revealed the microscopic nature of its thermosensitive properties, helping to explain for the wet experiments why adjusting the degree of hydroxybutyl substitution to alter the system's hydrophobicity allows for precise and predictable control over the gelation temperature. Furthermore, the analysis of this molecular dynamics simulation indicates that the injection of the gel and the loading of drugs/cells must be carried out at low temperatures (e.g., 4°C), as the molecular network is most open and extended at this temperature, ensuring optimal flowability and encapsulation efficiency.

Molecular Dynamics Simulation of the Interaction between L-HBC and Antimicrobial Peptides

To deeply understand the potential of L-HBC as a delivery vehicle for antimicrobial peptides (AMPs) and to uncover the key molecular mechanisms governing its controlled release, we employed all-atom molecular dynamics (MD) simulations to study the interaction between L-HBC and AMPs. The core objective of this study is to elucidate, at the atomic level, the binding modes, driving forces, and the influence of environmental factors, such as temperature, on this interaction. Through these simulations, we aim to quantify binding affinity, identify critical amino acid residues and polymer segments, and observe the conformational changes of the AMP during its association or dissociation with L-HBC. This information is vital for optimizing the material's design and predicting its biological function.

As with any rigorous computational study, our simulations follow a standardized, multi-stage workflow to ensure the reliability and physical significance of the results. This workflow begins with a meticulous "system preprocessing and equilibration" phase. The objective of this phase is to construct an initial model containing L-HBC, AMPs, water, and ions, and subsequently bring it to a stress-free thermodynamic equilibrium through steps including energy minimization, gradual heating, and ensemble relaxation. Building upon this equilibrated foundation, we then conduct an extended "production simulation," during which the trajectories of all atoms are recorded. These trajectories form the basis for all subsequent quantitative analysis, providing a dynamic and detailed microscopic picture of the L-HBC and AMP interaction.

Preprocessing

Construction of Pexiganan Topology File

To simulate the microscopic mechanism of the affinity between Pexiganan and the six-chain polysaccharide aggregate, it was necessary to obtain the structure file for Pexiganan. After an unsuccessful search for Pexiganan-related information in the RCSB Protein Data Bank (RCSB PDB), we predicted its structure based on its sequence information using AlphaFold3 and generated a ".pdb" file (converted from CIF). We then attempted to generate a topology file for Pexiganan using GROMACS. However, whether choosing the AMBER or CHARMM force field, it was difficult to maintain compatibility in GROMACS with our original six-chain polysaccharide aggregate (GAFF force field)—this was mainly evident in the construction of the force field and topology files. Therefore, we again used the Sobtop software to generate the topology file for Pexiganan and used Packmol to package the entire system. For this, the final equilibrated structure of the six-chain polysaccharide aggregate from the (37°C, 8 ns) simulation was used as the L-HBC part of the initial system for the interaction study. We then randomly placed 10 Pexiganan chains in a concentric spherical shell centered on the L-HBC system, with an inner radius of 120 Å and an outer radius of 150 Å. A schematic diagram of the entire system in PyMOL is shown below:

Pre-equilibration

The force field file for this system can be obtained by merging the force field files of L-HBC and Pexiganan. The pre-equilibration of this system was carried out using the same method as before and will not be detailed here again.

Production Simulation

The principle of the production simulation was also identical to the previous one. However, considering computational power and time constraints, we only conducted a 4 ns simulation, which is actually sufficient for investigating the trend of hydrogen bond changes. The simulation process is shown in the video below:

Video 5. The video of showing molecular dynamics simulation Process of Six Chain Polysaccharide and Pexiganan

Results and Analysis

To elucidate the source of affinity between the polysaccharide and the peptide chain from a molecular level, we conducted a detailed hydrogen bond analysis on the 4 ns molecular dynamics simulation trajectory. As a key non-covalent interaction that governs biomolecular recognition, binding, and complex stability, quantifying the formation process and dynamic changes of hydrogen bonds is the most direct and core evidence for revealing the binding strength and mechanism between the two.

We first counted the total number of hydrogen bonds formed between the polysaccharide and peptide molecules as a function of time during the simulation, with the results shown in the figure below:

As shown in the figure, shortly after the simulation began, the number of hydrogen bonds between the two rapidly increased and soon entered a dynamic equilibrium state. Throughout the entire simulation, the system consistently maintained a very substantial number of intermolecular hydrogen bonds, with the total number fluctuating around a stable average of 609.83 ± 17.78. This result reveals two key pieces of information:

• 1. The high average number of approximately 610 hydrogen bonds indicates that the polysaccharide and peptide chains form an extensive and dense interaction network at their contact interface.
• 2. The total number of hydrogen bonds remained at a high level throughout the entire simulation, showing no trend of decreasing or dissociating. This strongly confirms that the two have formed a structurally stable complex. The fluctuations around the average value reflect the dynamic nature of individual hydrogen bonds constantly breaking and reforming, yet the overall integrity of the hydrogen bond network is perfectly maintained.

Although the figure above confirms the existence of a stable interaction, it does not explain why, in an aqueous environment, these two types of molecules prefer to bind to each other rather than to water molecules. To investigate this fundamental thermodynamic driving force, we further analyzed the competitive relationship between "solute-solute" and "solute-water" hydrogen bonds, as shown in the figure below:

The dual Y-axis plot clearly contrasts the changes in the number of two types of hydrogen bonds: polysaccharide-peptide hydrogen bonds (red curve, left axis) and the hydrogen bonds formed by these two molecules with surrounding water molecules (blue curve, right axis). A very clear inverse correlation can be observed from the plot: as the red curve representing complex formation (polysaccharide-peptide H-bonds) rises and stabilizes, the blue curve representing the degree of molecular solvation (H-bonds with water) shows a corresponding significant downward trend.

Summary

This phenomenon reveals a classic Desolvation process, which is one of the core mechanisms of biomolecular recognition and binding. For the polysaccharide and peptide chains to form direct interactions, they must first break some of the hydrogen bonds formed with water molecules on their surfaces, in other words, "shake off" their respective solvation water layers. Therefore, we can conclude that the strong affinity between the polysaccharide and the peptide chain is not only derived from the energetic advantage of their direct interaction but is also largely driven by the favorable entropy gain from releasing surface water molecules. The system sacrifices some interactions with water in exchange for an overall more stable complex state with lower energy and higher entropy. This provides a microscopic theoretical basis for the wet lab experiments that combine antimicrobial peptides and hydrogels for treating diabetic wounds.

Discussion and Outlook

Summary and Analysis

This study has systematically elucidated the thermosensitive phase transition mechanism of hydroxybutyl chitosan (L-HBC) hydrogel at the atomic level and revealed the molecular basis of its high affinity for the antimicrobial peptide Pexiganan. The simulation results strongly confirm our core scientific hypothesis: the thermosensitivity of L-HBC originates from an entropy-driven dehydration process, where the macroscopic phase transition is co-driven by the disruption of the polymer-water hydrogen bond network and the enhancement of inter-polymer chain interactions upon heating. At low temperatures (4°C), the system is enthalpy-dominated, with polymer chains remaining extended due to full solvation by water molecules; at body temperature (37°C), the system shifts to being entropy-dominated, where polymer chains form compact aggregates through dehydration and collapse. Furthermore, we found that L-HBC and Pexiganan peptide chains can rapidly form a large and macroscopically stable hydrogen bond network, a binding process also accompanied by a significant desolvation effect. This provides a solid theoretical basis for the efficient immobilization of bioactive molecules by the hydrogel under mild physiological conditions. This research not only offers profound microscopic mechanistic insights for the design and application of this novel smart material, L-HBC, but also provides important theoretical guidance for the development of drug delivery and cell encapsulation systems based on natural polymers.

Further Discussion

Although this study has reached insightful conclusions, several limitations still exist. First, due to computational resource constraints, the simulation timescales (8 ns for phase transition studies, 4 ns for affinity studies) are relatively limited. In particular, the system at 37°C did not fully reach equilibrium, which might affect the precise determination of the final equilibrium conformation and thermodynamic properties. Second, the GAFF/UFF hybrid force field and MMFF94 charges used in this study are an efficient approximation. Their accuracy is inferior to that of high-precision quantum chemical parameterization for specific chemical groups, which may pose a challenge to the quantitative accuracy of the simulation results. Finally, the size and concentration of the simulated system are simplifications of a real hydrogel system. Its representativeness needs to be verified in the future by constructing larger-scale, longer-timescale models and combining them with experimental methods.

Outlook

During the process of molecular dynamics simulation, the biggest challenges encountered were the drawing of protein and material molecular structures, and the construction of topology and force field files for the polysaccharide-protein system. In practice, these tasks can be solved by fixed algorithms, but currently, there are no end-to-end automated processing tools. On one hand, the tedious file construction and the relatively high barrier to entry of simulation software deter many beginners. On the other hand, AI-designed biomacromolecules are becoming increasingly popular, leading to a growing demand for molecular dynamics simulation as a tool to explain microscopic mechanisms. Therefore, there is an urgent need for an automated tool that can assist users in conducting dynamics simulations and free beginners from the cumbersome task of file construction. We envision incorporating this functionality into AMPilot in the future. This would not only enhance the extensibility of AMPilot but also help transform MD from a static, open-loop computational tool into a dynamic, closed-loop autonomous discovery platform.

References

1. [1] O'Boyle, N. M.; Banck, M.; James, C. A.; Morley, C.; Vandermeersch, T.; Hutchison, G. R. Openbabel: An open chemical toolbox. J. Cheminform. 2011, 3, 33.
2. [2] Lu, T. Sobtop, Version 1.0(dev5). http://sobereva.com/soft/Sobtop (accessed on Sep 18, 2025).
3. [3] Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general amber force field. J. Comput. Chem. 2004, 25, 1157–1174.
4. [4] Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015, 1–2, 19–25.
5. [5] The PyMOL Molecular Graphics System, Version 2.5, Schrödinger, LLC.
6. [6] Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. PACKMOL: A Package for Building Initial Configurations for Molecular Dynamics Simulations. J. Comput. Chem. 2009, 30, 2157–2164.