# Team:UGM Indonesia/Engineering

<!DOCTYPE html> Engineering

# Engineering

Our main goal is to create a feasible on-off system of the cyanide production and degradation according to either the presence or absence of arabinose. This page displays three engineering cycles; i.e., cyanide production, cyanide degradation, and arabinose conversion. Afterwards, an engineering cycle about the Auviola on-off system is performed to associate the aforementioned three cycles. An engineering cycle about the bioreactor is also demonstrated as the implementation of our Auviola system.

#### HCN-regulating System

###### What are the properties of HCN production in C. violaceum?

The HCN synthase, or also known as glycine dehydrogenase, is an enzyme that facilitates the conversion of an amino acid glycine into cyanide.1 This enzyme is a heterotrimeric protein as encoded by a cluster of three genes called hcnABC operon. This enzyme belongs to the oxidoreductase class as it oxidizes the amine (CH-NH2) functional group into imine (C=NH) and consecutively cleaves the molecule into HCN and carbon dioxide (CO2).2

###### What are some concerns to be considered?

The resulted cyanide species were designed to dissolve the gold into the soluble coordination compound of (NaAu(CN)2) through the reaction below:

4Au + 8NaCN + O2 + 2H2O 4(NaAu(CN)2) + 4NaOH1

According to previous study, the gold leaching process requires HCN compounds with concentration of 600 mg/L.2 On the other hand, a wild-type of C. violaceum was previously studied to produce 0.71 - 16 micromoles/mL which nearly equals to 19 - 432 of HCN in 8 hours.3 Therefore, an optimization of produced HCN concentration should be conducted to achieve an effective gold leaching process.

###### What are the types of regulators tested in C. violaceum?

Two of several regulator genes have tested in C. violaceum are araC and tetR.4 The araC works dependably to arabinose level which acts as both an activator in the presence of arabinose and a repressor in the absence of arabinose.5 A PBAD is one of the inducible promoters presented in C. violaceum and regulated by the araC.6 On the other hand, tetR acts as a repressor and normally play a role in tetracycline resistance with PTET promoter.7

,8
###### Wetlab work: HCN Production

Our strategies for improving cyanide metabolism in the cells are to increase the accumulation of HCN, hydrolysis of HCN, and L-arabinose converter. In this study, the proposed design for increasing HCN production has been achieved. The principle behind this system is overexpression of hcnABC gene encoding cyanide synthase, an enzyme that is responsible for HCN synthesis, by addition of L-arabinose-induced pBAD-AraC regulon. The genetic circuit for this study is written in the Design page. At first, we were focusing on the first plasmid construction for the HCN production. To achieve this goal, several standard Biobrick parts as well as the hcnABC operon was required. Some standard parts that we use are listed below:

2. BBa_B0015 (double terminator)
3. BBa_B0034m (strong RBS)

The hcnABC operon as the main insert, were obtained through gene synthesis from Twist Bioscience. Subsequently, the plasmid map was created to determine the overall gene construction into the pSB1C3 vector. Benchling application system was used to design the circuit, particularly the HCN production plasmid (Kindly visit our Results page for our plasmid construct!).

We used AraC-pBAD derived from BBa_K1602055 to regulate the hcnABC operon. We modified the hcnABC operon by adding RBS before each gene (hcnA, hcnB, hcnC) and also double terminator B0015. The hcnB and hcnC genes have a high percentage of GC content so we faced some problems to synthesize it. Thus, we altered the hcnB gene sequence to reduce the GC content and to remove some illegal sites. We also altered some bases in hcnC gene to remove the illegal site. For a detailed description of the construction, please visit our Design page. The recombinant plasmid transformed was verified through colony PCR and the inserted parts were confirmed by sequencing analysis.

A Kinetic modelling for demonstrating an overexpression of cyanide was performed in Tellurium. This modelling used several parameters obtained from literature. In addition, genome-scale modelling was performed to optimize the metabolic flux of C. violaceum to achieve optimal cyanide.

###### Drylab work: Kinetic Modelling

Beside the wet lab, we also perform dry lab analysis to model the HCN production, HCN degradation, and L-Arabinose isomerase in engineered C. violaceum. We use kinetic modelling approach to perform this analysis.

Kinetic modelling of metabolic pathways is important for predicting the output of a metabolic process, and to understand what type of modification should be carried out if the output is not desirable. For our model, the complete equations are explained in the Model page. We collected most of our parameters from research papers and past iGEM team modelling.

###### Wetlab work: HCN Production

After we went through the design phase, we started to conduct our project at the Biotechnology Laboratory of Inter-University Center in Universitas Gadjah Mada. We ordered hcnA, hcnB, and hcnC sequences from Twist Bioscience. The sequences were referred to the genome of Chromobacterium violaceum ATCC 12472 (NCBI: GenBank: AE016825.1).

To meet the BioBrick standard, we altered some of the sequences i.e., hcnB and hcnC, to reduce the GC content and remove the illegal sites. However, Twist Bioscience reported a problem in synthesizing the hcnC gene, thus the hcnC synthesize failed. Thus, we isolated the hcnC gene using PCR from the genome of C. violaceum ATCC 12472 and removed an illegal site inside the sequence using overlap PCR (Phusion polymerase) (see Results page table 1, no. 7,8, and 9). The assembly of plasmid backbone (pSB1C3), AraC-pBAD (BBa_ K1602055), RBS (BBa_B0034m), hcnABC operon, and double terminator (BBa_B0015) were done using Gibson Assembly. The assembled constructs were then transformed into two strains of E. coli, BL21(DE3) and DH5α using standard heat shock protocol. The screening of transformats were selected on LB agar supplemented with Chloramphenicol for both strains.

###### Drylab work: Kinetic Modelling

The equations needed for modelling gold cyanidation are explained in our Model page. We were utilizing python in Jupyter with the Tellurium package for most of our model, except for Au and CN interaction.

Some of the kinetic constant cannot be found on any research paper or open source database that is available on the internet. Thus, we used the value of similar parameter(s).

###### Wetlab work: HCN Production

The resulting gene constructs were transformed into two chassisses and verified by colony PCR using VF2 VR primers. After the colony PCR and electrophoresis, the expected bands were incorrect in size and unobtainable. Our amplified insert using VF2 VR primers should be around 4600 bp, yet we merely obtain approximately 3000 bp band. On the upper side of some bands, ColA. There were thin bands of which the sizes are over 3000 bp, which indicated that the desirable construct was inserted to the chassis. Therefore, we subcultured the ColA colony to LB broth media supplemented with Chloramphenicol overnight for plasmid isolation (check the Results page for plasmid isolation result). For further confirmation, we send our isolated plasmid for sequencing. Now we are still waiting for the result of sequencing and we will update and report the result by the time we receive it.

###### Drylab work: Kinetic Modelling

The graph below shows concentration (M) vs time (second).

###### Wetlab work: HCN Production

It was such a challenge for us to construct the whole operon with the size over 4000 bp with high GC content. A more robust PCR method such as an adjustable KOD with high-fidelity PCR master mix kit is necessary to increase the PCR product and avoid point mutation that might happen in the amplified product.

###### Drylab work: Kinetic Modelling

We can learn that rhodanese was produced instantly, then the concentration of HCN decreased rapidly. The maximum speed of gold cyanidation only lasts briefly, less than 5 minutes. Biohydrometallurgy method is favored in mineral processing field due to their ability to extract low-concentrated ore. With this kind of trend, it is almost impossible to gain a feasible gold extraction yield. In lab experiments, most bioleaching processes managed to attain large mineral yield in laboratory experiments. This modelling might be not reflecting the real conditions of the bioleaching process.

###### Drylab work: Kinetic Modelling

We need to reconsider some of the parameters, especially those based on assumption. Rate of large organic molecule degradation is one of the main examples. There are 2 large organic molecules that are involved in this modelling, protein (enzyme) and monosaccharide. We need to separate the decay rate of protein and monosaccharide. SCN formation constant also needs to be changed. Due to larger scale use and relatively dynamic environment, we assumed that the collision rate would be significantly slower. Therefore, the constant also becomes lower. Rhodanese activates almost simultaneously after HCN hits its highest concentration, hence we need to set a starting point to delay rhodanese synthesis. The starting point will be based on L-arabinose concentration.

###### Drylab work: Kinetic Modelling

The parameter that we consider to be revised are:

ParameterDescriptionValueUnitsSource
\gamma_A1Rate of L-arabinose degradation3.46e x 10-5s-1[16]
pRhodanese synthesis starting point50% of initial L-arabinose concentration
kcat2Rate of thiocyanate (SCN) biosynthesis5.5 s-1
###### Drylab work: Kinetic Modelling

The results of modelling using those parameters are:

###### Drylab work: Kinetic Modelling

The HCN concentration reduction becomes less steep, thus gold cyanidation could be more efficient. It will prolong the almost optimum gold leaching period and surely increase the potential gold yield. Although we achieved desirable results through modelling, further laboratory experiments also needed to verify our model and obtain more realistic results. Large scale and pilot testing are also required for industrial application.

#### Genome Scale Metabolic Model: Gene knockout for enhancing the flux of cyanide production pathway

###### How to enhance cyanide production?

Our idea was to perform the Genome Scale Model (GEMs) for a knockout prediction and flux analysis to optimize the cyanide production. The output for this model were production envelope which compares the metabolic map between the knockout mutant and wild type as well as to illustrate the flux uplift to the cyanide production pathway. We use BiGG Model database (http://bigg.ucsd.edu ),17KEGG metabolic pathway (https://www.genome.jp/kegg/pathway/),18and Escher (https://escher.github.io/#/),19to confirm and simulate the pathway.

To perform the knockout prediction model, we use stochastic and deterministic approaches. Data resource that we used was the SBML model with XML language and Escher map which used JSON language. We carried out the engineering, simulation, and model alteration process using Jupyter Notebook and Jupyter widget which were provided in the Escher module. In general, the GEMs implementation used a particular module in Python which related to its task such as Escher, Cameo, and Cobra. These modules have a good documentation which enables us to understand its usage in either class or the methods.

Tay et al., (2013) once performed a comparative proteome analysis through cyanogenesis in C. violaceum and established the overview of cyanogenesis pathway.20On the other hand, Benarjee et al. (2019) also did the GEMs to C. violaceum, yet merely to describe and identify NAD+ recycling.21In our project, we updated the model and metabolic map from Benarjee et al. (2018) with the adaptation of cyanogenesis pathway that have been published by Tay et al. (2013).20

,21

After updating the model and metabolic map, we try to knockout the prediction gene using OptGene and OPTKnock.22

,23

These approaches were recommended to us by our instructor, Mr. Matin Nuhamunada. The results of this model were expected to illustrate the metabolic flux shiftment to cyanide production through gene knockout process.

We used a model and map file from Banerjee et al. (2019) with some modification. Curation is still needed to match the reactions and metabolites IDs between the map and the model. We found that there are 32 reactions with some IDs mismatch between the Escher map and the model, three ‘double’ reactions, and two metabolites which should be curated. After that, we started to fix the IDs in the model. Following that, we had our discussion with our instructor, Mr. Matin Nuhamunada who helped us to evaluate and fix the model. In addition to the model, we also modified the metabolic map. We manually added some reactions in the cyanide production pathway to the metabolic map using the Escher module, because at the beginning it still lacks the cyanide production pathway. For further detail, please visit Model: Genome Scale Metabolic Model .

We performed a phenotypic phase plane to show us a trade-off between our target (in this case cyanide production) and the biomass of C. violaceum. We also implemented the OptGene and OptKnock approaches from Cameo to predict the gene knockouts for enhancing the cyanide production. Based on the OptGene and OptKnock gene knockout result, we then analysed the reaction flux of the metabolic maps (Escher map) between the wild type and knockout mutant of C. violaceum, especially in cyanide production pathway. For further detail of our test result, please visit Model: Genome Scale Metabolic Model .

We found some difficulties due to the limited solution from OptGene and OptKnock. Some examples from the knockout prediction, especially which mentioned in https://cameo.bio/05-predict-gene-knockout-strategies.htmlwere merely discussed and model the exchange reaction as the target,22

,23

in fact we did not use exchange reaction in our model. Thus, we did the consultation with our instructor and they suggested us to use the exchange reaction because it can represent the knockout prediction for the accumulation of targeted extracellular metabolic. However, in our knockout model we solely try the cyanide production reaction (glycine:acceptor oxidoreductase) as the main target. Therefore, it is necessary to improve and add the cyanide exchange reaction to our model and transform the target in the knockout prediction into the cyanide exchange reaction.

#### Pyrite Dissolution Bioreactor

There are many possibilities of pyrite dissolution enzyme inhibition mechanism. The bioreactor main variables, such as µmax, KS, and KP are determined by fitting the experimental data with semi-batch bioreactor equation and model. This variable will also vary depending on: bioreactor operating condition, such as pH and temperature; substrate, such as different type of gold ore; and computational method, such as Lineweaver-Burk or linear regression.

Determining the inhibition mechanism of pyrite dissolution bioreactor will be useful for scale-up, developing process control (like PID), and further research for bioleaching of another metals using archaebacterial.

In this case, we modelled and calculated the results using python via jupyter notebook.

For bioreactor modelling, we wrote ODE based on the relation of X (cell count), P (product concentration), and S (substrate concentration) with V (total volume) by applying Monod model and dilution rate:

\begin{align*} \frac{d(XV)}{dt} = V\frac{dX}{dt} + X\frac{dV}{dt} = V\frac{dX}{dt} + F(t)X \\ \frac{dX}{dt} = - \frac{F(t)}{V}X + r_c(X,S) \\ \end{align*} \begin{align} \frac{d(PV)}{dt} = V\frac{dP}{dt} + P\frac{dV}{dt} = V\frac{dP}{dt} + F(t)P \\ \frac{dP}{dt} = - \frac{F(t)}{V}P + r_P(X,S) \\ \end{align} \begin{align*} \frac{d(SV)}{dt} = V\frac{dS}{dt} + S\frac{dV}{dt} = V\frac{dS}{dt} + F(t)S \\ \frac{dS}{dt} = \frac{F(t)}{V}(S_f - S) - \frac{1}{Y_{X/S}}r_g(X,S) \\ \end{align*}

The formula derivation from basic ODE correlation of X, P, and S with V is described in Bioreactor page.

The specific cell growth rate expressed as follows:

$\mu_{S} = \mu_{max}\frac{S}{K_S + S}$

However, that equation doesn’t include the inhibitor variable. Based on ferrous iron oxidation in T. Ferrooxidans periplasm, we hypothesized that the cell regulates itself to inhibit the oxidation process with competitive inhibition. Thus, the µs equation with competitive inhibition would be:

$\mu_{S} = \mu_{max}\frac{S}{K_S(1+ \frac{S_0-S}{K_P})+S}\\$

XCell concentration (grams/liter)
PProduct concentration (grams/liter)
SSubstrate concentration (grams/liter)
$$S_0$$      Feed substrate concentration (grams/liter)
$$r_c$$ Rate of cell or biomass growth (grams/liter/hr)
$$r_p$$ Rate of product formation (grams/liter/hr)
$$Y_{P/X}$$ Product yield coefficient
$$Y_{S/X}$$ Cell growth coefficient
F(t)Flowrate (liter/hr)
$$\mu_{S}$$Cell specific growth rate
$$\mu_{max}$$Maximum cell specific growth rate
$$K_S$$Half saturation constant
$$K_P$$Inhibition constant

The regulation of T. Ferrooxidans cell as competitive inhibitor assumption gives the following value for bioreactor-related variables:

Kinetic VariableFIELD2
$$\mu_{max}$$ 1.25 hr^{-1}
$$K_S$$ 0.125 gram/liter
$$K_P$$ 0.02 gram/liter
$$Y_{S/X}$$ 0.324 gram/gram
$$Y_{P/X}$$ 0.4348 gram/gram
$$S_0$$ 606.6225 gram/liter
Bioreactor Operating Condition
pH2.3
Temperature25oC

Suzuki et al, 1989

###### Cell-Regulated Competitive Inhibition

In this model, it is assumed that the operating conditions are maintained from the beginning until the end of the process batch. Then, the reaction was modelled via Jupyter Notebook.

These are the results of cell concentration (X), product concentration (P), substrate concentration (S), and total volume (V) in 40 hours batch time:

VariableInitial (time = 0 hr)End of Batch (time = 40 hr)Unit
X0.051,95E+10Gram/liter
P08.463.576.051Gram/liter
S6.066.225583.712.039Gram/liter
V2527Liter

To further extract the product, the residence time could be extended to be more than 40 hours. However, it is not efficient. The 40 hours residence time is considered to be moderate time for the gold bioleaching step, but it is still time consuming.

There is another component that potentially be the most potent inhibitor agent: Fe3+, the product itself. The product being the inhibitor itself is present in a lot of other enzymatic reactions. This hypothesis could also be true for pyrite dissolutions.

The second theory is by assuming Fe3+, the product itself to be the one who inhibits the iron oxidation process with a competitive inhibition mechanism. This theory gives the following value for bioreactor-related variables:

Kinetic VariableFIELD2
$$\mu_{max}$$  1.25 hr^{-1}
$$K_S$$  0.048 gram/liter
$$K_P$$  0.06 gram/liter
$$Y_{S/X}$$  0.324 gram/gram
$$Y_{P/X}$$  0.4348 gram/gram
$$S_0$$  606.6225 gram/liter
Bioreactor Operating Condition
pH1.6
Temperature30oC

Jones et al, 1983

In this model, it is assumed that the operating conditions are maintained from the beginning until the end of the process batch. Then, the reaction was modelled via Jupyter Notebook.

These are the results of cell concentration (X), product concentration (P), substrate concentration (S), and total volume (V) in just 4 hours batch time, 36 hours shorter than the previous mechanism:

VariableInitial (time = 0 hr)End of Batch (time = 4 hr)Unit
X0.051,96E+10Gram/liter
P08.508.121.549Gram/liter
S6.066.225267.506.643Gram/liter
V2525.2Liter

We managed to achieve 4 hours of pyrite dissolution process with this mechanism. Still, there are another inhibition mechanism by Fe3+ that is worth modeling.

There are 3 types of enzyme inhibition: competitive, non-competitive, and uncompetitive. According to one research, Fe3+ also has potential to inhibit the iron oxidation in T. Ferrooxidans by non-competitive inhibition.

The third theory is by assuming Fe3+ to inhibit the reaction by non-competitive inhibition mechanism. For non-competitive inhibition, the value of cell specific growth rate would be described as:

$\mu_{S} = \mu_{max}\frac{S}{K_S+S(1+ \frac{S_0-S}{K_P})}\\$

This theory gives the following value for bioreactor-related variables:

Kinetic VariableFIELD2
$$\mu_{max}$$   µmax1.33 hr^{-1}
$$K_S$$   KS0.04 gram/liter
$$K_P$$   KP0.06 gram/liter
$$Y_{S/X}$$   YS/X0.324 gram/gram
$$Y_{P/X}$$   YP/X0.4348 gram/gram
$$S_0$$   S0606.6225 gram/liter
Bioreactor Operating Condition
pH1.6
Temperature30 oC

Jones et al, 1983

###### Fe3+ Non-Competitive Inhibition

In this model, it is assumed that the operating conditions are maintained from the beginning until the end of the process batch. Then, the reaction was modelled via Jupyter Notebook.

These are the results of cell concentration (X), product concentration (P), substrate concentration (S), and total volume (V) in 17000 hours:

VariableInitial (time = 0 hr)End of Batch (time = 17000 hr)Unit
X0.051,97E+10Gram/liter
P08.545.373.735Gram/liter
S6.066.2250.0307269178Gram/liter
V25875Liter

From the non-competitive hypothesis, the residence time increased to 1700 hours. The volume also expanded 35 times.

Thus, the best inhibition mechanism to determine bioreactor for artisanal & small scale gold mining (ASGM) is competitive inhibition by Fe3+.

#### References

1. KEGG, ENZYME: 1.4.99.5 [Online] https://www.genome.jp/dbget-bin/www_bget?ec:1.4.99.5 [accessed on July 28, 2021 at 09:06 WIT]
2. Kianinia, Y., Khalesi, M.R., Abdollahy, M., Hefter, G., Senanayake, G., Hnedkovsky L., Darban, A.K., Shahbazi, M., 2018, Predicting Cyanide Consumption in Gold Leaching: A Kinetic and Thermodynamic Modeling Approach, Minerals, vol. 8, no. 3, pp. 110. doi:10.3390/min8030110
3. Michaels, R. and Corpe, W.A., 1965, Cyanide Formation by Chromobacterium violaceum, Journal of Bacteriology, vol. 89, no. 1, pp. 106-112.
4. Brazilian National Genome Project Consortium, 2003, The complete genome sequence of Chromobacterium violaceum reveals remarkable and exploitable bacterial adaptability, Proceedings of the National Academy of Sciences of the United States of America, vol. 100, no. 20, pp. 11660-11665.
5. Lobell, R.B., Schleif, R.F., 1990, DNA looping and unlooping by AraC protein, Science, vol. 250, no. 4980, pp. 528-532.
6. Schleif, R., 2003, AraC protein: a love-hate relationship, Bioessays, vol. 25, pp. 274-282.
7. Brazilian National Genome Project Consortium, 2003, The complete genome sequence of Chromobacterium violaceum reveals remarkable and exploitable bacterial adaptability, Proceedings of the National Academy of Sciences of the United States of America, vol. 100, no. 20, pp. 11660–11665. doi:10.1073/pnas.1832124100
8. Cuthbertson, L., Nodwell, J.R., 2013, The TetR family of regulators, Microbiology and molecular biology reviews: MMBR, vol. 77, no. 3, pp. 440–475. doi:10.1128/MMBR.00018-13
9. Anonymous, 2011, Registry of Standard Biological Parts - pBAD Strong Characterization, Last Updated 22 Oct 2009) [Online] http://partsregistry.org/Part:BBa_K206000:Characterization [accessed on August 21, 2021 at 08:34 WIT]
10. Schleif, R 1969, Induction of the L-arabinose operon, Journal of Molecular Biology, Vol 46, pp. 197-199. doi: 10.1016/0022-2836(69)90066-7
11. Hendrickson, W.., Schleif, R., 1984, Regulation of the Escherichia coli L-Arabinose Operon Studied by Gel Electrophoresis DNA Binding Assay, Journal of Molecular Biology, Vol 174, no. 3, pp. 611-628.
12. Koostra, A.M.J., Mosier, N.S., Scott E.L., Beeftink, H.H., Sanders, J.P.M., 2009 , Differential effects of mineral and organic acids on the kinetics of arabinose degradation under lignocellulose pretreatment conditions, Biochemical Engineering Journal, Vol 43, pp. 92-97.
13. Bionumbers Logo, Key Numbers for Cell Biologists [Online] http://bionumbers.hms.harvard.edu/KeyNumbers.aspx [accessed on October 22, 2021]
14. Mintel, R. and Westley, J., 1966, The Rhodanese Reaction : Mechanism of Thiosulfate Binding, The Journal of Biological Chemistry, vol. 241, no. 14, pp. 3386 – 3389.
15. Seo, J.M., 2013, Characterization of an L-Arabinose Isomerase from Bacillus Thermoglucosidasius for D-Tagatose Production, Bioscience, Biotechnology, and Biochemistry, vol. 77, no. 2, pp. 385-388.
16. Megerle, J.A., Fritz, G., Gerland, U., Jung, K., Radler, J.O., 2008, Timing and Dynamics of Single Cell Gene Expression in the Arabinose Utilization System, Biophysics Journal, vol. 95, no. 4, pp. 2013-2115. doi: 10.1529/biophysj.107.127191
17. King, Z.A., Lu, J.S., Dräger, A., Miller, P.C., Federowicz, S., Lerman, J.A., Ebrahim, A., Palsson, B.O., and Lewis, N.E., 2016, BiGG Models: A platform for integrating, standardizing, and sharing genome-scale models, Nucleic Acids Research, vol. 44, no. D1, pp. D515-D522. doi:10.1093/nar/gkv1049
18. KEGG, Metabolic pathways - Reference pathway [Online] https://www.genome.jp/kegg/pathway/ [accessed on October 22, 2021]
19. The Regents of the University of California, 2019, ESCHER: Build, share, and embed visualizations of metabolic pathways [Online] https://escher.github.io/#/ [accessed on October 22, 2021]
20. Tay, S., Natarajan, G., Rahim, M., Tan, H., Chung, M., Ting, Y. et al., 2013, Enhancing gold recovery from electronic waste via lixiviant metabolic engineering in Chromobacterium violaceum, Scientific Reports, vol. 3, no. 1.
21. Banerjee, D., Raghunathan, A., 2019, Constraints-based analysis identifies NAD+ recycling through metabolic reprogramming in antibiotic resistant Chromobacterium violaceum, PLOS ONE, vol. 14, no. 1, pp. e0210008.
22. Patil, K. R., Rocha, I., Förster, J., Nielsen, J., 2005, Evolutionary programming as a platform for in silico metabolic engineering, BMC Bioinformatics, vol. 6, no. 308. doi:10.1186/1471-2105-6-308
23. Burgard, A.P., Pharkya, P., Maranas, C.D., 2003, OptKnock: A Bilevel Programming Framework for Identifying Gene Knockout Strategies for Microbial Strain Optimization, Biotechnology and Bioengineering, vol. 84, no. 6, pp. 647-657.