Escherichia coli FtsZ molecular dynamics simulations

Abstract

Earlier molecular dynamics studies of the FtsZ protein revealed that the protein has high intrinsic flexibility which the crystal structures cannot reveal. However, the input structure in these simulation studies was based on the available crystal structure data and therefore, the effect of the C-terminal Intrinsically Disordered Region (IDR) of FtsZ could not be observed in any of these studies. Recent investigations have revealed that the C-terminal IDR is crucial for FtsZ assembly in vitro and Z ring formation in vivo. Therefore, in this study, we simulated FtsZ with the IDR. Simulations of the FtsZ monomer in different nucleotide bound forms (without nucleotide, GTP, GDP) were performed. In the conformations of FtsZ monomer with GTP, GTP binds variably with the protein. Such a variable interaction with the monomer has not been observed in any previous simulation studies of FtsZ and not observed in crystal structures. We found that central helix bends towards the C-terminal domain in the GTP bound form, hence, making way for polymerization. A nucleotide dependent shift/rotation of the C-terminal domain was observed in simulation time averaged structures.

Communicated by Ramaswamy H. Sarma

Keywords:

• FtsZ
• Intrinsically Disordered Peptide
• Gromacs Molecular Dynamics

1. Introduction

Many crystal structures of FtsZ are available in the PDB (the Protein Data Bank). For example, Pseudomonas aeruginosa and Methanocaldococcus jannaschii FtsZ crystal structures (PDB ID: 2VAW (Oliva et al., 2007) and 1W5A/1W5B (Oliva et al., 2004)). The protein has a conserved globular core which is approximately 300 amino acid residues. This core comprises of the ordered N- and C-terminal domains. GTP binds to the nucleotide binding site present in the N-terminal domain. The protein also contains a short N-terminal intrinsically disordered peptide (IDP) and a much longer C-terminal intrinsically disordered region (IDR). The N-terminal IDP is present in the crystal structure of P. aeruginosa FtsZ, 2VAW. However, the C-terminal IDR structure has not yet been resolved by X-ray crystallization. This sequence is not a conserved protein sequence, and its length is found to vary considerably across different bacterial species (Vaughan et al., 2004). However, the extreme C-terminus consists of a conserved ∼ 11 residue C-terminal helix which interacts with certain regulatory proteins like ZipA. The crystal structure of the extreme C-terminal residues 375–383 of Escherichia coli FtsZ co-crystallized with ZipA is available (PDB ID: 1F47 (Mosyak et al., 2000)).

In a previous study, FtsZ crystals obtained from different species were compared (Oliva et al., 2007). The root mean square deviation (RMSD) between the available structures and the angle of rotation of the C-terminal domain was calculated. The RMSD or difference in rotation angles were found to be insignificant to be indicative of the possibility of different FtsZ conformations in these crystal structures. The conclusion drawn from the comparison was that there could be a possibility of inter-species differences in the crystal structures, however, the existence of nucleotide dependent different conformations of FtsZ appeared to be unlikely.

By using electron microscopy, the FtsZ protofilaments have also been visualized (Lu et al., 2000). It was observed that FtsZ protofilaments may be straight or bent depending on nucleotide binding. In vitro polymerization studies support a more ‘straight’ conformation of the GTP bound protofilament and ‘curved’ GDP bound protofilaments (Lu et al., 2000).

Molecular dynamics simulations support similar ideas. While significant conformational changes in the monomer as such, had not observed in simulations (Hsin et al., 2012; Natarajan & Senapati, 2013), the studies definitely revealed a high intrinsic flexibility of the protein. In their FtsZ monomer simulations (Natarajan & Senapati, 2013) reported a ‘twisted bending’ of the core helix and the possibility that the ‘H6–H7’ region (top portion of the central helix) acts a switch which recognizes the nucleotide binding state of the protein. Principal component analysis (PCA) revealed that in the GTP bound simulations of the protein, the top portion of the central helix moves upwards (towards the C-terminal domain of the top subunit in the longitudinal protofilament) and the HC2, HC3 of C-terminal domain move downwards to make protofilament longitudinal contacts. Even though their simulations revealed important nucleotide dependent structural changes in the monomer, a significant conformational change in the monomer was not observed. RMSD between the GTP and GDP forms were very low, comparable to the effects of X-ray crystallization resolution. In the FtsZ dimer simulations (Hsin et al., 2012), the nucleotide dependent change in protofilament curvature was attributed only to the loss of monomer-monomer contact.

Recent studies shed light the role of the C-terminal linker (CTL) (i.e., the C-terminal IDR minus the extreme C-terminal conserved helix) in Z-ring assembly (Buske et al., 2015; Buske & Levin, 2013). Removal of the CTL, in Bacillus subtilis leads to a filamentous phenotype in which Z-ring assembly does not occur. In vitro, the removal of CTL resulted in a much lower GTPase activity (2.5 fold lower GTPase) and an increase in FtsZ critical concentration (by 1.5 fold) and inhibited protofilament formation. It was thus concluded that intrinsic flexibility is crucial for FtsZ function, and replacing the region with ‘rigid’ alpha-helix hinders FtsZ assembly in vitro and cell division in vivo. It was found that ‘FtsZ CTLH’ (FtsZ in which the CTL had been replaced by rigid helical repeats) could form only small oligomers in vitro. CTLs of varying lengths (25–100 residues) were functional in vivo.

Since CTL is required for FtsZ assembly in vitro (in the absence of regulatory proteins), therefore, it appears to be likely that the CTL might be playing an important role during the formation of protofilaments. Therefore, in the present study, we investigated the actual role of the FtsZ IDR by applying the molecular dynamics simulation method. Simulations were performed using a homology model for the E. coli FtsZ in which the C-terminal IDR has been included. This wild-type monomer was simulated without nucleotide, with GTP or with GDP (guanosine diphosphate). In all the simulations, it was observed that the homology model retains its overall fold of the core domains and that the IDR is very flexible. Differential flexibilities of the IDR were observed in the wild-type monomer simulations with GDP and GTP which re-affirms the role in protofilament formation as found in (Buske et al., 2015; Buske & Levin, 2013). In agreement with the simulations of Natarajan and Senapati (2013), simulations also reveal a significant bending of the central helix in GTP simulations and C-terminal domain rotation.

2. Methods

2.1. Homology model of the E. coli FtsZ protein (wild-type, WT)

For constructing a model of the WT E. coli FtsZ protein, the gene: ftsZ (EC2860050_0094) was identified (NCBI accession is AQDL01000001.1, locus tag is EC2860050_0094). The primary sequence was obtained using the ExPASy Translate tool (Gasteiger et al., 2003). The obtained primary sequence was aligned with primary sequences of FtsZ proteins from other bacterial species, for which crystal structures are already available in the Protein Data Bank. BLAST-p (Altschul et al., 1990) alignment scores are provided in Table 1.

Table 1. Alignment of E. coli FtsZ protein (obtained by translating Gene: ftsZ (EC2860050_0094)) with various subject sequences using the Standard Protein Blast algorithm.

Species and Transcript ID	Max. score	Total score	Query cover (%)	E-value	Ident. (%)	Positives (%)	No. of matches
P. aeruginosa PAM18_4498	392	392	93	4e − 138	57	69	1
M. tuberculosis MTCTRI2_2186	253	253	93	3e − 84	45	63	1
B. subtilis BEST7613_2297	279	298	98	8e − 94	47	64	2
T. maritima Tmari_0838	254	269	82	6e − 85	48	66	2
A. aeolicus aq_525	231	231	77	6e − 76	46	63	1
S. aureus SA21334_0269	270	285	81	9e − 91	51	72	2
S. epidermidis SEVCU045_1562	266	266	91	3e − 89	45	64	1
H. volcanii ADE02725	230	244	81	2e − 75	45	65	2

The P. aeruginosa FtsZ gave the largest total score, query cover and highest E-value. Therefore, the crystal structure 2VAW of P. aeruginosa FtsZ was selected as the primary template for homology modelling. E. coli FtsZ protein sequence is 383 amino acids. In all available FtsZ crystal structures including 2VAW, structural information is available for the ordered residues of the globular core only. In 2VAW, the short disordered peptide at the N-terminal is also present in addition. Therefore, this structure was used as a template for residues 2–316 of the E. coli FtsZ protein.

However, the IDR at the C-terminal is not resolved in any of the crystal structures. Since, residues 2–316 of E. coli FtsZ align with residues 2–316 of the P. aeruginosa FtsZ sequence, therefore, it may be assumed that residues 317–383 in E. coli FtsZ correspond to the C-terminal IDR (including the extreme C-terminal conserved helix).

On analysing the primary sequence (PrDOS probability described later) it was observed that the residues of the IDR have the highest disorder probability. Therefore, during molecular dynamics simulation, high RMSF (root mean square fluctuation) in this region was expected. Since this region is expected to have exceptionally high flexibility, the primary expectation from the template-based modelling approach for this region was only to predict the local secondary structures in this region. We selected templates for this region by aligning short sequences from the IDR with other proteins in the PDB database. PDB Extended Blast-P search (E. value cut off up to 10) for E. coli FtsZ query sequence gave an MP1-p14 complex (PDB accession: 1SKO, molecule: Mitogen-activated protein kinase kinase 1 interacting protein 1, chain A: 143 amino acids (Lunin et al., 2004)). A 41 residue sequence in this protein, from 21 to 61 (crystal residues), aligned with residues 342–382 in E. coli FtsZ. 48% positive alignment was observed for this sequence. For building homology model, residues 21–62 of the 1SKO structure were used as template for residues 342–383 of E. coli FtsZ. To identify suitable templates for residues 317–341, shorter peptides from the C-terminal flexible/disordered region were used as query sequences. For the query sequence 325–360 in E. coli FtsZ, one of the results was residues 172–207 in Marinibacillus marinus Adenylate kinase protein, PDB accession 3FB4 (Davlieva & Shamoo, 2009). For residues, 317–324 also, we used neighbouring residues from the 3FB4 sequence. Therefore, including additional 11 neighbouring residues, residues 161–207 of the 3FB4 structure were used as template for residues 317–360 of E. coli FtsZ (48% positive alignment, included 3 gaps in the match). The C-terminal residues 367–383 of the E. coli FtsZ protein has already been crystallized with the protein ZipA (Mosyak et al., 2000). It was used as a template for the corresponding residues.

The templates used for constructing E. coli FtsZ homology model are listed in Table 2. Modeller (Fiser et al., 2000; Marti-Renom et al., 2000; Sali & Blundell, 1993) program was used for building the homology model. One model was generated. The model generated includes residues Phenylalanine (residue no. 2)—Aspartate (residue no. 383) (Figure 1).

Figure 1. E. coli FtsZ model built using templates given in Table 2.

Table 2. Template crystal structures for constructing E. coli FtsZ homology model.

S. No.	PDB	Protein	Aligned residues from E. coli FtsZ
1	2VAW (Oliva et al., 2007) (2–316)	P. aeruginosa FtsZ	2–316
2	1SKO (Lunin et al., 2004) (A:21–62)	Homo sapiens Mp1-p14 kinase	342–383
3	3FB4 (Davlieva & Shamoo, 2009) (A:161–207)	M. marinus Adenylate kinase	317–360
4	1F47 (Mosyak et al., 2000) (A:1–17)	ZipA-FtsZ E. coli	367–383

The E. coli FtsZ homology model and the coordinates for nucleotides GDP and GTP are provided in the supplementary information, and the AMBER 94 forcefield files are also provided. Nucleotides were allowed to bind freely by placing inside the nucleotide binding site or close to it in the VMD (n.d.a, n.d.b) software tool by performing translations and rotations of the nucleotides. Initial position of the nucleotides is provided in Figure 2. Two independent repeat simulations were performed with GDP. Five independent repeat simulations were performed with GTP using two different starting orientations of the nucleotide.

Figure 2. Position of the nucleotides in NPT equilibrated structures: GDP in 2VAW reference (blue), GDP in 2 simulations (red), 2 different orientations of GTP in 5 simulations (green).

2.2. Molecular dynamics simulations

Atomistic molecular dynamics simulations were performed using GROMACS version 4.6.5 (Hess et al., 2008). Simulations were performed in the AMBER94 force field (DePaul et al., 2010; Sorin & Pande, 2005). AMBER parameters for nucleotides developed by Meagher et al. (2003) were used. The protein molecule was centred in the cubic box, at a minimum distance of 1.5 nm from box edges. Solvent (water) molecules were added in the coordinate file. SPC-E (Single Point Charge Extended) water model configuration was used (Berendsen et al., 1987). Mg²⁺ and Cl⁻ ions were added to neutralize the simulation box and at a minimum concentration of up to 10 mM (any of which was higher). Energy minimization was performed using the steepest descent minimization algorithm until the maximum force on any atom in the system was less than 1000.0 kJ/mol/nm. A step size of 0.01 nm was used during energy minimization. A cut-off distance of 1 nm was used for generating the neighbour list and this list was updated at every step. The electrostatic forces were calculated using Particle-Mesh Ewald method (PME) (Darden et al., 1993). A cut-off of 1.0 nm was used for calculating electrostatic and Van der Waals forces. Periodic boundary conditions were used. A short 100 ps NVT equilibration was performed. During equilibration, the protein molecule was restrained. Leap-frog integrator MD simulation algorithm (Hockney et al., 1974) was implemented with a time-step of 2 fs. All the bonds were constrained by the LINear Constraints Solver (LINCS) constraint-algorithm (Hess et al., 1997). Neighbour list was updated every 10 fs. A distance cut-off of 1 nm was used for calculating both electrostatic and van der Waals forces. The two groups, that is, protein and non-protein (solvent, ligand and ions) were coupled with the modified Berendsen thermostat (Berendsen et al., 1984) set at 300 K. The time constant for temperature coupling was set at 0.1 ps. Long range dispersive corrections were applied for both energy and pressure. The output coordinates were saved either every 2 ps or every 5 ps (in .xtc compressed trajectory format). A short 100 ps NPT equilibration was performed using the Parrinello-Rahman barostat (Nos’e & Klein, 1983; Parrinello & Rahman, 1981) algorithm, with time constant of 2 ps, was applied to maintain pressure at a constant value of 1 bar. For the production run a 100 ns MD simulation in an NPT ensemble was implemented.

2.3. PCA analysis methodology

PCA analysis was carried out using Gromacs commands “covar” and “anaeig” (Borkotoky & Murali, 2018). The coordinates of the protein molecule can be generated from the output trajectory and saved in a .gro file format (example: 100 ns_Prot.gro). For the covariance analysis, the trajectory must be fit on the ordered protein (residues 12 to 311 Main chain & beta carbon atoms). Then the entire protein Main chain & beta carbon can be used for covariance analysis. The filtered trajectory along an Eigen vectors using the Gromacs “anaeig” command. The filtered trajectory can then be visualized.

For example:

mpirun -np 1 gmx_mpi covar -f 100ns_Prot.gro -s Struct_mass.pdb -n index.ndx

mpirun -np 1 gmx_mpi anaeig -v eigenvec.trr -eig eigenval.xvg -f 100ns_Prot.gro -s Struct_mass.pdb -filt PCA1.gro -first 1 -last 1 -n index.ndx

In the current study, since the idea is to find the large, concerted motions between the core protein domains and the IDP, the results of the PCA analysis are found to be the same even if the C-alpha only atoms are used in the analysis.

3. Results and discussion

3.1. Protein DisOrder prediction system (PrDOS) prediction of disordered residues in the E. coli FtsZ protein sequence

Algorithms such as PrDos predict disordered regions in proteins (Ishida & Kinoshita, 2007). PrDos predictions are based on (a) local amino acid sequence information and (b) alignments with homologues for which structures have already been determined. Combination (based on weighted average) of these two independent predictions gives the probability of disorder. The threshold probability above which a residue is said to be disordered is 0.5. For the E. coli FtsZ protein sequence, the stretch of residues 319–368 was predicted as a disordered region in the protein (2 state prediction results in Figure 3a). The short N-terminal peptide of ∼ 8 to 10 residues, is found to be in a helical conformation in some FtsZ molecules (e.g., 1W5A), in a loop (e.g., 2VAW) or could not be resolved in some structures (e.g., 2R6R). According to PrDOS predictions these residues, residues 1–8, also have disorder probabilities greater than 0.5.

Figure 3. (a) PrDos 2 state prediction (False positive rate 5%) Red: Disordered residues, Black: Ordered residues. (b) Ramachandran distribution map of I) model phi and psi angles II) intrinsically disordered region in the protein.

3.2. Model assessment

The RMSD of the model and the 2VAW template based on backbone alignment was 3.3 Å. The IDR was modelled using the templates 3FB4 residues 161–207 for residues 317–360 of E. coli FtsZ. The Modeller predicted the residues 317–330 in either coil or turn conformations, 331–341 in an α-helix, 342–357 in a β-hairpin loop and 358–362 are in either coil or turn conformations, 363–366 are in a 3₁₀ helix, 367–371 are in a turn conformation, 372–381 are in an α-helix and 382 in a coil. The model was assessed using ProQ (Wallner & Elofsson, 2003) which generates LG and Maxsub scores based on atom-atom contacts, residue-residue contacts and solvent accessibility. LG score for the model was 5.396, that is, a very good model.

The model (ordered core as well as the intrinsically disordered regions) exhibits the Φ, Ψ torsion angles in accordance with the Ramachandran distribution (Figure 3b). Apart from a 10 to 12 residues on the extreme right of the chart, the rest of the residues are all within the allowed regions for the alanine dipeptide (Hollingsworth & Karplus, 2010). In the IDP sequence also conforms to the Ramachandran distribution of Φ, Ψ torsion angles.

3.3. MD visualization

Elaborate visualizations of the simulation trajectories were performed in VMD. To show the evolution of the protein in simulation, the structure of the protein at every 20 ns in the wild-type monomer simulation WT-FtsZ mono-GDP (1) is provided in Figure 4a–f. Structure of the protein at a single time frame from well equilibrated simulation provides considerable information about the simulation. The conformational flexibility of the loop regions which are lost during averaging of structures may be observed from such structures. In some of the simulations we observed a large change in curvature of the central helix of FtsZ. The extent of bending was found to differ between the independent repeat simulations of the same system. An average structure will not show the extent to which such variations occur. Therefore, for the discussion, we have used representative structures obtained from well equilibrated simulations. We use structures from the same time frame 93.9 ns from 100 ns MD simulation (Figure 5a–i).

Figure 4. (a–f) Structure of the protein in the wild-type monomer simulation WT-FtsZ mono-GDP (1) at every 20 ns from t = 0 ns to t = 100 ns during MD simulation: (a) t = 0 ns, (b) t = 20 ns (c) t = 40 ns, (d) t = 60 ns, (e) t = 80 ns, (f) t = 100 ns. The frames were aligned for the backbone atoms residues 12–311.

Figure 5. (a, b) FtsZ simulations without nucleotide, frame considered at 93.9 ns simulation time. (c,d) FtsZ simulations with GDP. (e-i) FtsZ simulations with GTP. VMD ‘new cartoon’ representation is used for the protein and ‘secondary structure’ colour scheme is used. ‘CPK model’ representation is used for nucleotide and coloured based on atom name; and ‘CPK model’ for Mg2+ ion, coloured in black. The bending of the central helix is shown using a red line drawn through the helix.

It may be observed that the structure of the globular core is stable. However, the structure of the C-terminal IDR is not stable and its structure varies during simulations (Figure 4a–f), and between different simulations also (Figure 5a–i), consistent with its intrinsic disorder.

In the two independent repeat simulations with GDP, labelled as WT-FtsZ mono-GDP (1) and (2), GDP is located well inside the nucleotide binding site (Figure 5c, d). Frames at 20 ns interval for the simulation WT-FtsZ mono-GDP (1) given in Figure 4a–f. It may be observed that after t = 20 ns, there is no further significant change in the location of GDP apart from a small change in the phosphate tail between t = 80 ns (Figure 4e) to t = 100 ns (Figure 4f). However, a very high variability is observed both in position and orientation of the GTP nucleotide. For example, in Figure 5f, GTP guanine group points towards the nucleotide binding site and its phosphate tail points outwards but in Figure 5g it is oriented in the opposite direction. The position of GTP was observed to fluctuate continuously during the simulations suggesting that the GTP binding in the monomer form is unstable, probably held by non-specific interactions with the more flexible loop residues of the nucleotide binding site.

The secondary structure of the sequence HL1–S3—(loop between S3 and H2A) (residue nos. 48–72) and S5–H4 loop varied considerably. These residues/regions are labelled in Figure 5a, the secondary structure elements are labelled in Figure 5b. The residues of HL1, form a 3₁₀ helix (Figure 5b, c), α-helix (Figure 5d). It may be observed that in the simulation WT-FtsZ mono-GTP (3) (Figure 5g, i), the short β-strand S3 converted from β-strand to a β-bridge (formed by a single residue). Residues between S3 and H2A form α-helix (Figure 5e) or a 3₁₀ helix (Figure 5c) in some simulations or are present in a loop conformation in others (Figure 5g, h). Similar secondary structure changes may be seen in the S5–H4 loop which forms a 3₁₀ helix in some simulations (Figure 5c–e, g–i). It is likely that the varying conformations of these residues of the nucleotide binding site is essential for allowing nucleotide binding and release. The different conformations observed in the independent repeat simulations with GTP, and at different time points in the same simulation, suggests that these residues GTP interacts with the monomer through many non-specific bonds at the active site.

The curvature of the central helix in the representative figures is shown using a red line along the helix. In the simulations with GTP, the residues 177–188, that is, the N-terminal residues of the central helix bends towards the C-terminal domain, thus opening the nucleotide binding site considerably. This was observed in at least four out of the five simulations, WT-FtsZ mono-GTP (2), (3), (4) and (5) (Figure 5f–i). In the simulation WT-FtsZ mono-GTP (1), the bend was observed for a short duration much before the 93.9 ns frame. In the simulations, WT-FtsZ mono-GTP (2) and WT-FtsZ mono-GTP (4) (Figure 5f, h) a break in the central helix is observed which divides it into two. This is similar to the kinked helix reported in the Staphylococcus aureus FtsZ crystal structure apo form (Matsui et al., 2012). However, the overall structure of the globular core observed in the S. aureus is quite different in which the two β-planes appear to merge and form a single large plane.

On the other hand, in simulations with GDP (Figure 5d) the central helix often bends towards the N-terminal domain (Figure 5d). There are no major distortions in the helix in GDP simulations (changes in secondary structure or break in helix).

It may be observed that residues of the HC2-SC2 loop (245–257), have varying conformations. In Figure 5a, b, only one 3₁₀ helix is present, in Figure 5c–e, g, i, two 3₁₀ helices are observed, in Figure 5f, h an α-helix and a 3₁₀ helix are present.

3.4. Root mean square deviation (RMSD)

To study the stability of the model during molecular dynamics simulations, the RMSD of backbone atoms was calculated. RMSD measures the structural stability of the protein during the course of simulation and may be used to test if there are large structural changes in the protein leading to high RMSD. RMSD was calculated for the backbone atoms of residues 12–311, that is, the globular core. The structure at time, t = 0, that is, after npt-equilibration was used as reference for fitting. The trajectories were fit over backbone atoms of residues 12–311. RMSD was calculated at every 20 ps.

RMSD for simulations is provided in Figure 6a–c. A steady increase in RMSD during the initial 20 ns period was observed. This corresponds to the equilibration/stabilization phase of the simulation. After this equilibration phase, structural drifts in the core globular domain minimize, and a more or less stable RMSD is attained. RMSD fluctuates within ∼0.5 Å. This range is quite similar to RMSD reported in previous studies (Hsin et al., 2012; Natarajan & Senapati, 2013) in which crystal structures were used in simulations.

Figure 6. (a-c) RMSD of backbone atoms of residues 12–311 of the globular tubulin like core of FtsZ in (a) WT-FtsZ mono-empty (1) (red) and WT-FtsZ mono-empty (2) (green). (b) WT-FtsZ mono-GDP (1) (red) and WT-FtsZ mono-GDP (2) (green).(c) WT-FtsZ mono-GTP (1) (red), WT-FtsZ mono-GTP (2) (green), WT-FtsZ mono-GTP (3) (blue), WT-FtsZ mono-GTP (4) (cyan), WT-FtsZ mono-GTP (5) (magenta). Note : For RMSD calculation, the simulation trajectory frames were aligned on backbone atoms of residues 12–311 of the structure obtained after respective NPT equilibrations.

The important conclusions that may be drawn from the RMSD plots are (1) the globular core is stable and does not undergo any unusual structural changes and (2) the presence of CTL does not affect the structure or the stability of the globular core.

3.5. Root mean square fluctuation (RMSF)

Average RMSF per residue was calculated for backbone atoms of residues 2–383 from 20–100 ns that is, after equilibration (Figure 7). The high RMSF of the N-terminal helix H0 is consistent with its intrinsic disorder predicted by PrDOS (Section 3.1). The amino acid residues of the loop HL1 have high flexibility. It may be observed that these residues 48–50 may take conformation of loop (Figure 4a, f), alpha helix (Figure 4d) and 3₁₀ helix (Figure 4b, c, e, h–i). The residues of the loop between S3 and H2A, 59–72, may have high flexibility. The residues of helix HL2, loop HL2-H5 and few N-terminal residues of the central helix H5 have high flexibility (residues 165–175). This is consistent with the changes in secondary structure (α-helix to 3₁₀ helix) and the bending observed in this region of the protein. The residues in this region form contacts with HC3 of the upper subunit during polymerization [M. jannaschii FtsZ dimer 1W5A, 1W5B]. The bending of the central helix in simulations with GTP will make the central helix more accessible for making protofilament contacts. Consistent with the intrinsic disorder of the IDR, RMSF for residue no. 320 onwards were much higher than the rest of the protein. Two to three major peaks are observed in this region (1) between 321 and 341, (2) near 361 and (3) the extreme C-terminal tail which has the maximum RMSF. In between the first two peaks of the IDR there is a sharp decline in RMSF. Such an RMSF profile indicates that there is some loose structural order due to the presence of certain less flexible residues within the IDR. Further, it may be observed that the maximum RMSF for residues in this region in simulations of FtsZ with GTP is less than 0.55 nm but in simulations with GDP it may go up to upto1.64 nm (0.8 nm in the other).

Figure 7. RMSF of protein backbone atoms and averaged per residue from 20–100 ns WT-FtsZ mono-empty (red); WT-FtsZ mono-GDP (blue) WT-FtsZ mono-GTP (green).

Since many interactions with regulatory proteins occur at the CTL, a highly flexible CTL would increase interaction of the GDP bound monomer with regulatory proteins. In contrast, a lower flexibility CTL, would decrease interaction with regulatory proteins, in turn allowing polymerization of the GTP form.

3.6. Hydrogen bonds between nucleotide and protein

For calculating hydrogen bonds the distance cut-off of 3.5 Å and the angle cut-off 25° was used. Hydrogen bonds were calculated during 0–100 ns MD simulation. ‘Hydrogen bond occupancy’ is defined as the percentage time during which the hydrogen bond is present in the simulation. Hydrogen bonds which have occupancy greater than 60% are given in Table 3. In the simulation WT-FtsZ mono-GDP (1), GDP forms hydrogen bonds with Glu137, Arg141, Gly106 and Thr131. In the simulation WT-FtsZ mono-GDP (2), GDP forms 2 hydrogen bonds only with Glu137, Arg141. In the simulation WT-FtsZ mono-GTP (1), GTP forms hydrogen bonds with Arg141 and Lys140. In the simulation WT-FtsZ mono-GTP (2), GTP forms hydrogen bonds with Arg141, Lys140 and Glu137. In the other three simulations, GTP forms hydrogen bonds with Arg141 only.

Table 3. Hydrogen bonds between nucleotide and protein.

S.No.	Simulation	Hydrogen bond with nucleotide	Occupancy (%)
1	WT-FtsZ mono-GDP (1)	GLU137-Side	95
		ARG141-Side	71
2	WT-FtsZ mono-GDP (2)	GLU137-Side	84
		ARG141-Side	66
3	WT-FtsZ mono-GTP (1)	ARG141-Side	75
		LYS140-Side	68
4	WT-FtsZ mono-GTP (2)	ARG141-Side	63
		LYS140-Side	61
		GLU137-Side	61
5	WT-FtsZ mono-GTP (3)	ARG141-Side	86
6	WT-FtsZ mono-GTP (4)	ARG141-Side	82
7	WT-FtsZ mono-GTP (5)	ARG141-Side	68

The residues Arg141, Lys140 and Glu137 belong to the S5-H4 loop. Since both GTP and GDP bind to the same S5-H4 loop residues and GTP can bind from various positions from the outside by forming non-specific hydrogen bonds (Figure 5e, h). Such diverse interactions of GTP would allow the release of GDP on introduction of GTP. Further, the bending of the central helix on the introduction of GTP would also enable the release of bound GDP.

In response to peer review query, a 50 ns simulation of FtsZ with bound GDP and another 50 ns simulation of FtsZ with bound GTP was performed to find the most stable hydrogen bonds that persist for more than 75% of the simulation duration. These bonds were found to be with GLU137 and ARG141 and were formed with both the nucleotides.

3.7. Principal component analysis

Principal component analysis was performed to identify significant concerted motions in different regions of the protein and understand protein dynamics. Covariance analysis was performed for the Cα atoms of the monomer protein. Filtered trajectories for the Cα atoms were generated for the first Eigen vector (Figure 8). The PCA contribution along the first Eigen vector are in the order of 32 and 44% in the simulations with GTP and GDP, respectively.

Figure 8. (a-i) Principal motion along the first eigen vector for (a,b) WT-FtsZ-empty (1) and (2), (c,d) WT-FtsZ-GDP (1) and (2) and (e–i) WT-FtsZ-GTP (1), (2), (3), (4) and (5). The position of atoms at different simulation timesteps is shown using the ‘CPK’ representation in VMD. The ‘trajectory-timestep’ colour scheme in VMD has been used with the ‘RGB’ colour scale such that the first frame (t = 100 ps MD simulation time) is in red and the frame at t = 100,000 ps is in blue.

It may be observed that the Cα atoms of the helix H0, HL1, the S3–H2A loop, the S5–H4 loop, helix HL2, the HL2–H5 loop, the N-terminal residues of the central helix H5, the loop SC1–HC2, helix HC2, loop HC2–SC2, HC3 and the loop SC3–SC4 have higher displacements among the Cα atoms of the residues 1–312, that is, the N-terminal helix H0 and the ordered residues of the N- and C-terminal.

The higher displacements of H0 is expected due to its intrinsic disorder. The nucleotide binding site residues HL1, the S3–H2A loop, the S5–H4 loop, helix HL2, the HL2–H5 loop, the N-terminal residues of the central helix H5 have high displacements. The secondary structure changes were described in Section 3.3 and 3.5.

Rotation of the C-terminal domain may be observed in Figure 8a, b (downwards) and e, g, h (upwards). In GDP bound simulations, the magnitude is less (compared to the two empty- and three GTP bound simulations Figure 8e, g, h). Between the three GTP bound simulations, in which rotation is observed, the domain motion is towards the front in (e), towards the nucleotide binding site in (g) and more perpendicular in (h). The displacement of the Cα atoms of the C-terminal flexible or disordered region are much higher than the rest of the protein.

The average Cα coordinates for the ordered residues were compared after alignment of N-terminal ordered residues (12–164) (Figure 9). It was observed that the C-terminal domain is rotated or shifted more towards the left (towards the nucleotide binding site) in simulations with GTP. In contrast, in the simulations with GDP, the domain is shifted towards the right.

Figure 9. Alignment of average Cα coordinates after superpositioning Cα atoms of N-terminal residues 12–164: WT-FtsZ mono-empty simulations (blue), WT-FtsZ mono-GDP simulations (red), WT-FtsZ GTP simulations (magenta and green).

4. Discussion and conclusions

The FtsZ monomer was simulated without nucleotide, with GTP and with GDP. The E. coli FtsZ protein was simulated. The structure for the protein was constructed by homology modelling. It also includes the C-terminal IDR. The protein model was found to be stable during simulations. The structural stability of the ordered domains was found to be comparable to what was reported in previous studies. The IDR has high structural variability. Residues with much lower flexibility are present within the IDR as seen from the higher RMSF and higher PCA contribution along the first Eigen vector. Such residues with low RMSF provide some loose 3-D form to the IDR. In the simulations with GDP, the IDR is more flexible than in simulations with GTP or without nucleotide. It was observed that GTP binds from the outside of the nucleotide binding site, non-specifically to the long chain residues of the S5-H4 loop (Arg141, Lys140 and Glu137). A significant opening of the nucleotide binding site was observed in the simulations with GTP which is caused by the bending of the core helix. Such a bending of the central helix, would facilitate nucleation polymerization to occur. The structure of nucleus being the GTP bound FtsZ monomer with its open active site at which other GTP bound monomers may add easily.

Principal component analysis suggested concerted motion of the Cα atoms of helix H0, the S3–H2A loop, the S5–H4 loop, the helix HL2, the HL2–H5 loop, the N-terminal residues of the central helix H5, the loop SC1–HC2, the helix HC2, the helix HC3, the loop SC3–SC4 and the IDR. In the simulations of the Aquifex aeolicus FtsZ monomer (Natarajan & Senapati, 2013), concerted motion was observed for residues in the HL2–H5 loop, the S3–H2A loop, the H5–HC1 loop, the helix HC2 and HC3. However, it was observed that in the GTP/GDP bound simulations, HC2 and HC3 move downwards (diagonally opposite from the nucleotide binding site in the GTP bound state, that is, towards the adjacent monomer of the protofilament). In our simulations, we observed that in all the simulations with GTP, the domain shifts forward, upwards or diagonally towards the nucleotide binding site (Figure 8e, g, h). Therefore, our observation suggests that the GTP bound FtsZ favours polymerization by the bending of the central helix, allowing contacts to be formed at the N-terminal domain but the C-terminal domain shifts upwards in the monomer.

An earlier study proposed that the CTL has an important role in transforming FtsZ from an inactive to an active state (Buske et al., 2015; Buske & Levin, 2013) by providing structural flexibility to the N-terminal core and in inter-subunit interaction. The differential flexibilities observed in the presence of GTP and GDP in this study also supports a regulatory role for the IDR. The absence of nucleotide or the presence of GTP higher up in the nucleotide binding site, results in a relatively lower flexible IDR. In contrast, in simulations with GDP, the flexibility is much higher. Our simulations suggest that variable/non-specific GTP binding to the monomer is a key feature of FtsZ dynamics. In the presence of a well bound stable nucleotide, the central helix may be stabilized by interactions with nucleotide but in the presence of GTP, the monomer first allows polymerization to occur at the nucleotide binding site. Contact maps from the average Cα coordinates were generated to find contacts between residues of the IDR (316 to 382) and residues 164–310, HL2 to SC4. Usually 2–3 contacts are formed between residues in 341–361 and (1) in HL2-H5 (2) in HC2, HC2-SC2 loop. For example, in the simulation WT-FtsZ mono-GTP (1), 2 contacts are observed, between 347–350 and 169–171 (in HL2-H5), and between 354–356 and 245–248 (in HC2-SC2 loop) (Table 4). Consistent with this, a low RMSF region is present between residues 341–361 in all simulations. By performing sequence analysis in VMD, it was observed that in residues 341–361, there are 10 non polar residues (47%); in 164–177, there are 6 non-polar residues (43%) and in 243–253, there are 6 non-polar residues (54%). Since this region forms contacts with the HL2-H5 loop and the ordered C-terminal domain, it appears that the C-terminal IDR may be directly involved in bending of the core helix and rotation of the C-terminal domain through hydrophobic interactions. In Figure 9, average structures show a large difference in residues of the HC2–SC2 region (243–253). PCA analysis also shows significant displacements of these residues (Figure 8a, c, f, i). A twisted bending of the central helix in the nucleotide free simulation of A. aeolicus FtsZ monomer was observed in the simulation without nucleotide (Natarajan & Senapati, 2013). Super positioning of average Cα coordinates of their simulations showed that in the same simulation the C-terminal domain shifts towards the nucleotide binding site (similar to the GTP bound simulations of the E. coli FtsZ). Therefore, it appears that hydrophobic interaction between residues of HL2–H5 (164–177) and the helix HC2, HC2–SC2 causes a rotation in the nucleotide free simulation of A. aeolicus FtsZ (Natarajan & Senapati, 2013). However, in the presence of GTP/GDP they do not observe such bending of the central helix, similar to our GDP bound simulations of E. coli FtsZ.

Table 4. Contacts between the IDR and the globular core (contact maps provided in SI).

Simulation	Contacts with HL2-H5	Contacts with C-terminal domain
WT-FtsZ mono-empty (1)	351–355 and 168–171	339–342 and 228–233
WT-FtsZ mono-empty (2)		342–346 and 245–249
WT-FtsZ mono-GDP (1)		353–358 and 246–252
WT-FtsZ mono-GDP (2)	349–354 and 169–172	355–357 and 246–248
WT-FtsZ mono-GTP (1)	347–350 and 169–171	354–356 and 245–248
WT-FtsZ mono-GTP (2)	350–355 and 167–172	355–358 and 247–249
WT-FtsZ mono-GTP (3)	352–354 and 168–173	347 and 238–240
WT-FtsZ mono-GTP (4)	348–352 and 168–173
WT-FtsZ mono-GTP (5)	352–355 and 168–171	353–356 and 247–253

In our simulations, in the nucleotide free form, the hydrophobic interaction between HL2-H5 and HC2-SC2 is limited by steric interactions with the IDR. In the presence of GDP, the central helix is stabilized by the nucleotide. Due to its higher charge, GTP is not able to bind stably in the nucleotide binding site. Since it forms hydrogen bonds from the outside or higher up in the nucleotide binding site with the long chain charged residues of the S5-H4 loop (Arg141, Lys140 and Glu137), it may lead to a differential charge distribution in the nucleotide binding site which allows the IDR to pull the central helix downwards even in the presence of GTP.

Our simulations support recent experiments in which the presence of IDR was found to be essential for FtsZ polymerization (Buske et al., 2015; Buske & Levin, 2013). In the experiments, it was found that in the absence of the C-terminal linker, FtsZ critical concentration increases, however, GTPase activity reduces from 3.82 to 1.08 GTP FtsZ⁻¹min⁻¹, protofilaments were observed after a long delay (15 min after addition of GTP) (Buske & Levin, 2013). Our simulations support that in the absence of the IDR, few nucleus FtsZ would be available to allow polymerization. When the IDR residues are scrambled, then the distribution of non-polar residues is affected, due to which a lower GTPase activity compared to the wild-type is seen in the experiments.

In conclusion, we may say that the IDR is not an open chain of amino acids, it has some loose structural order within it to allow its interaction with the globular core. We observed that the C-terminal IDR forms hydrophobic contacts with the HL2-H5 residues and the HC2, HC2-SC2. Through such interactions with the globular core, it is able to bend the core helix in the GTP bound FtsZ monomer, favouring polymerization.

This study was conducted in 2014 and is among the first attempts in the investigation of intrinsically disordered protein sequences using the molecular dynamics technique. The work demonstrates that the AMBER94 forcefield can represent the dynamics of the protein in conjunction with its IDP sequence successfully. Subsequent work done by scientists in the field also shows that most modern forcefields, including AMBER99SB, correspond well with the NMR data with the CHARMM22* forcefield producing the most accurate results (Carballo-Pacheco & Strodel, 2017).

In the work done by Paul et al. (2022), 100 ns simulations of protein Bcl2 with a protein ligand BaxBH3 demonstrate how the modifications to the IDR sequence in turn tune ligand binding. Further, the effect of mutations in the IDR and at the protein active site can be captured using both MD simulations and/or molecular docking approaches (Jakhmola et al., 2022). An essential requirement for molecular docking is the accurate conformational sampling of the active site amino acid residues. This work and Paul et al. (2022), illustrate that MD simulations prove to be very useful in conformational sampling of proteins with IDR sequences.

Acknowledgements

We would like to thank the Marie Curie Actions program for funding the project, the Minerva supercomputer resource in Warwick University.

Disclosure Statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

Marie Skłodowska-Curie actions under the Seventh Framework Programme (FP7) for Research and Technological Development (EC Contract 316630).