What Binds Cationic Photosensitizers Better: Brownian Dynamics Reveals Key Interaction Sites on Spike Proteins of SARS-CoV, MERS-CoV, and SARS-CoV-2

We compared the electrostatic properties of the spike proteins (S-proteins) of three coronaviruses, SARS-CoV, MERS-CoV, and SARS-CoV-2, and their interactions with photosensitizers (PSs), octacationic octakis(cholinyl)zinc phthalocyanine (Zn-PcChol8+) and monocationic methylene blue (MB). We found a major common PS binding site at the connection of the S-protein stalk and head. The molecules of Zn-PcChol8+ and MB also form electrostatic encounter complexes with large area of negative electrostatic potential at the head of the S-protein of SARS-CoV-2, between fusion protein and heptad repeat 1 domain. The top of the SARS-CoV spike head demonstrates a notable area of electrostatic contacts with Zn-PcChol8+ and MB that corresponds to the N-terminal domain. The S-protein protomers of SARS-CoV-2 in “open” and “closed” conformations demonstrate different ability to attract PS molecules. In contrast with Zn-PcChol8+, MB possesses the ability to penetrate inside the pocket formed as a result of SARS-CoV-2 receptor binding domain transition into the “open” state. The existence of binding site for cationic PSs common to the S-proteins of SARS-CoV, SARS-CoV-2, and MERS-CoV creates prospects for the wide use of this type of PSs to combat the spread of coronaviruses.


Introduction
The 21st century has been followed by three outbreaks of highly pathogenic coronaviruses [1]. In 2003, the Severe Acute Respiratory Syndrome coronavirus (SARS-CoV) caused an outbreak of atypical pneumonia in 30 countries, more than 8000 people were infected, and 812 died [2]. Between April 2012 and December 2019, the Middle East Respiratory Syndrome coronavirus (MERS-CoV) infected 2499 people in 27 countries and 858 people died [3,4]. Finally, the Severe Acute Respiratory Syndrome type 2 coronavirus (SARS-CoV-2) caused the COVID-19 pandemic, which has already caused 4,265,903 deaths of 200,840,180 infected (at 9 August 2021) [5] and led to global social problems and economic losses. that such an insert is located at a distance of 10-13 nm from the RBD domain [29], it also potentiates the binding of the S-protein to the negatively charged ACE-2 receptor. Thus, electrostatic interactions are involved in the most important initial stages of the replication cycle of coronaviruses-recognition, binding to host cell receptors, and priming.
Computational virology tools contribute greatly to the understanding of viral structure, infectivity and pathogenesis, and design of antiviral drugs. Computational mutagenesis revealed key mutations that affect the electrostatic properties of coronavirus spikes and RBD binding interface [30,31]. Molecular Dynamics (MD) and Monte Carlo simulations shed light on the molecular interactions of the coronavirus RBD domain with host receptors and antibodies [32], and factors that influence the fusion process [33], as well as contributed to the design of antiviral peptides [34].
Recently, we analyzed the electrostatic interactions of the photosensitizer (PS) octacationic Zn-PcChol 8+ with the spike protein of SARS-CoV-2 using BD software ProKSim and proved PS effectiveness in the photodynamic inactivation of SARS-CoV-2 in vitro [35]. In the absence of irradiation, Zn-PcChol 8+ exhibits the lack of antiviral activity towards avian influenza A viruses [36] as well as towards SARS-CoV-2 [35]. The fact that Zn-PcChol 8+ itself does not influence SARS-CoV-2 infectivity coincides with the location of the main binding site of Zn-PcChol 8+ at the junction of the head and the stalk of SARS-CoV-2 S-protein [35] separated from the sites involved in the interactions with the host cells. In contrast, a well-known antimicrobial compound methylene blue (MB), even nonphotoactivated, shows in vitro activity against SARS-CoV-2 (strain IHUMI-3) at submicromolar concentrations [37]. This is apparently due to the ability of MB to inhibit interactions of the S-protein RBD with ACE2 thus blocking SARS-CoV-2 binding and entry [38]. Under irradiation, Zn-PcChol 8+ causes the loss of SARS-CoV-2 infectivity. With a minimal studied concentration of Zn-PcChol 8+ 1 µM and a dose of 692 nm LED light 3.75 J/cm 2 , this PS completely inactivated SARS-CoV-2 with the initial titer 5.00 lgTCID50/mL [35]. The efficiency of photodynamic inactivation of SARS-CoV-2 with Zn-PcChol 8+ is similar to that for another enveloped virus. Indeed, Zn-PcChol 8+ (2 µM) together with white light of 30 W halogen lamp (12 J/cm 2 ) completely inactivated avian influenza A virus H5N8 with initial titer 7.125 lgTCID50/mL [36]. The second compound, MB, in the concentrations of 1.0-10.0 µg/mL (the range about 3-30 µM) with a continuous laser irradiation with wavelength λ = 662 nm (16 J/cm 2 ) fully protected Vero E6 cells infected with 4 lgTCID 50 of SARS-CoV-2 and partly protected from 5 lgTCID 50 of SARS-CoV-2 [39]. The measured MB photodynamic activity against SARS-CoV-2 seems to be quite low taking into account that MB toxicity to SARS-CoV-2 without irradiation estimated in the same study was IC 50 0.22 µg/mL. Consequently, Zn-PcChol 8+ and MB differ in antiviral action probably due to different location of their binding sites. In this study, we compare the electrostatic properties of the S-proteins (the key viral surface proteins) of three coronaviruses, SARS-CoV, MERS-CoV, and SARS-CoV-2, and their interactions with antivirals Zn-PcChol 8+ and MB.

Protein 3D Models
The structural model of SARS-CoV-2 S-protein was adopted from [40]. The models of S-proteins of SARS-CoV and MERS-CoV were based on cryo-EM structures from the Protein Data Bank (PDB) with IDs 6NB3 and 5X58, respectively. The unresolved amino acid residues of the "head" domains of S-proteins were reconstructed using the i-TASSER software [41]. The "stalks" of S-proteins containing HR2, CP, and TM domains (1223-1353 for MERS-CoV and 1105-1255 for SARS-CoV) were absent in initial cryo-EM structures. We predicted their secondary structure on the basis of UniProt sequences A0A140AYW5 for MERS-CoV and P59594 for SARS-CoV using Jpred4 server [42]. Then we designed 3D structures of unresolved "stalk" parts according to the predicted secondary structure using Modeller 9.19 [43] with coiled-coil crystal structure template (PDB ID: 2WPQ) as it was previously done in [44]. The model of the Zn-PcChol 8+ molecule was adopted from [45]. The model of the MB molecule was created using ATB web-service [46]. We used PyMOL [47] for visualization of molecular structures.

Brownian Dynamics Simulations
We performed rigid-body simulations of Brownian diffusion and long-range electrostatic interaction of PS molecules with coronavirus spike proteins using BD software "ProKSim" (Protein Kinetics Simulator, [48,49]). In this approach, spike protein was represented as a low dielectric area (ε = 2) with spatially fixed partial charges. Partial charges on the S-proteins were assigned in accordance with the CHARMM27 force field [50,51] using Gromacs 5.1.4. [52]. Electrostatic field of molecules was calculated using Poisson-Boltzmann formalism [53], as described in detail earlier [54]. The water solvent (dielectric constant ε = 80) with ions was described implicitly. Electrostatic cutoff radius was 3.5 nm. Ionic strength was 100 mM.
Each studied S-protein molecule was placed in a virtual reaction volume with mirror boundary conditions and dimensions of 30 × 30 × 30 nm in such a way that the entire molecule was inside the reaction volume, except its TM and CP domains. In each BD simulation, one PS molecule was initially randomly placed in the reaction volume. The BD simulation then continued until the attracting electrostatic energy reached a predetermined threshold and the resulting structure of the electrostatic encounter complex was saved for further analysis. For Zn-PcChol 8+ the threshold was 8 kT and for MB-2 kT. Twenty thousand independent BD simulations with various initial positions of PS molecules were performed for each PS-spike protein system. Contacts between PS and particular amino acids of S-protein were characterized with a custom Python script, which identifies and counts residues of the S-protein that are within 5 A from the PS molecule throughout the obtained ensemble of structures. Contact probability for particular amino acid residue was calculated as the average fraction of contacts over the entire ensemble of structures. The contact probabilities for all amino acid residues were then visualized on the S-protein surface by color, and also on the S-protein primary sequences. S-protein stalks of SARS-CoV and SARS-CoV-2 show three large areas of negative electrostatic potential, with the largest negatively charged area located at the stalk and head connection, whereas the MERS-CoV S-protein stalk demonstrates only two major areas of negative potential, one at the stalk and head connection and the other at the cytoplasmic domain. The electric charge is highly heterogeneously distributed on the surface of the S-protein head. Nevertheless, areas with a pronounced negative potential are observed on the head of each of the considered coronaviruses. The noteworthy negatively charged area on the NTD-related surface of the SARS-CoV S-protein head is generated by aspartate and glutamate residues Asp15, Asp17, Asp23, Asp24, Asp134, Asp154, Asp243, and Glu131. There are no similar negatively charged regions on the NTD of SARS-CoV-2 and MERS-CoV, which is explained by the absence of some corresponding negatively charged amino acid residues. Instead, the head of MERS-CoV S-protein has a large area of negative electrostatic potential at the linker subdomain between RBD and FP (amino acid residues 676-686, 708-723, and 879-942). SARS-CoV-2 S-protein has the similarly located negatively charged area constituted by amino acid residues of linker between FP and HR1 domains. S-protein stalks of SARS-CoV and SARS-CoV-2 show three large areas of negative electrostatic potential, with the largest negatively charged area located at the stalk and head connection, whereas the MERS-CoV S-protein stalk demonstrates only two major areas of negative potential, one at the stalk and head connection and the other at the cytoplasmic domain. The electric charge is highly heterogeneously distributed on the surface of the S-protein head. Nevertheless, areas with a pronounced negative potential are observed on the head of each of the considered coronaviruses. The noteworthy negatively charged area on the NTD-related surface of the SARS-CoV S-protein head is generated by aspartate and glutamate residues Asp15, Asp17, Asp23, Asp24, Asp134, Asp154, Asp243, and Glu131. There are no similar negatively charged regions on the NTD of SARS-CoV-2 and MERS-CoV, which is explained by the absence of some corresponding negatively charged amino acid residues. Instead, the head of MERS-CoV S-protein has a large area of negative electrostatic potential at the linker subdomain between RBD and FP (amino acid residues 676-686, 708-723, and 879-942). SARS-CoV-2 S-protein has the similarly located negatively charged area constituted by amino acid residues of linker between FP and HR1 domains.

Electrostatic Potential Fields of Coronavirus Spike Proteins and Photosensitizer Molecules
The PS molecules MB and Zn-PcChol8+ demonstrate only positive electrostatic potential surface, with +1 at MB and +8 total charge generated by choline substitutes in Zn-PcChol8+. The PS molecules MB and Zn-PcChol 8+ demonstrate only positive electrostatic potential surface, with +1 at MB and +8 total charge generated by choline substitutes in Zn-PcChol 8+ .

Binding of Photosensitizer Molecules to Coronavirus S-Proteins
We analyzed the binding of octacationic and monocationic PSs to the S-proteins of three coronaviruses, SARS-CoV, SARS-CoV-2, and MERS-CoV. Figures 2-4 show amino acid residues involved in contacts with PSs on the surface of S-proteins (Figures 2 and 3) and in primary sequences with a domain structure (Figure 4), colored depending on the probability of PSs binding. For the SARS-CoV-2 S-protein, amino acid residues involved in contacts with PSs are given for each of the three spike protomers, one of them in the "open" state ( Figure 4(c1,d1)) and two others in the "closed" state ( Figure 4(c2,c3,d2,d3)). Figure 2 shows the areas on the S-protein surface with which Zn-PcChol 8+ molecules form encounter complexes with an electrostatic attraction energy exceeding 8 kT. As MB possesses lower total electric charge than Zn-PcChol 8+ , it has lower values of electrostatic energy when interacting with S-proteins. Therefore, we were able to detect encounter complexes of MB with S-proteins only with a lower threshold (2 kT) of electrostatic attraction energy. In general, the binding sites of this PS resemble those for Zn-PcChol 8+ , although there are some differences. The most prominent attractive areas at the junction of the stalk and the head of S-proteins are present both for Zn-PcChol 8+ and MB. three coronaviruses, SARS-CoV, SARS-CoV-2, and MERS-CoV. Figures 2-4 show amino acid residues involved in contacts with PSs on the surface of S-proteins (Figures 2 and 3) and in primary sequences with a domain structure (Figure 4), colored depending on the probability of PSs binding. For the SARS-CoV-2 S-protein, amino acid residues involved in contacts with PSs are given for each of the three spike protomers, one of them in the "open" state ( Figure 4(c1,d1)) and two others in the "closed" state ( Figure 4(c2,c3,d2,d3)).  As previously found for the interaction of Zn-PcChol 8+ with SARS-CoV-2 [35], the molecules of this PS formed electrostatic encounter complexes with certain areas of the S-protein surface in two other investigated coronaviruses. The only one common binding area for S-protein of all three coronaviruses is the vast area of negative electrostatic potential at the junction of the stem and the head. This area can also be clearly seen in Figure 4 at the linker between HR1 and HR2 domains for all three S-proteins. Note that there is almost no difference in PS binding to this area for the protomers in the "open" and the "closed" states.
Similar to the binding of Zn-PcChol 8+ to SARS-CoV-2 S-protein [35], the molecules of this PS formed electrostatic encounter complexes with the largest area of negative electrostatic potential at the head of the S-protein of MERS-CoV (Figure 2c, upper panel) just before the FP domain (Figure 4e). Unlike SARS-CoV-2 and MERS-CoV, this negatively charged region is absent on the S-protein of SARS-CoV. However, the top of the SARS-CoV spike head demonstrates a notable area of electrostatic contacts with Zn-PcChol 8+ (Figure 2a) and MB (Figure 3a). This area corresponds to the NTD (Figure 4a,b).    Figure 2 shows the areas on the S-protein surface with which Zn-PcChol8+ molecules form encounter complexes with an electrostatic attraction energy exceeding 8 kT. As MB possesses lower total electric charge than Zn-PcChol8+, it has lower values of electrostatic energy when interacting with S-proteins. Therefore, we were able to detect encounter complexes of MB with S-proteins only with a lower threshold (2 kT) of electrostatic attraction energy. In general, the binding sites of this PS resemble those for Zn-PcChol8+,    Figure 2 shows the areas on the S-protein surface with which Zn-PcChol8+ molecules form encounter complexes with an electrostatic attraction energy exceeding 8 kT. As MB possesses lower total electric charge than Zn-PcChol8+, it has lower values of electrostatic energy when interacting with S-proteins. Therefore, we were able to detect encounter complexes of MB with S-proteins only with a lower threshold (2 kT) of electrostatic attraction energy. In general, the binding sites of this PS resemble those for Zn-PcChol8+, The largest area of negative electrostatic potential at the head of the S-protein of SARS-CoV-2 also attracts MB molecules (Figure 3b) as it was for Zn-PcChol 8+ [35]. The attractive ability of the similar area of MERS-CoV S-protein for MB is almost totally reduced (Figure 3c). The S-protein protomers of SARS-CoV-2 in "open" and "closed" conformations demonstrate a different ability to attract PS molecules (Figure 4c,d). The mentioned above largest area of negative electrostatic potential at the head of the S-protein protomers in the "closed" conformation attracts Zn-PcChol 8+ (Figure 4(c2,c3)) and MB (Figure 4(d2,d3)). This area is located just after the FP domain. The S-protein protomer in the "open" conformation ( Figure 4(c1,d1)) hardly attracts PS molecules in this region. To understand the nature of this difference, we analyzed the spatial conformation of distinct amino acid residues generating electrostatic potential in this area for both "open" and "closed" states. Arg847 of the Sprotein protomer in the "open" state is exposed to solution and shields negatively charged residues Asp830, Asp839, Asp843, Asp848 of the same protomer, therefore, positively charged PS molecules hardly form electrostatically favorable complexes with this area (Figure 4(c1,d1) and Figure 5a,b). In the "closed" state, Arg847 forms hydrogen bonds with D574, D568, and D586 of the adjacent S-protein protomer. Thus, Arg847 is not exposed into solution and does not impose such a serious shielding effect on the negative electrostatic field of this area (Figure 5c,d) and PS binding (Figure 4(c2,c3,d2,d3)).
The largest area of negative electrostatic potential at the head of the S-protein of SARS-CoV-2 also attracts MB molecules (Figure 3b) as it was for Zn-PcChol8+ [35]. The attractive ability of the similar area of MERS-CoV S-protein for MB is almost totally reduced (Figure 3c). The S-protein protomers of SARS-CoV-2 in "open" and "closed" conformations demonstrate a different ability to attract PS molecules (Figure 4c,d). The mentioned above largest area of negative electrostatic potential at the head of the S-protein protomers in the "closed" conformation attracts Zn-PcChol8+ (Figure 4(c2,c3)) and MB (Figure 4(d2,d3)). This area is located just after the FP domain. The S-protein protomer in the "open" conformation ( Figure 4(c1,d1)) hardly attracts PS molecules in this region. To understand the nature of this difference, we analyzed the spatial conformation of distinct amino acid residues generating electrostatic potential in this area for both "open" and "closed" states. Arg847 of the S-protein protomer in the "open" state is exposed to solution and shields negatively charged residues Asp830, Asp839, Asp843, Asp848 of the same protomer, therefore, positively charged PS molecules hardly form electrostatically favorable complexes with this area (Figures 4(c1,d1) and 5a,b). In the "closed" state, Arg847 forms hydrogen bonds with D574, D568, and D586 of the adjacent S-protein protomer. Thus, Arg847 is not exposed into solution and does not impose such a serious shielding effect on the negative electrostatic field of this area (Figure 5c,d) and PS binding ( Figure  4(c2,c3,d2,d3)).  In contrast with Zn-PcChol8+, MB possesses the ability to penetrate inside the pocket formed as a result of SARS-CoV-2 RBD domain transition into the "open" state ( Figure  3b). This MB binding site possesses a complex labyrinthine geometry that is composed by amino acid residues of NTD, RBD, and HR1 domains belonging to protomers in the "open" and "closed" states ( Figure 4d). In contrast with Zn-PcChol 8+ , MB possesses the ability to penetrate inside the pocket formed as a result of SARS-CoV-2 RBD domain transition into the "open" state ( Figure 3b). This MB binding site possesses a complex labyrinthine geometry that is composed by amino acid residues of NTD, RBD, and HR1 domains belonging to protomers in the "open" and "closed" states ( Figure 4d).

Discussion
In our study, we used BD to simulate the binding of water-soluble cationic PSs to coronavirus S-proteins that are viral structures mostly protruded into a water environment. The S-proteins of all three studied coronaviruses have a significant negative total charge formed as a result of the excess of the number of negatively charged amino acid residues (Asp + Glu) over positively charged ones (Arg + Lys). This electric charge is highly heterogeneously distributed on the S-protein surface and there are only a few areas of pronounced negative electrostatic potential (colored in red in Figure 1). The largest of them, common to all three coronaviruses, is located at the connection of the S-protein stalk and the head on the linker between the HR1 and HR2 domains, adjacent to HR2. Consequently, it attracts the majority of Zn-PcChol 8+ and MB molecules and constitutes the main binding site for both cationic dyes located at a distance of about 10 nm from the viral membrane. Zn-PcChol 8+ and MB are both the Type II PSs, which produce singlet oxygen under red light irradiation. Earlier we proposed [35] that the location of PS in this binding site promotes effective oxidative damage to both the S-protein and the viral membrane. In vitro studies confirmed Zn-PcChol 8+ to be highly effective in SARS-CoV-2 photodynamic inactivation [35].
The photodynamic activity of MB towards MERS-CoV [55] and SARS-CoV-2 [39] is a well-known experimental fact. However, as it has been shown recently for SARS-CoV-2, MB antiviral activity occurs even in the absence of light [37,39]. Unlike MB, Zn-PcChol 8+ in the concentrations up to 5 µM has no antiviral activity without irradiation [35]. The mechanism of the "dark" MB antiviral effect on the infectivity of SARS-CoV-2 is elusive. One of the possible explanations comes from our simulations. We found that MB molecules are able to bind to the internal cavity inside the SARS-CoV-2 S-protein head (the inset in Figure 3), which opens for binding as the result of RBD transition into the "open" conformation. The binding of MB to the RBD domain may affect interaction of the S-protein with the ACE2 receptors on the host cell, thus, reducing virus infectivity. In contrast to MB, the much larger Zn-PcChol 8+ molecules did not appear in this cavity in our simulations. This agrees with our finding that Zn-PcChol 8+ lacks dark antiviral activity.
We also analyzed the binding of positively charged dyes to other S-protein key sites, such as the cleavage site and FP. It turned out that in all studied S-proteins in the loci between RBD and FP, Zn-PcChol 8+ and MB do not bind directly to the cleavage sites, apparently due to the presence of positively charged arginine residues in these sites. However, we cannot completely exclude the effect of cationic dyes on S-protein cleavage in SARS-CoV-2 and MERS-CoV, since one can see several binding sites between RBD and FP (Figure 4), although not directly at the cleavage sites. In SARS-CoV, subdomains between RBD and FP do not bind Zn-PcChol 8+ and MB, so, it is unlikely that cationic dyes would affect S-protein cleavage in this coronavirus. The facts that Zn-PcChol 8+ and MB bind to the region just after FP in SARS-CoV-2 and that the notable Zn-PcChol 8+ binding site in MERS-CoV S-protein appears just before FP, as well as the localization of the main binding sites adjacent to HR2 (Figure 4) mean that cationic dyes can potentially affect the process of fusion, but this hypothesis needs experimental proof.
The ability to recognize a variety of host cell receptors greatly impacts virus infectivity. Along with RBD binding to main protein receptors, NTDs of MERS-CoV can accommodate sialosides at neutral pH [56], while NTD of SARS-CoV uses only ACE2 receptor to bind to host cells [57]. In coronaviruses, NTDs of S-proteins have β-sandwich structure modified by various variable loops that affect recognition and binding to the host cell [25]. NTDs of SARS-CoV and SARS-CoV-2 show much lower (53.5%) homology compared to full-length S-proteins [58]. The total charge of the NTDs also varies among the studied coronaviruses, from −4 in SARS-CoV to +2 in SARS-CoV-2. According to the aligned NTD sequences, the key negatively charged residues that make up the prominent region on the head of the SARS-CoV S-protein in both SARS-CoV-2 and MERS-CoV are largely replaced by neutral amino acid residues. The substitution of the negative Asp23 in SARS-CoV with Arg21 or Lys27 in SARS-CoV-2 and MERS-CoV, respectively, is especially intriguing. The heterogeneous distribution of positively and negatively charged amino acid residues in the NTDs of SARS-CoV-2 and MERS-CoV results in the lack of the noticeable areas of negative electrostatic potential. Both Zn-PcChol 8+ and MB were rare guests in the SARS-CoV-2 NTD and did not bind to MERS-CoV NTD. In contrast, in each of the NTDs of SARS-CoV S-protein trimer, negatively charged amino acid residues form condensed areas of negative potential, attracting Zn-PcChol 8+ and MB molecules (Figures 2a and 3a). At the same time, the presence of such negatively charged regions can prevent SARS-CoV NTD contacts with negatively charged cellular receptors. This is in line with the experimental fact that NTD of SARS-CoV does not bind to sialic acid [57,58].
Finally, we found several PS binding sites in different S-protein domains involved in the initial processes of viral replication cycle-the primary recognition of cell receptors, binding to the main receptor, cleavage, and membrane fusion. The domain structure determines the selective binding of dye molecules and, thus, the nearby targets that can be damaged with singlet oxygen when illuminating the bound PS. In all three studied S-proteins, for both cationic dyes with photodynamic activity, Zn-PcChol 8+ and MB, we found a major binding site at the connection of the S-protein stalk and the head adjacent to the HR2 domain. The existence of such main binding site for cationic PSs, common to the S-proteins of SARS-CoV, SARS-CoV-2, and MERS-CoV, creates prospects for the wide use of PSs as photodesinfectants to combat the spread of coronaviruses. Our results can be useful in studying the initial electrostatic interactions of other cationic antiseptics with surface proteins of coronaviruses [59]. Brownian dynamics can help in studying the interaction of potential cationic antiviral drugs with negatively charged receptors of host cells that inhibit viral binding. A rational choice between the affinity for the receptors of host cells and viral structures can significantly increase the effectiveness of antiviral drugs, the mechanism of action of which involves electrostatic interactions.