The conventional degradation of halocarbons consists of two distinct steps: 

• Conversion of halocarbons to an intermediate via Cytochrome P450CAM under anaerobic conditions.
• Subsequent conversion of that intermediate to oxidised products via the action of Toluene Dioxygenase under aerobic conditions.

To further our understanding of microbial halocarbon degradation pathways, a kinetic enzyme model based on a system of ordinary differential equations was built. With the help of this model, the kinetics of the system were studied, the time scale of the enzymatic degradation steps was estimated, and the number of bacteria in the system was modelled as a function of time. A conceptual understanding of the kinetics of microbial halocarbon degradation played a critical role in experimental design to engineer the chassis (Pseudomonas putida), and to determine requirements for the design of the bioreactor.

INTRODUCTION

GENERAL SCHEMATIC OF ENZYME SYSTEM

A BASIS FOR OUR MODEL

SKELETON OF OUR MODEL

THE DIFFERENTIAL EQUATIONS

FINAL DIFFERENTIAL EQUATIONS

Results

Graphs

Please scroll the table if you cannot see the entire graph.

Oxygen concentration as a function of time The amount of oxygen in the system falls off approximately linearly with time. According to this, it can be expected that after 1hr of oxygen insertion, the oxygen should fall off to acceptable levels. This result helps in the switching between hypoxic and aerobic phases.	Substrate concentration as a function of time The concentration of the substrate decays exponentially which fits the results by Wackett et al.[1]
Concentration of Intermediate as a function of time The concentration of the Intermediate is predicted to increase exponentially time. Here, however, it is approximately linear	Concentration of Product as a function of time The concentration of the Product increases exponentially with time, as predicted.
Number of bacteria as a function of time The number of bacteria (N) increases logistically and then decays due to the product toxicity factor. Thus, a culture with an OD of 0.7 would be expected to stay alive for nearly 2 hours. This set the time scale for the assays. This graph is consistent with the experimental data which indicates that the bacteria die off in a span of 2hrs.	Compilation of all the graphs Shows the results of all the graphs in a collated format. The graph helps in inferring the relevent quantities simultaneously. The legend from the graph is given below.

Values of Parameters

MOLECULAR DOCKING AND DYNAMICAL SIMULATIONS

INTRODUCTION

Molecular Docking is a powerful tool used to study and make predictions regarding ligand-protein binding, the pose and site, and the relative strength of binding. A molecular docking software performs a search based on an algorithm that evaluates the system's energy due to electrostatic and Van der Waals interactions for different ligand conformations until it is minimised. Based on the predicted pose, an affinity score is calculated. Given a list of docking candidates, this score can be used to rank them. It is used extensively in drug design applications to screen potential drug candidates.

Due to the varying needs based on the variety of available ligands and receptors, several different software packages have been developed. Hence, a crucial part of the docking process is understanding the capabilities of different available software, their options, and deciding the protocol to use. 

Molecular Docking, however, has a broader scope than just predicting binding poses and scoring ligands; nowadays, many software packages are capable of predicting binding poses while accounting for side-chain flexibility, predicting active sites, and also determining binding and docking positions in the presence of water molecules. This simulates the natural interactions better and provides better predictive power. Some of the capabilities have been used to understand the enzymes' behaviour and generate qualitative and quantitative predictions crucial for this project. 

A thorough literature review was conducted to understand the features offered by the various docking software, including Autodock 4, Autodock Vina, GOLD, Flexx, and Glide, amongst others. Based on the current literature[15], the Autodock suite of tools was found to be the most suitable due to the high accuracy of pose prediction and better accessibility and reproducibility because it is freely available. Autodock is also one of the most highly cited docking software applications and is used by Professor Nagasuma Chandra's lab, which also helped with understanding the different applications and benefits of Autodock. 

The protocol for different experiments for the enzyme was decided based on several considerations with reference to published and cited protocols. The protocol broadly follows [15], with modifications made when necessary. The protocols relevant to the experiments conducted are given below.

Protocol

1. The enzyme database entries were found on the RSCB and Uniprot databases and cross-verified using BLAST search. Based on parameters such as resolution and presence of co-crystallised ligands, a suitable PDB file was selected with the best possible resolution and the minimal number of ligands co-crystallised (with necessary cofactors, redox factor, supporting ions and natural ligands). 
2. Suitable ligands were identified based on common Halocarbons used in the refrigerant industry, and availability for experimentation were selected. The SDF files were obtained from PubChem. The files were converted to PDB format using Pymol. In addition, Camphor and Toluene were also selected as a ligand for control and comparison. 
3. Ligand and Receptor coordinate preparation was carried out as listed in the protocol (2.1). Note that the Fe atom in Heme was assigned a Partial charge of +3 based on the literature on Cytochrome P450 enzymes and their mechanism. Similarly, in Toluene Dioxygenase, the Fe atoms (both in ferredoxin and ions) were assigned a charge of +2 and the S atoms a partial charge of -2. (refer to Protocol 2.1 Step 8)
4. Active sites were obtained from the literature [17] and were visualised in visualisation software. The list of possible active sites for Reductive Dehalogenation reactions in Cytochrome P450cam includes F87, Y96, F98, T101, T185, L244, V247, V295, D297, I395, and V396. 
5. Active sites were then predicted from AutoLigand (Protocol 2.3). Docking was carried out (Protocol 2.2) with Camphor and Toluene, and the docked and crystallographic positions were visualised and compared along with the active site and active site predictions. The comparison is used to verify the protocol.
6. Following verification, docking was carried out for all the other ligands of interest. (Note: Hydrated dockings were avoided as the protein and ligand interact via hydrophobic interactions)
7. All the obtained results were visualised on a suitable visualisation software the results were tabulated. 

Results:

Docking with Cytochrome P450cam (Camphor-5-monooxygenase)

The results obtained from docking the various ligands are presented below in a tabular format. Docking with camphor is also shown as a verification of the protocol.

Image of the enzyme Cytochrome P450cam →	Video Displaying Docking with Camphor and Active Site Pocket
Affinity of Camphor = -7.0 kcal/mol

Table with Docking Results

Ligand	Classification	Structure	Docking Position	Afffinity	Predicted Product
R-22	HCFC			-3.2	\( \ce{CO2} \) through \( \ce{COF2} \) intermediate
R-12	CFC			-3.7	Unknown
R-123	CFC			-4.7
R-110	Halocarbon			-5.0
R-112a	CFC			-5.2
R-113	CFC			-5.3
R-141b	HCFC			-3.8	Acetic acid (under aerobic conditions)
R-134a	HFC			-4.0	Trifluoroacetic acid (under aerobic conditions), and 1,1-difluoroethane (under anaerobic conditions)
R-124	HCFC			-4.7
R-142b	HCFC			-3.7	Unknown
TCE	Halocarbon			-3.8	Trichloroethanol and Trichloroacetic acid

Observations:  We observe that CFC's and CFC like halocarbons have highest negative affinity. Meanwhile HCFC's have lower affinity and HFC's are even lower. The lower affinity for HFCs and Methane based halocarbons may explain their reluctance to undergo the enzymatic reaction.

Docking with Toluene Dioxygenase (Benzene 1,2-dioxygenase)

Image of the Enzyme Toluene Dioxygenase	Image displaying Docking with Toluene and Active Site Pocket
Affinity of Toluene = -5.7 kcal/mol

Table with Docking Results

Initial Substrate	Ligand Structure	Docking Position	Afffinity	Predicted Product
R-110			-4.5	Glyoxylic acid/Oxalic Acid
R-112a			-4.6	Glyoxylic acid/Oxalic Acid
R-113			-4.6	Glyoxylic acid/Oxalic Acid
R-123			-4.0	Glyoxylic acid/Oxalic Acid
R-124			-4.0	Glyoxylic acid/Oxalic Acid
R-134a			-3.5	Glycolic acid
TCE			-3.9	Glyoxylic acid/Oxalic Acid

Observations  We observe that all the end-products are as expected and can be utilised in metabolism of the bacteria.

Inferences:

The docking images show that the predicted active site matches the active site pocket in the literature. In addition, the binding site of camphor and our halocarbon substrates is the same. An essential aspect of this binding site is allowing interaction with the heme group. The reversible reduction and oxidation of the iron atom in the Heme from Fe(III) to Fe(II) is essential for the proposed mechanism, as elaborated on the design page. 

Similarly, the dockings and active site for Toluene Dioxygenase are closely matched with that for the native ligand, Toluene. In addition, the docked positions are in close proximity with the \( \ce{Fe^2+} \) ion, which reversibly oxidises to \( \ce{Fe^3+} \) thereby enabling the reactions to occur as per the proposed mechanisms.  Hence, the docking experiments were essential in providing critical insights into our mechanism and helped provide credibility and better predictions for the reaction products.  

The affinity results also help us interpret the feasibility of our proposed reactions which are essential for proof of concept. Literature provides evidence for reductive dehalogenation reactions for substrates like Hexachloroethane and a few CFCs. The affinities are reflective of the binding energies of the ligand and protein in the active site. 

More negative values indicate a greater strength of binding; hence, the affinities are indicative of the feasibility of the reaction. R-110, i.e., Hexachloroethane, has a binding affinity of -5 kcal/mol, while some of the other refrigerants we are interested in degrading (CFCs and HCFCs) have a range from 4 kcal/mol to 5kcal/mol. This indicates that our proposed reactions are suitable from an energetic basis. 

Similarly, for Toluene Dioxygenase, TCE (trichloroethylene) has been shown to degrade experimentally. The binding affinity for the same was found to be -3.9 kcal/mol from docking experiments. Most of the other substrates for the enzyme (after Cytochrome P450cam) have similar, if not higher binding affinities indicating a high likelihood of degradation by the enzyme. 

The docking experiments proved to be an engineering success by taking inputs from the WetLab team and providing critical information for the proposed reactions and substrates. 

 Taking our Docking Results Further

One of the challenges faced in enzyme kinetic modelling was the unavailability of specific parameters since the reactions are not well studied. One critical parameter was the Dissociation constant, \( k_d \), which describes how well the ligand binds to the enzyme. An important result from statistical modelling of the ligand-receptor system gives us the equation:

Now, from the Molecular Docking results, we obtain the Binding energy \( E_{binding} \), and hence we require the solvation energy in order to estimate the Dissociation constant. For obtaining the value of \( E_{solvation} \), we have modelled the dynamics of our substrates in water using Molecular Dynamics simulations. 

A Physical Interlude[18]

Here a brief derivation of Equations 1 and 2 will be shown. Professor Sumantra Sarkar, Department of Physics, IISc, has been instrumental in guiding us through modelling the Ligand-Receptor model. 

Consider the system with a fixed temperature \( T \) and let the energy fluctuate about an average energy \( <E> \). Under these assumptions, we can apply the Boltzmann Distribution for probabilities of a microstate with energy \( E_i \) as

Considering our system of multiple ligands and a receptor, we have two possible states: ligand-bound and ligand-unbound. Consider the system's energy in the bound state to be \( E_1 \) and in the unbound state to be \( E_2 \). Hence we have

Now consider the system where there are L ligands and N possible states a ligand can occupy in the solution. In the unbound case, all the ligands are in the solution, and that gives a multiplicity of \( m_{ub} = \frac{N!}{(N-L)!L!} \), and the energy is \( LE_{sol} \). When one of the ligands is bound to the receptor, the number of ligands in the solution becomes \( N-1 \). This gives a multiplicity of of \( m_b = \frac{N!}{(N-L+1)!(L-1)!} \), and the energy is \( (L-1)E_{sol} + E_{bind} \). Applying the Boltzmann formulation along with the multiplicities gives us:

As \( N>>L, N-L+1 \approx N \). Hence, this gives us

Now, assume that the ligand has some volume \( V_0 \) and the total volume of the container here is \( V \). Then, \( \frac{L}{N} = \frac{L}{NV_0 \times c_0} \) where \( c_0 = 1/V_0 \) is the concentration of 1 molecule per ligand volume. Also, \( L/NV_0 = L/V = c \) where \( c \) is the concentration of the enzyme. Hence, this would give us \( \frac{L}{N} = \frac{c}{c_0} \) and

Hence, we have derived the equations we are going to use here. 

Molecular Dynamics Simulation

The solvation energies for our molecules of interest were unavailable on existing databases; hence, prediction by simulation was necessary. Out of several methods employed in the literature, Molecular Dynamics simulations are commonly used with accurate predictions. 

To perform the simulations, we used the package GROMACS[19], along with the tutorials, which were helpful in performing the simulations. It uses Free Energy Perturbations (or alchemical free energy computation)[20] to evaluate the difference in free energy between two states. In order to perform the simulations, the Jupyter notebook was followed with the GROMACS configuration file and the topology file obtained from the Virtual Chemistry database[21].  The topology files used Generalized Amber Force Field (GAFF) and the Optimised Potential for Liquid Simulations for All Atoms (OPLS/AA). GAFF, OPLS/AA along with SPCE(a water model) yield good results as is well characterised in literature.

In order to perform successful molecular dynamics simulations in water, the topol.top files need to be edited. The files available from online sources or those obtained by using 'pdb2pqr'. In particular, the following changes are neccessary for the topol files to be compatible with solvation experiments: necessary:

• Inclusion of the following code block at the top of the file 

  ;#define _FF_OPLS
                  #define _FF_OPLSAA
                    #include "oplsaa.ff/forcefield.itp"

• Commenting out the [ defaults ] section 
• Inclusion of the following code block at the end of the document above the [ system ] section. 

•                 ; Include SPC water topology
                    #include "oplsaa.ff/spc.itp"
                  

In addition, file name changes must be made in the Jupyter notebook depending on the naming of the files. The topology file should be saved as topol.top. 

The molecular dynamics simulation was run for Hexachloroethane and Camphor. Running the computation for all the ligands was prohibitive because of the high runtime and computational power required. Based on the simulations and docking with Cytochrome P450cam, the results obtained are tabulated below:

Our predictions based on molecular simulation methods give us a close match, hence verifying the capability of the methods and protocols. The simulations require intensive computational power. For better predictions, the program can be run with higher number of steps which would give more consistent results. Experiments run may yield slightly different results due to the fact that the calculations are perturbative and the approximations of statistical mechanics may not hold in every case. Regardless, the method helps in contributing not only to WetLab, but can also help in providing better constants needed for future implementations.

Active Side modification

The active site of the enzyme is the region where the substrate binds and undergoes a chemical reaction. The active site is usually a pocket in the enzyme, which is lined with amino acids that are involved in the chemical reaction. The residues around the active site define the specificity and the reactivity of the enzyme. Altering the residues in this region directly influences the interactions made with the substarte which can alter the specificty of the enzyme. The shape and size of the pocket determine what substrate the enzyme can accomodate and also determine the orientation in which the substrate binds. This in turn influences contacts with Catalytic residues, which may either increase or decrease the efficiency with which the reaction catalysis occurs.

We modelled active site modifications in-silico in order to analyse if engineering the protein can help facilitate the reactions required. We performed this analysis for Cytochome P450cam, which is involved in the first catalytic step. The active sites of the enzyme were reviewed and existing literature on active site modifications was studied. Based on that, it was shown that modification of the residues 87 (F -> W), 247 (V -> L) and 96 (Y -> F) can help in increasing the efficiency of the reaction [17]. Experimental evidence showed that F87W-V247L was the most active mutant, however, it was specualted that the Y96F mutation may also help in increasing the efficiency of the reaction further. We have carried out the three amino acid modifications and tested the change in binding affinity via Molecular Docking.

In order to carry out the modification, the sequences were edited and the structure prediction was carried out by AlphaFold2[26] and the predicted structure was used for docking of Camphor, Hexachloroethane (R-110), and HFC-134a as per the protocol described above. The results are tabulated below along with images of the enzyme and docked ligand.

Image of Overlap of Modified Protein (green) with Original Protein (pink)	Image of Overlap of Modified Protein Active Site (green) with Original Protein Active Site(pink) with labelled amino acid modifications

AlphaFold Predictions Validations

Observations: From the above table, we can observe that active site modification did indeed reduce the dissociation constant for Hexachloroethane as predicted in [17] and the activity is about 2.5 times higher as predicted. However, the same increase in Affinity is not observed for HFC-134a. With this small change, the ratio of the Dissociation constant after modification is \( e^{-\beta/10} = 1.06\) which is not a significant change.

From the observations, we can infer that the active side modification does help in the case of hexachloroethane and possibly other halocarbons, but it did not appear to work for HFC-134a, a very popular refrigerant. Hence, while active site modifications may help in increasing reactivity, properties like activation energy reduction and incorporation of stronger redox partners like Spinach ferredoxin [27]

Bioreactor and Financial Modelling

REFERENCES