Modelling
When facing an engineering challenge, models are our best tool to simulate and predict the outcome of our designs before even they becoe a physichal reality. Year by year modelling is gaining more importance for the synthetic biology community, since they allow us to speed up the development process by guiding our design decissions.
In this page you will find a brief summary of how biological modelling has been utilized within 4C_Fuels Project. However, we have also applied modelling for simulating an integrated bioprocess design. Further details about bioprocess modelling can be found at the Implementation and Proof Of Concept
Genome-Scale Modelling of PCC11801
During our project, we thought it might be useful to generate a Genome-Scale Model (GEM) of our chassis: Synechococcus elongatus PCC11801. After performing a comparison between both genomes and annotated features, we realized that most of the central metabolic reactions of both strains were almost identical. Main genetic differences are due to the presence of environmental-stress protective genes, where PCC11801 has a richer pool of genes dedicated to cope with metal, high temperatures or xenobiotics. Only small mutations in the ATPase and other central energy metabolism genes are found. Despite these differences does not change the general stoichiometric configuration of the organism metabolic network, previous research has stated that they are critical for conferring the strain a fast-growth phenotype. Then, we thought that most metabolic network should should be conserved across our chassis and other more conventional laboratory strains.
To achieve the task of adapting a model from another model organism to our chassis, we employed a previously existing model for Synechococcus elongatus PCC7942, the iJB785 Model. Thanks to bioinformatic tools, following a gap-filling approach, we have been able to create a genome-scale model of PCC11801 based on the available iJB785 model of the closely related PCC7942 strain. This new draft genome model contained only five reactions less that were present in the reference strain but not in PCC11801. This confirmed our initial guess that their metabolic network is highly conserved, and how the enhanced metabolic performance of PCC11801 could be derived from a distinctive network regulation, as well as mutations leading to higher active enzymes involved in carbon fixation and energy metabolism.
During our project, we thought it might be useful to generate a Genome-Scale Model (GEM) of our chassis: Synechococcus elongatus PCC11801. After performing a comparison between both genomes and annotated features, we realized that most of the central metabolic reactions of both strains were almost identical. Main genetic differences are due to the presence of environmental-stress protective genes, where PCC11801 has a richer pool of genes dedicated to cope with metal, high temperatures or xenobiotics. Only small mutations in the ATPase and other central energy metabolism genes are found. Despite these differences does not change the general stoichiometric configuration of the organism metabolic network, previous research has stated that they are critical for conferring the strain a fast-growth phenotype. Then, we thought that most metabolic network should should be conserved across our chassis and other more conventional laboratory strains.
To achieve the task of adapting a model from another model organism to our chassis, we employed a previously existing model for Synechococcus elongatus PCC7942, the iJB785 Model. Thanks to bioinformatic tools, following a gap-filling approach, we have been able to create a genome-scale model of PCC11801 based on the available iJB785 model of the closely related PCC7942 strain. This new draft genome model contained only five reactions less that were present in the reference strain but not in PCC11801. This confirmed our initial guess that their metabolic network is highly conserved, and how the enhanced metabolic performance of PCC11801 could be derived from a distinctive network regulation, as well as mutations leading to higher active enzymes involved in carbon fixation and energy metabolism.
Assesing Metabolic Feasibility of our Design
During our project, we thought it might be useful to generate a Genome-Scale Model (GEM) of our chassis: Synechococcus elongatus PCC11801. After performing a comparison between both genomes and annotated features, we realized that most of the central metabolic reactions of both strains were almost identical. Main genetic differences are due to the presence of environmental-stress protective genes, where PCC11801 has a richer pool of genes dedicated to cope with metal, high temperatures or xenobiotics. Only small mutations in the ATPase and other central energy metabolism genes are found. Despite these differences does not change the general stoichiometric configuration of the organism metabolic network, previous research has stated that they are critical for conferring the strain a fast-growth phenotype. Then, we thought that most metabolic network should should be conserved across our chassis and other more conventional laboratory strains.
To achieve the task of adapting a model from another model organism to our chassis, we employed a previously existing model for Synechococcus elongatus PCC7942, the iJB785 Model. Thanks to bioinformatic tools, following a gap-filling approach, we have been able to create a genome-scale model of PCC11801 based on the available iJB785 model of the closely related PCC7942 strain. This new draft genome model contained only five reactions less that were present in the reference strain but not in PCC11801. This confirmed our initial guess that their metabolic network is highly conserved, and how the enhanced metabolic performance of PCC11801 could be derived from a distinctive network regulation, as well as mutations leading to higher active enzymes involved in carbon fixation and energy metabolism.
Collaboration with MiamiU_OH Team
In the beginning of the summer, we met with the MiamiU_OH iGEM team, which is when we realized that not only both of our projects shared a similar motivation, but also how we have planned to address the scientific challenges in a very similar way. After speaking about our project ideas, we found how part of our metabolic design could be also used to enhance carbon capture and then we moved forward to think about collaborating together. MiamiU_OH iGEM team Team offered us the possibility to test our designed metabolic strategies, utilizing our adapted model and the reaction schemes we were planning to implement. SInce they were already evaluating via FBA their own strategies, and one of our metabolic designs could offer a benefit for their goal too, we thought it was a great idea. Then, we prepared all the required information and submitted it to the MiamiU_OH iGEM team team, who kindly studied the metabolic feasibility of our design and remitted it back to us.
The process of FBA and results received from MiamiU_OH iGEM team are explained below. However, to properly understand them, we encourage you to read first our Metabolic Engineering Design page. The general strategy consisted of evaluating simultaneously the performance of the organism in terms of product generation and growth rate. A total of 4 models were used. The adapted PCC11801 model (WT model). One model which only includes the reactions allowing for butanol synthesis (BUT model). A second model including the reactions for the phosphoketolase (PK) bypass and MCG carbon fixation pathway (MAL model). And a last one where the 3 pathways (butanol synthesis, PK bypass, MCG pathway) areimplemented (BUTMAL model). It is important to note that a sink reaction was also added for n-butanol in the BUT and MALBUT models. For n-butanol and acetyl-CoA respectively. These reactions were included for the verification of a nonzero flux for all added reactions increasing acetyl-CoA production. For growth assessment and subsequent analysis, this sink reaction was removed.
Verification of Models
First, all models must be able to exhibit proper growth estimations and to perform their expected functions based on their previously tested legitimacy in vivo. This means that some metabolic fluxes have to be fixed during the simulation in order to obtain a numeric result for the flux distribution in the network. Since we did not have time to obtain experimental data of carbon flux or product generation in our strain, all the calculated fluxes would be a mere theoretical estimation. However a comparison of fluxes among the different models can be performed, using the relative increment or reduction as a proof of its feasibility for enhancing or producing n-butanol. This way, the MALBUT should hypothetically increase butanol production compared with the MAL model, due to an increased pool of acetyl-CoA which will flow towards the butanol synthesis pathway. To determine that all added pathways are functional, the introduced metabolic reactions must be active, and the cell must be able to grow to some extent. To evaluate this, n-butanol and acetyl-CoA production reactions for the BUT and MAL models, respectively, were bound to have non zero fluxes. The model is forced to look for solutions where the synthetic pathways are active. Then, biomass production was optimized for all models; a nonzero solved optimized value would indicate a nonzero growth rate and therefore viable cells with active added pathways.
The WT model had a growth rate of 0.054 mgDW/h. This value corresponds with the baseline growth rate used for comparing the burden imposed by each artificial pathway. Then, the three adapted models were tested for biomass productivity. To do so, a sink reaction was added for the target metabolite, and a non-zero flux constraint was added for the optimization. This way, all the pathways are forced to be active.
The MAL model had a growth rate of 0.0470 mgDW/h
This result demonstrates that during PK and MCG pathways expression growth is not compromised. The reduced growth rate is a consequence of considering the “export” of acetyl-CoA in a sink reaction instead of its utilization.
The BUT model had a growth rate of 0.0423 mgDW/h.
This results demonstrates that growth is reduced due to carbon channeling towards a secreted product. However, the observed reduction in growth rate will derive from an excessive consumption of acetyl-CoA, whose formation in the WT may be limited.
The MALBUT model had a growth rate of 0.054 mgDW/h.
This growth is close to WT, indicating that the generated acetyl-CoA by phosphoketolase (PK) bypass and MCG carbon fixation pathway can efficiently supply the n-butanol biosynthesis, without compromising other metabolic functions. This allows the production of n-butanol while maintaining a higher growth rate
After seeing these results, we considered that our metabolic design could be viable and decided to focus our efforts within its biological implementation. However, now we have the resources to further explore the performance of our metabolic design.
Eventually, the MiamiU_OH iGEM team team performed a carbon fixation assessment of PK - MCG bypasses combined, to do so the models were optimized for both biomass production and the flux of the carbon fixing enzyme Rubisco. After assessing cellular viability RubisCO flux was assessed. The obtained results demonstrated how the combination of PK bypass and MCG cycle showed no observable improvement in growth rate and even slightly reduced the flux across RuBisCO. These results can confirm that both pathways actually require to be coupled with an acetyl-CoA consuming reaction in order to offer the expected carbon fixation improvement results.
Bioprocess Modelling
In order to evaluate the implementation of phototrophic microorganisms as a biomanufacturing platform, we decided to propose and evaluate an industrial bioprocess. In order to achieve this goal we have employed a simplified model to design and simulate the operations required for the industrial purification of the produced n-butanol.
Within this section you will find a more detailed explanation about the fundamental principles behind the mathematical models employed for downstream process design.
Designing an adsorption operation
At the industrial scale, adsorption operations are performed in multiple ways, however, the utilization of semi continuously operated fixed-beds are the most utilized ones. A fixed bed is a vessel, usually built in stainless steel whose interior is filled with an adsorbent material.
In order to design a fixed bed adsorption operation we have decided to outline the basic design parameters without entering into a deeply detailed design. Then, the design parameters to calculate belong to two groups: fixed-bed dimensions (length, diameter, required adsorbent quantity etc…) and operational parameters (adsorption flow and temperature, desorption flow and temperature, adsorption time and desorption time).
There exists multiple adsorption models which allow the simulation of adsorption operations. One of the simplest ones is based on the solute movement theory.
This approach considers that the equilibrium between the fluid phase and the adsorbent is reached instantaneously. Then product adsorption is determined by conditions of the fluid fed to the system, which can be easily determined with the adsorption isotherm of the fixed-bed filling. Then, the adsorption will be controlled either by the type of adsorbent, temperature, and composition of the stream from which the product has to be recovered. Briefly, solute movement theory established that as the fixed bed starts adsorbing the product, the first “layers” of the bed become saturated, and the more distant ones do not have adsorbed product. Then, a “saturation wave” can be imagined to be moving across the fixed bed, separating the “saturated bed zone” from the “regenerated bed zone”. The speed of this wave can be calculated via the material balance across a differential element of the fixed bed, leading to the equation 1.
A similar approach can be used for bed hot air bed heating during desorption, The hot air flow across the bed will generate a temperature profile, then a thermal wave will be also produced, which will separate the regions of the bed where the bed filling still holding the product of interest adsorbed, and those which has been already desorbed. Then an energy balance can be performed in order to determine the required time for column heating.
This equation is a simplified version which does not consider the heat requirements of the fixed-bed walls or heat losses in the environment. The “temperature wave” speed can be used to estimate the required time for bed heating.
However, for designing the fixed bed It is important to also consider that the application of solute movement theory is limited to certain conditions, such as the operation of the fixed bed in plug flow conditions. and neglectable mass and heat transfer limitations.
Considering these aspects and after performing preliminary simulations, we discovered that a more accurate and simpler approach could be used. We analyzed the available n-butanol adsorption isotherms data. Then, we realized that for our particular case, n-butanol adsorption over silicalite could be simulated taking an even simpler approach. After comparing the simulation of the gas-phase adsorption results employing these concepts with a simpler approach, we decided to utilize a different modelling strategy for the next steps.
Based on former research results, we have confirmed that silicalite particles with diameters < 0.4 mm achieve fast enough mass transfer to consider that for any reasonable fluid flow during adsorption, column saturation waves establish as an abrupt front. In addition, silicalite adsorption capacity is significantly reduced with temperature, to the extent of being possible to consider that at 130 ºC the amount of adsorbed product is neglectable. Likesiwse adsorption isotherm is highly pronounced, allowing to achieve maximum saturation at concentrations as low as 1E-5 mol*m3.
This way, we decided to base our adsorption simulation in an already studied setup for n-butanol adsorption-desorption, maintaining those parameters that allow us to scale-up the design. Within adsorption operations, keeping the geometrical design of the adsorption equipment allows the scale up of the system.
Considering these and the formerly mentioned aspects, a modular system composed of multiple adsorption columns of fixed dimensions will be utilized. Each column will have 0.15 m length and 0.0105m inner diameter. This design allows having moderate pressure drops across the column while column heating for adsorbent drying can be easily achieved. Each one of these fixed-bed modules has capacity for 3.88 g of adsorbent, considering the silicalite particles density, porosity and the voids volume within the fixed bed system.
The amount of required adsorbent will determine the system dimensions in terms of the number of parallel modules to utilize.
Then, to calculate the required amount of adsorbent, it's capacity should be calculated. To do so Langmuir adsorption isotherm for n-butanol can be utilized.
Next step will be to define the flow and concentration required. While concentration ill depend on the product concentration within the culture media, flow will depend on the productivity of the system. Performing a mass balance for stationary conditions (detailed below) concentration and flow for the adsorption operation can be obtained.
Then, since adsorbent load in each adsorption module is fixed, the number of beds required will only depend on a fixed adsorption time. Taking into consideration that the system process should be operated continuously, two different adsorption units will be used. When the first is in the adsorption phase, the second one is under desorption, product recovery and regeneration cycle. The easiest way to achieve this is considering an adsorption time equal to the time required for drying, desorption and cooling.
Available scientific literature has shown that the optimal drying conditions of the fixed bed are 220 min at 80ºC with an air flow of 12 L*h-1. The optimal thermal desorption conditions maintain the air flow of 12 L*h-1 while the temperature is increased to 130ºC for 34 minutes. With this operation conditions the product can be efficently desorbed. Then a total desorption time of 270 minutes is considered.
Once total desorption time, product concentration in the feed stream, adsorbent load and flow are known, the required number of parallel modules can be easily calculated as follows:
Other downstream modelling calculations
Besides the design of fixed-bed adsorption systems, other calculations are required to study the purification system. The most relevant ones are explained below.
● Gas-phase n-butanol concentration. In order to calculate the amount of n-butanol evaporated during reactor aeration, it is necessary to determine the n-butanol concentration within the gas-phase. To do so, vapor-liquid equilibrium must be considered. Since n-butanol is highly diluted ,its concentration in the gas phase can be calculated using Raoult's law. (Equation 4), and after that refer it as molar fraction considering the pressure of the system (atmospheric in the normal case).
● Biorreactor Model. In order to determine the flow of culture that will exit the reactor during the continuous cultivation process, a simplified non-structured, non-segregated model of the system is employed. Biomass production is modelled using a logistic equation. With this approach the influence of light and carbon availability are not considered. Then we have utilized the specific growth rate and maximum biomass accumulation considering the available literature data for PCC11801 considering a CO2 concentration of 400 ppm and a light intensity of 400 μE, corresponding with sunny day conditions.
Carbon uptake has also been considered, in order to determine the required amount of CO2. To do so, an empirical carbon uptake coefficient is considered, which relates the CO2 fixation with the biomass production. Then, the following carbon uptake rate equation can be defined. Carbon uptake coefficient has been obtained from Kuan, D et. all publication.
This equation was utilized in combination with the biomass generation one to determine the range of viable aeration rates. However, due to the high variability in mass transfer depending on the utilized photobioreactor system, both biomass carbon uptake can vary widely. Then this equations has only been utilized to estimate the range of possible aeration rates
Eventually, productivity within the bioreactor has been referred to as a constant value. Under the assumption that during continuous photobioreactor operation, a well established stationary amount of cells will yield a constant production rate. Then, butanol production can be utilized to estimate the required exit flow of the reactor to achieve a constant n-butanol concentration.
To do so, a mass balance can be performed around the reactor, considering the generation of n-butanol and it’s extraction in the bioreactor discharge stream.
This modelling consideration of n-butanol production is not directly correlated with the biomass concentration. Instead, it defines constant productivity in order to consider the achievable productivity within the system in an independent manner to relative photobiocatalyst or cell concentration. Then , the formerly mentioned biomass production equation could eventually be used to calculate the amount of cells to recalculate within the photobioreactor to achieve a stationary culture density. For Synechococcus spp. and Synechocystis ssp a correlation of 0.6 g*L-1 of biomass per each OD730 unit can be utilized (Curtis, W. R et al.).
● Energy consumption calculations. To calculate the energy consumption of the auxiliary equipment required for photobioreactor and downstream process operation, different approaches have been used.
For pumps, the pressure drop of each element in the liquid circuit is considered. Bernoulli equation has been used to estimate required fluid pressurization in terms of kinetic energy and level differences. For pressure drop within the fixed bed Ergun’s equation considering the fluid superficial velocity is used. After that pumping power is calculated, being corrected with an efficiency factor of 60%. Then required pumping power has been calculated considering the time of adsorption stage.
Air compression has been modelled as isentropic compression with a isentropic coefficient of 1,4 for air and an compressor efficiency of 80%. This efficiency factor considers the energy transfer from electric power towards the pressurized fluid. Then required pumping power has been calculated considering the time of adsorption stage.
Aeration pressure drop has been estimated considering the hydrostatic pressure of the water column within the reactor. When necessary fixed-bed pressure drop has been estimated with the Ergun equation. Isentropic compression with an 80% efficiency has been also utilized. To these values, a 20% excess of pressure drop has been considered as friction losses in the pipes. Likewise, aeration consumption has been considered. Required aeration flow has been established as an common value within phototrophic cultivation installations: 0.03 vvm ( an air volume of 0.03 times the total volume of culture per minute).
Heat consumption has been calculated considering the specific heat of each material or components of the process stream. Likewise, for the condensation steps, the vaporization enthalpy of n-butanol has been used. In this case, a condenser energy efficiency of 0,33 has been considered, in order to account for the requirement of refrigeration for an adequate condensation within the heat exchanger. In the case of electric resistances, an efficiency factor for heat transport of 0,4 has been considered.
For centrifuges energy consumption, the bibliographic reference of 1.2 kWh*m-3 of liquid feed has been considered.
To read more about how all of these modelling concepts have been applied for downstream design, visit our Proof of Concept page.
Genome Scale Modelling
Jaiswal, D., Sengupta, A., Sohoni, S., Sengupta, S., Phadnavis, A.G., Pakrasi, H.B., Wangikar, P.P., 2018. Genome Features and Biochemical Characteristics of a Robust, Fast Growing and Naturally Transformable Cyanobacterium Synechococcus elongatus PCC 11801 Isolated from India. Sci. Rep. 8. https://doi.org/10.1038/S41598-018-34872-Z
Jaiswal, D., Sengupta, A., Sengupta, S., Madhu, S., Pakrasi, H.B., Wangikar, P.P., 2020. A Novel Cyanobacterium Synechococcus elongatus PCC 11802 has Distinct Genomic and Metabolomic Characteristics Compared to its Neighbor PCC 11801. Sci. Reports 2020 101 10, 1–15. https://doi.org/10.1038/s41598-019-57051-0
Schroeder, W.L., Saha, R., 2020. OptFill: A Tool for Infeasible Cycle-Free Gapfilling of Stoichiometric Metabolic Models. iScience 23, 100783. https://doi.org/10.1016/J.ISCI.2019.100783
Mendoza, S.N., Olivier, B.G., Molenaar, D., Teusink, B., 2019. A systematic assessment of current genome-scale metabolic reconstruction tools. Genome Biol. 20, 1–20. https://doi.org/10.1186/S13059-019-1769-1/FIGURES/7
García-Jiménez, B., García, J.L., Nogales, J., 2018. FLYCOP: metabolic modeling-based analysis and engineering microbial communities. Bioinformatics 34, i954–i963. https://doi.org/10.1093/BIOINFORMATICS/BTY561
Qian, X., Kim, M.K., Kumaraswamy, G.K., Agarwal, A., Lun, D.S., Dismukes, G.C., 2017. Flux balance analysis of photoautotrophic metabolism: Uncovering new biological details of subsystems involved in cyanobacterial photosynthesis. Biochim. Biophys. Acta - Bioenerg. 1858, 276–287. https://doi.org/10.1016/J.BBABIO.2016.12.007
Broddrick, J.T., Rubin, B.E., Welkie, D.G., Du, N., Mih, N., Diamond, S., Lee, J.J., Golden, S.S., Palsson, B.O., 2016. Unique attributes of cyanobacterial metabolism revealed by improved genome-scale metabolic modeling and essential gene analysis. Proc. Natl. Acad. Sci. U. S. A. 113, E8344–E8353. https://doi.org/10.1073/PNAS.1613446113/-/DCSUPPLEMENTAL
Superior, E.P., Manuel Muñoz López, F., Nogales, J., David, E., León, S., Centro, G., De Biotecnología, N., De Biología, D., Sistemas, D.E., 2021. Mejora en la generación automática de modelos metabólicos (GEM) de organismos relevantes en la industria alimentaria mediante gap filling.
Bioprocess Modelling
Seader, J. D., Henley, E. J., & Roper, D. Keith. (2011). Separation Process Principles with Applications using Process Simulators.
Kuan, D., Duff, S., Posarac, D., & Bi, X. (2015). Growth optimization of Synechococcus elongatus PCC7942 in lab flasks and a 2-D photobioreactor. The Canadian Journal of Chemical Engineering, 93(4), 640–647. https://doi.org/10.1002/CJCE.22154
Myers, J. A., Curtis, B. S., & Curtis, W. R. (2013). Improving accuracy of cell and chromophore concentration measurements using optical density. BMC Biophysics 2013 6:1, 6(1), 1–16. https://doi.org/10.1186/2046-1682-6-4
Cousin Saint Remi, J., Baron, G., Denayer, J., 2012. Adsorptive separations for the recovery and purification of biobutanol. Adsorpt. 2012 185 18, 367–373. https://doi.org/10.1007/S10450-012-9415-1
Incropera, F.P., DeWitt, D.P., Bergman, T.L., Lavine, A.S., 2017. Incropera’s principles of heat and mass transfer. Wiley 1000.
García Rodríguez, Á., García Rodríguez, Á., 2017. Recuperación de biobutanol para la producción de combustibles mediante ciclos de adsorción-desorción.