Team:Manchester/fabAmodeltest

Figure 1. Two colour image of βHACdH homodimer with labelled N- and C-Termini (arrows)

Introduction
Our model of the fatty acid biosynthesis pathways highlighted several key enzymes that could be altered to produce palm oil, including Delta 9, Delta 12 and β-hydroxy acyl-ACP dehydrase (βHACdH), which agrees with previous publications on these enzymes. Delta 9 and Delta 12 are involved in linear pathways, meaning that there are obvious reactants and products, so the overexpression of them can be measured directly via LC-MS and other techniques. However, βHACdH which is expressed from fabA is a β-Hydroxydecanoyl Thiol Ester Dehydrase involved in a cyclic pathway^[1] specifically the conversion of β-hydroxy acyl-ACP to Enol acyl-ACP as part of the fatty acid biosynthetic pathway^[2]. Therefore characterising any specific change due to βHACdH overexpression will be challenging, as any products will automatically be involved in the proceeding stage of the cycle and it would therefore be difficult to determine the effect of overexpression. An alternative to measure the overexpression of βHACdH, would be through the addition of His-tags to either the N-terminal and or the C-terminal of βHACdH. However, depending on the structure of βHACdH, the addition of His-tags could potentially interfere with expression, protein folding, enzymatic functions and interactions^{[3][4][5][6][7]}. These problems are particularly applicable to βHACdH, as it forms a homodimer, as shown by Leesong et al., 1996^[1]. Therefore, the addition of His-tag could potentially interfere with the interaction domain and thus the formation of a homodimer, which would be consistent with several reports^[4][5][6]. To address this issue, we decided to perform a molecular dynamics simulation using the GROMACS software package^[8] on a structure of βHACdH, determined by X-ray crystallography by Leesong et al., 1996^[1]. This would allow for the trajectories of the N- and C-Termini of βHACdH over the course of the simulation to be studied, therefore allowing us to identify which terminal would be more suited for His-tag addition.

Installing GROMACS
The Groningen Machine for Chemical Simulation (GROMACS) Version 4.5.5 [8], was installed on a MacBook Pro 2011 model, operating OS X 10.8.3, with a 4 GB 1333 MHz DDR3 memory, 2.3 GHz Intel Core i5 processor and a 320 GB SATA disk drive. Prior to installing GROMACS, “command line tools” were installed within Xcode, Version 4.6.1. GROMACS was installed to single precision, with the source file downloaded from www.gromacs.org/Downloads. In addition to GROMACS, the FFTW library, Version 3.3.3^[9] was installed, with the source file downloaded from www.fftw.org/download.html.

Figure 2. GROMACS Molecular Dynamics Simulation Workflow. The commands used at each step are included, with a brief description of their function

Behind the scenes of the Molecular Dynamics Simulations
For the simulations of βHACdH, a dimeric structure of βHACdH (PDB ID: 1MKB) derived by X-ray Diffraction, with a resolution of 2 Å from an E.Coli expression system by^[1] was used. The molecular dynamics simulation was performed within the GROMACS package using version 4.5.5 [8], with the AMBER99SB force field^[10], the transferrable intermolecular potential 3P (TIP3P) water model^[11] and the original crystal waters, from the X-Ray Crystallography with periodic boundary conditions. The methods for generating both topologies and parameters for G16bP are as described above. For all cases of the Protein in water, a protocol derived from the “Lysozyme in Water” GROMACS tutorial (by J. Lemkul, Department of Biochemistry, Virginia Tech) and optimized for βHACdH was used. Briefly, a cubic cell of 2 nm in diameter with the protein centered was used and filled with the generic single point charge 216 (SPC216) water configuration^[12]. The system was neutralized to a salt concentration of 150 mM by adding Na+ and Cl-. Energy Minimization (EM) was performed using the Steepest Descent Algorithm^[13] with a tolerance of 1000 KJ mol-1 nm-1, to remove any steric clashes or inappropriate geometries. In the simulation, the long-range electrostatic interactions were modeled using the Particle-Mesh- Eswald (PME) method ^[14][15] and the Linear Constraint Solver (LINCS) algorithm^[16] to preserve chemical bond lengths. Temperature and pressure coupling were performed independently, using a modified Berendsen thermostat^[17] at constant temperature of 300 K at a time constant of 0.1 ps and the Parinello-Rahman barostat algorithm^[18] at a constant pressure of 1 bar, for a time constant of 2 ps and V-rescale, respectively. The Molecular Dynamics simulations trajectories were analysed with the GROMACS analysis tools^[19], to output structural molecular dynamics trajectories and PyMOL (The PyMOL Molecular Graphics System, Education-Use-Only, Version 1.3 Schrödinger, LLC) used to visualize and create both still structures and videos. All calculations of progression plots from the simulations were produced using the GROMACS analysis tools with output files viewed in the Microsoft Excel 2011.

Getting a working GROMACS simulation for βHACdH
To study the motions of the N- and C-Termini of βHACdH with molecular dynamics, we first had to get an optimized system preparation protocol, which we based on the “Lysozyme in Water” GROMACS tutorial by ^[20]. Firstly, we built a simulation box of 2nm around βHACdH, which was sufficient to satisfy the minimum image convention and the simulation cut-off schemes without adding excess solvent. βHACdH had to then be solvated within this box by a solvent configuration compatible with the solvent model applied to the protein, in this case the SPC216^[12], which is compatible with the TIP3 water model^[11]. The system is then neutralized through the addition of ions to a molarity of 150 mM, therefore allowing the system to reach a neutral state. Our results indicate that βHACdH is well solvated in SPC216 water and neutralized to a molarity of 150 mM in a cubic 2 nm cell.

Figure 3. Establishing, solvating and neutralising βHACdH, for simulation preparation. A. The establishment of a 2nm simulation cube around the dimerised structure of βHACdH. B. The solvation of the βHACdH structure by the addition of SPC216 water molecule as a solvent. C. Neutralisation of the solvated system by the addition of Na+ and Cl- ions. Na+ and Cl- ions are represented by Blue and Green spheres, respectively and SPC216 water is represented by Cyan coloured molecules.

Once we had prepared the system, it was relaxed by energy minimization. This ensures that there are no steric clashes or inappropriate geometries that exist, by applying the steepest descent algorithm^[13]. With the minimized structure, the system of solvent and ions around the protein are equilibrated to become orientated about the protein solute at the same temperature, by an isothermal-isochoric ensemble^[17]. Pressure is then applied using an isobaric-isochoric ensemble, thereby ensuring that the system reaches a proper density^[18]. Our energy minimized structure shows a decrease from -6.75E+05, to a maximum energy plateau for the system of -1.05 E+06 after 1238 ps of minimization time. The temperature of the system quickly reaches the target value of 300 K remaining stable, with an average temperature of 299.82 K, the equivalent of 27 °C over the 100 ps equilibration. Over the course of the 100 ps equilibration stage, both the pressure and density of the system averages 0 bar and 1015.19 Kg m-3. This is close to the experimental value of 0 bar and 1000 Kg m-3 and the equivalent of Earth’s atmospheric pressure at sea level and the density of water. The pressure fluctuations are consistent with the applied isothermal-isobaric ensemble and is suggestive of compatible molecular dynamics conditions for simulations of βHACdH.

Figure 4. Graphs of the Energy Minimisation, Temperature, Pressure and Density equilibration of the βHACdH simulation system prior to simulation.

To His-tag or not to His-tag
Now that βHACdH is ready to undergo simulations, we ran a simulation for 1 ns under the notion that we would be able to visualize motions around the N- and C-Terminals during the course of the simulation and therefore determine, which terminal would be more appropriate to add His-tags to. The conclusion for our simulation was that the N-Terminal is ideal for the addition of His-tags for several reasons. Firstly, the C-terminal is localized in close vicinity to the interaction domain of the βHACdH homodimer, therefore the addition of His-tags could possibly interfere with the dimer interaction^[6]. Our model also shows that the C-terminal is more dynamic compared to the N-terminal and there are several times in the simulation that the C-terminal interacts with the dimerization domain and may interfere with the folding and function of the protein^[4][5][6]. Therefore we concluded that the N-terminal would be ideal to add the His-tags to, as the N-terminal is less dynamic and will be less likely to interfere with folding and the protein function. With this in mind, our experimental team began to design the His-tagged FabA BioBrick (BBa_K1027003) to express a N-Terminal His-tagged βHACdH to use in the characterisation of βHACdH overexpression.

Figure 5. Overlay of structures from 1 ns βHACdH simulation. Images of overlaid from the following respective time points: 0 ps, 250 ps, 500 ps, 750 ps and 1000 ps with the following colours indicating each individual image: Green, Blue, Purple, Orange and Grey, respectively. Both the N-Terminal and C-Terminal, are specified (Dotted Box), with a zoom in on each respective terminal at an angle appropriate to visualise the positions of the terminals.

References
[1]. Leesong, M., Henderson, B.S., Gillig, J.R., Schwad, J.M. and Smith, J.L. 1996. Structure of a dehydratase-isomerase from the bacterial bathways for biosynthesis of unsaturated fatty acids> two catalytic activities in one active site. Structure 4:253-264.
[2]. Xu, P., Gu, G., Wang, L., Bower, A.G., Collins, C.H. and Koffas, M.A. 2013. Modular optimization of multi-gene pathways for fatty acids production in E.coli. Nature Communications. 4: (1409): 1-8. DOI: 10.1038/ncomms2425
[3]. Carson, M., Johnson, D.H., McDonald, H., Brouillette C. and Delucas, L.J. (2007) His-tag impact on structure. Acta Crystallogr D Biol Crystallogr. 63(3):295-301.
[4]. Chant, A., Kraemer-Pecore, C.M., Watkin, R. and Kneale G.G. 2005. Attachment of a histidine tag to the minimal zinc finger protein of the Aspergillus nidulans gene regulatory protein AreA causes a conformational change at the DNA-binding site. Protein Expr Purif. 39(2):152-9.
[5]. Freydank, A.C., Brandt, W. and Dräger, B. 2008. Protein structure modeling indicates hexahistidine-tag interference with enzyme activity. Proteins. 72(1):173-83. doi: 10.1002/prot.21905.
[6]. Klose, J., Wendt, N., Kubald, S., Krause, E., Fechner, K., Beyermann, M., Bienert, M., Rudolph, R. and Rothemund, S. 2004. Hexa-histidin tag position influences disulfide structure but not binding behavior of in vitro folded N-terminal domain of rat corticotropin-releasing factor receptor type 2a. Protein Sci. 13(9):2470-5.
[7]. Woestenenk, E.A., Hammarström M., Van Den Berg, S., Härd, T. and Berglund, H. 2004 His tag effect on solubility of human proteins produced in Escherichia coli: a comparison between four expression vectors. J Struct Funct Genomics.
5(3):217-29. [8]. Pronk, S., Páll, S., Schulz, R., Larsson, P., Bjelkmar, P., Apostolov, R., Shirts, M.R., Smith, J.C., Kasson, P.M., Van der Spoel, D., Hess, B., Lindahl, E., 2013. GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics: 1–10.
[9]. Frigo, M., Johnson, S.G., 2005. The Design and Implementation of FFTW3. Proceedings of the IEEE 93: 216–231.
[10]. Hornak, V., Abel, R., Okur, A., 2006. Comparison of multiple Amber force fields and development of improved protein backbone parameters. Proteins: Structure, Function and Bioinformatics 65: 712–725.
[11]. Jorgensen, W.L., Chandrasekhar, J., Madura, J.D., Impey, R.W., Klein, M.L., 1983. Comparison of simple potential functions for simulating liquid water. The Journal of Chemical Physics 79: 926.
[12]. Berendsen, H., 1981. Interaction models for water in relation to protein hydration. Intermolecular Forces: 331–338
[13]. Kiwiel, K., Murty, K., 1996. Convergence of the steepest descent method for minimizing quasiconvex functions. Journal of Optimization Theory and Applications 89: 221–226.
[14]. Darden, T., York, D., Pedersen, L., 1993. Particle mesh Ewald: An N⋅log(N) method for Ewald sums in large systems. The Journal of Chemical Physics 98: 10089.
[15]. Karttunen, M., Rottler, J., Vattulainen, I., Sagui, C., 2008. CHAPTER 2 Electrostatics in Biomolecular Simulations  : Where Are We Now and Where Are We Heading  ? Current Topics in Membranes 60: 49–89.
[16]. Hess, B., Bekker, H., Berendsen, H.J.C., Fraaije, J.G.E.M., 1997. LINCS: A linear constraint solver for molecular simulations. Journal of Computational Chemistry 18: 1463–1472.
[17]. Bussi, G., Donadio, D., Parrinello, M., 2007. Canonical sampling through velocity rescaling. The Journal of Chemical Physics 126: 014101.
[18]. Parrinello, M., 1981. Polymorphic transitions in single crystals: A new molecular dynamics method. Journal of Applied Physics 52: 7182.
[19]. Lindahl, E., Hess, B., Spoel, D. Van Der, 2001. GROMACS 3.0: a package for molecular simulation and trajectory analysis. Molecular Modeling Annual: 306–317.
[20] J. Lemkul, Department of Biochemistry, Virginia Tech. http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/lysozyme/. [Last Accessed: 27/09/2013]