Conformational ensemble of the full-length SARS-CoV-2 nucleocapsid (N) protein based on molecular simulations and SAXS data

The nucleocapsid protein of the SARS-CoV-2 virus comprises two RNA-binding domains and three regions that are intrinsically disordered. While the structures of the RNA-binding domains have been solved using protein crystallography and NMR, current knowledge of the conformations of the full-length nucleocapsid protein is rather limited. To fill in this knowledge gap, we combined coarse-grained molecular simulations with data from small-angle X-ray scattering (SAXS) experiments using the ensemble refinement of SAXS (EROS) method. Our results show that the dimer of the full-length nucleocapsid protein exhibits large conformational fluctuations with its radius of gyration ranging from about 4 to 8 nm. The RNA-binding domains do not make direct contacts. The disordered region that links these two domains comprises a hydrophobic α-helix which makes frequent and nonspecific contacts with the RNA-binding domains. Each of the intrinsically disordered regions adopts conformations that are locally compact, yet on average, much more extended than Gaussian chains of equivalent lengths. We offer a detailed picture of the conformational ensemble of the nucleocapsid protein dimer under near-physiological conditions, which will be important for understanding the nucleocapsid assembly process.


Introduction
Recent advances in cryo-EM microscopy have allowed us to structurally characterize very large proteins and their complexes [1,2], and structure prediction by AlphaFold has reached previously unseen accuracy [3]. Together with classical macromolecular crystallography and protein NMR, these advances could create an impression that mostif not allprotein structures can be rather easily solved or accurately predicted. Yet one class of proteins is outstanding in this respect, the intrinsically disordered proteins (IDPs) -especially proteins where wellfolded domains are connected by long, flexible, intrinsically disordered polypeptide segments. Such proteins are usually too big for NMR analysis and too flexible for crystallographic or cryo-EM analysis [4].
SARS-CoV-2 encodes for more than 20 proteins including 16 nonstructural proteins (nsp1-16), spike glycoprotein (S), nucleocapsid protein (N), membrane protein (M), envelope protein (E) and several accessory proteins [5]. Most of these are well folded and were structurally characterized using cryo-EM or macromolecular crystallography. However, the N protein is composed of two globular domains and three intrinsically disordered regions (IDR1-3, Fig. 1). The structured domains are named NTD (N-terminal domain) and CTD (C-terminal domain). The NTD has a right hand-like fold composed of a β-sheet core with a large protruding central loop that resembles a finger ( Fig. 1) [6]. Importantly, it binds both single-stranded and double-stranded RNA [7]. The CTD also binds RNA with a preference for viral genomic intergenic transcriptional regulatory sequences [8]. It is composed of eight helices and two β-strands and its shape reminds of the letter C ( Fig. 1) [8,9]. Two C-shaped monomers form a dimer which is responsible for the dimerization of the N protein [8][9][10][11]. As outlined above, SARS-CoV-2 nucleocapsid NTD and CTD are well structurally characterized. However, the current knowledge about conformations of the full-length N protein is more limited. Importantly, small angle X-ray scattering (SAXS) experiments have shown that the full-length N protein in solution is a dimer with a radius of gyration of 5.9 nm [12].
SAXS is a powerful method to characterize macromolecules in solution [13,14]. Despite the loss of structural information due to orientational averaging of the scattering signal, SAXS provides useful information about the shapes and dimensions of macromolecules in solution. In particular, standard analysis of SAXS data permits the determination of low-resolution structures of proteins in the form of molecular envelopes, where protein shapes are represented by a set of dummy atoms [15]. However, this approach is only valid for proteins with a stable ternary structure that can be represented as a rigid body, which the N protein is certainly not. SAXS analysis of IDPs requires different approaches that are based on generating an ensemble of protein conformations [16][17][18][19][20][21][22]. Useful and appropriate tools for generating conformational ensembles of biomolecules are molecular dynamics (MD) simulations. However, despite steady developments in this field [23][24][25], all-atom MD simulations of multi-domain and partially disordered proteins are still challengingnot only because of the large sizes of such proteins but also because of the long time scales on which large conformational fluctuations occur. Coarse-grained molecular simulations, on the other hand, provide means to efficiently sample conformational ensembles of large, multi-domain and dynamic proteins and their complexes [4,20,26]. However, the predictive power of coarsegrained simulations certainly lags behind that of all-atom MD simulations. One way to overcome these challenges is to use experimental data to guide or bias coarse-grained simulations.
Recently, original data are more and more reused within the scientific community and many studies are published as open-access articles. Usually the data are deposited in open-access databases such as PDB (Protein Data Bank) or EMDB (Electron Microscopy Data Bank). SASBDB (Small Angle Scattering Biological Data Bank) and PED (Protein Ensemble Database) [27,28] are also available but not yet broadly established and their use is not mandatory (unlike PDB) by most journals. Nevertheless, professors Tengchuan Jin and Hongliang He kindly provided us with the SAXS data they recently acquired on the SARS-CoV-2 N protein [12]. Here, we combined coarse-grained simulations of the N protein dimer with the data from SAXS experiments using the ensemble refinement of SAXS (EROS) method [18]. Our results provide a detailed picture of the conformational ensemble of the SARS-CoV-2 N protein dimer.

Coarse-grained simulations
We used the coarse-grained model for multi-domain IDPs introduced by Kim and Hummer [29]. Replica exchange Monte Carlo (REMC) simulations of the coarse-grained model were performed with replicas at 24 different temperatures T i below and above the room temperature T = 300 K. Specifically, we used T/T i = 0.45, 0.49, 0.52, 0.55, …, 1.09, 1.12, 1.15. Following the simulation methods of Kim and Hummer [29], the basic MC steps were rotations and translations of each of the rigid domains. For the flexible linkers and termini, in addition to local MC moves on each of the residue beads, crank-shaft moves were employed to enhance sampling. The REMC simulation consisted of 1.1 × 10 8 MC cycles, where the initial 10 7 MC cycles were used for equilibration and the subsequent 10 8 MC cycles for data collection. The protein conformations were saved every 5000 MC cycles at room temperature, T i = T, which resulted in an ensemble of N = 20,000 conformations for further analysis. Selected conformations were visualized using VMD [30].

Ensemble refinement of SAXS
We computed scattering intensity profiles of the recorded conformations using an algorithm co-developed with the EROS method [18] taking the default value of the electron density of the protein hydration shell (0.03 e/Å 3 ). Then the ensemble-averaged scattering intensity profile was computed as where the index k = 1, …, N labels the conformations in the ensemble, w k is the statistical weight assigned to conformation k, I k (q) denotes the scattering intensity profile of the k-th conformation, and q denotes the momentum transfer. (We use the notation q = 4πsin(θ/2)/λ, where θ is the scattering angle and λ is the X-ray wavelength). The discrepancy between the experimental SAXS data, I exp (q), and the ensembleaveraged scattering intensity profile, I sim (q), was quantified by where i = 1, …, N q labels points in the SAXS dataset, q i are the values of At first, all conformations in the simulation ensemble were taken to be equally relevant, which implied w k = w k (0) = 1/N for k = 1, …, N and resulted in χ 2 = 1.9 ( Fig. 2A). The agreement between the simulation ensemble and the experimental SAXS data was thus good but still could be improved. We therefore refined the conformational ensemble using the EROS method [18], which prevents over-fitting and keeps the refined ensemble (with weights w k ) as close as possible to the reference ensemble (with weights w k (0) = 1/N). Following the original EROS method, we introduced a function F = χ 2 -θS, where θ is a parameter that expresses the confidence in the reference ensemble whereas is the relative entropy, which equals to the negative Kullback-Leiber divergence. We minimized F with respect to the set of weights w k using a steepest descent method within the log-weights formulation Their R g and D max values are indicated. E) Cartoon representation of a dimer structure that best fits the experimental SAXS data (χ 2 = 1.12). The colour code is as in

Contacts between residues
Each residue bead was assigned a van der Waals radius as introduced by Kim and Hummer [29]. We assumed that a pair of beads, i and j, with van der Waals radii σ i and σ j , respectively, were in contact if their distance r ij was smaller than 1.5 × σ ij , where σ ij = (σ i + σ j )/2 [33]. According to this criterion, we calculated a map of residue contacts for each conformation obtained in the REMC simulations. The relative weight of a given contact was then taken as a sum of statistical weights w k of the conformations in which this contact was present.

Results
To efficiently sample conformations of the nucleocapsid protein in solution, we used an implicit-solvent, coarse-grained model introduced by Kim and Hummer [29]. This model is equipped with a transferable energy function and devised to simulate conformational ensembles of large, dynamic, multi-domain proteins. It has been successfully applied to systems ranging from membrane-protein trafficking machineries [34,35] to lipid kinases in dynamic complexes with regulatory proteins [36,37] to hijacking of host protein by viruses [38,39] and to protein complexes controlling the biogenesis of autophagosomes [40]. In the framework of this model, amino acid residues are represented as spherical beads centered at the α-carbon atoms. The interactions between the residue beads are described by amino-acid dependent pair potentials and Debye-Hückel-type electrostatics. Folded protein domains are treated as rigid bodies whereas inter-domain linkers and flexible termini (IDR1-3, Fig. 1) are represented as polymers of amino acid beads with bending, stretching and torsion potentials. Here, the rigid domains were the NTD (PDB entry 6YI3) [7], the dimer of the CTD (PDB entry 7DE1) [8] and the hydrophobic α-helix within the polypeptide segment joining the NTD and CTD (PDB entry 7PKU) [41] (Fig. 1).
We performed extensive replica exchange Monte Carlo simulations of the coarse-grained model (see Materials and methods). From the simulations we obtained an ensemble of N = 20,000 structures of the N protein dimer at room temperature (T = 300 K). We computed the scattering intensity profile for each of the simulation structures using an algorithm co-developed with the EROS method [18]. The ensembleaveraged scattering intensity profile was found to be in good agreement (χ 2 = 1.9) with the SAXS data published by Zeng et al. [12] ( Fig. 2A). We also performed analogous simulations of the monomeric form of the N protein (Fig. 1B), and computed the scattering intensity profile corresponding to the ensemble of the simulation structures of the N protein monomer (Fig. 2B). This scattering intensity profile is distinctly inconsistent with the experimental SAXS data (χ 2 = 52), which confirms that monomers of the N protein are very unlikely to exist in solution at the concentrations used in the SAXS experiments. (To our knowledge, the N protein has never been reported to be in the monomeric form). This result illustrates how SAXS and computer simulations can be used together to determine the multimeric state of a macromolecule.
Next, we refined the conformational ensemble of the N protein dimer using the EROS approach, which employs the maximum-entropy method to minimally modify the statistical weights of the simulation structures in an attempt to match the conformational ensemble to the experimental SAXS data (see Materials and methods). The refined ensemble was very close to the original ensemble and in perfect agreement with the SAXS data (χ 2 = 1, Fig. 2C). We also employed the minimum ensemble method [34] and picked out four structures of the N Fig. 3. Analysis of the radius of gyration (R g ) and maximum dimension (D max ) of the N protein dimer. Histograms of A) R g and B) D max . The lines in gray and black correspond to the ensemble of simulation structures before and after refinement, respectively. C-E) Snapshots of selected conformations with C) R g = 5 nm D) R g = 6 nm and E) R g = 7 nm, as indicated by the arrows in panel A. The colour code is as in Fig. 1. protein dimer that jointly fit the SAXS data perfectly well (χ 2 = 1, Fig. 2D). This set of structures gives us a glimpse into conformational fluctuations of the N protein dimer. We also identified a single structure that best fits the SAXS data (χ 2 = 1.12). We note that this structure (Fig. 2E) involves very few inter-domain contacts, thus it cannot be stable in solution: it represents only a single "snapshot" of the flexible N dimer.
To quantify the degree of conformational fluctuations of the N protein dimer, we computed the radius of gyration (R g ) and the maximum dimension (D max ) of each of the simulation structures. The maximum dimension is defined as the maximum distance between any two points of the protein. The resulting histograms of R g and D max together with three representative structures are shown in Fig. 3. We note that in the simulation ensemble (i.e. before refinement), structures with an R g between 5.3 and 6.3 nm and a D max between about 18 and 20 nm are slightly over-represented, whereas more compact structures with an R g smaller than 5 nm and a D max smaller than 18 nm are slightly underrepresented. Taken together, the histograms in Fig. 3 show that the N protein dimer exhibits large conformational fluctuations with an R g ranging from about 4 to over 8 nm and a D max between 12 and 32 nm.
To characterize the structural ensemble of the dimer, we identified contacts between pairs of amino acid residues and the relative frequency of their occurrence. To this end, we determined contacts between pairs of amino-acid beads in each of the simulation structures using a distance criterion [33]. The relative frequency f of a given contact in the refined ensemble was then taken as a sum of statistical weights of the conformations in which this contact was present (see Materials and methods). The maps of the relative frequency of the intra-and interchain contacts are shown in Fig. 4A and B, respectively. The intrachain contacts are the contacts formed within each of the polypeptide chains of the dimer. The inter-chain contacts, on the other hand, are those between the two polypeptide chains of the dimer. The points in red indicate frequent contacts with f > 0.1. The points in orange correspond to transient contacts with 0.01 < f < 0.1. The points in yellow indicate rare contacts with 0.001 < f < 0.01. Finally, the points in gray show contacts within rigid domains that do not evolve in the course of the REMC simulations.
The contact map in Fig. 4A shows that the frequent contacts can be grouped in two categories: i) contacts formed within each of the three IDRs and ii) contacts formed by either the NTD or the CTD with the hydrophobic α-helix in the IDR2. The contact map in Fig. 4B shows that there are rather few contacts between the two chains and that these contacts are mostly rare or transient.
To characterize the level of extension or compactness of the IDRs, we determined their end-to-end distances in each of the simulation structures. Fig. 5 shows the resulting histograms of the end-to-end distances. The end-to-end distances of IDR1, IDR2 and IDR3 are found to vary from 1 to 9 nm, from 1 to 13 nm, and from 0.5 to 11.5 nm, respectively. The average and root-mean-square values of the end-to-end distances are summarized in Table 1.
It is interesting to compare the root-mean-square end-to-end distances of the IDRs with those of freely-jointed chains, which are also often called ideal chains or Gaussian chains. The root-mean-square endto-end distance of a Gaussian chain of n beads is equal to the square root of n times the length of the bond between consecutive beads in the chain. As in the coarse-grained model for multi-domain IDPs introduced by Kim and Hummer, here we take the chain beads to be centered at the α-carbon atoms and the length of the pseudo-bond between consecutive α-carbon atoms to be 0.38 nm. The results presented in Table 1 show clearly that each of the three IDRs exhibits a root-mean-square end-toend distance larger than the Gaussian chains of the corresponding length. This observation implies that although the conformations of the IDRs are locally compact (see the contact map in Fig. 4A), they are clearly more extended than those of the Gaussian chains overall. Interestingly, the N-and C-terminal tails (IDR1 and IDR3) behave somewhat more like Gaussian chains compared to the middle linker (IDR2).

Discussion
The thermodynamic state of an IDP is an ensemble of rapidly interconverting conformations. Important examples of IDPs are large, multi-domain proteins in which several autonomously folded domains are held together by intrinsically disordered regions [42]. Their structural analysis is difficult because of their large sizes and dynamic nature [4]. A notable example is the nucleocapsid protein of SARS-CoV-2 (Fig. 1). The question we posed in this study was about the degree of conformational fluctuations of this flexible protein complex.
We performed extensive REMC simulations of the N protein dimer and refined the resulting ensemble of conformations using the available SAXS data (Fig. 2). Our results show that the N protein dimer exhibits large conformational fluctuations in solution. Its radius of gyration varies between about 4 and 8 nm with an average of 5.9 nm (Fig. 3A). The maximum dimension exhibits large fluctuations in a range of 12 to 32 nm (Fig. 3B). The NTD does not make direct contacts with the CTD (Fig. 4); however, the hydrophobic α-helix within the IDR2 makes frequent and transient contacts with both the NTD and the CTD. Frequent contacts within each of the three IDRs indicate that these regions can adopt conformations that are locally compact. Yet, our analysis of the distances between termini of the IDRs (Fig. 5 and Table 1) shows that each of the three IDRs exhibits conformations that are more extended than Gaussian chains of corresponding lengths. Taken together, our results provide a detailed picture of the conformational ensemble of the SARS-CoV-2 N protein dimer under near-physiological conditions.
The structure and flexibility of the SARS-CoV-2 N protein is likely to be important for the assembly of the nucleocapsid. Interestingly, both the NTD and CTD can bind RNA, yet it is evident from the contact maps in Fig. 4 that there are practically no direct contacts between the NTD and the CTD. Thus our structural analysis suggests that the NTD and CTD do not directly cooperate in RNA binding. However, the hydrophobic α-helix within the IDR2 makes multiple intra-and inter-chain contacts with both the NTD and CTD. This α-helix makes most frequent contacts with the NTD of the same chain (Fig. 4A) and least frequent contacts with the NTD of the other chain in the N dimer (Fig. 4B). Frequent contacts are observed also within each of the disordered segments IDR1-3. The presence of these contacts indicates that the IDRs adopt conformations that are locally compact. It remains to be shown how IDRs influence other interesting properties of the N protein such as its ability to phase separate with RNA [43][44][45][46][47]. However, we note that the conformations adopted by the N protein could be modified by RNA or protein binding [41,48].

Author contributions
EB and BR planned the research, BR performed the simulations and analyzed the results, EB and BR wrote the manuscript.

Declaration of Competing Interest
The authors declare no conflicts of interests.