Protein glycosylation is an essential aspect of post-translational modification, with 50-70% 36 of human proteins been estimated to be glycosylated to some degree ). These 37 glycans play significant roles in numerous biological processes, such as cellular signaling, 38 recognition and adhesion, protein folding, structural stability, and immune system 39 interactions (Dennis et al., 1999; Dwek, 1996 The HIV-1 Envelope Glycoprotein (Env) is a heterodimeric trimer composed of proteins 55 gp120 and gp41, and is responsible for the molecular recognition of the host receptor and fusion 56 into host target cells. A number of obstacles hinder traditional vaccine design methods in case of 57 HIV-1 -namely its remarkable sequence diversity, conformational plasticity, dramatic shifts in 58 position and number of glycans (Zhang et al., 2004) , and extremely dense glycosylation making up 59 to approximately half the mass of the entire Env molecule (Lasky et al., 1986) . As a result, there 60 has been only limited success in eliciting broadly neutralizing antibodies (bNAbs) to Env vaccine 61 immunogens (Bradley et landscape of the glycoprotein, long and often multiple trajectories need to be run, preferably with 99 enhanced sampling techniques (Yang et al., 2017) . 100
We had previously established a high throughput pipeline to robustly build atomistic 101 models of glycans, sample the glycan conformational space, construct the glycan network 102 topology, and extract molecular level description of the viral glycan shield (Berndsen et al., 2019) . 103
We had also set up a protocol to experimentally validate the models by generating synthetic cryo-104 EM datasets and comparing the resulting maps in a quantitative manner. Here, we have employed 105 a similar quantitative comparison using computational and experimental single-particle cryo-EM 106 data and show that the dynamics of high-mannose and native glycosylation remains practically 107 indistinguishable around the glycan stem b-mannose (BMA) residue. This provides further 108 rationale towards the need for computational modeling of complete glycans in a natively 109 glycosylated shield. Presently this seems to be the only way to elucidate the ramifications of such 110 complex sugars on the shield architecture. Therefore, we have built a glycosylated Env having 111 both oligomannose and complex glycans (Behrens et al., 2016; Cao et al., 2017 ) based on the 112 information provided by mass spectrometry (MS) (Behrens et al., 2016) . To our knowledge, this is 113 the first time that a complete computational model of Env SOSIP was generated to include 114 processed glycoforms, thus obtaining results which are immunologically more relevant. 115
Using these atomistic models, we have employed graph theory to capture the glycan shield 116 topological network, pinpoint potential interaction pathways, and identify concerted behavior of 117 the glycans. The glycan-glycan and glycan-protein interactions influence the behavior of the shield 118 as a whole and can affect distant sites in the glycan network (Behrens et al., 2016; Ferreira et al., 119 2018; Stewart-Jones et al., 2016). Analyses of various network attributes, such as relative centrality 120 of different glycans and critical subnetwork properties, have aided in detailed examination of the 121 native glycan shield in the context of bNAb interactions. Important global and local feature 122 differences due to native-like glycosylation pattern come to light as compared to all high-123 oligomannose glycans, arising from the presence of charged sialic acids at the tips and fucose at 124 the base. Both the simulated models and cryo-EM maps identify a change in projection angle 125 relative to the peptide surface for fucosylated glycans. We have also been successful in quantifying 126 the glycan shield based on the number of sugar heavy atoms encountered over the antigen surface 127 that can be implemented to define regions of relative vulnerability and regions where the glycan 128 shield blocks access. Due to the rapid yield rate of this method, it can be applied to a large number 129 of diverse HIV-1 sequences or can be seamlessly extended to model glycan shields of other viruses 130 in a rapid manner. ) Table S1 . The structure schematic for each of the selected glycan species have 153 been shown in Figure 1A . 154
In the static pictures as obtained from fully glycosylated structures (Gristick et al., 2016; 156 Stewart-Jones et al., 2016) and even single snapshots from computational models, a significant 157 fraction of the protein surface appears to be exposed, with each glycan taking up a particular 158 conformation ( Figure 1B) . However, due to their flexible and dynamic nature, these glycans 159
are not confined to a particular conformation, instead sampling a large volume in space. A 160 previous simulation study suggested (Guvench et al., 2011 ) that the root mean squared 161 deviations (RMSD) for the carbohydrate regions are more than 4 times larger than that of the 162 underlying protein loops. As a result of this conformational variability, the cumulative effect of 163 the glycans over time is like that of a cloud of glycan atoms that shield the underlying surface 164 from any approaching protein probe (an antibody for example) as illustrated in Figure 1C . trimer, like all HIV-1 Envs, has a highly dense glycosylation pattern and each glycan samples a 256 particular region in space (SI Figure S3A-B) , such that neighboring glycans can occupy 257 overlapping regions ( Figure 1D ). The fraction of volume overlap gives a measure of the 258 interaction probability between the constituent glycans (SI Figure S3C) . Within a protomer there 259 are three main regions of overlap -the apex, the gp41 base, and the glycan-dense central high-260 mannose patch (HMP). Inter-protomer overlaps are primarily due to V1 and V2 loop glycans near 261 the trimer apex. Figure 4A shows the obtained network of BG505 native glycosylation unfolded 262 and laid out in two dimensions for the ease of visualization. This is a force-directed layout, which 263 uses attractive forces between adjacent nodes and repulsive forces between distant nodes to reach 264 an optimum distribution of the node points in space (Fruchterman and Reingold, 1991) . For 265 conceptualization, the network relative to the BG505 structure is provided in SI Figure S3D . . We calculated the shortest paths to all glycans from the N332 glycan site using the 288 Floyd-Warshall algorithm (Floyd, 1962; Warshall, 1962) , where the glycan distances were edge-289 weighted ( Figure 4B ). Our calculation shows the most probable pathway over which each glycan 290 feels the influence of the glycan at site N332, providing useful insights into how these inter-glycan 291 communications can occur. 292
To elucidate the relative influence of each of the glycans on the network, we calculated the 293 eigenvector centrality, or eigencentrality, of the nodes to measure its connectivity to the network. 294
A large relative value indicates that the node is well connected to the network and thus is 295 "centrally" located, whereas a low relative value indicates that the node is on the periphery of the 296 network. Operationally, the eigencentrality is the eigenvector associated with the largest 297 eigenvalue of the adjacency matrix, that in position (i,j) reports the overlap of the ensemble of 298 glycan at sites i and j. Value of that particular matrix element will be zero if the glycans in position 299 i and j do not physically interact (details given in Methods section). The normalized eigencentrality 300 of the glycans are projected on the network as a colormap in Figure 4A . 
In order to interact and neutralize the Env protein, antibodies look for breaches in the shield 354 where glycan densities are sparse. This information of relative sparsity is encoded within the 355 glycan network, in terms of edge clustering. We examine this clustering of glycan interactions 356 within the network by defining communities using a modularity maximization approach that 357 divides the network into sub-modules or groups (see Methods). Communities have high density of 358 edges between the nodes within a module and comparatively sparse edges connecting different 359 modules (Figure 5D inset) . This analysis identifies five distinct communities within the BG505 360 native glycan network (Figure 5D and E) . roughly correspond to our modules 1, 2, (3+4) and 5. However, in that study, mannose-5 glycoform 371 was used at all sites, and due to the smaller length of these glycans, some of the interactions, 372 including inter-protomer interactions were not observed. Remarkably, in that study, junctions 373 between microdomains were found to indicate regions of relative vulnerability. The communities 374 we identified, also demarcate the regions where the glycan shield can be penetrated. Broadly 375 neutralizing antibodies whose binding epitopes are known target these community boundaries 376 ( Figure 5D) . Thus, glycan community dynamics can help to determine susceptible regions of the 377 glycan shield, and can be further used for guiding immunological studies. 378
To capture the impact of the glycan shield acting as an immunological barrier, we 380 quantified the sugar barrier over the Env protein surface using the ensemble of structures that have 381 been generated. Towards this goal, we have defined the glycan encounter factor (GEF) as the 382 number of glycan heavy atoms encountered by an external probe approaching the protein surface. 383
This parameter is calculated at each residue on the surface of the protein within a probe of diameter 384 6Å calculated using our ensemble. The highest value obtained is 14 for the native model, and is 385 located in the HMP. Based on existing structures, the main interaction points between Env and 386 bNAbs are often hairpin-like loop regions. Even large-scale atomistic simulations suggest that the 387 first line of contact between Env and an antibody is through a loop region (Schmidt et al., 2013) . 388
Accordingly, for our analysis, we have used a probe size of 6Å, which is the typical diameter of a 389 hairpin loop. At each residue present on the surface of the protein, the approaching probe was 390 considered in three directions ( Figure 6A) : perpendicular to the surface z, and then the x and y 391 directions spanning the plane parallel to the surface. The geometric mean of the x, y and z values 392 provide the GEF per residue on the Env surface. This value will go to zero when the glycan 393 encounter factor is zero from any one of the three cardinal directions. Thus, for any point on the 394 surface which has a dense glycan covering, such as D in Figure 6B , has a high glycan encounter 395 factor value, versus a point such as R where the glycan covering is sparse, which will have a low 396
Previous evidence suggests that 70% of the Env ECD surface area is covered by 398 glycans (Pancera et al., 2014) . Based on this, we determined the lower cut-off of GEF below which 399 we can define glycan holes. Figure 6C shows GEFs mapped onto the trimer protein surface. The 400 calculated solvent-exposed surface area of the protein part of our modeled BG505 structure 401 without considering glycans is 86,055 Å 2 . Excluding the exposed region at the base of the soluble 402 SOSIP, this reduces to 79,672 Å 2 . Varying the lower cut-off of GEF, we find that a cut-off of 1.5 403 GEF leads to 30% of the surface to be exposed. Regions of BG505 surface having a GEF less than 404 1.5 is colored in black in Figure 6D . Deletion of 197 does not only affect the V1/V2 loop region glycans. The difference in adjacency 443 matrix due to this perturbation is illustrated in SI Figure S7A . There is a slight modification of the 444 community structure within the network as well, with glycan 392 leaving module 4 (red) and 445 entering module 3 (yellow) instead (SI Figure S7B) . Interactions such as those between glycans 446 N88 -N625, N133 -N363, and interprotomer N156 -N185e are reduced more than 10%. On the 447 other hand, glycan 276, which was originally interacting strongly with 197 now forms new 448 interactions with N295, N332 and N411 glycans (SI Figure S7C, D) . In summary, HIV and the human immune system are at war constantly, and the virus uses 521 the Env glycan shield to mask the human immune surveillance. One battlefront in this war is the 522 glycan shield: while the virus evolves to develop resilience, the immune response counteracts by 523 looking for vulnerabilities. As efforts are underway to aid the immune system overcome this race 524 by conditioning it with engineered immunogens, there is a need to quantify the resistance and 525 vulnerability of glycan shield in a more quantitative manner. This is that first time that the native 526 glycan network and shielding has been spatially quantified. Our derivation of the Glycan 527 Encounter Factor measures the relative barrier over the Env protein surface, and can potentially 528 aid to distinguish subtle differences on the shield due to variations in the glycosylation or even 529 protein sequences. We are therefore armed with a set of validated in silico tools with which to help 530 the anti-HIV war efforts and guide immunogen design. 531 532
Further information and requests for resources and reagents should be directed to and will be 535 fulfilled by the Lead Contacts S. Gnanakaran (gnana@lanl.gov) and Andrew B. Ward 536 (andrew@scripps.edu). 537
High-throughput conformational modeling 539
The ensemble of BG505 glycoprotein 3D conformations were built in atomistic detail by Ligand-free SOSIP structure with PDB accession ID 4ZMJ (Kwon et al., 2015) was used 547 as the base template for gp120, that also guided the three-fold symmetry of the trimer. 548 5CEZ (Garces et al., 2015) was used as the base template for gp41 and since it has the least 549 number of missing residues among the available structures. we devised the following protocol. We aimed to construct 10 distinct protein models for 554 these flexible regions as starting scaffolds for glycan building. In order to model these, 34 555 different BG505 SOSIP experimental structures were selected below a resolution cut-off 556 of 4Å (see SI Table S2 ). Open-like, chimeric and synthetic models were not considered. 557
The set of 34 were structurally aligned using PyMOL (Schrödinger, 2015) and their residue-558 wise backbone RMSF was calculated (SI Figure S9A) . The base-level of 0.1 was used as 559 a threshold in the normalized RMSF measure. The unstructured residues (coils and turns) 560
were identified on 4ZMJ and 5CEZ using the STRIDE program (Frishman and Argos, 561 1995) . All the identified regions with RMSF higher than the threshold (SI Figure S9A missing residues in many of the PDBs (SI Figure S9B) . Wherever enough models were 567 available, 10 distinct conformations were randomly selected as templates. If this number 568 10 could not be reached for some regions (for example, in loops V2 and V4) due to missing 569 residues (SI Figure S9B) , those residues were modeled ab initio. All known disulphide 570 bonds were kept as additional restraints. We made sure these short-listed regions for ab 571 initio modeling were less than MODELLER-prescribed 14 residues in length. This is 572 because, modeling of such long loops are usually considered to be error-prone (Fiser et al., 573 2000) . The loop region 546-569 in gp41 is present only in the structure 5CEZ. We did not 574 model this region ab initio, and used 5CEZ as the only template for all models. The 10 575 models thus built were further tested for proper stereochemistry, based on MODELLER 576 optimization scores and PROCHECK scores (Laskowski et al., 1993) is the average atomic position of heavy atom k in glycan n. K is the total number 606 of heavy atoms in the glycan. It is 127 for mannose-9, and varies depending on the type of glycan. 607
The ensemble for each model contains 1000 snapshots, making M = 1000 snapshots for each of 608 the two models. The standard deviations were obtained by dividing the 1000 snapshots into 4 sets 609 of M=250, and calculating the four sets of RMSF values (SI Figure S5) . 610
The 1000 structures of the total ensemble are built from 10 initial starting protein 611 conformations, as described above. The main difference between these 10 conformations are the 612 variations in the loop regions, due to the missing residues in the templates. In order to remove the 613 effects of the loop fluctuations and consider the RMSF coming from the glycans alone, the reduced 614 RMSF values were also calculated in each of these 10 sub-models, and their average and s.d. 615 calculated (SI Figure S5B) . The RMSF difference between the models (SI Figure S5E) were 616 obtained by subtracting the reduced RMSF values of all-man9 from native model (native minus 617 all-man9). Since, the average reduced RMSF value is ~4Å, only those RMSF difference values are 618 statistically significant which are above 0.2 Å, which corresponds to a p-value of 0.05, rejecting 619 those values within the null hypothesis. Ensemble of individual glycans were converted to atomic 620 densities at 1 Å gridsize in order to calculate the volume. Volume calculation was done using 621 CHIMERA package (Pettersen et al., 2004) . 622
The inter-glycan overlap is calculated as the total fraction of heavy atoms from the two 624 glycans that come within 5Å of each other. Let us consider the example of mannose-9 to illustrate 625 the parameter of overlap. A single mannose-9 glycan has 127 heavy atoms. Since our ensemble is 626 composed of 1000 possible structures, there are effectively 127,000 heavy atoms per ensemble of 627 mannose-9 at one position. The fraction of the total number of heavy atoms from two neighboring 628 ensembles that come within contact distance defines the overlap fraction. Since mannose-9 is the 629 most commonly occurring glycoform in our system, we have used it as our reference for 630 normalization of the overlap probability. An overlap greater than or equal to 50% of heavy atoms 631 from two neighboring mannose-9 glycans is assigned as 1. This overlap matrix is used to define 632 the adjacency matrix for our network analysis. Each glycan functions as a node of the graph 633 ( Figure 4A inset) , and two nodes are connected by an edge if there is at least 5% overlap as per 634 our overlap definition given above. The edge length is inversely proportional to the overlap value, 635
i.e., the larger the overlap, the closer two nodes (glycans) are in the graph. Only those glycans 636 from the neighboring protomers are considered, that have an inter-protomer edge. All graph theory 637 and network analyses were performed using Python (Rossum, 1995) Where HH is the fraction of edges in module i; ai is the fraction of edges with at least one end in 661 module i. It is calculated by a greedy heuristic, beginning with the trivial system of each node 662 being a cluster, and merging two clusters that will increase modularity by the largest value, 663 stopping when any further merge would decrease the modularity. This approach is known to work 664 well for small networks similar to our system. The calculations were implemented through a 665 standard algorithm (Blondel et al., 2008) present on the surface of the protein, the approaching probe is considered in three directions x, y 1100 and z. (B) Any point on the surface which has a dense glycan covering, such as D, has a high 1101 glycan encounter factor value, versus a point such as R where the glycan covering is sparse, which 1102 will have a low GEF. (C) Representation of GEF on Env surface, given by a colormap. Regions 1103 of BG505 surface having a GEF less than 0.14 is colored in black (D). Typical known antibody 1104 epitopes are indicated by arrows. BG505-specific GH and COT epitope region demarcated by 1105 yellow dashed circle. 1106 Footprint of glycan N197 on Env surface, colored by glycan encounter factor contributed by N197 1110
alone. (C) Glycan encounter factor over the PG9 and (D) PGT145 epitope residues. GEF decreases 1111
for both the epitopes due to deletion of N197 glycan. 1112