RIP-MD: a tool to study residue interaction networks in protein molecular dynamics

Input and pre-processing steps

As previously mentioned, RIP-MD takes as input either a dynamic or a static protein structure (Fig. 1A). Structural analyses within RIP-MD are handled by MD analysis (Michaud-Agrawal et al., 2011; Gowers et al., 2016), a python library to manipulate MD trajectories. The first action in this step is to delete heteroatoms maintaining only protein atoms. If the user wants to add missing atoms, such as hydrogens, these are added based on internal coordinates with PDB2PQR (Dolinsky et al., 2007, 2004). After all these steps, additional parameters such as (partial) charges, Lennard–Jones parameters and secondary structure are either defined or calculated. In the case of MD trajectories, partial charges are assigned from the topology file, that is, PSF file, employed to perform the simulation. Charges for static structures and Lennard–Jones parameters are assigned from the parameters files of the CHARMM force field (MacKerell et al., 1998). Secondary structure and solvent accessibility for each AA are defined using DSSP (Joosten et al., 2011; Kabsch & Sander, 1983). For further details about input files and their format, please refer to the RIP-MD user manual.

Interactions defined in RIP-MD

Residue interaction networks in protein molecular dynamics defines several types of interactions between the AAs in a protein structure: C_α contacts, HBs, SBs, disulfide bonds, cation-π, π–π, Arg–Arg, Coulomb, and vdW contacts. Parameters defining each interaction, together with their mathematical formulation is explained in detail in Text S1 and resumed in Table 1. Importantly, users may employ the provided by-default parameters defining each of these interactions or set them accordingly to his/her needs.

RIP-MD versions and availability

Residue interaction networks in protein molecular dynamics is available free of charge at http://www.dlab.cl/ripmd in three options:

• Standalone program: This version is thought for users who want to take advantage of high-performance computing architectures to perform analysis of very large systems, such very long MD trajectories that are impossible to handle otherwise.

• Visual molecular dynamics (Humphrey, Dalke & Schulten, 1996) plugin: This version benefits from the graphical interface provided by VMD. This plugin performs a system call to execute the standalone RIP-MD program.

• Webserver: This form of RIP-MD is for those users who want to analyze a single PDB structure without installing the software locally. In this web-server, users first upload their selected structure and then, after few steps, the results can be easily downloaded.

It is important to note that results generated by the web-server and by the stand-alone version of the program are compatible with the VMD plugin, making possible to load and display these results into the VMD graphical user interface. The installation guide for both the stand-alone and plug-in versions, together with the user manual, is available at the RIP-MD webpage. In addition, users can also download from the webpage a virtual machine with RIP-MD already available which only requires the installation of all proprietary software that cannot be distributed by us.

Other scripts employed to generate the figures shown in this article are also available at http://www.dlab.cl/ripmd or directly at https://github.com/networkbiolab/supp_script_ripmd.

Output files

Residue interaction networks in protein molecular dynamics provides two types of output: interaction graphs and Pearson correlation plots.

Comparison with other methods

A comparison with two other approaches, Carma (Glykos, 2006) and MD-TASK (Brown et al., 2017), used on our first example (the MD2 pocket closure, see next section) is explained in detail in Text S1. These methods were chosen due to their availability, they are easy to install and they consider both inter and intra chain interactions, even thought both methods only generate a single type of interactions between residues (C_α and C_β contacts for Carma and MD-TASK, respectively). Interested readers can also find in Text S1 a brief description of other methods or libraries that can be used to generate RINs for MD simulations, even if, as Carma and MD-TASK, they use MDs instead of non-covalent interactions in their RIN definitions.

MD2: an example of a conformational change

MD2 is a soluble protein part of the innate immune response in humans. Upon binding of lipopolysaccharides belonging to Gram-bacteria, MD2 triggers an immune response via the interaction with Toll-Like receptors 4 (Dziarski et al., 2001; Ohto et al., 2007). As shown by MD simulations (DeMarco & Woods, 2011; Garate & Oostenbrink, 2013), MD2 possesses a hydrophobic cavity that rapidly closes upon ligand removal (Fig. 2). The later is a clear example of a conformational change switching the RIN, suitable to be visualized with RIP-MD. To do so, we analyzed MD simulations of apo-MD2 obtained from (Garate & Oostenbrink, 2013). Briefly, a 20 ns MD of a solvated MD2 was divided into three windows (Fig. 3) which characterize the closing process. The closing event along the MD simulation was projected onto the collective variable described in Eq. (1). (1) $$D_0=\frac{1}{N}{\displaystyle \sum _{n=1}^N\sqrt{{({\text{COMβ}}_n-{\text{COM}}_{\text{MD}2})}^2}}$$ where n stands for each of the 10 β strands of MD2, and COMβ_n and COM_MD2 are the centers of mass of strand n and MD2, respectively. As reported in (Garate & Oostenbrink, 2013), D₀ unambiguously differentiates the three well defined stages of the closing process. We used RIP-MD to obtain RINs for both HBs and vdW interactions between AAs, employing default parameters (see Table 1 and Supplementary Material). For each window, Pearson correlations were calculated for the number of HBs and vdW interactions for any given pair of AAs.

Residue interaction networks obtained during the closing process of MD2 (Fig. 3) for both HBs and vdW interactions between AAs, reflecting absolute Pearson correlations values |r| ≥ 0.5 are shown in Fig. 4. There is a higher amount of correlated pairs for both HBs and vdW interactions at the initial stage, as the closing process continues the number of correlated pairs decay (see Fig. 4). The latter reveals an initial concerted action that triggers the closing event, but is lost after closure, a clear indication of the unspecific nature of non-polar interactions which are dominant within the hydrophobic cavity of MD2. To further quantify these changes, Table 2 presents the number of correlated pairs at each window normalized by the amount of correlated pairs of the first window. Overall, a monotonic decrease upon window increment, not withstanding any correlation threshold, is observed. Slight deviations from the latter regarding HBs occur at Pearson values equal or above 0.7 (third and fourth columns of Table 2). This increase is due to the formation of stable H-bonds upon MD2 closure (see Table 2). On the other hand, the number of vdW interactions in Table 2 is higher than that of H-bonds, as the latter are specific interactions depending on both the relative distance and orientations of the participant residues. Finally, while most of the pairs of residues are uncorrelated, the β strands tend to be coordinated with respect to the formation and break of both types of interactions, showing how all secondary structure elements are kept in the closure of MD2.

Interestingly, this behavior is not expected for other conformational changes such as folding, in which the formation of polar interactions will lead to the appearance of highly correlated pairs, for example, the formation of an alpha-helix. In this way, RIP-MD served to quantify and visualize the counter-intuitive idea that a large conformational change can indeed lead to lower correlations, a consequence of the nature of the interactions that dominate a given structural process; in this case, non-polar VdWs interactions among residues within the MD2 cavity.

Gap-junction channel: an example for “inter-monomeric” and “inter-molecular” interactions in a large system

GJCs are intercellular hydrophilic channels connecting the cytoplasm of two adjacent cells (Villanelo et al., 2017). GJCs allow the exchange of water, ions and small molecules of up to 1 kDa (Söhl & Willecke, 2004; Araya-Secchi et al., 2014). GJCs are formed by the extracellular docking of two HCs (Fig. 2C), where each HC is formed by six CXs monomers (Fig. 2D). In this example, we studied a 20 ns all-atom MD simulation of a complete GJC formed by the human CX26 (Escalona et al., 2016). This molecular system contains 12 identical CX monomers in total, each one comprising 226 AAs. Due to the high computational cost of studying vdW and Coulombic interactions in such a large system, these were disregarded in this example. In addition, due to the covalent nature of disulfide bridges these were not considered in the following analyses. RIP-MD was run with default parameters, keeping only those appearing in at least 75% of the simulation time.

Table 3 exhibits the interactions found by RIP-MD over the human CX26 GJC simulation. All interaction types are found to be highly stable over the simulation with an overall persistence of 95–100%, indicating that this molecular structure remains structurally stable at least over the 20 ns of the MD simulation. Since C_α contacts describe spatial relationships between AAs with a plenty of available methods to study their interaction (Sethi et al., 2009), only the other types of interactions calculated by RIP-MD are further analyzed.

Table 4 shows those interactions occurring at least during the 75% of the MD simulation for all the following cases: interactions between AAs of the same chain for each of the two HCs (Intra-chain), between chains of the same HC (Inter-chain) and between AAs of different HCs (Inter-HCs). It should be noted that HBs and SBs are interactions appearing mainly between different chains, while π–π and Arg–Arg interactions are both intra-monomeric interactions. Interestingly, only HBs appear at the interface of the two docked CXs (Fig. 5). The role of these non-bonded interactions in the maintenance of the GJC quaternary structure were also described in the original paper reporting the CX26 crystal structure (Maeda et al., 2009).

The information shown in the bottom section of Table 4 was then compared with the original static structure of the GJC after energy minimization (model based on PDB ID 2ZW3, see Araya-Secchi et al. (2014) for details). The top section of Table 4 resumes the number of interactions present in this structure. The first thing to notice when comparing the dynamic and static versions of the network is the lower number of interactions in the dynamic network. This reduction can be due to their low stability at physiological temperatures, the replacement of intra-H-bonds with H-bonds established with solvent molecules or a poor description of a given interaction by the force-field, for example, π–π interactions. Notoriously, the number inter-CXs HB interactions increases from 9 in the static structure to 26 in the dynamic network, highlighting the relevance of HBs to maintain the GJC complex. It is also very relevant that most of the inter-CXs HBs observed in the non-dynamic RIN are absent in the dynamic version, with only two HBs appearing in the entire simulation (Table 5). Regarding the other types of interactions, most of the intra and inter-chain HBs, SBs, Cation-π, and π–π interactions appear in the dynamic and non-dynamic RINs.