Varying molecular interactions explain aspects of crowder-dependent enzyme function of a viral protease

Biochemical processes in cells, including enzyme-catalyzed reactions, occur in crowded conditions with various background macromolecules occupying up to 40% of cytoplasm’s volume. Viral enzymes in the host cell also encounter such crowded conditions as they often function at the endoplasmic reticulum membranes. We focus on an enzyme encoded by the hepatitis C virus, the NS3/4A protease, which is crucial for viral replication. We have previously found experimentally that synthetic crowders, polyethylene glycol (PEG) and branched polysucrose (Ficoll), differently affect the kinetic parameters of peptide hydrolysis catalyzed by NS3/4A. To gain understanding of the reasons for such behavior, we perform atomistic molecular dynamics simulations of NS3/4A in the presence of either PEG or Ficoll crowders and with and without the peptide substrates. We find that both crowder types make nanosecond long contacts with the protease and slow down its diffusion. However, they also affect the enzyme structural dynamics; crowders induce functionally relevant helical structures in the disordered parts of the protease cofactor, NS4A, with the PEG effect being more pronounced. Overall, PEG interactions with NS3/4A are slightly stronger but Ficoll forms more hydrogen bonds with NS3. The crowders also interact with substrates; we find that the substrate diffusion is reduced much more in the presence of PEG than Ficoll. However, contrary to NS3, the substrate interacts more strongly with Ficoll than with PEG crowders, with the substrate diffusion being similar to crowder diffusion. Importantly, crowders also affect the substrate-enzyme interactions. We observe that both PEG and Ficoll enhance the presence of substrates near the active site, especially near catalytic H57 but Ficoll crowders increase substrate binding more than PEG molecules.

a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 [10]. PEG, dextran, and Ficoll crowders also enhanced the enzymatic activity of adenylate kinase (AK3L1) [11]. Dextran was shown to induce the activity of the ribonuclease T1 most probably because the crowder promoted the folding of this protein [12]. We have also found that the kinetic parameters of peptide hydrolysis catalyzed by trypsin [13] and HIV-1 protease [14] were affected by PEG and BSA crowders. Our experiments on another protease, NS3/4A, encoded by the hepatitis C virus (HCV), showed different effects of Ficoll and PEG, as well as BSA on the proteolytic reaction [15]. These observed differences are the main motivation for the work presented here. Since enzymes are important drug targets, it becomes crucial to understand how and why the reaction environment modulates the enzymatic rates and equilibria. In general, different factors contribute to the observed changes in crowding-induced catalytic activity such as increased concentration of the substrate due to decreased space, induced folding of the enzyme to its functional form, induced oligomerization of the enzyme to an inactive form, decreased diffusion of the reactive molecules, occupation of the substrate binding or active site by crowders, and the change of the conformation or flexibility of the binding or dimerization site.
Many of the above effects cannot be elucidated using experimental methods alone. Therefore, the impact of crowding on different solute properties and stages of enzymatic reactions has also been investigated using computational approaches [16][17][18][19]. However, apart from a few atomistic models [4,20], most simulations applied simplified crowder models that primarily mimicked only the effect of volume exclusion due to the impenetrability of molecules, e.g., [21][22][23]. One reason for using such simplified crowder models is computational efficiency, but the choice of such crowders also reflects the assumption that the artificial crowders mainly introduce the volume exclusion effect, and a more detailed representation of molecular-level interactions with crowders may not be needed. However, there are now many studies that describe crowder-type dependent effects on protein structure [24][25][26][27], interaction with substrates [28,29], and enzymatic function [15,30,31].
We focus here on the NS3/4A protein complex which is a bifunctional heterodimeric enzyme composed of the serine protease and RNA helicase domains that can function independently [32]. The NS3/4A activity to process the HCV polyprotein is necessary for viral replication. HCV is responsible for a chronic illness affecting nearly 200 million people and this enzyme is the main target in HCV therapy. However, the emerging drug-resistant genotypes and subtypes of this virus necessitate the design of new inhibitors. The NS3 protease domain is a 180-residue protein, with the chymotrypsin-like fold dominated by β-sheets (Fig 1A). The active site contains a catalytic Ser139-His57-Asp81 triad (with the residues referred to as S139, H57, and D81 in the remainder of the manuscript). The structure is stabilized by a Zinc ion, which is at the opposite side of the active site. Zinc is coordinated by three cysteines (C97, C99 and C145) and one histidine (H149 (Fig 1B) that are conserved in the HCV strains [32] and mutations of any of the cysteines into alanines significantly decreased protease activity [33]. For efficient catalysis, the NS3 protease domain requires a 54 amino acid long peptide co-factor, called NS4A, but only its central part, embedded in the NS3 core, is sufficient for the in vitro activity of the protease. The N-and C-terminal parts of NS4A, which are not bound with NS3, are believed to be intrinsically disordered in solution. In the host cell, the NS3/4A complex associates with the endoplasmic reticulum (ER) membrane where the HCV RNA is replicated. The disordered N-terminal tail of the NS4A cofactor can be induced to form a helix that inserts into the ER membrane to anchor the NS3/4A complex [34,35].
Using fluorescence spectroscopy, we have found that PEG 6000 and Ficoll 400 differently affect (in a concentration-dependent manner) the initial and maximum velocities as well as the turnover number for the peptide hydrolysis catalyzed by the NS3/4A protease [15]. The kinetic parameters decreased in the presence of PEG and increased with Ficoll as compared to dilute buffer solutions. Ficoll also increased by 40% the inhibition constant K i of telaprevir, the clinically-approved drug targeting the NS3/4A enzyme [15]. We have previously explained how PEG crowders may promote folding of the NS3/4A unstructured N-terminal fragment possibly enabling membrane anchoring [35]. However, it remains unclear how the two different crowders may impact enzyme function in an opposite manner. We hypothesize that there are specific interactions between the crowders, NS3/4A, and the substrate that differ between PEG and Ficoll, and the goal of this work is to not just explain the experimental data for NS3/4A but gain a broader understanding of how PEG and Ficoll may interact with proteins at a molecular level.
To probe these questions, we use extensive atomistic molecular dynamics simulations of NS3/4A with and without peptide substrates in the absence or presence of PEG and Ficoll-like crowders. We analyze how the crowders impact NS3/4A structure and dynamics, diffusion of the enzyme and the substrate, and how the crowders may specifically impact the substrate binding near the active site. Based on this analysis we explain how PEG and Ficoll may differently impact the kinetic parameters of the NS3/4A-catalyzed reaction and provide more general insights into how PEG and Ficoll crowders interact with peptides and proteins.

Preparation of solute structures
As in our previous studies [15,35], the NS3/4A crystal structure (PDB ID: 4JMY [36]) was used as the starting point for the MD simulations. This structure contains the coordinates of the full protease domain of NS3, residues 21-32 of the NS4A cofactor, the zinc ion critical for maintaining the active structure of the protease [37] and a fragment of the peptide substrate. To  match the 1b genotype of HCV, the D30E, L36V, G66A, A87K, M94L, S147F, V150A, I170V,  A181S, and S182P mutations were introduced using Chimera [38] with the Dunbrack rotamer library [39]. The missing N-and C-terminal residues of the NS4A cofactor were added with MODELLER [40], ver. 9.19., using a template-free modeling method. Five models of complete NS4A were generated, and, based on the Dope [41] scoring function, the lowest-scoring conformation was selected. Residues were assigned standard protonation states at pH 7, with histidine neutral and protonated at Nδ1 atom. Zinc-coordinating cysteines were modified with a specific cysteine-zinc patch to reflect deprotonation and polarized charges with coordinating cysteines modeled with a -1e charge and the Zn ion with a +2e charge [42,43]. The goal of these parameters was simply to maintain Zn coordination, whereas charge transfer and quantum corrections would need to be considered to achieve a more realistic description of zinc binding equilibria [44,45]. The net charge of the entire NS3/4A-zinc complex was +2e.
The peptide substrates with the sequence Ac-DEDEEAASK-NH 2 , corresponding to the sequence of the peptide used in experimental assays [15], were built with CHIMERA [38]. The substrate models were subjected to 2,000 steps of minimization and 20 copies of 200 ns MD simulations in NAMD [46]. Next, the substrate trajectories were clustered with ProDy [47] and the five most probable conformations were used as starting structures surrounding NS3/4A. The NS3/4A protease and substrates were parameterized using the CHARMM36m force field [48] which is suitable for simulating proteins with unstructured fragments.

Preparation of crowder models
All-atom 28-mer PEG crowder models (Fig 2) were prepared as described elsewhere [15,35]. Starting conformations were taken from MD simulations of a single PEG polymer in explicit water performed by Lee et al. [49]. The CHARMM36 CgenFF force field was used with a modification of the torsion angles based on the work of Leonard et al. [50]. Simulations of PEG using these parameters has been found to result in good agreement with experimental observables in another study [51].
All-atom Ficoll models were built to match the mass of the PEG crowders (cf. S1 Text, S1 and S2 Figs for additional details) so that crowders with similar sizes can be compared. Our Ficoll-like molecule (referred subsequently as simply 'Ficoll' here) consists of four sucroses and three glycerol linkers, forming the smallest possible branched Ficoll-like polysucrose (Fig 2). Sucrose and glycerol parameters were assigned based on the CHARMM36 force field parameters for carbohydrates [52,53]. Additional patches were constructed for sucrose-glycerol connections, with parameters derived by analogy to already parameterized polysaccharides in the CHARMM force field, e.g., using parameters for isomaltulose and melezitose (cf. S1 Text, S3 Fig). To obtain the Ficoll starting conformation for the simulations of the crowded systems, one tetra-sucrose was energy-minimized and an MD simulation in explicit TIP3P water was carried out over 200 ns. The starting Ficoll crowder conformations were randomly chosen from the 20% of the MD snapshots with the lowest R g values.

Simulation systems
Six types of systems were simulated: NS3/4A protease with and without 10 peptide substrates, NS3/4A surrounded by 130 PEG polymers, with and without the substrates, NS3/4A surrounded by 110 Ficoll polymers, with and without the substrates (see S1 Table). All simulations were performed with explicit TIP3P water molecules and 2 Clions to neutralize the Zn 2+ ion coordinated by NS3. Na + ions were added to neutralize the simulations with 10 substrates (each substrate has a net charge of -4e). Additional Na + and Clions were added to achieve an ionic strength of~20 mM. The MMTSB Toolset [54] was used to solvate the systems and add ions. Initial snapshots of the simulation systems are shown in Fig 3. Before solvating the systems, 130 copies of PEG or 110 copies of Ficoll molecules were added in random locations around NS3/4A. Three different sets of crowder positions were generated. The numbers of crowder molecules were chosen to match the mass concentration of crowders; the resulting volume fractions of crowders varied between systems, ranging from 11 to 15%. To avoid clashes, an additional step was carried out for the protease-PEG systems. They were minimized in vacuum before adding water molecules. Considering the high degree of crowding, substrates were added manually using the VMD [55] tools for the systems comprised of both crowders and substrates.

Molecular dynamics simulations
The simulations were performed for the systems listed in S1 Table. Each system was first energy-minimized over 3,000 steps using the steepest descent algorithm. Second, water molecules and ions were thermalized by increasing the temperature from 10 to 310 K in steps of 10 K. Each step was simulated for 25 ps, with 10 kcal/mol/Å 2 harmonic constraints on the positions of the solute atoms. Third, the restraints on the solute were decreased in six 25-ps steps with harmonic force constants gradually decreasing from 10 to 5, 2, 1, 0.1, and 0 kcal/mol/Å 2 . These thermalization and equilibration steps were performed in the NVT ensemble, with a time step of 1 fs. Finally, an additional equilibration step of 25 ps was performed in the NPT ensemble with a time step of 2 fs using the SHAKE algorithm [56] to constrain bonds made by hydrogen atoms. For each system, a 500 ns production simulation was performed either using NAMD 2.11 [46] or OpenMM [57]. A temperature of 310 K and a pressure of 1 atm were maintained with the Langevin thermostat and Langevin piston methods [58]. Two friction coefficients in the Langevin thermostat were used: 1 ps -1 , the value used in our previous simulations [15,35], and 0.01 ps -1 . The lower value gives more accurate estimates of kinetic properties such as diffusion [59,60]. All trajectories were used to evaluate thermodynamic properties but only the trajectories obtained with a value of 0.01 ps -1 were used to analyze diffusive properties and contact life-times. To calculate the electrostatic interactions, the Ewald summation [61] method was used, with a cut-off of 12 Å. To avoid aggregation of the crowders, interactions between the solute and explicit water molecules were increased by a factor of 1.09 because it was found that this modification ensures realistic diffusion of the crowders [35,60].
Multiple replicates of each system were simulated over at least 500 ns. The NS3/4A enzyme was stable in all simulations as indicated by Cα root-mean square deviation (RMSD) values generally between 1 and 2 Å throughout the trajectories (S4 Fig). Replicate-averaged RMSD values increased significantly during the first tens of nanoseconds, but they changed slowly afterwards (S5 Fig). Therefore, we decided to discard the first 50 ns as equilibration from all analyses. To compare the same simulation intervals, all trajectories were analyzed over 50-500 ns, even though some simulations were somewhat longer (S1 Table).

Trajectory analysis
RMSD and RMSF values were calculated for the Cα atoms, after NS3/4A superposition onto the reference crystal structure as a starting point. Translational diffusion was calculated from the mean-square displacements (MSD) of the centers of mass of the molecules. Diffusion coefficients, D 0 , were calculated from the slopes of linear fits to MSD(τ) = h(r(t+ τ) − r(t)) 2 i according to the Einstein equation: where r(t) and r(t+τ) are the vector positions of the center of mass of a molecule in time t and t+τ (τ is the time lag). D 0 values were corrected for the artifacts resulting from using the periodic boundary conditions (PBC) [62] and reduced viscosity of the TIP3P water model (as the ratio of the shear viscosity of TIP3P equal to 0.308 cP and bulk water-0.89 cP) giving: The correction D PBC was calculated as: where T is the temperature set to 310 K, k B is the Boltzmann constant, B is equal to 2.837, R h is the hydrodynamic radius of a molecule (calculated with HYDROPRO [63]), and L is the length of the simulation box. η = η w (1+2.5φ) describes the shear viscosity of the solvent as the viscosity of the TIP3P water scaled with a factor depending on φ, i.e., the volume fraction of crowders. Rotational diffusion of the NS3/4A protease was determined using the method introduced by Case et al. [64]. First, a trajectory of a 1,000 unity vectors pointing in random directions was merged with a trajectory of the centered protease. Second, the vectors were rotated along with the protease rotations by fitting the protease in each frame to a reference structure. To calculate the rotational correlation times, the average correlation function for the vectors was then fitted to a double-exponential function: Obtaining the τ 1 and τ 2 times describing rotations in longer and shorter time-scale, with a being the weight corresponding to the shorter rotational correlation time. Rotational diffusion coefficients were related to calculated correlation times as: where the 3.08/8.9 is the scaling factor included to account for the reduced viscosity of the TIP3P water model (see above) and the overall t was calculated according to: Correlation analysis was used to estimate the time scales of contacts between the molecules and of the dynamics of active site residues. For contact analysis, a contact function was used to describe whether a contact between the molecules present at time t 0 is also present at time t 0 +Δt: where N is the total number of trajectory frames, Δt is the k-th time interval and N p is the number of molecule pairs, and the summation runs over j trajectory frames and i pairs of molecules. The function δ i (t) takes the value of 1 when the distance between the i-th pair of molecules at time t is smaller than 5 Å, and 0 when the molecules are further away from each other. To obtain estimated times of contact survivals, the contact functions were fitted to double exponential functions according to Eq 4, but including an additional constant to allow for contacts persisting beyond the length of the simulations. The above analyses, along with radii of gyration and radial distribution functions, were conducted using the MMTSB Toolset [54], CHARMM version 46b2 [65], and in-house scripts, in particular to automate efficient data processing. The MMTSB Toolset was also used to cluster substrate binding poses. VMD was used to prepare the structural figures and analyze the hydrogen bonds, with the maximum donor-acceptor distance set to 3.5 Å and the donorhydrogen-acceptor angle criteria set to 120 degrees. Secondary structures were analyzed using the Stride algorithm [66]. Plots were generated using scripted Gnuplot, version 5.2 [67].

Results and discussion
The analysis of the simulations is organized by first describing the structural stability, conformational dynamics, and diffusive properties of NS3/4A as a function of crowders. We then describe how different crowders interact with NS3/4A and finally how crowders affect the substrates and substrate-enzyme interactions. Finally, we analyze the interaction of water molecules with NS3/4A and the crowders. All of the results are then interpreted in the context of the experimentally observed functional differences of the enzyme in the presence of different crowders.

Structural stability of NS3/4A in the presence of crowders
The structure of the NS3/4A enzyme remained stable throughout the simulations with average Cα RMSD values of about 1.6 Å for NS3 and the 13 residues of NS4A resolved in the crystal structure (S4 Fig and S2 Table). RMSD values appeared to be slightly larger in the presence of Ficoll or PEG crowders compared to water; in the presence of substrates, RMSD values were smaller with Ficoll and larger with PEG. However, in all cases, the differences in RMSD values were on the order of the standard errors and therefore not considered significant.
Average radius of gyration values R g of 16.2 Å for NS3 (S3 Table) were slightly larger than the value of 15.8 Å calculated from the experimental coordinates based on 4JMY. Again, there were only slight differences in the absence or presence of crowders, but it appears that the distribution of R g values is shifted to slightly larger values in the presence of NS4A has flexible N-and C-termini which result in larger R g values and larger fluctuations (S6C Fig). Average R g values of NS4A appear to be smaller with Ficoll and larger with PEG compared to the system without crowders (S3 Table); with substrates the opposite trend is found, i.e., NS4A is less compact in the presence of Ficoll and more compact with PEG ( S6D  Fig). However, all of the differences are well within the standard errors and therefore not considered significant.
The presence of crowders is generally believed to favor more compact conformations, especially for more dynamic structural elements. However, for the system studied here, the effect of the crowders on the NS3/4A structure is very modest if any, including for the extended NS4A termini.

Conformational dynamics of NS3/4A in the presence of crowders
Since NS3 remained overall stable, no major conformational changes were observed. The remaining conformational dynamics characterized via root-mean square fluctuations (RMSF) around the trajectory-averaged structures ( This suggests that both crowders perturb the conformational sampling of NS3, Ficoll more than PEG, but with the detailed effects of the two types of crowders being qualitatively different. Previous work has suggested that altered conformational fluctuations in the presence of crowders may correspond to different overall stabilities of the folded state in the presence of the crowders [68]. Some of these effects are retained when substrates are introduced into the systems (S7B Fig), but together with substrates, PEG has a greater tendency to lead to increased fluctuations, for example in the entire N-terminal helix (residues 12-23), and near residue 89, whereas Ficoll with substrates results in more similar fluctuations than without Ficoll. The NS3 N-terminal helix (Fig 1A, magenta) is essential for membrane association of the NS3/4A complex [34] and both PEG and Ficoll seem to increase its conformational dynamics.
The cofactor NS4A has flexible N-and C-termini giving rise to greater RMSF values (S9 Fig). Crowders somewhat reduce the fluctuations, but without substrates the effect is only statistically significant for the N-terminus. When substrates are present, the fluctuations in the NS4A termini are reduced more in the presence of crowders, with no clear differences between PEG and Ficoll (S9B Fig).
We previously described the ability of PEG crowders to induce helical structures in the Nterminus of PEG [35]. Such helix formation prepares for membrane anchoring of NS3 and thus we have previously concluded that cellular crowding may assist with this process [35]. Based on the simulations presented here, we again find a greater tendency towards helix formation in the N-terminus of NS4A in the presence of PEG crowders (S10 Fig, S4 Table).
Without substrates, Ficoll does not appear to induce helical conformations to the same extent. However, when substrates are present, both PEG and Ficoll induce significantly more helical structures in NS4A compared to dilute solutions, and PEG induces more helical structures than Ficoll (S11 Fig, S4 Table). Similar conclusions are also found when analyzing the φ/ψ peptide backbone torsion angle distributions in the N-and C-termini of NS4A (S12 and S13 Figs). This confirms an overall role of crowding in promoting functionally relevant secondary structure formation in NS4A, but also indicates that this observation may depend on the choice of the crowder molecule used to mimic cellular environments.

Self-diffusion of NS3/4A in the presence of crowders
Translational and rotational diffusion was analyzed from the trajectories with the reduced friction coefficients as described in the Methods section. Translational diffusion coefficients were obtained from linear fits to the mean-square displacement curves shown in S14 Fig followed by corrections for periodic boundary artifacts and reduced viscosity of the TIP3P water model. The results are given in Table 1. As expected, the presence of crowders significantly reduces the diffusion rates of NS3/4A by about one third from diffusion in dilute aqueous solvent. Interestingly, diffusion is reduced more in the presence of PEG than in the presence of Ficoll. When substrate molecules are present, diffusion is reduced further, both without and with crowders and, again, diffusion is slower with PEG crowders than with Ficoll crowders when substrate is present. The substrate molecules therefore act themselves as additional crowders with substrate and crowders apparently affecting translational diffusion of NS3/4A in an independent, additive manner.
Rotational diffusion coefficients were obtained from double-exponential fits to rotational correlation functions shown in S15 Fig. The resulting rotational time scales are given in S5 Table. Overall rotational diffusion rates, corrected for the faster viscosity of TIP3P water, are reported in Table 2. Different from translational diffusion, the presence of crowders alone appears to affect rotational diffusion rates of NS3/4A less. Diffusion actually appears to be slightly higher with PEG, although this result is based on only two simulations (S1 Table). Moreover, this is a result of slightly faster rotation on fast time scales (3.4 ns with PEG vs. 3.7 ns with water) whereas rotation on slower time scales is retarded with PEG (28.2 ns with PEG vs. 22.7 ns with water) according to the double-exponential fits (S5 Table). Double-exponential fits resulted in significantly better fits to the data than single-exponential fits as indicated by lower χ 2 values (S5 Table). This suggests that rotational diffusion involves multiple time With Ficoll crowders, rotational diffusion of NS3/4A is slowed down more. The presence of substrates reduces rotational diffusion slightly without crowders, but there is a more significant reduction when either PEG or Ficoll crowders are present along with substrates.
The suggestion that translational and rotational diffusion of NS3/4A may be affected differently in the presence of crowders, led us to further look for evidence of anomalous diffusion, i. e., a variation of diffusion rates over different time scales. This could result, for example from a cage effect where crowders do not hinder short-term motion and rotational motions in place but present molecular obstacles when diffusing over longer distances. Based on plots of log (MSD/t) vs. time (S16 Fig) we found that the presence of the crowders does indeed give rise to anomalous diffusion since log(MSD/t) decreases significantly from a peak at around 50 ps towards longer time scales. This effect appears to be stronger with PEG than with Ficoll, both without and with substrate (S16 Fig).
Diffusion of the crowders themselves is also significantly different. Translational diffusion of PEG is significantly slower than that of Ficoll ( Table 1). This is easily explained by the larger average size of PEG (S3 Table) compared to Ficoll, but PEG again appears to rotate faster than Ficoll ( Table 2), potentially indicating differences in crowder-crowder interactions. Finally, substrate diffusion is reduced upon crowding and follows crowder diffusion rates, both for translational and rotational diffusion. For example, substrate diffusion is reduced significantly more in the presence of PEG than in the presence of Ficoll. This suggests that substrates may spend a significant time near crowders and effectively move together. This will be discussed in more detail below.

Interaction of NS3 with crowders
The PEG and Ficoll crowders are present at high concentrations. Based on radial distribution functions (Fig 4), the crowders overall do not have a strong preference for interactions with NS3 since crowder atom densities are higher away from the surface of NS3 than near it. As described previously [35], a peak at 4 Å for NS3-PEG interactions indicates specific interactions when PEG is near NS3. This is largely due to interactions involving PEG's oxygen atoms (Fig 4B). Ficoll mostly lacks distinct features in the RDF, although the number of oxygen atoms (in hydroxyls) that could be involved in interactions with the protein is larger in Ficoll than in PEG (Fig 4B).
To quantify NS3-crowder interactions, we calculated crowder contacts per amino acid residue as well as the average number of crowders in contact with NS3 ( Table 3). Perhaps Results from double-exponential fits to combined correlation functions from replicate trajectories with reduced friction coefficient (see S1 Table), corrected for the TIP3P water model. Standard errors of the mean given in parentheses were estimated from variations between fits to individual correlation functions. surprisingly, there were more than five PEGs and more than three Ficoll molecules in contact with NS3 at any given time. However, this is a consequence of the dense molecular environment with more than 100 crowder molecules (S1 Table). Interestingly, the number of PEGs in contact slightly decreases and the number of Ficolls in contact increases when the substrate is added. Crowder contacts per residue, about 0.2-0.3, tell a similar story ( Table 3). In addition to simply counting crowder contacts, we also specifically analyzed hydrogen bonding interactions as those are likely to be stronger and potentially longer-lasting contacts. The number of NS3-crowder hydrogen bonds are a small fraction of the total number of contacts, but interestingly, Ficoll formed significantly more hydrogen bonds with NS3 than PEG, even though PEG overall interacted more strongly with NS3 based on the radial distribution functions and total number of contacts. This may be explained by the fact that PEG only has ether oxygens acting as hydrogen bond acceptors, whereas Ficoll has more oxygens and its hydroxyl groups can act both as hydrogen bond acceptors and donors.
To further compare NS3-crowder interactions between PEG and Ficoll, we broke down interactions according to NS3 amino acid types (Fig 5). Both crowders interacted with all types of amino acids, with about half of the contacts with charged or polar residues. However, there are significant differences between PEG and Ficoll with respect to interactions with acidic residues (Ficoll interacts more strongly) and the classic hydrophobic residues (PEG interacts more strongly). When projecting crowder interactions onto the NS3 surface (Fig 6), we find broad coverage across most of the surface, but with preferences for certain parts of the structure. As may be expected, both crowders interacted more with the parts of the NS3 structure that are protruding most from the overall global shape, such as some of the solventexposed β-turns. In more detail, there are differences between the parts of NS3 in contact with PEG or Ficoll, presumably reflecting the differences in amino acid preferences. Many of the regions in frequent contact with the crowders correspond to the parts of the NS3 structure, where the presence of the crowders alters conformational fluctuations. This leads to the conclusion that crowder contacts directly modulate NS3 conformational dynamics where frequent contacts are formed. , and for substrate heavy atoms (C). Dashed lines in A and B indicate crowder distributions in the presence of substrates. Distances from the crowder or substrate atoms to the closest NS3 atom were counted and normalized by the total available volume at a given distance from the NS3 surface. For the simulations with only substrates (C), the distribution functions were scaled by a factor of 2 to account for the difference in substrate molarity between the systems without and with crowders (see S1 Table). The dotted line shows the results obtained without scaling. Distribution functions were averaged over all trajectories for the systems with the crowders but without substrates. The shaded areas indicate standard errors. https://doi.org/10.1371/journal.pcbi.1011054.g004

PLOS COMPUTATIONAL BIOLOGY
Finally, NS3-crowder contact life times were analyzed based on double-exponential fits to contact survival correlation functions (see Methods). Again, double-exponentials fit the data significantly better than single-exponential fits (S6 Table). We find that contacts persisted on two life times with roughly equal weight (S6 Table). Short-lived contacts lasted for about 0.2-0.3 ns, with no statistically significant differences between PEG and Ficoll. Longer contacts lasted for about 6.5 ns with PEG and, significantly shorter, for 4.5 ns, with Ficoll. Contacts on different time scales likely reflect different types of contacts with different kinetic barriers to dissociation once contacts are formed. The observation that NS3 forms frequent contacts with multiple crowders and lasting on time scales relevant for diffusion corresponds to the overall slow-down in diffusion of NS3 in the presence of the crowders described above.
The overall picture from the analysis so far is that crowders do not favor interactions with the enzyme, but interactions lasting as long as several nanoseconds occur nevertheless simply due to the highly concentrated environment. Several crowders are in contact with NS3 at any given time and PEG interactions are slightly stronger, with more specific interactions and longer contact life times compared to Ficoll. However, Ficoll forms more hydrogen bonding interactions with NS3 than PEG and has a stronger preference for acidic residues than PEG. The different amino acid preferences lead to different surface interaction patterns explaining different effects of PEG and Ficoll on the internal dynamics of NS3. Further consequences on substrate binding and enzyme function will be discussed in the following.

Interaction of substrates with NS3 and crowders
Peptide substrates were present in part of the simulations, without and with crowders. This allows us to analyze substrate dynamics and NS3-substrate interactions in the absence and presence of crowders as well as substrate-crowder interactions. In contrast to the crowders, the substrate peptides interact preferentially with NS3/4A (Fig 4C). We note that the

PLOS COMPUTATIONAL BIOLOGY
Crowding effects on enzyme function concentration of substrates in terms of molarity is about two-fold higher in the control simulations without crowders than in the simulations with crowders (S1 Table) in order to achieve a comparable number of substrates interacting with NS3/4A with and without crowders. To compare relative preferences of NS3-substrate interactions relative to bulk solvent with and without crowders, the distribution function without crowders was therefore scaled by a factor of two in Fig 4C. We find that the presence of crowders increases substrate interactions with NS3/4A and between crowders; Ficoll crowders increase substrate binding more than PEG crowders. On average, there are about two substrates bound to NS3 at any given time ( Table 3) and substrate contacts per residue are almost the same as crowder contacts per residue. Substrate interactions projected onto the surface of NS3 (Fig 6) show strong preferences for certain parts of NS3, mostly near the active site. The substrate is highly acidic and substrate interactions appear to be guided by overall electrostatic attraction to surface patches of NS3/4A with a high concentration of basic amino acids (Fig 6).
The substrates themselves sample mostly extended conformations with small percentages of helical structure (S4 Table). Crowders do not significantly change the secondary structure sampling of the substrates (S4 Table and S17 Fig). However, in the presence of the crowders, the substrates appear to be slightly more compact as average radii of gyration are reduced by 1-2% (S3 Table).
Substrate contacts with crowders are characterized in Fig 7. As for NS3, the peptide substrates do not interact preferentially with the crowders. In contrast to NS3, the substrate interacts more strongly with Ficoll than with PEG. This is a result of Ficoll preferring interactions with acidic amino acids (Fig 7B) relative to PEG. As a result, there is, on average, almost 1 contact with Ficoll per substrate residue, but only 0.3 contacts per residue with PEG (Table 3). Nevertheless, the number of crowder molecules in contact anywhere with a given substrate is similar, around one crowder per substrate on average ( Table 3). Substrate-crowder life times are somewhat shorter than NS3-crowder life times, with short contact times of 0.1 ns and long contact times of around 1.5 ns (S7 Table). This suggests a picture where substrates mostly travel attached to a crowder. Substrates detach from crowders and exchange interactions frequently, but the contacts still remain long enough for substrate diffusion to be dominated by crowder diffusion as discussed above.

Crowding effects on NS3 active site
We now focus on the specific effects near the NS3 active site with the ultimate goal of explaining the experimental data that showed enhancement of NS3/4A activity in the presence of Ficoll crowders but reduced activity with PEG crowders [15]. According to radial distribution functions, the crowders themselves do not interact strongly with the active site residues, although Ficoll interacts more than PEG (Fig 8). At 10 Å distance, crowder atom densities range from less than 0.0005 Å -3 near S139 to 0.0012 Å -3 near H57 which is much less than the average crowder density of 0.0035 Å -3 with respect to any residue of NS3 (Fig 4). This may be expected since the active site residues lie at the bottom of a surface pocket (Fig 1A-1C) where non-specific binding with the crowders is less likely.
Substrate molecules are expected to bind near the active site for cleavage, but, without crowders, substrate atom densities near the active site are similar to crowder atom densities (Fig 8). However, both crowders significantly enhance the presence of substrates near the active site, especially near H57. At the same time, Ficoll crowder densities near the active site are reduced when substrate is present.
The geometry of the active site residues is important for enzyme function. D81 is expected to coordinate with H57, and H57 would attach the hydroxyl of S139 in the presence of substrate in the first step of the peptide cleavage reaction [69]. Therefore, we analyzed whether the presence of a crowder and/or a substrate at close distance to the active site would affect the conformational sampling and interactions between the active site residues.
We used RMSD of the active site residues (57, 81, 139) to compare the overall geometry of the active site relative to the crystal structure during the simulations. Probabilities of finding different RMSD values are shown in Fig 9. As the structure of the enzyme fluctuates, there is a distribution of conformations ranging from structures very similar to the crystal structure (<1 Å) that are likely reflecting a conformation close to the catalytically competent state to structures that deviate more significantly (up to about 3 Å). The most likely structures seen in the simulations deviate by about 2.2 Å. When crowders are present, close contacts of PEG increase the sampling of crystal-like structures, whereas close Ficoll contacts actually make the sampling of such conformations less likely (Fig 9A). When substrates are present as well, neither crowder enhances the sampling of crystal-like conformations when coming in contact with the active site (Fig 9A). The functionally more relevant scenario involves cases where substrate molecules approach the active site. Without crowders, crystal-like structures are not sampled as frequently as in pure water when substrate molecules come close to the active site (Fig 9B).
The additional presence of PEG does not change that. However, the addition of Ficoll restores the more frequent sampling of crystal-like active site conformations when substrate approaches the active site seen in pure water (Fig 9B). This suggests that PEG and Ficoll crowders act differently in promoting catalytically relevant active site geometries along with substrate binding.
To analyze the effect of crowders and substrate on the NS3 active site in more detail, we calculated the distribution of side chain χ 1 torsion angles and inter-residue distances between active site residues. χ 1 torsion angles as a function of the closest crowder or substrate distance from the active site are shown in S18-S20 Figs. In addition, we compare the distribution of χ 1 angles between dilute solvent without crowders or substrates with the distributions when a PEG, Ficoll, or substrate molecule is within 5 Å of any of the active site residues (S21 Fig). The sampling of side chain orientations is broad, covering different rotameric states and not just the conformation found in the crystal structure (PDB ID: 4JMY) (S21 Fig). This may be expected since the residues are solvent-exposed and the crystallographic conformation may be further constrained by the presence of a tightly bound inhibitory peptide. When crowders interact without substrates, the side chain orientations appear to change somewhat, for example the rotameric state with χ 1 = 60˚for H57 seen also in the crystal structure became less populated when Ficoll crowders interacted but more populated when PEG crowders were close. However, most of these effects disappear when substrates are present along with the crowders (S21 Fig). On the other hand, close substrate contacts themselves also appear to modulate the torsional sampling, especially when a substrate is in very close direct contact (<3 Å). More specifically, χ 1 = 60˚for H57 and χ 1 = -150˚for D81 become preferred states similar to the crystal structures when substrate is in very close contact (<3Å) in the presence of Ficoll (S20 Fig, lower panels). For S139, close substrate contacts promote non-crystal values of χ 1 = -150m ore in the presence of Ficoll than in the presence of PEG, although the crystal-like configuration with χ 1 = -60˚remains the preferred state with both crowders (S20 and S21 Figs).
We further analyzed side chain center distances between active site residues. As the side chain orientations fluctuated, side chain distances also displayed significant dynamics between close contacts of 4 Å to larger distances of 11 Å with broad sampling around preferred distances depending on residue pair (S22-S24 Figs). The results are summarized for close crowder contacts without substrates in S25 Fig. We find that with PEG crowders near the active site, close distances between all active site residues were most likely 4.5 Å for H57-D81, 5 Å for H57-S139, and 8 Å for D81-S139 (S25 Fig). In contrast, Ficoll near the active site led to the opposite effect of active site residues moving further away from each other as the most likely distances increased to 8 Å for H57-D81, 7.5 Å for H57-S139, and 9.5 Å for D81-S139 (S25 Fig). Without crowders the distance distributions were found to be intermediate between the distributions obtained with PEG and Ficoll (S25 Fig). However, as for the torsion angle sample, these effects mostly disappeared when substrates were present along with the crowders since distance distributions with or without crowders became largely indistinguishable when considering the statistical uncertainties (S25 Fig). Close substrate contacts by themselves were also correlated with shorter distances becoming more probable when PEG was present ( S25  Fig). With Ficoll, the same effect was observed, but required very close substrate distances to the active site (<3 Å, S24 Fig). These findings suggest that crowder-induced close substrate interactions may facilitate an induced-fit mechanism of substrate binding in the active site, especially with Ficoll crowders.
We proceeded to look at two key atomistic distances, H57:Nδ-H . . . Oδ1/2-D81 and H57-Nε . . . HOγ-S139, in more detail (S26-S28 Figs). Again, there are significant fluctuations of these distances ranging from hydrogen-bond distances of 2 Å to values beyond 10 Å as a reflection of significant dynamics (S26-S28 Figs). These distances need to be close for the cleavage reaction to proceed. Comparing distance distributions again for the case where crowders or substrates are in close contact with the active site (S29 Fig), we find similar conclusions as for the side chain distances, namely PEG by itself promotes shorter distances whereas Ficoll by itself increases these catalytically important distances. Substrate contacts by themselves without crowders do not strongly promote closer distances since hydrogen-bond distances of less than <3Å are observed much less frequently than longer distances (S29 Fig, bottom panels). However, close H57-Nε . . . HOγ-S139 hydrogen bond distances, necessary to start the cleavage reaction, become relatively more favorable than intermediate distances of 4-6 Å when Ficoll crowders are present (S28 and S29 Figs, right lower panel). In fact, with PEG crowders and substrate, sampling of close distances (<2.5 Å) of this hydrogen bond almost vanishes (S28 Fig), whereas Ficoll enhances this interaction when a substrate is close to the active site.
To complement the conformational sampling analysis of active site geometries, we also characterized side chain fluctuations from orientational correlation function decays (S30 Fig). Double-exponential fits to the correlation functions (S8 Table) reveal fast motions on time scales of around 1 ns and slow motions on time scales between about 50 and 500 ns. The slower time scales presumably correspond to transitions between different rotamer states. We find that the presence of Ficoll crowders appears to accelerate the slow motions over dilute solvent dynamics by 20-50% when substrates are not present (S8 Table), whereas PEG crowders may have slowed down those motions for H57 (by 10%) and S139 (by 11%) while increasing the slow motion for D81 (by 40%) (S8 Table). However, these variations may lie within the uncertainties of our analysis given the estimated errors of the correlation functions (S30 Fig). The presence of substrates by themselves, accelerated the side chain motions significantly, about twofold (S8 Table). When both crowders and substrates are present, faster motions were retained for H57, but mostly slowed down again for D81 and S139 with varying effects depending on the crowder type and residue (S8 Table) that may be difficult to distinguish with certainty given our data. This analysis may suggest that the crowders can influence the dynamics of functionally relevant residues to different degrees with possible consequences for enzyme catalysis. However, much longer simulations are likely needed to come to firmer conclusions.

Crowding effects on substrate binding modes
Finally, we looked at the binding modes of substrates near the active site. Experimental structures of NS3/4A provide partial evidence of where the substrate is expected to bind NS3, mainly based on the binding mode of inhibitory peptides [70,71]. According to these studies, the substrate binds in a mostly extended conformation with the N-terminal part of the substrate interacting with D168 of NS3 and pointing towards the NS3/4A active site (comprised of H57, D81, and S139). We analyzed the preferred binding modes from the simulations based on clustering analysis of substrates found near the active site (Fig 10 and S9-S10 Tables). Substrates interacted in a variety of ways. Only some conformations were close to the exact conformation of the inhibitory peptide bound in the crystal structure (PDB ID: 4JMY). Without crowders, we find that bound substrates were generally located in the surface depressions above and to the left of the active site in Fig 10B, interacting with key residues known from previous studies to be involved in substrate binding (such as Asp168 at the N-terminal side of the substrate and Gln41, Ser42, Arg109, and Lys136 for the C-terminal part of the substrate near the cleavage site [71,72]). However, some of the bound conformations appeared to be relatively far from the surface (Fig 10F). In the presence of Ficoll crowders, substrate binding appeared to become more focused towards the inhibitory peptide binding site in the crystal structure, and mostly involved linearly extended substrates across the active site (Fig 10C-10G). In contrast, the presence of PEG appeared to have an opposite effect of more diffuse binding of the substrate pointing in different directions but without clear examples of the substrate following the binding mode along the cleft above the active site suggested by the inhibitory peptide bound in the crystal structure. Additional quantitative analysis of individual clusters focused on end-to-end distances of substrates to measure their compactness, alignment of substrates with respect to the substrate fragment in the 4JMY structure to assess how similar the substrate orientation was to the presumed binding configuration, and the question to what degree different substrates contributed to a given cluster and how long substrates remained in contact near the active site (S9-S11 Tables). Accordingly, we found that substrates interacting near the active site are mostly compact in the presence of PEG crowders (S10) whereas substrates are significantly mostly extended in the most populated clusters in the presence of Ficoll crowders (S11 Table). At the same time, substrates in the most populated clusters in the presence of PEG are oriented mostly orthogonal to the substrate fragment present in 4JMY, as indicated by near-zero cosθ values (S10 Table). In contrast, substrates in the presence of Ficoll were highly aligned with the substrate fragment (S11 Table), although in some cases the direction of the peptide was opposite as indicated by negative cosθ values. Many of the most populated clusters involve predominantly a single substrate remaining in contact near the active site for longer times (up to 10 ns). However, different clusters capture different substrate interactions with similar features suggesting exchange of substrates and that the findings are not dominated by a single substrate remaining bound near the active site for a very long time.
Along with the differences in substrate binding modes, Fig 10 also shows the preferences of crowder molecules to interact with NS3 residues. It can be seen that Ficoll interacts more prominently on the right side of the active site (near D81 and above), whereas PEG crowders interact much less in that area, presumably because of Ficoll's preference for acidic residues. Based on this analysis, we speculate that Ficoll binding to NS3 guides substrates to more easily find the catalytically relevant binding mode whereas PEG binding to NS3 may actually have the opposite effect of discouraging catalytically competent substrate binding.

Conclusions
Understanding enzyme function in the presence of crowders is important for understanding enzyme function in vivo. Proteins are the main type of crowder molecules encountered in biology but experiments often use artificial crowders to focus on generic crowder properties and avoid the complexity of biological environments. The present study was motivated by the observation that the enzyme kinetics of NS3/4A responded very differently to two different types of artificial crowders, PEG and Ficoll. Although K M values in the absence and presence of crowders were found to be similar, catalytic rates decreased with PEG and significantly increased with Ficoll compared to the enzyme in buffer solution [15]. This suggested that these commonly used crowders may affect enzyme function via crowder-specific interactions rather than the more generic volume exclusion effect typically assumed to be dominant for these types of crowders [1].
The extensive simulations of NS3/4A in the presence of PEG and Ficoll as well as peptide substrates that are presented here offer molecular-level insights into how PEG and Ficoll may interact differently with NS3 and the substrates, and how they may lead to different enzyme behavior. As expected, neither PEG and Ficoll interact strongly with NS3. However, interactions do occur with enough frequency for NS3 to be in contact with several crowder molecules at any given time. Moreover, NS3-crowder interactions are relatively long-lived, up to several nanoseconds. An interesting finding is that PEG and Ficoll have different preferences for different amino acid types. Ficoll clearly prefers interactions with acidic residues, whereas PEG is more likely to interact with amino acids that have hydrophobic aliphatic side chains. As a result, crowders decorate different parts of the NS3 surface according to their amino acid preferences.
The crowders also interact with the substrate peptides, especially Ficoll because the substrates are highly acidic. In fact, based on the analysis presented here, substrates rarely diffuse on their own but are in contact with a crowder at any given time. This suggests that substrate diffusion is essentially slaved to crowder diffusion. The crowders considered here were relatively small, and therefore the effective substrate diffusion was reduced more moderately in the presence of crowders, but one expects this effect to be significantly more pronounced when much larger crowders such as PEG 6000 or Ficoll 400 are used as common in experiments. In that case, reduced substrate diffusion could become a rate-limiting step for enzyme reactions. On the other hand, we found that the presence of crowders enhanced substrate binding to NS3, perhaps by delivering substrates to NS3 and then hindering substrates from leaving once bound to NS3. Moreover, crowders appeared to modulate the substrate binding modes, essentially by binding to certain NS3 surface regions thereby focusing substrate binding on the remaining patches that are unoccupied by crowders. With PEG crowders this led to more diffuse binding, while Ficoll may focus substrate binding towards a catalytically competent pose. In addition, Ficoll crowders also appeared to facilitate the formation of catalysis-relevant conformations of active site residue when a substrate molecule was closely interacting. Taken together, this may explain the experimental observation of increased NS3/4A activity in the presence of Ficoll compared to buffer and when PEG is present.
While the study highlights specific effects of PEG and Ficoll due to different interactions with proteins and peptides at the molecular level, it does not fully interpret the experimental data. Substrate binding near the active site was observed, but it is unclear whether any of the bound substrates fully reached a catalytically competent binding pose. Probably much more extensive simulations with enhanced sampling techniques could overcome this challenge. However, it is also not clear from experiment how exactly cleavable substrates would bind to the active enzyme. Another limitation is that the crowder models used here may compare to PEG 600 or 2000, but not the much larger PEG 6000 and Ficoll 400 crowders used in the experiments. Simulations with much larger crowders would be computationally prohibitive. Another complication is that the higher-order structure of Ficoll is not well-characterized. Finally, the simulations presented here only address classical aspects of biomolecular interactions and dynamics but do not directly describe the peptide cleavage reaction, which requires quantum-mechanical treatments [69]. Studies of enzyme catalysis via QM/MM calculations where crowder molecules are considered explicitly as part of the studied system would be interesting future work.
While the work here focuses on the specifics of the NS3/4A protease and compares the effects of two types of artificial crowders, the more general conclusions are that even artificial crowders may significantly impact functional studies under crowding conditions as a result of molecule-specific interactions with proteins. We hope that the insights gained here will motivate further studies of enzymes under crowded conditions and guide the interpretation of results obtained with different crowding agents. Fluctuations were calculated based on the Cα-Nε vector (H57, A/B), the Cα-Cγ vector (D81, C/D), and the Cα-Oγ vector (S139, E/F) after superposition of the overall NS3 structure to a common reference in water (red) or in the presence of Ficoll (blue) or PEG (green) crowder without (A/C/E) and with (B/D/F) substrates. The correlation function was calculated as the second-order Legendre polynomial from the inner dot product of the orientational vector at different time points. Solid lines show simulation averages, standard errors are indicated by the shaded area. Results from double exponential fits to the correlation functions are given in S8 Table. Only trajectories with the reduced friction were used in this analysis (S1 Table). (TIF) S1  Table. NS3-crowder contact life-times from single-and double-exponential fits to contact survival decays. (PDF) S7 Table. Substrate-crowder contact life-times from double-exponential fits to contact survival decays. (PDF) S8 Table. Double-exponential fits to the combined orientational correlation functions of the H57, D81, and S139 active site side chains. (PDF) S9 Table. Cluster analysis for substrate binding near the active site in simulations without crowders. (PDF) S10