{'fraction_non_alphanumeric': 0.0345370978332239, 'fraction_numerical': 0.029711096520026264, 'mean_word_length': 4.692580826013829, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Molecular dynamics simulations of the flexibility and inhibition of SARS-CoV-2 NSP 13 helicase

The helicase protein of the Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) is both a good potential drug target and very flexible. The flexibility, and therefore its function, could be reduced through knowledge of these motions and identification of allosteric pockets. Using molecular dynamics simulations with enhanced sampling, we determined key modes of motion and sites on the protein that are at the interface between flexible domains of the proteins. We developed an approach to map the principal components of motion onto the surface of a potential binding pocket to help in the identification of allosteric sites.

SARS-CoV-2 NSP13 helicase has 596 residues and is 99% homologous to SARS-CoV helicase. The X-ray structure from SARS-CoV (PDB ID 6JYT) shows five domains [20]. The three domains, IA, 1B, and 2A make up a triangular base with the N-terminal zinc binding domain (ZBD) at the top, and a stalk region connecting the base to the ZBD (Fig. 1). Fig. 1 shows the SWISS-MODEL [21] predicted structure based on 6JYT. Recently, Newman, et al., determined the X-ray structure of SARS-CoV-2 NSP13 helicase, for the apo protein and in complex with phosphate, adenosine 5 ′ -triphosphate, ATP, and other ligands [22]. The catalytic function of the helicase is to separate double stranded nucleic acids; the protein hydrolyzes nucleoside triphospate molecules (most often ATP), a chemical reaction that ultimately produces the mechanical energy necessary for the helicase to change conformations every step in the unwinding process [23]. The flexibility of the protein, key to its function, might make helicases a challenge for structural-based drug design [24], although a number of groups [25][26][27][28][29], have used computational docking and molecular dynamics (MD) simulations to identify potential SARS-CoV-2 inhibitors. The inhibitors are designed to bind to the ATP binding site of helicase [25][26][27][28], located between domains 1A, 1B, and 2A [25,30]. Another study designed inhibitors to bind to a potential allosteric site, involving resides of the ZBD, the stalk, and domain 1A [29]. Beck et al., used machine learning to find potential inhibitors without using structural information [31]. White, Lin, and Cheng demonstrated that two of their inhibitors inhibited helicase activity [26]. One simulation study, as part of a study of many SARS-Cov-2 proteins, using Markov state models and distributed computing, found a large domain motion [32].
For the SARS-CoV inhibitors [12][13][14][15][16][17][18][19] there is only a little known about the interactions between the molecules and the protein. The inhibitors bananin and its derivatives of Tanner and co-workers show a non-competitive mechanism [12] and mutagenesis experiments and molecular docking by Wang and co-workers suggest that the molecules interact with residues 257 and 259 on domain 1B [33]. The inhibitors of Adedeji and co-workers are also believed to be non-competitive, perhaps by inducing a structural change, although the binding site is not known [17]. Bismuth complexes may act as inhibitors by competing for the zinc binding sites [13]. Yu and co-workers reported inhibitors that modelling suggests bind to the ATP site [16]. For the other inhibitors, no binding sites have been proposed [14,15,18,19].
Inhibitor design could benefit from a better understanding of the flexibility of helicase, both with and without an inhibitor present. MD simulations, used to explore conformational space, can be aided using enhanced sampling methods. One general purpose, widely-used enhanced sampling approach is replica exchange [34,35]. In replica exchange (RE), simulations at a high temperature, which can more easily overcome high energy barriers, are run together with the simulations at the temperature of interest, with Monte Carlo swaps taken between the replicas. A limitation of RE is poor system size scaling. As the system size increased, the number of replicas required increases, making RE inefficient for large systems [36]. Both the utility and the limitations of RE are evident in a recent simulation of the SARS-CoV-2 spike protein [37]. This study required 46 replicas to span a temperature range of 310-350 K (and that was with many of the atoms restrained, effectively reducing the degrees-of-freedom of the system). Despite this, RE is a powerful method. Periole and Mark showed RE was 30-50 times more efficient than conventional MD in determining the melting curve for a small peptide in explicit solvent [38]. One method to increase the efficiency of RE is the Replica Exchange with Dynamical Scaling (REDS), which can reduce the number of replicas required by a factor of 8-10 [39,40].
In this work, we use computer simulations with the REDS method to understand the flexibility of SARS-CoV-2 NSP13 helicase and to predict alternative structures that the protein adopts. Simulations will be carried out for both the apo protein and the protein in complex with a potential inhibitor in the ATP site.

Methods
Examining the conformational flexibility of helicase presents a sampling challenge, as the time scales for the structural changes can be longer than typical simulation times, and a force field challenge, as the models must correctly weigh different, nearly energetically equal conformations [41]. Considerable effort has gone into the development of force fields for proteins [42,43]. The recently developed a99SB-disp model [44] has been shown to accurately reproduce the structure of a range of folded and disordered proteins [44], as well as a protein that contains both folded and disordered domains [45]. The a99SB-disp model is the result of many revisions of the Amber 99SB model [46], which was a revision of the Amber 94 force field [47]. The SB model revised the backbone φ/ψ dihedral terms, which corrected errors in the backbone secondary structure, particularly for glycine [46]. Best and Hummer revised the ψ dihedral terms to give the correct helix/coil populations, creating Amber 99SB* [48]. The side chain torsion potentials for four amino acids (isoleucine, leucine, aspartic acid, and asparagine) that showed disagreement with experimental data were reoptimized against quantum calculations, resulting in the a99SB*-ILDN model, also referred to as ff99SB*-ILDN [49]. To fix the helix propensities of charged residues, new charges from a restrained electrostatic potential (RESP) fit with constraints on the backbone charges to match neutral side chains, were assigned to charges for charged side chains, to give the a99SB*-ILDN-Q model [50]. The a99SB-disp model makes a number of changes to the a99SB*-ILDN-Q model, starting with enhanced dispersion interactions, through the Lennard-Jones potential, between the amide hydrogen and carbonyl oxygen backbone atoms. Side chain torsional terms were optimized against ab initio and experimental data. Side chain charges for the charged residues aspartate, glutamate, and arginine were re-parameterized, as was backbone torsional terms for glycine. Lastly, a new water model, a99SB-disp water with slightly stronger Lennard-Jones dispersion was developed to be used with the a99SB-disp model [44].
For the three zinc ions in the ZBD and the complexing cysteine and histidine residues, we used the Zinc Amber Force Field (ZAFF) [51]. The ZAFF force field gave different charges for the backbone atoms for cysteine residues. To make the ZAFF force field compatible with a99SB-disp and have the same charges for the N,H, C, and O backbone atoms, we modified the charges of the Cα and Hα atoms for the cysteine residues complexing a zinc ion, by adding 0.0077e to both atoms, so that charge is conserved. The rest of the side chain charges were from the ZAFF force field. The backbone charges for the histidine residues from ZAFF were the same as a99SB-disp.
Simulations were carried out both for the apo enzyme and with a potential inhibitor in the ATP binding site. This inhibitor was identified through a combination of docking using the FRED and HYBRID programs [52], and molecular dynamics simulations, as described in the Supplementary Information. The Generalized Amber Force Field (GAFF) [53] parameters were used for the potential inhibitor, using AM1-BCC charges [54]. The apo simulation used 36002 and the complex simulation used 35993 water molecules, using the a99SB-disp model, in a rhombic dodecahedron box. Charge neutrality was maintained by adding 9 chloride ions for the apo and 10 chloride ions for the complex, using the a99SB-disp force field for ions.
The simulations were setup using the tleap and antechamber programs of AMBER 16 [55,56]. Those coordinate and topology files were converted to GROMACS input files using the ACPYPE program [57]. Parameter sets for the a99SB-disp model were downloaded from GitHub at https://github.com/paulrobustelli/Force-Fields.
Configurational space sampling is enhanced using the REDS method [39,40], which combines conventional replicas with replicas that have a modified potential energy. The modified, or scaled replicas, can span a large range of energy and are placed between two conventional replicas, so that fewer replicas are required. Replica exchange works by having a replica at a temperature high enough to overcome energetic barriers. Choosing the high temperature is important to the efficiency of replica exchange [58][59][60][61]. An optimal high temperature of 366 K was found by Rosta and Hummer for a small protein in water (the l-repressor fragment  [61] and Heo and Feig used a temperature of 360 K for the prediction of SARS-CoV-2 protein structures (but not including helicase) [62], suggesting that 360 K is a reasonable high temperature. Replica exchange simulations for eight SARS-CoV-2 proteins, but not helicase, used a temperature range of 310 K-350 K, with a number of replicas from 20 to 60 [63]. Our REDS simulations had conventional replicas, using the plain, unscaled replicas, at temperatures of 310, 319, 331, 345, and 360 K, with scaled replicas at 314, 325, 338, and 352 K, for a total of 9 replicas. These temperatures were determined using the method of Patriksson and van der Spoel [64]. RE would require 42 replicas for this system. The simulations used the REDS method as implemented in the GROMACS simulation package [65]. Details of this implementation are given elsewhere [40] and the additions to the GROMACS source code were downloaded from Gitlab at https://gitlab.com/csumma_uno/grom acs.
Initial coordinates of the protein were taken from the SWISS-MODEL [21] predicted structure (P0DTD1 PRO_0000449630, Model 1 [66]) based on the 6JYT [20] structure from SARS-CoV. The helicases from the two viruses have a close homology, differing only by one residue (570 which is an ILE in SARS-CoV and VAL in SARS-CoV-2). After our simulations were set-up and running, crystal structures of SARS-CoV-2 NSP13 were published [22]. The apo structure (PDB ID 7NIO) has two chains, with a Cα root-mean-square deviation (RMSD) of 0.15 nm. The SWISS-MODEL structure has a Cα RMSD of 0.20 nm with Chain A and 0.19 nm with Chain E of the 7NIO structure. Both Chain A and Chain E have unresolved residues, mostly in domain 1B. Also after our simulations were running, White, et al., published a paper pointing out problems with the 6JYT structure, including cis peptide bonds in the loop near the ATP binding site [26]. The 6JYT structure has two chains. The SWISS-MODEL structures used Chain B, which does not have cis peptide bonds in this loop. Comparisons between the SWISS-MODEL and the 7NIO structures are shown in Fig. S2. The loop does have a different structure from the 7NIO structures, but the residues that have important interactions with the ligand (LYS 288, ASP 374, and GLU 375) are in similar positions. The SWISS-MODEL structure was protonated using the H++ server [67]. The structure was solvated and chloride ions were added to give a neutral system. This structure was minimized using steepest decent for 1000 steps, followed by constant molecule number, volume, and temperature (NVT) simulations for 100 ps and then 100 ps at constant N, pressure, and T (NPT) for each of the temperatures.
The simulations used a 0.002 ps time step. All bonds involving hydrogens were constrained using the LINCS [68] and Settle [69] algorithms. Long-ranged electrostatics were treated using Particle mesh Ewald [70] with a Fourier mesh spacing of 0.16 nm, a real-space cut-off of 1.2 nm, and fourth-order interpolation. Lennard-Jones interactions used a 1.2 nm cut-off, with long-ranged dispersion corrections for energy and pressure. Temperature was maintained using stochastic velocity rescaling [71], with separate protein and solvent thermostating and a relaxation time of 0.1 ps. Pressure was maintained, at 1 bar, using a Parrinello-Rahman barostat [72] with a relaxation time of 2 ps and a compressibility of 4.5x10 − 5 bar -1 . Exchanges between replicas were attempted between adjacent replicas every 100 time steps. The simulations were run for a total of 150 ns. The last 100 ns were used as production runs for data analysis. Root mean square deviations and root mean square fluctuations were determined using gmx rms and rmsf utilities in Gromacs [73], respectively. Clustering of structures based on the root-mean-square-deviation (RMSD) with a 0.25 nm cut-off used the gmx clustering utility in Gromacs [73]. The solvent accessible surface area (SASA) was calculated with using the sasa utility in Gromacs [73], with a probe radius of 0.14 nm. Principal component analysis and correlation matrix calculations used the gmx covar utility in Gromacs [73]. The correlation matrix is defined as where the angle brackets donate an ensemble average and Δr i is the distance of atom i from its average position [74]. Using a structure that resulted from the lengthy MD simulations, additional binding pockets were explored using PrankWeb [75], as web interface for the P2Rank ligand binding site prediction method [76]. Docking into these sites was done using the program Autodock Vina [77], using 3868 compounds from the SWEETLEAD [78] library. The SWEETLEAD library consists of approved pharmaceuticals, illicit drugs, and traditional medicinal herbs. Spatial data file (SDF) input files were converted to PDBQT (protein data bank format plus partial charges and atom types) input for Autodock using OpenBabel [79].

Results
The top scoring compounds from the docking analysis were simulated for 100 ns, of which three stayed in the ATP site for the entire simulation length. Full details are given in the Supplementary Information. The top scoring inhibitor FCID1 was selected for further simulation using enhanced sampling. The initial placement of FCID1 is shown in Fig. 2.
Molecular simulation results. For the simulation of the complex, some of the replicas give trajectories where the FCID1 molecule leaves the binding pocket. For the analysis of the results of the complex that follows, only the bound structures are used. Structures are saved every 0.01 ns, so the entire 150 ns simulation has 15,000 saved structures. Of those 11,373 have FCID1 in the ATP site. In two of the nine replicas, the inhibitor will eventually leave the binding site. In the plots below showing the time evolution of selected properties, the absence of the data corresponding to the unbound data creates small gaps.
The structures for both the apo protein and the complex both change considerably from the initial structure (Fig. 3). The use of replica exchange results in sudden in changes in the root-mean-square deviation (RMSD) as configurations from different temperatures are swapped. Structures with an RMSD over 1 nm are observed. The RMSDs from the initial structure for the five domains, average over the last 100 ns of the simulations, are given in Table 1. The ZBD and the stalk domains are the most rigid, as indicated by the small RMSD and the small standard deviation of the RMSD. The 1A, 1B, and 2A domains show large structural changes, showing that the conformational changes involve not only motion of the domains relative to each other but structural changes within these three domains. Domain 2A has the largest RMSD. The apo and the complex show similar structural changes.
The relative motion of the domains is determined by calculating the distance between the center of masses of the domains, using the Cα atoms (Fig. 4). The ZBD and stalk domains do not show any significant inter-domain flexibility, so that data is not shown. The distance between domains 1A and 1B is relatively constant over the simulation (Fig. 4A), For the apo-protein the 1A-1B distance tends to decrease with time, although it does sample some configurations near the end of time simulation where that distance is larger. The distance between domains 1A and 2A shows a larger change (Fig. 4C), fluctuating over distances of almost 1 nm. From the distribution of these distances (Fig. 4D), it can be seen that the complex fluctuates between two conformations, one with a 1A-2A distance around 3.1 nm and another around 3.7 nm. The apo protein shows a single broad distribution. The 1B-2A domains show the largest change in distance (Fig. 4E) The residues in the ATP binding site, particularly residues ASP374 and GLU375 have small RMSF values, around 0.2 nm. Other rigid regions are around residue 136, near the boundary between the Stalk and 1B domains, around residue PRO234, and around residue ASN388.
The structures from the last 100 ns of the simulations at 310K were clustered based on RMSD. The positions of the Cα atoms for residues 1-590, excluding the highly flexible C-terminus, were used. This resulted in 20 clusters for the apo protein and 11 for the complex. There is no overlap between the structures from the two separate simulations; no clusters contained structures form both the apo and complex simulations. Selected properties of the top 5 clusters are given in Table 2. The    top 5 clusters contained about 90% of the total structures. Also shown are properties for the two clusters for the apo protein that have the largest RMSD from the initial structure. For comparison, the properties of the SWISS-MODEL structure are given. For the apo structures, the most populated cluster contains the structures that are predominant over the last 20 ns of the simulation, suggesting that this is the structure that is most stable. For the complex, the structures that are prevalent in the last 20 ns are in clusters 1, 3 and 5. The solvent accessible surface area for the structures from the simulations are larger than the initial structure, indicating, along with the domain distances, that the simulation structures are more open. Fig. 6 shows the structures of selected clusters, including those with the biggest populations and one with the largest RSMD from the original structure (cluster 9). A significant structural change involves motion of domain 1B (shown in blue). This domain moves between a structure closer to domain 2A (in orange) and a structure closer to the ZBD domain (in red). The motion involves hinge regions around residue TYR149, at the end of the stalk domain, and VAL232. Both of these regions have small fluctuations, bracketing a region of high fluctuations (Fig. 6). Domain 2A moves by bending around residue GLY439. The Markov state simulations using distributed computing found a large motion between domains, which they identified as domains 1A and 2A [32]. Upon inspection of Fig. S10 from that publication, the authors have mis-identified domain 1B as domain 1A, so that the large motion is actually between domains 1B and 2A.
The FCID1 inhibitor interacts with residues on both the 1A and 2A domains, as initially docked (Fig. 2). As the distances between the two domains increases, by over 1 nm, the molecule stays with domain 1A, where it has the strongest interactions, but loses interactions with residues on domain 2A. Fig. 7 shows the distances between the FCID1 atoms and the protein which are present in the original structure. There is a strong interaction between the ASP374 side chain and the cyclopentyl ammonium nitrogen (identified as N20), which is maintained during the length of the simulation (Fig. 8A). Another strong interaction is between the GLU375 side chain and two nitrogens on the inhibitor, the amide nitrogen (N18) and pyrrole ring nitrogen (N19). This interaction is also maintained during the simulation. Fig. 7B shows the N18-E375 distance; the N19-E375 distance is fairly constant as well. The Nζ atom of LYS288 interacts with the oxygen atom on the inhibitor. This interaction is not constantly maintained during the simulation (Fig. 8B). In cluster 1, the LYS288 Nζ and O atoms are in direct contact. In clusters 2-5, the atoms are over 0.5 nm apart, separated by solvent. The GLY538 Cα atom, part of domain 2A, has a close contact with a carbon atom (C3) on the cyclopentyl group. This interaction is lost for many of the simulation structures, including in clusters 1, 3 and 4, as domain 2A moves away from domain 1A (Fig. 9A).
The correlation of motion between different parts of the protein as given by the correlation matrix is shown in Fig. 10 Table 2 Cluster rank, population fraction, RMSD from initial structure, distances between domains, and solvent accessible surface area (SASA). The large scale motion of the protein involves movement of domain 2A, relative to the other domains, including the neighboring domain 1A. A possible means to suppress the flexibility of the helicase might be designing an inhibitor to bind at the interface between domains 1A and 2A, in a region away from the both the ATP binding site and the nucleic acid binding channel. The nucleic acid binding channel is located between domain 1B on one side and domains 1A and 2A on the other [80]. Using the predominant structure of the apo protein from the clustering analysis, potential drug binding sites were identified using P2Rank [76]. Four, out of 17, binding sites were found at the interface between domains 1A and 2A, as shown in Fig. 10A. (Pockets 3 and 4 are close together but are separated by a narrow channel created by residues SER264, ASN268, and LYS462.) The pockets are numbered according to their rankings from PrRank, so that among the four sites in the region between domain 1A and 2A, pocket 1 is the highest ranked and pocket 4 is the lowest rank. Pocket 1 is the ATP binding site. In order to inhibit motion of domain 2A, the inhibitor would ideally be at a site with residues that are involved in domain motion but moving in opposite directions, as near a hinge or shear motion. The molecule would possibly inhibit the motion by interacting with the residues that need to move in opposite directions. The four binding sites can be assessed based on the motions of the residues in those binding sites using the results of the PCA. P2Rank returns the points of the solvent accessible surface (SAS) for each binding site. For a point j on the surface, the nearest residue is  found, which is labeled i j . The motion of that residue can be characterized from q α (i j ) the eigenvector of principal component α for residue i j .
The total motion of the binding site can be found from where N is the number of SAS points for a specific pocket. The angle between the direction of motion of one point on the surface relative the average direction is a way to determine if a principal component of motion is moving the whole surface together, or if parts of the surface move in different directions. The cosine of this angle, given by can be mapped onto the SAS. Of the four pockets, site 4 shows the largest spread d α (i) values, or most anisotropy, for the first PC, see Fig. 10 (B). Areas of the surface that move in the same direction as the average (positive d α (i)) of the pocket are shown in blue and regions that move in the opposite direction as shown in red. Pocket 4 has contact with SER264, ALA267, and ASN268 on domain 1A and LYS462 on domain 2A, which move in opposite directions in PC1. The pocket also makes contact with PHE437 and LYS460. Pocket 4 shows the most anisotropy for the second PC as well, although pockets 1 and 3 also show a high anisotropy, as shown in Fig. 10 (C). In pocket 4, LYS460 and LYS462 move in opposite directions from ASN 268 for PC2. Pocket 3 has contact with ASN265 and ARG442, which move in opposite directions. Pocket 1, the ATP site, the residues GLY287, PHE291, GLY400, and GLN537 move in opposition to residues SER377, MET378, and ASP 274. Pocket 2, in contact with GLN281, ASP401, ALA403, and MET429, has a small anisotropy for both the first and second PC, meaning all those residues move in the same direction in those modes.
The PC map analysis suggests that pockets 3 and 4 might be good targets for a non-competitive inhibitor, which could limit the flexibility of domain 2A. If an inhibitor could simultaneously make contact with residues among SER264, ALA267, and ASN268 on domain 1A and LYS460 and LYS462 on domain 2A, it could reduce the movement of both PC1 and PC2. In pocket 3, an inhibitor in contact with ASN265 and ARG442 might reduce movement of PC2. From the correlation matrix residues 267 and 268 have a negative correlation (around − 0.1) with Pocket 4, while promising in terms of inhibiting domain motion, is not ranked highly by PrRank (13th out of 17 potential pockets). The druggability of the pockets was assessed through virtual screening. Docking into pocket 4 gave weakly binding compounds; in contrast, docking into pocket 3, which is ranked fourth out of the 17 pockets, gave more promising scores (details are in the Supplementary Information). A top scoring compound made hydrogen bonds with ASN265, THR440, and ARG442, as well as hydrophobic contacts between benzyl group and ALA446, VAL449, and the methyl group on THR440. It therefore makes significant contacts with residues on both domains 1A and 2A.

Conclusions
The simulations revealed that the NSP13 helicase, both with and without an inhibitor at the ATP site, is highly flexible, with an RMSD from the original structure of over 0.6 nm (see Fig. 3 and Table 1). Some structures have RMSD values over 1.3 nm ( Table 1). The largest structural changes involve motion of domain 2A, relative to the other domains. Principal component (PC) analysis showed that about 43% of the fluctuations of the protein are captured in the first two principal components, both of which involve large amplitude motion of domain 2A. The flexibility of the helicase might be altered with an inhibitor to bind at the interface between domains 1A and 2A, in a region away from the ATP and the nucleic acid active sites. Four potential drug binding pockets between domains 1A and 2A were identified. We sought a binding site that could change the flexibility so we developed an approach that maps the motion of the atoms on the surface of each pocket onto the pocket's surface, as shown in Fig. 10. Maps were done using the motion as characterized either by the first or second PC. If a pocket is at the boundary between regions that tend to move in different directions, then the PC map would show negative and positive areas on the surface. The combination of the suitability of the pocket for drug binding and the indication from the PC map that it is involved domain motion, one pocket was identified as a promising allosteric site.
A number of studies have compared PCA from ligand bound and ligand free simulations, including those for helicase [81,82], to determine how does a ligand changes the principal components. We are looking at the reverse question, how does PCA point to a site where a ligand might influence the dynamics, and hence the PCAs. Other work has combined PCA and binding pocket characterization. The GRID/-CPCA method generates data for PCA for different ligand types in a pocket [83]. Other methods perform PCA on binding pocket grids from multiple MD structures and to analyze how a pocket changes shape [84][85][86][87]. The aim of these methods is to guide structure-based drug design and not necessarily to identify allosteric sites. Our approach differs in that it maps the relative motion of pocket surface atoms as determined from PCA onto the surface.
The proposed allosteric site is in a different location on the protein than some of the previously proposed allosteric sites for SARS-Cov and SARS-Cov-2 NSP13 helicase [29,33]. X-ray crystallographic fragment screening found three possible allosteric sites: between the stalk and ZBD domains, between the stalk and 1B domains, and between domains 1A and 2A [22]. The fragments in this third site (as determined from the 5RMB and 5RMF structures) makes contacts with residues GLN281, LYS430, PRO434, and PHY 437 on domain 1A and residues TYR457, ASP458, and LYS460 on domain 2A. These residues are in a similar region to one of the binding pockets identified from our simulations (pocket 3). Our results from the PCA analysis support the conclusion from Newman et al., [22] that this region could be a good target for allosteric inhibitor design.
Helicase NSP13 is part of a replication and transcription complex (RTC) together with other proteins, including NSP7, NSP8, and NSP12, as well as another copy of itself [80]. The helicase structure in that complex is likely to be different that it is in solution. However, in silico design has led to drugs that inhibit helicase ATPase activity in vitro. [26] Further studies would be needed to determine if molecules designed to bind to the allosteric site are effective inhibitors of helicase. While many structural details into the exact mechanism of replication remain yet to be discovers, the relevance of the range of motion of domain 2A observed in our studies provides some insight into how the helicase carries out its crucial role in the RTC.

Data and software accessibility
The REDS modifications to the GROMACS simulation package are available from Gitlab at https://gitlab.com/csumma_uno/gromacs. Input files for the simulations and structures from the cluster analysis are provided in associated content.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.