Investigating Intrinsically Disordered Proteins With Brownian Dynamics.

Intrinsically disordered proteins (IDPs) have recently become systems of great interest due to their involvement in modulating many biological processes and their aggregation being implicated in many diseases. Since IDPs do not have a stable, folded structure, however, they cannot be easily studied with experimental techniques. Hence, conducting a computational study of these systems can be helpful and be complementary with experimental work to elucidate their mechanisms. Thus, we have implemented the coarse-grained force field for proteins (COFFDROP) in Browndye 2.0 to study IDPs using Brownian dynamics (BD) simulations, which are often used to study large-scale motions with longer time scales and diffusion-limited molecular associations. Specifically, we have checked our COFFDROP implementation with eight naturally occurring IDPs and have investigated five (Glu-Lys)25 IDP sequence variants. From measuring the hydrodynamic radii of eight naturally occurring IDPs, we found the ideal scaling factor of 0.786 for non-bonded interactions. We have also measured the entanglement indices (average C α distances to the other chain) between two (Glu-Lys)25 IDP sequence variants, a property related to molecular association. We found that entanglement indices decrease for all possible pairs at excess salt concentration, which is consistent with long-range interactions of these IDP sequence variants getting weaker at increasing salt concentration.

1 Introduction

One of the main determinants of biological structure and function is the interaction of two or more molecules, especially protein molecules. Understanding the dynamics of these bimolecular interactions is important for the understanding of such cellular structures as the cytoskeleton (actin and tubulin, for example), ribosomes, chromosomes, and polymerases, as well as processes such as cell signaling and cell motility (Alberts et al., 2002; Pollard and Earnshaw, 2007). Furthermore, the encounter stages of such reactions, which are often the rate-limiting steps, are diffusion-limited (Elcock, 2004). Therefore, the use of Brownian dynamics (BD) is appropriate for such systems [see Huber and McCammon (2019) for a review]. For several decades, BD has found use in polymer and peptide simulations, simulations of enzyme-substrate reactions, and protein–protein association reactions. More recently BD has found use in studies of large-scale cytoplasm simulations, microtubule dynamics, assembly of protein complexes, retroviral infectivity, molecular motors, chromosome organization, the nuclear pore complex, synapses, and endocytosis. The previous version of the Browndye software package (Browndye 1.0), which was limited to two rigid bodies, has been used in enzyme kinetics and channeling (Huang et al., 2018), as well as protein-protein interactions (Grant et al., 2011).

The Browndye 2.0 software package, successor to the previous simulation package, consists of two simulation programs and about 38 auxiliary programs for processing data. Like the previous version, Browndye 2.0 can compute the second-order rate constants of the encounter of two bodies moving according to BD, compute the probabilities of the two bodies moving from one binding mode to another, and output the molecules’ trajectories. The main addition is the ability to model each molecule as a collection of large rigid cores with flexible connectors and loops. In its original two-rigid-body model, Browndye has functionality very similar to the packages SDA (Martinez et al., 2015), MacroDox (Northrup et al., 1993), and GeomBD (Roberts and Chang, 2016), and is intended primarily for simulations of large biological molecules like those three other packages. Its current limitations arise mainly from the structural rigidity approximations and the nature of the force computations between the molecules.

Using Browndye 2.0, we have investigated intrinsically disordered proteins (IDPs), which are proteins that do not have a stable, folded structure and instead take on various structures depending on their current tasks in modulating biological processes. Conducting a computational study of these systems will be critical to elucidate their mechanisms. Specifically, we have implemented the coarse-grained force field for proteins (COFFDROP) (Andrews and Elcock, 2014; Frembgen-Kesner et al., 2015) in Browndye 2.0 to study eight naturally occurring IDPs and five (Glu-Lys)25 IDP sequence variants. We have measured their structural properties, including radius of gyration (R ), interresidue distances (R ), and hydrodynamic radius (R ), and a property related to molecular association, namely the entanglement index (average C distance to the other chain).

2 Materials and Methods

2.1 Structure Preparation

The Alphafold Colab (Jumper et al., 2021) was used to prepare the starting structures for the five (Glu-Lys)25 IDP sequence variants and eight naturally occurring IDPs that were used in Frembgen-Kesner et al. (2015), which are Alzheimer amyloid β (1–40) (Aβ (1–40)) (Danielsson et al., 2002), suppressor of Mec1 lethality (Sml1) (Danielsson et al., 2008), Lotus japonicas intrinsically disordered protein 1 (LjIDP1) (Haaning et al., 2008), prothymosin α (ProTα) (Yi et al., 2007), abscisic acid stress ripening 1 (ASR1) (Goldgur et al., 2007), yeast nucleoporin 116 (Nup116) (Krishnan et al., 2008), α-synuclein (Uversky et al., 2001), and cystic fibrosis transmembrane conductance regulator regulatory region (CFTR R) (Baker, 2009). We have used Alphafold to prepare the starting structures for the naturally occurring IDPs since they have conditionally folded regions that have confident per-residue confidence scores (pLDDT) (above 70 in a range from 0 to 100), which are expected to be accurately predicted by Alphafold (Alderson et al., 2022). We have also used Alphafold for the five (Glu-Lys)25 IDP sequence variants since peptides composed of many Glu and Lys residues favor forming α-helical structures (Marqusee and Baldwin, 1987; Iqbalsyah and Doig, 2005; Meuzelaar et al., 2016; Wolny et al., 2017), and Alphafold had yielded α-helical structures for all five IDP sequence variants.

Figures 1, 2 show the amino acid sequences of these systems, respectively, and Table 1 summarizes the various characteristics of the systems obtained from the classification of intrinsically disordered ensemble regions (CIDER) program (Holehouse et al., 2015).

FIGURE 1

The five (Glu-Lys)25 IDP sequence variants used in the study. Glutamic acid (E) is colored in red for negative charge, and lysine (K) is colored in blue for positive charge. The labels for the sequence variants (sv) are from Das and Pappu (2013). The five sequence variants are the same ones tested in McCarty et al. (2019).

FIGURE 2

The eight naturally occurring IDPs used in the study. Glutamic acid (E) and aspartic acid (D) are colored in red for negative charge, and lysine (K) and arginine (R) are colored in blue for positive charge. The eight IDPs are from Frembgen-Kesner et al. (2015).

TABLE 1

Summary of classification of intrinsically disordered ensemble regions (CIDER) (Holehouse et al., 2015) results for the IDPs used in the study. NCPR denotes the net charge per residue, FCR denotes the fraction of charged residues, and κ denotes the measure of charge segregation from Das and Pappu (2013). Hydrophathy measures how hydrophobic the sequence is (0–9 with 0 being least hydrophobic and nine being most hydrophobic) (Kyte and Doolittle, 1982) and disorder measures the fraction of disorder promoting residues (Uversky, 2002). The categorization of each IDP is determined from the Das-Pappu phase diagram (Das and Pappu, 2013; Holehouse et al., 2015).

IDP	Length	NCPR	FCR	κ	Hydropathy	Disorder	Category
sv10	50	0.000	1.000	0.083	0.800	1.000	Strong polyampholytes
sv15	50	0.000	1.000	0.135	0.800	1.000	Strong polyampholytes
sv20	50	0.000	1.000	0.272	0.800	1.000	Strong polyampholytes
sv25	50	0.000	1.000	0.528	0.800	1.000	Strong polyampholytes
sv30	50	0.000	1.000	1.000	0.800	1.000	Strong polyampholytes
Aβ (1−40)	40	-0.075	0.225	0.211	4.558	0.600	Weak polyampholytes
Sml1	104	-0.048	0.221	0.143	3.712	0.635	Weak polyampholytes
LjIDP1	107	0.009	0.271	0.174	3.890	0.729	Janus sequences
ProTα	109	-0.394	0.578	0.424	2.507	0.881	Strong polyelectrolytes
ASR1	115	-0.017	0.383	0.100	3.326	0.809	Strong polyampholytes
Nup116	126	0.040	0.040	0.278	3.709	0.762	Weak polyampholytes
α-Synuclein	140	-0.064	0.279	0.172	4.097	0.729	Janus sequences
CFTR R	190	-0.026	0.289	0.285	3.743	0.679	Janus sequences

The protonation states were assigned using PROPKA 3 (Olsson et al., 2011; Søndergaard et al., 2011) at pH 7.0 for the five (Glu-Lys)25 IDP sequence variants and at appropriate pH’s for the eight IDPs as done in Frembgen-Kesner et al. (2015), which are listed in the Supplementary Material. PDB2PQR 3.4 (Jurrus et al., 2018; Unni et al., 2011; Dolinsky et al., 2007, 2004) was used to convert the PDB files to PQR format for the BD simulations. The temperature T was set to 298 K, and the dielectric constant was set to 78.4 for all systems. For the five (Glu-Lys)25 IDP sequence variants, the ionic concentration was set to NaCl 15 mM (reference concentration) or NaCl 125 mM (excess salt concentration) as done in Das and Pappu (2013) by setting the appropriate Debye length λ using Equation 1 where ϵ 0 is the permittivity of the free space, ϵ is the dielectric constant (of water in this case), k is the Boltzmann constant, T is the temperature (298 K in this case), e is the elementary charge, N is Avogadro’s constant, and C is the ionic strength in mol/m3 units. The Debye length λ was set to be 7.85 for NaCl 15 mM (reference concentration) and 2.72 for NaCl 125 mM (excess salt concentration).

2.2 Brownian Dynamics Simulations

The BD simulations were run using Browndye 2.0 (Huber and McCammon, 2010) with the spline-based potential coarse-grained force field for proteins (COFFDROP) (Andrews and Elcock, 2014; Frembgen-Kesner et al., 2015), which was newly implemented for Browndye 2.0. In COFFDROP, each amino acid is represented as a “bead” so that a protein sequence can be represented as a flexible “chain” composed of beads. In addition, since the scaling of non-bonded interactions improved COFFDROP’s ability to reproduce experimental results (Frembgen-Kesner et al., 2015), this feature was also implemented for Browndye 2.0. Moreover in Browndye 2.0, interactions can be computed less frequently, which is useful since computing these interactions take up most of the simulation time for longer chains. Finally in Browndye 2.0, a constant time step size can be set, and the recommended value is 0.05 ps for COFFDROP chains, unless bond constraints are used in which case a larger constant time step size is allowed, which can make the simulations run faster.

For the eight naturally occurring IDPs, the maximum number of BD simulation steps was set to 80,000,000, and a constant time step size of 0.05 ps was used (no bond constraints used). To calculate the hydrodynamic radius (R ) for each COFFDROP potential with a scaling factor (0.5–1.0 in intervals of 0.1) for non-bonded interactions, ten trajectories were run for each system and potential, and simulation snapshots were recorded every 200,000 steps. Hydrodynamic interactions were updated every 400 steps. The specific parameter values follow the parameter values from Frembgen-Kesner et al. (2015) since these IDPs were used to check the COFFDROP implementation for Browndye 2.0.

For the five (Glu-Lys)25 IDP sequence variants, the maximum number of BD simulation steps was set to 125,000,000, and by using bond constraints, a constant time step size of 0.2 ps was used (25 μs total). The autocorrelation functions of N-terminus C to C-terminus C distance, N-terminus C to middle C distance, and middle C to C-terminus C distance were measured and plotted with new Browndye 2.0 functions chain_atom_distances and autocor to check whether the simulation time was sufficiently long enough to obtain converged properties. As seen in Figure 3, the three autocorrelation functions converge, and the simulation time was regarded to be sufficiently long enough. The rest of the autocorrelation functions are included in the Supplementary Material. As seen from the Supplementary Material, the shortest simulation was 12 ns, whereas the longest simulation was 25 µs. To calculate structural properties such as radius of gyration (R ), ten trajectories were run for each system, and simulation snapshots were recorded every 100,000 steps. Hydrodynamic interactions were updated every 400 steps.

FIGURE 3

Autocorrelation functions of N-terminus C to C-terminus C distance, N-terminus C to middle C distance, and middle C to C-terminus C distance for (Glu-Lys)25 sv = 10 at NaCl 15 mM (reference concentration) with scaling factor 0.786. The autocorrelation functions from one of the ten trajectories are shown for clarity and the rest are included in the Supplementary Material.

3 Results

3.1 Eight Naturally Occurring IDPs

We first investigated the eight naturally occurring IDPs to see if Browndye 2.0 can reproduce the COFFDROP results in Frembgen-Kesner et al. (2015). In particular, we measured the hydrodynamic radius (R ) for each system and COFFDROP potential with a scaling factor (0.5–1.0 in intervals of 0.1) for non-bonded interactions. R is the radius of a hard-sphere that diffuses at the same rate as solute and is dependent on the size and hydration of protein. The Kirkwood definition (Kirkwood, 1996) was used to calculate R as stated in Equation 2 where r denotes pairwise distances between C of amino acids i and j, as done in Nygaard et al. (2017). R was calculated for each simulation snapshot (every 200,000 steps), and the final R value for each simulation was obtained by averaging the R values from the simulation. The average R values, along with standard error bars (95% confidence interval), from ten independent simulations, are plotted in Figure 4. To match up with the COFFDROP results that used the HYDROPRO program (Ortega et al., 2011), the average R values and standard error bars were multiplied by 1.186 and added by 1.03 as done in Nygaard et al. (2017). The R values are in good agreement with those in Frembgen-Kesner et al. (2015), which are marked as dashed lines with square markers in Figure 4, indicating that the COFFDROP implementation in Browndye 2.0 is reliable. The small discrepancies between the two results could be from the long-range electrostatic interactions being computed differently, i.e., Frembgen-Kesner et al. (2015) used a treecode algorithm (Li et al., 2009) that involves Taylor expansion to compute particle-cluster interactions, whereas this study used pairwise summations of potentials evaluated by a cubic spline using tabulated COFFDROP potential data. Except for ProTα, the ideal scaling factor for the naturally occurring IDPs is between 0.7 and 0.8, which allows the BD simulation results to match up with experimental values. We can consider ProTα to be an outlier among the eight naturally occurring IDPs since it substantially has more like-charged residues (i.e., positively charged residues aspartic acid (D) and glutamic acid (E)) as seen in Figure 2 and as noted in Frembgen-Kesner et al. (2015). The averaged ideal scaling factor, after leaving ProTα out as an outlier, is 0.786, which is slightly different from the scaling factor in Frembgen-Kesner et al. (2015) (0.825). This scaling factor was used for subsequent COFFDROP BD simulations of the five (Glu-Lys)25 IDP sequence variants. Finally, Figure 4 shows that R generally increases with sequence length, except for ProTα that is shorter than ASR1, Nup116, α-synuclein, and CFTR R.

FIGURE 4

Hydrodynamic radius (R ) values for each of the eight naturally occurring IDPs with different scaling factors for non-bonded interactions. The experimental R values are marked as dashed straight lines and correspond to the IDP with the same color in each graph. The approximate R values from Frembgen-Kesner et al. (2015) are marked as dashed lines with square markers and correspond to the IDP with the same color in each graph. The ideal scaling factor value would be where R matches with the experimental value.

3.2 Five (Glu-Lys)25 IDP Sequence Variants

We then investigated five (Glu-Lys)25 IDP sequence variants for rates of association, which were model IDP systems in Das and Pappu (2013), Sawle and Ghosh (2015), and McCarty et al. (2019). These block polymers of glutamate and lysine residues with different patterns serve as model IDPs since IDPs mostly consist of oppositely charged residues (i.e., they are polyampholytes) and do not have significant secondary structures.

We first measured the radius of gyration (R ), which serves as an indicator of protein structure compactness, i.e., the smaller the R , the tighter the packing of the protein is. R was calculated for each simulation snapshot (every 100,000 steps), and the final R value for each simulation was obtained by averaging the R values from the simulation. The average R values, along with standard error bars (95% confidence interval), from ten independent simulations, are plotted in Figure 5. As observed in Das and Pappu (2013), R generally decreases as κ, which represents the measure of charge segregation (Das and Pappu, 2013), increases. The R values are smaller than those from Das and Pappu (2013), all within the value for classical Flory random coils (∼18 ) and compact globules (∼11 ). The R values never reach near the value for self-avoiding random walks (∼28 ), which is expected for well-mixed sequence variants or those with low κ values. This is most likely attributed from using different force fields and potentially shows the limitation for the COFFDROP potential in modeling highly charged systems. However, when using the averaged ideal scaling factor for IDPs (0.786), the R values increase, show closer to expected R values, and its minimum R range match with that in Das and Pappu (2013). As κ → 1, the R values get closer to the value for compact globules (∼11 ) (Dima and Thirumalai, 2004). Finally, the R values increase as the salt concentration increases due to long-range interactions getting weaker, which is consistent with the results from Das and Pappu (2013). Overall, we were able to observe correct trends for R for the five (Glu-Lys)25 IDP sequence variants using the COFFDROP potential.

FIGURE 5

Radius of gyration (R ) values for each of the five (Glu-Lys)25 IDP sequence variants. Since each (Glu-Lys)25 IDP sequence variant has a different κ value, which represents the measure of charge segregation (Das and Pappu, 2013), κ was used as one of the axes of the graph.

We then measured the interresidue distances between residue i and residue j (R ) against residue separations |j − i|, which can characterize local concentrations of chain segments within the IDP (Das and Pappu, 2013). Specifically, the distance between residue i’s C and residue j’s C was measured. The scaling factor was set to the averaged ideal scaling factor of 0.786. R was calculated for each simulation snapshot (every 100,000 steps), and the final R value for each simulation was obtained by averaging the R values from the simulation. The average R values from ten independent simulations are plotted in Figure 6. The R values follow similar trends as observed in Das and Pappu (2013), and the concave upward parts show indications of long-range interactions between oppositely charged blocks. As observed for the R values, the R values were also smaller than those from Das and Pappu (2013), which could be attributed from using different force fields. Finally, the R values also increase as the salt concentration increases due to long-range interactions getting weaker, which is consistent with the results from Das and Pappu (2013). The effects of the salt concentration are the smallest for sv10, which has the most well-mixed sequence in comparison with the rest and can counterbalance electrostatic repulsions and attractions (Das and Pappu, 2013). Overall, we were also able to observe correct trends with R for the five (Glu-Lys)25 IDP sequence variants using the COFFDROP potential.

FIGURE 6

Interresidue distances between residue i and residue j (R ) against residue separations |j − i| for each of the five (Glu-Lys)25 IDP sequence variants.

Finally, we measured a property related to molecular association, namely the entanglement indices, or the average C distances to the other chain, between the five (Glu-Lys)25 IDP sequence variants. Since all possible pair combinations were tested, 15 simulations were run for each salt concentration (15 and 125 mM, respectively). The scaling factor was set to the averaged ideal scaling factor of 0.786. The pairwise simulations start with two IDP sequence variants oriented crosswise and translated 15 apart. All five IDP sequence variants have their middle C centered at (0.0, 0.0, 0.0). The middle C ’s of the two IDP sequence variants in the pairwise simulation were restrained to be less than 20 apart. The entanglement indices were measured using the entanglement_index, which is a new implementation in Browndye 2.0. The entanglement index was calculated for each simulation snapshot (every 100,000 steps), and the final entanglement index for each simulation was obtained by averaging the entanglement indices from the simulation. The average entanglement indices, along with standard error bars (95% confidence interval), from ten independent simulations are listed in Table 2. The entanglement indices were similar across all possible IDP sequence variant pairs, indicating that there is no direct relation between entanglement indices and charge segregation κ. This may be from all sequence variants having relatively similar R values or degrees of compactness, however, which could be from using the COFFDROP potential. All entanglement indices, with the exception for two pairs, increased, however, at excess salt concentration. This is consistent with long-range interactions getting weaker at increasing salt concentration, and published work demonstrating that long-range interactions accelerate protein-protein encounter for IDPs (Chu et al., 2012; Ganguly et al., 2013; Pang and Zhou, 2016; Tsai et al., 2016; Chu et al., 2017; Yang et al., 2019).

TABLE 2

Summary of entanglement index values of all possible pair combinations between the five (Glu-Lys)25 IDP sequence variants.

IDP #1	IDP #2	Entanglement index (A°) at 15 mM	Entanglement index (A°) at 125 mM
sv10	sv10	26.72 ± 0.98	26.38 ± 0.58
sv10	sv15	25.93 ± 0.63	27.12 ± 0.39
sv10	sv20	26.54 ± 0.98	26.70 ± 0.49
sv10	sv25	26.72 ± 0.71	26.92 ± 0.50
sv10	sv30	25.43 ± 1.61	26.70 ± 1.21
sv15	sv15	26.72 ± 0.53	27.68 ± 0.44
sv15	sv20	25.72 ± 0.69	27.29 ± 0.18
sv15	sv25	26.58 ± 0.75	26.87 ± 1.10
sv15	sv30	26.17 ± 0.28	27.98 ± 0.81
sv20	sv20	25.48 ± 0.31	25.78 ± 0.56
sv20	sv25	26.89 ± 0.83	25.85 ± 0.57
sv20	sv30	26.07 ± 0.31	26.38 ± 0.70
sv25	sv25	27.00 ± 0.25	27.44 ± 0.58
sv25	sv30	25.07 ± 3.05	27.15 ± 0.40
sv30	sv30	24.17 ± 2.38	26.59 ± 0.81

4 Discussion and Conclusion

We have presented our new COFFDROP force field implementation on Browndye 2.0 that enabled us to study IDPs computationally. BD simulations are ideal to study large-scale motions with longer time scales and diffusion-limited molecular associations, including the aggregation of IDPs. We have presented results that show that our COFFDROP implementation is reliable to study naturally occurring IDPs. We have also studied model (Glu-Lys)25 IDPs using our COFFDROP implementation and found that there is no relation between entanglement indices and how well the charges are mixed and segregated within the IDPs. However, this may be from the limitation of the COFFDROP potential in studying highly charged systems, which was also noted in Frembgen-Kesner et al. (2015). The COFFDROP potential was derived from MD simulations of all possible amino acid pairs (Andrews and Elcock, 2014; Frembgen-Kesner et al., 2015), but the simulations did not include salt so the COFFDROP potential may be limiting in modeling systems with strong charge-charge interactions.

For future work, we plan to implement a program to measure the rates of association with an appropriate reaction criterion as done in Ganguly et al. (2013), Liu et al. (2019). Then we plan to measure the rates of association of a highly positive IDP binding to a highly negative IDP (i.e., oppositely charged IDPs), an interaction that may be abundant in eukaryotes for regulation (e.g., cellular localization) (Borgia et al., 2018). We also plan to look at the rates of association between IDPs and folded proteins with secondary structures (Ruff et al., 2019). However, since COFFDROP is meant to model IDPs or systems without significant secondary or tertiary structures, the secondary structural elements would need to have constraints to have them fixed throughout the simulation, and the folded protein would be treated as a rigid body. The IDP would still be modeled as a flexible chain, and the scaling factor of 0.786 would be used for the simulation.