IDPConformerGenerator: A Flexible Software Suite for Sampling the Conformational Space of Disordered Protein States

The power of structural information for informing biological mechanisms is clear for stable folded macromolecules, but similar structure–function insight is more difficult to obtain for highly dynamic systems such as intrinsically disordered proteins (IDPs) which must be described as structural ensembles. Here, we present IDPConformerGenerator, a flexible, modular open-source software platform for generating large and diverse ensembles of disordered protein states that builds conformers that obey geometric, steric, and other physical restraints on the input sequence. IDPConformerGenerator samples backbone phi (φ), psi (ψ), and omega (ω) torsion angles of relevant sequence fragments from loops and secondary structure elements extracted from folded protein structures in the RCSB Protein Data Bank and builds side chains from robust Monte Carlo algorithms using expanded rotamer libraries. IDPConformerGenerator has many user-defined options enabling variable fractional sampling of secondary structures, supports Bayesian models for assessing the agreement of IDP ensembles for consistency with experimental data, and introduces a machine learning approach to transform between internal and Cartesian coordinates with reduced error. IDPConformerGenerator will facilitate the characterization of disordered proteins to ultimately provide structural insights into these states that have key biological functions.

S2 chains derived from protein coordinates. During the preparation of the torsion angle database, the first and last residues of each protein are discarded to reduce artifacts in torsion angles. To create this database, IDPConformerGenerator has a series of commands that can be used in a modular fashion. First, the PDB downloader command (`pdbdl`) downloads a list of PDB IDs from the RCSB databank. Both PDB and mmCIF format files are accepted. Listed PDB IDs have the form of ABCDX where ABCD is the PDB ID and X is the optional chain id. If given, only that chain is stored in disk. We recommend using a culled list of PDB IDs from the Dunbrack PISCES database (http://dunbrack.fccc.edu/lab/pisces) 1 . The PDB downloader client further processes each PDB file by removing known solvent and ligands, considering only elements composed of amino acids, considering only the alternative location 'A' or ' ' and selecting only the first model for multi-model structures. If no chain ID is provided, IDPConformerGenerator will save each chain of the PDB file separately, and treat them as individual structures in subsequent steps. In this step, the user can easily add or remove PDBs to feed the sampling process. PDB files can also be added manually to the set of downloaded PDBs.
After the PDB files have been downloaded, the secondary structure calculator client (`sscalc`) annotates the secondary structure information from the PDBs using third-party software. We have implemented parsers for DSSP 2,3 , but IDPConformerGenerator code structure allows for facile implementation of additional parsers. Because DSSP v2 and v3 do not consider polyproline type II helices, we adapted the DSSP-PPII Perl script 4 to work with DSSP v2 or v3 annotate poly-proline type II helices as "P". The output of the `sscalc` client produces a JSON file containing aligned primary sequences, residue numbers and DSSP annotations. Additionally, the input PDB files are split according to chain gaps, if existent, and resulting chains treated as individual PDBs. Here, segments shorter than 3 (configurable) are discarded. We suggest using the "reduce" flag (`-rd`), which directs IDPConformerGenerator to group DSSP codes into either Loops (L), Helices (H), or Extended strands (E). IDPConformerGenerator's definition of a loop incorporates DSSP's hydrogen-bonded turn (T), bend (S), beta-bridge (B), loops ('space'), poly-proline type II helix (P), and pi-helix (I); helices include DSSP's alpha-helix (H), and 3-10 helix (G); extended strands (E) are identical to DSSP's extended strands that participate in beta-ladder structures 5 .
Finally, the torsion angles for each PDB are calculated using the `torsions` command in IDPConformerGenerator. The torsion angle calculator client calculates phi, psi, and omega torsion angles from the coordinates. On top of all of IDPConformerGenerator's filters in previous steps, the integrity of the PDB files is further confirmed at this stage, most importantly, by ensuring input structures have all required backbone atoms and no gaps. The torsion angle calculator client stores all the successfully processed PDB IDs aligned with their secondary structure code, primary sequence, and per-residue torsion angles in a JSON file. This is the initial IDPConformerGenerator database that will be used during the building algorithm. We chose a JSON structure because the amount of data used in this step, despite being large, is still workable within a JSON file, and because we advocate for using text-editable files as much as possible in research software pipelines for the sake of traceability, editing, and easy reviewing.
Once the database file has been created, it can be used to build conformers of any protein sequence, and can be shared among collaborators. The database only needs to be recreated if different sets of structures are desired (e.g., including new structures or using different resolution cutoffs). Otherwise, users only need the primary sequence of the disordered protein of interest. If users want to bias secondary structure propensities by NMR chemical shifts, they need an NMR STAR v2 file 6 for CheSPI 7 or SHIFTY 8 format for d2D 9 . The outputs of secondary structure propensity predictions from CheSPI or d2D can be directly fed into the custom secondary structure sampling conversion `csssconv` client which will read the fractional secondary structure information and standardize them for each secondary structure code on a per-residue basis for the conformer building module. The standardized CSSS input file can also be edited in the `makecsss` module. If no chemical shift-based secondary structure propensity prediction exists, however, the user can specify secondary structure biases using `makecsss`, as well. The fact that the database preparation process is separated into three commands allows advanced users to manually alter PDB files, add or remove them from the JSON files, or otherwise intervene in the different steps, if needed.
Building conformers with IDPConformerGenerator. IDPConformerGenerator builds coordinates by sampling torsion angles from the database file. Therefore, the user has complete control over which PDBs will be used to drive conformational sampling. IDPConformerGenerator's `build` command-line interface requires the database file (created in the previous steps) and the IDP target sequence. Several details of the building process can be configured through the client's parameter flags. The most important will be explained here, and the complete list is available in the documentation pages, or via the `IDPConformerGenerator build -h` menu. Before starting building the actual coordinates, IDPConformerGenerator maps the input sequence to the respective positions in the input database. For this, the input protein sequence is split into fragments of different sizes (described below). Identified fragments allow overlap, i.e., abcde has three fragments of 3 residues: abc, bcd, and cde. Then, IDPConformerGenerator searches the input database for positions matching the fragments generated from the input protein. IDPConformerGenerator considers fragments preceding proline differently, i.e., abcP is treated differently than abcX, where X is not P 10,11 . Users can also specify a list of residues allowing substitutions to increase the likelihood of matching. For example, A → AV allows the input fragment GAG to match GAG or GVG sequences in the database. This process happens internally and no additional files are saved to disk.
Users specify the fragment sizes utilized to build conformers. Different from other tools to build disordered chains, IDPConformerGenerator allows the user to select whatever sizes for the fragments and the probabilities by which they are sampled during the building process. By default, IDPConformerGenerator uses fragment sizes 1, 2, 3, 4, and 5 with relative probabilities of 10%, 10%, 30%, 30%, and 20%, respectively. Importantly, the user can sample fragment torsion angles specifically for user-defined secondary structure annotations. For example, fragments can be searched for residues in the database annotated as "loop" and/or "sheet" and/or "helix" regions depending on the user choice, using the `-rd` flag in the `sscalc` client. With these options, fragments consist only of residues matching a single secondary structure annotation. Alternatively, the user can specify "--dany" to allow fragment matching to disregard the secondary structure annotation. Moreover, users can specify sampling specific secondary structure annotations for distinct parts of the input protein and in certain probabilities. This option allows, for example, conformers to directly sample alpha-helices for certain residue ranges of the sequence, matching experimental data. To benefit from this feature, users need to prepare the required CSSS files as noted in the previous section. These options allow a rich combination of sampling possibilities adaptable to the project requirements. S4 Generation of structures requires transformation from internal to Cartesian coordinates. The default approach used for this work is based on an internal database built from structures with resolutions better than 1.6Å created with the bgeo command line interface. The database is hierarchically organized with bond angles separated by type, by residue type and residue neighbors, placed in 10° bins. While building backbones, the angles are statistically sampled from these values stored in the bgeo database. The bond lengths are taken as fixed values for each backbone atom pair and residue type, and are calculated from the same structures.
IDPConformerGenerator also ships with the recently developed Int2Cart 12 recurrent neural network algorithm for transformation from internal to Cartesian coordinates, currently integrated as a flag `-bgeo-strategy int2cart` within IDPConformerGenerator. This deep learning model provides improved prediction of bond lengths and bond angles from backbone ω, φ, and ψ torsion angles and residue types of a protein, relative to using fixed bond lengths and bond angles, or the original IDPConformerGenerator statistical approach used for most ensemble calculation in this work. Running Int2Cart within IDPConformerGenerator requires CUDA compatible GPUs, with at least two good GPUs needed for optimal performance if it is utilized at each fragment incorporation step. More comprehensive benchmarks and optimization of Int2Cart are forthcoming, with the goal to exploit the improved geometry that Int2Cart provides with no speed compromises. In addition to being implemented within IDPConformerGenerator, the Int2Cart algorithm is also available as an individual Python package at https://github.com/THGLab/int2cart.
IDPConformerGenerator can create backbone-only or full sidechain-containing conformers using the recommended Monte Carlo Side Chain Entropy (MC-SCE) algorithm 13 or FASPR 14 . The rationale during the process is however the same: conformers are accepted if the final energy is below the threshold. Accepted conformers are saved to disk. The conformer building process is an embarrassingly parallel process, for which users can give a number of cores across which the conformers will be distributed: 1000 conformers in 10 CPUs results in each CPU building 100 conformers, simultaneously. Also, IDPConformerGenerator is deterministic, and the building process can be reproduced by giving the same input and the same seed number.
Utilities available and technical information. As IDPConformerGenerator has been designed with modularity in mind, certain functions could be used as Python libraries for scripting or integration into other programs. One example of this is the parser for PDB files. Since IDPConformerGenerator parses only the labels and coordinate data from PDB and PDBx/mmCIF files, its parser module is much faster at processing structure files than other Python tools, including the highly utilized BioPython, or alternatives such as pdb2sql. IDPConfGen is approximately 2.5x faster than BioPython and 1.4x faster than pdb2sql for processing PDB files and approximately 847x faster than BioPython for mmCIF files, testing the same 22,599 structures with variable protein lengths and complexities in both PDB and mmCIF formats (without mulit-processing). BioPython had 48 errors due to these PDB IDs having only mmCIF file formats, while IDPConformerGenerator had no errors. Note that pdb2sql does not handle the mmCIF file format. Therefore, IDPConformerGenerator's PDB parsing algorithm could be imported as a library for custom scripts and computational processing of structure files.
Other examples are tools for analysis of torsion angles and fractional secondary structure of ensembles and for analysis of the database to provide histograms of matching fragment sequences and to search PDB headers for keywords of interest. For a basic analysis of conformer ensembles, there are three plotting functions within the "torsions" and "sscalc" sub-clients. Within "torsions", given a folder/tarball of generated PDBs, a torsion angle distribution scatter plot can be plotted as well as fractional secondary structure information based on Ramachandran ranges for Phi/Psi angles. Fractional secondary structure based on DSSP codes can be plotted in the "sscalc" sub-client.
IDPConformerGenerator is written in Python, its source is open and can be found at Julie Forman-Kay's GitHub repository (https://github.com/julie-forman-kay-lab/IDPConformerGenerator). While IDPConformerGenerator is written all in Python, we incorporated a modification of the C++ code for the FASPR sidechain building tool in our repository to facilitate its usage (license compatible). For further details, please refer to the usage examples on the ReadTheDocs webpage (https://IDPConformerGenerator.readthedocs.io), as well as within the help documentation in IDPConformerGenerator, providing a guide to users of all levels from installation on multiple operating systems, to updating, and to using IDPConformerGenerator along with all of its functions.

Monte Carlo Side Chain Entropy (MC-SCE).
The original code for MC-SCE 13 was re-written and implemented in Python, incorporating elements of the IDPConformerGenerator code. For each side chain other than glycine and alanine, the candidate torsion angle rotamers and their prior probabilities are obtained from the Dunbrack library 1 with the specific backbone / values rounded to the closest tenth. For the N-terminal residue that does not have and the C-terminal residue that does not have , −180° was used as the undefined backbone torsion angle. The torsion angles ) and * were extended by a standard deviation, i.e. + ± were also considered as rotamer candidates in addition to + itself. Accordingly, probabilities for all rotamers from residues that only have ) were set to 1/3 of their original probabilities in the Dunbrack library, and those from residues having and * were set to 1/9 of the original probabilities. Torsion angles beyond * were not expanded.
At each side chain growing position, let { 0 } denote all candidate rotamers, the probability for each rotamer state to be selected was defined by: where (67+87) ( 0 ) was the prior probability for the rotamer 0 from the Dunbrack library scaled by 1/3 or 1/9, ( 0 )was the energy of the peptide structure with partially completed side chains up until the current building position, and = ) 3 D E where G was Boltzmann's constant and was selected as 300K. In this study only the 12-6 Lennard Jones term and a pairwise distance-based clash term were used for energy evaluation. The Lennard Jones term followed the implementation of Amber14 15 force field, and the clash term set the energy to ∞ if the distance for any pair of atoms in the structure was smaller than 0.6 times the sum of the van der Waals radii of the two atoms. We intentionally kept minimal terms in the energy function to avoid computational burden. If the clashes were inevitable, the algorithm selects the lowest energy conformation. Further details can be found in the MCSCE project (https://github.com/THGLab/MCSCE).
While it is possible to run MC-SCE through IDPConformerGenerator as a side chain method (`-scm`), for longer protein sequences we would advise running MC-SCE separately after the backbones have been successfully generated by IDPConformerGenerator. This is due to the current integration of MC-SCE within IDPConformerGenerator requiring MC-SCE to run immediately after each successful generation of a backbone. Although both programs are capable of Pythonic multiprocessing, running MC-SCE with multiple trials within IDPConformerGenerator will significantly stress the single worker responsible for calculating that conformer. For example, when requesting the typical 128 trials for the simple mode, if we specified 32 workers in IDPConformerGenerator, each of those workers would be independently generating backbones which then require 128 sidechain trials, leading to 4096 possible concurrent processes. In contrast, running MC-SCE with the same settings once all backbones have been generated would distribute the 128 trials among 32 workers, each worker running 4 calculations as opposed to 128.
Calculating Ensembles with FastFloppyTail. In order to use FastFloppyTail, input files are required to generate fragment libraries. We used the input files for Sic1 and alpha-synuclein provided by the authors of FastFloppyTail 16 . For the drkN SH3 domain, the sequence used by the authors of FastFloppyTail was different (one residue substitution) so we regenerated all the input files. We also created our own for inhibitor-2 and Tau since those were not modelled in the FastFloppyTail study. For proteins for which we needed to create input files, we used the suggested "best-reweighting scheme" to reweight the secondary structure predictions generated with RaptorX 17 and PSIPRED 18 . The disorder prediction was done using RaptorX, since this is the only predictor the authors found to predict alphasynuclein as disordered. We created fragment libraries for all sequences using these input files and the Quota protocol, with the times required to create these libraries provided in Supplementary Table 2C. For calculating the structures, we used the 3-mer fragment libraries, ignoring the 9-mer fragment libraries. We generated ensembles with FastFloppyTail for all sequences other than the drkN SH3 domain using the suggested approach using the sequence and the fragment library as inputs, without supplying a disorder prediction file so that the sequence is treated as fully disordered, since RaptorX predicted all sequences other than the drkN SH3 domain to be largely disordered. For the unfolded state of the drkN SH3 domain, we ran two calculations, one treating it as partially ordered and disordered as RaptorX predicts, using the default protocol for a sequence with more than 20% predicted order as 42% of residues were predicted to be ordered. This entails first running AbInitioVO and providing the lowest energy output to use in FastFloppyTail along with the disordered probability prediction from RaptorX so that any ordered regions are maintained. For the second set of calculations for the unfolded state of the drkN SH3 domain, we treated it as fully disordered. For these, we edited the disorder probability prediction file from the RaptorX server, changing all residues predicted to be ordered to disordered with probability 1.0. This edited disorder file was used to reweight the secondary structure predictions with the best-reweighting scheme and to then generate the fragment libraries. AbInitioVO was not utilized and we did not supply FastFloppyTail with the disorder prediction file so the sequence was treated as fully disordered.

Analysis of Conformational Ensembles and Experimental Data Sources.
To analyze the generated ensembles for fractional secondary structure information, torsion angle distribution, and tertiary contacts, well established Python libraries were used from MDAnalysis 19,20 and MDTraj 21 . For the analysis of pairwise RMSDs, the `diffusionmap` and `align` functions were used from the MDAnalysis `analysis` library. To calculate secondary structure information as DSSP codes, we used the `compute_dssp` function in MDTraj. Alpha, beta, and coil regions on the Ramachandran plot were S7 determined using the `Ramachandran` function in the MDAnalysis `dihedrals` library. Tertiary contacts (Cα-Cα) analysis was done with the `distances` function in the `analysis` library of MDAnalysis. Statistics performed for tertiary contacts analysis used the `stat` function from the `stats` library in SciPy.
RMSD values to experimental restraints and scores were calculated using ENSEMBLE 22 and X-EISD 23,24 . ENSEMBLE was modified to generate human-readable back-calculated data and to enable calculation of ENSEMBLE scores without optimizing. NMR and SAXS experimental data for comparison to ensembles were taken from supplemental information of papers or from our previous publications, and, for some chemical shifts, from the Biological Magnetic Resonance Data Bank (BMRB 25 ). For the drkN SH3 domain unfolded state and I-2, data were from Marsh & Forman-Kay (2012) 26 , with I-2 chemical shifts also in BMRB (ID 15179). For Sic1, data were from Gomes et al (2020) 27 , with chemical shifts also in the BMRB (ID 16657). For alpha-synuclein, RDC, 3 J-coupling, chemical shifts, and PRE restraints were from Ferrie et al (2020) 16 and SAXS restraints were from Pesce & Lindorff-Larsen (2021) 28 . For Tau chemical shifts, BMRB IDs: 17920 and 50701 were used and for SAXS, data were also from Pesce & Lindorff-Larsen (2021) 28 . Chemical shifts for the drkN SH3 domain unfolded state and I-2 were also used for d2D 9 and CheSPI 7 calculations.
Analysis of MD Ensembles. Conformations from the MD simulations of the unfolded state of the drkN SH3 domain, Sic1, and aSyn generated using the a99SB-disp 29 force field were extracted using taurenmd 30 . Every frame was extracted for the secondary structure profile analysis while every 29 th frame for drkN SH3 and aSyn and every 166 th frame for Sic1 were extracted to give at least 1000 conformers for the calculation of ENSEMBLE scores and RMSDs against experimental data using previously described approaches. Renumbering residues of Sic1 for agreement with experimental data formatting was done with pdb-tools 31 . S8 SUPPLEMENTARY FIGURES. All plots, in the main text and Supplement, were generated using Matplotlib 32 unless otherwise specified. Figure S1. Histogram of omega dihedral angle distributions for structures found within the IDPConformerGenerator database. PDBs from the 24,003 PDBID database were used, with sets of this with resolutions better than 1.5Å (~5,000 structures), better than 1.8Å (~16,000 structures) and better than 2.0Å (full set). The deviations from an omega torsion angle of 180° (trans) is plotted. values are centered around 0o to facilitate visualization of the distribution. Cis peptide bonds are ignored. Distributions for omega are shown for various secondary structures as defined by the reduced DSSP codes (loops, L+; helices, H+; and extended, E+) as well as for every protein in the database (whole database). Radius of gyration (R g ), end-to-end distance (R ee ), asphericity (A) and pairwise root-mean-squareddeviations of atomic positions (pwRMSDs), are shown as a function of secondary structure sampling parameters for ensembles generated with different secondary structure sampling, including loops (L+), loops and helices (L+H+), loops, helices and extended strands (L+H+E+), all torsion angles agnostic to secondary structure (ANY) and biased by d2D chemical shifts (CSSS), or with FastFloppyTail (FFT), for Sic1 (row 1), alpha-synuclein (row 2), and Tau (row 3). The error bars represent the standard deviations. For all but Tau, 1000-conformer ensembles were generated and analyzed. For Tau, variable number of conformers were generated and analyzed, with FFT, L+, L+H+, L+H+E+, ANY, CSSS having 500, 240, 200, 500, 350, 480 conformers, respectively. S14 Supplementary Figure S7. Pairwise RMSD distance matrices. Ensembles of the drkN SH3 domain unfolded state (A, B, C) and I-2 (D, E, F) with 1000 conformers each, generated using different secondary structure sampling methods, ANY and CSSS with either CheSPI or d2D to bias sampling. differences between C alpha -C alpha distance matrixes (P < 0.05 from a Mann-Whitney U test). S15 Supplementary Figure S8. Fractional secondary structure in ensembles of the drkN SH3 domain unfolded state. Analyses were performed on sets of 1000-conformer pools of the drkN SH3 domain unfolded state generated with different combinations of secondary structure sampling consisting of loops, helix, and extended strands. Orange indicates alpha-helix for DSSP (solid) and the alpha-region on the Ramachandran diagram (dashed). Blue indicates extended strand for DSSP (solid) and betaregion on the Ramachandran diagram (dashed). Black indicates coil/loop for DSSP (solid) and other regions on the Ramachandran diagram (dashed). S16 Supplementary Figure S9. Fractional secondary structure in calculated ensembles of the drkN SH3 domain unfolded state generated with IDPConformerGenerator using ANY and CSSS and FastFloppyTail. Analyses were with ensembles of 1000 conformers, a backbone energy threshold of 100 kJ, and side chains added on with MC-SCE with 128 trials per conformer in simple mode. Sequence substitutions were from the EDSSMat50 substitution matrix rows 5, 3, and 2. * is for the standard protocol which for this case treats the protein as a mixture of ordered and disordered, while the other is for a modified protocol in which the protein is considered to be fully disordered.  Figure S14. Analysis of tertiary contacts for Sic1 ensembles with varying secondary structure sampling preferences. (top row) C alpha -C alpha distance matrixes (lower) with deviations (upper) for 1000-conformer ensembles of Sic1 generated with the loops only and ANY flag for secondary structure, with substitutions from columns 5, 3, and 2 of the EDSSMat50 amino acid substitution matrix and with variable fragment lengths. (bottom row) Significant differences between the two C alpha -C alpha distance matrixes (lower) and deviations (upper), (P < 0.05 from a Mann-Whitney U test).

S22
Supplementary Figure S15. Fractional secondary structure in calculated ensembles of the drkN SH3 domain unfolded state generated with IDPConformerGenerator (using ANY and CSSS), δ2D predictions, and molecular dynamics. Analyses were with IDPConformerGenerator ensembles of 1000 conformers, a backbone energy threshold of 100 kJ, and side chains added on with MC-SCE with 128 trials per conformer in simple mode. 29977 conformers were used from the molecular dynamics simulation of drkN SH3 from Robustelli et al 29 .

S23
Supplementary Figure S16. Fractional secondary structure in calculated ensembles of Sic1 generated with IDPConformerGenerator (using ANY and CSSS), δ2D predictions, and molecular dynamics. Analyses were with IDPConformerGenerator ensembles of 1000 conformers, a backbone energy threshold of 100 kJ, and side chains added on with MC-SCE with 128 trials per conformer in simple mode. 33328 conformers were used from the molecular dynamics simulation of Sic1 from Robustelli et al 29 .

S24
Supplementary Figure S17. Fractional secondary structure in calculated ensembles of the alpha-Synuclien generated with IDPConformerGenerator (using ANY and CSSS), δ2D predictions, and molecular dynamics. Analyses were with IDPConformerGenerator ensembles of 1000 conformers, a backbone energy threshold of 100 kJ, and side chains added on with MC-SCE with 128 trials per conformer in simple mode. 29977 conformers were used from the molecular dynamics simulation of aSyn from Robustelli et al 29 .