Ensembler: Enabling High-Throughput Molecular Simulations at the Superfamily Scale

Abstract

The rapidly expanding body of available genomic and protein structural data provides a rich resource for understanding protein dynamics with biomolecular simulation. While computational infrastructure has grown rapidly, simulations on an omics scale are not yet widespread, primarily because software infrastructure to enable simulations at this scale has not kept pace. It should now be possible to study protein dynamics across entire (super)families, exploiting both available structural biology data and conformational similarities across homologous proteins. Here, we present a new tool for enabling high-throughput simulation in the genomics era. Ensembler takes any set of sequences--from a single sequence to an entire superfamily--and shepherds them through various stages of modeling and refinement to produce simulation-ready structures. This includes comparative modeling to all relevant PDB structures (which may span multiple conformational states of interest), reconstruction of missing loops, addition of missing atoms, culling of nearly identical structures, assignment of appropriate protonation states, solvation in explicit solvent, and refinement and filtering with molecular simulation to ensure stable simulation. The output of this pipeline is an ensemble of structures ready for subsequent molecular simulations using computer clusters, supercomputers, or distributed computing projects like Folding@home. Ensembler thus automates much of the time-consuming process of preparing protein models suitable for simulation, while allowing scalability up to entire superfamilies. A particular advantage of this approach can be found in the construction of kinetic models of conformational dynamics--such as Markov state models (MSMs)--which benefit from a diverse array of initial configurations that span the accessible conformational states to aid sampling. We demonstrate the power of this approach by constructing models for all catalytic domains in the human tyrosine kinase family, using all available kinase catalytic domain structures from any organism as structural templates. Ensembler is free and open source software licensed under the GNU General Public License (GPL) v2. It is compatible with Linux and OS X. The latest release can be installed via the conda package manager, and the latest source can be downloaded from https://github.com/choderalab/ensembler.

The rapidly expanding body of available genomic and protein structural data provides a rich resource for understanding protein dynamics with biomolecular simulation. While computational infrastructure has grown rapidly, simulations on an omics scale are not yet widespread, primarily because software infrastructure to enable simulations at this scale has not kept pace. It should now be possible to study protein dynamics across entire (super)families, exploiting both available structural biology data and conformational similarities across homologous proteins. Here, we present a new tool for enabling high-throughput simulation in the genomics era. Ensembler takes any set of sequences--from a single sequence to an entire superfamily--and shepherds them through various stages of modeling and refinement to produce simulation-ready structures. This includes comparative modeling to all relevant PDB structures (which may span multiple conformational states of interest), reconstruction of missing loops, addition of missing atoms, culling of nearly identical structures, assignment of appropriate protonation states, solvation in explicit solvent, and refinement and filtering with molecular simulation to ensure stable simulation. The output of this pipeline is an ensemble of structures ready for subsequent molecular simulations using computer clusters, supercomputers, or distributed computing projects like Folding@home. Ensembler thus automates much of the time-consuming process of preparing protein models suitable for simulation, while allowing scalability up to entire superfamilies. A particular advantage of this approach can be found in the construction of kinetic models of conformational dynamics--such as Markov state models (MSMs)--which benefit from a diverse array of initial configurations that span the accessible conformational states to aid sampling. We demonstrate the power of this approach by constructing models for all catalytic domains in the human tyrosine kinase family, using all available kinase catalytic domain structures from any organism as structural templates. Ensembler is free and open source software licensed under the GNU General Public License (GPL) v2. It is compatible with Linux and OS X. The latest release can be installed via the conda package manager, and the latest source can be

Fig 1. Diagrammatic representation of the stages of the Ensembler pipeline and illustrative statistics for modeling all human tyrosine kinase catalytic domains.

On the left, the various stages of the Ensembler pipeline are shown. The red labels indicate the corresponding text description provided for each stage in the Design and Implementation section. On the right, the number of viable models surviving each stage of the pipeline is shown for the 93 target TK domains and for two representative individual TK domains (SRC and ABL). Typical timings on a computer cluster (containing Intel Xeon E5-2665 2.4GHz hyperthreaded processors and NVIDIA GTX-680 or GTX-Titan GPUs) is reported to illustrate resource requirements per model for modeling the entire set of tyrosine kinases. Note that CPU-h denotes the number of hours consumed by the equivalent of a single CPU hyperthread and GPU-h on a single GPU--parallel execution via MPI reduces wall clock time nearly linearly.

http://dx.doi.org/10.1371/journal.pcbi.1004728.g001

Target selection and retrieval (A)

The first stage entails the selection of a set of target protein sequences--the sequences for which the user is interested in generating simulation-ready structural models. This may be a single sequence--such as a full-length protein or a construct representing a single domain--or a collection of sequences, such as a particular domain from an entire family of proteins. The output of this stage is a FASTA-formatted text file containing the desired target sequences with corresponding arbitrary identifiers.

The ensembler command-line tool allows targets to be selected from UniProt--a freely accessible resource for protein sequence and functional data (uniprot.org) [20]--via a UniProt search query. To retrieve target sequences from UniProt, the subcommand gather_targets is used with the --query flag followed by a UniProt query string conforming to the same syntax as the search function available on the UniProt website. For example, --query 'mnemonic:SRC_HUMAN' would select the full-length human Src sequence, while the query shown in Box 1 would select all human tyrosine protein kinases which have been reviewed by a human curator. In this way, the user may select a single protein, many proteins, or an entire superfamily from UniProt. The program outputs a FASTA file, setting the UniProt mnemonic (e.g. SRC_HUMAN) as the identifier for each target protein.

Box 1. Ensembler command-line functions used to select targets and templates

The commands retrieve target and template data by querying UniProt. The query string provided to the gather_targets command selects all human tyrosine protein kinases which have been reviewed by a curator, while the query string provided to the gather_templates command selects all reviewed protein kinases of any species. The --uniprot_domain_regex flag is used to select a subset of the domains belonging to the returned UniProt protein entries, by matching the domain annotations against a given regular expression. In this example, domains of type 'Protein kinase', 'Protein kinase 1', and 'Protein kinase 2' were selected, while excluding many other domain types such as 'Protein kinase; truncated', 'Protein kinase; inactive', 'SH2', 'SH3', etc. Target selection simply entails the selection of sequences corresponding to each matching UniProt domain. Template selection entails the selection of the sequences and structures of any PDB entries corresponding to the matching UniProt domains.

ensembler gather_targets --query 'family:"tyr protein kinase family" AND organism:"homo sapiens" AND reviewed:yes'

--uniprot_domain_regex '^Protein kinase(?!; truncated)(?!; inactive)'

ensembler gather_templates --gather_from uniprot --query 'domain:"Protein kinase" AND reviewed:yes'

--uniprot_domain_regex '^Protein kinase(?!; truncated)(?!; inactive)'

In many cases, it will be desirable to build models of an isolated protein domain, rather than the full-length protein. The gather_targets subcommand allows protein domains to be selected from UniProt data by passing a regular expression string to the --uniprot_domain_regex flag. For example, the above --query flag for selecting all human protein kinases returns UniProt entries with domain annotations including 'Protein kinase', 'Protein kinase 1', 'Protein kinase 2', 'Protein kinase; truncated', 'Protein kinase; inactive', 'SH2', 'SH3', etc. The regular expression shown in Box 1 selects only domains of the first three types. If the --uniprot_domain_regex flag is used, target identifiers are set with the form [UniProt mnemonic]_D[domain index], where the latter part represents a 0-based index for the domain--necessary because a single target protein may contain multiple domains of interest (e.g. JAK1_HUMAN_D0, JAK1_HUMAN_D1).

Target sequences can also be defined manually (or from another program) by providing a FASTA-formatted text file containing the desired target sequences with corresponding arbitrary identifiers.

Template selection and retrieval (B)

Ensembler uses comparative modeling to build models, and as such requires a set of structures to be used as templates. The second stage thus entails the selection of templates and storage of associated sequences, structures, and identifiers. These templates can be specified manually, or using the ensembler gather_templates subcommand to automatically select templates based on a search of the Protein Data Bank (PDB) or UniProt. A recommended approach is to select templates from UniProt which belong to the same protein family as the targets, establishing homology and some degree of sequence identity between targets and templates.

The ensembler gather_templates subcommand provides methods for selecting template structures from either UniProt or the PDB (http://www.rcsb.org/pdb), specified by the --gather_from flag. Both methods select templates at the level of PDB chains--a PDB structure containing multiple chains with identical sequence spans (e.g. for crystals with non-crystallographic symmetry giving rise to independent conformations of the protein within the asymmetric unit) would thus give rise to multiple template structures.

Selection of templates from the PDB simply requires passing a list of PDB IDs as a comma-separated string, e.g. --query 2H8H,1Y57. Specific PDB chain IDs can optionally also be selected via the --chainids flag. The program retrieves structures from the PDB server, as well as associated data from the SIFTS service (www.ebi.ac.uk/pdbe/docs/sifts) [21], which provides residue-level mappings between PDB and UniProt entries. The SIFTS data is used to extract template sequences, retaining only residues which are resolved and match the equivalent residue in the UniProt sequence--non-wildtype residues are thus removed from the template structures. Furthermore, PDB chains with less than a given percentage of resolved residues (default: 70%) are filtered out. Sequences are stored in a FASTA file, with identifiers of the form [UniProt mnemonic]_D[UniProt domain index]_[PDB ID]_[PDB chain ID], e.g. SRC_HUMAN_D0_2H8H_A. Matching residues then extracted from the original coordinate files and stored as PDB-format coordinate files.

Selection of templates from UniProt proceeds in a similar fashion as for target selection; the --query flag is used to select full-length proteins from UniProt, while the optional --uniprot_domain_regex flag allows selection of individual domains with a regular expression string (Box 1). The returned UniProt data for each protein includes a list of associated PDB chains and their residue spans, and this information is used to select template structures, using the same method as for template selection from the PDB. Only structures solved by X-ray crystallography or NMR are selected, thus excluding computer-generated models available from the PDB. If the --uniprot_domain_regex flag is used, then templates are truncated at the start and end of the domain sequence.

Templates can also be defined manually. Manual specification of templates simply requires storing the sequences and arbitrary identifiers in a FASTA file, and the structures as PDB-format coordinate files with filenames matching the identifiers in the sequence file. The structure residues must also match those in the sequence file.

Template refinement (C)

Unresolved template residues can optionally be modeled into template structures with the loopmodel subcommand, which employs a kinematic closure algorithm provided via the loopmodel tool of the Rosetta software suite [22, 23]. We expect that in certain cases, pre-building template loops with Rosetta loopmodel prior to the main modeling stage (with MODELLER) may result in improved model quality. Loop remodeling may fail for a small proportion of templates due to spatial constraints imposed by the original structure; the subsequent modeling step thus automatically uses the remodeled version of a template if available, but otherwise falls back to using the non-remodeled version. Furthermore, the Rosetta loopmodel program will not model missing residues at the termini of a structure--such residue spans are modeled in the subsequent stage.

Alignment and comparative modeling (D)

In the modeling stage, structural models of the target sequence are generated from the template structures, with the goal of modeling the target in a variety of conformations that could be significantly populated under equilibrium conditions.

Modeling is performed using the automodel function of the MODELLER software package [24, 25] to rapidly generate a single model of the target sequence from each template structure. MODELLER uses simulated annealing cycles along with a minimal forcefield and spatial restraints--generally Gaussian interatomic probability densities extracted from the template structure with database-derived statistics determining the distribution width--to rapidly generate candidate structures of the target sequence from the provided template sequence [24, 25].

While MODELLER's automodel function can generate its own alignments automatically, a standalone function was preferable for reasons of programming convenience. As such, we implemented pairwise alignment functionality using the BioPython pairwise2 module [26]--which uses a dynamic programming algorithm--with the PAM 250 scoring matrix of Gonnet et al. [27], though other choices of scoring matrices available within the module can be selected. The alignments are carried out with the align subcommand, prior to the modeling step which is carried out with the build_models subcommand. The align subcommand also writes a list of the sequence identities for each template to a text file, and this can be used to select models from a desired range of sequence identities. The build_models subcommand and all subsequent pipeline functions have a --template_seqid_cutoff flag which can be used to select only models with sequence identities greater than the given value. We also note that alternative approaches could be used for the alignment stage. For example, multiple sequence alignment algorithms [28], allow alignments to be guided using sequence data from across the entire protein family of interest, while (multiple) structural alignment algorithms such as MODELLER's salign routine [24, 25], PROMALS3D [29], and Expresso and 3DCoffee [30, 31], can additionally exploit structural data. Ensembler's modular architecture facilitates the implementation of alternative alignment approaches, and we plan to implement some of these in future versions, to allow exploration of the influence of different alignment methods on model quality.

Models are output as PDB-format coordinate files. To minimize file storage requirements, Ensembler uses the Python gzip library to apply compression to all sizeable text files from the modeling stage onwards. The restraints used by MODELLER could potentially be used in alternative additional refinement schemes, and Ensembler thus provides a flag (--write_modeller_restraints_file) for optionally saving these restraints to file. This option is turned off by default, as the restraint files are relatively large (e.g. ∼400 kB per model for protein kinase domain targets), and are not expected to be used by the majority of users.

At this time, the alignment and modeling functions cannot be used to model non-standard amino acids, though we plan to be able to provide this functionality in future versions. Note that the Ensembler functions for target and template selection only include standard amino acids which match the UniProt canonical isoform sequence, and thus any set of targets and templates selected this way should be compatible with the Ensembler alignment and modeling functions.

Filtering of nearly identical models

Because Ensembler treats individual chains from source PDB structures as individual templates, a number of models may be generated with very similar structures if these individual chains are nearly identical in conformation. For this reason, and also to allow users to select for high diversity if they so choose, Ensembler provides a way to filter out models that are very similar in RMSD. The cluster subcommand can thus be used to identify models which differ from other models in terms of RMSD distance by a user-specified cutoff. Clustering is performed using the regular spatial clustering algorithm [9], as implemented in the MSMBuilder Python library [18], which uses mdtraj [32] to calculate RMSD (for C[alpha] atoms only) with a fast quaternion characteristic polynomial (QCP) [33-35] implementation. A minimum distance cutoff (which defaults to 0.6 Å) is used to retain only a single model per cluster.

Refinement of models and filtering of poor models by simulation (E)

A number of refinement methods have been developed to help guide comparative modeling techniques toward more "native-like" and physically consistent conformations [36, 37]. Both short [37] and long [38] molecular dynamics simulations have been employed for this purpose. Here, we utilize short molecular dynamic simulations for two purposes: both to slightly relax the initial comparative models and to eliminate those comparative models that result in highly implausible conformations. This is especially critical here due to the inclusion of even very low sequence identity template structures. We stress that the limited refinement by molecular simulation here is primarily intended as initial relaxation and filtering stages, where implausible models might cause simulations to immediately fail, crash, or generate implausibly high energies or unstable dynamics. Exploration of conformational dynamics to derive MSMs, for example, will inevitably require orders of magnitude more simulation effort--very likely tens of microseconds to milliseconds of aggregate dynamics [8, 10].

Ensembler thus includes a refinement module, which uses short molecular dynamics simulations to refine the models built in the previous step. As well as improving model quality, this also prepares models for subsequent production MD simulation, including solvation with explicit water molecules, if desired.

Models are first subjected to energy minimization (using the L-BFGS algorithm [39], followed by a short molecular dynamics (MD) simulation with an implicit solvent representation. This is implemented using the OpenMM molecular simulation toolkit [2], chosen for its flexible Python API, and high performance GPU-acclerated simulation code. The simulation is run for a default of 100 ps, which in our example applications has been sufficient to filter out poor models (i.e. those with atomic overlaps unresolved by energy minimization, which result in an unstable simulation), as well as helping to relax model conformations. As discussed in the Results section, our example application of the Ensembler pipeline to the human tyrosine kinase family indicated that of the models which failed implicit solvent MD refinement, the vast majority failed within the first 1 ps of simulation.

The simulation protocol and default parameter values have been chosen to represent current "best practices" for the refinement simulations carried out here. As such, the simulation is performed using Langevin dynamics, with a default force field choice of Amber99SB-ILDN [40], along with a modified generalized Born solvent model [41] as implemented in the OpenMM package [2]. Any of the other force fields or implicit water models implemented in OpenMM can be specified using the --ff and --water_model flags respectively. The simulation length can also be controlled via the --simlength flag, and many other important simulation parameters can be controlled from either the API or CLI (via the --api_params flag). The default values are set as follows--timestep: 2 fs; temperature: 300 K; Langevin collision rate: 20 ps-1; pH (used by OpenMM for protonation state assignment): 7. We also draw attention to a recent paper which indicates that lower Langevin collision rates may result in faster phase space exploration [42].

For some studies, it may be useful to specify the protonation states of individual amino acids, rather than rely only on automatic protonation state assignment by OpenMM. The user can do this by listing the residue numbers and their protonation states in a configuration file (manual_overrides.yaml). The necessary formatting for the configuration file is specified in the software documentation, and a template file is written when initializing an Ensembler project. Protonation states are specified by naming the appropriate residue variant type in the force field, e.g. 'ASH' for an aspartic acid residue, as opposed to the aspartate base 'ASP'. Any residues which do not have specific protonation states listed in the configuration file will have protonation states assigned automatically by OpenMM. Note that Ensembler currently only supports residue definitions provided by the forcefield definition files--it does not yet have the ability to derive new forcefield parameters for uncommon amino acids, cofactors, or ions not provided by the forcefield.

Solvation and NPT equilibration (F)

While protein-only models may be sufficient for structural analysis or implicit solvent simulations, Ensembler also provides a stage for solvating models with explicit water and performing a round of explicit-solvent MD refinement/equilibration under isothermal-isobaric (NPT) conditions. The solvation step solvates each model for a given target with the same number of waters to facilitate the integration of data from multiple simulations, which is important for methods such as the construction of MSMs. The target number of waters is selected by first solvating each model with a specified padding distance (default: 10 Å), then taking a percentile value from the distribution (default: 68th percentile). This helps to prevent models with particularly long, extended loops--such as those arising from template structures with unresolved termini--from imposing very large box sizes on the entire set of models. The TIP3P water model [43] is used by default, but any of the other explicit water models available in OpenMM, such as TIP4P-Ew [44], can be specified using the --water_model flag. Models are resolvated with the target number of waters by first solvating with zero padding, then incrementally increasing the box size and resolvating until the target is exceeded, then finally deleting sufficient waters to match the target value. The explicit solvent MD simulation is also implemented using OpenMM, using the Amber99SB-ILDN force field [40] and TIP3P water [43] by default. The force field, water model, and simulation length can again be specified using the --ff, --water_model, and --simlength flags respectively. Further simulation parameters can be controlled via the API or via the CLI --api_params flag. Pressure control is performed with a Monte Carlo barostat as implemented in OpenMM, with a default pressure of 1 atm and a period of 50 timesteps. The remaining simulation parameters have default values set to the same as for the implicit solvent MD refinement.

Model validation with MolProbity score

Ensembler provides a function for validating model quality and filtering models using MolProbity [45, 46]--a widely used tool for validation of protein models, which provides a numerical score derived from features such as steric clashes between atoms, bond geometry, Ramachandran angles, sidechain rotamer outliers, backbone deviations, and the presence of cis-peptides. This function is accessed via the validate subcommand, which for a given target will output a text file containing a list of model IDs sorted by validation score. The optional --modeling_stage flag specifies which of the three main Ensembler modeling stages to validate--the initial comparative modeling stage, the implicit MD refinement stage, or the explicit MD refinement stage. If this flag is not used, Ensembler defaults to selecting the latest stage for which models have been generated. The output text file can be used to filter models based on validation score, for example by using the package_models subcommand. Protein model validation is a challenging problem and an active area of research for many groups, including the developers of Ensembler. We plan to implement further validation methods in future versions of Ensembler.

Packaging

Ensembler provides a packaging module which can be used to prepare models for subsequent downstream use, such as the use of distributed or cluster computing resources for the generation of MSMs [8-10]. The package_models subcommand currently provides functions (specified via the --package_for flag) for compressing models in preparation for data transfer, or for organizing them with the appropriate directory and file structure for production simulation on the distributed computing platform Folding@home [4]. The module could easily be extended to add methods for preparing models for other purposes. For example, production simulations could alternatively be run using Copernicus [5, 6]--a framework for performing parallel adaptive MD simulations-- or GPUGrid [7]--a distributing computing platform which relies on computational power voluntarily donated by the owners of nondedicated GPU-equipped computers.

An important use of the packaging stage is to filter models based on model quality. At the current time, the available filtering options are based on either target-template sequence identity or MolProbity validation score. The package_models subcommand includes optional flags for specifying a sequence identity cutoff (so that only models with a target-template sequence identity above the specified percentage are chosen), a MolProbity validation score cutoff (to choose only models with lower validation scores, which indicate better model quality), or a MolProbity validation score percentile (to choose only models with validation scores lower than the value at the given percentile).

Models can also be exported into trajectory files for the purpose of performing structural analyses across model ensembles using tools like MDTraj [32]. This is done using the mktraj subcommand, which writes model coordinates for a given target to a Gromacs [47, 48] XTC format trajectory (chosen for its wide usage and data compression). Each frame in the trajectory represents a single model, and models are sorted in descending order of target-template sequence identity. Also output for each target are a PDB coordinate file (for use as a topology input file) and a CSV file containing model IDs (in the same order as the frames in the trajectory file) and other data such as target-template sequence identity. Using the --modeling_stage flag, models can be selected from any of three Ensembler modeling stages--after the initial comparative modeling stage, after implicit MD refinement, or after explicit MD refinement. If this flag is not used, Ensembler defaults to selecting the latest stage for which models have been generated.

We stress that, despite evidence suggesting that there is a correspondence between solution-state dynamics and structural diversity of related template proteins [16], all models--especially those derived from low sequence identity templates--are not necessarily representative of conformations thermally accessible to the template proteins of interest. Care must be exercised in the use and analysis of these models.

Other features

Tracking provenance information.

To aid the user in tracking the provenance of each model, each pipeline function also outputs a metadata file, which helps to link data to the software version used to generate it (both Ensembler and its dependencies), and also provides timing and performance information, and other data such as hostname.

Rapidly modeling a single template.

For users interested in simply using Ensembler to rapidly generate a set of models for a single template sequence, Ensembler provides a command-line tool quickmodel, which performs the entire pipeline for a single target with a small number of templates. For larger numbers of models (such as entire protein families), modeling time is greatly reduced by using the main modeling pipeline, which is parallelized via MPI, distributing computation across each model (or across each template, in the case of the loop reconstruction code), and scaling (in a "pleasantly parallel" manner) up to the number of models generated.

Results

Modeling of all human tyrosine kinase catalytic domains

As a first application of Ensembler, we have built models for the human TK family. TKs (and protein kinases in general) play important roles in many cellular processes and are involved in a number of types of cancer [49]. For example, a translocation between the TK Abl1 and the pseudokinase Bcr is closely associated with chronic myelogenous leukemia [50], while mutations of Src are associated with colon, breast, prostate, lung, and pancreatic cancers [51]. Protein kinase domains are thought to have multiple accessible metastable conformation states, and much effort is directed at developing kinase inhibitor drugs which bind to and stabilize inactive conformations [52]. Kinases are thus a particularly interesting subject for study with MSM methods [53], and this approach stands to benefit greatly from the ability to exploit the full body of available genomic and structural data within the kinase family, e.g. by generating large numbers of starting configurations to be used in highly parallel MD simulation.

We selected all human TK domains annotated in UniProt as targets, and all available structures of protein kinase domains (of any species) as templates, using the commands shown in Box 1. This returned 93 target sequences and 4433 template structures, giving a total of 412269 target-template pairs. The templates were derived from 3028 individual PDB entries and encompassed 23 different species, with 3634 template structures from human kinase constructs.

The resultant models are available as part of a supplementary dataset which can be

Fig 2. Number of PDB structures available for each TK target.

Data is shown for each of the 93 TK kinase domains, sorted in order of the number of available PDB structures for each domain. The labels indicate the UniProt name for the target protein plus an index for the kinase domain (three of the selected proteins have two kinase domains). Each PDB chain is counted individually, and only chains which contain the target domain are counted.

http://dx.doi.org/10.1371/journal.pcbi.1004728.g002

Ensembler modeling statistics

Crystallographic structures of kinase catalytic domains generally contain a significant number of missing residues (median 11, mean 14, standard deviation 13, max 102) due to the high mobility of several loops (Fig 3, top), with a number of these missing spans being significant in length (median 5, mean 7, standard deviation 6, max 82; Fig 3, bottom). To reduce the reliance on the MODELLER rapid model construction stage to reconstruct very long unresolved loops, unresolved template residues were first remodeled using the loopmodel subcommand. Out of 3666 templates with one or more missing residues, 3134 were successfully remodeled by the Rosetta loop modeling stage (with success defined simply as program termination without error); most remodeling failures were attributable to unsatisfiable spatial constraints imposed by the original template structure. There was some correlation between remodeling failures and the number of missing residues (Fig 3, top); templates for which remodeling failed had a median of 20 missing residues, compared to a median of 14 missing residues for templates for which remodeling was successful (when excluding templates with no missing residues).

Fig 3. Distributions for the number of missing residues in the TK templates.

Upper: The number of missing residues per template, for all templates (blue) and for only those templates for which template remodeling with the loopmodel subcommand failed (red). Templates for which remodeling failed had a median of 20 missing residues, compared to a median of 14 missing residues for templates for which remodeling was successful. Lower: The number of residues in each missing loop, for all templates.

http://dx.doi.org/10.1371/journal.pcbi.1004728.g003

Following loop remodeling, the Ensembler pipeline was performed up to and including the implicit solvent MD refinement stage, which completed with 389067 (94%) surviving models across all TKs. To obtain statistics for the solvation stage without generating a sizeable amount of coordinate data (with solvated PDB coordinate files taking up about 0.9 MB each), the solvate subcommand was performed for two representative individual kinases (Src and Abl1).

The number of models which survived each stage are shown in Fig 1, indicating that the greatest attrition occurred during the modeling stage. The number of refined models for each target ranged from 4046 to 4289, with a median of 4185, mean of 4184, and standard deviation of 57. Fig 1 also indicates the typical timing achieved on a cluster for each stage, showing that the build_models and refine_implicit_md stages are by far the most compute-intensive.

The files generated for each model (up to and including the implicit solvent MD refinement stage) totaled ∼116 kB in size, totalling 0.5 GB per TK target or 42 GB for all 93 targets. The data generated per model breaks down as 39 kB for the output from the modeling stage (without saving MODELLER restraints files, which are about 397 kB per model) and 77 kB for the implicit solvent MD refinement stage.

Evaluation of model quality and utility

All tyrosine kinases.

To evaluate the variety of template sequence similarities relative to each target sequence, we calculated sequence identity distributions, as shown in Fig 4. This suggests an intuitive division into three categories, with 355712 models in the 0-35% sequence identity range, 51330 models in the 35-55% range, and 5227 models in the 55-100% range. We then computed the RMSD distributions for the models created for each target (relative to the model derived from the template with highest sequence identity) Fig 5, to assess the diversity of conformations captured by the modeling pipeline. Furthermore, to understand the influence of sequence identity on the conformational similarities of the resulting models, the RMSD distributions were stratified based on the three sequence identity categories described above. This analysis indicates that higher sequence identity templates result in models with lower RMSDs, while templates with remote sequence identities result in larger RMSDs on average, recapitulating the observation made years ago by Chothia and Lesk [54].

Fig 4. Template-target sequence identity distribution for human tyrosine kinase catalytic domains.

Sequence identities are calculated from all pairwise target-template alignments, where targets are human kinase catalytic domain sequences and templates are all kinase catalytic domains from any organism with structures in the PDB, as described in the text. A kernel density estimate of the target-template sequence identity probability density function is shown as a solid line with shaded region, while the corresponding cumulative distribution function is shown as a dashed line.

http://dx.doi.org/10.1371/journal.pcbi.1004728.g004

Fig 5. Distribution of RMSDs to all TK catalytic domain models relative to the model derived from the highest sequence identity template.

Distributions are built from data from all 93 TK domain targets. To better illustrate how conformational similarity depends on sequence identity, the lower plot illustrates the distributions as stratified into three sequence identity classes: high identity (55-100%), moderate identity (35-55%), and remote identity (0-35%). The plotted distributions have been smoothed using kernel density estimation.

http://dx.doi.org/10.1371/journal.pcbi.1004728.g005

We also used the ensembler validate subcommand to subject the refined models to analysis with MolProbity [45, 46]. The MolProbity scores varied from 0.92 to 4.80, with a median of 3.84, a mean of 3.22, and a standard deviation of 1.07. Lower numbers represent better quality models. When stratified by the same sequence identity ranges as above, the mean scores were as follows: 2.96 (55-100% sequence identity), 3.13 (35-55% sequence identity), 3.24 (0-35% sequence identity). This indicates that models with lower target-template sequence identities tend to be lower quality according to MolProbity analysis, as would be expected.

We also analyzed the potential energies of the models at the end of the implicit solvent MD refinement stage. These ranged from -14180 kT to -3160 kT, with a median of -9501 kT, mean of -9418 kT, and a standard deviation of 1198 kT (with a simulation temperature of 300 K). The distributions--stratified using the same sequence identity ranges again--are plotted in Fig 6, indicating that higher sequence identity templates tend to result in slightly lower energy models. Of the 4973 models which failed to complete the implicit refinement MD stage, all except 9 failed within the first 1 ps of simulation.

Fig 6. Distribution of final energies from implicit solvent MD refinement of TK catalytic domain models.

To illustrate how the energies are affected by sequence identity, the models are separated into three sequence identity classes: high identity (55-100%), moderate identity (35-55%), and remote identity (0-35%). The plotted distributions have been smoothed using kernel density estimation. Refinement simulations were carried out at the default temperature of 300 K.

http://dx.doi.org/10.1371/journal.pcbi.1004728.g006

Src and Abl1.

To provide a more detailed evaluation of the variety and utility of generated models, we have analyzed two specific TKs (Src and Abl1) in depth. Due to their importance in cancer, these kinases have been the subject of numerous detailed structural and simulation studies. In terms of structural data, a large number of crystal structures have been solved (with or without ligands such as nucleotide substrate mimetics or small-molecule inhibitors), revealing a variety of conformations accessible to these kinases. A recent large-scale MSM study has also studied the activation pathway of Src [53], while a separate study employed biased sampling techniques to dissect the role of conformational changes in selectivity and affinity of imatinib recognition of Abl [55].

Visualizing model structural diversity.

Fig 7 shows a superposition of a set of representative models of Src and Abl1. Models were first stratified into three ranges, based on the structure of the sequence identity distribution (Fig 4), then subjected to RMSD-based k-medoids clustering (using the msmbuilder clustering package [18]) to pick three representative models from each sequence identity range. Each model is colored and given a transparency based on the sequence identity between the target and template sequence. The figure gives an idea of the variance present in the generated models. High sequence identity models (in opaque blue) tend to be quite structurally similar, with some variation in loops or changes in domain orientation.

Fig 7. Superposition of clustered models of Src and Abl1.

Superposed renderings of nine models each for Src and Abl1, giving some indication the diversity of conformations generated by Ensembler. The models for each target were divided into three sequence identity ranges (as in Fig 5), and RMSD-based k-medoids clustering was performed (using the msmbuilder clustering package [18]) to select three clusters from each. The models shown are the centroids of each cluster. Models are colored and given transparency based on their sequence identity, so that high sequence identity models are blue and opaque, while lower sequence identity models are transparent and red.

http://dx.doi.org/10.1371/journal.pcbi.1004728.g007

The Abl1 renderings in Fig 7 indicate one high sequence identity model with a long unstructured region at one of the termini, which was unresolved in the original template structure. While such models are not necessarily incorrect or undesirable, it is important to be aware of the effects they may have on production simulations performed under periodic boundary conditions, as long unstructured termini can be prone to interact with a protein's periodic image. Lower sequence identity models (in transparent white or red) indicate much greater variation in all parts of the structure. We believe the mix of high and low sequence identity models to be particularly useful for methods such as MSM building, which require thorough sampling of the conformational landscape. The high sequence identity models could be considered to be the most likely to accurately represent true metastable states. Conversely, the lower sequence identity models could be expected to help push a simulation into regions of conformation space which might take intractably long to reach if starting a single metastable conformation.

Comparison with known biochemically relevant conformations.

To evaluate the models of Src and Abl1 in the context of the published structural biology literature on functionally relevant conformations, we have focused on two residue pair distances thought to be important order parameters for the regulation of protein kinase domain activity. We use the residue numbering schemes for chicken Src (commonly employed in the literature even in reference to human Src) [56, 57] and human Abl1 isoform A [58-60] respectively; see S1 Text for the exact sequences and numbering schemes used.

Fig 8 shows two structures of Src believed to represent inactive (PDB code: 2SRC) [56] and active (PDB code: 1Y57) [57] states. One notable feature which distinguishes the two structures is the transfer of an electrostatic interaction of E310 from R409 (in the inactive state) to K295 (in the active state), brought about by a rotation of the [alpha] C-helix. These three residues are also well conserved [61], and a number of experimental and simulation studies have suggested that this electrostatic switching process plays a role in a regulatory mechanism shared across the protein kinase family [53, 62, 63]. As such, we have plotted the distance between these two residue pairs for the Ensembler models for Src and Abl1, as well as Flt4 (Fig 9). The models all show strong coverage of regions in which either of the electrostatic interactions is fully formed (for models across all levels of target-template sequence identity), as well as a wide range of regions in-between (mainly models with low sequence identity). We thus expect that such a set of models, if used as starting configurations for highly parallel MD simulation, could greatly aid in sampling of functionally relevant conformational states. The Flt4 models (Fig 8c) are of particular note, as there are no available crystal structures of the kinase domain of this TK protein (which is involved in tumor angiogenesis and lymphangiogenesis [64]), yet the models generated here include structural motifs which are conserved and of known importance to other proteins of the same family.

Fig 8. Two structures of Src, indicating certain residues involved in activation.

In the inactive state, E310 forms a salt bridge with R409. During activation, the [alpha] C-helix (green) moves and rotates, orienting E310 towards the ATP-binding site and allowing it to instead form a salt bridge with K295. This positions K295 in the appropriate position for catalysis. Note that ANP (phosphoaminophosphonic acid-adenylate ester; an analog of ATP) is only physically present in the 2SRC structure. To aid visualization of the active site in 1Y57, it has been included in the rendering by structurally aligning the surrounding homologous protein residues.

http://dx.doi.org/10.1371/journal.pcbi.1004728.g008

Fig 9. Src, Abl1, and Flt4 models projected onto the distances between two conserved residue pairs, colored by sequence identity.

Two Src structures (PDB entries 1Y57 [57] and 2SRC [56]) are projected onto the plots for reference, representing active and inactive states respectively. These structures and the residue pairs analyzed here are depicted in Fig 8. Distances are measured between the center of masses of the three terminal sidechain heavy atoms of each residue. The atom names for these atoms, according to the PDB coordinate files for both reference structures, are--Lys: NZ, CD, CE (ethylamine); Glu: OE1, CD, OE2 (carboxylate); Arg: NH1, CZ, NH2 (part of guanidine).

http://dx.doi.org/10.1371/journal.pcbi.1004728.g009

Availability and Future Directions

Availability

The code for Ensembler is hosted on the collaborative open source software development platform GitHub (github.com/choderalab/ensembler). The latest release can be installed via the conda package manager for Python (conda.pydata.org), using the commands shown in Box 2. This will install all required dependencies, though a license must first be obtained for MODELLER. The optional dependencies Rosetta and MolProbity are not available through the conda package manager, and thus must be installed separately by the user according to the instructions for those packages. The latest source can be downloaded from the GitHub repository, which also contains up-to-date instructions for building and installing the code. Documentation can be found at ensembler.readthedocs.org.

Box 2. Ensembler installation using conda in a bash shell. A license for MODELLER (free for academic use) must first be obtained, and the license key should be used in place of the text "Modeller license key"

export KEY_MODELLER = "Modeller license key"

conda config -add channels https://conda.anaconda.org/omnia

conda config -add channels https://conda.anaconda.org/salilab

conda install ensembler

A supplementary dataset can also be downloaded from the Dryad Digital Repository (DOI: 10.5061/dryad.7fg32). This contains the TK models described in the section, general information on the targets and templates, plus a script and instructions for regenerating the same dataset.

Future directions

We recognize that the current version of Ensembler has a number of limitations that bound its domain of applicability: support for nonnatural amino acids is currently rudimentary and confined to those already appearing in the forcefield; cofactors cannot currently be automatically modeled in; ligands, cofactors, and nonnatural amino acids cannot yet be automatically parameterized; protonation state assignment is limited to selection of the most populated state based on the intrinsic pKa or user-specified overrides; the modeling of missing loops is rudimentary, relying on the subsequent dynamics for relaxation; there is not yet support for modeling of distinct domains from different templates, or the use of multiple templates to model a single domain. Nevertheless, there are a great number of use cases for this first version of an automated tool for simulation preparation at the superfamily scale. To expand this domain of applicability, there are a number of obvious additions and improvements which we plan to implement in future versions of Ensembler.

Template remodeling.

The lack of crystallographically-resolved regions of template structures presents a challenge to deriving structures from these templates by comparative modeling, especially in kinases, where loops are frequently unresolved. Improvements over the Rosetta-based strategy described here are likely possible, especially given the number of modeling failures observed in the template refinement stage (Fig 3). An alternative approach could be to re-refine complete-chain template structures to the experimentally-derived electron density or scattering data deposited in the RCSB using methods capable of exploiting the scattering data and crystallographic symmetry [65]. Even if definitive placement of these unresolved regions is impossible, plausible locations constrained by weak scattering data and strong steric exclusion of crystallographic neighbors may provide a great deal of useful information, especially when combined with forcefield priors [66].

Comparative modeling.

Comparative protein modeling can be approached in a number of different ways, with varying degrees of complexity. The comparative modeling stage of Ensembler currently uses MODELLER, but a number of excellent alternatives--such as RosettaCM [13] and the I-TASSER Suite [14]--can be added as user-selectable alternative choices. Additional options could be added to allow more expensive loop-modeling approaches to be employed to handle long insertions.

Protonation states.

Some amino acids can exist in different protonation states, depending on pH and on their local environment. These protonation states can have important effects on biological processes. For example, long timescale MD simulations have suggested that the conformation of the DFG motif of the TK Abl1--believed to be an important regulatory mechanism [67]--is controlled by protonation of the aspartate [68]. Currently, protonation states are assigned simply based on pH (a user-controllable parameter). At neutral pH, histidines have two protonation states which are approximately equally likely, and in this situation the selection is therefore made based on which state results in a better hydrogen bond. It would be highly desirable to instead use a method which assigns amino acid protonation states based on a rigorous assessment of the local environment. We thus plan to implement an interface and command-line function for assigning protonation states with MCCE2 [69-71], which uses electrostatics calculations combined with Monte Carlo sampling of side chain conformers to calculate pKa values.

Cofactors, structural ions, and ligands.

Many proteins require the presence of various types of non-protein atoms and molecules for proper function, such as metal ions (e.g. Mg+2), cofactors (e.g. ATP) or post-translational modifications (e.g. phosphorylation, methylation, glycosylation, etc.), and we thus plan for Ensembler to eventually have the capability to include such entities in the generated models. Binding sites for metal ions are frequently found in proteins, often playing a role in catalysis. For example, protein kinase domains contain two binding sites for divalent metal cations, and display significantly increased activity in the presence of Mg2+ [72], the divalent cation with highest concentration in mammalian cells. Metal ions are often not resolved in experimental structures of proteins, but by taking into account the full range of available structural data, it should be possible in many cases to include metal ions based on the structures of homologous proteins. We are careful to point out, however, that metal ion parameters in classical MD force fields have significant limitations, particularly in their interactions with proteins [73]. Cofactors and post-translational modifications are also often not fully resolved in experimental structures, and endogenous cofactors are frequently substituted with other molecules to facilitate experimental structural analysis. Future extensions to Ensembler could transfer cofactor and ion coordinates from homologous proteins in which these components are resolved.

Post-translationally modified amino acids and other molecules without forcefield parameters.

A major challenge in the preparation of simulations of proteins of interest is the wide variety of post-translational modifications possible that are often functionally or structurally relevant. Often, forcefields lack parameters for these residues, or for other cofactors or ligands that might be vital to probing the relevant structural dynamics of these systems. While tools such as Antechamber [74, 75] can rapidly generate small molecule parameters in an automated manner, the parameterization of polymeric residues or covalently attached cofactors is much more challenging. In addition, small molecule forcefields are generally tied to specific corresponding protein and nucleic acid forcefields, meaning that different procedures may be needed to generate consistent parameters.

Long insertions and deletions.

Another limitation with the present version of Ensembler involves the treatment of members of a protein family with especially long residue insertions or deletions. For example, the set of all human protein kinase domains listed in UniProt have a median length of 265 residues (mean 277) and a standard deviation of 45, yet the minimum and maximum lengths are 102 and 801 respectively. The latter value corresponds to the protein kinase domain of serine/threonine-kinase greatwall, which includes a long insertion between the two main lobes of the catalytic domain. In principle, such insertions could be excluded from the generated models, though a number of questions would arise as to how best to approach this.

Markov state model (MSM) construction and model utility.

We are actively utilizing Ensembler-generated models to seed the construction of Markov state models (MSMs) [8, 10]. While the observation that high sequence identity templates are likely to reflect accessible solution-phase conformations suggests that a number of these models occupy thermally accessible regions of configuration space [16], many models--especially those derived from very low sequence identity templates--are likely to be highly unrepresentative of conformations populated at equilibrium by the target protein. It is likely that even with hundreds of microseconds to milliseconds of aggregated dynamics, many of these poor quality models will remain trapped in inaccessible and irrelevant regions of configuration space. Standard approaches to MSM construction now employ an ergodic trimming step [18, 19] to prune away disconnected minor regions of configuration space, and this step is expected to be essential in the successful construction of MSMs using Ensembler-derived models.

Conclusion

We believe Ensembler to be an important first step toward enabling computational modeling and simulation of proteins on the scale of entire protein families, and suggest that it could likely prove useful for tasks beyond its original aim of providing diverse starting configurations for MD simulations. The code is open source and has been developed with extensibility in mind, in order to facilitate its customization for a wide range of potential uses by the wider scientific community.

Supporting Information

S1 Text. Sequences and residue numbering schemes for Src and Abl1.

(DOCX)

Acknowledgments

The authors are grateful to Robert McGibbon (Stanford) and Arien S. Rustenburg (MSKCC) for many excellent software engineering suggestions. The authors thank Nicholas M. Levinson (University of Minnesota), Markus A. Seeliger (Stony Brook), Diwakar Shukla (Stanford), and Avner Schlessinger (Mount Sinai) for helpful scientific feedback on modeling kinases. The authors are grateful to Benjamin Webb and Andrej Sali (UCSF) for help with the MODELLER package, Peter Eastman and Vijay Pande (Stanford) for assistance with OpenMM, and Marilyn Gunner (CCNY) for assistance with MCCE2, as well as the anonymous referees for constructive feedback on this manuscript.

Author Contributions

Conceived and designed the experiments: DLP JDC. Performed the experiments: DLP JDC. Analyzed the data: DLP. Contributed reagents/materials/analysis tools: DLP PBG SMH KAB JDC. Wrote the paper: DLP PBG SMH KAB JDC.

References

Lee GM, Craik CS. Trapping moving targets with small molecules. Science. 2009;324:213. doi: 10.1126/science.1169378. pmid:19359579

Eastman P, Friedrichs MS, Chodera JD, Radmer RJ, Bruns CM, Ku JP, et al. OpenMM 4: A Reusable, Extensible, Hardware Independent Library for High Performance Molecular Simulation. J Chem Theory Comput. 2012;9(1):461-469. doi: 10.1021/ct300857j. pmid:23316124

Salomon-Ferrer R, Götz AW, Poole D, Grand SL, Walker RC. Routine microsecond molecular dynamics simulations with AMBER on GPUs. 2. Explicit solvent particle mesh Ewald. J Chem Theor Comput. 2013;9(9):3878-3888. doi: 10.1021/ct400314y

Shirts M, Pande VS. Screen Savers of the World Unite! Science. 2000 Dec;290(5498):1903-1904. pmid:17742054 doi: 10.1126/science.290.5498.1903

Pronk S, Larsson P, Pouya I, Bowman GR, Haque IS, Beauchamp K, et al. Copernicus: A New Paradigm for Parallel Adaptive Molecular Dynamics. In: Proceedings of 2011 International Conference for High Performance Computing, Networking, Storage and Analysis. SC'11. New York, NY, USA: ACM; 2011. p. 60:1-60:10

Pronk S, Pouya I, Lundborg M, Rotskoff G, Wesén B, Kasson PM, et al. Molecular Simulation Workflows as Parallel Algorithms: The Execution Engine of Copernicus, a Distributed High-Performance Computing Platform. Journal of Chemical Theory and Computation. 2015 Apr;. doi: 10.1021/acs.jctc.5b00234. pmid:26575558

Buch I, Harvey MJ, Giorgino T, Anderson DP, De Fabritiis G. High-Throughput All-Atom Molecular Dynamics Simulations Using Distributed Computing. Journal of Chemical Information and Modeling. 2010 Mar;50(3):397-403. doi: 10.1021/ci900455r. pmid:20199097

Pande VS, Beauchamp K, Bowman GR. Everything you wanted to know about Markov State Models but were afraid to ask. Methods. 2010;52:99-105. doi: 10.1016/j.ymeth.2010.06.002. pmid:20570730

Prinz JH, Wu H, Sarich M, Keller B, Fischbach M, Held M, et al. Markov models of molecular kinetics: Generation and validation. J Chem Phys. 2011;134:174105. doi: 10.1063/1.3565032. pmid:21548671

Chodera JD, Noé F. Markov state models of biomolecular conformational dynamics. Curr Opin Struct Biol. 2014;25:135-144. doi: 10.1016/j.sbi.2014.04.002. pmid:24836551

Moult J, Fidelis K, Kryshtafovych A, Schwede T, Tramontano A. Critical assessment of methods of protein structure prediction (CASP) x2014; round x. Proteins: Structure, Function, and Bioinformatics. 2014 Feb;82:1-6. doi: 10.1002/prot.24452

Baker D, Sali A. Protein Structure Prediction and Structural Genomics. Science. 2001 Oct;294(5540):93-96. doi: 10.1126/science.1065659. pmid:11588250

Song Y, DiMaio F, Wang RYR, Kim D, Miles C, Brunette T, et al. High-resolution comparative modeling with RosettaCM. Structure. 2013;21:1735-1742. doi: 10.1016/j.str.2013.08.005. pmid:24035711

Yang J, Yan R, Roy A, Xu D, Poisson J, Zhang Y. The I-TASSER Suite: Protein structure and function prediction. Nature Methods. 2015;12:7. doi: 10.1038/nmeth.3213. pmid:25549265

van der Kamp MW, Schaeffer RD, Jonsson AL, Scouras AD, Simms AM, Toofanny RD, et al. Dynameomics: A comprehensive database of protein dynamics. Structure. 2010;18:423-435. doi: 10.1016/j.str.2010.01.012. pmid:20399180

Friedland GD, Lakomek NA, Griesinger C, Meiler J, Kortemme T. A correspondence between solution-state dynamics of an individual protein and the sequence and conformational diversity of its family. PLoS Comput Biol. 2009;5:e1000393. doi: 10.1371/journal.pcbi.1000393. pmid:19478996

Weinkam P, Pons J, Sali A. Structure-based model of allostery predicts coupling between distant sites. Proceedings of the National Academy of Sciences of the United States of America. 2012 Mar;109(13):4875-4880. doi: 10.1073/pnas.1116274109. pmid:22403063

Beauchamp KA, Bowman GR, Lane TJ, Maibaum L, Haque IS, Pande VS. MSMBuilder2: Modeling Conformational Dynamics on the Picosecond to Millisecond Scale. Journal of Chemical Theory and Computation. 2011 Oct;7(10):3412-3419. doi: 10.1021/ct200463m. pmid:22125474

Scalco R, Caflisch A. Equilibrium distribution from distributed computing (simulations of protein folding). The Journal of Physical Chemistry B. 2011 May;115(19):6358-6365. doi: 10.1021/jp2014918. pmid:21517045

Consortium TU. UniProt: a hub for protein information. Nucleic Acids Research. 2015 Jan;43(D1):D204-D212. doi: 10.1093/nar/gku989

Velankar S, Dana JM, Jacobsen J, van Ginkel G, Gane PJ, Luo J, et al. SIFTS: Structure Integration with Function, Taxonomy and Sequences resource. Nucleic Acids Research. 2013 Jan;41(D1):D483-D489. doi: 10.1093/nar/gks1258. pmid:23203869

Qian B, Raman S, Das R, Bradley P, McCoy AJ, Read RJ, et al. High-resolution structure prediction and the crystallographic phase problem. Nature. 2007 Nov;450(7167):259-264. doi: 10.1038/nature06249. pmid:17934447

Wang C, Bradley P, Baker D. Protein-Protein Docking with Backbone Flexibility. Journal of Molecular Biology. 2007 Oct;373(2):503-519. pmid:17825317 doi: 10.1016/j.jmb.2007.07.050

Fiser A, Do RKG, Sali A. Modeling of loops in protein structures. Protein Science. 2000 Jan;9(9):1753-1773. doi: 10.1110/ps.9.9.1753. pmid:11045621

Sali A, Blundell TL. Comparative Protein Modelling by Satisfaction of Spatial Restraints. Journal of Molecular Biology. 1993 Dec;234(3):779-815. doi: 10.1006/jmbi.1993.1626. pmid:8254673

Cock PJA, Antao T, Chang JT, Chapman BA, Cox CJ, Dalke A, et al. Biopython: freely available Python tools for computational molecular biology and bioinformatics. Bioinformatics (Oxford, England). 2009 Jun;25(11):1422-1423. doi: 10.1093/bioinformatics/btp163

Gonnet GH, Cohen MA, Brenner SA. Exhaustive matching of the entire protein sequence database. Science. 1992;256:1443-1445. doi: 10.1126/science.1604319. pmid:1604319

Thompson JD, Linard B, Lecompte O, Poch O. A Comprehensive Benchmark Study of Multiple Sequence Alignment Methods: Current Challenges and Future Perspectives. PLoS ONE. 2011 Mar;6(3):e18093. doi: 10.1371/journal.pone.0018093. pmid:21483869

Pei J, Kim BH, Grishin NV. PROMALS3D: a tool for multiple protein sequence and structure alignments. Nucleic Acids Research. 2008 Apr;36(7):2295-2300. doi: 10.1093/nar/gkn072. pmid:18287115

Armougom F, Moretti S, Poirot O, Audic S, Dumas P, Schaeli B, et al. Expresso: automatic incorporation of structural information in multiple sequence alignments using 3D-Coffee. Nucleic Acids Research. 2006 Jul;34(Web Server issue):W604-608. doi: 10.1093/nar/gkl092. pmid:16845081

Poirot O, Suhre K, Abergel C, O'Toole E, Notredame C. 3DCoffee@igs: a web server for combining sequences and structures into a multiple sequence alignment. Nucleic Acids Research. 2004 Jul;32(Web Server issue):W37-40. doi: 10.1093/nar/gkh382. pmid:15215345

McGibbon RT, Beauchamp KA, Schwantes CR, Wang LP, Hernández CX, Harrigan MP, et al. MDTraj: a modern, open library for the analysis of molecular dynamics trajectories. bioRxiv. 2014;. doi: 10.1101/008896

Theobald DL. Rapid calculation of RMSDs using a quaternion-based characteristic polynomial. Acta Cryst A. 2005;61:478-480. doi: 10.1107/S0108767305015266

Liu P, Agrafiotis DK, Theobald DL. Fast determination of the optimal rotation matrix for macromolecular superpositions. J Comput Chem. 2010;31:1561-1563. doi: 10.1002/jcc.21439. pmid:20017124

Liu P, Agrafiotis DK, Theobald DL. Reply to comment on "Fast determination of the optimal rotation matrix for macromolecular superpositions". J Comput Chem. 2011;32:185-186. doi: 10.1002/jcc.21606

MacCallum JL, Pérez A, Schnieders MJ, Hua L, Jacobson MP, Dill KA. Assessment of protein structure refinement in CASP9. Proteins: Structure, Function, and Bioinformatics. 2011 Jan;79(S10):74-90. doi: 10.1002/prot.23131

Zhang Y. Protein structure prediction: when is it useful? Current Opinion in Structural Biology. 2009 Apr;19(2):145-155. doi: 10.1016/j.sbi.2009.02.005. pmid:19327982

Raval A, Piana S, Eastwood MP, Dror RO, Shaw DE. Refinement of protein structure homology models via long, all-atom molecular dynamics simulations. Proteins: Structure, Function, and Bioinformatics. 2012 Aug;80(8):2071-2079. doi: 10.1002/prot.24098

Liu DC, Nocedal J. On the Limited Memory BFGS Method For Large Scale Optimization. Mathematical Programming. 1989;45:503-528. doi: 10.1007/BF01589116

Lindorff-Larsen K, anad Kim Palmo SP, Maragakis P, Klepeis JL, Dror RO, Shaw DE. Improved side-chain torsion potentials for the Amber ff99SB protein force field. Proteins. 2010;78:1950-1958. doi: 10.1002/prot.22711. pmid:20408171

Onufriev A, Bashford D, Case DA. Exploring protein native states and large-scale conformational changes with a modified generalized Born model. Proteins. 2004;55(22):383-394. doi: 10.1002/prot.20033. pmid:15048829

Basconi JE, Shirts MR. Effects of Temperature Control Algorithms on Transport Properties and Kinetics in Molecular Dynamics Simulations. Journal of Chemical Theory and Computation. 2013 Jul;9(7):2887-2899. doi: 10.1021/ct400109a. pmid:26583973

Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML. Comparison of simple potential functions for simulating liquid water. Journal of Chemical Physics. 1983;79:926-935. doi: 10.1063/1.445869

Horn HW, Swope WC, Pitera JW, Madura JD, Dick TJ, Hura GL, et al. Development of an improved four-site water model for biomolecular simulations: TIP4P-Ew. The Journal of Chemical Physics. 2004 May;120(20):9665-9678. doi: 10.1063/1.1683075. pmid:15267980

Davis IW, Leaver-Fay A, Chen VB, Block JN, Kapral GJ, Wang X, et al. MolProbity: all-atom contacts and structure validation for proteins and nucleic acids. Nucleic Acids Research. 2007 Jul;35(suppl 2):W375-W383. doi: 10.1093/nar/gkm216. pmid:17452350

Chen VB, Arendall WB, Headd JJ, Keedy DA, Immormino RM, Kapral GJ, et al. MolProbity: all-atom structure validation for macromolecular crystallography. Acta Crystallographica Section D, Biological Crystallography. 2010 Jan;66(Pt 1):12-21. doi: 10.1107/S0907444909042073. pmid:20057044

Berendsen HJC, van der Spoel D, van Drunen R. GROMACS: A message-passing parallel molecular dynamics implementation. Comp Phys Comm. 1995 Sep;91(1-3):43-56. doi: 10.1016/0010-4655(95)00042-E

Lindahl E, Hess B, van der Spoel D. GROMACS 3.0: a package for molecular simulation and trajectory analysis. J Mol Model. 2001;7:306-317

Krause DS, Van Etten RA. Tyrosine Kinases as Targets for Cancer Therapy. New England Journal of Medicine. 2005 Jul;353(2):172-187. doi: 10.1056/NEJMra044389. pmid:16014887

Greuber EK, Smith-Pearson P, Wang J, Pendergast AM. Role of ABL family kinases in cancer: from leukaemia to solid tumours. Nature Reviews Cancer. 2013 Aug;13(8):559-571. doi: 10.1038/nrc3563. pmid:23842646

Kim LC, Song L, Haura EB. Src kinases as therapeutic targets for cancer. Nature Reviews Clinical Oncology. 2009 Oct;6(10):587-595. doi: 10.1038/nrclinonc.2009.129. pmid:19787002

Liu Y, Gray NS. Rational design of inhibitors that bind to inactive kinase conformations. Nature Chemical Biology. 2006 Jul;2(7):358-364. doi: 10.1038/nchembio799. pmid:16783341

Shukla D, Meng Y, Roux B, Pande VS. Activation pathway of Src kinase reveals intermediate states as targets for drug design. Nature Commun. 2014;5:3397. doi: 10.1038/ncomms4397

Chothia C, Lesk AM. The relation between the divergence of sequence and structure in proteins. EMBO J. 1986;5:823-826. pmid:3709526

Lin YL, Meng Y, Jiang W, Roux B. Explaining why Gleevec is a specific and potent inhibitor of Abl kinase. Proc Natl Acad Sci USA. 2013;110:1664. doi: 10.1073/pnas.1214330110. pmid:23319661

Xu W, Doshi A, Lei M, Eck MJ, Harrison SC. Crystal Structures of c-Src Reveal Features of Its Autoinhibitory Mechanism. Molecular Cell. 1999 May;3(5):629-638. doi: 10.1016/S1097-2765(00)80356-1. pmid:10360179

Cowan-Jacob SW, Fendrich G, Manley PW, Jahnke W, Fabbro D, Liebetanz J, et al. The Crystal Structure of a c-Src Complex in an Active Conformation Suggests Possible Steps in c-Src Activation. Structure. 2005 Jun;13(6):861-871. doi: 10.1016/j.str.2005.03.012. pmid:15939018

Young MA, Shah NP, Chao LH, Seeliger M, Milanov ZV, Biggs WH, et al. Structure of the Kinase Domain of an Imatinib-Resistant Abl Mutant in Complex with the Aurora Kinase Inhibitor VX-680. Cancer Research. 2006 Jan;66(2):1007-1014. doi: 10.1158/0008-5472.CAN-05-2788. pmid:16424036

Cowan-Jacob SW, Fendrich G, Floersheimer A, Furet P, Liebetanz J, Rummel G, et al. Structural biology contributions to the discovery of drugs to treat chronic myelogenous leukaemia. Acta Crystallographica Section D: Biological Crystallography. 2006 Dec;63(Pt 1):80-93. doi: 10.1107/s0907444906047287

Levinson NM, Kuchment O, Shen K, Young MA, Koldobskiy M, Karplus M, et al. A Src-Like Inactive Conformation in the Abl Tyrosine Kinase Domain. PLoS Biol. 2006 May;4(5):e144. doi: 10.1371/journal.pbio.0040144. pmid:16640460

Kannan N, Neuwald AF. Did Protein Kinase Regulatory Mechanisms Evolve Through Elaboration of a Simple Structural Component? Journal of Molecular Biology. 2005 Sep;351(5):956-972. pmid:16051269 doi: 10.1016/j.jmb.2005.06.057

Foda ZH, Shan Y, Kim ET, Shaw DE, Seeliger MA. A dynamically coupled allosteric network underlies binding cooperativity in Src kinase. Nature Communications. 2015 Jan;6:5939. doi: 10.1038/ncomms6939. pmid:25600932

Ozkirimli E, Yadav SS, Miller WT, Post CB. An electrostatic network and long-range regulation of Src kinases. Protein Science: A Publication of the Protein Society. 2008 Nov;17(11):1871-1880. doi: 10.1110/ps.037457.108

Su JL, Yang PC, Shih JY, Yang CY, Wei LH, Hsieh CY, et al. The VEGF-C/Flt-4 axis promotes invasion and metastasis of cancer cells. Cancer Cell. 2006 Mar;9(3):209-223. doi: 10.1016/j.ccr.2006.02.018. pmid:16530705

Fenn TD, Schnieders MJ, Brunger AT, Pande VS. Polarizable Atomic Multipole X-Ray Refinement: Hydration Geometry and Application to Macromolecules. Biophysical Journal. 2010 Jun;98(12):2984-2992. doi: 10.1016/j.bpj.2010.02.057. pmid:20550911

Schnieders MJ, Fenn TD, Pande VS. Polarizable Atomic Multipole X-Ray Refinement: Particle Mesh Ewald Electrostatics for Macromolecular Crystals. Journal of Chemical Theory and Computation. 2011 Apr;7(4):1141-1156. doi: 10.1021/ct100506d. pmid:26606362

Nagar B, Hantschel O, Young MA, Scheffzek K, Veach D, Bornmann W, et al. Structural Basis for the Autoinhibition of c-Abl Tyrosine Kinase. Cell. 2003 Mar;112(6):859-871. doi: 10.1016/S0092-8674(03)00194-6. pmid:12654251

Shan Y, Seeliger MA, Eastwood MP, Frank F, Xu H, Jensen MØ, et al. A conserved protonation-dependent switch controls drug binding in the Abl kinase. Proceedings of the National Academy of Sciences. 2009 Jan;106(1):139-144. doi: 10.1073/pnas.0811223106

Alexov EG, Gunner MR. Incorporating protein conformational flexibility into the calculation of pH-dependent protein properties. Biophys J. 1997;72:2075-2093. doi: 10.1016/S0006-3495(97)78851-9. pmid:9129810

Georgescu RE, Alexov EG, Gunner MR. Combining conformational flexibility and continuum electrostatics for calculating pKas in proteins. Biophys J. 2002;83:1731-1748. doi: 10.1016/S0006-3495(02)73940-4. pmid:12324397

Song Y, Mao J, Gunner MR. MCCE2: Improving protein pKa calculations with extensive sidechain rotamer sampling. J Comput Chem. 2009;30:2231-2247. doi: 10.1002/jcc.21222. pmid:19274707

Adams JA, Taylor SS. Divalent metal ions influence catalysis and active-site accessibility in the cAMP-dependent protein kinase. Protein Science. 1993 Dec;2(12):2177-2186. doi: 10.1002/pro.5560021217. pmid:8298463

Sousa SF, Fernandes RA, Ramos MJ. Molecular Dynamics Simulations: Difficulties, Solutions and Strategies for Treating Metalloenzymes. In: Paneth PaDDA, editor. Kinetics and Dynamics: From Nano- to Bio-Scale. vol. 12 of Challenges and Advances in Computational Chemistry and Physics. Berlin: Springer Science & Business Media; 2010. p. 530

Wang J, Wolf RM, Caldwell JW, Kollman PA, Case DA. Development and testing of a general AMBER force field. J Comput Chem. 2004;25:1157-1174. doi: 10.1002/jcc.20035. pmid:15116359

Wang J, Wang W, Kollman PA, Case DA. Automatic atom type and bond type perception in molecular mechanics calculations. J Mol Graph Model. 2006;25:247-260. doi: 10.1016/j.jmgm.2005.12.005