Computational protocol: IBiSA_Tools: A Computational Toolkit for Ion-Binding State Analysis in Molecular Dynamics Trajectories of Ion Channels

Similar protocols

Protocol publication

[…] IBiSA_tools consists of a C++ program and several Python scripts (version 2.6 or later is required). For visualization of the results, some R scripts are also included. The overview of this toolkit is provided in .First, a series of MD trajectory files and an initial structure file are needed as input (—fn-trr setting in the configuration file). Because the current version accepts only .trr format (GROMACS), trajectory files written in other formats have to be converted into .trr format with some other tools, e.g., VMD [] or MDAnalysis []. See the documentation for brief instructions on the file conversion. The C++ program named TraChan parses a trajectory and converts Cartesian coordinates of ions into the coordinates along the channel pore axis, which is defined by two sets of atoms as the bases; the line through the center of the first set to that of the other set is defined as the pore axis (; the settings—pore-axis-basis-from and—pore-axis-basis-to in the configuration file specify these sets of atoms). The output files record positions of ions along the pore axis and the radial distance from the pore axis. The former can be visualized with an attached R script (). The histogram of frequency of ion coordinates across the pore axis is analyzed by means of ion_histogram.py and an attached R script (). When the histogram shows a well-separated multimodal distribution like , boundaries of a series of ion-binding sites can be defined by the pore axis coordinates at the local minimum between every successive two peaks (the coordinates of boundaries are specified by the keyword—site-boundary in the configuration file). If the peaks cannot be separated by certain values on the pore axis, then the definition of the pore axis should be revised. Note that, there is no standard to determine whether a hypothetical peak is treated as separated two peaks or a broad single peak based on the density distribution. In such twilight cases, ion-binding sites can be defined in terms of physicochemical perspective. For example, in a K+ channel, each ion-binding sites is formed as an octahedral site by the backbone oxygen atoms of the SF. In addition, the boundary of ion-binding sites on the radial coordinate also can be defined (—site-max-radius setting in the configuration file). This setting is not important for usual cases because the channel pore is narrowed by the protein structure and the range of ion-binding sites on the radial coordinates is sterically defined by the protein structure. Details of the other settings are described in the software documentation. Next, program site_occupancy.py analyzes what ion-binding sites have ions in each snapshot. We define the ion-binding state of each snapshot as a set of ion-binding sites with ions. For example, when a series of ion-binding sites are named as 0-origin numbers from the extracellular side to the intracellular side () and K+ ions bind to the ion-binding sites 0, 2, and 4 at the same time, the ion-binding state of this snapshot is referred to as K:0:2:4. A transition from K:0:2:4 to K:2:4 means that the outermost ion is released from the pore into the extracellular fluid. Program analyze_ion_path.py shows what ion-binding sites were passed by each ion during each conduction process.To summarize ion-binding states observed in a trajectory, the network diagram called ion-binding state graph is generated (). A node represents an ion-binding state, and an edge between two nodes denotes transition between them. The probabilities of each node and edge are depicted as a size of the node and width of the edge, respectively. This graph can be drawn with standard network analysis tools, e.g., Cytoscape [] via an output file generated by analyze_site_state.py written in the standard graph format Graph Modeling Language (GML). In the ion-binding state graph, a simulation trajectory is expressed as a single path, and the latter can be decomposed into a series of cyclic paths, each of which starts from the most stable state and reaches the same state. When a cyclic path includes both ion association and dissociation processes, the path corresponds to an ion conduction event. A list of cyclic paths is obtained by means of program extract_cycles.py.To characterize patterns of ion conduction events, all the cyclic paths are compared and classified into some groups. Program convert_state_characters.py assigns a unique character to each ion-binding state. For example, the * character (asterisk) is assigned to the most stable ion-binding state, which is K:0:2:4. On the basis of this assignment, program cycle_to_sequence.py converts cyclic paths into sequences. The similarities between ion conduction events are assessed on the basis of the sequence alignment, which is performed by the dp_align.py program. This alignment task requires a matrix of similarities among ion-binding states (characters), which is generated by make_score_matrix.py. This program uses a simple similarity definition: the similarity score between the same states is 1.0, that between states with the same number of ions is 0.5, otherwise it is 0.0. Program align_similarity.py generates a similarity matrix of ion conduction events. This matrix can be easily analyzed by data analysis software; an attached R script performs hierarchical clustering and draws a dendrogram (). As shown in our previous reports, the clustering objectively categorizes ion conduction events into several different mechanisms [,]. […]

Pipeline specifications

Software tools GROMACS, VMD, MDAnalysis
Application Protein structure analysis
Diseases Genetic Diseases, Inborn