MOLECULAR DYNAMICS STUDY OF THE LYSOZYME-BASED DRUG DELIVERY NANOSYSTEMS LOADED WITH ANTIVIRAL DRUGS AND CYANINE DYES 1

Protein-based drug nanocarriers are increasingly recognized as promising candidates for effective drug delivery, owing to a multitude of beneficial advantages over synthetic materials including low cytotoxicity, biocompatibility, biodegradability, abundance, renewability, and high drug loading capacity mediated by diverse functional groups and interactions. In the present study the molecular dynamics simulation was employed to explore the stability of lysozyme-based drug delivery nanosystems functionalized by the antiviral drugs (favipiravir, molnupiravir, nirmatrelvir and ritonavir) and cyanine dyes (AK7-5, AK5-6, AK3-11). A series of 5 ns or 100 ns MD simulations for the top-scored docked drug-dye-protein complexes, obtained using the PatchDock server was performed at 310 K with GROMACS software using the CHARMM General Force Field. The MD results have been analyzed in terms of the parameters, such as the backbone root mean-square deviation, gyration radius, solvent accessible surface area, the root means square fluctuations. The analysis of calculated parameters for the studied systems enabled us to improve the previously acquired molecular docking data. Taken together, the results obtained indicate that Lz-F-AK3-11, Lz-R-AK75, Lz-R-AK56, Lz-N-AK75, Lz-N-AK3-11, and Lz-M-AK75 systems exhibit the highest stability among the examined dye-drug-protein systems and represent potential candidates for the targeted delivery of the explored antiviral agents.

Over the past few decades protein-based drug delivery nanosystems have emerged as a promising avenue in the field of pharmaceutical and biomedical research due to a multitude of their beneficial biochemical and biophysical advantages over synthetic materials such as: i) biocompatibility and biodegradability; ii) remarkable capacity to bind various drugs owing to the presence of numerous binding pockets and functional groups within proteins; iii) safeguarding the drugs from enzymatic breakdown and swift renal excretion; iv) ample availability of proteins sourced from nature v) the ability to attach surface ligands specific to target tissues; vi) efficient and cost-effective synthesis procedures [1][2][3][4][5][6][7], etc.A critical aim in the fabrication of efficient protein-based nanotherapeutics is their optimization, which demands both: i) a profound knowledge of the behavior of the designed nanocarrier in various conditions as well as ii) a comprehensive understanding of potential interactions within the system [8,9].Despite the increasing sophistication of experimental endeavors to design and optimize the structure and dynamics of the protein-based drug delivery nanosystems, such research inevitably encounters intrinsic and practical limitations since the systematic variations in the nanocarrier properties like the type of protein or visualizing agent, protein size, hydrophobicity, surface charge, etc. can be prohibitively time-consuming and expensive.Hence, relying solely on experimental research makes it challenging to elucidate the general biophysical and biochemical principles that connect the effectiveness of a potential protein-based drug delivery system with its composition.Given this, the methods of computer modeling are becoming widespread in the early-stage development and optimization of potential protein-based drug-delivery nanosystems [10][11][12][13][14].
In our previous work, we employed the multiple ligand simultaneous docking technique to evaluate the possibility of functionalization of protein-based nanoparticles by the antiviral drugs and cyanine dyes [15,16].More specifically, we evaluated the suitability of a series of cyanine dyes (four monomethines, six trimethines, seven pentamethines and two heptamethines) to serve as visualizing agents in the drug-protein-dye systems, encompassing four functionally significant proteins (cytochrome c, serum albumin, lysozyme and insulin) and four antiviral drugs (favipiravir (F), molnupiravir (M), nirmatrelvir (N) and ritonavir (R)) [15,16].We identified the ternary systems with the highest dye-protein surface shape complementarity for each group of the examined cyanine dyes and assessed the impact of cyanine dye structure on the stability of the drug-protein-dye complexes [15,16].In continuation of our previous work, in the present study the molecular dynamics simulation was performed to elucidate the stability of the most prospective drug delivery systems obtained by multiple ligand simultaneous docking approach [15,16].To this end, the present study aimed to characterize the effects of three cyanine dyes (AK7-5, AK5-6, AK3-11) and 4 antiviral drugs (favipiravir, molnupiravir, nirmatrelvir and ritonavir) on the structure and dynamics of lysozyme (Lz) in the drug-protein-dye systems.

MOLECULAR DYNAMICS SIMULATIONS
The top-scored docked drug-dye-protein complexes, obtained using the PatchDock server as described in [15,16], were considered as a starting structure for the molecular dynamic's simulations.First, the .pdb-file of docked system was corrected in PDB Editor [17] to add the missing chain ID for ligands.Next, this .pdb-filewas uploaded in the webbased graphical interface CHARMM-GUI [18] to obtain the input files for MD calculations.The .itp files of cyanines (AK7-5, AK5-6, AK3-11) and antiviral drugs (favipiravir, molnupiravir, nirmatrelvir and ritonavir) were constructed from the dye.mol2 files, using the CHARMM General Force Field, followed by replacing the dye/drug partial charges with those assigned by RESP ESP charge Derive Server [19].The drug-dye-protein complexes were solvated in a rectangular box with a minimum distance of 10 Å from the protein to the box edges and 0.15 M NaCl (neutralizing ions) were added to the systems.The TIP3P water model was used.The molecular dynamics simulations and analysis of the trajectories were performed using the GROMACS software (version 2023.3) with the CHARMM36m force field in the NPT ensemble with the time step for MD simulations 2 fs.The calculations were performed at a temperature of 310 K.The minimization and equilibration of the systems were carried out during 50000 and 125000 steps, respectively.The V-rescale thermostat provided the constant temperature conditions.The time interval for MD calculations was 5 ns for stability complex preview and 100 ns for the most stable systems.The molecular dynamics trajectories were corrected after the MD run, using the gmx trjconv GROMACS command [20,21].The GROMACS commands gmx rms, gmx gyrate, gmx rmsf were used to calculate the protein backbone root-mean-square deviation (RMSD), protein radius of gyration (Rg), root-mean-square fluctuations of the C-alpha atoms (RMSF).Visualization of the snapshots of the MD runs and analysis of the protein secondary structures, calculating distances between the dye/drug and protein centers of geometry and protein solvent-accessible surface area (SASA) were performed in VMD.

RESULTS AND DISCUSSION
Lysozyme, a globular functional protein containing 129 amino acids in the monomeric state, is currently widely employed in the development of drug-delivery systems, due to lysozyme ability to bind the drugs and facilitate their accumulation at target sites [22,23].Beyond its drug-delivery role, lysozyme exhibits various intrinsic medicinal properties, including anti-inflammatory, antiviral, immune-modulatory, antihistaminic, and antitumor activities [24,25].
In our previous studies we evaluated the possibility of creating the lysozyme-based nanosystems carrying antiviral drugs and cyanine dyes as visualizing agents [15,16].Varying both the structure of the cyanine dyes (mono-, tri-, pentaand heptamethines) and the antiviral agents (favipiravir, molnupiravir, nirmatrelvir and ritonavir), we identified three cyanines (AK7-5, AK5-6, AK3-11) as the most promising candidates for the fabrication of multicomponent drugdelivery nanosystems [15,16].To verify the stability of the obtained ternary complexes, the next logical step was to perform the molecular dynamics simulation for the obtained systems with the highest docking score.Initially, the 5 ns simulation was performed for all dye-lysozyme-drug complexes.The representative snapshots are shown in Figure 1.At the first stage of our study we investigated the overall stability of the systems during the simulation time through analyzing the time course evolution of the set of structural parameters such as the protein backbone root-meansquare deviation (RMSD), protein radius of gyration (R g ) and root-mean-square fluctuations of the C-alpha atoms time dependence of RMSD for the Lz-F-AK75 is characterized by the two periods: 1) slight fluctuations ~ 0.07 nm during the first 2.5 ns of simulation and 2) increase of RMSD values to 0.12-0.14nm during the next 0.5 ns followed by the fluctuations at this level during the remaining time of the simulation.For the Lz-F-AK56 system, a substantial increase of RMSD values to 0.15 nm was detected starting from ~4.6 ns.The Lz-F-AK75, Lz-N-AK56 and Lz-M-AK3-11 produced less stable trajectories in comparison with the other systems.The radius of gyration (R g ) serves as an indicator of the structural compactness of a protein.A lower R g indicates a more tightly packed polypeptide chain, while a higher value of the radius of gyration suggests a more open protein structure.Figure 3 shows the minor fluctuation of the radius of gyration for all examined systems.Notably, for the Lz-F-AK75 drug-protein-dye system the R g value was decreased from 1.44 to 1.41 nm at 2.3 ns of simulation.After 2.3 ns, the R g seems to be stabilized and equilibrated during the entire simulation period of 5 ns.The obtained results indicate the strongest fluctuations in the systems containing ritonavir.To determine the dynamic behavior of amino acid residues, the RMSF values of the Calpha atoms of lysozyme were calculated (Figure 3).The RMSF values for the majority of the protein residues ranged from 0.05 to 0.15 nm during the simulation.The exception was the residues 105-115 for all systems under study and the residues 45-54 for the ritonavir-and nirmatrelvir-containing systems.Taken together, the comparison of the calculated parameters such as the protein backbone root-mean-square deviation (RMSD), protein radius of gyration (Rg) root-mean-square fluctuations of the C-alpha atoms (RMSF) and solvent accessible surface area for the investigated systems allowed us to determine the most stable systems, viz.Lz-F-AK3-11, Lz-R-AK75, Lz-R-AK56, Lz-N-AK75, Lz-N-AK3-11 and Lz-M-AK75.The obtained results indicate that during the 5 ns MD simulation, the antiviral drugs and cyanines undergo the changes in their binding interactions in the docked complexes with the protein to adjust to the potential energy function from the MD force field resulting in the destabilization of the drug-dye-protein interactions and unbinding either the drug or the dye from the protein for some systems.Notably, to assess the efficacy of 5 ns MD simulations in determining the stability of the dye-protein-drug complexes, we conducted additional 100 ns MD simulations for the three chosen systems (Figure 5).As shown in Figure 5, the Lz-R-AK75 exhibits the highest stability among the three selected complexes, as both the cyanine dye and the antiviral drug remain bound to the lysozyme during the entire simulation time.On the contrary, in the case of Lz-F-AK75, the antiviral drug detached from the protein at the beginning of the simulation.Additionally, the stability Lz-F-AK56 complex is lower compared to Lz-R-AK75, as we observed the detachment of dye molecule in the middle of the 100 ns simulation, followed by the binding of AK56 to the Lz-ritonavir complex in the end of MD trajectory.The time course evolution of RMSD parameters confirmed the highest stability of the Lz-R-AK75 drugprotein-dye system (Figure 6).

CONCLUSIONS
In conclusion, in the present study we evaluated the potential stability of the drug-dye-lysozyme complexes functionalized with antiviral drugs (favipiravir, molnupiravir, nirmatrelvir and ritonavir) and the cyanine dyes (AK3-11, AK7-5 and AK5-6).The 5-ns and 100 ns MD simulations were performed to elucidate the stability of the top-scored docked drug-dye-protein complexes obtained using the molecular docking method in our previous work.It was found that the performed molecular dynamics simulations were highly effective in the improvement of the results obtained by the molecular docking studies.The analysis of calculated parameters, including the protein backbone root-mean-square deviation (RMSD), protein radius of gyration (Rg), root-mean-square fluctuations of the C-alpha atoms (RMSF), and solvent accessible surface area, for the studied systems enabled the identification of the most stable systems such as Lz-F-AK3-11, Lz-R-AK75, Lz-R-AK56, Lz-N-AK75, Lz-N-AK3-11, and Lz-M-AK75.

Figure 1 .
Figure 1.Representative snapshots of the drug-lysozyme-dye complexes.Cyanine dyes are colored green, and antiviral drugs are marked as blue.

Figure 2 .
Figure 2. Time course evolution of the root-mean-square deviation and the radius of gyration for drug-lysozyme-dye complexes As seen in Fig. 2, the calculated RMSD values for the potential drug-delivery nanosystems do not exceed 0.1 nm, except for short-time fluctuations of Lz-F-AK75, Lz-F-AK56, Lz-N-AK56 and Lz-M-AK3-11.More specifically, the

Figure 3 .
Figure 3.Time course evolution of the root-mean-square deviation for the drug-lysozyme-dye complexes To assess the alterations in the environment of protein residues during the simulation, we concentrated on the solvent-accessible surface area per residue SASA (Figure.4).It was observed that the SASA of the majority of complexes under study fluctuated within the range 71-75, while the SASA of Lz-M-AK3-11 and Lz-R-AK3-11 slightly increased during the simulation.The SASA of the Lz-N-AK56 complex increased from 72 to 78 during the initial 3 ns of simulation and maintained equilibrium during the rest of the simulation time.Taken together, the comparison of the calculated parameters such as the protein backbone root-mean-square deviation (RMSD), protein radius of gyration (Rg) root-mean-square fluctuations of the C-alpha atoms (RMSF) and solvent accessible surface area for the investigated systems allowed us to determine the most stable systems, viz.Lz-F-AK3-11, Lz-R-AK75, Lz-R-AK56, Lz-N-AK75, Lz-N-AK3-11 and Lz-M-AK75.The obtained results indicate that during the 5 ns MD simulation, the antiviral drugs and cyanines undergo the changes in their binding interactions in the docked complexes with the protein to adjust to the potential energy function from the MD force field resulting in the destabilization of the drug-dye-protein interactions and unbinding either the drug or the dye from the protein for some systems.Notably, to assess the efficacy of 5 ns MD simulations in determining the stability of the dye-protein-drug complexes, we conducted additional 100 ns MD simulations for the three chosen systems (Figure5).

Figure 4 .Figure 5 .
Figure 4. Time course evolution of the solvent-accessible surface area for the drug-lysozyme-dye complexes

Figure 6 .
Figure 6.Time course evolution of the MD parameter for the drug-lysozyme-dye complexes obtained during the 100 ns simulation